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.

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 Figure 1. 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 with respect to the smallest critical dimension of the structure is important for the choice of the fracture model [8]. For , LEFM is valid while for , 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 . The length of the fracture process zone is of the order of the characteristic length [11]: where is the modulus of elasticity, is the Poisson’s ratio, is the fracture energy and denotes the tensile strength of the material. Carpinteri [12] introduced a nondimensional brittleness number where 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 . If the process zone is small compared to , 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 Figure 3. 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 () and the dimension of the mass density (), there is indeed a length scale associated with the viscosity given by [14]: where is the longitudinal propagation velocity in 1D. However, the well-posedness of the IBVP is only guaranteed during the time span of the viscosity (). 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 Figure 4. 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), Figure 5(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].

(a)

(b)

(c)

(d)

(e)

(f)

(a)
(b)
(c)
(d)
(e)
(f)

Figure 5 

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.

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 energy . In the one-dimensional case of pure mode I-failure, the fracture energy is the energy that is dissipated during crack opening. It is related to the area of the cohesive crack; is the cohesive traction, is the crack opening and 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 the crack 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].

(a)

(b)

(c)

(a)
(b)
(c)

Figure 6 

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].

The smeared crack models can be further classified into rotating and fixed crack models. In the fixed 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).

(a)

(b)

(a)
(b)

Figure 7 

Principle of (a) rotating crack models, (b) fixed crack models.

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 Figure 8; 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.

(a)

(b)

(a)
(b)

Figure 8 

Mesh orientation bias for smeared crack models.

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-called interelement-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 -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.

(a)

(b)

(c)

(a)
(b)
(c)

Figure 9 

Element with (a) one weak discontinuity, (b) two weak discontinuities, and (c) one strong discontinuity.

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] where are the nodal parameters of the usual approximation, is the enrichment, and the superimposed “” indicates the enrichment on element level. The function is given by where is the step function acting on the crack line and is the number of nodes of element that belong to the domain ; see Figure 9(c). The discretization of the strain field can be written as where the superimposed denotes the symmetric part of a tensor, is a regularized Dirac delta function, and is a collocation function defined as with the thickness 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 with where is the elasticity tensor, is the B operator containing the derivatives of the shape functions and is the B operator of the enrichment that depends on the embedded element formulation (SOS and KOS or KSON). For the KSON formulation, , that will result in a non-symmetric stiffness matrix. For the SOS and KOS formulation, . The elemental enrichment can be condensed on the element level: that leads with the expression for to the system of equations with 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)

(b)

(a)
(b)

Figure 10 

(a) Piecewise constant crack opening in embedded elements and (b) linear crack opening for XFEM.

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 part and a discontinuous part : 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 (). 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 cracks with crack tips reads where is the set of nodes in the entire discretization, is the set of nodes around the crack tip, and is the set of nodes associated to elements completely cut by the crack; is the enrichment function of crack ; is the enrichment function for the crack tip ; and 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 with where denotes the normal vector to the crack surface and is a level set function; see Figure 12. The zero-isobar 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 vanishes in the neighboring element. The shifting also maintains the interpolatory character of the XFEM approximation or in other words: the nodal parameters remain the physical displacements.

Figure 12 

Representation of the crack surface by level sets and .

In LEFM, the crack tip enrichment is choosen according to the analytical solution: where 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: with . The sin-term ensures the discontinuity in the displacement field around the crack tip while 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)

(b)

(c)

(a)
(b)
(c)

Figure 13 

(a) Topological enrichment versus (b) and (c) geometrical enrichment.

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 is and , see Figure 12; is a level set orthogonal to (i.e., ) 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 where commonly the shape functions 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, Figure 15. 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 , blending domain , and the total (standard) domain .

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 or with the stiffness matrix where the superscripts , and indicate the “standard”, “step-enriched” and “tip-enriched” part of the discretization; is the vector with the nodal parameters, is the external force vector with and with body forces and applied tractions . The stiffness matrix is where , is the B operator defined by: 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 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 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 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)

(b)

(c) effective crack length

(d)

(a)
(b)
(c) effective crack length
(d)

Figure 16 

(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.

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” terms and (or and 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 Figure 18 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 is extended to the phantom domain . Then, the displacement in the real domain can be interpolated by using the degrees of freedom for the nodes in the phantom domain . 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]: where is the signed distance measured from the crack, , , and are nodes belonging to , , and , respectively. 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, and . 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: where denotes the meshfree shape function and 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 by . Usually, radial kernel functions are employed where is expressed in terms of . The kernel functions are chosen to have compact support, that is, if , where is the “interpolation" radius often called dilation parameter. The kernel function can be expressed in terms of material coordinates (Lagrangian kernel) or spatial coordinates (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: where is the polynomial basis. Usually, a linear polynomial basis is choosen, given in 2D by

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 Figure 19: 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 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)

(b)

(a)
(b)

Figure 20 

(a) Undesired introduced discontinuities by the visibility method, (b) Difficulties with the visibility method for concave boundaries and kinks.

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 Figure 20 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 the undesired 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 .

(a)

(b)

(a)
(b)

Figure 21 

The diffraction method.

The idea of the diffraction method is not only to treat the crack as opaque but also to evaluate the length of the ray 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 with ; they suggested .

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 with as a cutoff value, the boundary can be considered as convex, otherwise non-convex. The cutoff value of corresponds to the wedge angle of in Figure 22.

Figure 22 

Domain of influence near a wedge-shaped non-convex boundary. The boundary is enforced if .

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: 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: where are the material coordinates, denotes the continuous displacement, and denotes the discontinuous part.

The crack is modelled by a set of discrete cracks as shown in Figure 23. 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.

(a)

(b)

(a)
(b)

Figure 23 

Schematic on the right shows a crack model for the crack on the left.

The approximation of the displacement field in the cracking particles method is given by where is the total set of nodes in the model and is the set of cracked nodes; 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 field and hence the enrichment parameters . It can be shown from (33) that 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 , ) 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 , the crack prevention zone is given through the angle : For branching cracks, it is often observed that the orientation of neighboring cracks severely differ, see node in Figure 25. Therefore, cracks are also allowed in the prevention zone when the following criterion is fulfilled: with and . 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)

(b)

(a)
(b)

Figure 24 

(a) Spurious cracking and (b) improved crack pattern.

Figure 25 

Crack prevention: nodes and are located in the crack prevention zone; node 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 as , where is a postulated quantity known as the surface (dissipation) energy of the crack system. Thermodynamic equilibrium with respect to crack growth then requires that , where 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 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 may be obtained, equal to for the case of plane stress. We now assume that the surface energy to be , 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: . According to the classical theory, the existing crack will not propagate for and will do so for .

Despite improvements made by subsequent authors to account for nonlinearity and inelasticity, the original form of the fracture criterion remained, that is, . 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 (), 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 (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 where denotes the -D Hausdorff measure. The function represents the energy required to create an infinitesimal crack at the point 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 -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 in which is the set of kinematically admissible displacement fields. The total energy of the body is the sum of these two terms: The problem then is to find the mapping such that at some given time , the crack geometry is the one that minimizes and furthermore contains all for (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, where is a function obeying the constraints , , and . 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 where the positive and negative parts of the energies are defined by and are ramp functions defined as . 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 are , , and . To maintain numerical stability, a small positive constant 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 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 where and are two vector fields, and denotes the vector of the normal to the crack surface. The family of surfaces, enveloping both vector fields and are described by a temperature-like function when the following conditions hold 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.