Phase-field modelling of brittle fracture in thin shell elements based on the MITC4+ approach

We present a phase field based MITC4+ shell element formulation to simulate fracture propagation in thin shell structures. The employed MITC4+ approach renders the element shear- and membrane- locking free, hence providing high-fidelity fracture simulations in planar and curved topologies. To capture the mechanical response under bending-dominated fracture, a crack-driving force description based on the maximum strain energy density through the shell-thickness is considered. Several numerical examples simulating fracture in flat and curved shell structures are presented, and the accuracy of the proposed formulation is examined by comparing the predicted critical fracture loads against analytical estimates.


Introduction
Thin shell structures find numerous applications in a wide range of industries within the aerospace, automotive, and construction sectors. Thin composite laminates in particular are being deployed in aircraft structures and comprise the chassis of automotive vehicles. Hence, providing highfidelity simulation of damage processes pertinent to thinshells becomes vital for estimating their critical load bearing capacities while at the same time reducing the number of high-cost experimental tests that are required for certification purposes.
Numerical simulation of evolving damage in thin shelllike structures is often performed using Reissner-Mindlin shell elements which allow efficient modelling of both in-plane (membrane) and out-of-plane (bending) deformations at a reduced computational cost. Especially when using an explicit time-integration scheme, shell elements do not penalize the stable time-increment even when the thickness is extremely small [51]. This makes Mindlin shells an ideal candidate for modelling computationally complex fracture problems involving, e.g., impact driven damage scenarios.
Damage modelling methods can be broadly categorized onto two types, i.e., Discrete or Smeared/Diffuse. In discrete methods, a crack is treated either explicitly as a geometrical entity or implicitly as a discontinuity in the displacement field. In diffuse methods, the crack is smeared over the surrounding domain and the stress degradation effects are incorporated by means of a damage variable embedded directly into the constitutive formulations.
Discrete crack approaches primarily rely on modifying an existing finite element mesh in the locations where crack propagates, see, e.g., the robust remeshing algorithms developed by Ingraffea and Saouma [33], Bouchard et al. [18], [19], Rethore et al. [47], Shahani and Fasakhodi [49]. The extended finite element method (XFEM), first introduced in Belytschko and Black [13], see, also [22], eliminates the need of expensive mesh-updating algorithms for tracking crack paths by decoupling the crack topology from the underlying finite-element mesh. The XFEM models cracks by introducing a set of additional (enriched) degrees of freedom and corresponding discontinuous basis functions. Over the past fifteen years, the method has evolved onto the industrial standard for resolving crack-tip stress singularities without the requirement of very fine discretizations. However, the XFEM is not free from computational complexities pertinent to the the number of additional DOFs; furthermore, it relies on the definition of ad-hoc assumptions vis-a-vis the stress field at the crack-tip. The extension of XFEM to 3-D problems is not straightforward and poses challenges in specifying the crack propagation increment in 3-D [25].
Cohesive Zone Modelling (CZM) is a discrete method [10,23,32] that simulates fracture propagation by redistributing the stresses ahead of the crack-tip over a finite fracture process zone (FPZ). The constitutive behaviour of the FPZ is defined on the basis of a traction-separation law. With the exception of the Cohesive Segments Method (CSM) [46], CZM relies on the pre-definition of the crack surfaces. Hence, it cannot predict arbitrary crack propagation scenarios and is mostly applied in cases where crack path is known a-priori, e.g., in composite delamination.
Diffuse damage modelling approaches such as the Phasefield method (PFM) [20,27] and the thick level set method [40], overcome these challenges and have been proven robust in treating complex crack patterns, e.g., branching, merging, and curvilinear crack paths. The PFM emerged from the step-changing works of Francfort and Marigo [27], Bourdin et al. [20] and has garnered much attention in the past 10 years. The main advantage of the PFM is that the crack initiation location and crack-paths do not need to be predefined, but naturally emerge from the solution of a PDE that is derived on the basis of energy-minimisation principles and solved over the entire computational domain. The PFM relies on replacing the sharp crack edges with a diffusive crack interface represented by the phase field and hence resolves difficulties of numerically tracking discontinuities in the displacement field during crack propagation. To this point, the PFM has been extended to treat brittle fracture [37,38,41], ductile fracture [4,17], hydraulic fracture [26,31,44,53], and has also been applied within material-point method (MPM) [34] and virtual-element method (VEM) setting [1].
Despite the significant advantages provided by shell elements in resolving three dimensional surfaces in a robust and efficient manner, there have been only limited efforts to apply the PFM for simulating shell damage problems; a detailed review is provided in [54]. The PFM has been used to model thin-shell fractures based on the Kirchoff-Love shell theory [7,35,52]. Kiendl et al. [35] adopted higher order smooth basis functions (NURBS), whereas Amiri et al. [7] employed maximum entropy meshfree approximations based on C 1 continuous basis functions. Reinoso et al. [45] extended the PFM for brittle fracture in large-deformation solid shell elements based on enhanced assumed strain (EAS) formulations.
An important challenge to address when using thin Mindlin shell elements is that they display membrane and transverse shear locking [36], which significantly affects the evolution the simulated crack path. Transverse shear locking occurs purely due to the displacement-based interpolation that is also used for the calculation of strains. This leads to a significant over-prediction of the bending stiffness and an under-prediction of the transverse deformations which may become lower than the theoretical estimates by orders of magnitude [24]. In addition, when the shell elements are curved or become overly distorted during nonlinear deformation, spurious coupling may occur between membrane and transverse shear strains; this also increases the element stiffness and leads to membrane locking [36]. Since in thin shells the membrane stiffness can be significantly larger than the bending stiffness, membrane locking leads to the exclusion of the desired bending modes from the overall element response [21].
To this point, several approaches have been proposed to alleviate locking in shell elements. Selective/reduced integration schemes have been employed [14,15,55], that however result in spurious zero energy modes necessitating additional hourglass stabilization techniques. More notably, the precise prediction of crack paths using elements based on reduced integration necessitates an even finer mesh discretisation in the critical regions which adds up to the computational complexity. The assumed strain approach based on the Mixed Interpolation of Tensorial Components (MITC) formulation proposed in the works of Dvorkin and Bathe [24], Bathe and Dvorkin [12], Bathe [11], and more recently the MITC4+ approach proposed by Ko et al. [36] has been successful in alleviating both transverse shear and membrane locking issues and also pass all basic patch tests in an optimal convergence behaviour for both uniform and distorted meshes.
In this work, we extend the phase-field modelling framework to simulate brittle fracture in MITC4+ based thin Mindlin shell elements, wherein damage initiates and evolves due to coupled membrane/bending deformations. We restrict our implementation to thin 4-noded shell elements subjected to small strain deformations; however, the approach is general and can be straight-forwardly extended to higher order shell elements. We use the proposed formulation to examine the post-fracture response of 3D surfaces and establish its accuracy by comparing against analytically predicted critical fracture loads.
The paper is structured as follows: In Sect. 2, the geometrical and kinematic considerations for the Mindlin shell element based on small-strain theory and coupled bending/membrane deformations are discussed. This is followed by a brief review of MITC4/MITC4+ formulations in Sect. 2.3. In Sect. 3 the combined constitutive relations extending brittle phase-field theory to MITC4+ shells are proposed, followed by numerical validations in Sect. 4.

Geometrical considerations
Point of departure for the formulation presented herein is the Reissner-Mindlin degenerated 4-node shell element [21]. The element comprises 6 local degrees of freedom (DOF), i.e., 3 translations and 3 rotations, as shown in Fig. 1.
The vector of the local nodal DOF at each node i is (Fig. 1b). The translational DOF, i.e., [u i , v i , w i ] are defined with respect to the global coordinate system x yz. The rotational DOF, i.e., [α i , β i , γ i ] are aligned with the local shell vectors, i.e., V 1i , V 2i , and V 3i , respectively. The vector V 3i is normal to the shell midsurface; the coplanar vectors V 1i , and V 2i are perpendicular to V 3i .
The coordinates of any arbitrary point x within the shell element are expressed in terms of the mid-surface nodal coordinates according to Eq. (1) where, t i is the shell thickness, and N i and x i = [x i y i z i ] T are the shape functions and coordinate vector for mid-surface nodes, respectively. Furthermore, ζ is the parametric coordinate along the thickness direction (ζ ∈ [−1, 1]), see, also, Fig. 1a.

Kinematics
The displacement at any point P lying above or below the shell mid-surface ( Fig. 1a) is derived with respect to the midsurface according to Eq. (2) [21].
where μ i contains the direction cosines of the shell vectors V 1i and V 2i and assumes the following form [Eq. (3)] The strain tensor [ε] xyz in the global cartesian system is defined according to Eq. (4) below.

Remark 1
The drilling DOF γ i have no stiffness associated with them. Hence, when coplanar elements share a common structural node, the drilling rotation about the shell normal V 3i at that node is not resisted and the system matrix becomes singular. On the contrary when not all elements surrounding a structural node are coplanar, the normal rotation of any element at the shared node has a component which gets resisted by the bending stiffness of adjacent elements. This means that in flat-shell geometries, the drilling rotation DOFs γ i can be omitted from the list of overall structural DOFs. However when the shell geometry is curved, any such suppression of γ i would lead to an over-constrained model and unwarranted stiffening of the structure [21]. Keeping this in view, in this work all 6 DOFs [u i , v i , w i , α i , β i , γ i ] are retained at nodes which are shared by non-coplanar elements; they are however omitted for nodes shared by coplanar elements.
To conveniently describe the kinematics of the shell element, the following coordinate systems are introduced (Fig. 2 Here, g i = x ,ζ i are the tangent vectors to the shell-surface at any arbitrary point with a position vector x, where ζ i ∈ {ξ, η, ζ } represents the parametric directions.

MITC4/MITC4+ formulations
In this section, the modified formulations for the transverse shear strain components based on the MITC4+ approach [24,36] are briefly presented. The 4-noded flat shell element shown in Fig. 2  In the original MITC4 formulations [24], the transverse shear strains ε st and ε rt are assumed constant along the edges perpendicular to the r and s axes, respectively (Fig. 3a). Furthermore, instead of using the displacement based interpolations shown in Eq. (4), the transverse shear strain components at any arbitrary point inside the element are interpolated based on the strain values at a pre-defined set of tying points {A, B, C, D} (Fig. 3a) using Eq. (5).
The transverse shear strains at these tying points, i.e., {ε st }, are calculated using the standard approach in Eq. (4) ε (T P) rt = (ε rt ) at TP using DI ε (T P) st = (ε st ) at TP using DI (6) where T P ∈ {A, B, C, D} denotes the tying points, and D I denotes the direct displacement-based interpolation analogous to Eq. (4). Similarly, in the MITC4+ formulations the membrane strain components {ε rr , ε ss , ε rs } are interpolated using the strains at the membrane tying points {A, B, C, D, E} shown in Fig. 3b. The detailed expressions are omitted herein and can be found in [36].

Coordinate transformations
To formulate the local element matrices and the constitutive relations, the strain tensor in Eq. (4) must be transformed into the shell-aligned local coordinate system [1, 2, 3] using the strain-transformation matrix T T T ε according to Eq. (7) [ε] 123 = ε 11 ε 22 ε 33 γ 12 γ 23 γ 13 A general definition forT T T ε involving strain-transformation between any two arbitrary coordinate systems is provided in "Appendix B" for completeness.
The assumed strains introduced in Eq. (5) are defined in the convected coordinate system [r , s, t], whereas the strains in Eq. (7) are expressed with respect to the shell-aligned local system [1,2,3]. Hence, to impose the MITC4+ modification, the shell-aligned local strains [ε] 123 must be first transformed into the convective strains [ε] rst . Due to the planar geometry of the 4-noded Mindlin shell elements, the in-plane directions for both coordinate systems [r , s] and [1, 2] are co-planar, but rotated with respect to each other. The rotation for transverse shear strains γ 13 , γ 23 T into the convected coordinates where In Eq. (9), α and β are the angles between the r and V 1 axes and s and V 1 axes respectively.
The in-plane convective strain components [ε rr , ε ss , γ rs ] is derived according to Eq. (10) where [ε] 123 is provided in Eq. (7). The transformation matrix T T T ε is directly derived from T T T ε in "Appendix B" using (a) (b) Fig. 3 Location of tying points used for assumption of a transverse-shear strains [24], b membrane strains within MITC4+ approach [36] only the elements of the 1st, 2nd, and 4th rows of T T T ε that correspond to the in-plane strain components [ε rr , ε ss , γ rs ].
After performing the MITC4+ modifications on the convective transverse shear strains {γ st , γ rt } and in-plane membrane strains {ε rr , ε ss , γ rs }, the total convected strain tensor [ε] rst is transformed back into the shell-aligned local coordinate system [ε] 123 by applying the inverse of the linear transformations shown in Eqs. (8)- (10).
The overall shell-aligned local strain tensor can then be expressed according to Eq. (11).
In the MITC4+ shell element, plane-stress assumptions hold, i.e. the out-of-plane tensile stress σ 33 = 0 in the shell-aligned local coordinate system [1,2,3]. Hence, the expression for the out-of-plane tensile strain ε 33 is derived according to Eq. (12) where ν is the material Poisson's ratio. with, is the global vector of DOF and the expressions for the direction cosines

Constitutive phase-field model
Griffith's theory of brittle fracture [30] derives from the assumption that the total potential energy of a fractured solid is additively decomposed into the bulk strain energy depending on the elastic deformations and the crack surface energy [Eq. (14)] In Eq. (14), and also Fig. 4, u is the displacement vector at any arbitrary point within the domain Ω, b and t represent the body forces within Ω and the surface-traction forces on the external boundary ∂Ω respectively, Γ c is the internal discontinuous boundary, ψ e is the elastic energy density and G c is the critical fracture energy density. The linearised strain tensor ε(u) is In the variational phase-field formulation, the sharp crack surface energy term in Eq. (14) is replaced by the regularized volume integral of a diffuse crack term shown in Eq. (16), i.e., where φ ∈ [0, 1] is the phase-field variable. For a quadratic fracture surface energy approximation introduced in Ambrosio and Tortorelli [5,6], the phase-field function γ (φ, ∇ φ) assumes the following form, i.e., where l o is the length-scale parameter controlling the width of phase-field diffusion zone. Using the functional definition of Eq. (17) it is straight-forward to show that φ = 0 and φ = 1 correspond to the fully-cracked and fully-intact states of the material, respectively. As a crack evolves, the elastic strain energy and induced stresses of the solid must decrease to compensate for the fracture energy required to generate new crack surfaces. This degradation mechanism is achieved by means of a degradation function g(φ) ∈ [0, 1] so that the elastic strain energy becomes Combining Eqs. (14)- (18), the following expression for the regularized potential energy of a cracked solid is obtained with u i , b i and t i as the vector components of displacement u, body-force b and surface traction force t respectively. Equation (19) corresponds to the phase-field model with an isotropic energy split; this however results also in cracks evolving under pure compression.
To address the issue of non-physical crack evolution under pure compression, phase-field models based on an anisotropic energy-splitting have been proposed, see, e.g., [3,8,38]. In the current work, we employ the spectral decomposition of the strain tensor as introduced in Miehe et al. [38] to facilitate comparisons with published results. To effectively impose plane-stress assumptions and calculate the in-plane and out-of-plane contributions of the strain energy density accurately, an additional 2-D strain tensor [ε] comprising only in-plane strain components [ε 11 , ε 22 , ε 12 ] is defined, i.e., The effective Cauchy stress vector is defined accordingly as Remark 2 To effectively impose the plane-stress assumption after damage has initiated, the in-plane membrane stress components [σ 11 , σ 22 , τ 12 ] T and their corresponding contributions to the total strain energy density must be calculated based on the 2-D strain tensor ε in Eq. (20), whereas the out-of-plane components [τ 23 , τ 13 ] T and their strain-energy contributions must be calculated using the complete 3-D strain tensor ε in Eq. (11). In addition, the out-of-plane tensile stress σ 33 can be explicitly set to zero to achieve optimal convergence characteristics and ensure that the plane-stress assumptions hold even after damage initiation.
According to Eqs. (24) and (25), only the positive tensile parts of the strain energy density and the Cauchy stress tensor, resepctively are multiplied by the degradation function g(φ). In this work, we employ the quadratic degradation function originally introduced in Pham and Marigo [42] and Miehe et al. [39], i.e., where the parameter η r was first defined in Ambrosio and Tortorelli [5] and denotes the residual stiffness to prevent ill-conditioning of system matrices when damage has fully propagated.
In the standard Mindlin shell theory, the transverse shear stresses along the shell thickness are not constant; rather they follow a parabolic distribution. To account for this effect, the transverse shear strains in Eq. (28) are scaled by a factor of 5/6 as also highlihgted in Cook et al. [21].
Based on the in-plane and out-of-plane contributions given in Eqs. (24) and (27), the overall tensile and compressive components of the total strain energy density can be given as in Eq. (30).
and hence, the expression for the total potential energy in Eq. (19) can be modified to naturally suppress crack growth in the regions under pure compression.
The strong form of the governing linear momentum and phase-field evolution equations are henceforth obtained by minimizing the total potential energy in Eq. (31) with respect to the field variables {u, φ}.
where the boundary conditions satisfy, with n i , i ∈ {1, 2, 3, . . . r } being the outward pointing normal vectors at the crack boundary.
To facilitate crack-irreversibility, a history variable (also referred to as crack-driving force D) proposed by [38], based on maximum strain energy density throughout the deformation history is adopted in the current formulations. The expression for D can be given as: Using the history variable to impose crack irreversibility produces acceptable and accurate results in cyclic loading scenarios. It must be emphasized however that it also disrupts the original variational formulation, see also [28,29] for alternative techniques to impose crack irreversibility.

Effective material tangent operator
The undamaged material elastic constitutive law for homogeneous materials is expressed in the local shell-aligned coordinate system [1,2,3] as where E = E/(1 − ν 2 ) with E and ν as Young's modulus and Poison's ratio respectively, and G = 0.5G/(1 + ν) is the shear modulus of the material [21]. Equation (35) is derived on the basis of a plane-stress state and indicates that the in-plane components of the elastic Cauchy stress [σ 11 , σ 22  To achieve optimal convergence rates even with the modified stress definitions in Eqs. (25) and (28), plane-stress assumptions must hold even when the material is undergoing damage. To achieve this, we consider a split of the damaged tangent stiffness matrix C d into its corresponding components as shown in Eqs. (36) and (39), which are based on in-plane {σ I P , ε } or out-of-plane {σ O P , ε} stresses and strains, respectively. where The in-plane material tangent operator [C C C d ] I P can also be represented as the 4 × 4 tensor shown in Eq. (38).
The out-of-plane component of material tangent operator can be similarly given as Eqs. (39) and (40).

Crack driving force variation along shell-thickness
The 3-D kinematics of Mindlin shell elements are defined with respect to the kinematics of the mid-surface. Furthermore, damage evolution as manifested by the evolution of the phase field is obtained only at the mid-surface nodes as a 2-D field. Hence, achieving an accurate and realistic stress degradation along the thickness becomes a challenging task [see, e.g., 35]. Driven by the observation that, especially in thin shell structures, the crack propagation through all thickness layers is often sudden and brutal, we employ a maximum through the thickness driving force rule to control the evolution of the phase field. Within this setting, the crack driving forces D at each integration point [see, Eq. (34) and Fig. 5 for the case of 3 thickness layers] is calculated based on the 3-D stress state within the shell volume. The crack-driving force at all thickness integration points corresponding to a particular mid-surface location is then set equal to the maximum of driving forces prevalent at those integration points. Finally, the phase-field evolution equation is integrated at each Gauss-point over the entire shell-element volume, thus causing phase-field (or damage) to evolve based on the max crack-driving force description.
This assumption captures the physical cracking phenomena through shell thickness very well and leads to highly accurate critical fracture strength predictions, especially during bending dominated failure scenarios, as would be shown in the benchmark numerical examples.

Discretization and solution procedure
The coupled strong-form evolution equation (32) are discretized via a Galerkin approximation. The test S and weighting W function spaces for the displacement field are defined as The corresponding spaces for the phase field are Multiplying the strong form equation (32), integrating by parts and performing the necessary algebraic manipulations eventually leads to the the following convenient nodal residual form for the equilibrium equation at node i, and the phase-field evolution equation at node i respectively. In Eqs. (46) and (47), V is the element volume, N i is the 2-D shape function and B u i is the straindisplacement matrix as expressed in Eq. (4), and T T T rot , T T T ε are the rotation and strain-transformation tensors defined in Eqs. (13) and (7), respectively, which facilitate the calculation of the internal forces F F F u int in the local shell coordinate system [1, 2, 3] and their subsequent rotation into global [x, y, z] system. The explicit expressions for N i and B u i can be obtained from [21]. Furthermore, B φ i is defined with respect to shelllocal system [1, 2, 3] as shown in Eq. (48).

Remark 3 In practice, the components of B
φ i can be effectively obtained by choosing the relevant components of locally transformed strain-displacement tensor T T T ε B u i . Since in Mindlin shell theory, the kinematics of the shellelement is represented using 2-D shape functions at the mid-surface, N i, 3 can be effectively set as zero.
Assembling the contributions from each element shown in Eqs. (46) and (47) into the overall residual vectors R R R u and R R R φ , the solution {u, φ} to the combined system of equations (32) can be obtained by setting R R R u → 0 and R R R φ → 0.
In the current work, the solution is obtained using the staggered or alternating minimization approach based on [38]. To ensure accuracy of the obtained solution, either both equations must be solved using staggered iterations [2] or the analysis must be solved using small incremental steps [38].

Integration procedure
For the MITC4+ shell element analyzed in the current work, a full-integration technique is employed with 4 Gauss integration points defined at each parametric thickness layer within the element. The integral expressions in Eqs. (46) and (47) are expressed in terms of parametric coordinates [ξ, η, ζ ] according to Eq. (49) (49) where I is evaluated at each integration point through the shell-volume and the definition for Jacobian [J] is provided in "Appendix A". The in-plane integration over {ξ, η} within each thickness layer ζ is performed using the Gaussintegration rule, where i ∈ {1, 2, 3, 4} are the in-plane integration points and w i ∈ {1, 1, 1, 1} are the weights associated with each of these points. The out-of-plane integration for all thickness layers is performed using the Simpson's rule, which can be expressed as in Eq. (51) for any integrand I .

Numerical examples
In all the test cases examined in this Section, a displacement controlled analysis has been employed. Unless explicitly stated, a one-pass staggered (alternating minimization) approach with a very small time-increment size ( 1.e −06 − 1.e −05 ) has been used for the solution of the coupled displacement-phase-field problem, and the residual stiffness η r is set to 0.

Notched square plate subjected to in-plane tension
The standard benchmark of the notched square plate shown in Fig. 6 under tension is examined herein. The material properties considered are E = 210 GPa, ν = 0.3, and G c = 0.0027 kN/mm. The mesh-size is h e = 0.0025 mm in the central strip where the crack is expected to propagate and the length scale parameter is l o = 0.0075 mm. A displacement control analysis is performed with an equilibrium tolerance of tol u = 1.e −08 .
It is interesting to note that the length-scale parameter l 0 adopted by Miehe et al. [38] is twice the size of l 0 used by Borden et al. [16]. This implies that the formulation detailed in [38] requires the minimum value of l 0 to be at-least twice  [16]. Indeed both the definitions of l 0 are equivalent, and one must be careful while appropriately choosing the value of l 0 when comparing results from the two formulations. The current work uses the formulations from [16], and hence the definition l 0 ≥ h e consistently hold for all the numerical simulations performed in this paper.
The resulting crack-path and load-displacement response are shown in Figs. 7 and 8

Notched square plate subjected to in-plane shear
The square plate specimen examined in Sect. 4.1 is subjected to horizontal in-plane tractions. Due to the nature of the loading and boundary conditions in this case (Fig. 9),

1-D beam subjected to transverse bending
A simply-supported rectangular plate subjected to a uniformly distributed pressure over the entire top face is considered as shown in Fig. 12. The aim of this example is to verify the proposed formulation predictions under bending-dominated fracture scenarios. The material and fracture properties are E = 1.e10 MPa, ν = 0, G c = 3 N/mm, and l o = 0.01 mm. The mesh is refined with h e = 0.003 mm in the entire mid-span of the plate where the crack propagation is expected. The thickness of the beam t = 0.01 mm is very small in comparison to the other two plate-dimensions (l = 8 mm and w = 1 mm) so that the effects of trans- verse shear and membrane locking on the critical fracture characteristics can be monitored. The vertical displacement is monitored at the centre-node of the plate, and the total applied distributed load is analysed with tol u = 1.e −06 . The crack initiates at the plate's mid-span which is also the location of maximum transverse deformation u z , as shown in Fig. 13. The load-displacement response is shown in Fig. 14 where a brittle fracture response under pure bending is indeed recovered. Since the Poisson's ratio is null, the transverse bending stiffness and the critical fracture loads should be identical to those predicted by the classical Euler/ Bernoulli beam theory. According to the Euler/ Bernoulli beam theory, the analytical elastic stiffness/length of the beam is established in Eq. (52) as where δ is the maximum transverse deformation obtained at the centre-span, E is the Young's modulus, I = wt 3 /12 is the area moment of inertia for the beam, and P = F/l is the total distributed applied load/length on the beam with units in N/mm, wherein F is the total applied load in N. For the current case, the analytical elastic stiffness of the beam can be calculated using Eq. (52) as k = P/δ ≈ 15.625 N/mm 2 . The slope of the predicted elastic loaddisplacement response in Fig. 14b (k = 0.06249/0.004 = 15.6225 N/mm 2 ) is in close agreement with this analytical estimate.
Considering the case of isotropic phase field fracture, i.e., fracture initiating both at tension and compression, the critical fracture load of the beam can be evaluated as where M cr is the critical bending moment required for crack initiation M cr = σ cr wt 2 /6 ( 5 4 ) and σ cr is the critical fracture stress. Based on derivations in [16], the critical fracture stress can be evaluated as Eq. (55).
For the given material and fracture properties, the critical stress in Eq. (55) is σ cr = 3.9775 × 10 5 N/mm 2 . This can be inserted into Eq. (53) to obtain the critical fracture load P cr = 0.8286 N/mm. Comparing the load-displacement responses in Fig. 14a, it is evident that the maximum crack-driving force description through thickness (detailed in Sect. 3.2) produces good agreement with the analytical fracture force estimated by Eq. (53) for the isotropic phase-field model. This reinstates the validity of the assumption that in thin shells, all transverse thickness layers at a given location would fracture simultaneously as soon as the crack is initiated in any one of these layers. Hence to incorporate this effect, the material stiffness degradation at that shell location must start as soon as the crack-driving force in any one of the associated thickness layers attains a critical limit. Such a description of crackdriving force D enables a 3-D description of crack topology and stress-degradation effects, albeit using a 2-D phase-field, refer to Sect. 3.2 for details.
Solving the phase-field evolution Eq. (32) using the spectral split proposed in [38] and with the same crack-driving force definition (Fig. 5) results in the load-displacement response in Fig. 14b. The corresponding critical fracture load is higher than the one provided by the isotropic model as the in this case material degradation occurs only on the part of the shell undergoing tension. The accuracy of the predicted critical force for the spectral-split case [38] is verified against the analytical estimates and XFEM results in Sect. 4.4.

Rectangular plate with a through crack subjected to pure bending moments
The rectangular plate specimen with a through crack shown in Fig. 15 is subjected to pure bending moments on its opposite edges and the accuracy of predicted peak moments are compared with the corresponding analytical values obtained using the stress-intensity factors in [50]. This example has been examined previously in Rouzegar and Mirzaei [48], where a comparison between SIFs obtained with XFEM and the analytical SIFs was performed. In this example, we attempt a comparison between the critical fracture loads predicted by the proposed phase-field model and the analytical formulations provided in Sih et al. [50]. The material properties are E = 210,000 MPa and ν = 0.33. The rotational increment Δθ X is monitored at the top-right corner node, and the plate is analysed with respect to varying sizes of Δθ X until the peak critical bending-moment is converged. An equilibrium tolerance of tol u = 1.e −06 is used in each case. According to [50], the analytical expression for the critical stress-intensity factor (SIF) for a centrallycracked plate with infinite width and subjected to remotely applied pure bending moment is evaluated as where K c is the equivalent critical SIF, t is the plate thickness, M 0,crit is the critical bending moment and a is half-length of the central crack. The analytical value of critical SIF for this example is provided in [48] as Assuming plane-stress conditions, the corresponding critical energy release rate G c is Substituting the value of K 1c from Eq. (57) into (56) and considering the edge length l = 70 mm, the critical bending moment/edge-length is derived as In our phase-field simulations, the mesh is refined in the central region with the element size h e = 0.25 mm where the crack is expected to propagate. The length-scale parameter and residual stiffness are chosen as l 0 = 0.25 mm and η r = 1.0e −3 , respectively. In the original variational formulation proposed by Bourdin et al. [20], it was shown that the fracture energy is overestimated depending on the size of finite element discretization. To compensate for this amplification, an effective critical energy release rate was proposed for the purpose of phase-field simulations, see also [43].  Fig. 16. The resulting crack topology is shown in Fig. 17. The crack originates simultaneously at both notch-tips and propagates horizontally towards the ends of the plate.
The results of a convergence analysis on the rotation increment Δθ X are also shown in Fig. 16. The converged value for the critical moment/length in Fig. 16 is M 0,P F M = 10.83 N mm/mm, which is in close agreement with the analytical bending moment/length derived in Eq. (59). This example further establishes the validity of assumptions made in Sect. 3.2 for the phase-field model based on anisotropic spectral strain decomposition, and verifies the accuracy of the proposed phase-field formulations in characterising realistic bending-dominated fracture scenarios.

Simply supported plate subjected to bi-directional bending loads
To demonstrate cracking phenomena under bi-directional bending loads, a simply supported plate with a uniformly distributed surface load is examined. The material and fracture properties are E = 1.9e 5 MPa, ν = 0.3, l o = 0.01 mm, G c = 0.295 N/mm, and the boundary conditions are as shown in Fig. 18. The mesh is refined along the plate's diagonals with h e = 0.005 mm. Only a quarter section of the plate is analyzed due to symmetry. The quarter-section is simply supported on the outer edges of the plate, whereas the internal shared edges are subjected to symmetric boundary conditions. A uniformly distributed load is applied over the entire top face until complete fracture of the plate, and the vertical displacement is monitored at the centre node of the plate. The analysis is run until a convergence tolerance of tol u = 1.e −06 is reached. The crack-path is shown in Fig. 19 which is consistent with the results reported previously in [9,35]. The loaddisplacement curve is illustrated in Fig. 20.

Cylinder with/without spherical closing cap subjected to uniform pressure loads
A cylindrical shell geometry with small axial notches placed on diametrically opposite ends and a uniformly applied pressure load on its inner surface is considered. Owing to the problem symmetry across the xy and xz planes, only the quarter part of the full cylinder is analyzed as shown in Fig. 21.
To examine the robustness of the approach, two different cases are considered, i.e. with and without a spherical cap at the two ends of the cylindrical shell. The latter is expected to give rise to crack branching at the spherical cap. The material and fracture properties are E = 7.0e 4 MPa, ν = 0.3, l o = 0.125 mm, G c = 1.5 N/mm. The mesh is refined with the size h e = 0.1 mm in all the cylindrical and spherical cap regions where the crack is expected to propagate. A displacement controlled analysis is performed with an equilibrium tolerance of tol u = 1.e −05 . For the cylinder specimen without spherical cap (Fig. 21a), the vertical circular arc BC is fixed along the x and z directions, whereas symmetric boundary conditions are imposed on horizontal edges AB, CD, and AD. The specimen with the spherical cap (Fig. 21b) is subjected to symmetric boundary conditions on all free edges, i.e. the vertical circular arc AD towards the notch is subjected to y-symmetric and horizontal edges AB, BC and CD are subjected to z-symmetric boundary conditions. The example demonstrates the capability of the proposed phase-field formulation to simulate damage in thin curved geometries that would otherwise demonstrate significant membrane as well as transverse shear locking.
The responses between the total applied pressure load and the displacement-norm measured at the notch-tip are compared in Fig. 22 for both the uncapped and capped specimens.
The crack-path at increasing load-increments for the uncapped and capped cylinders are shown in Figs. 23 and 24 , respectively. In the former case, the crack initiates at the notch-tip and propagates along the longitudinal direction of the shell. In the latter, the specimen demonstrates a similar response (Fig. 24), however, in this case the crack initiates at a slightly lower critical fracture load (Fig. 22). Over the spherical cap region, the crack first propagates linearly, but subsequently splits into two symmetric crack branches; these further evolve simultaneously.

Assymetric hyperboloid subjected to uniform internal pressure
To further demonstrate the robustness of proposed formulations in analysing curved shell problems, an assymetric hyperboloid geometry is considered which is subjected to a uniform internal pressure applied in the direction normal to its surface. The thin-shell assumptions apply as the thickness of the geometry t = 0.1 mm is significantly smaller than the other dimensions of the tower. A notch is introduced at the mid-height along the longitudinal direction of the shell. Due to the model symmetry only half part of the complete model as shown in Fig. 25 is analysed. To reduce the effect of bending at the boundary, the hyperboloid geometry is supported by an elastic shell structure, displayed as ABFE in Fig. 25 in which the evolution of phase-field (or damage) is restricted. Furthermore, the translational DOFs at the bottom-most part of the elastic base-support is completely fixed (u x = u y = u z = 0) while the rotational DOFs are kept free. For x + u 2 y + u 2 z measured at the notch-tip the curved side-edges BC and AD, z-symmetric boundary conditions are imposed whereas the top-edge CD is unrestrained. The internal distributed load is applied only on the hyperboloid region EFCD in the direction of outwardpointing normals to its surface. The elastic support ABFE is unloaded. The radial displacement is monitored at the bottom notch-tip shown by P in Fig. 25, and tol u = 1.e −05 . The crack initiates at the bottom notch-tip P as shown in Fig. 26, and propagates vertically downwards followed by a second branch that initiates at the top notch-tip Q. The two cracks propagate simultaneously and crack-branching is eventually observed at the bottom crack due to the shell-curvature at which point the shell loses all bearing capacity. The response between the vertical z-displacement at the bottom notch-tip P and the total applied load is shown in Fig. 27.

Conclusion
A phase-field driven shell element formulation is presented for brittle fracture in Reissner-Mindlin shells. We employ an MITC4+ approach to alleviate shear and membrane locking. Our method is based on the assumption of a maximum through the thickness crack driving force definition. Considering an anisotrpic split for damage evolution, we impose the plane stress assumptions directly on the tangent constitutive matrix; this approach has been found to provide optimum convergence rates. The accuracy of the proposed model is demonstrated by a set of illustrative numerical examples. Our solutions are verified against the analytical estimates both in the isotropic and anisotropic phase field case. The validity of the proposed model is further established by obtaining realistic and accurate fracture predictions in curved shell geometries, which display significant membrane and transverse shear locking due to the coupling of membrane and bending deformations. Whereas the proposed model highlights the capabilities brittle fracture phase field modelling to harness the advantages of MITC4+ formulations, research should be directed to account for more complex responses as,e .g., the case of finite strain ductile fracture. In the near future, we aim to extend the capabilities of the proposed phase-field model in simulating diverse anisotropic fracture scenarios. where where x = [x, y, z] is the position vector of any arbitrary point within the shell element, {ξ, η, ζ } are the shell parametric coordinates, t i is the shell thickness and {l 3i , m 3i , n 3i } are the direction cosines of normal vector V 3i to the shell midsurface at any node i.