# Computational Methods for Fracture in Brittle and Quasi-Brittle Solids: State-of-the-Art Review and Future Perspectives

Authors: Timon Rabczuk
Journal: ISRN Applied Mathematics
 (2013)

---

## Abstract

An overview of computational methods to model fracture in brittle and quasi-brittle materials is given. The overview focuses on continuum models for fracture. First, numerical difficulties related to modelling fracture for quasi-brittle materials will be discussed. Different techniques to eliminate or circumvent those difficulties will be described subsequently. In that context, regularization techniques such as nonlocal models, gradient enhanced models, viscous models, cohesive zone models, and smeared crack models will be discussed. The main focus of this paper will be on computational methods for discrete fracture (discrete cracks). Element erosion technques, inter-element separation methods, the embedded finite element method (EFEM), the extended finite element method (XFEM), meshfree methods (MMs), boundary elements (BEMs), isogeometric analysis, and the variational approach to fracture will be reviewed elucidating advantages and drawbacks of each approach. As tracking the crack path is of major concern in computational methods that preserve crack path continuity, one section will discuss different crack tracking techniques. Finally, cracking criteria will be reviewed before the paper ends with future research perspectives.


---

## Body

## 1. Introduction

The prediction of material failure is of major importance in engineering and materials Science. In the past two decades, a huge effort was made to develop novel, efficient, and accurate computational methods for fracture, and an enormous progress was made. This paper will give an overview on those advances. The paper focuses on continuum models for brittle and quasi-brittle fracture. It will not discuss issues related to ductile fracture with plastic deformation prior to localization; see for example, the excellent work of the group of Prof. Li [1–7].There are three failure modes see Figure1. Mode-I failure is related to crack opening, mode-II failure (sliding) is a pure shear failure mode, and mode-III failure (tearing) can be considered as out-of-plane shearing. In many applications, materials will fail due to a mixed mode failure. Fracture in brittle materials can be modelled by linear elastic fracture mechanics (LEFMs). Brittle materials are characterized by linear elastic material behavior in the bulk and a small fracture process zone. However, LEFM cannot be used when the failure process zone is of the order of the size of the structure. The relative size of the fracture process zone (l) with respect to the smallest critical dimension of the structure (D) is important for the choice of the fracture model [8]. For D/l>100, LEFM is valid while for 5<D/l<100, a (nonlinear) quasi-brittle fracture approach, should be used instead; see also the contributions by [9, 10]. Gdoutos [8] recommends the use of nonlocal damage models for D/l<5. The length of the fracture process zone is of the order of the characteristic length [11]:
(1)lch=EGfft2orlch=EGf(1-ϑ2)ft2,
where E is the modulus of elasticity, ϑ is the Poisson’s ratio, Gf is the fracture energy and ft denotes the tensile strength of the material. Carpinteri [12] introduced a nondimensional brittleness number sE=Gf/(ftb) where b is a geometric measure. He tested pre-notched beams under 3-point bending. The distance between the notch and the upper boundary of the beam is b. If the process zone is small compared to b, the failure is brittle and LEFM is applicable. Otherwise, in the absence of plasticity, we may call the failure quasi-brittle. The stress-strain curve for quasi-brittle materials such as concrete, ceramics, and rock is characterized by a drop of the stress after the peak stress is reached. This softening can be explained by damage mechanics due to the initiation, growth, and coalescences of microcracks that reduce the effective cross-sectional area, Figure 2. For more details, the reader is referred to numerous articles and books on damage mechanics [13–26].Figure 1
Different failure modes.Figure 2
The equivalence principe of continuum damage mechanics.The use of strain-softening material models leads to an ill-posed boundary value problem (BVP) [27–31]. The solution of ill-posed BVPs with computational methods leads to certain difficulties. Bažant and Belytschko [32] showed that the deformations localize in a set of measure zero; in 1D the localization will occur in a point, in 2D in a line, and in 3D in a surface. They further showed that the energy dissipation approaches zero with increasing refinement of the discretization yielding physically meaningless results.The difficulties that emerge with the material instabilities due to the ill-posedness of the BVP can be avoided by regularization techniques such as higher order continuum models, gradient based models, and polar theories (the most well-known theory is probably the Cosserat continuum), nonlocal models, viscous models or cohesive zone models. These regularization techniques introduce a characteristic length into the discretization. In other words, the local character of the deformation is lost, and the fracture is “smeared” over a certain domain involving several elements. Therefore, a very fine mesh has to be used in order to resolve the crack. As the crack introduces a jump in the displacement field, methods that smear the crack over a certain region are not capable of representing the correct crack kinematics. Moreover, the different length scales (of the structure and the characteristic length) may significantly increase the computational cost. A lot of effort has been devoted to develop methods that are capable of capturing the crack and take advantage of regularization techniques that can be coupled to these strong discontinuity approaches. Modeling fracture with strong-discontinuity approaches require two key ingredients:(1)
a method that is capable of capturing the crack kinematics and(2)
a cracking criterion that determines the orientation and the length of the crack. Computational methods that describe the topology of the crack as continuous surface require also methods to describe the crack surface and track the crack path. The paper is structured as follows. In the next section, popular regularization techniques are summarized including gradient enhanced models, non-local models, viscous models, cohesive zone models, and smeared crack models. Section 3 gives an overview of the state-of-the-art computational methods for (discrete) fracture of brittle and quasi-brittle solids. Section 4 is devoted to crack tracking procedures as it is a major challenge for computational methods preserving crack path continuity. Section 5 discusses cracking criterion that determine the transition from the continuum to the discontinuum. The article finishes with future research directions and some concluding remarks.


## 2. Regularization

### 2.1. Gradient-Enhanced Models

Gradient-enhanced models, or briefly called gradient models, are typically described by differential equations that contain higher order spatial derivatives. The coefficients multiplying the terms of different order have different physical dimensions. It is possible to deduce the characteristic length from their ratio. Gradient models can be incorporated into plasticity format, [33, 34], damage format [22, 35] or in a combination of both [15]. Since they require higher order spatial derivatives, higher order elements need to be used. Gradient models are sometimes referred to as weakly non-local models, see Bažant and Jirásek [14] though some recent models [22] are also classified as strongly non-local.


### 2.2. Nonlocal Models

Strongly non-local models are models of the integral type, see Figure3. In non-local formulations, the stress at a given material point does not only depend on that local strain but also on the strain of the neighboring points. This effect is obtained by averaging certain internal variables in a given neighborhood, Figure 3. The size of this neighborhood is a material parameter, that depends also on the geometry of the structure of interest. In other words, changing the dimensions of the problem requires recalibration of the parameters of the non-local model. They can be obtained by an inverse analysis. Bažant and Jirásek [14] give an excellent overview on non-local models of the integral types and physical motivations (see also Bažant [36]) as well as suggestions for calibrating material parameters.Figure 3
Principle of nonlocal constitutive models with the typical bell-shaped domain of influence.


### 2.3. Viscous Models

The introduction of a viscosity can also restore the well-posedness of the BVP or initial BVP (IBVP). It can be regarded as introducing higher order time derivatives, similar to the gradient models. Considering the dimensions of the viscosityη (kg/(ms)), the dimensions of the Young’s modulus E (kg/(ms2)) and the dimension of the mass density ϱ(kg/m3), there is indeed a length scale lch associated with the viscosity given by [14]:
(2)lch=ηE/ϱ,
where c=E/ϱ is the longitudinal propagation velocity in 1D. However, the well-posedness of the IBVP is only guaranteed during the time span of the viscosity (t0=lch/c). For visco-plastic models, Etse and Willam [37] have shown that after discretization, hyperbolicity of the linearized momentum equation can be lost if a critical time step is exceeded.The introduction of a viscosity can sometimes be physically motivated. The strain rate effect and the corresponding dynamic strength, increase, for example, can be captured by viscous damage models [38].


### 2.4.  Cohesive Zone Models

Cohesive zone models (CZMs) also restore the well-posedness of the (I)BVP. In contrast to the models described before, CZMs can be combined with computational methods that maintain the local character of the crack. In cohesive cracks, a traction-separation model is applied across the crack surface that links the cohesive traction transmitted by the discontinuity surface to the displacement jump, characterized by the separation vector. CZMs go back to the 60’s and were originally developed from Dugdale [39] and Barenblatt [40]. They were applied in metal plasticity to take into account friction along neighboring grains. Hillerborg et al. [41] extended this concept to model crack growth in concrete. They called their approach fictitious crack model. The main difference between the fictitious crack model of Hillerborg and the CZM by Dugdale [39] and Barenblatt [40] is that crack initiation and propagation is not restricted along a predetermined path but cracks can initiate anywhere in the structure.The principal idea of cohesive cracks is shown in Figure4. The cohesive model is used in the so-called process zone, sometimes also called cohesive zone. The complex stress state around the crack tip is lumped into a single surface. The first approach by Hillerborg et al. [41] was limited to mode-I fracture. Meanwhile, many models have been developed that are able to handle mixed mode fracture and other complex phenomena including irreversible deformations, stress triaxiality, and rate dependence [42–47]. Some CZMs include fine scale effects obtained from micromechanical models or atomistic simulations [48–51]. Note that the methods developed in [52–54] are different to conventional cohesive zone models as the interphase zone is a finite size zone and therefore enables the simulation of complex mixed-mode fracture and thermomechanical fracture. They also allow the combination of the Cauchy-Born rule with cohesive strength.Figure 4
Principle of cohesive crack models.CZMs can be classified into initially rigid models (sometimes also called extrinsic models), Figure5(d), and initially elastic models, Figure 5(e). While the initially rigid models are characterized by a monotonic decrease in the cohesive traction, initially elastic models exhibit initially a positive slope, Figure 5(e). As the cohesive traction increases in initially elastic models, the stresses in the neighboring areas will increase as well leading to spurious cracking. Cohesive cracking will extend to a finite zone with infinitely close cracks with infinitely small crack openings, which means that localized cracking is in fact forbidden [55].Typical one-dimensional cohesive models (a) The model by Dugdale1 for fracture across grain boundaries in metal, (b) a modified Dugdale model, (c) linear decaying cohesive model, (d) exponential decaying model, (e) model with positive initial stiffness for brittle fracture, and (f) non-admissible cohesive model.
(a)
(b)
(c)
(d)
(e)
(f)Initially rigid models cause numerical difficulties when elastic unloading occurs at an early stage such that the stiffness for the unloading case tends to infinity. It is believed that initially rigid models are better suited particularly in the context of dynamic fracture. A good overview and a discussion on initially elastic and initially rigid models are given, for example, in Papoulia et al. [56] and references therein.An important material parameter of the CZM is the fracture energyGf. In the one-dimensional case of pure mode I-failure, the fracture energy Gf=∫0wmax‍t(w)dw is the energy that is dissipated during crack opening. It is related to the area of the cohesive crack; t is the cohesive traction, w is the crack opening and wmax is the crack opening at which no cohesive forces are transmitted; see also Figure 4.


### 2.5.  Smeared Crack Models

The main idea of smeared crack models is to spread the energy release along the width of the localization band usually within a single element so that it be objective. This is achieved by calibrating the width of the band such that the dissipated energy is the correct one. This introduces a characteristic length into the discretization that depends on the size of the elements.A milestone in smeared crack models is thecrack band model by Bažant and Oh [57]. It is based on the observation that the process zone localizes in a single element, the so-called crack band. Therefore, Bažant and Oh [57] modified the constitutive model after localization so that the correct energy is dissipated. Since the width of the numerically resolved fracture process zone depends on the element size and tends to zero if the mesh is refined, the crack band model cannot be considered as real localization limiter. It only partially regularizes the problem in the sense that the global characteristics do not exhibit spurious mesh dependence. When the element size exceeds the size of the process zone, the stress-strain curve has to be adjusted horizontally as shown in Figure 6(b). This makes the stress-strain curve in the post-localization region steeper, and the appropriate scaling ratio achieves the dissipation energy remains the same. The crack band model was successfully used for mode-I fracture. The extension to mixed mode failure and complex three-dimensional stress states is difficult. Besides, exact scaling leads to a system of non-linear equations that can only be solved iteratively. Moreover, when the elements are too large (larger than the physical process zone), the adjusted post-localization diagram will develop a snap back, see Figure 6(c). In turn, when the elements are much smaller than the process zone, the adjusted post-localization diagram becomes less steep, and the band of cracking elements converges to a set of measure zero. Thus, the limit corresponds to the solution of a cohesive crack model [58].Adjustment of the stress-strain curve in the crack band model, (a) master curve and decomposition of the strains into prelocalizaton and post-localization parts, (b) horizontal scaling for an element larger than the physical process zone, and (c) snap back for too large elements reproduced from Bažant and Jirásek [14].
(a)
(b)
(c)The smeared crack models can be further classified into rotating and fixed crack models. In thefixed crack model, the orientation of the crack is not allowed to change. Fixed crack models can be again classified into single fixed crack models and multiple fixed crack models. In the latter case, more than one crack initiates in an element. Often, a second crack is only allowed to be initiated perpendicular to the first crack (this should avoid the initiation of too many cracks), see Figure 7(b).Principle of (a) rotating crack models, (b) fixed crack models.
(a)
(b)The orientation of the crack is typically chosen to be perpendicular to the direction of the principal tensile strain or principal tensile stress. Since the direction of the principal tensile strain changes during loading, fixed crack models result in too stiff system responses. Fixing the crack orientation leads to stress locking, meaning stresses are transmitted even over wide open cracks mainly caused due to shear stresses generated by a rotation of the principal strain axes after the crack initiation. Therefore,rotating crack models have been introduced where the crack is rotated with the principal strain axis; see Figure 7(a). However, Jirásek and Zimmermann [59] have shown that even rotating crack models suffer from stress locking caused by the poor kinematic representation of the discontinuous displacement field around the macroscopic crack. Unless the direction of the macroscopic crack that is represented by a band of cracked elements, is parallel to the element sides, the directions of the maximum principal strain determined from the finite element interpolation at individual material points deviate from the normal to the macroscopic crack. The lateral principal stress has a nonzero projection on the crack normal, which generates cohesive forces acting across the crack even at very late stages of the cracking process when the crack should be completely stress-free. Therefore, Jirásek and Zimmermann [60] proposed to combine rotating crack models with scalar damage models. Since in scalar damage models, all stress components tend to zero when the crack opens widely, no spurious stress transfer occurs any more.Some models combine the rotating and fixed smeared crack approach [61]. Up to a certain crack opening, the crack is allowed to rotate. Then, the crack orientation is fixed. Other interesting contributions can be found, for example, in [17, 61–65].Smeared crack models in finite element analysis often exhibit a so-called mesh alignment sensitivity or mesh orientation bias, see Figure8; meaning the orientation of the smeared crack depends on the orientation of the discretization. When smeared crack models are used in the context of meshfree methods with isotropic (circular) support, the “mesh” orientation bias can be alleviated.Mesh orientation bias for smeared crack models.
(a)
(b)


## 2.1. Gradient-Enhanced Models

Gradient-enhanced models, or briefly called gradient models, are typically described by differential equations that contain higher order spatial derivatives. The coefficients multiplying the terms of different order have different physical dimensions. It is possible to deduce the characteristic length from their ratio. Gradient models can be incorporated into plasticity format, [33, 34], damage format [22, 35] or in a combination of both [15]. Since they require higher order spatial derivatives, higher order elements need to be used. Gradient models are sometimes referred to as weakly non-local models, see Bažant and Jirásek [14] though some recent models [22] are also classified as strongly non-local.


## 2.2. Nonlocal Models

Strongly non-local models are models of the integral type, see Figure3. In non-local formulations, the stress at a given material point does not only depend on that local strain but also on the strain of the neighboring points. This effect is obtained by averaging certain internal variables in a given neighborhood, Figure 3. The size of this neighborhood is a material parameter, that depends also on the geometry of the structure of interest. In other words, changing the dimensions of the problem requires recalibration of the parameters of the non-local model. They can be obtained by an inverse analysis. Bažant and Jirásek [14] give an excellent overview on non-local models of the integral types and physical motivations (see also Bažant [36]) as well as suggestions for calibrating material parameters.Figure 3
Principle of nonlocal constitutive models with the typical bell-shaped domain of influence.


## 2.3. Viscous Models

The introduction of a viscosity can also restore the well-posedness of the BVP or initial BVP (IBVP). It can be regarded as introducing higher order time derivatives, similar to the gradient models. Considering the dimensions of the viscosityη (kg/(ms)), the dimensions of the Young’s modulus E (kg/(ms2)) and the dimension of the mass density ϱ(kg/m3), there is indeed a length scale lch associated with the viscosity given by [14]:
(2)lch=ηE/ϱ,
where c=E/ϱ is the longitudinal propagation velocity in 1D. However, the well-posedness of the IBVP is only guaranteed during the time span of the viscosity (t0=lch/c). For visco-plastic models, Etse and Willam [37] have shown that after discretization, hyperbolicity of the linearized momentum equation can be lost if a critical time step is exceeded.The introduction of a viscosity can sometimes be physically motivated. The strain rate effect and the corresponding dynamic strength, increase, for example, can be captured by viscous damage models [38].


## 2.4.  Cohesive Zone Models

Cohesive zone models (CZMs) also restore the well-posedness of the (I)BVP. In contrast to the models described before, CZMs can be combined with computational methods that maintain the local character of the crack. In cohesive cracks, a traction-separation model is applied across the crack surface that links the cohesive traction transmitted by the discontinuity surface to the displacement jump, characterized by the separation vector. CZMs go back to the 60’s and were originally developed from Dugdale [39] and Barenblatt [40]. They were applied in metal plasticity to take into account friction along neighboring grains. Hillerborg et al. [41] extended this concept to model crack growth in concrete. They called their approach fictitious crack model. The main difference between the fictitious crack model of Hillerborg and the CZM by Dugdale [39] and Barenblatt [40] is that crack initiation and propagation is not restricted along a predetermined path but cracks can initiate anywhere in the structure.The principal idea of cohesive cracks is shown in Figure4. The cohesive model is used in the so-called process zone, sometimes also called cohesive zone. The complex stress state around the crack tip is lumped into a single surface. The first approach by Hillerborg et al. [41] was limited to mode-I fracture. Meanwhile, many models have been developed that are able to handle mixed mode fracture and other complex phenomena including irreversible deformations, stress triaxiality, and rate dependence [42–47]. Some CZMs include fine scale effects obtained from micromechanical models or atomistic simulations [48–51]. Note that the methods developed in [52–54] are different to conventional cohesive zone models as the interphase zone is a finite size zone and therefore enables the simulation of complex mixed-mode fracture and thermomechanical fracture. They also allow the combination of the Cauchy-Born rule with cohesive strength.Figure 4
Principle of cohesive crack models.CZMs can be classified into initially rigid models (sometimes also called extrinsic models), Figure5(d), and initially elastic models, Figure 5(e). While the initially rigid models are characterized by a monotonic decrease in the cohesive traction, initially elastic models exhibit initially a positive slope, Figure 5(e). As the cohesive traction increases in initially elastic models, the stresses in the neighboring areas will increase as well leading to spurious cracking. Cohesive cracking will extend to a finite zone with infinitely close cracks with infinitely small crack openings, which means that localized cracking is in fact forbidden [55].Typical one-dimensional cohesive models (a) The model by Dugdale1 for fracture across grain boundaries in metal, (b) a modified Dugdale model, (c) linear decaying cohesive model, (d) exponential decaying model, (e) model with positive initial stiffness for brittle fracture, and (f) non-admissible cohesive model.
(a)
(b)
(c)
(d)
(e)
(f)Initially rigid models cause numerical difficulties when elastic unloading occurs at an early stage such that the stiffness for the unloading case tends to infinity. It is believed that initially rigid models are better suited particularly in the context of dynamic fracture. A good overview and a discussion on initially elastic and initially rigid models are given, for example, in Papoulia et al. [56] and references therein.An important material parameter of the CZM is the fracture energyGf. In the one-dimensional case of pure mode I-failure, the fracture energy Gf=∫0wmax‍t(w)dw is the energy that is dissipated during crack opening. It is related to the area of the cohesive crack; t is the cohesive traction, w is the crack opening and wmax is the crack opening at which no cohesive forces are transmitted; see also Figure 4.


## 2.5.  Smeared Crack Models

The main idea of smeared crack models is to spread the energy release along the width of the localization band usually within a single element so that it be objective. This is achieved by calibrating the width of the band such that the dissipated energy is the correct one. This introduces a characteristic length into the discretization that depends on the size of the elements.A milestone in smeared crack models is thecrack band model by Bažant and Oh [57]. It is based on the observation that the process zone localizes in a single element, the so-called crack band. Therefore, Bažant and Oh [57] modified the constitutive model after localization so that the correct energy is dissipated. Since the width of the numerically resolved fracture process zone depends on the element size and tends to zero if the mesh is refined, the crack band model cannot be considered as real localization limiter. It only partially regularizes the problem in the sense that the global characteristics do not exhibit spurious mesh dependence. When the element size exceeds the size of the process zone, the stress-strain curve has to be adjusted horizontally as shown in Figure 6(b). This makes the stress-strain curve in the post-localization region steeper, and the appropriate scaling ratio achieves the dissipation energy remains the same. The crack band model was successfully used for mode-I fracture. The extension to mixed mode failure and complex three-dimensional stress states is difficult. Besides, exact scaling leads to a system of non-linear equations that can only be solved iteratively. Moreover, when the elements are too large (larger than the physical process zone), the adjusted post-localization diagram will develop a snap back, see Figure 6(c). In turn, when the elements are much smaller than the process zone, the adjusted post-localization diagram becomes less steep, and the band of cracking elements converges to a set of measure zero. Thus, the limit corresponds to the solution of a cohesive crack model [58].Adjustment of the stress-strain curve in the crack band model, (a) master curve and decomposition of the strains into prelocalizaton and post-localization parts, (b) horizontal scaling for an element larger than the physical process zone, and (c) snap back for too large elements reproduced from Bažant and Jirásek [14].
(a)
(b)
(c)The smeared crack models can be further classified into rotating and fixed crack models. In thefixed crack model, the orientation of the crack is not allowed to change. Fixed crack models can be again classified into single fixed crack models and multiple fixed crack models. In the latter case, more than one crack initiates in an element. Often, a second crack is only allowed to be initiated perpendicular to the first crack (this should avoid the initiation of too many cracks), see Figure 7(b).Principle of (a) rotating crack models, (b) fixed crack models.
(a)
(b)The orientation of the crack is typically chosen to be perpendicular to the direction of the principal tensile strain or principal tensile stress. Since the direction of the principal tensile strain changes during loading, fixed crack models result in too stiff system responses. Fixing the crack orientation leads to stress locking, meaning stresses are transmitted even over wide open cracks mainly caused due to shear stresses generated by a rotation of the principal strain axes after the crack initiation. Therefore,rotating crack models have been introduced where the crack is rotated with the principal strain axis; see Figure 7(a). However, Jirásek and Zimmermann [59] have shown that even rotating crack models suffer from stress locking caused by the poor kinematic representation of the discontinuous displacement field around the macroscopic crack. Unless the direction of the macroscopic crack that is represented by a band of cracked elements, is parallel to the element sides, the directions of the maximum principal strain determined from the finite element interpolation at individual material points deviate from the normal to the macroscopic crack. The lateral principal stress has a nonzero projection on the crack normal, which generates cohesive forces acting across the crack even at very late stages of the cracking process when the crack should be completely stress-free. Therefore, Jirásek and Zimmermann [60] proposed to combine rotating crack models with scalar damage models. Since in scalar damage models, all stress components tend to zero when the crack opens widely, no spurious stress transfer occurs any more.Some models combine the rotating and fixed smeared crack approach [61]. Up to a certain crack opening, the crack is allowed to rotate. Then, the crack orientation is fixed. Other interesting contributions can be found, for example, in [17, 61–65].Smeared crack models in finite element analysis often exhibit a so-called mesh alignment sensitivity or mesh orientation bias, see Figure8; meaning the orientation of the smeared crack depends on the orientation of the discretization. When smeared crack models are used in the context of meshfree methods with isotropic (circular) support, the “mesh” orientation bias can be alleviated.Mesh orientation bias for smeared crack models.
(a)
(b)


## 3. Computational Methods for Fracture

### 3.1. Element Erosion

The easiest way to deal with discrete fracture is the element erosion algorithm [66–69]. Element erosion algorithms do not require any representation of the crack’s topology. Most approaches do not really delete the elements but set the stresses in the “deleted” elements to zero. Approaches that indeed remove elements often replace the deleted elements with rigid masses. In order to account for the correct energy dissipation in the post-localization domain, concepts from smeared crack models are applied. Song et al. [70] compared different computational methods for certain benchmark problems in dynamic fracture. They report extreme mesh sensitivity for element erosion algorithms and conclude they are not well suited for dynamic fracture. However, due to their simplicity, they are incorporated in commercial software packages such as LS DYNA.For brittle fracture, Pandolfi and Ortiz [71] derived an interesting eigenerosion approach. Their idea is based on the eigenfracture approach by Schmidt et al. [72], where the eigendeformation is restricted in a binary sense. The element can be either intact or fractured. In the first case, the element behavior is elastic. Otherwise, the element is eroded. They showed that in contrast to many other element erosion approaches [66, 67, 73, 74], the eigenerosion approach converges to Griffith fracture with increasing mesh refinement.


### 3.2. Interelement-Separation Methods

Standard finite elements have difficulties to capture the crack kinematics since they use continuous trial functions that are not particularly well adapted for solutions with discontinuous displacement fields. One of the first models capable of modelling cracks within the FEM are the so-calledinterelement-separation methods [46, 47, 75–80]. Interelement seperation methods were mostly applied to dynamic fracture with cohesive cracks. In the approaches by Xu and Needlemen [75, 76], cohesive surfaces were introduced at the beginning of the simulation. In contrast, Camacho and Ortiz [46] adaptively introduced cohesive surfaces at element edges when a certain cracking criteria is met. In interelement-separation methods, cracks are only allowed to develop along existing interelement edges. This endows the method with comparative simplicity, but can result in an overestimate of the fracture energy when the actual crack paths are not coincident with element edges. The results depend not only on the mesh size (and form of the chosen element), but also on the mesh bias that can be compensated only by remeshing that is particularly cumbersome in 3D. It has been noted that the solutions sometimes depend on mesh refinement; see Falk et al. [81]. This sensitivity has been mollified by adding randomness to the strength, as in Zhou and Molinairi [80] and Espinosa et al. [82]. For dynamic problems, interelement-separation methods are often used with rigid cohesive zone models. Special care has to be taken in order to guarantee time continuity. Time continuity requires temporal convergence when the time step tends to zero. Violating time continuity might lead to unstable or oscillatory responses, and the results severely depend on the chosen time step [56, 83].Many interesting problems have been studied by interelement-separation methods [47, 75, 79, 82–84]. They are most popular for dynamic fracture and problems involving numerous cracks with complex fracture patterns.Special crack tip elements for linear fracture problems were developed in the early 70’s [85–88]. They account for the crack tip singularity and were applied to stationary cracks. Recently, such approaches were incorporated in the smoothed finite element method by Liu et al. [89, 90]. These approaches were extended to anisotropic materials [91], interface cracks [92], and plastic fracture mechanics [93]. In [94, 95], these approaches were applied to crack propagation problems in two dimensions. Adaptive remeshing techniques were employed. As the crack is aligned to the crack boundary, these methods can also be classified as interelement-separation methods. One drawback is that only the r-term is accounted for in the crack tip element but not the discontinuity behind the crack tip.


### 3.3. Embedded Elements (EFEM)

In 1987, Ortiz et al. [77] modified the approximation of the strain field to capture weak discontinuities in finite elements to improve the resolution of shear bands. Therefore, they enriched the strain field to obtain the kinematic relations shown in Figure 9(a). Based on this idea, Belytschko et al. [96]  allowed two parallel weak discontinuity lines in a single element, Figure 9(b), so that the element was able to contain a band of localized strain. Dvorkin et al. [97] were the first who developed a method that was able to deal successfully with strong discontinuities in finite elements, the class of embedded elements (EFEM). The name comes from the fact that the localization zone is embedded in a single element; see Figure 9(c). This way, crack growth can be modeled without remeshing. This class of methods are much more flexible than schemes that allow discontinuities only at element interfaces, and it eliminates the need for continuous remeshing. Several different versions of the EFEM were developed. All techniques had in common that the approximation of the displacement field was enriched with additional parameters to capture the jump in the displacement field in a single element. These parameters, in the following called enrichment, were inherent to the cracked element and the jump in the displacement field depends only on this enrichment. The different methods were derived from various principles, starting from the extended principle of virtual work over Hellinger Reissner to Hu-Washizu variational principle, using an enhanced assumed strain (EAS) or B-bar format. One major advantage of the EFEM (compared to XFEM discussed in the next section) is that the additional unknowns could be condensed on the element level, so that discontinuities could be captured with very small changes of the existing code.Element with (a) one weak discontinuity, (b) two weak discontinuities, and (c) one strong discontinuity.
(a)
(b)
(c)The first version of the EFEM is often called statical optimal symmetric (SOS) since traction continuity is fulfilled, but it is not possible to capture the correct crack kinematics. The correct relative rigid body motion of the element is not guaranteed and it has been shown by Jirásek [98] that SOS formulations lead to stress locking. The kinematical optimal symmetric (KOS) version by Lotfi and Shing [99] ensures the correct crack kinematics but violates the traction continuity condition. Consequently, the criteria for the onset of localization written in terms of stresses in the bulk are no longer equivalent to the same criteria written in terms of the tractions on the discontinuity area; for example, for the Rankine criterion, the normal traction at the onset of localization should be equal to the tensile strength and the shear traction should be zero. This cannot be properly reproduced by the KOS formulation. The kinematical and statical optimal nonsymmetric (KSON) version of embedded elements (see Dvorkin et al. [97], Klisinski et al. [100]) guarantees traction continuity and the appropriate crack kinematics but leads to a non-symmetric stiffness matrix with all its disadvantages with respect to solving the linearized system of equations.The approximation of the displacement field in the EFEM is given by [101]
(3)uh(X)=∑I∈𝒮NI(X)u~J+ℳs(e)(X)⟦u^I(e)(X)⟧,
where u~J are the nodal parameters of the usual approximation, u^  is the enrichment, and the superimposed “(e)” indicates the enrichment on element level. The function ℳs(e) is given by
(4)ℳs(e)(X)={0∀(e)∉SHs(e)-ρ(e)∀(e)∈Sρ(e)=∑I=1Ne+NI+(X),
where Hs is the step function acting on the crack line S and Ne+ is the number of nodes of element (e) that belong to the domain Ω0+;  see Figure 9(c). The discretization of the strain field can be written as
(5)ϵh(X)=∑I∈𝒮‍(∇0NI(X)⊗u~I)S-(∇ρ(e)⊗⟦u^I(e)(X)⟧)S+ηs(e)k(⟦u^I(e)(X)⟧⊗n)S,
where the superimposed S denotes the symmetric part of a tensor, ηs(e)/k is a regularized Dirac delta function, and ηs(e) is a collocation function defined as
(6)ηs(e)={1∀X∈Sek0∀X∉Sek,
with the thickness k of the localization band. Considering the equilibrium equation in elastostatics and with the trial and test function of the structure of (3), the discrete equations in matrix form can be written as
(7)[Ku~u~(e)Ku~u^(e)Ku^u~(e)Ku^u^(e)]{u~(e)⟦u^I(e)⟧}={FIext0},
with
(8)Ku~u~(e)=∫Ω0BTCBdΩ0,Ku~u^(e)=∫Ω0BTCB^dΩ0,Ku^u~(e)=∫Ω0B^*TCBdΩ0,Ku^u^(e)=∫Ω0B^*TCB^dΩ0,
where C is the elasticity tensor, B is the B operator containing the derivatives of the shape functions and B^ is the B operator of the enrichment that depends on the embedded element formulation (SOS and KOS or KSON). For the KSON formulation, B^*≠B^, that will result in a non-symmetric stiffness matrix. For the SOS and KOS formulation, B^*=B^. The elemental enrichment ⟦u^I(e)⟧ can be condensed on the element level:
(9)⟦u^I(e)⟧=-[Ku^u~(e)]-1Ku^u^(e)u~(e),
that leads with the expression for ⟦u^I(e)⟧ to the system of equations
(10)K¯u=f,
with
(11)K¯=Ku~u~-Ku~u^Ku^u^-1Ku^u~.
The crack in embedded elements is modelled by a set of crack segments with piecewise constant crack opening, see Figure 10 though recent versions of the EFEM can handle linear crack opening [102]. The EFEM has been applied to many interesting problems [103–106] including piezoelectric materials [107], geomaterials [108], coupled formulations [109], finite strains [110], crack branching [111], and dynamic fracture [112], to name a few.(a) Piecewise constant crack opening in embedded elements and (b) linear crack opening for XFEM.
(a)
(b)Oliver et al. [113] compared embedded elements with the extended finite element method (that will be discussed in the next section) (XFEM) for several selected examples and showed that embedded elements can be as accurate as XFEM. However, examples in LEFM that could be compared to analytical solutions are missing.Methods that model the crack as continuous surface and therefore contain a “physical” crack tip are expected to be more accurate. However, the requirement to represent the crack surface is a drawback of those methods, in particular for complex crack patterns (e.g., in three-dimensions or for branching cracks). EFEM in principle does not require crack path continuity [114, 115]. However, some authors report that the elements then become sensitive to the orientation of the discontinuity surface with respect to the nodes, and in some unfavorable situations the response of the element ceases to be unique [116].


### 3.4. Extended Finite Element Method (XFEM)

A very flexible method that can handle linear and nonlinear crack openings is the extended finite element method (XFEM) developed by Belytschko et al. [11, 117, 118]. This method employs the local partition of unity concept [119–121] and introduces additional nodal parameters for the elements cut by the crack. Hence, the additional unknowns cannot be condensed on the element level as in EFEM. The displacement discontinuity depends only on the additional nodal parameters. The principle idea of XFEM is shown in Figure 11, simplified for a crack in one dimension. Though XFEM has originally been developed for fracture, it has been extended to numerous applications including two-phase flow [122, 123], fluid-structure interaction [124, 125], biomechanics [126], inverse problems [127, 128], multifield problems [129–132], among others. XFEM has been incorporated meanwhile into commercial software (XFEM) and has become one of the most popular methods for fracture. The generalized finite element method (GFEM) [120, 121, 133–139] is similar to XFEM.Figure 11
Principle of XFEM for strong discontinuities in 1D.The basic idea of XFEM is to decompose the displacement field into a continuous partucont and a discontinuous part udisc:
(12)uh(X)=ucont(X)+udisc(X).
The continuous part is the standard FE interpolation and additional information is introduced into the FE interpolation through the local partition of unity approach by adding an enrichment (udisc). Therefore, additional degrees of freedom are introduced for the discretization of the test- and trial functions. Those additional degrees of freedom are multiplied with the so-called enriched shape functions. They are the product of the “standard” element shape functions and enrichment functions that contain information about the solution. Note that the solution does not need to be known exactly. For example, a crack results in a discontinuous displacement field. Therefore, a discontinuous function is used as enrichment function. In linear elastic fracture mechanics (LEFM), the analytical near crack tip solution is known as well and can be included in the XFEM approximation. The approximation of the displacement field for nc cracks with mt crack tips reads
(13)uh(X)=∑I∈𝒮NI(X)uI+∑N=1nc∑I∈𝒮cNI(X)ψI(N)aI(N)+∑M=1mt∑I∈𝒮tNI(X)∑K=1NKϕKI(M)bKI(M),
where 𝒮 is the set of nodes in the entire discretization, 𝒮t is the set of nodes around the crack tip, and 𝒮c is the set of nodes associated to elements completely cut by the crack; ψI(N) is the enrichment function of crack N; ϕKI(M) is the enrichment function for the crack tip M; aI and bKI are additional degrees of freedom that needs to be solved for.For elements completely cut by the crack, the jump function is often chosen to cause a discontinuous displacement field(14)ψI(N)=sign[f(N)(X)]-sign[f(N)(XI)],
with
(15)f(N)(X)=sign[n·(X-X(N))]min(X-X(N))︸,XN∈Γc(N)
where n denotes the normal vector to the crack surface and f(N)(X) is a level set function; see Figure 12. The zero-isobar f=0 describes the position of the crack surface. Often a shifting is employed, see the second term on the RHS of (14) and Figure 11, and ensures that the enriched shape function NI(x)ψI(N) vanishes in the neighboring element. The shifting also maintains the interpolatory character of the XFEM approximation or in other words: the nodal parameters uI remain the physical displacements.Figure 12
Representation of the crack surface by level setsf and g.In LEFM, the crack tip enrichment is choosen according to the analytical solution:(16)ϕKI(M)=[rsin(θ2),rcos(θ2),rsin(θ2),sinθ,rcos(θ2),sinθ],
where r and θ are explained in Figure 12. Early approaches used a so-called topological enrichment, Figure 13(a). However, a topological enrichment leads to a deterioration of the convergence rate as the domain with the accurate tip enrichment shrinks when the mesh is refined. Therefore, it is important to use a geometrical enrichment, Figures 13(b) and 13(c). Note that a geometrical enrichment can drastically increase the condition number of the system matrix and requires the use of pre-conditioners [140] or other techniques [121, 141]. For non-linear materials and cohesive cracks, a tip enrichment should be used that avoids the crack tip singularity:
(17)ϕI(M)=rssin(θ2),
with s≥1. The sin-term ensures the discontinuity in the displacement field around the crack tip while rs guarantees that the crack closes at the crack tip. Note that it is standard to omit the crack tip enrichment for cohesive cracks. Omitting the crack tip enrichment tremendously simplifies the integration (only polynomial terms need to be integrated), the representation of the crack surface and the entire formulation. In that case, the crack closes at the element edge, see Figure 17, though more advanced techniques have been developed to close the crack inside an element [142–144] (without tip enrichment).(a) Topological enrichment versus (b) and (c) geometrical enrichment.
(a)
(b)
(c)XFEM has been mainly applied to problems involving crack growth of a single crack or a few cracks. There are only a few contributions dealing with fracture problems that involve the growth of numerous (several hundreds of) cracks [143, 144]. Methodologies have been developed in XFEM to account for complex crack patterns such as branching cracks and crack coalescence (within a single element) [145, 146]. However, the formulation becomes cumbersome with increasing number of cracks and crack branches. Moreover, those methods suffer (in practice) from the absence of reliable crack branching criteria.


#### 3.4.1. Representing the Crack Surface in XFEM

XFEM is a method that ensures crack path continuity. Enforcing crack path continuity is especially challenging in three dimensions. The crack path can be represented explicitly, usually with piecewise straight/planar crack segments [147, 148] or with other techniques such as level sets [149, 150] that trace the crack path in a more elegant way. In many publications, XFEM is closely related to level sets though any other technique can be used to describe the crack’s topology.Level sets are particularly attractive when a crack tip enrichment is used as geometry quantities can be calculated through the level sets. For example, the distance of a material point to the crack tip isr=f2+g2 and θ=arctan(g/f), see Figure 12; g is a level set orthogonal to f (i.e., ∇f·∇g=0) and is needed to uniquely identify the position of a material point with respect to the crack surface. Often, the level-set function is discretized
(18)fh(X)=∑I∈𝒮NI(X)fI;gh(X)=∑I∈𝒮NI(X)gI,
where commonly the shape functions NI(X) of the mechanical variables are used. Therefore, the enrichment can be expressed in terms of nodal values of both level set functions. Prabel et al. [151] used a (finer) finite difference scheme for the discretization of the level set function. However, the implementation effort is higher when different meshes are employed to discretize the level set and the mechanical variables. On the other hand, errors in the crack’s geometry will be introduced for low-order finite elements and curved cracks, Figure 14. As in most applications, a straight crack segment is attached in front of an existing crack front; a low-order finite elements should be sufficient for most crack propagation problems.Figure 14
Curved crack and its representation with level set functions when the level set is discretized with linear triangles.


#### 3.4.2. Integration in XFEM

Numerical integration in XFEM requires special attention, in particular, around the crack tip when employing non-polynomial enrichment functions. Laborde et al. [152] reported that several thousand Gauss-points are required to adequately reduce integration error in a “brute-force” approach. For step-enriched (or Heaviside-enriched) elements, a simple sub-triangulation is sufficient. Ventura [153] proposed an integration scheme for piece-wise straight cracks that avoid a sub-triangulation. Another alternative is to increase the number of integration points [151]. For singularities, Ventura et al. [154], Gracie et al. [155], Bordas et al. [156], Cheng and Fries [157] proposed efficient integration schemes. One of the most efficient schemes is the “almost polar integration” by Nagarajan and Mukherjee [158]. The almost polar integration is based on mapping a rectangle onto a triangle and avoids integration of the crack tip singularity. It was exploited for the first time in the context of XFEM by Laborde et al. [152]  and simultaneously by Béchet et al. [159].


#### 3.4.3. Blending Elements

Another major concern in XFEM are the so-called blending elements, Figure15. Blending elements are the elements adjacent to the enriched elements. It was shown by Chessa et al. [160], Stazi et al. [161], and Fries [162] that those blending elements result in the deterioration of the convergence rate in XFEM. When only a step enrichment is employed, a simple shifting can avoid problems due to blending. However, any non-constant enrichment function causes difficulties. Numerous approaches have been developed to eliminate drawbacks associated to blending [154, 160, 163]. A very simple and efficient approach was developed by Fries [162]. The idea is based on the modification of the enrichment function. Therefore, he enriched the nodes in the blending elements. He guaranteed that the enrichment in the enriched elements remains unchanged but vanishes in the “standard” elements.Figure 15
XFEM discretization with enriched domainΩenr, blending domain Ωblnd,  and the total (standard) domain Ωstd.


#### 3.4.4. Application to Fracture

Applying XFEM to linear problems in elastostatics, the final discrete equations for a single crack are obtained by substituting the test and trial functions into the weak form of the equilibrium equation yielding(19)[KIJuuKIJuaKIJKubKIJauKIJaaKIJKabKIJKbuKIJKbaKIJKbb]{uJaJbJK}={fIextfIextfIKext}
or
(20)KIJdJ=fIext,
with the stiffness matrix
(21)KIJ=[KIJuuKIJuaKIJKubKIJauKIJaaKIJKabKIJKbuKIJKbaKIJKbb],
where the superscripts u, a  and b indicate the “standard”, “step-enriched” and “tip-enriched” part of the discretization; dJ={uJaJbKJ}T is the vector with the nodal parameters, fIext={fIufIafKIb}T is the external force vector with  fKIb={f1Ibf2Ibf3Ibf4Ib} and
(22)fIu=∫ΩNI(X)bdΩ+∫ΓtNI(X)t-dΓ,fIa=∫ΩNI(X)ψIbdΩ+∫ΓtNI(X)ψIt-dΓ,fKIb=∫Ω‍NI(X)(ϕKI(X)-ϕKI(XI))bdΩ+∫ΓtNI(X)(ϕKI(X)-ϕKI(XI))t-dΓ
with body forces b and applied tractions t-. The stiffness matrix is
(23)KIJij=∫Ω‍(B¯Ii)TCB-JjdΩ,i,j=u,a,b,
where B¯Ii, i=u,a,b is the B operator defined by:
(24)BIu=∇NI(X),BIa=∇NI(X)ψI,BIb=∇[NI(X)(ϕKIl(X)-ϕKIl(XI)]).
For cohesive cracks without tip enrichment, the terms associated to the tip enrichment are omitted facilitating the implementation of the method.Although the XFEM approximation is capable of representing crack geometries that are independent of element boundaries, it relies on the interaction between the mesh and the crack geometry to determine the sets of enriched nodes [164]. This leads to particular crack configurations that cannot be accurately be represented by (13). Such cases are shown in Figure 16. As the crack size approaches the local nodal spacing, the set 𝒮c of nodes for the Heaviside or step enrichment is empty, Figure 16(b). Moreover, node 1 for the cracking case in Figure 16(a) or nodes 1 to 4 for the cracking case in Figure 16(b), respectively, contain two branch enrichments. Thus, the standard approximation gets difficulties with this crack configuration since the discontinuous function rsin(θ/2) extends too far. This problem always arises when one or more nodal supports contain the entire crack geometry. Usually, this kind of problem arises whenever a crack nucleates. Similar difficulties occur for approximations without any crack tip enrichment; see Figure 16(d). In order to close the crack within a single element, the set 𝒮c is empty as well. One solution is to refine the mesh locally such that the characteristic element size falls below that of the crack. An admissible crack configuration is shown in Figure 17.(a) and (b) Crack length that approach the local element size cannot be accurately represented by the standard XFEM approximation. Dots denote single enriched nodes and squares denote double (in our case, the node will contain the enrichment of two-crack tips) enriched nodes; (c) the dashed line shows the effective crack length; (d) even if no crack tip enrichment is used, in order to close the crack within a single element, no nodes have to be enriched with a step function.
(a)
(b)
(c)
effective crack length
(d)Figure 17
Admissible crack representation.More details on XFEM and its applications can be found in the excellent review papers by Karihaloo and Xiao [165], Fries and Belytschko [166], Belytschko et al. [167]  or the book from Mohammadi [168].


### 3.5. Phantom Node Method

An alternative to the standard XFEM for strong discontinuities was proposed by A. Hansbo and P. Hansbo [169]. They do not model the crack with additional degrees of freedom but by overlapping elements; see Figure 18. In 1D, it was shown by Song et al.  [170] that the crack kinematics of the A. Hansbo and P. Hansbo [169] version of XFEM can be derived from “standard” XFEM. However, the Hansbo version of XFEM has some advantages over standard XFEM.(i)
It avoids the “mixed” termsKua and Kau (or Mua and Mau in dynamics), (19), facilitating the implementation.(ii)
It leads to a better conditioned system matrix.(iii)
Standard mass lumping schemes such as the row-sum technique can be used. Note that, efficient mass lumping schemes have meanwhile been proposed for XFEM [171, 172] (and XEFG [173]).(iv)
Based on the full interpolation basis of the overlapping elements, the Hansbo-XFEM can straightforwardly integrate enhanced techniques that are more difficult to integrate in an enriched XFEM formulation [174].(v)
The implementation of the Hansbo XFEM is simpler. It has already been implemented in the commercial software package ABAQUS. In ABAQUS, it has to be used in combination with cohesive zone models. The ABAQUS implementation of the Hansbo-XFEM can handle numerous cracks (crack propagation) in statics. In contrast, XFEM with tip enrichment has also been incorporated in ABAQUS. However, a crack propagation algorithm for the tip-enriched XFEM is yet not available in ABAQUS. A plugin-in (Morfeo) that can handle three-dimensional crack growth in LEFM was recently developed by the company Cenaero.On the downside, we notice. (i)
It is difficult to apply the Hansbo XFEM to other problems besides cracks.(ii)
It is also difficult to incorporate a tip enrichment for the Hansbo-XFEM. A crack tip element has been developed for problems in 2D [175]. However, the extension of such a model to 3D is cumbersome.The Hansbo version of XFEM was subsequently developed in a static setting by Mergheim et al. [176], Areias et al. [148], Mergheim and Steinmann [177], and in dynamics by Song et al. [170] who called their approach phantom node method. The phantom node method is well suited for non-linear problems with cohesive cracks and problems with numerous cracks, in particular when the method needs to be implemented in an existing finite element code. It is less accurate for problems in LEFM due to the absence of a tip enrichment. Subsequently, the basic idea of the phantom node method is briefly summarized.Figure 18
The principle of the phantom node method in which the hatched area is integrated to build the discrete momentum equation; the solid circles represent real nodes and the empty ones phantom nodes.Consider a body that is cracked as shown in Figure18 and the corresponding finite element discretization. There are elements cut by the crack. To have a set of full interpolation bases, the part of the cracked elements which belongs in the real domainΩ0 is extended to the phantom domainΩp. Then, the displacement in the real domain Ω0 can be interpolated by using the degrees of freedom for the nodes in the phantom domain Ωp. The nodes are called the phantom nodes and marked by empty circles in Figure 18. The approximation of the displacement field is then given by [170]:
(25)uh(X,t)=∑I∈{𝒲0+,𝒲P-}uI(t)NI(X)H(f(X))+∑J∈{𝒲0-,𝒲P+}uJ(t)NJ(X)H(-f(X)),
where f(X) is the signed distance measured from the crack, 𝒲0+, 𝒲0-, 𝒲P+ and 𝒲P- are nodes belonging to Ω0+, Ω0-, ΩP+ and ΩP-, respectively. H(x) is the Heaviside function. As can be seen from Figure 18, cracked elements have both real nodes and phantom nodes. The jump in the displacement field is realized by simply integrating only over the area from the side of the real nodes up to the crack, that is, Ω0+ and Ω0-. We note that there are certain similarities to the visibility method used originally in meshfree methods [178–181].


### 3.6. Cohesive Segment Method and Cracking Node Method

Remmers et al. [182, 183], proposed the cohesive segment method that has similarities to XFEM. The crack is represented by a set of overlapping cohesive segments that cross three entire elements. In contrast to XFEM, no crack path continuity is required in the Cohesive Segment Method. A similar approach was developed by Song and Belytschko [184] in the Cracking Node Method. In the Cracking Node Method, the crack is required to directly pass through the node of an element. Also the cracking node method does not require crack path continuity.


### 3.7. Meshfree Methods

Thanks to the absence of a mesh, meshfree methods (MMs) offer another alternative to model fracture. They are particularly well suited for dynamic fracture and large deformations. Moreover, adaptive h-refinement can be easily implemented in a meshfree context [185–187]. Though all methods discussed so far are capable of capturing discrete cracks, a certain refinement around the crack front is still needed for accuracy and efficiency reasons [188]. The meshfree approximation is similar to the FEM formulation:
(26)uh(X)=∑J∈𝒮NJ(X)uJ,
where NJ(X) denotes the meshfree shape function and uJ are the nodal parameters. Most meshfree methods are not interpolatory complicating the imposition of essential boundary conditions. Popular meshfree methods include the smoothed particle hydrodynamics (SPH) method [189], the element-free galerkin (EFG) method [190], the reproducing kernel particle method (RKPM) [191], the point interpolation method (PIM) [192], the meshless local-petrov galerkin method (MLPG) [193] or the Finite Point Method [194], among others.Most meshfree shape functions depend on a weighting or kernel function, which is denoted byWJ(X). Usually, radial kernel functions are employed where WJ(X)=W(r0) is expressed in terms of r0=∥X-XJ∥. The kernel functions are chosen to have compact support, that is, WJ(X)=0  if  ∥X-XJ∥>hJ, where hJ is the “interpolation" radius often called dilation parameter. The kernel function can be expressed in terms of material coordinates X (Lagrangian kernel) or spatial coordinates x (Eulerian kernel). Formulations based on Lagrangian kernels are more efficient as the shape functions can be constructed at the beginning of the simulation while Eulerian kernels require updating during the course of the simulation. More details on meshfree methods can be found in the literature, for example, in several review articles and books [195–202]. Before discussing different possibilities to incorporate fracture into MMs, the shape functions of the EFG-method will be given as it is helpful for the explanations in the subsequent sections. They are constructed as follows:
(27)NJ(X)=pT(X)·A(X)-1·D(XJ),A(X)=∑Jp(XJ)pT(XJ)W(r-J;h),D(XJ)=∑Jp(XJ)W(r-J;hJ),
where p is the polynomial basis. Usually, a linear polynomial basis is choosen, given in 2D by
(28)p=[1,X,Y].


#### 3.7.1.  “Natural” Fracture

Due to the absence of a mesh, fracture in meshfree methods can occur naturally [38, 203, 204]. In contrast to element erosion, there is no need to erode particles. Fracture occurs when particles move outside the domain of influence of neighboring particles. Libersky et al. [203] presented impressive results of dynamic fracture and fragmentation predicting fragment distributions quite accurately. Dilts [205] presented an algorithm to detect fracture surfaces within the SPH method in 2D. He extended the approach later to three dimensions [206] and solved impressive impact problems. However, care has to be taken as fracture can occur artificially as shown for example, by Rabczuk et al. [207]. This is mostly the case when Eulerian kernels are employed. Formulations based on Lagrangian kernels do not suffer from numerical fracture. However, as the neighbor list does not change during the course of the simulation, fracture cannot be modeled “naturally” any more. Moreover, the Lagrangian kernel has to be consistently updated for problems involving large deformations. A very efficient algorithm to treat dynamic fracture and fragmentation in meshfree methods based on Lagrangian kernels was proposed by Maurel and Combescure [208] and Combescure et al. [209]. In their work, connectivities between neighboring particles were broken when their distance exceeds a given threshold. They solved complex problems including fluid-structure interaction and fragmentation [210, 211]. Another effective and simple method to treat fracture in meshfree methods is the material point method (MPM) [212]. The advantage of the MPM over most SPH methods for fracture is their computational efficiency [213]. While those methods are quite attractive to model dynamic fracture and fragmentation involving an enormous amount of cracks, they are not well suited for static fracture with few cracks due to their less accurate representation of the crack kinematics.


#### 3.7.2. Visibility Method

The visibility method is the first approach that introduces a discrete crack into the meshfree discretization. In the visibility method, the crack boundary is considered to be opaque. Thus, the displacement discontinuity is modeled by excluding the particles on the opposite side of the crack in the approximation of the displacement field, see Figure19:
(29)uh(X)=∑I∈𝒮+‍NI(X)uI.
This is identical to setting the shape function across the crack to zero as shown in Figure 19. The jump in the displacement is then computed by
(30)[[u(X)]]=(∑I∈𝒮+NI(X)uI)-(∑I∈𝒮-NI(X)uI).
Difficulties arise for particles close to the crack tip since undesired interior discontinuities occur (Figure 20) due to the abrupt cut of the shape function, see Figure 19. Nevertheless, Krysl and Belytschko [214] showed convergence for the visibility method. For linear complete EFG shape functions, they even showed that the convergence rate is not affected by the discontinuity. The length and size of the (undesired) discontinuities depend on the nodal spacing near the crack boundary. If the nodal spacing approaches zero, the length of the discontinuities tends to zero. With this argument and the theory of nonconforming finite elements, convergence of the discontinuous displacement field can be shown.Figure 19
Principle of the visibility, diffraction, and transparency methods with corresponding shape functions, from Belytschko et al. [200].(a) Undesired introduced discontinuities by the visibility method, (b) Difficulties with the visibility method for concave boundaries and kinks.
(a)
(b)It should also be noted that the visibility criterion leads to discontinuities in shape functions near nonconvex boundaries such as kinks, crack edges, and holes, as shown in Figure20 in two dimensions.An efficient implementation of the visibility method in 2D is given, for example, in [215].


### 3.8. The Diffraction Method

The diffraction method is an improvement of the visibility method. It eliminates theundesired interior discontinuities; see Figure 21 (see also Figure 19). The diffraction method is also suitable for non-convex crack boundaries. It is motivated by the way light diffracts around a sharp corner, but the equations used in constructing the domain of influence and the weight function bear almost no relationship to the equation of diffraction. The method is only applicable to radial basis kernel functions with a single parameter hI0.The diffraction method.
(a)
(b)The idea of the diffraction method is not only  to treat the crack as opaque but also to evaluate the length of the rayh0 by a path which passes around the corner of the discontinuity. A typical weight function is shown in Figure 19. It should be noted that the shape function of the diffraction method is quite complex with several areas of rapidly varying derivatives that complicates quadrature of the discrete Galerkin form [200]. Moreover, the extension of the diffraction method into three dimensions is complex. Nevertheless, several authors have presented implementations of the diffraction method in 3D [216, 217].


### 3.9. The Transparency Method

The transparency method was developed as an alternative to the diffraction method by Organ et al. [178]. The transparency method is easier extendable into three dimensions than the diffraction method. In the transparency method, the crack is made transparent near the crack tip. The degree of transparency is related to the distance from the crack tip to the point of intersection; see Figure 19.An additional requirement is usually imposed for particles close to the crack. Since the angle between the crack and the ray from the node to the crack tip is small, a sharp gradient in the weight function across the line ahead of the crack is introduced. In order to reduce this effect, Organ et al. [178] imposed that all nodes have a minimum distance from the crack surface, that is, the normal distance to the crack surface must be larger than  γh  with 0<γ<1; they suggested γ=1/4.


### 3.10. The “See Through” and “Continuous Line” Method

The “see-through” method was proposed by Terry [218] for constructing continuous approximations near nonconvex boundaries. Therefore, the boundary was considered as completely transparent such that the discontinuity is removed. Though the “see-through” method works well for capturing features such as interior holes, it is not well suited to model cracks.In the continuous line method from Krysl and Belytschko [214] and Duarte and Oden [219], the crack is completely transparent at the crack tip. In other words, particles with a partially cut domain of influence can see through the crack. This drastically shortens the crack. If no special techniques are introduced, the crack also does not close at the crack tip that leads to inaccurate solutions. Crack closure at the crack tip can be enforced with Lagrange multiplier or by decreasing the domain of influence of nodes close to the crack tip.Belytschko and Fleming [220] suggest to combine different methods depending on the convexity of the crack boundary for example, to use the visibility for convex boundaries, and other methods for non-convex crack boundaries. They suggest a criterion based on the angle of the wedge that can be written in terms of the surface normal; see Figure 22. When  nA·nB≤β with β=0° as a cutoff value, the boundary can be considered as convex, otherwise non-convex. The cutoff value of β=0° corresponds to the wedge angle of  ω=90° in Figure 22.Figure 22
Domain of influence near a wedge-shaped non-convex boundary. The boundary is enforced ifnA·nB≤β.


### 3.11. Extended Meshfree Methods

Enrichment in meshfree methods was used before the development of XFEM. A good overview paper on enriched meshfree methods in the context of LEFM is given by Fleming et al. [221], see also the review articles by Belytschko et al. [200], Nguyen et al. [201], and Huerta et al. [202]. The enrichment schemes are classified into intrinsic enrichment and extrinsic enrichment. An intrinsic enrichment introduces enrichment functions as additional functions to the polynomial basis compare to (28) as the following:
(31)p=[1,X,Y,rsin(θ2),rcos(θ2),rsin(θ2),sinθ,rcos(θ2),sinθ].
This concept was later also used in the intrinsic XFEM [222]. The major drawback of this approach is its higher computational cost as the enrichment needs to be employed everywhere to avoid artificial discontinuities. Techniques and blending enriched approximations to “standard” approximations have been used but they introduce additional complexity. Duflot and Nguyen-Dang [223] proposed an alternative intrinsic enrichment by enriching the kernel functions. Their approach does not require blending techniques and is computationally cheaper. However, it is also less accurate. The extrinsic enrichment can be categorized into extrinsic MLS enrichment and extrinsic PU enrichment. The extrinsic MLS enrichment introduces enrichment parameters that can be interpreted as stress-intensity factors (SIFs), at least when the enrichment is introduced globally. Hence, the SIFs can be obtained directly from the numerical analysis without computing the J-integral or interaction integral. However, it was shown by Zamani et al. [224, 225] that such approaches are less accurate. The extrinsic PU enrichment is comparable to an XFEM enrichment. It is believed that it is the best choice in terms of accuracy, efficiency, and flexibility.An extended element-free Galerkin (XEFG) method based on vector level sets was first proposed by Ventura et al. [226] in the context of LEFM. It was later extended to nonlinear materials and cohesive cracks by Rabczuk and Zi [227], Rabczuk et al. [228]. In contrast to XFEM, due to the heavily overlapping shape functions, a crack tip enrichment has to be employed to ensure crack closure at the crack tip though alternative methods have been developed that avoid the use of tip enrichments [229, 230]. However, those approaches are cumbersome in three dimensions.


### 3.12. Cracking Particles Method

While all previous meshfree approaches for fracture require crack path continuity, the cracking particles method [231–235] does not represent the crack as continuous surface. This makes the method particularly useful for 3D applications including complex crack patterns such as crack branching and crack coalescence. In contrast to methods that enforce crack path continuity, crack branching and crack coalescence are a (more) natural outcome of the analysis in the cracking particles method.To model cracks in the cracking particles method, the displacement is decomposed into continuous and discontinuous parts:(32)u(X)=ucont(X)+uenr(X),
where X are the material coordinates, ucont denotes the continuous displacement, and uenr denotes the discontinuous part.The crack is modelled by a set of discrete cracks as shown in Figure23. These discrete cracks are restricted to lie on the particles, that is, each crack plane (or line in 2D) always passes through a particle. Since the crack geometry is described by the set of cracked particles, a representation for the geometry of the crack is not needed.Schematic on the right shows a crack model for the crack on the left.
(a)
(b)The approximation of the displacement field in the cracking particles method is given by(33)uh(X)=∑I∈𝒩NI^(X)uI+∑I∈𝒩cNI(X)s(fI(X))qI,
where 𝒩 is the total set of nodes in the model and 𝒩c  is the set of cracked nodes;  s(fI(X)) is the step function. In general, different shape functions can be used for the continuous and discontinuous part of the displacement field.The jump in the displacement across the crack depends only on the discontinuous part of the displacement fielduenr(X) and hence the enrichment parameters qI. It can be shown from (33) that
(34)⟦u(X)⟧=2∑I∈𝒩cNI(X)qI.
The factor 2 comes from the fact that the step function is used that is −1 on one side of the crack and +1 on the other side of the crack.The cracking particles method might suffer from spurious crack patterns [236, 237]. However, they can be eliminated by a set of simple crack injection rules [184]. Firstly, one has to distinguish between crack propagation and crack initiation. Since the crack is not considered as continuous surface, crack propagation is assumed when no cracked node is detected in the vicinity (in a circle [sphere in 3D] with radius r=αh, α≥1.5) of a sampling point, where cracking is detected. For propagating cracks, spurious cracks as shown in Figure 24 are avoided by introducing a prevention zone, Figure 25. For example, for particle L, the crack prevention zone is given through the angle α:
(35)|n·x-xL∥x-xL∥|<cos(Π2-α).
For branching cracks, it is often observed that the orientation of neighboring cracks severely differ, see node K in Figure 25. Therefore, cracks are also allowed in the prevention zone when the following criterion is fulfilled:
(36)nL·nK<cosβ,
with α=2Π/9 and β=Π/3. Similar criteria are used within the cracking particles method in H. X. Wang and S. X. Wang [238]. Zhang [237] reported that cracking rules can be avoided by the use of rotating crack segments.(a) Spurious cracking and (b) improved crack pattern.
(a)
(b)Figure 25
Crack prevention: nodesJ and K are located in the crack prevention zone; node I is located in the crack propagation zone.Recently, the cracking particles method has been presented without enrichment to model fracture in continua [235] and structures [239]. This results in less degrees of freedom and facilitates diagonalizing the mass matrix when the method is used in explicit dynamics.


### 3.13. Boundary Element Method

In the Boundary Element Method (BEM), the weak form is formulated in boundary integral form. It reduces therefore the dimension of the problem, for example, a three-dimensional problem reduces to two dimensions. However, the BEM is only applicable to problems where Green’s functions can be computed restricting the application range of the method. Due to the dimension-reduction, the remeshing procedure is fairly simple as the crack is part of the boundary. Interesting applications of the BEM to fracture can be found, for example, in [240–257].


### 3.14. Isogeometric Analysis for Fracture

Isogeometric analysis (IGA) [258] was developed to unify CAD and CAE. Verhoosel et al. [130] were the first who incorporated fracture in isogeometric analysis based on T-splines. Therefore, they used the concept of knot insertion to introduce discontinuities in the displacement approximation. This method suffers from similar problems as interelement-separation methods. Moreover, for complex crack patterns with slightly deviating angles between the crack branches, it is expected that the mesh will be extremely distorted. De Luycker et al. [259], Ghorashi et al. [260] and recently Tambat and Subbarayan [261] combined XFEM with isogeometric analysis (XIGA) allowing crack growth without remeshing also in the context of IGA.


### 3.15. The Variational Approach to Fracture

Variational methods in fracture mechanics are a relatively new development in the field. The underlying theory was laid down by Francfort and Marigo [262], who proposed the idea that cracks should propagate along a path of least energy. The goal of variational methods is to circumvent several important weaknesses of the classical theory emanating from the work of Griffith in 1920.In classical fracture mechanics, crack propagation is assumed to occur under thermodynamic equilibrium. In particular for crack formation occurring under constant load, the external work done by the applied loading is equal to twice the energy in the bulk, so that the total energy of the system may be expressed asU=-WL+UE+US=-UE+US, where US is a postulated quantity known as the surface (dissipation) energy of the crack system. Thermodynamic equilibrium with respect to crack growth then requires that dU/da=0, where a is the crack length. The previous expression provides the critical condition for fracture. The simplest case to consider is an infinite solid with a straight crack of length 2a and subjected to far-end tractions of magnitude σ directed perpendicular to the crack. For this problem, the exact solutions for the stresses and strains are available due to Inglis, so that an analytic form of UE may be obtained, equal to πa2σL2/E for the case of plane stress. We now assume that the surface energy to be US=4γa,  where γ is the surface energy per unit area, which is considered a material property. Enforcement of thermodynamic equilibrium yields the expression for the critical applied stress required for crack propagation: σfa=2Eγ/π. According to the classical theory, the existing crack will not propagate for σ<σf and will do so for σ=σf.Despite improvements made by subsequent authors to account for nonlinearity and inelasticity, the original form of the fracture criterion remained, that is,σfa=constant. As pointed out in [262], this inverse proportionality of the critical stress to the square root of the initial crack length means that classical fracture mechanics is unable to predict crack initiation, since for a body without an initial crack (a=0), σf becomes infinite. Other important weaknesses of the classical theory listed in [262] are its inadequacy for predicting the direction of the crack path, and its inability to handle crack jumps (brutal cracking). The former arises from the fact that Griffith’s criterion gives only one constraint, whereas the description of the crack tip propagation requires a number of functions equal to the dimensionality of the problem. The latter weakness stems from the classical theory being silent on what happens when σ>σf (unstable cracking). Griffith’s criterion is also incapable of predicting when a crack will branch, so that additional criteria are required in order to evaluate the possibility of this kind of response.In order to address the above mentioned weaknesses of classical fracture mechanics, Francfort and Marigo [262] proposed an alternative formulation of brittle fracture based on energy minimization concepts. For a body denoted by Ω- in D-dimensional space, energy (similar to the surface energy concept of Griffith) is assigned to the crack system Γ∈Ω- and expressed as
(37)Es(Γ)=∫Γk(x)dℋN-1(x),
where ℋN-1 denotes the (N-1)-D Hausdorff measure. The function k(x) represents the energy required to create an infinitesimal crack at the point x  and is strictly positive and bounded away from zero. The above definition is rather accommodating in terms of crack geometries, and in addition to the usual “sharp” crack surfaces of classical fracture mechanics, it also admits entities like point clouds and edge cracks. However, cracks having N-D Hausdorff measure (termed as “fat” cracks in [262]) are not allowed, for example a crack with computable area in a 2D or volume in 3D, since these will require an infinite amount of energy to propagate. The bulk energy is defined as
(38)Ed(Γ,U)=infv∈𝒞(Γ,U)∫Ω∖ΓW(x,ϵ(v,(x)))dx,
in which 𝒞(Γ,U) is the set of kinematically admissible displacement fields. The total energy of the body is the sum of these two terms:
(39)E(Γ,U)=Ed(Γ,U)+Es(Γ).
The problem then is to find the mapping Γ(t) such that at some given time t>0, the crack geometry Γ(t) is the one that minimizes E(U,Γ(t)) and furthermore contains all Γ(s) for s<t (condition of irreversibility). The authors warn that the use of global energy minimization to drive the crack evolution is not based on thermodynamic principles, but is adapted as “… a rather convenient postulate which provides useful insight into a variety of behaviors” [262]. They show also that this new formulation is free of the earlier identified weaknesses of the classical theory.Numerical experiments on the proposed formulation were carried out by Bourdin et al. [263], who pointed to the similarity of the new formulation to models of image segmentation obtained by minimization of the Mumford-Shah functional. This meant that the inapplicability of standard numerical methods, owing to the fact that the formulation allowed for jump sets of the displacement field (representing the cracks) whose locus was unknown a priori. Two methods of solution were demonstrated, based on the mathematical concept of Γ-convergence. In one of these methods, a second field was included in the energy functional in addition to the primary displacement field. This auxiliary scalar field acted as a regularizer for the jump sets of the displacement field, and took on values of 0 on the crack and 1 away from it. The regularizing variable was referred to by later researchers as the phase field, and it allowed for the treatment of the global energy minimization as a standard variational problem for which classical FEM is up to the task, albeit with the restriction that the characteristic length of the mesh should tend to zero faster than the characteristic length of the regularization so as not to overestimate the surface energy of the crack. Due to the inclusion of a second field, a coupled system of equations must now be solved consisting of the original equilibrium/linear momentum equations and the evolution PDE of the phase field.A phase field model for mode III dynamic fracture was devised by Karma et al. [264], which made use of a phase field evolution equation based on the standard two-minimum Ginzburg-Landau form, that is,
(40)τ∂tϕ(x,t)=Dϕ∇2ϕ-ddϕ(14ϕ2(1-ϕ2))-μ2g′(ϕ)(ϵ:ϵ-ϵc2),
where g is a function obeying the constraints g(0)=0, g(1)=1, and g′(0)=g′(1)=0. The KKL phase field model was subsequently utilized by Hakim and Karma [265] to analyze the laws of quasistatic crack tip motion. Among other things, they showed that for kinked cracks in anisotropic media, force-balance gives predictions that are significantly different from those obtained using the principle of maximum energy release rate, and that the predictions obtained from force-balance hold even when fracture is modeled as irreversible. It should be mentioned that up to this point, phase field implementations featured isotropic damage wherein fracture occurred both in tension and compression. This sometimes resulted in unphysical response, such as sample interpenetration in the crack branching simulation in Bourdin et al. [263].Thermodynamically consistent phase field models of fracture were developed by Miehe et al. [266, 267], along with corresponding incremental variational principles and numerical implementation within a finite element framework. Their implementation of the phase field was also slightly different from that of previous authors whereas Bourdin et al. [263] and subsequent papers utilized the convention that the phase field took values of 0 at cracks and 1 at unbroken states (hence a pseudo material density function) and Miehe et al. reversed the convention by assigning to the phase field values of 1 and 0 for fully cracked and fully unbroken states respectively. As a consequence, the functional associated with the phase field may now be seen to represent the crack surface itself. They then derived the evolution equation for the phase field based on the assumption that the solution is negative exponential in nature. Furthermore, they extended the application of phase fields to fracture simulations involving viscous overforce response via a time-regularized three-field formulation, and showed that the coupled system (linear momentum + phase field evolution) could be solved either monolithically or via a robust staggered-update scheme [266, 267].An important addition of Miehe et al. [266, 267] to the existing theory was the modelling of anisotropic degradation by using an additive decomposition of the stored energy of the undamaged solid, to come up with an anisotropic energy function of the form
(41)Ψ(ϵ,ϕ)=[g(ϕ)+k]Ψ0+(ϵ)+Ψ0-(ϵ),
where the positive and negative parts of the energies are defined by
(42)Ψ0±(ϵ)=λ2〈ϵ1+ϵ2+ϵ3〉±2+μ(〈ϵ1〉±2+〈ϵ2〉±2+〈ϵ3〉±2),
and 〈x〉± are ramp functions defined as 〈x〉±=(1/2)[|x|±x]. The energy decomposition allows for the case where fracture occurs in tension only as opposed to both tension and compression, and effectively fixes the sample inter-penetration issue encountered by Bourdin in crack branching simulations. Following the reversed convention on the phase field, the constraints on g(ϕ) are g(0)=1, g(1)=0, and g′(1)=0. To maintain numerical stability, a small positive constant k is included in the formulation so that the material retains some residual stiffness even when it attains a fully damaged state.One challenge with the use of phase field models for fracture is the computational expense associated with mesh size requirements, since the use of a mesh having a characteristic length that is not small enough compared to the crack regularization parameter yields erroneous results with regard to the energy. Based on numerical experiments utilizing the stationary phase field equation, Miehe et al. [266, 267] found that the characteristic length of the mesh should be no more than half the value of the regularization parameter, at least presumably in the vicinity of the regularized crack. In an effort to reduce computational expense, Kuhn and Müller [268, 269] (see also their work in [270]) introduced specially engineered FEM shape functions of an exponential nature to discretize the phase field. These shape functions have the form
(43)N1e(ξ,δ)=1-exp(-δ(1+ξ)/4)-1exp(δ/2)-1,N2e(ξ,δ)=exp(-δ(1+ξ)/4)-1exp(δ/2)-1,
in 1D for the 2-node line element. 2D shape functions for the 4-node quadrilateral are obtained as a tensor products of the above, with certain tweaks. The shape functions are parametrized by δ which is the ratio between the element size and phase field regularization parameter. Their results showed that with the special element shape functions, accurate prediction of surface energy associated with the phase field is possible with a much lower level of refinement compared to standard FEM shape functions. The drawback to their approach is that some information regarding crack orientation is necessary for the proper construction of the exponential shape functions, so that numerical examples were confined to cases, where the general direction of the crack path is known a priori.


## 3.1. Element Erosion

The easiest way to deal with discrete fracture is the element erosion algorithm [66–69]. Element erosion algorithms do not require any representation of the crack’s topology. Most approaches do not really delete the elements but set the stresses in the “deleted” elements to zero. Approaches that indeed remove elements often replace the deleted elements with rigid masses. In order to account for the correct energy dissipation in the post-localization domain, concepts from smeared crack models are applied. Song et al. [70] compared different computational methods for certain benchmark problems in dynamic fracture. They report extreme mesh sensitivity for element erosion algorithms and conclude they are not well suited for dynamic fracture. However, due to their simplicity, they are incorporated in commercial software packages such as LS DYNA.For brittle fracture, Pandolfi and Ortiz [71] derived an interesting eigenerosion approach. Their idea is based on the eigenfracture approach by Schmidt et al. [72], where the eigendeformation is restricted in a binary sense. The element can be either intact or fractured. In the first case, the element behavior is elastic. Otherwise, the element is eroded. They showed that in contrast to many other element erosion approaches [66, 67, 73, 74], the eigenerosion approach converges to Griffith fracture with increasing mesh refinement.


## 3.2. Interelement-Separation Methods

Standard finite elements have difficulties to capture the crack kinematics since they use continuous trial functions that are not particularly well adapted for solutions with discontinuous displacement fields. One of the first models capable of modelling cracks within the FEM are the so-calledinterelement-separation methods [46, 47, 75–80]. Interelement seperation methods were mostly applied to dynamic fracture with cohesive cracks. In the approaches by Xu and Needlemen [75, 76], cohesive surfaces were introduced at the beginning of the simulation. In contrast, Camacho and Ortiz [46] adaptively introduced cohesive surfaces at element edges when a certain cracking criteria is met. In interelement-separation methods, cracks are only allowed to develop along existing interelement edges. This endows the method with comparative simplicity, but can result in an overestimate of the fracture energy when the actual crack paths are not coincident with element edges. The results depend not only on the mesh size (and form of the chosen element), but also on the mesh bias that can be compensated only by remeshing that is particularly cumbersome in 3D. It has been noted that the solutions sometimes depend on mesh refinement; see Falk et al. [81]. This sensitivity has been mollified by adding randomness to the strength, as in Zhou and Molinairi [80] and Espinosa et al. [82]. For dynamic problems, interelement-separation methods are often used with rigid cohesive zone models. Special care has to be taken in order to guarantee time continuity. Time continuity requires temporal convergence when the time step tends to zero. Violating time continuity might lead to unstable or oscillatory responses, and the results severely depend on the chosen time step [56, 83].Many interesting problems have been studied by interelement-separation methods [47, 75, 79, 82–84]. They are most popular for dynamic fracture and problems involving numerous cracks with complex fracture patterns.Special crack tip elements for linear fracture problems were developed in the early 70’s [85–88]. They account for the crack tip singularity and were applied to stationary cracks. Recently, such approaches were incorporated in the smoothed finite element method by Liu et al. [89, 90]. These approaches were extended to anisotropic materials [91], interface cracks [92], and plastic fracture mechanics [93]. In [94, 95], these approaches were applied to crack propagation problems in two dimensions. Adaptive remeshing techniques were employed. As the crack is aligned to the crack boundary, these methods can also be classified as interelement-separation methods. One drawback is that only the r-term is accounted for in the crack tip element but not the discontinuity behind the crack tip.


## 3.3. Embedded Elements (EFEM)

In 1987, Ortiz et al. [77] modified the approximation of the strain field to capture weak discontinuities in finite elements to improve the resolution of shear bands. Therefore, they enriched the strain field to obtain the kinematic relations shown in Figure 9(a). Based on this idea, Belytschko et al. [96]  allowed two parallel weak discontinuity lines in a single element, Figure 9(b), so that the element was able to contain a band of localized strain. Dvorkin et al. [97] were the first who developed a method that was able to deal successfully with strong discontinuities in finite elements, the class of embedded elements (EFEM). The name comes from the fact that the localization zone is embedded in a single element; see Figure 9(c). This way, crack growth can be modeled without remeshing. This class of methods are much more flexible than schemes that allow discontinuities only at element interfaces, and it eliminates the need for continuous remeshing. Several different versions of the EFEM were developed. All techniques had in common that the approximation of the displacement field was enriched with additional parameters to capture the jump in the displacement field in a single element. These parameters, in the following called enrichment, were inherent to the cracked element and the jump in the displacement field depends only on this enrichment. The different methods were derived from various principles, starting from the extended principle of virtual work over Hellinger Reissner to Hu-Washizu variational principle, using an enhanced assumed strain (EAS) or B-bar format. One major advantage of the EFEM (compared to XFEM discussed in the next section) is that the additional unknowns could be condensed on the element level, so that discontinuities could be captured with very small changes of the existing code.Element with (a) one weak discontinuity, (b) two weak discontinuities, and (c) one strong discontinuity.
(a)
(b)
(c)The first version of the EFEM is often called statical optimal symmetric (SOS) since traction continuity is fulfilled, but it is not possible to capture the correct crack kinematics. The correct relative rigid body motion of the element is not guaranteed and it has been shown by Jirásek [98] that SOS formulations lead to stress locking. The kinematical optimal symmetric (KOS) version by Lotfi and Shing [99] ensures the correct crack kinematics but violates the traction continuity condition. Consequently, the criteria for the onset of localization written in terms of stresses in the bulk are no longer equivalent to the same criteria written in terms of the tractions on the discontinuity area; for example, for the Rankine criterion, the normal traction at the onset of localization should be equal to the tensile strength and the shear traction should be zero. This cannot be properly reproduced by the KOS formulation. The kinematical and statical optimal nonsymmetric (KSON) version of embedded elements (see Dvorkin et al. [97], Klisinski et al. [100]) guarantees traction continuity and the appropriate crack kinematics but leads to a non-symmetric stiffness matrix with all its disadvantages with respect to solving the linearized system of equations.The approximation of the displacement field in the EFEM is given by [101]
(3)uh(X)=∑I∈𝒮NI(X)u~J+ℳs(e)(X)⟦u^I(e)(X)⟧,
where u~J are the nodal parameters of the usual approximation, u^  is the enrichment, and the superimposed “(e)” indicates the enrichment on element level. The function ℳs(e) is given by
(4)ℳs(e)(X)={0∀(e)∉SHs(e)-ρ(e)∀(e)∈Sρ(e)=∑I=1Ne+NI+(X),
where Hs is the step function acting on the crack line S and Ne+ is the number of nodes of element (e) that belong to the domain Ω0+;  see Figure 9(c). The discretization of the strain field can be written as
(5)ϵh(X)=∑I∈𝒮‍(∇0NI(X)⊗u~I)S-(∇ρ(e)⊗⟦u^I(e)(X)⟧)S+ηs(e)k(⟦u^I(e)(X)⟧⊗n)S,
where the superimposed S denotes the symmetric part of a tensor, ηs(e)/k is a regularized Dirac delta function, and ηs(e) is a collocation function defined as
(6)ηs(e)={1∀X∈Sek0∀X∉Sek,
with the thickness k of the localization band. Considering the equilibrium equation in elastostatics and with the trial and test function of the structure of (3), the discrete equations in matrix form can be written as
(7)[Ku~u~(e)Ku~u^(e)Ku^u~(e)Ku^u^(e)]{u~(e)⟦u^I(e)⟧}={FIext0},
with
(8)Ku~u~(e)=∫Ω0BTCBdΩ0,Ku~u^(e)=∫Ω0BTCB^dΩ0,Ku^u~(e)=∫Ω0B^*TCBdΩ0,Ku^u^(e)=∫Ω0B^*TCB^dΩ0,
where C is the elasticity tensor, B is the B operator containing the derivatives of the shape functions and B^ is the B operator of the enrichment that depends on the embedded element formulation (SOS and KOS or KSON). For the KSON formulation, B^*≠B^, that will result in a non-symmetric stiffness matrix. For the SOS and KOS formulation, B^*=B^. The elemental enrichment ⟦u^I(e)⟧ can be condensed on the element level:
(9)⟦u^I(e)⟧=-[Ku^u~(e)]-1Ku^u^(e)u~(e),
that leads with the expression for ⟦u^I(e)⟧ to the system of equations
(10)K¯u=f,
with
(11)K¯=Ku~u~-Ku~u^Ku^u^-1Ku^u~.
The crack in embedded elements is modelled by a set of crack segments with piecewise constant crack opening, see Figure 10 though recent versions of the EFEM can handle linear crack opening [102]. The EFEM has been applied to many interesting problems [103–106] including piezoelectric materials [107], geomaterials [108], coupled formulations [109], finite strains [110], crack branching [111], and dynamic fracture [112], to name a few.(a) Piecewise constant crack opening in embedded elements and (b) linear crack opening for XFEM.
(a)
(b)Oliver et al. [113] compared embedded elements with the extended finite element method (that will be discussed in the next section) (XFEM) for several selected examples and showed that embedded elements can be as accurate as XFEM. However, examples in LEFM that could be compared to analytical solutions are missing.Methods that model the crack as continuous surface and therefore contain a “physical” crack tip are expected to be more accurate. However, the requirement to represent the crack surface is a drawback of those methods, in particular for complex crack patterns (e.g., in three-dimensions or for branching cracks). EFEM in principle does not require crack path continuity [114, 115]. However, some authors report that the elements then become sensitive to the orientation of the discontinuity surface with respect to the nodes, and in some unfavorable situations the response of the element ceases to be unique [116].


## 3.4. Extended Finite Element Method (XFEM)

A very flexible method that can handle linear and nonlinear crack openings is the extended finite element method (XFEM) developed by Belytschko et al. [11, 117, 118]. This method employs the local partition of unity concept [119–121] and introduces additional nodal parameters for the elements cut by the crack. Hence, the additional unknowns cannot be condensed on the element level as in EFEM. The displacement discontinuity depends only on the additional nodal parameters. The principle idea of XFEM is shown in Figure 11, simplified for a crack in one dimension. Though XFEM has originally been developed for fracture, it has been extended to numerous applications including two-phase flow [122, 123], fluid-structure interaction [124, 125], biomechanics [126], inverse problems [127, 128], multifield problems [129–132], among others. XFEM has been incorporated meanwhile into commercial software (XFEM) and has become one of the most popular methods for fracture. The generalized finite element method (GFEM) [120, 121, 133–139] is similar to XFEM.Figure 11
Principle of XFEM for strong discontinuities in 1D.The basic idea of XFEM is to decompose the displacement field into a continuous partucont and a discontinuous part udisc:
(12)uh(X)=ucont(X)+udisc(X).
The continuous part is the standard FE interpolation and additional information is introduced into the FE interpolation through the local partition of unity approach by adding an enrichment (udisc). Therefore, additional degrees of freedom are introduced for the discretization of the test- and trial functions. Those additional degrees of freedom are multiplied with the so-called enriched shape functions. They are the product of the “standard” element shape functions and enrichment functions that contain information about the solution. Note that the solution does not need to be known exactly. For example, a crack results in a discontinuous displacement field. Therefore, a discontinuous function is used as enrichment function. In linear elastic fracture mechanics (LEFM), the analytical near crack tip solution is known as well and can be included in the XFEM approximation. The approximation of the displacement field for nc cracks with mt crack tips reads
(13)uh(X)=∑I∈𝒮NI(X)uI+∑N=1nc∑I∈𝒮cNI(X)ψI(N)aI(N)+∑M=1mt∑I∈𝒮tNI(X)∑K=1NKϕKI(M)bKI(M),
where 𝒮 is the set of nodes in the entire discretization, 𝒮t is the set of nodes around the crack tip, and 𝒮c is the set of nodes associated to elements completely cut by the crack; ψI(N) is the enrichment function of crack N; ϕKI(M) is the enrichment function for the crack tip M; aI and bKI are additional degrees of freedom that needs to be solved for.For elements completely cut by the crack, the jump function is often chosen to cause a discontinuous displacement field(14)ψI(N)=sign[f(N)(X)]-sign[f(N)(XI)],
with
(15)f(N)(X)=sign[n·(X-X(N))]min(X-X(N))︸,XN∈Γc(N)
where n denotes the normal vector to the crack surface and f(N)(X) is a level set function; see Figure 12. The zero-isobar f=0 describes the position of the crack surface. Often a shifting is employed, see the second term on the RHS of (14) and Figure 11, and ensures that the enriched shape function NI(x)ψI(N) vanishes in the neighboring element. The shifting also maintains the interpolatory character of the XFEM approximation or in other words: the nodal parameters uI remain the physical displacements.Figure 12
Representation of the crack surface by level setsf and g.In LEFM, the crack tip enrichment is choosen according to the analytical solution:(16)ϕKI(M)=[rsin(θ2),rcos(θ2),rsin(θ2),sinθ,rcos(θ2),sinθ],
where r and θ are explained in Figure 12. Early approaches used a so-called topological enrichment, Figure 13(a). However, a topological enrichment leads to a deterioration of the convergence rate as the domain with the accurate tip enrichment shrinks when the mesh is refined. Therefore, it is important to use a geometrical enrichment, Figures 13(b) and 13(c). Note that a geometrical enrichment can drastically increase the condition number of the system matrix and requires the use of pre-conditioners [140] or other techniques [121, 141]. For non-linear materials and cohesive cracks, a tip enrichment should be used that avoids the crack tip singularity:
(17)ϕI(M)=rssin(θ2),
with s≥1. The sin-term ensures the discontinuity in the displacement field around the crack tip while rs guarantees that the crack closes at the crack tip. Note that it is standard to omit the crack tip enrichment for cohesive cracks. Omitting the crack tip enrichment tremendously simplifies the integration (only polynomial terms need to be integrated), the representation of the crack surface and the entire formulation. In that case, the crack closes at the element edge, see Figure 17, though more advanced techniques have been developed to close the crack inside an element [142–144] (without tip enrichment).(a) Topological enrichment versus (b) and (c) geometrical enrichment.
(a)
(b)
(c)XFEM has been mainly applied to problems involving crack growth of a single crack or a few cracks. There are only a few contributions dealing with fracture problems that involve the growth of numerous (several hundreds of) cracks [143, 144]. Methodologies have been developed in XFEM to account for complex crack patterns such as branching cracks and crack coalescence (within a single element) [145, 146]. However, the formulation becomes cumbersome with increasing number of cracks and crack branches. Moreover, those methods suffer (in practice) from the absence of reliable crack branching criteria.


### 3.4.1. Representing the Crack Surface in XFEM

XFEM is a method that ensures crack path continuity. Enforcing crack path continuity is especially challenging in three dimensions. The crack path can be represented explicitly, usually with piecewise straight/planar crack segments [147, 148] or with other techniques such as level sets [149, 150] that trace the crack path in a more elegant way. In many publications, XFEM is closely related to level sets though any other technique can be used to describe the crack’s topology.Level sets are particularly attractive when a crack tip enrichment is used as geometry quantities can be calculated through the level sets. For example, the distance of a material point to the crack tip isr=f2+g2 and θ=arctan(g/f), see Figure 12; g is a level set orthogonal to f (i.e., ∇f·∇g=0) and is needed to uniquely identify the position of a material point with respect to the crack surface. Often, the level-set function is discretized
(18)fh(X)=∑I∈𝒮NI(X)fI;gh(X)=∑I∈𝒮NI(X)gI,
where commonly the shape functions NI(X) of the mechanical variables are used. Therefore, the enrichment can be expressed in terms of nodal values of both level set functions. Prabel et al. [151] used a (finer) finite difference scheme for the discretization of the level set function. However, the implementation effort is higher when different meshes are employed to discretize the level set and the mechanical variables. On the other hand, errors in the crack’s geometry will be introduced for low-order finite elements and curved cracks, Figure 14. As in most applications, a straight crack segment is attached in front of an existing crack front; a low-order finite elements should be sufficient for most crack propagation problems.Figure 14
Curved crack and its representation with level set functions when the level set is discretized with linear triangles.


### 3.4.2. Integration in XFEM

Numerical integration in XFEM requires special attention, in particular, around the crack tip when employing non-polynomial enrichment functions. Laborde et al. [152] reported that several thousand Gauss-points are required to adequately reduce integration error in a “brute-force” approach. For step-enriched (or Heaviside-enriched) elements, a simple sub-triangulation is sufficient. Ventura [153] proposed an integration scheme for piece-wise straight cracks that avoid a sub-triangulation. Another alternative is to increase the number of integration points [151]. For singularities, Ventura et al. [154], Gracie et al. [155], Bordas et al. [156], Cheng and Fries [157] proposed efficient integration schemes. One of the most efficient schemes is the “almost polar integration” by Nagarajan and Mukherjee [158]. The almost polar integration is based on mapping a rectangle onto a triangle and avoids integration of the crack tip singularity. It was exploited for the first time in the context of XFEM by Laborde et al. [152]  and simultaneously by Béchet et al. [159].


### 3.4.3. Blending Elements

Another major concern in XFEM are the so-called blending elements, Figure15. Blending elements are the elements adjacent to the enriched elements. It was shown by Chessa et al. [160], Stazi et al. [161], and Fries [162] that those blending elements result in the deterioration of the convergence rate in XFEM. When only a step enrichment is employed, a simple shifting can avoid problems due to blending. However, any non-constant enrichment function causes difficulties. Numerous approaches have been developed to eliminate drawbacks associated to blending [154, 160, 163]. A very simple and efficient approach was developed by Fries [162]. The idea is based on the modification of the enrichment function. Therefore, he enriched the nodes in the blending elements. He guaranteed that the enrichment in the enriched elements remains unchanged but vanishes in the “standard” elements.Figure 15
XFEM discretization with enriched domainΩenr, blending domain Ωblnd,  and the total (standard) domain Ωstd.


### 3.4.4. Application to Fracture

Applying XFEM to linear problems in elastostatics, the final discrete equations for a single crack are obtained by substituting the test and trial functions into the weak form of the equilibrium equation yielding(19)[KIJuuKIJuaKIJKubKIJauKIJaaKIJKabKIJKbuKIJKbaKIJKbb]{uJaJbJK}={fIextfIextfIKext}
or
(20)KIJdJ=fIext,
with the stiffness matrix
(21)KIJ=[KIJuuKIJuaKIJKubKIJauKIJaaKIJKabKIJKbuKIJKbaKIJKbb],
where the superscripts u, a  and b indicate the “standard”, “step-enriched” and “tip-enriched” part of the discretization; dJ={uJaJbKJ}T is the vector with the nodal parameters, fIext={fIufIafKIb}T is the external force vector with  fKIb={f1Ibf2Ibf3Ibf4Ib} and
(22)fIu=∫ΩNI(X)bdΩ+∫ΓtNI(X)t-dΓ,fIa=∫ΩNI(X)ψIbdΩ+∫ΓtNI(X)ψIt-dΓ,fKIb=∫Ω‍NI(X)(ϕKI(X)-ϕKI(XI))bdΩ+∫ΓtNI(X)(ϕKI(X)-ϕKI(XI))t-dΓ
with body forces b and applied tractions t-. The stiffness matrix is
(23)KIJij=∫Ω‍(B¯Ii)TCB-JjdΩ,i,j=u,a,b,
where B¯Ii, i=u,a,b is the B operator defined by:
(24)BIu=∇NI(X),BIa=∇NI(X)ψI,BIb=∇[NI(X)(ϕKIl(X)-ϕKIl(XI)]).
For cohesive cracks without tip enrichment, the terms associated to the tip enrichment are omitted facilitating the implementation of the method.Although the XFEM approximation is capable of representing crack geometries that are independent of element boundaries, it relies on the interaction between the mesh and the crack geometry to determine the sets of enriched nodes [164]. This leads to particular crack configurations that cannot be accurately be represented by (13). Such cases are shown in Figure 16. As the crack size approaches the local nodal spacing, the set 𝒮c of nodes for the Heaviside or step enrichment is empty, Figure 16(b). Moreover, node 1 for the cracking case in Figure 16(a) or nodes 1 to 4 for the cracking case in Figure 16(b), respectively, contain two branch enrichments. Thus, the standard approximation gets difficulties with this crack configuration since the discontinuous function rsin(θ/2) extends too far. This problem always arises when one or more nodal supports contain the entire crack geometry. Usually, this kind of problem arises whenever a crack nucleates. Similar difficulties occur for approximations without any crack tip enrichment; see Figure 16(d). In order to close the crack within a single element, the set 𝒮c is empty as well. One solution is to refine the mesh locally such that the characteristic element size falls below that of the crack. An admissible crack configuration is shown in Figure 17.(a) and (b) Crack length that approach the local element size cannot be accurately represented by the standard XFEM approximation. Dots denote single enriched nodes and squares denote double (in our case, the node will contain the enrichment of two-crack tips) enriched nodes; (c) the dashed line shows the effective crack length; (d) even if no crack tip enrichment is used, in order to close the crack within a single element, no nodes have to be enriched with a step function.
(a)
(b)
(c)
effective crack length
(d)Figure 17
Admissible crack representation.More details on XFEM and its applications can be found in the excellent review papers by Karihaloo and Xiao [165], Fries and Belytschko [166], Belytschko et al. [167]  or the book from Mohammadi [168].


## 3.4.1. Representing the Crack Surface in XFEM

XFEM is a method that ensures crack path continuity. Enforcing crack path continuity is especially challenging in three dimensions. The crack path can be represented explicitly, usually with piecewise straight/planar crack segments [147, 148] or with other techniques such as level sets [149, 150] that trace the crack path in a more elegant way. In many publications, XFEM is closely related to level sets though any other technique can be used to describe the crack’s topology.Level sets are particularly attractive when a crack tip enrichment is used as geometry quantities can be calculated through the level sets. For example, the distance of a material point to the crack tip isr=f2+g2 and θ=arctan(g/f), see Figure 12; g is a level set orthogonal to f (i.e., ∇f·∇g=0) and is needed to uniquely identify the position of a material point with respect to the crack surface. Often, the level-set function is discretized
(18)fh(X)=∑I∈𝒮NI(X)fI;gh(X)=∑I∈𝒮NI(X)gI,
where commonly the shape functions NI(X) of the mechanical variables are used. Therefore, the enrichment can be expressed in terms of nodal values of both level set functions. Prabel et al. [151] used a (finer) finite difference scheme for the discretization of the level set function. However, the implementation effort is higher when different meshes are employed to discretize the level set and the mechanical variables. On the other hand, errors in the crack’s geometry will be introduced for low-order finite elements and curved cracks, Figure 14. As in most applications, a straight crack segment is attached in front of an existing crack front; a low-order finite elements should be sufficient for most crack propagation problems.Figure 14
Curved crack and its representation with level set functions when the level set is discretized with linear triangles.


## 3.4.2. Integration in XFEM

Numerical integration in XFEM requires special attention, in particular, around the crack tip when employing non-polynomial enrichment functions. Laborde et al. [152] reported that several thousand Gauss-points are required to adequately reduce integration error in a “brute-force” approach. For step-enriched (or Heaviside-enriched) elements, a simple sub-triangulation is sufficient. Ventura [153] proposed an integration scheme for piece-wise straight cracks that avoid a sub-triangulation. Another alternative is to increase the number of integration points [151]. For singularities, Ventura et al. [154], Gracie et al. [155], Bordas et al. [156], Cheng and Fries [157] proposed efficient integration schemes. One of the most efficient schemes is the “almost polar integration” by Nagarajan and Mukherjee [158]. The almost polar integration is based on mapping a rectangle onto a triangle and avoids integration of the crack tip singularity. It was exploited for the first time in the context of XFEM by Laborde et al. [152]  and simultaneously by Béchet et al. [159].


## 3.4.3. Blending Elements

Another major concern in XFEM are the so-called blending elements, Figure15. Blending elements are the elements adjacent to the enriched elements. It was shown by Chessa et al. [160], Stazi et al. [161], and Fries [162] that those blending elements result in the deterioration of the convergence rate in XFEM. When only a step enrichment is employed, a simple shifting can avoid problems due to blending. However, any non-constant enrichment function causes difficulties. Numerous approaches have been developed to eliminate drawbacks associated to blending [154, 160, 163]. A very simple and efficient approach was developed by Fries [162]. The idea is based on the modification of the enrichment function. Therefore, he enriched the nodes in the blending elements. He guaranteed that the enrichment in the enriched elements remains unchanged but vanishes in the “standard” elements.Figure 15
XFEM discretization with enriched domainΩenr, blending domain Ωblnd,  and the total (standard) domain Ωstd.


## 3.4.4. Application to Fracture

Applying XFEM to linear problems in elastostatics, the final discrete equations for a single crack are obtained by substituting the test and trial functions into the weak form of the equilibrium equation yielding(19)[KIJuuKIJuaKIJKubKIJauKIJaaKIJKabKIJKbuKIJKbaKIJKbb]{uJaJbJK}={fIextfIextfIKext}
or
(20)KIJdJ=fIext,
with the stiffness matrix
(21)KIJ=[KIJuuKIJuaKIJKubKIJauKIJaaKIJKabKIJKbuKIJKbaKIJKbb],
where the superscripts u, a  and b indicate the “standard”, “step-enriched” and “tip-enriched” part of the discretization; dJ={uJaJbKJ}T is the vector with the nodal parameters, fIext={fIufIafKIb}T is the external force vector with  fKIb={f1Ibf2Ibf3Ibf4Ib} and
(22)fIu=∫ΩNI(X)bdΩ+∫ΓtNI(X)t-dΓ,fIa=∫ΩNI(X)ψIbdΩ+∫ΓtNI(X)ψIt-dΓ,fKIb=∫Ω‍NI(X)(ϕKI(X)-ϕKI(XI))bdΩ+∫ΓtNI(X)(ϕKI(X)-ϕKI(XI))t-dΓ
with body forces b and applied tractions t-. The stiffness matrix is
(23)KIJij=∫Ω‍(B¯Ii)TCB-JjdΩ,i,j=u,a,b,
where B¯Ii, i=u,a,b is the B operator defined by:
(24)BIu=∇NI(X),BIa=∇NI(X)ψI,BIb=∇[NI(X)(ϕKIl(X)-ϕKIl(XI)]).
For cohesive cracks without tip enrichment, the terms associated to the tip enrichment are omitted facilitating the implementation of the method.Although the XFEM approximation is capable of representing crack geometries that are independent of element boundaries, it relies on the interaction between the mesh and the crack geometry to determine the sets of enriched nodes [164]. This leads to particular crack configurations that cannot be accurately be represented by (13). Such cases are shown in Figure 16. As the crack size approaches the local nodal spacing, the set 𝒮c of nodes for the Heaviside or step enrichment is empty, Figure 16(b). Moreover, node 1 for the cracking case in Figure 16(a) or nodes 1 to 4 for the cracking case in Figure 16(b), respectively, contain two branch enrichments. Thus, the standard approximation gets difficulties with this crack configuration since the discontinuous function rsin(θ/2) extends too far. This problem always arises when one or more nodal supports contain the entire crack geometry. Usually, this kind of problem arises whenever a crack nucleates. Similar difficulties occur for approximations without any crack tip enrichment; see Figure 16(d). In order to close the crack within a single element, the set 𝒮c is empty as well. One solution is to refine the mesh locally such that the characteristic element size falls below that of the crack. An admissible crack configuration is shown in Figure 17.(a) and (b) Crack length that approach the local element size cannot be accurately represented by the standard XFEM approximation. Dots denote single enriched nodes and squares denote double (in our case, the node will contain the enrichment of two-crack tips) enriched nodes; (c) the dashed line shows the effective crack length; (d) even if no crack tip enrichment is used, in order to close the crack within a single element, no nodes have to be enriched with a step function.
(a)
(b)
(c)
effective crack length
(d)Figure 17
Admissible crack representation.More details on XFEM and its applications can be found in the excellent review papers by Karihaloo and Xiao [165], Fries and Belytschko [166], Belytschko et al. [167]  or the book from Mohammadi [168].


## 3.5. Phantom Node Method

An alternative to the standard XFEM for strong discontinuities was proposed by A. Hansbo and P. Hansbo [169]. They do not model the crack with additional degrees of freedom but by overlapping elements; see Figure 18. In 1D, it was shown by Song et al.  [170] that the crack kinematics of the A. Hansbo and P. Hansbo [169] version of XFEM can be derived from “standard” XFEM. However, the Hansbo version of XFEM has some advantages over standard XFEM.(i)
It avoids the “mixed” termsKua and Kau (or Mua and Mau in dynamics), (19), facilitating the implementation.(ii)
It leads to a better conditioned system matrix.(iii)
Standard mass lumping schemes such as the row-sum technique can be used. Note that, efficient mass lumping schemes have meanwhile been proposed for XFEM [171, 172] (and XEFG [173]).(iv)
Based on the full interpolation basis of the overlapping elements, the Hansbo-XFEM can straightforwardly integrate enhanced techniques that are more difficult to integrate in an enriched XFEM formulation [174].(v)
The implementation of the Hansbo XFEM is simpler. It has already been implemented in the commercial software package ABAQUS. In ABAQUS, it has to be used in combination with cohesive zone models. The ABAQUS implementation of the Hansbo-XFEM can handle numerous cracks (crack propagation) in statics. In contrast, XFEM with tip enrichment has also been incorporated in ABAQUS. However, a crack propagation algorithm for the tip-enriched XFEM is yet not available in ABAQUS. A plugin-in (Morfeo) that can handle three-dimensional crack growth in LEFM was recently developed by the company Cenaero.On the downside, we notice. (i)
It is difficult to apply the Hansbo XFEM to other problems besides cracks.(ii)
It is also difficult to incorporate a tip enrichment for the Hansbo-XFEM. A crack tip element has been developed for problems in 2D [175]. However, the extension of such a model to 3D is cumbersome.The Hansbo version of XFEM was subsequently developed in a static setting by Mergheim et al. [176], Areias et al. [148], Mergheim and Steinmann [177], and in dynamics by Song et al. [170] who called their approach phantom node method. The phantom node method is well suited for non-linear problems with cohesive cracks and problems with numerous cracks, in particular when the method needs to be implemented in an existing finite element code. It is less accurate for problems in LEFM due to the absence of a tip enrichment. Subsequently, the basic idea of the phantom node method is briefly summarized.Figure 18
The principle of the phantom node method in which the hatched area is integrated to build the discrete momentum equation; the solid circles represent real nodes and the empty ones phantom nodes.Consider a body that is cracked as shown in Figure18 and the corresponding finite element discretization. There are elements cut by the crack. To have a set of full interpolation bases, the part of the cracked elements which belongs in the real domainΩ0 is extended to the phantom domainΩp. Then, the displacement in the real domain Ω0 can be interpolated by using the degrees of freedom for the nodes in the phantom domain Ωp. The nodes are called the phantom nodes and marked by empty circles in Figure 18. The approximation of the displacement field is then given by [170]:
(25)uh(X,t)=∑I∈{𝒲0+,𝒲P-}uI(t)NI(X)H(f(X))+∑J∈{𝒲0-,𝒲P+}uJ(t)NJ(X)H(-f(X)),
where f(X) is the signed distance measured from the crack, 𝒲0+, 𝒲0-, 𝒲P+ and 𝒲P- are nodes belonging to Ω0+, Ω0-, ΩP+ and ΩP-, respectively. H(x) is the Heaviside function. As can be seen from Figure 18, cracked elements have both real nodes and phantom nodes. The jump in the displacement field is realized by simply integrating only over the area from the side of the real nodes up to the crack, that is, Ω0+ and Ω0-. We note that there are certain similarities to the visibility method used originally in meshfree methods [178–181].


## 3.6. Cohesive Segment Method and Cracking Node Method

Remmers et al. [182, 183], proposed the cohesive segment method that has similarities to XFEM. The crack is represented by a set of overlapping cohesive segments that cross three entire elements. In contrast to XFEM, no crack path continuity is required in the Cohesive Segment Method. A similar approach was developed by Song and Belytschko [184] in the Cracking Node Method. In the Cracking Node Method, the crack is required to directly pass through the node of an element. Also the cracking node method does not require crack path continuity.


## 3.7. Meshfree Methods

Thanks to the absence of a mesh, meshfree methods (MMs) offer another alternative to model fracture. They are particularly well suited for dynamic fracture and large deformations. Moreover, adaptive h-refinement can be easily implemented in a meshfree context [185–187]. Though all methods discussed so far are capable of capturing discrete cracks, a certain refinement around the crack front is still needed for accuracy and efficiency reasons [188]. The meshfree approximation is similar to the FEM formulation:
(26)uh(X)=∑J∈𝒮NJ(X)uJ,
where NJ(X) denotes the meshfree shape function and uJ are the nodal parameters. Most meshfree methods are not interpolatory complicating the imposition of essential boundary conditions. Popular meshfree methods include the smoothed particle hydrodynamics (SPH) method [189], the element-free galerkin (EFG) method [190], the reproducing kernel particle method (RKPM) [191], the point interpolation method (PIM) [192], the meshless local-petrov galerkin method (MLPG) [193] or the Finite Point Method [194], among others.Most meshfree shape functions depend on a weighting or kernel function, which is denoted byWJ(X). Usually, radial kernel functions are employed where WJ(X)=W(r0) is expressed in terms of r0=∥X-XJ∥. The kernel functions are chosen to have compact support, that is, WJ(X)=0  if  ∥X-XJ∥>hJ, where hJ is the “interpolation" radius often called dilation parameter. The kernel function can be expressed in terms of material coordinates X (Lagrangian kernel) or spatial coordinates x (Eulerian kernel). Formulations based on Lagrangian kernels are more efficient as the shape functions can be constructed at the beginning of the simulation while Eulerian kernels require updating during the course of the simulation. More details on meshfree methods can be found in the literature, for example, in several review articles and books [195–202]. Before discussing different possibilities to incorporate fracture into MMs, the shape functions of the EFG-method will be given as it is helpful for the explanations in the subsequent sections. They are constructed as follows:
(27)NJ(X)=pT(X)·A(X)-1·D(XJ),A(X)=∑Jp(XJ)pT(XJ)W(r-J;h),D(XJ)=∑Jp(XJ)W(r-J;hJ),
where p is the polynomial basis. Usually, a linear polynomial basis is choosen, given in 2D by
(28)p=[1,X,Y].


### 3.7.1.  “Natural” Fracture

Due to the absence of a mesh, fracture in meshfree methods can occur naturally [38, 203, 204]. In contrast to element erosion, there is no need to erode particles. Fracture occurs when particles move outside the domain of influence of neighboring particles. Libersky et al. [203] presented impressive results of dynamic fracture and fragmentation predicting fragment distributions quite accurately. Dilts [205] presented an algorithm to detect fracture surfaces within the SPH method in 2D. He extended the approach later to three dimensions [206] and solved impressive impact problems. However, care has to be taken as fracture can occur artificially as shown for example, by Rabczuk et al. [207]. This is mostly the case when Eulerian kernels are employed. Formulations based on Lagrangian kernels do not suffer from numerical fracture. However, as the neighbor list does not change during the course of the simulation, fracture cannot be modeled “naturally” any more. Moreover, the Lagrangian kernel has to be consistently updated for problems involving large deformations. A very efficient algorithm to treat dynamic fracture and fragmentation in meshfree methods based on Lagrangian kernels was proposed by Maurel and Combescure [208] and Combescure et al. [209]. In their work, connectivities between neighboring particles were broken when their distance exceeds a given threshold. They solved complex problems including fluid-structure interaction and fragmentation [210, 211]. Another effective and simple method to treat fracture in meshfree methods is the material point method (MPM) [212]. The advantage of the MPM over most SPH methods for fracture is their computational efficiency [213]. While those methods are quite attractive to model dynamic fracture and fragmentation involving an enormous amount of cracks, they are not well suited for static fracture with few cracks due to their less accurate representation of the crack kinematics.


### 3.7.2. Visibility Method

The visibility method is the first approach that introduces a discrete crack into the meshfree discretization. In the visibility method, the crack boundary is considered to be opaque. Thus, the displacement discontinuity is modeled by excluding the particles on the opposite side of the crack in the approximation of the displacement field, see Figure19:
(29)uh(X)=∑I∈𝒮+‍NI(X)uI.
This is identical to setting the shape function across the crack to zero as shown in Figure 19. The jump in the displacement is then computed by
(30)[[u(X)]]=(∑I∈𝒮+NI(X)uI)-(∑I∈𝒮-NI(X)uI).
Difficulties arise for particles close to the crack tip since undesired interior discontinuities occur (Figure 20) due to the abrupt cut of the shape function, see Figure 19. Nevertheless, Krysl and Belytschko [214] showed convergence for the visibility method. For linear complete EFG shape functions, they even showed that the convergence rate is not affected by the discontinuity. The length and size of the (undesired) discontinuities depend on the nodal spacing near the crack boundary. If the nodal spacing approaches zero, the length of the discontinuities tends to zero. With this argument and the theory of nonconforming finite elements, convergence of the discontinuous displacement field can be shown.Figure 19
Principle of the visibility, diffraction, and transparency methods with corresponding shape functions, from Belytschko et al. [200].(a) Undesired introduced discontinuities by the visibility method, (b) Difficulties with the visibility method for concave boundaries and kinks.
(a)
(b)It should also be noted that the visibility criterion leads to discontinuities in shape functions near nonconvex boundaries such as kinks, crack edges, and holes, as shown in Figure20 in two dimensions.An efficient implementation of the visibility method in 2D is given, for example, in [215].


## 3.7.1.  “Natural” Fracture

Due to the absence of a mesh, fracture in meshfree methods can occur naturally [38, 203, 204]. In contrast to element erosion, there is no need to erode particles. Fracture occurs when particles move outside the domain of influence of neighboring particles. Libersky et al. [203] presented impressive results of dynamic fracture and fragmentation predicting fragment distributions quite accurately. Dilts [205] presented an algorithm to detect fracture surfaces within the SPH method in 2D. He extended the approach later to three dimensions [206] and solved impressive impact problems. However, care has to be taken as fracture can occur artificially as shown for example, by Rabczuk et al. [207]. This is mostly the case when Eulerian kernels are employed. Formulations based on Lagrangian kernels do not suffer from numerical fracture. However, as the neighbor list does not change during the course of the simulation, fracture cannot be modeled “naturally” any more. Moreover, the Lagrangian kernel has to be consistently updated for problems involving large deformations. A very efficient algorithm to treat dynamic fracture and fragmentation in meshfree methods based on Lagrangian kernels was proposed by Maurel and Combescure [208] and Combescure et al. [209]. In their work, connectivities between neighboring particles were broken when their distance exceeds a given threshold. They solved complex problems including fluid-structure interaction and fragmentation [210, 211]. Another effective and simple method to treat fracture in meshfree methods is the material point method (MPM) [212]. The advantage of the MPM over most SPH methods for fracture is their computational efficiency [213]. While those methods are quite attractive to model dynamic fracture and fragmentation involving an enormous amount of cracks, they are not well suited for static fracture with few cracks due to their less accurate representation of the crack kinematics.


## 3.7.2. Visibility Method

The visibility method is the first approach that introduces a discrete crack into the meshfree discretization. In the visibility method, the crack boundary is considered to be opaque. Thus, the displacement discontinuity is modeled by excluding the particles on the opposite side of the crack in the approximation of the displacement field, see Figure19:
(29)uh(X)=∑I∈𝒮+‍NI(X)uI.
This is identical to setting the shape function across the crack to zero as shown in Figure 19. The jump in the displacement is then computed by
(30)[[u(X)]]=(∑I∈𝒮+NI(X)uI)-(∑I∈𝒮-NI(X)uI).
Difficulties arise for particles close to the crack tip since undesired interior discontinuities occur (Figure 20) due to the abrupt cut of the shape function, see Figure 19. Nevertheless, Krysl and Belytschko [214] showed convergence for the visibility method. For linear complete EFG shape functions, they even showed that the convergence rate is not affected by the discontinuity. The length and size of the (undesired) discontinuities depend on the nodal spacing near the crack boundary. If the nodal spacing approaches zero, the length of the discontinuities tends to zero. With this argument and the theory of nonconforming finite elements, convergence of the discontinuous displacement field can be shown.Figure 19
Principle of the visibility, diffraction, and transparency methods with corresponding shape functions, from Belytschko et al. [200].(a) Undesired introduced discontinuities by the visibility method, (b) Difficulties with the visibility method for concave boundaries and kinks.
(a)
(b)It should also be noted that the visibility criterion leads to discontinuities in shape functions near nonconvex boundaries such as kinks, crack edges, and holes, as shown in Figure20 in two dimensions.An efficient implementation of the visibility method in 2D is given, for example, in [215].


## 3.8. The Diffraction Method

The diffraction method is an improvement of the visibility method. It eliminates theundesired interior discontinuities; see Figure 21 (see also Figure 19). The diffraction method is also suitable for non-convex crack boundaries. It is motivated by the way light diffracts around a sharp corner, but the equations used in constructing the domain of influence and the weight function bear almost no relationship to the equation of diffraction. The method is only applicable to radial basis kernel functions with a single parameter hI0.The diffraction method.
(a)
(b)The idea of the diffraction method is not only  to treat the crack as opaque but also to evaluate the length of the rayh0 by a path which passes around the corner of the discontinuity. A typical weight function is shown in Figure 19. It should be noted that the shape function of the diffraction method is quite complex with several areas of rapidly varying derivatives that complicates quadrature of the discrete Galerkin form [200]. Moreover, the extension of the diffraction method into three dimensions is complex. Nevertheless, several authors have presented implementations of the diffraction method in 3D [216, 217].


## 3.9. The Transparency Method

The transparency method was developed as an alternative to the diffraction method by Organ et al. [178]. The transparency method is easier extendable into three dimensions than the diffraction method. In the transparency method, the crack is made transparent near the crack tip. The degree of transparency is related to the distance from the crack tip to the point of intersection; see Figure 19.An additional requirement is usually imposed for particles close to the crack. Since the angle between the crack and the ray from the node to the crack tip is small, a sharp gradient in the weight function across the line ahead of the crack is introduced. In order to reduce this effect, Organ et al. [178] imposed that all nodes have a minimum distance from the crack surface, that is, the normal distance to the crack surface must be larger than  γh  with 0<γ<1; they suggested γ=1/4.


## 3.10. The “See Through” and “Continuous Line” Method

The “see-through” method was proposed by Terry [218] for constructing continuous approximations near nonconvex boundaries. Therefore, the boundary was considered as completely transparent such that the discontinuity is removed. Though the “see-through” method works well for capturing features such as interior holes, it is not well suited to model cracks.In the continuous line method from Krysl and Belytschko [214] and Duarte and Oden [219], the crack is completely transparent at the crack tip. In other words, particles with a partially cut domain of influence can see through the crack. This drastically shortens the crack. If no special techniques are introduced, the crack also does not close at the crack tip that leads to inaccurate solutions. Crack closure at the crack tip can be enforced with Lagrange multiplier or by decreasing the domain of influence of nodes close to the crack tip.Belytschko and Fleming [220] suggest to combine different methods depending on the convexity of the crack boundary for example, to use the visibility for convex boundaries, and other methods for non-convex crack boundaries. They suggest a criterion based on the angle of the wedge that can be written in terms of the surface normal; see Figure 22. When  nA·nB≤β with β=0° as a cutoff value, the boundary can be considered as convex, otherwise non-convex. The cutoff value of β=0° corresponds to the wedge angle of  ω=90° in Figure 22.Figure 22
Domain of influence near a wedge-shaped non-convex boundary. The boundary is enforced ifnA·nB≤β.


## 3.11. Extended Meshfree Methods

Enrichment in meshfree methods was used before the development of XFEM. A good overview paper on enriched meshfree methods in the context of LEFM is given by Fleming et al. [221], see also the review articles by Belytschko et al. [200], Nguyen et al. [201], and Huerta et al. [202]. The enrichment schemes are classified into intrinsic enrichment and extrinsic enrichment. An intrinsic enrichment introduces enrichment functions as additional functions to the polynomial basis compare to (28) as the following:
(31)p=[1,X,Y,rsin(θ2),rcos(θ2),rsin(θ2),sinθ,rcos(θ2),sinθ].
This concept was later also used in the intrinsic XFEM [222]. The major drawback of this approach is its higher computational cost as the enrichment needs to be employed everywhere to avoid artificial discontinuities. Techniques and blending enriched approximations to “standard” approximations have been used but they introduce additional complexity. Duflot and Nguyen-Dang [223] proposed an alternative intrinsic enrichment by enriching the kernel functions. Their approach does not require blending techniques and is computationally cheaper. However, it is also less accurate. The extrinsic enrichment can be categorized into extrinsic MLS enrichment and extrinsic PU enrichment. The extrinsic MLS enrichment introduces enrichment parameters that can be interpreted as stress-intensity factors (SIFs), at least when the enrichment is introduced globally. Hence, the SIFs can be obtained directly from the numerical analysis without computing the J-integral or interaction integral. However, it was shown by Zamani et al. [224, 225] that such approaches are less accurate. The extrinsic PU enrichment is comparable to an XFEM enrichment. It is believed that it is the best choice in terms of accuracy, efficiency, and flexibility.An extended element-free Galerkin (XEFG) method based on vector level sets was first proposed by Ventura et al. [226] in the context of LEFM. It was later extended to nonlinear materials and cohesive cracks by Rabczuk and Zi [227], Rabczuk et al. [228]. In contrast to XFEM, due to the heavily overlapping shape functions, a crack tip enrichment has to be employed to ensure crack closure at the crack tip though alternative methods have been developed that avoid the use of tip enrichments [229, 230]. However, those approaches are cumbersome in three dimensions.


## 3.12. Cracking Particles Method

While all previous meshfree approaches for fracture require crack path continuity, the cracking particles method [231–235] does not represent the crack as continuous surface. This makes the method particularly useful for 3D applications including complex crack patterns such as crack branching and crack coalescence. In contrast to methods that enforce crack path continuity, crack branching and crack coalescence are a (more) natural outcome of the analysis in the cracking particles method.To model cracks in the cracking particles method, the displacement is decomposed into continuous and discontinuous parts:(32)u(X)=ucont(X)+uenr(X),
where X are the material coordinates, ucont denotes the continuous displacement, and uenr denotes the discontinuous part.The crack is modelled by a set of discrete cracks as shown in Figure23. These discrete cracks are restricted to lie on the particles, that is, each crack plane (or line in 2D) always passes through a particle. Since the crack geometry is described by the set of cracked particles, a representation for the geometry of the crack is not needed.Schematic on the right shows a crack model for the crack on the left.
(a)
(b)The approximation of the displacement field in the cracking particles method is given by(33)uh(X)=∑I∈𝒩NI^(X)uI+∑I∈𝒩cNI(X)s(fI(X))qI,
where 𝒩 is the total set of nodes in the model and 𝒩c  is the set of cracked nodes;  s(fI(X)) is the step function. In general, different shape functions can be used for the continuous and discontinuous part of the displacement field.The jump in the displacement across the crack depends only on the discontinuous part of the displacement fielduenr(X) and hence the enrichment parameters qI. It can be shown from (33) that
(34)⟦u(X)⟧=2∑I∈𝒩cNI(X)qI.
The factor 2 comes from the fact that the step function is used that is −1 on one side of the crack and +1 on the other side of the crack.The cracking particles method might suffer from spurious crack patterns [236, 237]. However, they can be eliminated by a set of simple crack injection rules [184]. Firstly, one has to distinguish between crack propagation and crack initiation. Since the crack is not considered as continuous surface, crack propagation is assumed when no cracked node is detected in the vicinity (in a circle [sphere in 3D] with radius r=αh, α≥1.5) of a sampling point, where cracking is detected. For propagating cracks, spurious cracks as shown in Figure 24 are avoided by introducing a prevention zone, Figure 25. For example, for particle L, the crack prevention zone is given through the angle α:
(35)|n·x-xL∥x-xL∥|<cos(Π2-α).
For branching cracks, it is often observed that the orientation of neighboring cracks severely differ, see node K in Figure 25. Therefore, cracks are also allowed in the prevention zone when the following criterion is fulfilled:
(36)nL·nK<cosβ,
with α=2Π/9 and β=Π/3. Similar criteria are used within the cracking particles method in H. X. Wang and S. X. Wang [238]. Zhang [237] reported that cracking rules can be avoided by the use of rotating crack segments.(a) Spurious cracking and (b) improved crack pattern.
(a)
(b)Figure 25
Crack prevention: nodesJ and K are located in the crack prevention zone; node I is located in the crack propagation zone.Recently, the cracking particles method has been presented without enrichment to model fracture in continua [235] and structures [239]. This results in less degrees of freedom and facilitates diagonalizing the mass matrix when the method is used in explicit dynamics.


## 3.13. Boundary Element Method

In the Boundary Element Method (BEM), the weak form is formulated in boundary integral form. It reduces therefore the dimension of the problem, for example, a three-dimensional problem reduces to two dimensions. However, the BEM is only applicable to problems where Green’s functions can be computed restricting the application range of the method. Due to the dimension-reduction, the remeshing procedure is fairly simple as the crack is part of the boundary. Interesting applications of the BEM to fracture can be found, for example, in [240–257].


## 3.14. Isogeometric Analysis for Fracture

Isogeometric analysis (IGA) [258] was developed to unify CAD and CAE. Verhoosel et al. [130] were the first who incorporated fracture in isogeometric analysis based on T-splines. Therefore, they used the concept of knot insertion to introduce discontinuities in the displacement approximation. This method suffers from similar problems as interelement-separation methods. Moreover, for complex crack patterns with slightly deviating angles between the crack branches, it is expected that the mesh will be extremely distorted. De Luycker et al. [259], Ghorashi et al. [260] and recently Tambat and Subbarayan [261] combined XFEM with isogeometric analysis (XIGA) allowing crack growth without remeshing also in the context of IGA.


## 3.15. The Variational Approach to Fracture

Variational methods in fracture mechanics are a relatively new development in the field. The underlying theory was laid down by Francfort and Marigo [262], who proposed the idea that cracks should propagate along a path of least energy. The goal of variational methods is to circumvent several important weaknesses of the classical theory emanating from the work of Griffith in 1920.In classical fracture mechanics, crack propagation is assumed to occur under thermodynamic equilibrium. In particular for crack formation occurring under constant load, the external work done by the applied loading is equal to twice the energy in the bulk, so that the total energy of the system may be expressed asU=-WL+UE+US=-UE+US, where US is a postulated quantity known as the surface (dissipation) energy of the crack system. Thermodynamic equilibrium with respect to crack growth then requires that dU/da=0, where a is the crack length. The previous expression provides the critical condition for fracture. The simplest case to consider is an infinite solid with a straight crack of length 2a and subjected to far-end tractions of magnitude σ directed perpendicular to the crack. For this problem, the exact solutions for the stresses and strains are available due to Inglis, so that an analytic form of UE may be obtained, equal to πa2σL2/E for the case of plane stress. We now assume that the surface energy to be US=4γa,  where γ is the surface energy per unit area, which is considered a material property. Enforcement of thermodynamic equilibrium yields the expression for the critical applied stress required for crack propagation: σfa=2Eγ/π. According to the classical theory, the existing crack will not propagate for σ<σf and will do so for σ=σf.Despite improvements made by subsequent authors to account for nonlinearity and inelasticity, the original form of the fracture criterion remained, that is,σfa=constant. As pointed out in [262], this inverse proportionality of the critical stress to the square root of the initial crack length means that classical fracture mechanics is unable to predict crack initiation, since for a body without an initial crack (a=0), σf becomes infinite. Other important weaknesses of the classical theory listed in [262] are its inadequacy for predicting the direction of the crack path, and its inability to handle crack jumps (brutal cracking). The former arises from the fact that Griffith’s criterion gives only one constraint, whereas the description of the crack tip propagation requires a number of functions equal to the dimensionality of the problem. The latter weakness stems from the classical theory being silent on what happens when σ>σf (unstable cracking). Griffith’s criterion is also incapable of predicting when a crack will branch, so that additional criteria are required in order to evaluate the possibility of this kind of response.In order to address the above mentioned weaknesses of classical fracture mechanics, Francfort and Marigo [262] proposed an alternative formulation of brittle fracture based on energy minimization concepts. For a body denoted by Ω- in D-dimensional space, energy (similar to the surface energy concept of Griffith) is assigned to the crack system Γ∈Ω- and expressed as
(37)Es(Γ)=∫Γk(x)dℋN-1(x),
where ℋN-1 denotes the (N-1)-D Hausdorff measure. The function k(x) represents the energy required to create an infinitesimal crack at the point x  and is strictly positive and bounded away from zero. The above definition is rather accommodating in terms of crack geometries, and in addition to the usual “sharp” crack surfaces of classical fracture mechanics, it also admits entities like point clouds and edge cracks. However, cracks having N-D Hausdorff measure (termed as “fat” cracks in [262]) are not allowed, for example a crack with computable area in a 2D or volume in 3D, since these will require an infinite amount of energy to propagate. The bulk energy is defined as
(38)Ed(Γ,U)=infv∈𝒞(Γ,U)∫Ω∖ΓW(x,ϵ(v,(x)))dx,
in which 𝒞(Γ,U) is the set of kinematically admissible displacement fields. The total energy of the body is the sum of these two terms:
(39)E(Γ,U)=Ed(Γ,U)+Es(Γ).
The problem then is to find the mapping Γ(t) such that at some given time t>0, the crack geometry Γ(t) is the one that minimizes E(U,Γ(t)) and furthermore contains all Γ(s) for s<t (condition of irreversibility). The authors warn that the use of global energy minimization to drive the crack evolution is not based on thermodynamic principles, but is adapted as “… a rather convenient postulate which provides useful insight into a variety of behaviors” [262]. They show also that this new formulation is free of the earlier identified weaknesses of the classical theory.Numerical experiments on the proposed formulation were carried out by Bourdin et al. [263], who pointed to the similarity of the new formulation to models of image segmentation obtained by minimization of the Mumford-Shah functional. This meant that the inapplicability of standard numerical methods, owing to the fact that the formulation allowed for jump sets of the displacement field (representing the cracks) whose locus was unknown a priori. Two methods of solution were demonstrated, based on the mathematical concept of Γ-convergence. In one of these methods, a second field was included in the energy functional in addition to the primary displacement field. This auxiliary scalar field acted as a regularizer for the jump sets of the displacement field, and took on values of 0 on the crack and 1 away from it. The regularizing variable was referred to by later researchers as the phase field, and it allowed for the treatment of the global energy minimization as a standard variational problem for which classical FEM is up to the task, albeit with the restriction that the characteristic length of the mesh should tend to zero faster than the characteristic length of the regularization so as not to overestimate the surface energy of the crack. Due to the inclusion of a second field, a coupled system of equations must now be solved consisting of the original equilibrium/linear momentum equations and the evolution PDE of the phase field.A phase field model for mode III dynamic fracture was devised by Karma et al. [264], which made use of a phase field evolution equation based on the standard two-minimum Ginzburg-Landau form, that is,
(40)τ∂tϕ(x,t)=Dϕ∇2ϕ-ddϕ(14ϕ2(1-ϕ2))-μ2g′(ϕ)(ϵ:ϵ-ϵc2),
where g is a function obeying the constraints g(0)=0, g(1)=1, and g′(0)=g′(1)=0. The KKL phase field model was subsequently utilized by Hakim and Karma [265] to analyze the laws of quasistatic crack tip motion. Among other things, they showed that for kinked cracks in anisotropic media, force-balance gives predictions that are significantly different from those obtained using the principle of maximum energy release rate, and that the predictions obtained from force-balance hold even when fracture is modeled as irreversible. It should be mentioned that up to this point, phase field implementations featured isotropic damage wherein fracture occurred both in tension and compression. This sometimes resulted in unphysical response, such as sample interpenetration in the crack branching simulation in Bourdin et al. [263].Thermodynamically consistent phase field models of fracture were developed by Miehe et al. [266, 267], along with corresponding incremental variational principles and numerical implementation within a finite element framework. Their implementation of the phase field was also slightly different from that of previous authors whereas Bourdin et al. [263] and subsequent papers utilized the convention that the phase field took values of 0 at cracks and 1 at unbroken states (hence a pseudo material density function) and Miehe et al. reversed the convention by assigning to the phase field values of 1 and 0 for fully cracked and fully unbroken states respectively. As a consequence, the functional associated with the phase field may now be seen to represent the crack surface itself. They then derived the evolution equation for the phase field based on the assumption that the solution is negative exponential in nature. Furthermore, they extended the application of phase fields to fracture simulations involving viscous overforce response via a time-regularized three-field formulation, and showed that the coupled system (linear momentum + phase field evolution) could be solved either monolithically or via a robust staggered-update scheme [266, 267].An important addition of Miehe et al. [266, 267] to the existing theory was the modelling of anisotropic degradation by using an additive decomposition of the stored energy of the undamaged solid, to come up with an anisotropic energy function of the form
(41)Ψ(ϵ,ϕ)=[g(ϕ)+k]Ψ0+(ϵ)+Ψ0-(ϵ),
where the positive and negative parts of the energies are defined by
(42)Ψ0±(ϵ)=λ2〈ϵ1+ϵ2+ϵ3〉±2+μ(〈ϵ1〉±2+〈ϵ2〉±2+〈ϵ3〉±2),
and 〈x〉± are ramp functions defined as 〈x〉±=(1/2)[|x|±x]. The energy decomposition allows for the case where fracture occurs in tension only as opposed to both tension and compression, and effectively fixes the sample inter-penetration issue encountered by Bourdin in crack branching simulations. Following the reversed convention on the phase field, the constraints on g(ϕ) are g(0)=1, g(1)=0, and g′(1)=0. To maintain numerical stability, a small positive constant k is included in the formulation so that the material retains some residual stiffness even when it attains a fully damaged state.One challenge with the use of phase field models for fracture is the computational expense associated with mesh size requirements, since the use of a mesh having a characteristic length that is not small enough compared to the crack regularization parameter yields erroneous results with regard to the energy. Based on numerical experiments utilizing the stationary phase field equation, Miehe et al. [266, 267] found that the characteristic length of the mesh should be no more than half the value of the regularization parameter, at least presumably in the vicinity of the regularized crack. In an effort to reduce computational expense, Kuhn and Müller [268, 269] (see also their work in [270]) introduced specially engineered FEM shape functions of an exponential nature to discretize the phase field. These shape functions have the form
(43)N1e(ξ,δ)=1-exp(-δ(1+ξ)/4)-1exp(δ/2)-1,N2e(ξ,δ)=exp(-δ(1+ξ)/4)-1exp(δ/2)-1,
in 1D for the 2-node line element. 2D shape functions for the 4-node quadrilateral are obtained as a tensor products of the above, with certain tweaks. The shape functions are parametrized by δ which is the ratio between the element size and phase field regularization parameter. Their results showed that with the special element shape functions, accurate prediction of surface energy associated with the phase field is possible with a much lower level of refinement compared to standard FEM shape functions. The drawback to their approach is that some information regarding crack orientation is necessary for the proper construction of the exponential shape functions, so that numerical examples were confined to cases, where the general direction of the crack path is known a priori.


## 4. Tracking the Crack Path

Methods that ensure crack path continuity such as XFEM require the representation of the crack surface and algorithms to track the crack path. While those tasks are relatively easy to implement in two dimensions, their implementation in 3D is challenging.The topology of the crack surface is commonly represented either explicitly by piece-wise planar crack segments or implicitly by level set functions. Meshfree methods and the phantom node method usually use the former method while many XFEM-implementations are based on an implicit representation of the crack surface using level sets.There are three major approaches to track the evolving crack surface:(i)
local methods,(ii)
global methods,(iii)
level set method. A review article of those methods including implementation details in the context of partition of unity methods is given by Rabczuk et al. [271].


### 4.1. Local Methods

With local crack tracking algorithms, the alignment of the crack surfaces is enforced with respect to its neighborhood. Local crack tracking algorithms are usually characterized by recursively “cutting” elements (or background cells in meshfree methods) and are especially effective in three dimensions. Local methods have difficulties to ensure a continuous crack surface in 3D. Jäger et al. [272] report, for this reason, they are not well suited to model curved cracks. However, adaptive refinement and smoothing techniques allow also the representation of curved surfaces (modelled by piece-wise planar crack segments) [186, 232]. Local methods were pursued in XFEM, for example, by Gasser and Holzapfel [273, 274] and Areias and Belytschko [147] or in the context of meshfree methods by Rabczuk and Belytschko [186, 232], Krysl and Belytschko [275], among others.


### 4.2. Global Methods

In global methods, a linear heat conduction problem is solved each load step for the mechanical problem that needs to fulfill the condition(44)a·n0=b·n0=0,
where a and b are two vector fields, and n0 denotes the vector of the normal to the crack surface. The family of surfaces, enveloping both vector fields a and b are described by a temperature-like function T(X) when the following conditions hold
(45)a·∇0T=∇0T·a=0inΩ0,b·∇0T=∇0T·b=0inΩ0.
Equation (45) can be rephrased as an anisotropic heat conduction problem [101, 276, 277]. The drawback of this method is that the heat conduction problem has to be solved at every time step making global methods computationally expensive.Feist and Hofstetter [114] proposed the domain crack tracking algorithm, a slight variant of the global method. In their algorithm, the scalar field T is not constructed globally but only for a subset of the domain that is affected by the crack. Jäger et al. [272] report that global algorithm are well suited to track curved cracks but less suitable for straight cracks.


### 4.3. Level Set Method

The level set method was originally developed to track interfaces that propagate orthogonal to their surface. To update the interface, the Hamilton-Jacobi equation was solved with respect to the level set. It makes the method particularly attractive for problems in fluid mechanics. However, the original level set method is not well suited to track crack surfaces for three reasons as noted, for example, by Ventura et al. [278] and Duflot [279]. (i)
The zero isobar of the level set must be updated behind the crack front to account for the fact once a material point is cracked, it remains cracked.(ii)
The crack surface is an open surface that extends during the crack propagation. It requires the introduction of another level set function (orthogonal to the level set function describing the crack surface) in order to uniquely determine the position of a material point with respect to the crack surface. When the crack propagates, this level set function needs to be updated as well. Moreover, both level set functions must be periodically reinitialized to the signed-distance functions to preserve stability [280, 281].(iii)
The level set functions are not updated with the speed of the interface in the normal direction, and hence the Hamilton-Jacobi equation cannot be used. Instead, the level set function propagates with the speed of the crack front.Crack propagation with level sets can be modelled by different techniques that can be classified into four groups [279]. In the first group, the level set is updated by the solution of differential equations. The second group is defined on algebraic relations between the coordinates of a given point, the coordinates of the crack front, and the crack advance vector, see for example, Stolarska et al. [282]. This method is not easily extendable into three dimensions. The Vector level set method [226, 278] is defined in terms of geometric transformations. In the vector level set method, the distance to the crack surface is stored in addition to the signed distance function. This facilitates the implementation since there is no need to solve a PDE to update the level set. The last class of methods are based on algebraic and trigonometric equations involving the initial value of the level set functions and the crack advance vector. Some of them also require the description of the crack front. They are well suited for three-dimensional applications. Duflot [279] suggested a mixed method defined with differential and non-differential equations. He also gives an excellent state-of-the-art overview on level set techniques for cracks.A very efficient and elegant crack propagation and track cracking algorithm in the context of XFEM was recently presented by Fries and Baydoun [283] and in the context of meshfree methods by Zhuang et al. [284, 285].


## 4.1. Local Methods

With local crack tracking algorithms, the alignment of the crack surfaces is enforced with respect to its neighborhood. Local crack tracking algorithms are usually characterized by recursively “cutting” elements (or background cells in meshfree methods) and are especially effective in three dimensions. Local methods have difficulties to ensure a continuous crack surface in 3D. Jäger et al. [272] report, for this reason, they are not well suited to model curved cracks. However, adaptive refinement and smoothing techniques allow also the representation of curved surfaces (modelled by piece-wise planar crack segments) [186, 232]. Local methods were pursued in XFEM, for example, by Gasser and Holzapfel [273, 274] and Areias and Belytschko [147] or in the context of meshfree methods by Rabczuk and Belytschko [186, 232], Krysl and Belytschko [275], among others.


## 4.2. Global Methods

In global methods, a linear heat conduction problem is solved each load step for the mechanical problem that needs to fulfill the condition(44)a·n0=b·n0=0,
where a and b are two vector fields, and n0 denotes the vector of the normal to the crack surface. The family of surfaces, enveloping both vector fields a and b are described by a temperature-like function T(X) when the following conditions hold
(45)a·∇0T=∇0T·a=0inΩ0,b·∇0T=∇0T·b=0inΩ0.
Equation (45) can be rephrased as an anisotropic heat conduction problem [101, 276, 277]. The drawback of this method is that the heat conduction problem has to be solved at every time step making global methods computationally expensive.Feist and Hofstetter [114] proposed the domain crack tracking algorithm, a slight variant of the global method. In their algorithm, the scalar field T is not constructed globally but only for a subset of the domain that is affected by the crack. Jäger et al. [272] report that global algorithm are well suited to track curved cracks but less suitable for straight cracks.


## 4.3. Level Set Method

The level set method was originally developed to track interfaces that propagate orthogonal to their surface. To update the interface, the Hamilton-Jacobi equation was solved with respect to the level set. It makes the method particularly attractive for problems in fluid mechanics. However, the original level set method is not well suited to track crack surfaces for three reasons as noted, for example, by Ventura et al. [278] and Duflot [279]. (i)
The zero isobar of the level set must be updated behind the crack front to account for the fact once a material point is cracked, it remains cracked.(ii)
The crack surface is an open surface that extends during the crack propagation. It requires the introduction of another level set function (orthogonal to the level set function describing the crack surface) in order to uniquely determine the position of a material point with respect to the crack surface. When the crack propagates, this level set function needs to be updated as well. Moreover, both level set functions must be periodically reinitialized to the signed-distance functions to preserve stability [280, 281].(iii)
The level set functions are not updated with the speed of the interface in the normal direction, and hence the Hamilton-Jacobi equation cannot be used. Instead, the level set function propagates with the speed of the crack front.Crack propagation with level sets can be modelled by different techniques that can be classified into four groups [279]. In the first group, the level set is updated by the solution of differential equations. The second group is defined on algebraic relations between the coordinates of a given point, the coordinates of the crack front, and the crack advance vector, see for example, Stolarska et al. [282]. This method is not easily extendable into three dimensions. The Vector level set method [226, 278] is defined in terms of geometric transformations. In the vector level set method, the distance to the crack surface is stored in addition to the signed distance function. This facilitates the implementation since there is no need to solve a PDE to update the level set. The last class of methods are based on algebraic and trigonometric equations involving the initial value of the level set functions and the crack advance vector. Some of them also require the description of the crack front. They are well suited for three-dimensional applications. Duflot [279] suggested a mixed method defined with differential and non-differential equations. He also gives an excellent state-of-the-art overview on level set techniques for cracks.A very efficient and elegant crack propagation and track cracking algorithm in the context of XFEM was recently presented by Fries and Baydoun [283] and in the context of meshfree methods by Zhuang et al. [284, 285].


## 5. Fracture Criteria

The fracture criterion is needed to determine whether or not a crack will propagate or nucleate. Moreover, the fracture criterion should provide the orientation and the “length” of the crack as well as whether or not cracks branch or join. Methods that ensure crack path continuity need to distinguish between crack nucleation and crack propagation. Detecting branching cracks in dynamic problems is particularly difficult in those methods. The fracture criterion is often met at several quadrature points in front of the crack tip, and reliable criteria to branch cracks in practice are still missing. The best results are obtained when the crack branches are known in advance.In most applications, the “length” of the crack is controlled. Usually, it is correlated to the underlying discretization. To the best knowledge of the author, all of the methods assume a straight/planar extension of an existing crack surface.In the following section, different cracking criteria to obtain the orientation of a newly (nucleated or propagated) crack segment are reviewed.


### 5.1. Fracture Mechanics-Based Criteria

Besides approaches based on configurational forces [286–289] there are four major cracking criteria in LEFM.(i)
Maximum hoop stress criterion or maximum principal stress criterion.(ii)
Minimum strain energy density criterion, Sih [290].(iii)
Maximum energy release rate criterion, Wu [291].(iv)
The zeroKII criterion (Vanishing in-plane SIF (KII) in shear mode for infinitesimally small crack extension), Goldstein and Salganik [292].The first two criteria predict the direction of the crack trajectory from the stress state prior to the crack extension. The last two criteria require stress analysis for virtually extended cracks in various directions to find the appropriate crack-growth direction. Duflot and Hung [223] report for medium mixed-mode problems, all criteria yield almost identical results. However, according to Shen and Stephansson [293], only the maximum energy release rate criterion allows to predict secondary cracks in compressed specimens. The crack is propagated in an angle of θc from the crack tip. Note that LEFM can only deal with crack propagation but not with crack nucleation. It also can not handle crack branching.In the maximum hoop stress or maximum principal stress criterion, the maximum circumferential stressσθθ, often called hoop stress, in the polar coordinate system around the crack tip corresponds to the maximum principal stress and is given for a crack propagating with constant velocity vc by:
(46)σθθ=KI2πrfhI(θ,vc)+KII2πrfhII(θ,vc),
where the functions fhI and fhII represent the angular variation of the stress for different values of the crack-tip speed vc. When the maximum hoop stress is larger equal a critical hoop stress σθθc then the crack is propagated in the direction perpendicular to the maximum hoop stress. For pure mode I fracture, σθθc is given by
(47)σθθc=KIc2πr,
with the fracture toughness KIc.In LEFM, the local direction of the crack growth is determined by the condition that the local shear stress is zero that leads to the condition:(48)KIsinθc+KII(3cosθc-1)=0,
that results in the crack propagation angle
(49)θc=2arctan(KI-KI2+8KII24KII).The minimum strain energy density criterion is based on a critical strain-energy-density factorSrc. The basic assumption is that fracture occurs when the interior minimum of the strain-energy-density factor S reached Src. The strain energy density factor S represents the strength of the elastic energy field in the vicinity of the crack tip which is singular of the order r-1. The factor Src is assumed to be a material parameter and can be used as a measure of the fracture toughness under mixed mode conditions. For pure mode I fracture, Src can be directly expressed in terms of KIc:
(50)Src=(κ-1)π2(KIc)28G,
where G denotes the shear modulus and κ is the kosolov constant. For general mixed mode fracture, Src is:
(51)Src=a11KI2+2a12KIKII+a22KII2∀θ=θ0,
with -π<θ0<π and
(52)a11=(1+cosθ)(κ-cosθ)16G,a12=(2cosθ-κ+1)sinθ16G,a22=(κ+1)(1-cosθ)+(1+cosθ)(3cosθ-1)16G.
The parameters aij are obtained from the strain-energy density function. For further details see Sih [290]. Note that the minimum strain energy-density criterion is a local criterion since fracture occurs when the energy density in a volume element near the crack tip reaches a critical value while the classical fracture mechanics theory is based on global energy balance.In the maximum energy release rate criterion, the crack propagates in the direction defined by the angleα=αc(σ),  where αc satisfies the condition:
(53)∂𝒢∂α|α=αc=0,∂2𝒢∂α2|α=αc≤0,
where α is defined in Figure 26 and 𝒢 denotes the energy release rate. Crack propagation occurs at a stress state σ=σc when 𝒢(σc,αc(σc))=𝒢g, where 𝒢g is a material parameter. According to Wu [291], the energy release rate can be decomposed into an antiplane (antiplane shear load) part 𝒢A and a plane (plane load) part 𝒢P:
(54)𝒢(σ,α)=𝒢P(σ,α)+𝒢A(σ,α).
For the antiplane shear load, 𝒢A is given by
(55)𝒢A(σ,α)=π2σ232(1-α1+αα).
For plane strain conditions:
(56)𝒢P(σ,α)=1-ν2σ232(KI2+KII2).Figure 26
Crack propagation with the energy release rate criterion.


### 5.2. Rankine Criterion

For a Rankine material, a crack is introduced when the principal tensile stress reaches the uniaxial tensile strength. The crack is initiated perpendicular to the direction of the principal tensile stress. Usually, some kind of smoothing technique is applied that either averages the crack normal or the stress tensor, [176, 228, 273, 294, 295]. This is done in order to improve the reliability of the computed stresses in front of the crack tip. For example, Wells and Sluys [294] and Mergheim et al. [176] use an averaged stress tensor of the form:
(57)σm=∫ΩcwσdΩ,
where Ωc is a certain domain around the crack tip and  w is a weighting function:
(58)w=1(2π)3/2l3exp(-r22l2),
where l determines how fast the weight function decays from the crack tip. Gasser and Holzapfel [273] and Rabczuk et al. [228] smooth the crack normal instead of the stress tensor. The Rankine criterion is applicable to brittle materials and works well for mode I-fracture.


### 5.3. Loss of Material Stability

Fracture is caused by a material instability. A classical definition of material stability is based on the so-called Legendre-Hadamard condition, which establishes that for any non zero vectorsn and h the following pointwise inequality must hold
(59)e=minn,h(n⊗h):A:(n⊗h)>0,
where A=Ct+σ⊗δ,Ct is the constitutive tangent operator and e is the so-called ellipticity indicator. For a linear comparison solid, if non-propagating singular surfaces do not occur, the Legendre-Hadamard condition is identical to the strong-ellipticity condition, used for example, in Simo et al. [58]. In case (59) is no longer fulfilled, the material loses stability and n defines the direction of propagation, and h is the polarization of the wave. This condition ensures that the speed of propagating waves in a solid remains real. Equality in expression (59) is the necessary condition for stationary waves. Belytschko et al. [31] give a textbook description for obtaining condition (59) by means of a stability analysis of the momentum equation when a perturbation of the form u=hexp(iωt+kn·x) is applied. The Legendre-Hadamard condition is occasionally called strong ellipticity of the constitutive relation (see Marsden and Hughes [296]). In the dynamical case, the Legendre-Hadamard condition implies the hyperbolicity of the IBVP. For a rate-dependent material, we refer to it as loss of material stability. The reader is referred to Šilhavý [25] or Ogden [297] for details about the concepts mentioned above.Loss of hyperbolicity (or loss of ellipticity or loss of material stability, resp.) is determined by minimizinge with respect to n and h; if e is negative for any combination of  n and h, the material has lost stability at that material point. In 2D, e can be expressed as a function of the two angles θ and ϕ, where n=[cos(θ)sin(θ)] and h=[cos(ϕ)sin(ϕ)]. Therefore, based on expression in (59), let us define for a given material point of a solid at a given time,
(60)Q=n·A·n,
where Q is the acoustic tensor with components Qik=njnlAijkl. We say that a material point is stable whenever the minimum eigenvalue of Q is strictly positive, and unstable otherwise. It can be shown that this condition is equivalent to condition (59).One difficulty is that the analysis of the acoustic tensor for isotropic materials generally will yield two directionsn from which one direction has to be chosen. At times the loss of material stability criterion becomes ambiguous; an example of the minimum eigenvector of Q as a function of the angle θ is shown in Figure 28; simplified for a two-dimensional example. For uniaxial tension, two angles are obtained, but they correspond to the same plane. In other stress states, four angles corresponding to two planes are sometimes obtained. A typical eigenvalue landscape as function of the crack angles in three dimension is shown in Figure 27 for one material point and isotropic J2-plasticity with strain softening.Figure 27
Minimum eigenvalue of the acoustic tensor as a function of the two angles for one material point.Minimum eigenvalue of the acoustic tensorQ versus different angles of the orientation of  n  for (a) uniaxial tension (b) an arbitrarily directed load.
(a)
(b)In [232], we chose the direction of the maximum displacement gradient by maximizing
(61)g=max︸l(nlT·(∇u·hl)),l=1,2,
where the normals nl correspond to minima of Q, (60); h is the corresponding eigenvector of Q. Oliver et al. [298] have shown that an anisotropic continuum model will lead to a unique angle θ in the material stability analysis.Since the normaln=[cosαcosφ,cosαsinφ,sinα] depends on two angles in three dimensions, the procedure of finding the minimum eigenvalue of Q can become computationally expensive. A bisection method can reduce computational cost. Another efficient way to compute the eigenvalues of the acoustic tensor is given by Ortiz et al. [77].


### 5.4. Rank-One Stability Criterion

Belytschko et al. [299] suggest to check the rank-one stability criterion
(62)F˙:A:F˙≤0∀F˙
before carrying out a material stability analysis at every quadrature point. They point out that it is a stricter but computational cheaper criterion than the loss-of material stability criterion. However, the rank-one stability criterion does not provide the orientation of the crack.


### 5.5. Energy Criteria

The orientation of the crack can be obtained by global energy minimization. For different orientations of the crack, the global energy is computed, and the crack is propagated in the direction where the global energy has its minimum. This global fracture criterion can be expressed as:(63)ℐ=∫Ω‍W(u)dΩ+∫Γ‍𝒢tdΓ,
where 𝒢t is the surface toughness that characterizes the surface energy necessary to create the new crack surface.Such criteria were used in the context of XFEM by Meschke and Dummerstorf [300], Dummerstorf and Meschke [301]. It was presented in a two-dimensional setting for cohesive and cohesionless cracks. As pointed out by the authors, an extension to 3D is cumbersome. Miehe and Gürses [289] presented a similar method in the context of an r-adaptive interelement-separation method in two dimension. The crack orientation as well as the crack length is introduced as additional unknowns in the variational formulation.


### 5.6. Computation of the Crack “Length”

As mentioned above, the easiest way is to control the “length” of the crack. In the context of XFEM, Belytschko et al. [302] suggested to take advantage of the level set method to obtain the new position of the crack tip. The key assumption is that the hyperbolicity indicator e=h·Q·h must vanish at the crack tip:
(64)∂e∂t+vc·∇e=0withvc=vcs,
where vc is the crack speed, s  gives its direction which must fulfill the condition n·s=0. To obtain the crack length, (64) has to be solved for vc. Due to the hyperbolic character of the differential equation, the standard Bubnov Galerkin formulation leads to instabilities and stabilization techniques such as upwinding, Petrov-Galerkin, Galerkin least squares, streamline upwind Petrov Galerkin, or finite increment calculus Stabilization is needed as pointed out for example, in [285].


## 5.1. Fracture Mechanics-Based Criteria

Besides approaches based on configurational forces [286–289] there are four major cracking criteria in LEFM.(i)
Maximum hoop stress criterion or maximum principal stress criterion.(ii)
Minimum strain energy density criterion, Sih [290].(iii)
Maximum energy release rate criterion, Wu [291].(iv)
The zeroKII criterion (Vanishing in-plane SIF (KII) in shear mode for infinitesimally small crack extension), Goldstein and Salganik [292].The first two criteria predict the direction of the crack trajectory from the stress state prior to the crack extension. The last two criteria require stress analysis for virtually extended cracks in various directions to find the appropriate crack-growth direction. Duflot and Hung [223] report for medium mixed-mode problems, all criteria yield almost identical results. However, according to Shen and Stephansson [293], only the maximum energy release rate criterion allows to predict secondary cracks in compressed specimens. The crack is propagated in an angle of θc from the crack tip. Note that LEFM can only deal with crack propagation but not with crack nucleation. It also can not handle crack branching.In the maximum hoop stress or maximum principal stress criterion, the maximum circumferential stressσθθ, often called hoop stress, in the polar coordinate system around the crack tip corresponds to the maximum principal stress and is given for a crack propagating with constant velocity vc by:
(46)σθθ=KI2πrfhI(θ,vc)+KII2πrfhII(θ,vc),
where the functions fhI and fhII represent the angular variation of the stress for different values of the crack-tip speed vc. When the maximum hoop stress is larger equal a critical hoop stress σθθc then the crack is propagated in the direction perpendicular to the maximum hoop stress. For pure mode I fracture, σθθc is given by
(47)σθθc=KIc2πr,
with the fracture toughness KIc.In LEFM, the local direction of the crack growth is determined by the condition that the local shear stress is zero that leads to the condition:(48)KIsinθc+KII(3cosθc-1)=0,
that results in the crack propagation angle
(49)θc=2arctan(KI-KI2+8KII24KII).The minimum strain energy density criterion is based on a critical strain-energy-density factorSrc. The basic assumption is that fracture occurs when the interior minimum of the strain-energy-density factor S reached Src. The strain energy density factor S represents the strength of the elastic energy field in the vicinity of the crack tip which is singular of the order r-1. The factor Src is assumed to be a material parameter and can be used as a measure of the fracture toughness under mixed mode conditions. For pure mode I fracture, Src can be directly expressed in terms of KIc:
(50)Src=(κ-1)π2(KIc)28G,
where G denotes the shear modulus and κ is the kosolov constant. For general mixed mode fracture, Src is:
(51)Src=a11KI2+2a12KIKII+a22KII2∀θ=θ0,
with -π<θ0<π and
(52)a11=(1+cosθ)(κ-cosθ)16G,a12=(2cosθ-κ+1)sinθ16G,a22=(κ+1)(1-cosθ)+(1+cosθ)(3cosθ-1)16G.
The parameters aij are obtained from the strain-energy density function. For further details see Sih [290]. Note that the minimum strain energy-density criterion is a local criterion since fracture occurs when the energy density in a volume element near the crack tip reaches a critical value while the classical fracture mechanics theory is based on global energy balance.In the maximum energy release rate criterion, the crack propagates in the direction defined by the angleα=αc(σ),  where αc satisfies the condition:
(53)∂𝒢∂α|α=αc=0,∂2𝒢∂α2|α=αc≤0,
where α is defined in Figure 26 and 𝒢 denotes the energy release rate. Crack propagation occurs at a stress state σ=σc when 𝒢(σc,αc(σc))=𝒢g, where 𝒢g is a material parameter. According to Wu [291], the energy release rate can be decomposed into an antiplane (antiplane shear load) part 𝒢A and a plane (plane load) part 𝒢P:
(54)𝒢(σ,α)=𝒢P(σ,α)+𝒢A(σ,α).
For the antiplane shear load, 𝒢A is given by
(55)𝒢A(σ,α)=π2σ232(1-α1+αα).
For plane strain conditions:
(56)𝒢P(σ,α)=1-ν2σ232(KI2+KII2).Figure 26
Crack propagation with the energy release rate criterion.


## 5.2. Rankine Criterion

For a Rankine material, a crack is introduced when the principal tensile stress reaches the uniaxial tensile strength. The crack is initiated perpendicular to the direction of the principal tensile stress. Usually, some kind of smoothing technique is applied that either averages the crack normal or the stress tensor, [176, 228, 273, 294, 295]. This is done in order to improve the reliability of the computed stresses in front of the crack tip. For example, Wells and Sluys [294] and Mergheim et al. [176] use an averaged stress tensor of the form:
(57)σm=∫ΩcwσdΩ,
where Ωc is a certain domain around the crack tip and  w is a weighting function:
(58)w=1(2π)3/2l3exp(-r22l2),
where l determines how fast the weight function decays from the crack tip. Gasser and Holzapfel [273] and Rabczuk et al. [228] smooth the crack normal instead of the stress tensor. The Rankine criterion is applicable to brittle materials and works well for mode I-fracture.


## 5.3. Loss of Material Stability

Fracture is caused by a material instability. A classical definition of material stability is based on the so-called Legendre-Hadamard condition, which establishes that for any non zero vectorsn and h the following pointwise inequality must hold
(59)e=minn,h(n⊗h):A:(n⊗h)>0,
where A=Ct+σ⊗δ,Ct is the constitutive tangent operator and e is the so-called ellipticity indicator. For a linear comparison solid, if non-propagating singular surfaces do not occur, the Legendre-Hadamard condition is identical to the strong-ellipticity condition, used for example, in Simo et al. [58]. In case (59) is no longer fulfilled, the material loses stability and n defines the direction of propagation, and h is the polarization of the wave. This condition ensures that the speed of propagating waves in a solid remains real. Equality in expression (59) is the necessary condition for stationary waves. Belytschko et al. [31] give a textbook description for obtaining condition (59) by means of a stability analysis of the momentum equation when a perturbation of the form u=hexp(iωt+kn·x) is applied. The Legendre-Hadamard condition is occasionally called strong ellipticity of the constitutive relation (see Marsden and Hughes [296]). In the dynamical case, the Legendre-Hadamard condition implies the hyperbolicity of the IBVP. For a rate-dependent material, we refer to it as loss of material stability. The reader is referred to Šilhavý [25] or Ogden [297] for details about the concepts mentioned above.Loss of hyperbolicity (or loss of ellipticity or loss of material stability, resp.) is determined by minimizinge with respect to n and h; if e is negative for any combination of  n and h, the material has lost stability at that material point. In 2D, e can be expressed as a function of the two angles θ and ϕ, where n=[cos(θ)sin(θ)] and h=[cos(ϕ)sin(ϕ)]. Therefore, based on expression in (59), let us define for a given material point of a solid at a given time,
(60)Q=n·A·n,
where Q is the acoustic tensor with components Qik=njnlAijkl. We say that a material point is stable whenever the minimum eigenvalue of Q is strictly positive, and unstable otherwise. It can be shown that this condition is equivalent to condition (59).One difficulty is that the analysis of the acoustic tensor for isotropic materials generally will yield two directionsn from which one direction has to be chosen. At times the loss of material stability criterion becomes ambiguous; an example of the minimum eigenvector of Q as a function of the angle θ is shown in Figure 28; simplified for a two-dimensional example. For uniaxial tension, two angles are obtained, but they correspond to the same plane. In other stress states, four angles corresponding to two planes are sometimes obtained. A typical eigenvalue landscape as function of the crack angles in three dimension is shown in Figure 27 for one material point and isotropic J2-plasticity with strain softening.Figure 27
Minimum eigenvalue of the acoustic tensor as a function of the two angles for one material point.Minimum eigenvalue of the acoustic tensorQ versus different angles of the orientation of  n  for (a) uniaxial tension (b) an arbitrarily directed load.
(a)
(b)In [232], we chose the direction of the maximum displacement gradient by maximizing
(61)g=max︸l(nlT·(∇u·hl)),l=1,2,
where the normals nl correspond to minima of Q, (60); h is the corresponding eigenvector of Q. Oliver et al. [298] have shown that an anisotropic continuum model will lead to a unique angle θ in the material stability analysis.Since the normaln=[cosαcosφ,cosαsinφ,sinα] depends on two angles in three dimensions, the procedure of finding the minimum eigenvalue of Q can become computationally expensive. A bisection method can reduce computational cost. Another efficient way to compute the eigenvalues of the acoustic tensor is given by Ortiz et al. [77].


## 5.4. Rank-One Stability Criterion

Belytschko et al. [299] suggest to check the rank-one stability criterion
(62)F˙:A:F˙≤0∀F˙
before carrying out a material stability analysis at every quadrature point. They point out that it is a stricter but computational cheaper criterion than the loss-of material stability criterion. However, the rank-one stability criterion does not provide the orientation of the crack.


## 5.5. Energy Criteria

The orientation of the crack can be obtained by global energy minimization. For different orientations of the crack, the global energy is computed, and the crack is propagated in the direction where the global energy has its minimum. This global fracture criterion can be expressed as:(63)ℐ=∫Ω‍W(u)dΩ+∫Γ‍𝒢tdΓ,
where 𝒢t is the surface toughness that characterizes the surface energy necessary to create the new crack surface.Such criteria were used in the context of XFEM by Meschke and Dummerstorf [300], Dummerstorf and Meschke [301]. It was presented in a two-dimensional setting for cohesive and cohesionless cracks. As pointed out by the authors, an extension to 3D is cumbersome. Miehe and Gürses [289] presented a similar method in the context of an r-adaptive interelement-separation method in two dimension. The crack orientation as well as the crack length is introduced as additional unknowns in the variational formulation.


## 5.6. Computation of the Crack “Length”

As mentioned above, the easiest way is to control the “length” of the crack. In the context of XFEM, Belytschko et al. [302] suggested to take advantage of the level set method to obtain the new position of the crack tip. The key assumption is that the hyperbolicity indicator e=h·Q·h must vanish at the crack tip:
(64)∂e∂t+vc·∇e=0withvc=vcs,
where vc is the crack speed, s  gives its direction which must fulfill the condition n·s=0. To obtain the crack length, (64) has to be solved for vc. Due to the hyperbolic character of the differential equation, the standard Bubnov Galerkin formulation leads to instabilities and stabilization techniques such as upwinding, Petrov-Galerkin, Galerkin least squares, streamline upwind Petrov Galerkin, or finite increment calculus Stabilization is needed as pointed out for example, in [285].


## 6. Future Perspectives and Conclusion

Numerous computational methods for fracture have been developed in the past two decades. Advances in partition of unity methods have provided effective and reliable tools to analyze fracture for numerous applications. Improvements were made particularly concerning the accuracy (through enrichment) and the efficiency (no remeshing). Some of those methods are already available in commercial CAE software packages and can be used for commercial applications as pointed out previously. Most of the partition of unity methods for fracture are well suited for static fracture when a moderate number of cracks occur. Many applications (of partition of unity methods) focus on crack propagation, mainly for propagation of single crack surfaces without branching cracks and crack coalescence. While the majority of the computational methods discussed above are capable of handling complex phenomena such as branching cracks, reliable criteria to branch a crack in practical FE-simulations are still missing. Those fracture criteria can be categorized in local criteria and global criteria. Local criteria are based on the stress state in the vicinity of the crack tip while global criteria are based on energy minimization. Both criteria should in principle predict complex fracture pattern naturally. For example, it is known that the ellipticity indicatore should be zero at the crack tip, and if one could find the locations where e=0, the positions of the new crack tip are obtained naturally. Besides some efforts to seaming-less model the transition from continuum to discontinuum, the combination of the fracture criterion and the computational method remains a challenge. A very interesting approach is the variational approach to fracture [262] as the crack path is the natural outcome of the analysis. It also circumvents the implementation of complex crack tracking algorithms and the need to describe the topology of the crack surface. Moreover, the analysis of coupled problems seems to be simpler as the enrichment in XFEM can become quite cumbersome. However, most implementations of the (already computational expensive) variational approach to fracture require a high resolution of the crack (with several elements) increasing computational cost further. Combining the advantages of partition of unity methods for cracks (such as XFEM) and the variational approach to fracture does not seem to be easy, but it could lead to a new efficient and reliable pathway to model complex fracture patterns.While—as stated above—the majority of the applications of computational methods for fracture focus on crack propagation problems, far less studies are concerned with crack nucleation. Methods as XFEM seem not well suited for dynamic fracture and fragmentation involving the nucleation, growth, and coalescence of an enormous number of cracks. An alternative pathway to dynamic fracture offer methods such as the distinct element method (DEM) [303, 304], discontinuous deformation analysis (DDA) [305, 306], the particle flow code (PFC) [307] or peridynamics [308]; in brief discrete methods (not discussed in this paper) that are not based on continuum approaches.The exact prediction of complex fracture patterns is nearly impossible for many applications due to its stochastic nature. The onset of crack nucleation and in fact also the direction in crack propagation depends on numerous factors such as the micro-structure of the material (imperfections, voids, microcracks, etc.) or loading conditions, that are not known exactly. The author has made the experience that results are particularly sensitive with respect to boundary conditions for example. Most of the novel computational methods for fracture are based on deterministic approaches. There are far less contributions on statistical computational methods for fracture [309–314]. The development of efficient and reliable stochastic computational methods is one of the key challenges in the future. Two pathways can be followed: (i)
deterministic methods in a stochastic setting and(ii)
full stochastic methods.The first approach seems simpler as no modifications of the existing computational methods are needed. Sampling methods belong to the first category while variance reduction methods for example belong to the second category. One major challenge will be the design ofcomputational efficient methods.The key objective of computational methods is their application to “real-world” problems, for example, in order to support the design of new products. The choice of the method depends on the application. Most effort in the past was devoted to develop new methods. However, there has been minimal effort to assess those methods. In view of the growing number of computational methods for fracture, it would be of high practical relevance to provide a framework to quantitatively assess the quality of existing computational methods with respect to their accuracy, reliability, and robustness for a specific application. For example, computational methods (and constitutive models) have a number of uncertain input parameters such as the dilation parameter in meshfree methods. It is of utmost importance to quantitatively determine the sensitivity of the uncertain input parameters with respect to a quantity of interest (the output). It will also help to quantify the predictive capabilities of the computational methods.One major application of computational methods for fracture is their use in the design of new materials. Obtaining a fundamental understanding of how materials fail was/is a main research direction in materials science. In failure, the response of a structure is driven by fine scale features (nano- or micro-structure of the material). For an accurate prediction of material failure, it is important to account for fine scale features in the process zone. Therefore, an important future research direction is the development of multiscale methods for fracture. While numerous multiscale methods (see e.g., [316, 317]) were developed for intact materials, far fewer methods are applicable for fracture simulations though the effort is rapidly increasing.Multiscale methods can be categorized into hierarchical, semiconcurrent, and concurrent methods [318]. In hierarchical multiscale methods, information is passed from the fine-scale to the coarse scale but not vice versa. Computational homogenization [319] is a classical hierarchical upscaling technique. However, hierarchical multiscale approaches seem less suitable to model fracture. One basic assumption for the application of homogenization theories is the existence of disparate length scales [320]: ℒCr≪ℒRVE≪ℒSpec, where ℒCr,  ℒRVE and  ℒSpec are the crack length, the representative volume element (RVE) and specimen-size, respectively. For problems involving fracture, the first condition is violated as ℒCr is of the order of ℒRVE. Moreover, periodic boundary conditions (PBCs) often used at the fine-scale, cannot be used when a crack touches a boundary as the displacement jump in that boundary violates the PBC. Nevertheless, fine-scale features can also be incorporated “hierarchically” through enrichment. For example, Gracie et al. [321, 322], Belytschko and Gracie [323] developed special enrichment functions for dislocations.The basic idea of semiconcurrent multiscale methods is illustrated in Figure29(b). In semi-concurrent multiscale methods, information is passed from the fine-scale to the coarse-scale and vice versa. Semi-concurrent multiscale methods are of the same computational efficiency as concurrent multiscale methods [324]. The key advantage of semi-concurrent multiscale methods over concurrent multiscale methods is their flexibility, that is, their ability to couple two different software packages, for example, MD software to FE software. Parallelization is generally simple as well. A popular semi-concurrent multiscale method is the FE2 [325] originally developed for intact materials. Kouznetsova [319] was the first who extended this method to problems involving material failure, see also Kouznetsova et al. [326] or recent contribution by Nguyen et al. [327], Verhoosel et al. [328], and Belytschko et al. [299].Schematic of a (a) hierarchical, (b) semiconcurrent, and (c) concurrent multiscale methods [315].
(a)
(b)
(c)Numerous concurrent multiscale methods [329–332] have been developed that can be classified into “interface” coupling methods and “Handshake” coupling methods. Interface coupling methods seem to be less effective for dynamic applications as avoiding spurious wave reflections at the “artificial” interface seem to be more problematic. Some of the concurrent multiscale methods have been extended to modeling fracture [318, 333–335].Future challenges will lie in the development of efficient methods to transfer length scales when fracture occurs, to adaptively choose the discretization and the model based on error estimation, and to bridge disparate time scale. To overcome the high computational cost will be another challenge, in particular, when these methods will be applied to “real-world” applications, for example, in computational materials design.



---

*Source: 849231-2013-03-20.xml*