A phase-field formulation for dynamic cohesive fracture

A phase-field formulation for dynamic cohesive fracture


We extend a phase-field/gradient damage formulation for cohesive fracture to the dynamic case. The model is characterized by a regularized fracture energy that is linear in the damage field, as well as non-polynomial degradation functions. Two categories of degradation functions are examined, and a process to derive a given degradation function based on a local stress-strain response in the cohesive zone is presented. The resulting model is characterized by a linear elastic regime prior to the onset of damage, and controlled strain-softening thereafter. The governing equations are derived according to macro- and microforce balance theories, naturally accounting for the irreversible nature of the fracture process by introducing suitable constraints for the kinetics of the underlying microstructural changes. The model is complemented by an efficient staggered solution scheme based on an augmented Lagrangian method. Numerical examples demonstrate that the proposed model is a robust and effective method for simulating cohesive crack propagation, with particular emphasis on dynamic fracture.

phase-field models; cohesive fracture; dynamic fracture; gradient models; damage

1 Introduction

In recent years diffuse crack approaches have become increasingly popular for predicting crack propagation in engineering materials. This is largely due to their potential for capturing the evolution of complex crack patterns without the need for algorithmic crack-front tracking methods. There has been a proliferation of modeling techniques developed around a variational formulation of, most commonly, a Griffith-type description of fracture. In this manuscript, we take an alternative approach and develop a phase-field model of dynamic fracture that approximates a cohesive type of response. Critically, our model extends the work of Lorentz et al. (2011, 2012) to the dynamic regime and to a broad set of local cohesive behaviors. In contrast to regularized formulations based on Griffith-type descriptions of fracture, the effective macroscopic fracture parameters can be held fixed as the regularization length scale is decreased. We employ Lagrange multipliers to restrict the evolution of the damage to satisfy physically-motivated constraints. The robustness of the resulting formulation is demonstrated for multi-dimensional quasi-static and dynamic problems in fracture and fragmentation.

Phase-field models for fracture rely on the regularization of a sharp crack surface by means of an elliptical functional. In Ambrosio and Tortorelli (1990), a phase-field approximation to discontinuous fields was proposed, inspired by the work on image segmentation by Mumford and Shah (1989). In Francfort and Marigo (1998) and Bourdin et al. (2000, 2008) a phase-field approximation was proposed to a variational formulation of brittle fracture based on the same potential. In these variational frameworks, which have largely become the standard in the phase-field for fracture community, sharp crack topologies are regularized by their diffusive counterparts and governed by a scalar continuous field variable. More recently, alternative derivations have been proposed in Miehe et al. (2010a, b) and da Silva et al. (2013), which have broadened the scope of phase-field modeling for fracture. With respect to dynamic fracture, phase-field models have been considered in Bourdin et al. (2011), Borden et al. (2012), Hofacker and Miehe (2013), Schlüter et al. (2014) and Li et al. (2016), just to name a few. The vast majority of these works have been based on a regularization of a Griffith type description of fracture, in which the stress field at the crack tip is singular.

In fact, it bears emphasis that most materials are not perfectly brittle in the Griffith sense, but display some ductility after reaching their ultimate strength. Often, there exists a process zone in the vicinity of the crack tip. This region is typically characterized by small-scale yielding, micro-crack initiation, growth and coalescence, as graphically depicted in Figure 1. In those cases in which the fracture process zone is sufficiently small compared to the representative size of the structure, the use of a Griffith description of fracture is justified. In the more general case, cohesive forces in the fracture process zone must be taken into account, and the utility of a Griffith model diminishes.


[width=.46]fpz-eps-converted-to.pdf Traction-freeFracture process zoneIntact

Figure 1: Schematic representation of a fracture process zone.

In addition to this restriction, there are clear disadvantages to phase-field models based on a Griffith description of fracture. It is well known that the regularization length is intrinsically tied to the material properties, e.g. see Borden et al. (2012), which poses a significant limitation for many important crack propagation problems. In addition, such formulations typically exhibit nonlinear behavior immediately at the onset of loading. In the quest to develop a remedy to these issues, researchers have adopted a number of different perspectives. In Pham and Marigo (2013) a modification to the phase-field potential was proposed which results in a purely elastic response up to the onset of damage. In Miehe et al. (2015a, b) the same modification was considered and complemented by a generalization of the crack driving forces, allowing for a large spectrum of multi-field problems to be studied in the phase-field framework. In particular, it is characterized by a threshold energy in the elastic regime, after which stress-based failure criteria are defined to govern crack growth in the strain-softening regime. The same motivation lead Borden et al. (2016) to employ a cubic degradation function instead of the commonly adopted quadratic functional. This function essentially provides a stress-strain response prior to crack growth that more closely approximates linear elasticity. Despite the quasi-linear behavior, however, there is no means of calibrating the effective cohesive response to different materials. Therefore we view this approach as somewhat limited in use in the context of cohesive fracture problems. In Tanné et al. (2018) it was shown that both crack nucleation and propagation can be predicted for a large variety of materials by carefully choosing the length scale in a standard phase-field formulation. However, for materials with large process zone sizes relative to the geometry, this strategy results in regularized damage bands that are extremely wide.

An extension of the variational formulation of brittle fracture to cohesive fracture has been considered in Bourdin et al. (2008), but the development of a phase-field approximation is widely regarded as non-trivial. Several recent efforts have been made to establish such a formulation. In Lorentz and Godard (2011) a gradient damage model was proposed that decouples the regularization length from the macroscopic fracture parameters. The same formulation was demonstrated to converge to a cohesive zone model in the limit of vanishing regularization length in Lorentz et al. (2011). However, its extension to the multidimensional case remains, from a conceptual standpoint, a non-trivial challenge. In Verhoosel and de Borst (2013), an attempt was made to construct a phase-field model for cohesive fracture by casting the cohesive zone approach, pioneered by Dugdale (1960) and Barenblatt (1962), in an energetic framework. The approach relies on the construction of an auxiliary field to capture the displacement jump across the regularized fracture surface.

With particular regard to the gradient damage model from Lorentz and Godard (2011), we caution against establishing too strong a connection between regularized descriptions of cohesive fracture and the cohesive zone concept. In a conventional cohesive zone model, tractions are transmitted across a two-dimensional surface, which is embedded in a three-dimensional continuum. The relevant kinematic quantities are the crack opening displacements in the normal and shear directions. However, a kinematic quantity that represents the stretching of the fracture plane, i.e. the strain component parallel to the crack, is lacking in classical cohesive zone models. Generally speaking, these in-plane components cannot be ignored. For this reason, the cohesive zone model should be regarded as a strictly uniaxial model, as argued in Bažant (2002). In fact, a number of approaches have recently been proposed to overcome this limitation of cohesive zone models, such as the finite band method from Huespe et al. (2009) and the cohesive band model from Remmers et al. (2013), just to mention a few. It is well known that conventional zero thickness cohesive zone models are very limited in use, despite their conceptual simplicity. The aforementioned works make an attempt to include stress triaxiality effects in the constitutive response, which play an important role in the context of the ductile failure of materials. While the approach proposed in this paper embraces a diffuse description of fracture, the formulation is seen to be more in line with these contributions than with conventional cohesive zone formulations.

In this contribution, we explore the fundamental properties of dynamic cohesive fracture in a phase-field setting. The particular choice of the crack surface density functional will be motivated and various degradation functions are derived and examined. We examine the issue of convergence of these types of models as the regularization length scale is reduced, while maintaining sufficient numerical resolution. While some evidence of mesh insensitivity has previously been demonstrated in Lorentz et al. (2012), a formal numerical validation and convergence study has yet to be provided. Finally, we examine various options for enforcing the irreversibility of the damage field as fracture progresses. Several benchmark problems in quasi-static and dynamic fracture are considered to demonstrate the robustness and efficacy of the models and methods.

The paper is structured as follows. In Section 2 a phase-field/gradient damage formulation is presented that is particularly well-suited for cohesive fracture mechanics. Special attention is given to the elastic and dissipation functionals, after which the governing equations are derived using macro- and microforce balance theories. The section concludes with a discussion on the crack driving forces and the stiffness degradation function. Section 3 provides some analysis for the model described in Section 2. Some of the fundamental properties of the method are illustrated by means of a one-dimensional analysis. We then present a process for deriving degradation functions that give rise to particular stress-strain behaviors in the process zone, with an emphasis on linear decay. Finally, we discuss some aspects concerning the regularization of crack topologies and the issue of -convergence. Section 4 discusses the finite element implementation of the proposed formulation. The augmented Lagrangian method for enforcing irreversibility is reviewed and the coupled problem is formulated by means of the Galerkin method. The numerical formulation is completed by the definition of a robust solution scheme based on operator splits. Section 5 provides a series of numerical experiments under quasi-static and dynamic loading conditions. Finally, a summary and concluding remarks are presented in Section 6.

2 Cohesive fracture in a phase-field setting

We now define a phase-field/gradient damage formulation for dynamic cohesive fracture in elastic solids under small-strains and iso-thermal conditions. The theory is set up around materials for which a Griffith’s description of fracture is no longer justified due to the size of the fracture process zone. The governing equations are obtained following macro- and microforce balance theories, which have been commonly applied in the development of phase-field theories. The theory developed here incorporates irreversibility by introducing suitable constraints on the fracture evolution.

2.1 General considerations


[width=.85]domain.pdf (a)(b)

Figure 2: (a) Sketch of a body, , with an internal discontinuity . (b) A regularized representation of the internal discontinuity.

Consider a body (with ) with external boundary with , and an internal discontinuity boundary as shown in Figure 2a. The state of the system is described by two independent variables, the vector displacement field and a scalar damage field . As in standard phase-field descriptions, the damage plays the role of approximating a given crack by its regularized counterpart , as shown in Figure 2b. Following continuum damage mechanics conventions, it takes values in , with away from the crack surface and inside the crack. Small deformations and deformation gradients are assumed. The infinitesimal strain tensor is defined as


We assume that the damage field acts only to degrade the tensile resistance of the body and that crack propagation is prohibited under compression. Following Miehe et al. (2010a), this is effected by employing a spectral decomposition of into positive and negative components, via


where are the principal strains and the principal strain directions. The positive and negative operators and are defined in accord with


The decomposition of the strain makes it possible to decompose the strain energy density into tensile and compressive contributions. We restrict attention to isotropic linear elasticity, and define these strain energy densities by


where and are the Lamé coefficients.

We assume that the total potential energy of a body consists of bulk and fracture contributions. In accordance with an energetic description of fractured bodies, its total free energy can be written as


where was introduced as the critical fracture energy per unit area and equals the unknown fracture surface. As in Bourdin et al. (2008), we approximate this functional by


where is the degrading elastic bulk energy and the crack surface density functional. To account for a tension-compression asymmetry, the elastic bulk energy is frequently decomposed as


where is a stiffness degradation function. The degradation function is assumed to be a monotonically decreasing function of damage, with (initial damage) and (complete loss of stiffness). Explicit forms of this function are provided in Section 2.3.

With respect to the energetic contribution from the crack surface, we propose the use of the following crack surface density functional:


in which the regularization length scale is introduced. We note that, in contrast to the commonly employed phase-field approximiation from Bourdin et al. (2008), the damage field enters the energy through a linear term. Similar formulations have been discussed, amongst others, in Frémond and Nedjar (1996), Pham et al. (2011), Miehe (2011) and Lorentz and Godard (2011); Lorentz et al. (2012), in the context of gradient damage mechanics. It will be demonstrated that such an approximation is central in the development of a phase-field description of cohesive fracture.

2.2 Evolution equations

We suppose the microstructural changes, governed by the physics of the underlying problem, act to enforce a crack irreversibility condition that can be expressed as the inequality


where the superposed dot denotes a differentiation with respect to time.

We follow the approach outlined in da Silva et al. (2013) based on the existence of a set of internal constraints. This leads to the definition of a macroscopic momentum balance


and a microscopic force balance


where denotes a combined strain energy density and fracture energy density, given by


and is a kinetic modulus. The right-hand side of (11) is essential as it embodies the irreversibility of damage evolution. It takes the form and is complemented by a reactive microforce to enforce that constraint. In particular, vanishes for and is given by the left-hand side of (11) for . Given a positive kinetic modulus , the damage only increases if the left-hand side of (11) is positive.

We now turn our focus to the kinetic modulus . In Miehe et al. (2010a) a similar term was employed, which was coined the viscous regularization, with the purpose of stabilizing the numerical treatment of the algebraic equations. While it is clearly demonstrated how such a regularization has an impact, it remains unclear, generally speaking, how to determine such a parameter in practice. Here, we simply consider the rate-independent case where . Taking the above considerations and definitions into account, and using (7), (8), and (12), we arrive at the following coupled system of evolution equations:


where the Cauchy stress tensor is given by


The strong form of the governing equations is complemented by the following boundary conditions


where and are the prescribed boundary displacements and surface tractions, respectively. Additionally, we supplement (13) with initial conditions


for both the displacement field and the damage field. As a final remark in this subsection, we note that the formulation for quasi-static fracture can be retrieved by simply omitting the inertial term in (13).

2.3 Damage initiation and degradation functions

The evolution equation for the damage, (13) allows for the construction of a threshold level of the tensile strain energy before the onset of damage. As a means to trivially satisfy this equation for tensile strain energies below the threshold, we employ a crack driving force given by


where we have introduced as the critical fracture energy per unit volume of material


where is the critical tensile strength and the Young’s modulus.

Alternative motivations for such functions can be found in the work of Miehe et al. (2015a, b), for example, where analogous expressions have been used to generalize the onset of damage to a broader class of thresholds. For the moment, we simply note that this function is used in place of the tensile strain energy in (13), arriving at the following system of coupled equations:


The utility of this crack driving force as a threshold for damage initiation is discussed further in Section 3.1.


[width=0.31]gl1-eps-converted-to.pdf (a){overpic}[width=0.31]gl3-eps-converted-to.pdf (b){overpic}[width=0.31]gl10-eps-converted-to.pdf (c){overpic}[width=0.3]gq3-eps-converted-to.pdf {overpic}[width=0.31]gq3-eps-converted-to.pdf (d){overpic}[width=0.31]gq10-eps-converted-to.pdf (e)

Figure 3: Influence of the constant on the quasi-linear degradation function (top) and the quasi-quadratic degradation function (bottom). For illustration purposes (a); (b,d) and (c,e) are considered.

The formulation is completed by the definition of a stiffness degradation function. In this paper, we are concerned with two particular functions, which hereafter we refer to as the quasi-linear and quasi-quadratic degradation functions.

The quasi-linear and quasi-quadratic degradation functions are both rational functions of the damage field, with the numerator being either a linear function of damage or a quadratic function of damage, respectively. Each of the functions has an associated upper bound on the regularization length scale that can be used in conjunction with it. The quasi-linear degradation function is given by


where is a constant. This function and the associated bounds on and the regularization length are derived in Section 3.2.

The quasi-quadratic degradation function follows the work of Lorentz et al. (2011) and Lorentz (2017), and is given by


where is a shape parameter. The dependence of both degradation functions on the constant and shape parameter is graphically depicted in Figure 3, for illustration purposes. We note that both degradation functions satisfy the condition that


In the present work, we require that


The justification for this choice is provided in Section 3.1. This leaves only the shape parameter to be defined in (21). In the remainder of this paper, unless otherwise noted, we set .

In general, different choices of the shape parameter will give rise to different fracture responses, and influence the load-deflection response of the structure as well as the fracture evolution. A thorough examination of this effect is beyond the scope of this manuscript, and we refer the reader to the discussions in Lorentz (2017) as to how multiple shape parameters can be employed to effectively calibrate the model to individual materials.

The constraints on the regularization length and shape parameter that go along with (20) and (21) ensure an increasing damage band width, and a decreasing stress. The latter essentially stems from the objectivity condition in cohesive zone modeling, as discussed in Bažant (2002). Finally, we note that the upper bound on the regularization length is closely related to the characteristic length of the fracture process zone, which was approximated as in Hillerborg et al. (1976).

3 Analysis & discussion

We now provide some analysis of the model described in Section 2. We begin by examining a simple one-dimensional problem. We then provide the derivation of the quasi-linear degradation function. Finally, some aspects concerning the regularization of crack topologies and the issue of -convergence are discussed.

3.1 Analytical solution of the one-dimensional rate-independent problem

To illustrate the fundamental properties of our formulation, we begin by studying a particular boundary value problem. Consider a one-dimensional bar subjected to a uniaxial tensile load. A symmetrical solution for the localization band is expected in a domain of interest . Without a loss of generality, this point of symmetry is taken at center of the bar, i.e. . We assume a non-negative strain-field and ignore inertial effects. Under these assumptions (19) simplifies to


with , where denotes Young’s modulus.

Provided that the strain and damage are monotonically increasing functions of time at every point in the domain, can be written as


We begin our analysis by considering the situation before the onset of any damage in the bar. In such a configuration, the damage everywhere, and the stress is spatially constant and below the critical stress, i.e. . Under these conditions (24) and (25) simplify to


In order for this equation to be satisfied, the two terms must balance. This is trivial, provided that


We note that, given the choice of in (23), this constraint is satisfied for both the quasi-linear, (20), and quasi-quadratic, (21), degradation functions. By contrast, the use of a simple quadratic degradation function of requires the use of a particular value of for this constraint to be satisfied.

We turn now to considering the onset of damage and the post-critical behavior. If we want damage to begin when , we can effect this by setting


which relates the critical fracture energy per unit volume to the tensile strength . At the onset of damage and beyond, (24) and (25) simplify to the following nonlinear ordinary differential equation:


A solution for the damage field can be found by making use of the fact that the stress across the bar is constant and multiplying (29) with , to obtain


Because of symmetry about , (30) is integrated from to for positive values of , and from to when is negative. Since we have specified the crack to be centered at , we require the damage field to have a maximum at . From this requirement, we obtain


where is given by


By definition, substitution of the homogeneous solution into (31) yields a zero damage gradient. For a fully developed crack, i.e. and , (31) reduces to


and by applying the boundary conditions and , we find


as the solution that satisfies the specified boundary conditions. We note that this regularization is different from what is traditionally used in the phase-field modeling of brittle fracture. The different ultimate damage distributions are compared in Figure 4. We refer the reader to Borden et al. (2012) for an analogous investigation in the context of phase-field models for brittle fracture.


[width=.95]nonlocal-eps-converted-to.pdf 111(a)(b)(c)

Figure 4: Sharp and diffuse crack topologies for a crack at . (a) Sharp crack. (b) Smeared crack modeled with a ‘standard’ phase-field approach. (c) Smeared crack modeled with the proposed formulation.

Given a closed-form analytical solution for the ultimate damage profile in a one-dimensional setting, we calculate the corresponding dissipated energy in the system


For a fully developed crack we can deduce from (30) that this differential equation reduces to:


with the ultimate damage distribution. This reduces (35) to


which effectively demonstrates that (8) indeed constitutes a crack surface density functional, in the phase-field for fracture sense. The same line of reasoning was presented in the gradient damage framework of Lorentz and Godard (2011).

3.2 Derivation of quasi-linear degradation function

Up to this point in this Section, the only conditions that have been placed on the degradation function concern its slope at the onset of damage. In this subsection, we demonstrate how a desired stress-strain behavior in the damage zone can be effected through a particular choice of the degradation function. Attention is focused on obtaining a relatively simple, linearly decaying traction-separation law in the damaged region. Such a response is common in standard cohesive methods for fracture, see e.g. Ortiz and Pandolfi (1999). Here, we show how such a response gives rise to the quasi-linear degradation function (20).

A general discussion on the requirements of a degradation function can be found in the work of Pham and Marigo (2013). Generally speaking, a valid degradation function is assumed to satisfy the following set of conditions:


We consider the local response (i.e. within the localization band) in a one-dimensional system, and assume that all fields are spatially uniform. In this setting, the stress-strain relationship simplifies to


The critical stress at which damage initiates can be associated with a critical strain through . Assume that the desired stress-strain relationship in the post-critical regime is given by the following linear decay


where denotes the final strain at which the stress vanishes. This contrasts to most other phase-field approaches in which the stress vanishes only in the limit as .

Combining (39) and (40), we obtain


As the elastic regime is overstepped in this simplified one-dimensional setting, (24) reduces to


Inserting the strain (41) into this local evolution equation, we obtain a nonlinear ordinary differential equation for the degradation function of the form


where are constants which depend on and . It is easy to show that the quasi-linear degradation function (20) is a generic solution to this equation, satisfying and .

As the critical strain is naturally accounted for by fixing the critical strength, we are left to determine the final strain , which is reached when . In that case, (42) becomes


which by means of inserting (20), reduces to


For this model to have a clear physical interpretation, it is obviously required that 1. This constraint can be rewritten in terms of an upper bound on the regularization length scale , given by


In case the damage has reached unity and the stress has dropped to zero, the “opening displacement” within the damage band can be approximated as


This quantity which has dimensions of displacement does not depend on the regularization length scale. This is consistent with the desire to asymptotically approach a cohesive model with a known maximum opening displacement. As the damage localization must occur over a smaller length as , the average strain over this localization zone must increase to compensate. To confirm that the post-critical stress-strain behavior matches (40), a series of numerical experiments are conducted in Section 5.

While the above derivation was employed to extract a linear decay of the stress-strain behavior in the localization band, more complex constitutive relationships can obviously be obtained by following a similar procedure. Finally, we note that the form of the resulting quasi-linear degradation function (20) differs markedly from the quasi-quadratic function (21) as . In particular, the slope does not vanish as complete damage is approached. An unfortunate consequence of this feature is that there is no means to constrain the damage from exceeding unity, and we will find the need to prevent the overshoot in our discrete formulation.

3.3 On the regularization of crack topologies

We now investigate the approximations of the fracture energy that are employed by the proposed cohesive model. Our investigation follows the work described in Miehe et al. (2010b) for the numerical investigation of error in regularized crack surface representations. In particular, we examine a dimensionless version of the evolution equation (19), given by


where is a prescribed form of the driving force.

We consider two separate boundary-value problems. In the first, the prescribed driving force is set to , and the Dirichlet constraint


is enforced on the damage field along the sharp crack surface. A standard approach is to generate a finite-element approximation to the solution of (48)-(49), and then examine the ability of the corresponding (dimensionless) fracture energy to capture the exact energy, as a function of regularization length and mesh spacing.

In the second boundary value problem, we replace the Dirichlet boundary condition (49) with a prescribed driving force that is sufficiently large in the vicinity of the sharp crack. In effect, the magnitude of the driving force is increased significantly (approximating a Dirac distribution) at quadrature points along the sharp crack surface to drive the regularization. By adopting such an approach, the regularized topologies are expected to be more representative of those obtained in fully coupled displacement-damage calculations.


[width=0.85]model_problem.pdf (a)(b)(c)(d)

Figure 5: (a) The model problem from the work of Miehe et al. (2010b). A damage field is approximated using a uniform 401x401 finite element mesh and a large, prescribed driving force along the sharp crack front. Results are shown for the regularized crack topologies and regularization lengths of (a) , (b) and (c) .

[width=0.31]g3-eps-converted-to.pdf (a){overpic}[width=0.31]g2-eps-converted-to.pdf (b)

Figure 6: Regularization results using (a) Dirichlet boundary conditions along the crack surface and (b) prescribing a large driving force along the crack surface.

In the following numerical example, the model problem from Miehe et al. (2010b) is revisited with the intent of providing some insight into the -convergence properties of the proposed method. Consider a two-dimensional continuum with an embedded crack surface from the left side to the center of the domain as depicted in Figure 5a. The regularization is driven by prescribing the crack driving forces to a value of at the quadrature points along . A similar approach was adopted by Borden et al. (2012) for modeling preexisting cracks in a continuous body. The computations are performed on uniform finite element meshes consisting of bilinear quadrilateral elements. The regularized crack functionals


are then calculated as a post-processing step of the boundary value problem (48)-(49).

In Figure 5b-d the computed regularized cracks are shown for different values of the regularization length. The results are in excellent agreement with the ones reported in Miehe et al. (2010b). The results for the numerical -convergence investigation are summarized in Figure 6. In the Dirichlet case, the calculated fracture energy approaches the theoretical estimate as the mesh is sufficiently refined. The results suggest that better results are obtained for smaller ratios of the regularization length to the mesh size. This is counter-intuitive, and the results when the problem is driven by a prescribed driving force, Figure 6b, suggest that they are simply an artifact of the Dirichlet boundary condition. The driving force results indicate that while the error between the sharp and regularized crack surfaces may eventually asymptote to zero, much more refined meshes will be required.

4 Finite element implementation

We now discuss the numerical discretization of the evolution equations (19). We begin by discussing existing strategies for enforcing the irreversibility constraint, (19), and then present the finite-element formulation. Finally, details of the staggered solution algorithm are provided.

4.1 Augmented Lagrangian implementation

A now standard approach to enforcing the irreversibility constraint is to replace the crack driving function in (19) with a monotonic version :


over the full temporal history . Miehe et al. (2010a) demonstrated that such a strategy leads to a thermodynamically consistent phase-field model for fracture. However, while such an approximation effectively enforces irreversibility, it is variationally inconsistent. An alternative approach using an augmented Lagrangian method was proposed in Wheeler et al. (2014), in which the solution to the evolution equations is obtained from a constrained minimization problem.

In the following, we briefly review the augmented Lagrangian method described in Wheeler et al. (2014). The formulation starts with the energy function from (6) which is minimized with respect to the unknown solution variables, which are the displacements and continuous damage variable . Differentiation with respect to and leads to the Euler-Lagrange equations. Given a previous solution we seek a solution


which is approximated as


where are the Lagrange multipliers and is a penalty kernel. An iterative scheme for solving this constrained optimization problem is described in Section 4.3.

We note that the first penalty term in (53) enforces the monotonicity constraint (19) on the damage field, while the second prevents the damage field from exceeding unity. While we have not found the second condition to be necessary with the quasi-quadratic degradation function, we have found it to be essential with the quasi-linear degradation function. Additional comments to this effect are provided in Section 5.

4.2 Galerkin finite element discretization

To approximate the solution to the energy functional from (53) using a finite element method, we employ a finite element discretization as follows. The admissible function spaces for the displacements, for the damage field, and for the Lagrange multiplier field are defined as


respectively, while the associated weighting spaces are given by


Next the differentiation of (53) is performed with respect to and . After multiplication by the appropriate weighting functions and integrating by parts, we obtain


as well as


Where the bilinear forms and were introduced for the displacement field and damage field problems, and denotes the inner product on .

Following a standard Galerkin method, the problem is recast in finite-dimensional subspaces and . Given any we find satisfying


and similarly satisfying


In this work, standard finite elements are used along with full integration. In particular, we employ bilinear quadrilateral elements to approximate the displacement and damage fields. The Lagrange multiplier is also approximated with bilinear quadrilateral elements.

We note that in the phase-field literature, it is common to introduce a small regularization parameter which provides a lower bound on the tensile stresses as the damage reaches unity, as in Miehe et al. (2010a, b). This is a form of numerical regularization that provides completely damaged elements with a small degree of stiffness as a means to circumvent a loss of ellipticity in the discrete equations. In this work, however, we have not found such a parameter to be necessary. We note that this is not necessarily the result of our cohesive formulation. Indeed, recent phase-field implementations based on Griffith models have reported similar findings, as in Borden et al. (2012). As a result, none of the calculations presented in Section 5 rely on any form of regularization beyond the aforementioned penalty kernel.

4.3 Staggered solution scheme

We now describe an efficient solution procedure based on a staggered solution scheme for the successive update of the damage and displacement fields. In the first part of the algorithm, an outer loop is considered on the update for the damage evolution equation and the explicitly calculated Lagrange multipliers, while the displacement field is held fixed. A new damage field is obtained once the algorithm converges within a predefined tolerance . The second part of the algorithm then holds this damage field fixed while considering the elastodynamics problem, which is advanced using an explicit algorithm based on the central difference time integration scheme. As is standard with explicit algorithms, we employ a lumped mass matrix.

We note that this approach is distinct from that of Wheeler et al. (2014), in which staggered iterations were used in an attempt to more closely capture the true equilibrium solution through an iterative procedure. In general, such a strategy is significantly more expensive per time or load step, and a detailed investigation into the optimal trade-off between accuracy and computational time has yet to be performed. The approach presented here, summarized in Algorithm 1, relies on using sufficiently small time steps to ensure accuracy. In general, we have not found the need to use time steps that are significantly below the stability limit as established by the usual CFL condition in explicit dynamics.

Choose = , , and let
     Solve the nonlinear damage equation using a Newton-type scheme and initial guess
     Update the Lagrange multipliers
Solve the elastodynamics problem
Algorithm 1 Staggered solution scheme in the time interval .

5 Numerical results and validation

We now present the results from a series of representative numerical experiments designed to demonstrate the ability of the proposed approach to capture pertinent aspects of cohesive fracture processes. The goal of these examples is threefold: (i) to numerically verify the convergence properties of the proposed formulation with respect to the regularization length as well as to provide insight into its interpretation in a multidimensional setting; (ii) to examine various options at imposing irreversibility, as discussed in Section 4.1; and (iii) to demonstrate the performance of the regularized cohesive model for benchmark problems in quasi-static and dynamic fracture.

In our calculations, we rely on meshes that are sufficiently refined to capture the variation in the damage field around crack surfaces. Unless specified otherwise, we use element sizes of as a rule of thumb. We also rely on the quasi-quadratic degradation function (21) with , except where indicated. For all two-dimensional numerical experiments that are considered in this Section, plane strain conditions are assumed to hold. Finally, all calculations employing the augmented Lagrangian method use a tolerance of in the outer loop. In our experience, some problem-dependent fine-tuning of the penalty kernel is required to obtain good convergence behavior in the augmented Lagrangian approach. In practice, we have employed penalty kernels in the range of , which is consistent with those reported in Wheeler et al. (2014). In all cases, the monotonicity penalization is employed, whereas the threshold penalization is only employed in conjunction with the quasi-linear degradation function.

5.1 A one-dimensional bar under tension

Consider a one-dimensional bar with a reduced cross sectional area in the middle, subjected to a tensile load applied to the right end, as shown in Figure 7. The material parameters are  MPa,  N/mm and  MPa. The bar has a length of mm and a nominal cross-sectional area of mm.

Numerically demonstrating convergence towards a cohesive zone formulation for a vanishing regularization length is a non-trivial task. To ensure the equilibrium solution is captured accurately at every load step, we depart slightly from the algorithm described in Section 4.3 and implement a monolithic solution scheme for this problem. A path-following constraint was introduced to control the loading process, as in Verhoosel et al. (2009). In a one-dimensional setting, a constant mesh size of was employed, providing sufficient resolution to resolve the damage profile for all of the regularization lengths under consideration.

The influence of the regularization length on the force-displacement curve and damage profile is shown in Figure 8. Here the reaction force , normalized by a peak force , is plotted as a function of the displacement at the right end of the bar, normalized by its value at the onset of damage. Prior to softening, the force-displacement curves superimpose, as the constitutive response is purely elastic. Once the critical stress in the bar is exceeded, damage is initiated, and strain softening begins. Despite significant differences in the corresponding damage distributions, the regularization length clearly has a negligible effect on the structural response. These results confirm similar ones made in Lorentz et al. (2011), in which it was demonstrated that the one-dimensional response of the proposed phase-field/gradient damage model converges toward that of a cohesive zone model in the limit of vanishing regularization length.


[width=.44]1d.pdf ()

Figure 7: Bar with reduced cross sectional area in the middle, subjected to a tensile load .

[width=.44]cpf_fd-eps-converted-to.pdf (a){overpic}[width=.44]cpf_d-eps-converted-to.pdf (b)

Figure 8: Effect of the regularization length on (a) the normalized force-displacement curve, and (b) the damage field using a constant shape parameter .

[width=.44]p_fd-eps-converted-to.pdf (a){overpic}[width=.44]p_d-eps-converted-to.pdf (b)

Figure 9: Effect of the cohesive shape parameter on (a) the force-displacement curve, and (b) the damage field when the regularization length is held fixed at .

[width=.44]gl_fd-eps-converted-to.pdf (a){overpic}[width=.44]gl_d-eps-converted-to.pdf (b)

Figure 10: Effect of the regularization length on (a) the normalized force-displacement curve, and (b) the damage field using quasi-linear degradation function (20).

Next we study the impact of the softening shape parameter on the constitutive response. For a fixed regularization length, an increasing value for this parameter leads to a more rapid decay of the stress in the early stages of the strain-softening process, as demonstrated in Figure 9. This behavior is also implied by the form of the quasi-quadratic degradation function (21). In addition, we note that the shape parameter has minimal influence on the ultimate damage distribution.

We now study the same one-dimensional problem, but using the quasi-linear degradation function (20). In this case, the upper bound on the regularization length is computed to be mm. We therefore consider regularization length scales of mm, mm and  mm. We recall that (20) does not possess any shape parameters as it is solely designed to approximate a bilinear stress-stress response upon strain-softening.



Figure 11: A quasi-linear stress-strain response at the center of the bar for the cohesive model when employing the quasi-linear degradation function (20) and .

The non-dimensional force-displacement curves for the quasi-linear degradation function are given in Figure 10. Similar to the results obtained with the quasi-quadratic degradation function, we see virtually no sensitivity to the regularization length scale. However, the constitutive response markedly differs from the one obtained through the quasi-quadratic degradation function. The force-displacement response is a lot softer, while the stresses rapidly drop to zero as . The findings are consistent with the arguments outlined in Section 3.2. In addition, Figure 11 demonstrates that the stress-strain response at the center of the bar effectively approximates a linear decay in the post-peak regime.

5.2 Single edge notched tests


[width=0.7]testing.pdf (a)(b)0.5 mm0.5 mm0.5 mm0.5 mm0.5 mm0.5 mm0.5 mm0.5 mmmmmm

Figure 12: Geometry and boundary conditions for the single-edge notched tests (a) in tension, and (b) in shear.

We now investigate two benchmark problems, which have become canonical in the phase-field for fracture literature. Consider a square plate with an initially horizontal edge crack extending to the middle of the specimen, as shown in Figure 12a. Following Miehe et al. (2010b), the material parameters are chosen to be GPa, , and N/mm. In the following, the cohesive strength is set to GPa. The corresponding estimate for the length of the process zone is mm. The computations are performed in a displacement driven context where the displacement increment is adjusted upon approaching the peak load to ensure accuracy. To accurately capture the evolution of the damage field, the mesh is refined along the anticipated crack path. In addition, the initial notch was modeled with a small radius mm to obtain a more realistic representation of the stresses in the vicinity of the initial crack tip.

First we apply a vertical displacement to the complete top edge, as demonstrated in Figure 12a. The fracture patterns for different values of the regularization length are shown in Figure 13. The blue and red contour levels indicate an intact and fully damaged material state, respectively. In contrast to results obtained for this problem using phase-field approximations of a Griffith model of fracture, we do not observe any evolution of the damage field away from the regularized fracture surface.

The same specimen is now subjected to a shear load, as depicted in Figure 12b. We note that the crack trajectories presented here are in excellent agreement with the ones reported in Miehe et al. (2010a) and display little sensitivity to the regularization length scale. For a more comprehensive investigation into the effect of a vanishing regularization length, we refer to the trapezoid problem investigated in the next subsection.


[width=0.75]mode1_results.pdf (a)(b)(c){overpic}[width=0.015]legend.png

Figure 13: Single edge notched test under tension for (a)  mm, (b)  mm and (c)  mm.

[width=0.75]mode2_results.pdf (a)(b)(c){overpic}[width=0.015]legend.png

Figure 14: Single edge notched test under shear for (a)  mm, (b)  mm and (c)  mm using shape parameter .

5.3 Trapezoid problem

We now examine a problem first proposed by Lorentz and Godard (2011) in which a particular geometric setup is designed to facilitate stable crack propagation without exhibiting any snap-back behavior. In addition to this, the dimensions of the virtual specimen, several meters in size, are significantly larger then the characteristic size of the fracture process zone in plain concrete. This separation of scales is less distinct in the single-edge notched specimens from Section 5.2. For these reasons, the trapezoid problem is viewed as being particularly well-suited for an investigation into the effects of a vanishing regularization length.


[width=0.28]trapezoid.pdf 4000 mm3000 mm9400 mm{overpic}[width=0.45]trap_results.pdf {overpic}[width=0.015]legend.png

Figure 15: (a) Geometry and boundary conditions for a trapezoidal-shaped fracture specimen, taken from Lorentz and Godard (2011). (b) Fracture pattern and the impact of the shape parameter on the effective process zone.

The material parameters are chosen to be representative of concrete, namely GPa, , N/mm and MPa. The corresponding size of the fracture process zone is mm. The analysis was performed for two different values of the cohesive shape parameter, namely and . The finite element mesh is refined along the expected crack path.


[width=0.44]trap1-eps-converted-to.pdf (a){overpic}[width=0.44]trap2-eps-converted-to.pdf (b)

Figure 16: Load-displacement curves for different values of the regularization length of the trapezoid problem. The shape parameter under consideration equals (a) , and (b) .

The force-displacement curves for this problem are shown in Figure 16, and are in excellent agreement with those reported in Lorentz and Godard (2011) and Lorentz et al. (2012). The results illustrate the significance of the shape parameter and its impact on the critical load the structure is able to sustain. An increasing value for can also considerably elongate the length of the fracture process zone, as indicated in Figure 15b. In this case, the stresses in the diffusive process zone are effectively smeared over a larger distance, which actually facilitates a convergence study with respect to the regularization length. For a value of , the differences in the force-displacement curves for the various regularization lengths are negligible. As the fracture process zone decreases in size, strain components in the plane parallel to the crack become more significant. As a result, there is greater separation in the force displacement curves for , as indicated by 16. Given these conditions, we can effectively demonstrate the insensitively of the regularization length to the constitutive response. We refer the reader to Lorentz et al. (2012) for a more comprehensive analysis of this problem.

5.4 Three-point bending test


[width=0.55]3pb.pdf 600 mmCMOD700 mm150 mm50 mm()(a){overpic}[width=0.51]3pb_result.png (b)

Figure 17: (a) Dimensions and boundary conditions for the three-point bending test, from Roesler et al. (2007). (b) Associated damage pattern for the three-point bending test in the deformed configuration with degradation function (21) and . The displacements have been scaled by a factor of 50 and areas of the model where have been removed in order to show a representation of the fractured geometry.

We now consider a classical benchmark problem in cohesive fracture mechanics, the three-point bending (TPB) test. Consider the problem indicated in Figure 17a. The single notched concrete specimen is pushed downwards with imposed displacements . Experimental results of interest for this problem have been reported in Roesler et al. (2007). Previous numerical results based on cohesive zone approaches have been reported for this problem by Roesler et al. (2007) and Kim and Duarte (2014), among others.

The material parameters are chosen as GPa, , MPa, and N/mm. The regularization length is set to mm. The length of the fracture process zone is mm, which is fairly large relative to the overall dimensions of the beam. The mesh is refined in the center of the beam where the crack is expected to propagate.


[width=0.44]3pb2-eps-converted-to.pdf (a){overpic}[width=0.44]3pb1-eps-converted-to.pdf (b)

Figure 18: Three-point bending test: reaction force versus crack mouth opening displacement (CMOD) using (a) the quasi-linear degradation function (b) the quasi-quadratic degradation function.

We examine the influence of the quasi-linear and quasi-quadratic degradation functions on the constitutive response. It is emphasized that the former function requires additional penalization to prevent the approximation to the damage field from exceeding unity, as discussed in Section 4.1.

Figure 18 shows the computed force-displacement curves for both degradation functions under consideration. We emphasize that all calculations employ the same values for the critical energy and cohesive strength. Nevertheless, the results shown in Figure 18 indicate that the peak load sustained by the structure is sensitive to the particular degradation function. The quasi-linear degradation function gives rise to the largest peak load. As the shape parameter in the quasi-quadratic function is increased, the peak load decreases. We note that among the results shown here, those obtained using the quasi-quadratic degradation function and compare most favorably to the experimental load-displacement results reported by Roesler et al. (2007).

5.5 Dynamic crack branching

We now examine a standard benchmark problem in dynamic fracture that has been extensively investigated in previous works, including Belytschko et al. (2003) and Borden et al. (2012), to name a few. Consider a notched rectangular plate of dimensions  mm dynamically loaded in tension, as shown in Figure 19. The domain contains an initial crack that is horizontal and spans from the left edge of the plate to its center. A tensile load of  MPa is applied to the top and bottom surfaces and held constant throughout the course of the simulation.


[width=0.52]branching.pdf 40 mm20 mm100 mm50 mmMPa

Figure 19: Geometry and boundary conditions for a dynamic crack branching problem.

The material parameters are set to  GPa, ,  N/mm, MPa, and  kg/m. The corresponding dilatational, shear and Rayleigh wave speeds are  m/s,  m/s and m/s, respectively. The characteristic length of the fracture process zone is approximately mm. The regularization length is set to mm. The computation is performed on a uniform mesh with spacing .


[width=0.95]branching_results.pdf (a)(c)(b)(d)

Figure 20: Evolution of the damage field with time for the crack branching problem: (a) s; (b) s; (c) s and (d) s. The displacements have been scaled by a factor of 50. In addition, areas of the model where have been removed in order to show a representation of the fractured geometry.

[width=0.44]branching_diss-eps-converted-to.pdf (a){overpic}[width=0.44]branching_elastic-eps-converted-to.pdf (b)

Figure 21: Evolution of (a) the dissipated energy and (b) the elastic strain energy through time for the Lagrange multiplier (LM) and history-based formulations for the crack branching problem.


Figure 22: Comparison of the damage distributions across line segment for the two different strategies for enforcing the irreversibility constraint.

The damage field at selected time steps is shown in Figure 20. The proposed model for cohesive fracture is clearly capable of capturing the bifurcation of a rapidly evolving cohesive crack. Although cracks are not tracked algorithmically in this approach, crack-tip velocities can be estimated by post-processing a fixed damage iso-contour. Our results indicate that the crack tip velocity remains below 60% of the Rayleigh wave speed at all times, agreeing with previous numerical investigations of this problem as well as the experimental observations of Ravi-Chandar and Knauss (1984).

This particular problem also lends itself to an investigation of the impact of the various strategies for enforcing the irreversibility constraint. In Figure 21, the evolution of the elastic strain and dissipated energies is displayed for both the history and Lagrange multiplier methods of imposing a monotonically increasing damage field. While the history formulation is computationally more efficient (requiring a single solve of the nonlinear damage equation per time step), it results in a greater degree of energy dissipation. In the calculations, this manifests itself mainly through a widening of the regularized crack surface as compared to the results obtained with the Lagrange multiplier implementation. This is clearly demonstrated by looking at a cross-section A-A across the regularized fracture surface, shown in Figure 22.

5.6 The Kalthoff-Winkler experiments

We now focus on a second benchmark problem in dynamic fracture, namely models of the experiments by Kalthoff and Winkler (1987). Consider a plate with two edge notches that is impacted by a projectile, as shown in Figure 23. The notches extend halfway through the plate width. Due to symmetry, only the upper half of the specimen is explicitly modeled in the simulations that follow.


[width=0.35]kalthoff.pdf mmmmmmmmmmmm

Figure 23: The Kalthoff problem: geometry and boundary conditions.

[width=0.9]kalthoff_results.pdf (a)(c)(b)(d)(e){overpic}[width=0.015]legend.png +

Figure 24: Left: damage fields for the Kalthoff problem at various times: (a) s; (b) s; (c) s and (d) s. The displacements are magnified by a factor of 3. In addition, areas of the model where have been removed in order to show a representation of the fractured geometry. (e) The resulting crack propagation occurs at an angle of approximately 63.

The experiments of Kalthoff and Winkler (1987) exhibited different fracture/damage behaviors of a maraging steel material as a function of the impact velocity. For relatively low impact velocities, brittle fracture was observed with propagation angles of approximately 70 from the original crack plane. For higher impact velocities, failure was observed to occur due to shear localization originating from shear band formation ahead of the notch. In this work, we study the brittle failure mode. In particular, we are concerned with the load case where the impact velocity of the projectile is m/s. Assuming that the projectile has the same impedence as the specimen, we apply only half of its velocity. Following Belytschko et al. (2003), we therefore consider a velocity boundary condition applied to the notched area with magnitude m/s, linearly ramped up from zero over a time interval of s.


[width=0.45]kalthoff_diss-eps-converted-to.pdf (a){overpic}[width=0.45]kalthoff_elastic-eps-converted-to.pdf (b)

Figure 25: Evolution of (a) the dissipated energy and (b) the elastic strain energy with time for the Lagrange multiplier and history-based formulations, for the Kalthoff problem.

The maraging steel used in the original experiments is modeled using the following material properties:  GPa, , N/mm and kg/m. This leads to dilatational, shear and Rayleigh wave speeds of m/s, m/s and m/s, respectively. The cohesive strength has a value of  GPa. To avoid the use of a fairly small regularization length, a frequently adopted strategy in the phase-field community is to lower the effective strength of the material to facilitate a computational analysis, as in Borden et al. (2012). In that work, the maximum uniaxial tensile strength was set to 1.07 GPa. We found that the use of an artificially small cohesive strength gives rise to secondary cracking at the bottom of the specimen, very similar to what was reported in Borden (2012) and Moreau et al. (2015). Here, we report simulation results using the “true” cohesive strength of  GPa, which gives rise to a fracture process zone with characteristic size mm. We conduct simulations using a regularization length of  mm, using a uniform mesh with mesh spacing .

The damage field at selected times during the simulation is shown in Figure 24a-d. Our results indicate a crack propagation with an average orientation of approximately , as shown in Figure 24e. The result compares favorably with experimental observations of an angle of roughly , as well as the result from the phase-field simulation of Borden et al. (2012). Figure 25 shows the evolution of the dissipation and elastic energy functionals through time for the various irreversibility implementations. The history implementation once again is seen to dissipate more energy than the Lagrange multiplier approach advocated in this work.

5.7 Fragmentation of a thick cylinder


[width=0.75]song.pdf (a)(b)160 mm300 mm

Figure 26: (a) Geometry and boundary conditions for the fragmentation of a thick cylinder. (b) one sample of a spatially random Young’s modulus field.

In order to demonstrate the capabilities of the proposed method for simulating fragmentation problems, we consider a thick cylinder subjected to an impulsive internal pressure, as shown in Figure 26a. The cylinder has inner and outer radii of 80 and 150 mm, respectively. The inner surface is subjected to a spatially constant pressure that evolves in time as , after an initial linear ramp to MPa over a rise time of s. This problem has been studied previously by Song and Belytschko (2009), using a cracking-node algorithm, and by Hirmand and Papoulia (2018) using a particular cohesive zone model, to mention a few.

The material parameters are GPa, kg/m, and . This leads to dilatational, shear and Rayleigh wave speeds of m/s, m/s and m/s, respectively. Given the critical strain of provided by Song and Belytschko (2009), the cohesive strength was set to GPa. The fracture energy was set to N/mm, while the regularization length scale was set to mm. The size of the process zone associated with this set of material parameters is mm.


[width=0.95]song_results.pdf (a)(b)(c)

Figure 27: Damage field for the thick-walled cylinder fragmentation problem on mesh 3 at times: (a) s; (b) s; and (c) s.

[width=0.95]song_results2.pdf (a)(b)(c)

Figure 28: Final crack patterns in the reference configuration for (a) mesh 1 (b) mesh 2 and (c) mesh 3.

[width=0.45]tc_diss-eps-converted-to.pdf (a){overpic}[width=0.45]tc_elas-eps-converted-to.pdf (b)

Figure 29: Evolution of the (a) dissipated energy and (b) elastic strain energy with time for the thick cylinder fragmentation problem for a sequence of meshes of increasing refinement.

To introduce a stochastic element to the simulations, we construct a spatially-varying Young’s modulus field according to


where is the mean value of the random field; the uncorrelated Gaussian random variable and the set of eigenvalues and basis functions associated with the Fredholm equation of the second kind, respectively; and is the number of basis functions. We construct the spatially varying Young’s modulus field shown in Figure 26b using basis functions. The resulting field possesses a mean of GPa, with a correlation length of mm, and a standard deviation of 5%. Additional details concerning the construction of such fields are provided in Shang and Yun (2013).

We now hold the regularization length and spatially varying Young’s modulus field fixed, and examine the fracture patterns predicted by our simulations over a sequence of increasingly refined meshes. Table 1 provides the number of elements for each mesh and the corresponding ratio of regularization length to element size.

Mesh 1 Mesh 2 Mesh 3
Number of elements 218400 400000 728000
ratio 4 6 8
Table 1: The sequence of meshes considered in the convergence study of the thick cylinder fragmentation problem.

Figure 27 shows the fragmentation process at different times for the most refined mesh, mesh 3. The fragmentation process begins around s with a large number of small cracks along the inner surface. However, with increasing time, some of the initial cracks arrest, and the remainder continue to propagate to the outer surface of the cylinder. The fragmentation process completes around s, when the multiple fragments move outward from the center of the domain with no further crack initiation and virtually no growth of the arrested cracks.

The final damage fields for each mesh are shown in Figure 28. As can be seen, the fragmentation patterns are very similar and, in fact, the number of large fragments does not appear to change upon mesh refinement. This in contrast to the results presented in Song and Belytschko (2009).

Finally, Figures 29a-b compares the dissipated and elastic strain energies for the uniformly refined sequence of meshes from Table 1. The dissipated energy shown in Figure 29a is over-predicted on the coarsest mesh 1, and appears to converge quickly as the mesh is refined.

6 Conclusions

In this work, we have extended a phase-field/gradient damage model, based on the work of Lorentz and Godard (2011), to the dynamic case. Such a formulation is particularly well-suited to operate within the context of cohesive fracture, i.e. when a Griffith description of fracture is insufficient. We demonstrated that the combination of an alternative phase-field approximation, supported by a suitable stiffness degradation function, effectively leads to a regularized description of cohesive fracture. Such a formulation is distinctly characterized by a linear elastic regime prior to the onset of damage and controlled strain-softening thereafter. The governing equations are derived according to macro- and microforce balance theories, naturally addressing the irreversibility of crack growth in time by introducing suitable constraints for the kinetics of the underlying microstructural changes.

Our work departs from that of Lorentz and Godard (2011) in a number of key respects. In particular, we have developed a new degradation function that recovers a linear stress-strain response in the localization zone. We have also employed a tension-compression decomposition of the strain, similar to what is employed in phase-field models of fracture. And we have explored options for enforcing the irreversibility of the damage field. Finally, we have demonstrated the ability of the proposed method to reproduce the results of both idealized benchmark and experimental problems in quasi-static and dynamic fracture.

The current manuscript gives rise to several areas for future investigation. First and foremost, we highlight that no calibration of the degradation and/or corresponding strain-softening parameters was conducted in this work. It is anticipated that such an approach, possibly supported by data-driven optimization methods, could further reduce any discrepancies between simulations and experimental observations. Secondly, we mention the extension of some of the basic elements of a cohesive fracture based phase-field description to an anisotropic setting. Such an extension is expected to be non-trivial, with particular regard to the calculation of the effective properties along given material directions.


This work was performed under a research grant from Sandia National Laboratories, to Duke University. That support is gratefully acknowledged. The work was also partially supported by the Laboratory Directed Research and Development program at Sandia National Laboratories.


  1. journal: Computer Methods in Applied Mechanics and Engineering


  1. Ambrosio, L., Tortorelli, V. M., 1990. Approximation of functional depending on jumps by elliptic functional via -convergence. Communications on Pure and Applied Mathematics 43 (8), 999–1036.
  2. Barenblatt, G. I., 1962. The mathematical theory of equilibrium cracks in brittle fracture. Advances in applied mechanics 7, 55–129.
  3. Bažant, Z. P., 2002. Concrete fracture models: testing and practice. Engineering Fracture Mechanics 69 (2), 165 – 205.
  4. Belytschko, T., Chen, H., Xu, J., Zi, G., 2003. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. International Journal for Numerical Methods in Engineering 58 (12), 1873–1905.
  5. Borden, M. J., 2012. Isogeometric analysis of phase-field models for dynamic brittle and ductile fracture. Ph.D. thesis.
  6. Borden, M. J., Hughes, T. J., Landis, C. M., Anvari, A., Lee, I. J., 2016. A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects. Computer Methods in Applied Mechanics and Engineering 312, 130 – 166, phase Field Approaches to Fracture.
  7. Borden, M. J., Verhoosel, C. V., Scott, M. A., Hughes, T. J., Landis, C. M., 2012. A phase-field description of dynamic brittle fracture. Computer Methods in Applied Mechanics and Engineering 217-220, 77 – 95.
  8. Bourdin, B., Francfort, G., Marigo, J.-J., 2000. Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids 48 (4), 797 – 826.
  9. Bourdin, B., Francfort, G. A., Marigo, J.-J., 2008. The variational approach to fracture. Journal of elasticity 91 (1-3), 5–148.
  10. Bourdin, B., Larsen, C. J., Richardson, C. L., Apr 2011. A time-discrete model for dynamic fracture based on crack regularization. International Journal of Fracture 168 (2), 133–143.
  11. da Silva, M. N., Duda, F. P., Fried, E., 2013. Sharp-crack limit of a phase-field model for brittle fracture. Journal of the Mechanics and Physics of Solids 61 (11), 2178 – 2195.
  12. Dugdale, D. S., 1960. Yielding of steel sheets containing slits. Journal of the Mechanics and Physics of Solids 8 (2), 100–104.
  13. Francfort, G., Marigo, J.-J., 1998. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids 46 (8), 1319 – 1342.
  14. Frémond, M., Nedjar, B., 1996. Damage, gradient of damage and principle of virtual power. International Journal of Solids and Structures 33 (8), 1083 – 1103.
  15. Hillerborg, A., Modéer, M., Petersson, P.-E., 1976. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research 6 (6), 773 – 781.
  16. Hirmand, M. R., Papoulia, K. D., 2018. A continuation method for rigid-cohesive fracture in a discontinuous galerkin finite element setting. International Journal for Numerical Methods in Engineering 115 (5), 627–650.
  17. Hofacker, M., Miehe, C., 2013. A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns. International Journal for Numerical Methods in Engineering 93 (3), 276–301.
  18. Huespe, A., Needleman, A., Oliver, J., Sánchez, P., 2009. A finite thickness band method for ductile fracture analysis. International Journal of Plasticity 25 (12), 2349 – 2365.
  19. Kalthoff, J., Winkler, S., 1987. Failure mode transition at high rates of shear loading. International Conference on Impact Loading and Dynamic Behavior of Materials 1, 185–195.
  20. Kim, J., Duarte, C. A., 2014. A new generalized finite element method for two-scale simulations of propagating cohesive fractures in 3-d. International Journal for Numerical Methods in Engineering 104 (13), 1139–1172.
  21. Li, T., Marigo, J.-J., Guilbaud, D., Potapov, S., 2016. Gradient damage modeling of brittle fracture in an explicit dynamics context. International Journal for Numerical Methods in Engineering 108 (11), 1381–1405, nme.5262.
  22. Lorentz, E., Jun 2017. A nonlocal damage model for plain concrete consistent with cohesive fracture. International Journal of Fracture.
  23. Lorentz, E., Cuvilliez, S., Kazymyrenko, K., 2011. Convergence of a gradient damage model toward a cohesive zone model. Comptes Rendus Mécanique 339 (1), 20 – 26.
  24. Lorentz, E., Cuvilliez, S., Kazymyrenko, K., 2012. Modelling large crack propagation: from gradient damage to cohesive zone models. International Journal of Fracture 178 (1), 85–95.
  25. Lorentz, E., Godard, V., 2011. Gradient damage models: Toward full-scale computations. Computer Methods in Applied Mechanics and Engineering 200 (21–22), 1927 – 1944.
  26. Miehe, C., 2011. A multi-field incremental variational framework for gradient-extended standard dissipative solids. Journal of the Mechanics and Physics of Solids 59 (4), 898 – 923.
  27. Miehe, C., Hofacker, M., Schänzel, L.-M., Aldakheel, F., 2015a. Phase field modeling of fracture in multi-physics problems. part ii. coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic–plastic solids. Computer Methods in Applied Mechanics and Engineering 294, 486 – 522.
  28. Miehe, C., Hofacker, M., Welschinger, F., 2010a. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering 199 (45–48), 2765 – 2778.
  29. Miehe, C., Schänzel, L.-M., Ulmer, H., 2015b. Phase field modeling of fracture in multi-physics problems. part i. balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids. Computer Methods in Applied Mechanics and Engineering 294, 449 – 485.
  30. Miehe, C., Welschinger, F., Hofacker, M., 2010b. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations. International Journal for Numerical Methods in Engineering 83 (10), 1273–1311.
  31. Moreau, K., Moës, N., Picart, D., Stainier, L., 2015. Explicit dynamics with a non-local damage model using the thick level set approach. International Journal for Numerical Methods in Engineering 102 (3-4), 808–838.
  32. Mumford, D., Shah, J., 1989. Optimal approximations by piecewise smooth functions and associated variational problems. Communications on Pure and Applied Mathematics 42 (5), 577–685.
  33. Ortiz, M., Pandolfi, A., 1999. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. International Journal for Numerical Methods in Engineering 44 (9), 1267–1282.
  34. Pham, K., Amor, H., Marigo, J.-J., Maurini, C., 2011. Gradient damage models and their use to approximate brittle fracture. International Journal of Damage Mechanics 20 (4), 618–652.
  35. Pham, K., Marigo, J.-J., Mar 2013. From the onset of damage to rupture: construc