Phase field modeling of brittle fracture in large-deformation solid shells with the efficient quasi-Newton solution and global–local approach

Available online at www.sciencedirect.com ScienceDirect Comput. Methods Appl. Mech. Engrg. 399 (2022) 115410 www.elsevier.com/locate/cma Phase field modeling of brittle fracture in large-deformation solid shells with the efficient quasi-Newton solution and global–local approach Z. Liua,b,∗, J. Reinosob, M. Paggia aIMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100, Lucca, Italy bDepartamento de Mecánica de Medios Continuos y Teorría de Estructuras, Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092, Seville, Spain Received 4 March 2022; received in revised form 13 July 2022; accepted 14 July 2022 Available online 4 August 2022 Abstract To efficiently predict the crack propagation in thin-walled structures, a global–local approach for phase field modeling using large-deformation solid shell finite elements considering the enhanced assumed strain (EAS) and the assumed natural strain (ANS) methods for the alleviation of locking effects is developed in this work. Aiming at tackling the poor convergence performance of standard Newton schemes, a quasi-Newton (QN) scheme is proposed for the solution of coupled governing equations stemming from the enhanced assumed strain shell formulation in a monolithic manner. The excellent convergence performance of this QN monolithic scheme for the multi-field shell formulation is demonstrated through several paradigmatic boundary value problems, including single edge notched tension and shear, fracture of cylindrical structure under mixed loading and fatigue induced crack growth. Compared with the popular alternating minimization (AM) or staggered solution scheme, it is also found that the QN monolithic solution scheme for the phase field modeling using enhanced strain shell formulation is very efficient without the loss of robustness, and significant computational gains are observed in all the numerical examples. In addition, to further reduce the computational cost in fracture modeling of large-scale thin-walled structures, a specific global–local phase field approach for solid shell elements in the 3D setting is proposed, in which the full displacement-phase field problem is considered at the local level, while addressing only the elastic problem at the global level. Its capability is demonstrated by the modeling of a cylindrical structure subjected to both static and fatigue cyclic loading conditions, which can be appealing to industrial applications. © 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Keywords: Global–local approach; Phase field fracture; Solid shell element; Quasi-Newton scheme 1. Introduction Thin-walled structures are widely employed in various industrial fields, and components of a large number of engineering systems are shell-like structures, such as photovoltaic modules [1,2], fiber-reinforced composite structures in the fields of astronautics and aeronautics [3,4], laminated glass structures of high-speed trains [5], ∗Corresponding author at: IMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100, Lucca, Italy. E-mail address: [email protected] (Z. Liu). https://doi.org/10.1016/j.cma.2022.115410 0045-7825/© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http: //creativecommons.org/licenses/by-nc-nd/4.0/). Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 among many other applications. The prediction of crack events in thin-walled structures has received a great deal of attention in the past decades due to its significant importance. As a consequence of complex interaction between kinematics, large deformation constitutive law and fracture description, computational fracture modeling of shelllike thin-walled structures is quite challenging and a great deal of research has been focused on the development of reliable numerical methods to deal with. With regard to computational methods for fracture and damage modeling of solids, several formulations have been developed based on different numerical techniques in the past decades [6–9]. Specifically for fracture modeling of thin-walled structures, the partition of unity concept has been proposed in [10], and most formulations assume the classic Kirchhoff–Love and Reissner–Mindlin shell models for the description of kinematics [11,12]. Besides, the so-called phantom node method was proposed in [13,14] to model arbitrary through-thickness cracks normal to the mid-surface of the element. Some other contributions [7,15] considered the use of non-local or gradient enhanced schemes to alleviate mesh dependence pathologies. Following this kind of scheme, the continuous kinematic formulation was incorporated into the solid shell formulation and embodied an extra degree of freedom in the center of the shell element in [16]. Though these techniques have been successfully applied for modeling of specific fracture events, their applicability for brittle fracture is generally limited. To overcome these limitations, the phase field approach proposed in [17,18] has emerged as a promising framework for brittle fracture modeling of solids. This methodology is established on the basis of the classical Griffith fracture criterion [19], in which failure occurs when the critical energy-release rate criterion is met. The essence of the phase field approach is to introduce a characteristic length scale within the variational formulation and idealize the sharp crack into a diffusive crack topology. Consequently, some complex fracture patterns such as crack branching and multiple crack coalescence can be captured through the finite element modeling in a natural way. This approach can also be categorized into the gradient damage methods since the discontinuities for the description of geometric and kinematic fields do not need to be considered explicitly due to the smeared representation of crack. It is worth mentioning that a thermo-dynamically consistent framework for phase field approach is established in [20,21], which provides an efficient solution scheme for treatment of diffusive crack. Further investigations exploit the extension of this approach to multi-physics [22,23], cracking in laminated structures [24] and dynamic analyses [25], demonstrating the great potential and advantages of this methodology in fracture modeling. Recently, the variational phase field methods have also been successfully extended to the prediction of fatigue crack growth. The J-integral law was combined with the phase field theory in [26] to capture the crack growth behavior under fatigue loadings. By the introduction of a fatigue function to degrade the material toughness, Carrara et al. [27] proposed a phase field framework to capture the fatigue behavior of brittle materials. This framework was recently extended to simulate the fatigue failure of a NiTi stent with the consideration of phase transformation [28], as well as hydrogen assisted fatigue [29]. Besides the fatigue analysis of linear elastic materials, the phase field methods have also been used to predict fatigue crack growth in elastic–plastic solids, see [30–32], among others. For the use of phase field approach combined with shell kinematics, first attempts have been proposed in [33] to model geometrically linear fracture problems, in which the developments are limited to standard finite elements according to Reissner–Mindlin theory. Another investigation in this regard is proposed in [34], which adopted Local Maximum-Entropy approximations and C1-continuous basis functions, making it hard to be adapted to standard finite element packages. The theory of performing tension–compression split through the thickness of shell and then integrating the energy contribution along the integration points has been proposed by Paul et al. [35] and Kiendl et al. [36], and later extended to the isogeometric analysis of multipatch structures [37,38]. Recently, this approach has been adopted by Kikis et al. [39] for phase field modeling of brittle fracture in Reissner–Mindlin shells and plates. A similar phase field formulation has also been developed by Pillai et al. [40] for brittle fracture modeling in a MITC4+ Reissner–Mindlin degenerated shell element. Alternative formulations have been proposed in combination of phase field approach with solid shell and Kirchhoff–Love formulation to model brittle and ductile fracture [41–44]. In the phase field research, great efforts have been devoted so far to the efficient solution of the coupled phase field-displacement governing equations. Due to the non-convexity of potential energy functional with respect to the displacement field and phase field, the Jacobian matrix in the Newton solution technique may become indefinite, leading to the convergence issues in the monolithic scheme, in which kinematic variables are solved simultaneously. To address this drawback, different solution methods have been developed, like line search algorithms [45], modified Newton method [46], and error-oriented Newton method [47]. Although good results have been obtained, the 2 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 convergence issues for the monolithic solution of phase field approach are still quite challenging as pointed out in [46]. On the other hand, a robust staggered scheme on the basis of alternating minimization has received much attention [48]. By fixing either the displacement field or the phase field, the functional will be convex with respect to the other kinematic variable. However, this staggered scheme is very computationally expensive, and convergence at critical loading increments requires a lot of iterations [45]. Besides, this approach is not unconditionally stable, and thus it is essential to apply very small increment step to ensure the accuracy of equilibrium solution [49]. To address the drawbacks of both staggered or AM and conventional Newton schemes, the robustness and performance of quasi-Newton methods in dealing with non-convex minimization problems have also been investigated, see [50,51]. The first attempt to study the feasibility of the quasi-Newton method in phase field analyses is proposed in [52] to solve the so-called unified phase field theory [53], showing the potential of this scheme in the context of phase field regularization of the cohesive zone model [54]. Recently, the quasi-Newton method is extended to the analysis of standard phase field formulation in the cases of different nature, which further emphasizes the power and promise of this scheme for the solution of coupled displacement-phase field problems [55]. To model brittle fracture in thin-walled structures, in this research, a phase field approach incorporated into the solid shell formulation with assumed strain methods to alleviate locking effects is developed. With regard to the solid shell formulations, the parametrization adopted in this work is based on the 6-parameter shell model proposed in [56]. The kinematics of solid shell element is described through parametrization of the top and bottom surfaces of the body with the displacement interpolation scheme, which features a similar discretization to the standard 8-node brick element [57,58]. Solid shell elements are widely used to discretize thin-walled structures owning to the several advantages over the conventional shell elements. For instance, only translational degrees of freedom are required to describe the kinematics, which circumvents the difficulties concerning nodal rotations in the element formulation, and consequently, it is very straightforward to incorporate any three-dimensional constitutive law and deliver more accurate strain and stress fields through the thickness direction [59,60]. To reduce effects of the locking pathologies, two numerical techniques, EAS [61] and ANS [62], are adopted in the solid shell formulation. In line with many shell formulations, the EAS method that considers the Hu–Washizu multi-field variational principle as departure point for its derivation is used herein to tackle volumetric and Poisson thickness locking effects, while the ANS method is applied to overcome transverse shear and trapezoidal deficiencies. This work aims at presenting a global–local phase field approach for large-deformation solid shells with the QN monolithic solution so as to provide a very efficient modeling framework for the prediction of local crack propagation in large-scale thin-walled structures. The two main innovative aspects of this work are: (i) the coupling of phase field approach for fracture modeling with solid shell formulation considering the EAS and ANS methods for the alleviation of locking effects, and its efficient and robust quasi-Newton monolithic implementation into the commercial finite element code ABAQUS, which was comprehensively evaluated by comparison with the standard staggered Newton scheme through different benchmark examples, including single edge notched tension and shear, fracture of cylindrical structure under mixed bending and tension, and computationally demanding fatigue induced crack growth; (ii) a specific global–local approach in the 3D setting tailored for phase field modeling with solid shell elements is proposed here to save computational cost in fracture modeling of large-scale thin-walled structures, and its capability was demonstrated by simulating the local crack growth of cylindrical structure under both static and fatigue cyclic loadings. 2. Phase field approach of brittle fracture In this section, a thermodynamically consistent phase field approach of brittle fracture, which is based on the Griffith and Irwin fracture theory recalling [20,21]. In the three-dimensional setting, let B0⊂R3and Bt⊂R3 denote the reference configuration and current configuration, and the corresponding position vectors are represented by Xand x, respectively. The body motion denoted by ϕ(X,t):B0×[0,t]→R3maps the material point X in the reference configuration onto its corresponding point xin the current configuration during the time interval [0,t]. Deformation gradient is defined as Fu:= ∂Xϕ(X,t), where ∂Xdenotes the partial derivative with respect to the reference position, and its determinant Ju=det[Fu] is the Jacobian of deformation. The essence of phase field approach to model brittle fracture is the regularization of sharp crack topology within a diffusive crack zone characterized by the scalar-valued function d:B0×[0,t]→[0,1]. This parameter smears out the sharp crack by a diffusive crack area of width l, as shown in Fig. 1. Note that the width of regularization region depends on the length scale l, which controls the transition between intact and broken parts. In this phase field modeling framework, 3 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 Fig. 1. Schematic representation of solid body with (a) sharp crack topology, and (b) phase field approximation of diffusive crack. Fig. 2. Definition of the cracking shell body in both the reference and current configurations with the solid shell concept and phase field approach for fracture modeling. dis defined as a smooth function of (X,t) in the reference configuration, and d=0 and d=1 represent the intact and cracked states, respectively. According to [63], the phase field equations are given by d−l2∆d=0 in B0and ∇Xd·n=0 in ∂B0(1) where ∆dis the Laplacian of phase field variable, and ∇Xdstands for the gradient of phase field in the reference configuration. The weak form reads ∫B0(1 lδd+l∇X(δd)· ∇Xd)dΩ=0,∀δd∈Vd(2) where Vo={δd∈H1(B0)|δd=0 on Γc}is the admission space of the phase field variable. To represent the spatial variation of total cracked surface, a crack surface density functional γ(d,∇Xd)is defined as Γc(d):= ∫B0 γ(d,∇Xd)dΩ(3) It should be pointed out that Γc(d) converges to the sharp crack surface when the value of characteristic length l approaches zero [17,18,64]. The crack surface density per unit volume of body γ(d,∇Xd)reads γ(d,∇Xd)=1 2ld2+l 2|∇Xd|2(4) In this regard, the crack energy integral defined on the sharp crack surface can be approximated by volume integrals, ∫Γc GcdΓ≈∫B0 Gcγ(d,∇Xd)dΩ(5) where Gcrepresents the fracture toughness of the body material based on the Griffith fracture theory. 4 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 3. Kinematics of solid shell and variational formulation 3.1. Solid shell formulation and kinematics In the solid shell concept, the position vector Xin the reference configuration, the position vector xin the current configuration and phase field variable dof any material point can be approximated by the corresponding vectors of material points on the bottom and top surfaces of shell element, and their representations take the form as follows X(ξ1, ξ2, ξ3)=1 2(1+ξ3)Xt(ξ1, ξ2)+1 2(1−ξ3)Xb(ξ1, ξ2)(6a) x(ξ1, ξ2, ξ3)=1 2(1+ξ3)xt(ξ1, ξ2)+1 2(1−ξ3)xb(ξ1, ξ2)(6b) d(ξ1, ξ2, ξ3)=1 2(1+ξ3)dt(ξ1, ξ2)+1 2(1−ξ3)db(ξ1, ξ2)(6c) where the parametric space is identified as: A:= {ξ=(ξ1, ξ2, ξ3)∈R3| −1≤ξi≤ +1;i=1,2,3}, the subscripts t and b represent top and bottom surfaces, respectively, and (ξ1, ξ2, ξ3)denote coordinates in three different directions of the parametric space, see Fig. 2. The kinematics of solid shell are described using the convective curvilinear coordinates, which is necessary due to the ANS interpolation for the transverse normal and shear strain components. The covariant tangent vectors Gi(ξ) in the reference configuration and gi(ξ) in the current configuration are defined as the partial derivatives of position vectors with respect to the convective coordinates ξi Gi(ξ):= ∂X(ξ) ∂ξi,gi(ξ):= ∂x(ξ) ∂ξi,i=1,2,3 (7) The contravariant basis vectors can be defined in a standard manner by Gi·Gj=δj iand gi·gj=δj i, and metric tensors are defined as G=Gi j Gi⊗Gj=Gi j Gi⊗Gj,g=gi j gi⊗gj=gi j gi⊗gj. In curvilinear setting, the deformation gradient Fuis defined as Fu=∂x ∂X=gi⊗Gi(8) where the Einstein summation convention on repeated indices is assumed. With the definition of metric tensor components Gi j =Gi·Gjand gi j =gi·gjin the reference and current configuration, the displacement derived Green–Lagrange strain tensor reads Eu:= 1 2[(Fu)TFu−I2]=1 2[gi j −Gi j ]Gi⊗Gj(9) where I2is the second-order identity tensor. The energetically conjugated second Piola–Kirchhoff stress is given by S=Si j Gi⊗Gj(10) where Si j represent the contravariant components. 3.2. Variational formulation In this section, the mixed Hu–Washizu variational principle incorporating the EAS and ANS methods to alleviate the locking pathologies is derived in combination with phase field model using solid shell element formulation. It is worth mentioning that the EAS method is adopted here to remedy volumetric and Poisson thickness locking, while the membrane and in-plane locking effects are alleviated by the ANS method [61,62]. In the large deformation setting with the consideration of enhancement based on EAS method, the Green–Lagrange strain is composed of displacement derived compatible strain Euand incompatible strain ˜ Efollowing the approach proposed in [56], and its complete form reads: E=Eu+˜ E. The multi-field Hu–Washizu variational functional of potential cracked body is given by Π(S,˜ E,u,d)=∫B0 g(d)Ψ(Eu,˜ E)dΩ+∫B0 Gcl 2(d2 l2+|∇Xd|2)dΩ−Πext (11) 5 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 where Ψrepresents the intact free energy function, g(d)=(1−d)2+Kstands for the monotonic degradation function, K≈0 is a positive parameter to ensure numerical stability in case of fully material degradation, and Πext identifies the external contribution. It is worth noting that the displacement u, the phase field variable d, and incompatible strain ˜ Eare the independent variables. Given the orthogonality condition between the second Piola–Kirchhoff stress and Green–Lagrange strain, the stress field is neglected here following the previous work [65]. The first variation of Eq. (11) with respect to the independent fields reads Ru(u, δu,˜ E,d)=∫B0 g(d)[∂Ψ ∂E:∂Eu ∂uδu]dΩ−δΠext(u)=0,∀δu∈Vu(12a) R˜ E(u,˜ E, δ ˜ E,d)=∫B0 g(d)[∂Ψ ∂E:δ˜ E]dΩ=0,∀δ˜ E∈V˜ E(12b) Rd(u,˜ E,o, δd)=∫B0 −2(1 −d)δdΨdΩ+∫B0 Gcl[1 l2dδd+ ∇Xd· ∇X(δd)]dΩ=0,∀δd∈Vd(12c) where Vu,V˜ Eand Vdare the admissible spaces of displacement, incompatible strain and phase field variable, respectively. 4. Finite element implementation 4.1. Finite element interpolation By means of isoparametric concept, the approximations of position vectors Xand x, displacement vector uand its variation δuon the solid shell element level are given by X=N˜ X(13a) x=N˜ x(13b) u=Nd, δu=Nδd(13c) where ˜ Xand ˜ xdenote the nodal position vectors in the reference and current configurations, respectively, and d represents the nodal displacement vector. The shape function matrix Nis defined as N=[N1,N2,N3,N4,N5,N6,N7,N8](14) where NI=diag [NI,NI,NI]with I =1,2,...,8, and its component NIis given by NI=1 8(1+ξ1 Iξ1)(1+ξ2 Iξ2)(1+ξ3 Iξ3)(15) with ξ1 I, ξ2 I, ξ3 I= ±1. In a similar manner, the phase field variable dand its variation δdare interpolated as d=Nd¯ d, δd=Ndδ¯ d(16) where ¯ drepresents the nodal phase field vector, and the shape function matrix Ndfor interpolation of phase field is given by Nd=[N1,N2,N3,N4,N5,N6,N7,N8](17) The material gradient of phase field ∇Xdand its variation ∇Xδdare interpolated as ∇Xd=G−T∇ξ¯ d=Bd(ξ)¯ d,∇Xδd=G−T∇ξδ¯ d=Bd(ξ)δ¯ d(18) where ∇ξrepresents the gradient with respect to the natural coordinates in the curvilinear setting, and Bdis the gradient interpolation matrix of phase field variable. The Green–Lagrange strain vector in this work is given by E=[E11,2E12,2E13,E22,2E23,E33]T. To overcome the deficiency of curvature thickness locking, the ANS interpolation proposed in [66] to modify transverse normal strain components E33 is adopted. This ANS method considers four collocation points defined in convective 6 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 Fig. 3. Position of collocation points in the element parametric space for the ANS method to alleviate the transverse shear locking and trapezoidal locking, respectively. coordinates ξCias ξC1=(−1,−1,0), ξC2=(1,−1,0), ξC3=(1,1,0), and ξC4=(−1,1,0), see Fig. 3. Besides, to prevent transverse shear locking in case of distorted element geometry, the ANS interpolation method proposed in [62] is also employed in this work. The four collocation points for the treatment of transverse shear strain components are ξA1=(0,−1,0), ξA2=(0,1,0), ξB1=(−1,0,0), and ξB2=(1,0,0), see Fig. 3. With respect to the ANS interpolations, the approximation of strain vector is expressed as Eu= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 2(g11 −G11) (g12 −G12) (1−ξ2)(gA1 13 −GA1 13 )+(1+ξ2)(gA2 13 −GA2 13 ) 1 2(g22 −G22) (1−ξ1)(gB1 23 −GB1 23 )+(1+ξ1)(gB2 23 −GB2 23 ) ∑4 i=1 1 4(1+ξ1 iξ1)(1+ξ2 iξ2)1 2(gCi 33 −GCi 33) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (19) where the superscripts A1, A2, B1, B2, and Ciwith i=1,2,3,4 denote the values at the corresponding collocation points. The approximation of virtual strain is given by δEu=Bδdwith B=[B1,B2,B3,B4,B5,B6,B7,B8](20) where the matrix for each node I considering ANS interpolation is defined as BI= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ NI,1(gT 1) NI,1(gT 2)+NI,2(gT 1) (1−ξ2)(NA1 I,1(gA1 3)T +NA1 I,3(gA1 1)T)+(1+ξ2)(NA2 I,1(gA2 3)T +NA2 I,3(gA2 1)T) NI,2(gT 2) (1−ξ1)(NB1 I,2(gB1 3)T +NB1 I,3(gB1 2)T)+(1+ξ1)(NB2 I,2(gB2 3)T +NB2 I,3(gB2 2)T) ∑4 i=1 1 4(1+ξ1 iξ1)(1+ξ2 iξ2)NI,3(gi 3)T ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (21) with NI,1, NI,2and NI,3defined as the derivatives of shape function NIwith respect to the parametric variables ξ1, ξ2, and ξ3, respectively. The enhancing strain field is interpolated by the operator M(ξ) that will be defined in Section 4.3 and it also depends on the number of enhancing modes ζtaken into account as pointed out in [58,67]. The interpolations of incompatible strain vector ˜ Eand its variation δ˜ Eare given by ˜ E≈M(ξ)ζ, δ ˜ E≈M(ξ)δζ(22) Inserting the discretization formulae into Eqs. (12), the discrete form of residual equations can be expressed as Rd(d, δd,¯ d,ζ)=∫B0 g(d)B(d)TSdΩ−Rd ext (23a) 7 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 Rζ(d,¯ d,ζ, δζ)=∫B0 g(d)M(ξ)TSdΩ(23b) Rd(d,d, δd,ζ)=∫B0 −2(1 −d)N(ξ)TΨ(d,ζ)dΩ+∫B0 Gcl[(Bd(ξ))T∇Xd+1 l2N(ξ)Td]dΩ(23c) Considering the irreversible growth of the fracture process, a history variable His introduced to modify the corresponding phase field residual vector Rd[68], which reads Rd=∫B0 −2(1 −d)N(ξ)THdΩ+∫B0 Gcl[(Bd(ξ))T∇Xd+1 l2N(ξ)Td]dΩ(24) In order to ensure the irreversibility of phase field variable, the history variable Hshould comply with the Kuhn–Tucker conditions Ψ+−H⩽0,˙ H⩾0,˙ H(Ψ+−H)=0 (25) where Ψ+is the decomposed strain energy density corresponding to tension according to the spectral split method and hybrid approach proposed in [21,48]. At the certain time t, the history variable can be written as H=max τ∈[0,t] Ψ+(τ) (26) 4.2. Approximation of the weak form and linearization For the multi-field problem, an iterative scheme is adopted to solve the set of nonlinear residual equations given by Eqs. (23). The consistent linearization of this system obtained from the concept of Gateaux directional derivative in matrix form can be written as ⎡ ⎣ kdd kdζ0 kζdkζζ 0 0 0 kdd ⎤ ⎦⎡ ⎣ ∆d ∆ζ ∆¯ d ⎤ ⎦=⎡ ⎣ Rd ext 0 0 ⎤ ⎦−⎡ ⎣ Rd Rζ Rd ⎤ ⎦(27) The different components of the stiffness matrix are defined as kdd =∫B0[(1 −d)2+K][Q+B(d)TCB(d)]dΩ(28a) kdζ=∫B0[(1 −d)2+K]B(d)TCM(ξ) dΩ(28b) kζd=∫B0[(1 −d)2+K]M(ξ)TCB(d) dΩ(28c) kζζ =∫B0[(1 −d)2+K]M(ξ)TCM(ξ) dΩ(28d) kdd =∫B0[Gc l+2H]N(ξ)TN(ξ)dΩ+∫B0 Gcl(Bd(ξ))TBd(ξ) dΩ(28e) The derivative of Bwith respect to din Eq. (28a) appears in order to compute ∆δEu:S, where ∆δEuis the linearized virtual strain tensor, and the matrix form is given by Q=∂B(d)T ∂dS=⎡ ⎢ ⎢ ⎢ ⎣ Q11 Q12 · · · Q18 Q21 Q22 · · · Q28 . . .. . ..... . . Q81 Q82 · · · Q88 ⎤ ⎥ ⎥ ⎥ ⎦(29) 8 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 where QIJ is defined as QIJ =diag [QIJ,QIJ,QIJ]for the combination of node I and J, and the scalar QIJ reads QIJ =ST ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ NI,1NJ,1 NI,1NJ,2+NI,2NJ,1 (1−ξ2)(NA1 I,1NA1 J,3+NA1 I,3NA1 J,1)+(1+ξ2)(NA2 I,1NA2 J,3+NA2 I,3NA2 J,1) NI,2NJ,2 (1−ξ1)(NB1 I,2NB1 J,3+NB1 I,3NB1 J,2)+(1+ξ1)(NB2 I,2NB2 J,3+NB2 I,3NB2 J,2) ∑4 i=1 1 4(1+ξ1 iξ1)(1+ξ2 iξ2)NI,3NJ,3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (30) where Sis the approximate stress field, which is organized in matrix form as S=[S11,S12,S13,S22,S23,S33]T. The isotropic Kirchhoff–Saint-Venant constitutive model is considered in this study, and the tangent material stiffness tensor Cin the convective curvilinear setting can be expressed as C=[λGi j Gkl +µ(Gik Gjl +Gil Gjk)]Gi⊗Gj⊗Gk⊗Gl(31) where λand µare the Lame parameters. The matrix form of stiffness Cis given by C= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ C1111 C1112 C1113 C1122 C1123 C1133 C1212 C1213 C1222 C1223 C1233 C1313 C1322 C1323 C1333 C2222 C2223 C2233 sym. C2323 C2333 C3333 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (32) Since inter-element continuity is not required, the enhanced strain term can be condensed out in the element level [56], and the system of equations that couples the kinematics with phase field is given by [k∗ dd 0 0 kdd ][ ∆d ∆d]=[Rd ext 0]−[Rd∗ Rd](33) where components of modified stiffness matrix and residual vector are defined as k∗ dd =kdd −kdζk−1 ζζ kζd(34a) Rd∗=Rd−kdζk−1 ζζ Rζ(34b) 4.3. Interpolation of the assumed strain field To alleviate volumetric and Poisson thickness lockings, the EAS method is incorporated into the solid shell formulation and briefly summarized in the following. It should be pointed out that the number of enhancing modes in the EAS method greatly influence the numerical efficiency of the formulation. In [67], seven enhancing modes are adopted to simulate the out-of-plane and membrane bending, and it is found that volumetric locking for incompressible analysis and Poisson thickness locking phenomena are effectively reduced, which is also comprehensively discussed in [58,69] with regard to the different enhancing schemes. Following this enhancing strategy, the interpolation matrix ˜ M(ξ) of the enhancing modes defined in the parametric space takes the form ˜ M(ξ)= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ξ10 0 0 0 0 0 0ξ20 0 0 0 0 0 0 ξ3ξ1ξ3ξ2ξ30 0 0 0 0 0 0 ξ1ξ2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (35) Given the consistency constraint of incompatible strain with respect to the stress field, the interpolation matrix M(ξ) defined in Eq. (22) is given by M(ξ)=[det J0 det J]T−T 0˜ M(ξ) (36) 9 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 Fig. 10. The contour plots of phase field values at the loading displacements 0.0085 mm, 0.010 mm, and 0.0135 mm during the crack propagation in the single edge notched shear. Fig. 11. The force versus displacement curves obtained from both the quasi-Newton monolithic and staggered schemes for the phase field modeling of single edge notched shear using solid shell formulation. Table 2 Computation times for different schemes in the single edge notched shear. AM 103inc. AM 2X103inc. AM 104inc. QN 102inc. CPU hours (h) 1.50 2.83 14.02 1.08 total number of iterations versus displacement curves obtained from both the quasi-Newton monolithic and staggered schemes for the phase field modeling of single edge notched shear test using solid shell element formulation are shown in Fig. 12. Difference of total number of iterations between the two solution schemes is smaller compared with the previous single edge notched tension, but remains significant in this case. To obtain the same result as the quasi-Newton monolithic solution, the required total number of time increments for the staggered scheme is still one order larger than the former scheme. The computational cost for the two different solution schemes is listed in Table 2. It can be concluded that the quasi-Newton monolithic scheme in the single edge notched shear is also much more efficient as its computational cost is around ten times smaller than that required in the staggered solution with 104number of increments. 16 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 Fig. 12. The total number of iterations versus loading displacement curves obtained from both the quasi-Newton monolithic and staggered schemes for the phase field modeling of single edge notched shear using solid shell formulation. Fig. 13. (a) Sketch of the cylindrical structure under mixed tension and bending and (b) finite element mesh of the specimen. 6.3. Cylindrical structure under mixed tension and bending The performance of the quasi-Newton monolithic scheme for phase field modeling using enhanced assumed strain shell formulation is further assessed in capturing fracture events of thin-walled structure with curved geometry under mixed tension and bending. In this section, we consider a quarter of cylindrical structure with an initial crack notch in the middle of the specimen, see Fig. 13(a). The dimensions of the specimen are as follows: the length is equal to 240 mm, the internal radius is 100 mm, the external radius is 100.1 mm, the width of the initial notch is 3 mm and its length is 6 mm. To comply with the symmetric boundary conditions, the translational degrees of freedom at the X = 0 side of the specimen are constrained in the X direction, and the Y = 0 side is constrained in the translational vertical direction Y. Besides, the translational degree of freedom in the Z direction at the back of the specimen is constrained, while the displacement condition δZis applied to the front edges and the displacement condition δY is applied to the notch edges, as shown in Fig. 13(a). The same material properties as in the previous cases are used here. Regarding the fracture properties, the phase field length scale is set to 4 mm and fracture toughness is 2.7 N/mm. The specimen is discretized using 8474 solid shell elements with refinement in the potential crack propagation region, see Fig. 13(b). The contour plots of phase field and displacement during the crack propagation are shown in Fig. 14, which agrees well with the fracture pattern predicted by an alternate formulation proposed in [42]. 17 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 Fig. 14. The contour plots of (a) phase field and (b) axial displacement at the loading displacements δZequal to 0.14 mm, 0.145 mm, and 0.155 mm during the crack propagation in the fracture modeling of cylindrical structure under mixed tension and bending. The force versus displacement curves obtained from both the quasi-Newton monolithic and staggered schemes with different increment sizes are shown in Fig. 15. As with the previous single edge notched tension and shear cases, reproducing the exact numerical result from the quasi-Newton monolithic solution requires a large number of time increments for the staggered scheme. For the phase field fracture modeling of cylindrical structure under mixed tension and bending, the resultant curve obtained from the staggered scheme with 104time increments agrees perfectly well with that solved by quasi-Newton monolithic scheme that requires only 100 time increments. The total number of iterations versus loading displacement curves obtained from both the quasi-Newton monolithic and staggered schemes are shown in Fig. 16. The same conclusion as the previous two cases can be drawn that the total number of iterations to obtain the same numerical result for the staggered scheme is 100 times larger than that required by the quasi-Newton monolithic scheme. Computation times for the two different schemes with different time increment sizes are shown in Table 3. It takes approximately 12 times longer using the staggered scheme with 104time increments than that solved by the quasi-Newton monolithic scheme. It should be pointed out that the quasi-Newton scheme is very robust and efficient for the phase field modeling with solid shell formulation considering both EAS and ANS compared with the conventional Newton solution schemes. 6.4. Phase field fatigue modeling using solid shell formulation In this section, the quasi-Newton monolithic scheme for phase field modeling with solid shell formulation is employed to address fatigue problems. The framework presented in [27] has been adopted and implemented into 18 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 Fig. 15. The force versus displacement curves obtained from both the quasi-Newton monolithic and staggered schemes for the phase field modeling of cylindrical structure under mixed tension and bending using solid shell formulation. Fig. 16. The total number of iterations versus displacement curves obtained from both the quasi-Newton monolithic and staggered schemes for the phase field modeling of cylindrical structure under mixed tension and bending using solid shell formulation. Table 3 Computation times for different schemes in the fracture modeling of cylindrical structure. AM 103inc. AM 2X103inc. AM 104inc. QN 102inc. CPU hours (h) 1.22 2.79 11.31 0.89 the phase field solid shell formulation. Firstly, a history variable ¯αthat governs the accumulation of fatigue is defined and the crack propagation is driven by a fatigue function f(¯α), which is employed to lower the fracture energy around the vicinity of crack. The governing equation for the phase field fracture driven by fatigue is given by Rd=∫B0 −2(1 −d)δdΨdΩ+∫B0 f(¯α)Gcl[1 l2dδd+ ∇Xd· ∇X(δd)]dΩ=0 (43) 19 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 Fig. 17. (a) Sketch of the single edge notched fatigue and (b) finite element mesh of specimen. Fig. 18. Cyclic fatigue loading with equal tension and compression for the single edge notched specimen. The history variable ¯αshould be independent of the unloading process and takes the form ¯α(t)=∫t 0 θ(α˙α)| ˙α|dτ(44) where θ(α˙α) is the Heaviside function and τis the pseudo time. Note that the fatigue history variable ¯αis a scalar variable depending on the loading history of material, which is defined as α=g(φ)Ψ. The function f(¯α) characterizing the relationship between the degradation of fracture energy and the number of cycles in the fatigue loading is expressed as f(¯α(t)) ={1 if ¯α(t)≤αT 2αT ¯α(t)+αTif ¯α(t)≥αT (45) where αTis the threshold value of the fatigue variable, below which the fracture energy will not be influenced by fatigue loading, and its expression is given by αT=GC 12l(46) 20 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 Fig. 19. The curves of crack length versus number of cycles obtained from both the quasi-Newton monolithic and staggered schemes for the phase field fracture modeling under fatigue loading. Table 4 Computation times for different schemes in the fatigue induced fracture modeling. (inc. per cycle) AM 16 AM 32 AM 80 AM 200 QN 16 CPU hours (h) 1.15 1.5 5.38 8.67 1.5 The sketch of the single edge notched fatigue is shown in Fig. 17(a). The geometry and dimensions of specimen are all the same as those of the previous single edge notched tension and shear except the shape of initial notch. To avoid contact between the upper and lower faces during compression, a V-shaped notch is adopted here, see Fig. 17(a). The finite element mesh with refinement in the potential crack region is shown in Fig. 17(b), and the total number of solid shell elements is 4878. The material and fracture properties are also the same as the previous study. The cyclic loading with equal tension and compression is applied to the upper surface of the specimen, and the amplitude and cycle are 0.002 mm and 0.002 s, respectively, see Fig. 18. The curves of crack length driven by fatigue loading versus number of cycles for quasi-Newton monolithic and staggered schemes with different time increment sizes are shown in Fig. 19. It is clear that extremely small time increment size is required for the staggered scheme so as to obtain the same result as that of the quasi-Newton monolithic solution, and convergence rate of the staggered scheme is quite low. The curves of total iterations versus number of cycles for the two solution schemes in the phase field fatigue simulation are shown in Fig. 20. Results show that the quasi-Newton solution requires almost the same number of total iterations as the staggered solution with 80 time increments per cycle. However, even the staggered scheme with 200 time increments per cycle cannot reproduce the quasi-Newton monolithic solution, and obviously more increments per cycle are required, which is computationally too expensive. The computation times for fatigue using both staggered and quasi-Newton monolithic schemes are listed in Table 4. Notably, the staggered solution with 200 increments per cycle requires 6 times more computation times than that in the quasi-Newton monolithic analysis. It can be concluded that quasi-Newton monolithic scheme shows the tremendous advantage over the staggered scheme in computational efficiency for fatigue problems. The contour plots of phase field values in fatigue induced fracture modeling with solid shell formulation solved by the quasi-Newton monolithic scheme are shown in Fig. 21. As can be clearly seen, the crack propagation starts from the 15th cycle of fatigue loading, and it proceeds with the increase of fatigue cycles. The single edge notched specimen for fatigue modeling is fully cracked at the 50th cycle. 21 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 Fig. 20. The curves of total iterations versus number of cycles obtained from both the quasi-Newton monolithic and staggered schemes for the phase field fracture modeling under fatigue loading. Fig. 21. The contour plots of phase field values in fatigue modeling using solid shell element solved by the quasi-Newton monolithic scheme. 6.5. Global–local fracture modeling of thin-walled structure with curved geometry In this section, the global–local phase field approach using solid shells is tested through the fracture modeling of cylindrical structure with an initial notch in the center as shown in Fig. 22(a). The dimensions of this specimen in the global level are the same as that in Section 6.3. As shown in Fig. 22(a), the global specimen is separated into 9 different blocks in total and the initial notch lies within the central block, which is chosen as the region of local model. To save computational cost, the global model is discretized with relatively coarse mesh using the uniform 22 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 Fig. 22. (a) Geometry and (b) finite element mesh of the global and local models in the phase field fracture submodeling of cylindrical structure. size of 3 mm, while to appropriately capture the crack propagation in the phase field modeling of local region, very refined mesh is used to discretize the potential crack area around the notch of local model. Specifically, the mesh size around the edge of local model is 3 mm, which is the same as that of global model, and characteristic mesh size along the potential crack path is 0.8 mm, see Fig. 22(b). The continuum solid shell is chosen for the elastic finite element analysis of global model, and phase field modeling for solid shell using quasi-Newton monolithic scheme is subsequently employed for the analysis of refined local model. The same material properties are adopted for the global and local models with the Young’s modulus equal to 210 GPa and the Poisson’s ratio equal to 0.3, and both the phase field length scale and fracture toughness of local model are the same as the previous study in Section 6.3, which are equal to 4 mm and 2.7 N/mm, respectively. The contour plots of phase field values of local model, displacement fields of both global and local models solved by quasi-Newton monolithic scheme at two different time points are shown in Fig. 23. The first time point is the time increment when the crack of local model starts to propagate under the drive of the boundary conditions interpolated from the displacement output of global model, while the second time point is the time increment when the local model with curved geometry is fully cracked. This global–local approach is very computationally efficient in fracture modeling of large thin-walled structures, since only the local region of interest needs to be modeled by phase field approach. In the global model, the degree of freedom corresponding to phase field variable is ignored, and thus significant computational cost can be saved. Besides, mesh in the local fracture region of interest can be refined with more flexibility, which is very appealing for realistic application of phase field approach in fracture modeling of large-scale structures. To further assess the global–local phase field approach using solid shells, the above case study is extended to the analysis of fatigue induced crack growth by applying the cyclic loading to the frontal end of global model. The cyclic loading is characterized by the amplitude of 0.05 mm and the cycle of 0.01 s, as shown in Fig. 24. The time duration is 1 s, and thus there are 100 cycles in total. The phase field plots of local model in the fatigue submodeling solved by quasi-Newton monolithic scheme are shown in Fig. 25. It should be pointed out that the boundary conditions of local model are obtained by the interpolation of the displacement output from the global finite element analysis, which is subjected to cyclic loading shown in Fig. 24. In this fatigue submodeling case, the crack of local model starts to propagate from the 40th cycle and after 28 cycles of fatigue loading, the local model is fully cracked along the expected crack path, which further demonstrates the capability of this global–local phase field approach for fracture modeling. 23 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 Fig. 23. The contour plot of global and local displacements as well as phase field variable in the fracture submodeling of cylindrical structure with solid shell element formulation. Fig. 24. Cyclic fatigue loading with equal tension and compression for the global model of the cylindrical structure. 7. Conclusions In this work, a global–local phase field approach using the enhanced assumed strain solid shell formulation with the efficient quasi-Newton monolithic solution is proposed for the numerical prediction of crack events in thin-walled structures at finite deformations. The EAS method for the treatment of volumetric and Poisson thickness locking was employed through the postulation of multi-field variational framework, and the ANS method to alleviate the transverse shear and trapezoidal locking effects was also taken into account in the shell formulation. With the aim of improving the computational efficiency without the sacrifice of robustness, an efficient quasi-Newton solution scheme is implemented within the phase field solid shell formulation in a monolithic manner, and its performance is demonstrated by comparison with the popular staggered Newton scheme through different paradigmatic boundary value problems, including single edge notched tension and shear, fracture of cylindrical structure under mixed bending and tension, and computationally demanding fatigue induced crack growth. Significant computational gains of the quasi-Newton monolithic scheme can be observed in all the numerical 24 Z. Liu, J. Reinoso and M. Paggi Computer Methods in Applied Mechanics and Engineering 399 (2022) 115410 Fig. 25. The contour plots of phase field values of local model in fatigue submodeling solved by quasi-Newton monolithic scheme. examples. It is also found that this scheme is very robust, capable of solving different benchmark problems of varying complexity without convergence issues frequently occurring in the conventional monolithic Newton solution. Besides, in case that crack growth in the local region of interest has minor effect on the global stiffness, a specific global–local approach in the 3D setting tailored for the phase field modeling with solid shells is proposed to save computational cost. Coarse mesh can be adopted for the global model while the mesh of local model can be more refined to deal with the appropriate displacement-phase field problem. The capability of this approach was demonstrated by simulating the crack growth of cylindrical thin-walled structure under both static and fatigue cyclic loadings, which could open up new possibilities in tackling realistic industrial problems concerning fracture events in large-scale thin-walled structures. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability Data will be made available on request. Acknowledgments The authors acknowledge the funding received from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 861061 – Project NEWFRAC - New strategies for multifield fracture problems across scales in heterogeneous systems for energy, health and transport. JR is grateful to the funding of Consejer´ ıa de Econom´ ıa y Conocimiento of the Junta de Andaluc´ ıa (Spain) under the grant P20-0595 and Spanish Ministerio de Ciencia, Innovaci´ on y Universidades the under the grant PID2019-109723GB-I00. 25