A combined XFEM phase-field computational model for crack growth without remeshing
Alba Muixí, Onofre Marco, Antonio Rodríguez-Ferran, Sonia Fernández-Méndez
AA combined XFEM phase-field computational model for crack growthwithout remeshing
Alba Muix´ı, Onofre Marco, Antonio Rodr´ıguez-Ferran, Sonia Fern´andez-M´endez ∗ Laboratori de C`alcul Num`eric (LaC`aN)Universitat Polit`ecnica de CatalunyaBarcelonaTech, Barcelona, Spain
Abstract
This paper presents an adaptive strategy for phase-field simulations with transition to fracture. Thephase-field equations are solved only in small subdomains around crack tips to determine propagation,while an XFEM discretization is used in the rest of the domain to represent sharp cracks, enabling touse a coarser discretization and therefore reducing the computational cost. Crack-tip subdomains moveas cracks propagate in a fully automatic process. The same computational mesh is used during all thesimulation, with an h -refined approximation in the elements in the crack-tip subdomains. Continuityof the displacement between the refined subdomains and the XFEM region is imposed in weak formvia Nitsche’s method. The robustness of the strategy is shown for some numerical examples in 2D and3D, including branching and coalescence tests. Keywords:
Phase-field modeling · Brittle fracture · Crack propagation · Continuous-discontinuousmodels · Adaptive refinement · Nitsche’s method · XFEM
Typical approaches to model brittle or quasi-brittle fracture can be classified in discontinuous and con-tinuous models, depending on the way cracks are described.In discontinuous models , cracks are represented as discontinuities in the displacement field ( sharpcracks ), embedded in an elastic bulk. Specific mechanical criteria – typically based on linear elasticfracture mechanics (LEFM) concepts and based on crack-tip information – are needed to model crackinception, propagation and branching. A source of difficulty is the singularity of the strain and stressfields in classical linear elasticity at crack tips. For this reason, stress intensity factors (SIF) are oftenused [19].From a computational viewpoint, two different strategies may be used to describe the displacementjump associated to a sharp discontinuity: node and edge duplication, in the spirit of the original cohesivezone model (CZM), see [24], or the extended finite element method (XFEM) and related techniques,where the crack path is not constrained to follow element edges and therefore remeshing is avoided [2, 20].Together with a certain mesh-dependence of the crack path (for CZM) or geometrical complexity (forXFEM), discontinuous models offer various advantages, namely: i) the explicit representation of thecrack facilitates the modeling of crack-surface physics; ii) the sharp discontinuity avoids any spuriousinteraction between the two faces of the crack; iii) coarse discretizations may be used in the wake of thecrack tip. ∗ Corresponding author. E-mail: [email protected] a r X i v : . [ c s . C E ] J u l ontinuous models , usually phase-field or gradient-damage models, assume continuous displacementfields and represent cracks as damaged regions that have lost their load-carrying capacity ( diffuse cracks ).Information on the evolution of cracks is incorporated into the field equations; initiation, propagation,branching and coalescence of cracks are automatically handled by the model.The main drawback of these models is their high computational cost. Very fine meshes are neededto capture narrow damage bands that resemble sharp cracks. For cases in which the crack path is notknown a priori, one may use an adaptive strategy to automatically refine the discretization [23, 22, 21, 8].Continuous models do not explicitly identify a crack surface within damage bands and spurious interactionremains between physically separated parts. Several works adressing this issue have been presented whichcombine continuous and discontinuous approaches [29, 7].In continuous-discontinuous models both sharp and diffuse representations of cracks are overlapped.The discontinuous model is used to represent material separation while the continuous one is used todetermine the evolution of cracks. Critical issues are related to the introduction of discontinuities: whenand where (switching and locating criteria, respectively). Still, these models have a high computationalcost from solving the continuous model in the whole domain.Recent approches consider a phase-field problem only in moving subdomains containing crack tips toimprove the efficiency of the simulations. Giovanardi et al [10] propose a two-scale strategy. A globalsolution for the displacement field is obtained with XFEM. Then, the propagation is determined by solvinga phase-field problem in subdomains containing the crack tips with a finer mesh. The strategy is appliedto propagation examples in 2D, with no branching nor coalescence.Patil et al [25] present a similar approach: phase-field is used only in small crack-tip regions, andcracks are described by an XFEM discretization in the whole domain. In the so-called crack-tip regions,the mesh is adaptively refined with a uniform submesh. The multiscale FEM is used to couple thefiner discretization with the coarser one, by means of transition elements and the construction of multi-scale shape functions relating the degrees of freedom from the two levels of refinement. The strategy issuccessfully applied to branching and coalescence of cracks in 2D.With the same philosophy as the above-mentioned proposals, we present a model combining phase-field in small subdomains around crack tips, and a discontinuous model in the rest of the domain. Here,the two descriptions of cracks are not overlapped. The phase-field model is used to handle the propagationand coalescence of cracks, while the discontinuous model explicitly describes cracks and enables to use acoarser discretization in almost the whole domain. The presented methodology is stated in 2D and 3D.As a simplification, we use the term crack tips to refer both to crack tips, in 2D, and crack fronts, in 3D.The small phase-field subdomains consist of a set of elements for which the discretization is uniformly h -refined to capture the solution, following the refinement strategy proposed in [21], and move with cracktips as cracks propagate in a fully automatic process. The mesh is fixed during all the simulation. As thediffuse cracks evolve, refined elements are nested near crack tips. At the same time, refined elements inthe wake of the crack transition into an XFEM coarser approximation.On the interface between continuous and discontinuous subdomains, continuity of the displacementfield is imposed in weak form by means of Nitsche’s method. Imposing continuity in weak form avoidsthe definition of transition elements and having to constrain the degrees of freedom from the richer part.This leads to a very local refinement, which is handled by the method for any refinement factor.The paper is structured as follows. First, in Section 2, we recall the discontinuous and the phase-fieldmodels. Section 3 describes the idea behind the combined computational model. In Section 4 we givethe FE formulation, together with some implementation details. Sections 5, 6 and 7 are devoted to theevolution of the crack-tip subdomains and the transition to fracture as cracks propagate. In Section 8,we present some representative numerical examples in 2D and a fully 3D test to show the capabilities ofthe strategy. The conclusions in Section 9 close the paper.2 XFEM and phase-field computational models
In this paper, we combine a discontinuous model and a phase-field model to simulate crack propagationin the quasi-static regime. In the discontinuous model, cracks are described by discontinuous displace-ments fields. In phase-field, cracks are approximated by diffuse bands with a loss of stiffness, and thedisplacement field is continuous across them. In this section we give a brief overview of both approaches.As a representative phase-field model, we use here the hybrid model proposed by Ambati et al [1].Throughout the paper, we consider a linear elastic body with a traction-free crack, occupying a domainΩ ⊂ R n sd , with n sd = 2 ,
3, and under the assumption of small deformations, as illustrated in Figure 1.We restrict to linear isotropic materials.
If no body forces are applied and inertia effects are neglected, the balance of linear momentum for thebody in Ω leads to the system ∇ · σ = in Ω C , u = u D on Γ D , σ · n = t N on Γ N , σ · n = on Γ + C ∪ Γ − C , (1a)(1b)(1c)(1d)where Ω C = Ω \ Γ C is the cracked domain. The relation between the displacement field u and the stresstensor σ is given by the constitutive equation σ = ∂ Ψ ( ε ) ∂ ε , (2)where ε is the small strain tensor and Ψ is the elastic energy density. For linear isotropic materials,Ψ = ( ε : C : ε ) / σ = C : ε , with C a fourth order positive definite tensor depending on the Lam´eparameters.The Dirichlet and Neumann boundaries are denoted by Γ D and Γ N , respectively, and satisfy ∂ Ω =Γ D ∪ Γ N and Γ D ∩ Γ N = ∅ . Prescribed displacements and tractions are u D and t N , and n is the exteriorunit normal.Figure 1: Body with a sharp crack Γ C (left) and a diffuse phase-field representation of the crack withdamage variable d (right). The profile of the displacement field along the dotted section shows how thediscontinuity from the sharp crack is approximated by a continuous function with a steep variation inphase-field models. 3quation (1d) imposes traction-free conditions on the crack faces, Γ + C and Γ − C .Numerical techniques to solve the system in (1) need to account for the discontinuity of displacementsacross the crack. The most popular method used to this intent is XFEM [2, 20], which enriches thediscretization space to properly represent discontinuous displacements with no need to adapt the meshto the crack geometry. XFEM has demonstrated its applicability and computational efficiency in manyapplications [6], but it presents some limitations for crack simulations.It can be analytically shown that stresses around crack tips vary according to r − / , where r is thedistance to the crack tip [28]. This stress singularity is physically unrealistic and, from our experience,it leads to difficulties in the definition of a proper XFEM enrichment around crack tips. More precisely,to ensure convergence under uniform mesh refinement, the area were the crack-tip enrichment is appliedshould be kept fixed; but this strategy leads to severe ill-conditioning. On the other hand, if the area forthe crack-tip enrichment is reduced with the element size, asymptotic convergence cannot be observed[15]. Thus, from a practical point of view, crack-tip enrichment is not straightforward. An overview ofthe advances of the method to solve problems in fracture mechanics can be found in Sukumar et al [27].However, the main limitation of the discontinuous approach is that crack inception and propagationare not described by the model. Different criteria have been proposed to determine the direction andvelocity of propagation using local information of the solution at crack tips, see for instance [5, 26, 4],but a robust theory is not well-established yet. This is not the case for phase-field computational models,which handle crack evolution in a natural way. In phase-field models, cracks are represented by means of the so-called phase-field variable or damagefield , which is denoted by d . The damage takes value 1 on the fracture path, value 0 in intact parts ofthe material and smoothly varies between the two values, as sketched in Figure 1. The damage field isintroduced as an additional unknown into the system of equations. Thus, evolution of cracks is determinedby the model itself.The hybrid phase-field model by Ambati et al [1] is characterized by incorporating a tension-compressionsplitting, while maintaining a linear equilibrium equation within a staggered scheme to solve the system.The governing equations are ∇ · σ = with σ = g ( d ) ∂ Ψ ( ε ) ∂ ε , − l ∆ d + d = 2 lG C (1 − d ) H + ,g ( d ) := (cid:40) (1 − d ) where Ψ +0 ≥ Ψ − , G C is the energy release rate of the material. The length parameter l is related to the width of thediffuse cracks and it is usually taken very small to mimic sharp cracks. This entails the need of very finemeshes around cracks to properly approximate the solution.The stress tensor in the equilibrium equation (3a) is degraded by the function g ( d ). In fully brokenparts of the material, where d = 1, this implies a complete loss of stiffness. In parts of the material with d = 0, the linear elastic stress-strain constitutive relation (2) is recovered.The model considers an additive decomposition of Ψ , namely, Ψ = Ψ +0 + Ψ − with Ψ +0 and Ψ − thetensile and compressive parts of Ψ , respectively. Here, we assume the tension-compression splitting byMiehe et al [17, 18] which is based on the spectral decomposition of ε . Denoting the principal strains by { ε i } i =1 ,...,n sd and the principal strain directions by { d i } i =1 ,...,n sd ,Ψ ± ( ε ) = 12 λ (cid:104) tr( ε ) (cid:105) ± + µ tr (cid:0) ε ± (cid:1) , (4)4here ε ± = (cid:80) n sd i =1 (cid:104) ε i (cid:105) ± d i ⊗ d i and (cid:104)(cid:12)(cid:105) ± = ( (cid:12) ± | (cid:12) | ) / d is modeled by equation (3b). The history field variable, H + , isdefined as H + ( x , t ) = max τ ∈ [0 ,t ] Ψ +0 (cid:0) ε ( x , τ ) (cid:1) , (5)and therefore the propagation of d is only caused by the tensile Ψ +0 . The definition of H + also enforcesirreversibility of fracture [17, 18].The condition in equation (3c) recovers the original stiffness of the material under compression, inorder to prevent interpenetration of crack faces. Alternatively, the splitting of Ψ may be included intothe stress-strain equation in (3a), but then the equilibrium equation becomes nonlinear.The system in (3) is solved with boundary conditions σ · n = t N on Γ N , u = u D on Γ D , ∇ d · n = 0 on ∂ Ω . (6)It is important noting that the value of l determines the element size required around cracks. Thus,differently to other applications, adaptive refinement can be directly driven by the crack growth, witha refinement element size that is known a priori . This is the idea exploited in [21], where elementsalong the crack are automatically refined considering a refined reference element, and keeping the originalbackground finite element mesh in the whole computation. The weaknesses of the adaptive phase-fieldmodel are the inability of phase-field to explicitly describe crack opening, which is necessary in someapplications, and, more importantly, the computational cost associated to the refinement in all elementsalong the crack path. Thus, the next natural step is derefining the elements in the wake of the crack tipsand replacing the diffuse band by a sharp crack, where no variation of the damage is expected, by meansof an XFEM enrichment.In the next section, XFEM and the adaptive phase-field model in [21] are combined in a new methodthat has the computational efficiency of XFEM, and also the ability of phase-field models to handle crackevolution. A novel computational model combining a phase-field model and an XFEM discontinuous approach isproposed. The phase-field equations are solved in small subdomains around crack tips, Ω tips , and thediscontinuous model is used in the rest of the domain, Ω xfem . The corresponding diffuse and sharpdefinitions of cracks are not overlapped and we denote by Γ the interface between subdomains, this is,Ω = Ω tips ∪ Ω xfem , Ω tips ∩ Ω xfem = ∅ , Γ = Ω tips ∩ Ω xfem . The model overcomes the limitations of both approaches. The phase-field model drives the crackpropagation and deals with branching and merging of cracks, whereas the crack is explicitly described inalmost all the domain by a sharp discontinuity. Thus, the high spatial resolution needed around crackpaths in standard phase-field models is now only necessary in small neighborhoods of the crack tips, withthe consequent important saving in computational cost.The partition of the domain is based on the computational mesh. Crack-tip subdomains Ω tips aredefined as the set of elements close to crack tips and tips of notches in the domain, where crack growth isexpected. The rest of elements are Ω xfem . The initial definition of Ω tips is crucial to detect crack inception,since the damage field is solved only in this part of the domain.As crack tips advance during the computation, the domains Ω tips and Ω xfem are repeatedly updated,refining elements in the nose of the cracks and derefining elements in their wake. When derefining, diffusecracks are replaced by XFEM sharp traction-free cracks.5 xfem ⌦ tips Figure 2: Scheme of a discretization in two consecutive load steps. In Ω tips , elements are uniformly h -refined and cracks are represented by diffuse bands. The FE discretization in Ω xfem is the one corre-sponding to the initial mesh, with an XFEM Heaviside enrichment to introduce sharp cracks, representedwith black line. The interface Γ is depicted in red.Following the concept in [21], elements in Ω tips are automatically identified and uniformly refined toproperly approximate the phase-field solution, now only around crack tips. The refinement is based on arefined reference element, keeping the original background mesh during the whole computation.Figure 2 illustrates the discretizations for two consecutive load steps. The elements in Ω tips are h -refined with a uniform submesh with m n sd subelements. The refinement factor m is chosen a prioridepending on the length-scale parameter l of the phase-field model. As the crack evolves, the smallsubdomains move with the crack tips and the diffuse band becomes a sharp discontinuity in the wake ofthe crack.The main challenge is gluing the two models on the interface Γ, with different discretizations andrepresentations of cracks on each side.For the equilibrium equation, continuity of the displacement is imposed in weak form by means ofNitsche’s method. The Nitsche’s formulation, recalled in Section 4, is the one considered in [21], but nowwith an XFEM enriched approximation in Ω xfem . In addition, when imposing continuity, a small portionof the interface around the intersections with the crack is discarded.Numerical experiments show that unrealistic displacement fields are obtained if the continuity isenforced on the whole interface Γ. The reason is that the description of the crack on both sides of theinterface is different. From Ω xfem the crack is sharp and the material retains all the stiffness, whereas onΩ tips we have a smooth representation of the crack with a material degradation given by g ( d ), leading toalmost fully damaged material close to the crack path. Thus, for parts of Γ which are crossed by cracks,the difference of stiffness between the two subdomains is quite significant.Continuity is then enforced on a slightly cropped interface defined as (cid:98) Γ = { x ∈ Γ : d ( x ) < . } . (7)Figure 3 illustrates the concept. A simple test comparing the results obtained imposing continuity on Γand on the cropped interface, (cid:98) Γ, is presented in Section 8.1.It is worth noting that imposing continuity in weak form, by means of Nitsche’s method, handles thenon-conformal approximations in a natural way. That is, refined discretizations are directly attached tonon-refined discretizations, with a standard finite element approximation or with XFEM enrichment. Notransition elements are required, avoiding the spreading of the refinement or the enrichment, and withouthaving to deal with hanging nodes.The other important issue for gluing the computational models is setting proper boundary conditionsfor the damage variable d , since we solve for it only in Ω tips . In standard phase-field models, the boundary6igure 3: Continuity of displacements is imposed in weak form on the interface (cid:98) Γ, in blue (left). Cracksdo not necessarily intersect Γ perpendicularly. Dirichlet boundary conditions for the damage field d areimposed on Γ dD , i.e. on the side intersected by the crack and on the adjacent sides, in green (right).conditions for the damage equation are usually taken as homogeneous Neumann on ∂ Ω. However, wecannot assume that cracks intersect the interface Γ with perpendicularity. With this in mind, Dirichletboundary conditions are imposed on the faces that have been crossed by a crack. More precisely, theDirichlet boundary, Γ dD , consists of the element sides (in 2D, faces in 3D) on Γ that are cut by the crack,i.e. with d > .
95 at some node, together with their neighboring faces on Γ. See an example in Figure3 (right). When transitioning elements intersected by the crack from Ω tips to Ω xfem , the damage field onΓ dD is saved to be used as Dirichlet value for the damage equation. That is, the damage is assumed tobe frozen in the wake of crack tips. On the rest of the boundary, ∂ Ω tips \ Γ dD , homogeneous Neumannconditions are imposed not to favor any direction of crack growth. In this section, we present the formulation, commenting also on some implementation details. The fracturemodel is solved in a staggered scheme. Thus, the numerical formulations for the equilibrium and for thedamage equation are independent.In the equilibrium equation, the discretization has to account for discontinuous displacements acrosscracks in Ω xfem . Also, continuity of displacements between Ω tips and Ω xfem has to be handled by theformulation. On the other hand, the damage equation is solved with the standard finite element methodwith refined approximation in Ω tips . The equilibrium equation is solved in Ω, with different non-conformal approximations for the displacementfield in the partition given by Ω tips and Ω xfem . Thus, the displacement is assumed to be in the space[ V (Ω)] n sd with V (Ω) = (cid:8) v ∈ L (Ω) : v | Ω tips ∈ H (Ω tips ) , v | Ω xfem ∈ H (Ω xfem \ Γ c ) (cid:9) , which includes functions that are discontinuous across the interface Γ, and across the sharp crack Γ c ; thesecond one, handled by the XFEM enrichment.Our choice to weakly impose continuity on the interface Γ is Nitsche’s method. The formulation isbased on the Nitsche’s formulation for linear elasticity applied to interface problems by Hansbo [13], seealso [21].In Nitsche’s method, the weak form of classical finite element approach is modified in order to imposecontinuity of the solution on the interface, while keeping a symmetric and coercive bilinear form, andalso imposing equilibrium of tractions. Differently to Lagrange multipliers, no additional unknowns areintroduced to the system. 7he weak form reads: find u ∈ [ V (Ω)] n sd such that u = u D on Γ D and (cid:90) Ω ∇ v : σ ( u ) d V − (cid:90) (cid:98) Γ (cid:74) v ⊗ n (cid:75) : { σ ( u ) } d s − (cid:90) (cid:98) Γ { σ ( v ) } : (cid:74) u ⊗ n (cid:75) d s + β E (cid:90) (cid:98) Γ (cid:74) u ⊗ n (cid:75) : (cid:74) v ⊗ n (cid:75) d s − (cid:90) Γ N v · t N d s = 0 , (8)for all v ∈ [ V (Ω)] n sd such that v = on Γ D , where β E is the Nitsche’s stabilization parameter, and (cid:98) Γis almost the whole interface Γ, discarding small portions around the intersections with the cracks, asdefined in (7) and illustrated in Figure 3 (left).The mean and jump operators are defined as {(cid:12)} = ( (cid:12) t + (cid:12) x ) and (cid:74) (cid:12) n (cid:75) = (cid:12) t n t + (cid:12) x n x =( (cid:12) t − (cid:12) x ) n t , with the lower indices, t and x , indicating the values on (cid:98) Γ from Ω tips and Ω xfem , respectively,and n t , n x the corresponding unit exterior normals. A derivation of this weak form can be found in [21].The differences here are just the aproximation spaces to be glued, now including sharp cracks via XFEM,and the small cropping of the interface, (cid:98) Γ.Even though it has been omitted for simplicity, note that the stress tensor σ in (8) depends also onthe damage field d in Ω tips .The stability of the formulation depends on the value of β E . This parameter has to be large enoughto ensure coercivity of the bilinear form, thus leading to a discrete system with a positive-definite matrix.To obtain optimal orders of convergence, this parameter can be taken of the form β E = α E E ( h/m ) − , (9)with E the Young’s modulus, h the element size in the background mesh, m the refinement factor forrefined elements, and α E a large enough positive constant independent of other numerical and materialparameters. A lower bound for its value can be found by solving an eigenvalue problem [11]. However,the formulation is very robust on the parameter and the value of α E can be easily tuned by numericalexperimentation [21]. The computation is based on a fixed background mesh during the whole simulation, with elements { K i } n el i =1 .That is, Ω = n el (cid:91) i =1 K i , K i ∩ K j = ∅ if i (cid:54) = j. As cracks evolve, in every load step and iteration of the staggered scheme, the elements in the backgroundmesh are identified as refined, i.e. in Ω tips , if they are close to crack tips or tips of initial notches, or inΩ xfem , otherwise.In Ω xfem , the XFEM approximation is based on a standard finite element space, that is V h (Ω xfem ) = { v ∈ H (Ω xfem ) : v | K i ∈ P p for K i ⊆ Ω xfem } , with nodal basis { N i } i ∈ I xfem , where I xfem is the set of the indices of the nodes in Ω xfem and P p is the spaceof polynomials up to degree p . The XFEM approximation space in Ω xfem in the case of a unique crack isthen V h xfem = < N i > i ∈ I x ⊗ < HN i > i ∈ I enr , where I enr is the set of indices of enriched nodes, i.e. the nodes of the elements which are cut by thecrack, and H is the Heaviside function taking, for instance, values 1 and − V h xfem = < N i > i ∈ I x ⊗ < H N i > i ∈ I ⊗ . . . ⊗ < H n cracks N i > i ∈ I n cracksenr , p = 2 and refinement factor m = 2. Nodes arerepresented by gray dots and integration points, by black crosses. The refined element has 5 basis functionsand 6 integration points.where n cracks is the total number of crack branches [20]. Here, no asymptotic tip enrichment is necessarysince Ω xfem does not contain crack tips.A proper numerical quadrature must be defined in the elements intersected by the crack to accountfor the discontinuity of the XFEM approximation. That is, a numerical quadrature must be defined ineach portion of cut elements, see for instance [16, 12].On other hand, in Ω tips , a uniformly h -refined discretization is considered to capture the sharp variationof the phase-field solution. These elements are mapped into a uniformly refined reference element. Theresulting approximation space is equivalent to a standard finite element approximation on a finer mesh,non-conformal on the interface Γ. The chosen implementation, based on a refined reference element, issimpler since no mesh generation is needed, the assembly and the integration can be done as usual, andthe definition of the refined reference element is done as a simple preprocess. This is a viable optionin this application since all refined elements have the same refinement factor, which is known from thebeginning. The refinement factor is computed from the phase-field length, l , the characteristic elementsize of the background mesh, h , and the degree of approximation p . It can be taken as m = (cid:24) ahlp (cid:25) . In our experience, a reasonable value for the factor a is between 3 and 5. The approximation space inΩ tips is then defined as V h tips = { v ∈ H (Ω tips ) : v | K i ∈ P p ref ( K i ) for K i ⊆ Ω tips } , where P p ref ( K i ) is the refined space in each element. The element is split into m n sd uniform subelements,and the refined approximation space in the element is P p ref ( K i ) = { v ∈ H ( K i ) : v | K ij ∈ P p ( K ij ) , j = 1 ...m n sd } , with { K ij } m nsd j =1 denoting the subelements in the element K i . To illustrate the idea, a 1D refined elementis shown in Figure 4.Finally, the approximation space to discretize the Nitsche’s weak form (8) is the combined space V h (Ω) = { v ∈ H (Ω) : v | Ω xfem ∈ V h xfem , v | Ω tips ∈ V h tips } . The damage equation is solved only in Ω tips , this is, in the small subdomains containing crack tips andtips of notches, where damage evolution is expected. As commented in Section 3, and illustrated in Figure3 (right), Dirichlet boundary conditions are imposed on Γ dD , defined as a portion of the boundary around9he intersection with cracks. Homogeneous Neumann boundary conditions are imposed on the rest of theboundary.The weak form is then: find d ∈ H (Ω tips ) such that d = d D on Γ dD and (cid:90) Ω tips (cid:18) G C l + 2 H + (cid:19) vd d V + (cid:90) Ω tips G C l ∇ v · ∇ d d V = (cid:90) Ω tips v H + d V, (10)for all v ∈ H (Ω tips ) such that v = 0 on Γ dD .The value set on the Dirichlet boundary, d D , is the value saved from the damage solution in previousiterations, assuming that the damage is frozen on the wake of the crack.The weak form is discretized with the refined approximation space V h tips . Numerical simulations take an incremental process in load steps. At every load step, the system is solvedwith a staggered scheme, for which the equations are solved alternately until convergence. For load step n , we iterate over i. Computing the displacement field u in Ω by solving the weak form (8), with the current damage field d in Ω tips , and with boundary data u D = u nD on Γ D and t N = t nN on Γ N . Recall that traction-freeconditions on the crack are imposed by means of the XFEM enrichment in Ω xfem , and a refineddiscretization is considered in Ω tips . ii. Updating the history field H + in Ω tips . iii. Computing the damage field d in Ω tips by solving the weak form (10). iv. Updating the partition , Ω tips and Ω xfem , and the geometrical description of sharp cracks, as detailedin Sections 4.5, 5.1 and 6.As a convergence criterion we check if the error of the damage field in the Euclidean norm is lowerthan a certain tolerance.Note that crack-tip subdomains are updated at every staggered iteration. Since we are modelingbrittle fracture, cracks can significantly grow at a single load step and the phase-field subdomain has tobe accordingly modified to allow the propagation. The geometry of the domain is defined by the background mesh (
X, T ), with X the nodal coordinatesand T the connectivity matrix. This mesh is kept fixed during all the computation.In the preprocess, information of the faces is computed from the connectivity matrix. Faces arenumbered and, for each one of them, we store the number of the elements sharing the face and the localnumber of the face in each element. In 3D, the rotation is also stored, that is, which node in the face fromthe second element corresponds to the first node of the face from the first element. This information onfaces is used to compute the integrals over ˆΓ in the equilibrium weak form (8).Information of the partition { Ω tips , Ω xfem } is given by three vectors storing, respectively, the numbersof the elements in each subdomain and the numbers of the faces composing the interface Γ. Also, weneed to save the list of nodes on Γ dD and the value of the damage field on these nodes from the previousiteration, to be able to impose the boundary conditions in (10). During the computation, information onthe partition is updated at every staggered iteration according to the propagation of cracks by applyingthe criteria in Section 5.1.A geometrical description of sharp cracks in Ω xfem is needed to introduce the discontinuities into theXFEM discretization. Here, we store the crack path together with a nodal vector indicating whether each10nriched node is in the left or right side of the crack. Some details on the construction of the sharp crackcan be found in Section 6.In addition, to ensure continuity of the approximation in the refined region Ω tips , we define an auxiliarymesh ( X ref , T ref ). This new mesh is constructed by mappings of the refined reference element into theelements in Ω tips , and it is modified as the list of elements in Ω tips is updated. When an element is addedto Ω tips , the nodes in the refined element are added to X ref , avoiding repetitions, and the correspondingrow is added to T ref . When an element is supressed, nodes not belonging to other refined elements areremoved from X ref . T ref is modified according to the new nodal numbering and the row corresponding tothe element is removed. The use of a refined mesh is optional, since the computations are directly donewith the refined reference element, but it clearly facilitates the assembly for the elements in Ω tips . The criteria to decide if an element should be in Ω tips or in Ω xfem are based on i) the maximum valueof the damage reached in the element, and ii) the distance to elements containing crack tips (or tips ofnotches). These criteria are applied at every iteration of the staggered algorithm.In this section, we describe the criteria and also the algorithm to identify the elements that containcrack tips in 2D. The conversions from a phase-field to a sharp representation of cracks and vice versa arediscussed in Sections 6 and 7, respectively. The refining criterion is applied to elements in Ω xfem which are adjacent to the interface Γ, in order todecide if they should be transferred to Ω tips . Following the idea in [21], an element is refined if the damagefield reaches a threshold value d ∗ at one of its nodes. Here, since the damage is computed only in Ω tips ,an element in Ω xfem is refined if any of its nodes on the interface reaches the threshold value d ∗ and it isclose to some element containing a crack tip, with a threshold distance δ ∗ . This last consideration ensuresthat elements in the wake of cracks stay in Ω xfem .Then, to decide if some of the elements in Ω tips have to switch to the discontinuous approximation,the distance criterion with threshold value δ ∗ is also applied. More precisely, if the minimum distanceof an element to all the elements containing a crack tip (or notch tip) is larger than δ ∗ , the element istransferred to Ω xfem .Distances between elements are computed center-to-center. In practice, values of d ∗ between 0 . . δ ∗ is around 10 l . Here, we focus on quadrilateral elements in 2D. The logic of the algorithm is easily extendable to triangularelements. As a simplification, we assume a smeared crack can only intersect the same element edge once.Elements containing crack tips are identified using the number of sides which are intersected by thesmeared band and the area of the band inside the element. We assume a side is intersected if it reaches d > .
95 at some of its integration points. The area of the band is approximated as the sum of the weightsof integration points in the element with d > . tips , we count how many sides are intersected by the crack. The consideredconfigurations are summarized in Figure 5.In elements with only 1 intersected side, we mark the element as containing a crack tip if the area ofthe band is larger than a threshold A ∗ ; on the other case, the crack is assumed to be tangent to the side.For 2 intersected sides, the element contains a crack tip if the two intersected sides are adjacent, theshared node satisfies d > .
95 and the area of the band is larger than A ∗ . All other cases are discarded.11igure 5: Detail of a refined element with possible crack paths. Cases are classified depending on thenumber of intersected sides. Red crosses indicate integration points on the sides with d > . A ∗ = hl/ In elements in the wake of crack tips, no variation of crack paths is expected and the phase-field damageband is replaced by a sharp representation. This implies a significant reduction in the computational costof simulations. The transition is done for fully degraded material, so no additional energetic considerationsare needed.Sharp cracks are defined by the union of elemental contributions. Once a cut element transitions fromΩ tips to Ω xfem , the crack path is identified within the phase-field diffuse band and is then added to theexisting sharp crack. The process is illustrated in Figure 6.For both 2D and 3D problems, we search for intersections of the diffuse cracks along the edges of themesh of transitioning elements. In the 1-dimensional searchs, intersections are taken as the middle pointsof the nodes (in the refined discretization) with d > .
98. In 2D, this leads to piecewise linear cracks,with a segment in each cut element. In the event of crack branching, a proper piecewise representation is12onsidered in the element containing the branching point, as shown in Section 8.3. In 3D, crack surfacesare constructed as the tringular facets defined by the intersection points on edges. Note that these facetscan have very different sizes. This is not a problem since the surface is used only to define the integrationsubdomains in the element.This algorithm is enough to show the capabilities of the strategy. In our numerical examples, it leadsto similar results to the ones obtained by a plain phase-field model. More sophisticated techniques canbe introduced at this step if more accuracy is needed or crack paths are more challenging, such as ahigher-degree curved representation of the crack path in each cut element, given by more than 2 pointsin 2D [12], the medial-axis algorithm [29] or the optimization-based approach in [7].It is important for the different representations along a crack (sharp and diffuse) to match on theinterface Γ with enough accuracy, to avoid the creation of unphysical corners. In our experience, level-set representations of the crack in the coarse background mesh are not accurate enough to fulfill thisrequirement.Regarding the displacement field, since it is computed in the staggered scheme in Section 4.4, there isno need to compute its projection into the new approximation space.
When an element in Ω xfem is selected to be refined, nodal values for the damage have to be defined. If theelement is refined for the first time, nodal values of the damage can be set to zero except for the nodeson the interface, where the damage is known. If the element is intersected by a sharp crack, the damagefield in the element has to approximate the crack. This is the case in cracks coalescence.In experiments with merging of cracks, we may have elements cut by a sharp crack which are ap-proached by a phase-field crack tip. In this case, the criteria in Section 5.1 lead to refining again theelements, transitioning back to the continuous representation. The coalescence is then handled by thephase-field model. Once the cracks have merged, the elements transition to the discontinuous again withthe corresponding update in the involved sharp cracks.In the current implementation we store the damage field for elements that transition to discontinuousand recover its value in case these elements need to be refined again. Another option would be to constructthe damage band from scratch using the sharp representation, for instance, defining the correspondinghistory field variable, H + , following Borden et al [3].Figure 6: Sketch of transition to sharp crack. The diffuse crack is replaced by a sharp crack in the elementselected to be derefined, marked with red square. In 2D the crack is composed by linear segments, definedby the intersection of the diffuse crack with the elements sides.13 Numerical experiments
In this section, we aim to demonstrate the robustness of the proposed approach to efficiently simulatefracture processes. During the computations the background mesh is kept fixed and the discretization isautomatically modified to account for the different representations of cracks.First, a simple test is presented to numerically validate the definition of ˆΓ to impose continuityof displacements between the subdomains Ω tips and Ω xfem . Then, several benchmark tests are visited,covering a wide range of cases in fracture simulations. Examples in 2D prove the applicability of thestrategy to branching and coalescence of cracks. The approach is also applied to a fully 3D setting.Plane strain conditions are assumed in 2D. In all examples, we solve for degree of approximation p = 1. The tolerance for convergence in the staggered scheme is fixed to 10 − . The Nitsche’s parameterfor the equilibrium equation is α E = 100 and the refinement of elements is triggered by the thresholdvalue d ∗ = 0 .
2. A numerical study on the influence of α E and d ∗ in adaptive phase-field simulations canbe found in [21].Preexisting cracks in the domain are first introduced as diffuse phase-field bands by solving the damageequation with an initial history field, H +0 , as described in Borden et al [3]. Then, before initializing thesimulation, the procedure of transition to discontinuous is applied to replace the diffuse phase-field bandby a sharp crack where needed, according to the respective switching criterion. In this test, we analyse the effect of imposing continuity of displacements on the whole interface Γ and onthe cropped interface ˆΓ (defined in (7) as the part of Γ where the material is not significantly degraded)on a simple tension test with a fixed crack.Consider a square plate in [ − . , . × [ − . , .
5] mm with a horizontal crack at midheight crossingthe whole piece. The piece is clamped on its bottom face and has imposed displacements (0 , u D ) on itstop face, with u D = 10 − mm. The material parameters are E = 20 GPa, ν = 0 . G C = 10 − kN/mm.The crack is represented by a sharp crack for x < x >
0, witha length-scale parameter l = 0 .
012 mm. The domain is covered by a uniform quadrilateral mesh with12 ×
15 elements and Ω tips is taken as { x ≥ , | y | ≤ .
075 mm } . The rest of the domain is Ω xfem . Elementsare refined with a uniform submesh of refinement factor m = 15. The discretization is shown in Figure 7.Figure 8 shows the vertical displacement field u y when the equilibrium equation is solved imposingcontinuity of displacements on ˆΓ. As expected, the upper half moves rigidly, with a constant verticaldisplacement equal to u D , whereas the lower half does not move and has zero vertical displacement. InFigure 7: Test on continuity.
Initial discretization and approximation of the crack.14ontinuity on ˆΓ Continuity on ΓFigure 8:
Test on continuity.
Vertical displacement field u y imposing continuity on ˆΓ and on Γ.Figure 9: L-shaped test.
Geometry and boundary conditions. Dimensions in mm.Ω xfem we obtain the expected discontinuity and in Ω tips the displacement is continuous and abruptly variesbetween these two values.We repeat the experiment imposing continuity on the whole interface Γ. The vertical displacementsexhibit an unrealistic pattern near the gluing of the two representations of the crack, due to the differentmaterial stifnesses, see Figure 8. Consequently, continuity will be imposed on ˆΓ in all examples.
Consider an L-shaped panel with boundary conditions as illustrated in Figure 9. Following Ambati et al[1], the material parameters are E = 25 .
84 GPa, ν = 0 .
18 and G C = 8 . · − kN/mm. The phase-fieldlength-scale parameter is l = 2 . u D = 10 − mm.The background mesh is a uniform quadrilateral mesh with element size h = 10 mm, and elements inΩ tips are refined with refinement factor m = 20. The subdomain Ω tips initially consists of the 3 elements inthe corner of the piece, where crack inception is expected. The distance of derefinement in the switchingcriterion is taken as δ ∗ = 2 h . Recall that a correct initial definition of Ω tips , including crack tips andnotches of the domain, is essential since the damage field is computed only in this part of the domain.As a reference solution, we run the simulation with a plain phase-field model with adaptivity, following15F PF - XFEM . mm . mm . mm Figure 10:
L-shaped test.
Damage field at different load steps for PF and PF-XFEM. The sharp crack inPF-XFEM is plotted in white. Zoom into [ − , × [ − ,
60] mm . Displacement [mm] F o r c e [ k N ] PFPF - XFEM
Figure 11:
L-shaped test.
Load-displacement curve for PF and PF-XFEM.
Displacement [mm] % n D O F PFPF - XFEM
Figure 12:
L-shaped test.
Evolution of nDOF in the equilibrium equation for PF and PF-XFEM. Thepercentage is computed with respect to the nDOF of the initial discretization.16he refinement strategy proposed in [21]. The equations of the hybrid phase-field model are solved in thewhole domain. Elements along cracks are dynamically refined, but no transition to a sharp representationis considered. Thus, more elements become refined as the crack propagates. On the contrary, in thenewly proposed approach, only elements near the crack tip are refined. Thus, the size of Ω tips does notincrease constantly. Elements which are far enough from the crack tip transition to an XFEM coarserrepresentation, and the damage band is replaced by a sharp crack. Remeshing is avoided in both cases.In what follows, we abbreviate these approaches by PF and PF-XFEM, respectively.The crack path obtained with the two strategies at different load steps is depicted in Figure 10. Weobtain similar results, with a slightly faster propagation for PF-XFEM at load step u D = 0 .
26 mm. Thiscan be explained by the spurious transmission of forces across phase-field cracks, while sharp cracks arecompletely traction-free. In both cases, the discretization is updated according to crack growth.The slightly faster propagation when introducing sharp cracks is also observed in the load-displacementcurves in Figure 11. After the peak, when the crack starts propagating, we observe a steeper descent of thecurve, meaning the crack is slightly longer. As the simulations evolve, the difference between approachesdiminishes. In the PF curve, we observe the characteristic loss of stiffness prior to the peak of thesemodels; solving the damage equation in the whole domain with a quadratic degradation function g ( d )propagates the damage, softening the piece in the simulation. In PF-XFEM the behavior of the materialbefore cracking is closer to linear elasticity because the damage is computed only in small subdomains.For PF-XFEM, Ω tips takes a very small part of the domain. This leads to a reduction of the numberof degrees of freedom (nDOF) for the equilibrium equation. The reduction is substantial even though inPF the mesh is refined only in a narrow band containing the crack. Figure 12 shows the evolution ofnDOF with respect to its initial value. For PF, the number increases up to 800%, while for PF-XFEMnDOF the value increases only up to 150%. This example was first proposed in Muix´ı et al [22] and is revisited here to test the capability of thestrategy to reproduce bifurcations, as well as to handle multiple disconnected subdomains in Ω tips .We consider a square plate in [ − , mm , precracked at midheight and with boundary conditionsas illustrated in Figure 13. The piece is fixed on its right edge and has imposed vertical displacements onits top and bottom edges, with f ( x ) = u D ( x − /
8. The material parameters are E = 20 GPa, ν = 0 . G C = 8 . · − kN/mm. We take l = 0 .
01 mm and load increments of ∆ u D = 5 · − mm. Thedomain is discretized into a uniform quadrilateral mesh of 45 ×
45 elements. Elements are refined withrefinement factor m = 15. The threshold distance in the switching criterion is taken as δ ∗ = 3 h . Theinitial Ω tips consists of the elements containing the preexisting crack and the rest of the domain is part ofΩ xfem .The obtained crack at different stages is shown in Figure 14. The crack propagates horizontally andthen branches abruptly, maintaining the symmetry of the solution in the whole simulation. When thebranching occurs, we start having two crack tips and Ω tips contains the elements near both of them. InΩ xfem we have two separate sharp cracks, each one of them contributing to the XFEM discretization withan independent Heaviside function.The evolution of the crack at the load step when it branches, u D = 0 . g ( d ) in (3c), but when transitioningto the XFEM representation the crack faces intersect as shown in Figure 16. This issue can be tackledby implementing contact conditions into the XFEM discretization, see for instance [14, 9].17igure 13: Branching test.
Geometry and boundary conditions. Dimensions in mm. u D = 0 .
030 mm u D = 0 .
055 mm u D = 0 .
080 mmFigure 14:
Branching test.
Evolution of the crack at different load steps. The sharp crack is plotted inwhite. s = 0 s = 20 s = 40Figure 15: Branching test.
Detail of the crack path at some staggered iterations for imposed displacement u D = 0 . , . × [ − . , .
4] mm .18igure 16: Branching test.
Deformed piece at imposed displacement u D = 0 .
08 mm. Crack faces,highlighted in red, show a slight interpenetration. Zoom into [0 . , . × [ − . , .
2] mm .Figure 17: Multiple cracks test.
Geometry and boundary conditions.
This example, proposed in [21], exhibits the robustness of the proposed strategy to simulate coalescenceof cracks.Consider a square plate in [0 , mm , with 6 initial cracks and subjected to biaxial tension as depictedin Figure 17. The tips of the initial cracks are in Table 1. Material parameters are E = 20 GPa, ν = 0 . G C = 10 − kN/mm. We use l = 0 .
012 mm and take increments of ∆ u D = 5 · − mm. Thebackground mesh is a uniform quadrilateral mesh with 47 ×
47 elements and the refinement factor is m = 17. The distance to switch to discontinuous is δ ∗ = 3 h .The obtained crack pattern is plotted in Figure 18. The initial cracks propagate, coalescing betweenthem, until the piece is broken into 4 independent pieces. There are two mergings of cracks in theprocess. Since cracks grow abruptly, we need to plot the crack at staggered iterations to see the detailof how coalescence is handled by the method. Figure 19 shows some representative staggered iterationsfor imposed displacement u D = 0 .
262 mm. In the figure, we can see how the sharp crack in the right isreplaced by its diffuse representation as the crack tip approaches and the cut elements are refined. Then,after the merging has occured, all the elements in this subdomain of Ω tips transition back to Ω xfem .19able 1:
Multiple cracks test.
Tip coordinates for the initial cracks in the domain [0 , mm .Crack P (mm) P (mm)1 (0 . , .
5) (0 . , . , .
1) (1 , . . , .
5) (1 . , . . , .
9) (0 . , . . , .
5) (0 . , . . , .
5) (1 . , . u D = 0 . u D = 0 . u D = 0 . u D = 0 .
019 mm u D = 0 .
027 mmFigure 18:
Multiple cracks test.
Damage field at several imposed displacements.Now, we compare the results with the ones obtained by a PF approach. The final partition of thepiece for both approaches is shown in Figure 20. The damage field corresponds to the PF simulation,and the final sharp crack for PF-XFEM is plotted in black. The crack patterns are very similar also inthis more complex scenario. The corresponding load-displacement curves for the horizontal and verticalloads, F x and F y , in Figure 21, again indicate a loss of stiffness for PF in the precracking regime and asteeper descent of the curves for PF-XFEM. The evolution of nDOF is plotted in Figure 22, showing asubstantial decrease in nDOF. For PF-XFEM we observe a reduction of nDOF as crack tips disappear,while in PF the nDOF always increases. 20 = 30 s = 60 s = 67 s = 68Figure 19: Multiple cracks test.
Crack pattern at different staggered iterations for imposed displacement u D = 0 .
262 mm. Zoom into [0 , . × [0 . , .
1] mm .Figure 20: Multiple cracks test.
Comparison with PF solution.
Displacement [mm] F o r c e [ k N ] F x PF F y PFF x PF-XFEM F y PF-XFEM
Figure 21:
Multiple cracks test.
Load-displacement curve.21
Displacement [mm] % n D O F PFPF - XFEM
Figure 22:
Multiple cracks test.
Evolution of the number of degrees of freedom.Figure 23:
Twisting test.
Front and top view of the initial geometry of the piece. Lenghts are in mm.
In this example we test the performance of the method in a fully 3D example. This test was alsoproposed in [21]. We consider a square section beam in [0 , × [0 , × [0 ,
25] mm , which is clampedon { x = 0 mm } and has imposed displacements along x on { x = 125 mm } . The beam has two initialnotches on { y = 0 mm } and { y = 25 mm } , which are inclined with opposite angles, as shown in Figure23. The parameters are E = 32 GPa, ν = 0 . G C = 1 . · − kN/mm and l = 2 mm. The load stepstake increments ∆ u D = 5 · − mm. We use a uniform hexahedral mesh with element size h = 3 .
125 mmand refinement factor m = 7. The distance in the switching criterion is δ ∗ = 2 h . The initial notches aremodeled by diffuse cracks and Ω tips is the union of elements that contain them.In this case, the initial cracks coalesce, completely splitting the beam at a single load step, as canbe seen in the load-displacement curve in Figure 24. Figure 25 shows the final geometry at load step u D = 0 . Displacement [mm] F o r c e [ k N ] Figure 24:
Twisting test.
Load-displacement curve. The piece completely breaks at u D = 0 . Twisting test.
Geometry of the piece at imposed displacement u D = 0 . Conclusions
We propose a novel method to simulate fracture which is based on combining a phase-field model insmall subdomains around crack tips, Ω tips , and a discontinuous model in the rest of the domain, Ω xfem .The approach overcomes the limitations of both continuous and discontinuous models. Propagation isdescribed by the phase-field solution in Ω tips , while an XFEM approximation explicitly describes the crackopening and enables to use a coarser discretization in almost the whole domain. In all examples, a correctdefinition of the initial partition in Ω tips and Ω xfem is crucial to detect crack inception.Computationally, the same background mesh is used during the whole simulation. Refined elementsare nested in Ω tips to capture the phase-field solution and the sharp cracks are introduced via XFEM. Thediscretization is automatically updated as cracks propagate and remeshing is avoided. Nitche’s methodis used to impose continuity of displacements in weak form to mantain a very local refinement with notransitioning elements.The robustness of the strategy has been proved in 2D and 3D, including scenarios with branching andcoalescence. The obtained results are comparable to the ones from a plain phase-field approach, but withan important reduction in the number of degrees of freedom. Also, the response of the pieces in the linearelastic regime is more realistic and transmission of forces across cracks is prevented. The methodology isan efficient alternative for fracture simulation.Further investigation would include enhancing the method by incorporating contact conditions intothe XFEM formulation and by exploring other alternatives for the geometrical description of sharp cracks.
Acknowledgements
This work was supported by the Ag`encia de Gesti´o d’Ajuts Universitaris i de Recerca training grantFI-DGR 2017, the DAFOH2 project (Ministerio de Ciencia e Innovaci´on, MTM2013-46313-R) and theDepartament d’Innovaci´o, Universitats i Empresa, Generalitat de Catalunya (2017-SGR-1278).
References [1] M. Ambati, T. Gerasimov, and L. De Lorenzis. A review on phase-field models of brittle fractureand a new fast hybrid formulation.
Comput Mech , 55:383–405, 2015.[2] T. Belytschko and T. Black. Elastic crack growth in finite element with minimal remeshing.
Int JNumer Methods Eng , 45(5):601–620, 1999.[3] M. J. Borden, C. V. Verhoosel, M. A. Scott, T. J. R. Hughes, and C. M. Landis. A phase-fielddescription of dynamic brittle fracture.
Comput Methods Appl Mech Eng , 217-220:77–95, 2012.[4] K. J. Chang. On the maximum strain criterion – a new approach to the angled crack problem.
EngFract Mech , 14(1):107–124, 1981.[5] F. Erdogan and G. C. Sih. On the crack extension in plates under plane loading and transverse shear.
J Basic Eng , 85(4):519–525, 1963.[6] T.P. Fries and T. Belytschko. The extended/generalized finite element method: An overview of themethod and its applications.
Int J Numer Methods Eng , 84(3):253–304, 2010.[7] R. Geelen, Y. Liu, J. Dolbow, and A. Rodr´ıguez-Ferran. An optimization-based phase-field methodfor continuous-discontinuous crack propagation.
Int J Numer Methods Eng , 116(1):1–20, 2018.[8] R. Geelen, J. Plews, M. Tupek, and J. Dolbow. An extended/generalized phase-field finite elementmethod for crack growth with global-local enrichment.
Int J Numer Methods Eng , 121(11):2534–2557,2020. 249] E. Giner, M. Tur, J. Taranc´on, and F. Fuenmayor. Crack face contact in X-FEM using a segment-to-segment approach.
Int J Numer Methods Eng , 82(11):1424–1449, 2009.[10] B. Giovanardi, A. Scotti, and L. Formaggia. A hybrid XFEM-Phase field (
Xfield ) method for crackpropagation in brittle elastic materials.
Comput Methods Appl Mech Eng , 320:396–420, 2017.[11] M. Griebel and M. A. Schweitzer. A Particle-Partition of Unity Method. Part V: Boundary Condi-tions. In
Stefan Hildebrandt and Hermann Karcher (eds) Geometric Analysis and Nonlinear PartialDifferential Equations, Springer, Berlin , pages 519–542, 2003.[12] C. G¨urkan, E. Sala-Lardies, M. Kronbichler, and S. Fern´andez-M´endez. eXtended HybridizableDiscontinous Galerkin (X-HDG) for void problems.
J Sci Comput , 66:1313–1333, 2016.[13] P. Hansbo. Nitsche’s method for interface problems in computational mechanics.
GAMM-Mitteilungen , 28(2):183–206, 2005.[14] T. Y. Kim, J. Dolbow, and T. Laursen. A mortared finite element method for frictional contact onarbitrary interfaces.
Comput Mech , 39:223–235, 2007.[15] P. Laborde, J. Pommier, Y. Renard, and M. Sala¨un. High-order extended finite element method forcracked domains.
Int J Numer Methods Eng , 64(3):354–381, 2005.[16] O. Marco, R. Sevilla, Y. Zhang, J. J. R´odenas, and M. Tur. Exact 3D boundary representation infinite element analysis based on Cartesian grids independent of the geometry.
Int J Numer MethodsEng , 103(6):445–468, 2015.[17] C. Miehe, M. Hofacker, and F. Welschinger. A phase-field model for rate-independent crack prop-agation: robust algorithmic implementation based on operator splits.
Comput Methods Appl MechEng , 199(45):2765–2778, 2010.[18] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models offracture: variational principles and multi-field FE implementations.
Int J Numer Methods Eng ,83(10):1273–1311, 2010.[19] N. Mo¨es and T. Belytschko. Extended finite element method for cohesive crack growth.
Eng FractMech , 69:813–833, 2002.[20] N. Mo¨es, J. Dolbow, and T. Belytschko. A finite element method for crack growth without remeshing.
Int J Numer Methods Eng , 46(1):131–150, 1999.[21] A. Muix´ı, S. Fern´andez-M´endez, and A. Rodr´ıguez-Ferran. Adaptive refinement for phase-field modelsof brittle fracture based on Nitsche’s method.
Comput Mech , 66:69–85, 2020.[22] A. Muix´ı, A. Rodr´ıguez-Ferran, and S. Fern´andez-M´endez. A Hybridizable Discontinuous Galerkinphase-field model for brittle fracture with adaptive refinement.
Int J Numer Methods Eng ,121(6):1147–1169, 2020.[23] S. Nagaraja, M. Elhaddad, M. Ambati, S. Kollmannsberger, L. De Lorenzis, and E. Rank. Phase-field modeling of brittle fracture with multi-level hp -FEM and the finite cell method. Comput Mech ,63:1283–1300, 2019.[24] M. Ortiz and A. Pandolfi. Finite-deformation irreversible cohesive elements for three-dimensionalcrack-propagation analysis.
Int J Numer Methods Eng , 44(9):1267–1282, 1999.[25] R. U. Patil, B. K. Mishra, and I. V. Singh. A local moving extended phase field method (LMXPFM)for failure analysis of brittle fracture.
Comput Methods Appl Mech Eng , 342:674–709, 2018.2526] G. C. Sih. Strain-energy-density factor applied to mixed mode crack problems.
Int J Fract , 10(3):305–321, 1974.[27] N. Sukumar, J. Dolbow, and N. Mo¨es. Extended finite element method in computational fracturemechanics: a retrospective examination.
Int J Fract , 196:189–206, 2015.[28] C. T. Sun and Z. H. Jin. Fracture Mechanics (1st edition).
Academic Press , 2012.[29] E. Tamayo-Mas and A. Rodr´ıguez-Ferran. A medial-axis-based model for propagating cracks in aregularised bulk.