# Numerical investigation of dynamic brittle fracture via gradient damage models

Tianyi Li, Jean-Jacques Marigo, Daniel Guilbaud, Serguei Potapov

► To cite this version:

Tianyi Li, Jean-Jacques Marigo, Daniel Guilbaud, Serguei Potapov. Numerical investigation of dynamic brittle fracture via gradient damage models. Advanced Modeling and Simulation in Engineering Sciences, 2016, 3, pp.26. 10.1186/s40323-016-0080-x. hal-01344553v3

HAL Id: hal-01344553

https://hal.science/hal-01344553v3

Submitted on 30 Aug 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

Open Access

<figure><img src="image_3.png" /></figure>

CrossMark

# Numerical investigation of dynamic brittle fracture via gradient damage models

Tianyi Li<fnref n="1" />* , Jean-Jacques Marigo², Daniel Guilbaud<fnref n="1" />,³ and Serguei Potapov<fnref n="1" />,⁴

**Abstract**

**Background:** Gradient damage models can be acknowledged as a unified framework of dynamic brittle fracture. As a phase-field approach to fracture, they are gaining popularity over the last few years in the computational mechanics community. This paper concentrates on a better understanding of these models. We will highlight their properties during the initiation and propagation phases of defect evolution.

**Methods:** The variational ingredients of the dynamic gradient damage model are recalled. Temporal discretization based on the Newmark-$\beta$ scheme is performed. Several energy release rates in gradient damage models are introduced to bridge the link from damage to fracture.

**Results and discussion:** An antiplane tearing numerical experiment is considered. It is found that the phase-field crack tip is governed by the asymptotic Griffith's law. In the absence of unstable crack propagation, the dynamic gradient damage model converges to the quasi-static one. The defect evolution is in quantitative accordance with the linear elastic fracture mechanics predictions.

**Conclusion:** These numerical experiments provide a justification of the dynamic gradient damage model along with its current implementation, when it is used as a phase-field model for complex real-world dynamic fracture problems.

**Keywords:** Dynamic brittle fracture, Gradient damage models, Griffith's theory, Quasi-static limits

**Background**

Contrary to a sharp interface description of cracks, in the gradient damage approach the introduction of a continuous phase field regularizes displacement discontinuities which are now replaced by strain localizations within a finite band. As can be seen from Fig. 1, cracks are tracked with the help of a scalar damage field $0 \le \alpha \le 1$ which introduces a smooth transition between the undamaged part of the structure where $\alpha = 0$, and the crack where $\alpha = 1$.

The internal length $\ell$ induced naturally through dimensional analysis within the gradient damage model [1] admits several interpretations. From a geometric point of view, it controls the width of the damage band such that a sharp description of cracks can be retrieved in the limit $\ell \to 0$. Meanwhile, it turns out that in such process the total energy in the gradient damage model actually converges in a certain sense to the Griffith functional defined in the variational approach to fracture [2]. The gradient damage model can thus be

<figure><img src="image_4.png" /></figure>

<footnote>*Correspondence:</footnote>

<footnote>tianyi.li@polytechnique.edu</footnote>

<footnote n="1">¹Institute of Mechanical Sciences and Industrial Applications, UMR EDF-CNRS-CEA-ENSTA ParisTech 9219, Université Paris-Saclay, 7 boulevard Gaspard Monge, 91120 Palaiseau, France</footnote>

<footnote>Full list of author information is available at the end of the article</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 1** Phase field description of a discrete crack $\Gamma \subset \Omega$ by a damage field $0 \le \alpha \le 1$. The internal length $\ell$ of the gradient damage model geometrically controls the width of the damage band</figcaption></figure>

regarded as an elliptic regularization of the previous sharp-interface variational fracture model and the internal length $\ell$ serves as a purely numerical parameter which should be as small as possible. This is the interpretation undertaken among others by [3–5] where the gradient damage model is applied to drying, thin films debonding, and other fracture mechanics problems.

On the other hand, the gradient damage model can also be acknowledged as a genuine model per se of brittle fracture, where the internal length $\ell$ is then interpreted as a material parameter which contributes to the fracture or damage behavior of materials. This interpretation presents several advantages from a physical point of view. First of all, this additional length parameter could be related to the maximal stress that the material can sustain and hence introduces additional experimentally validated size effects which are not present in the Griffith model of fracture mechanics, see the work of [6–9] among others. Secondly, the tension-compression asymmetry phenomenon as observed for brittle materials can be easily formulated directly in the gradient damage model. The resulting sharp interface fracture model as $\ell \to 0$ remains unclear and inversely the elliptic regularization of the variational approach to fracture that actually accounts for unilateral contact between crack lips is still considered as a difficult task both from the physical and mathematical point of view, see [10]. Nevertheless, these tension-compression asymmetry formulations as summarized for instance in [11] constitute an improvement of the original gradient damage model [6] and can be regarded as an approximation of the actual non-interpenetration condition.

In this paper we will adhere to the second viewpoint and interpret the gradient damage model as a phase-field approach to dynamic brittle fracture. The formulation of dynamic gradient damage models that extends the original quasi-static ones [1] is sketched in [12]. The governing equations derived from the variational principles resemble those of other phase-field models originated from the computational mechanics community [8, 13–15], with a particular choice of damage constitutive law. These models settle down a unified and coherent numerical framework covering the onset and the space-time propagation of cracks with possible complex topologies and have been successfully applied to study various real-world dynamic fracture problems. Meanwhile, in the quasi-static setting, further physical insights into the gradient damage model seen as a phase-field model of fracture are given among others in [16, 17] where a comparison with the cohesive zone model and the Griffith’s linear elastic theory is conducted. In dynamics, more well-designed numerical experiments should be performed to carry out such verification, see [18] for instance for an investigation of the phase-field crack speeds for plane problems. The main
objective of this contribution is thus to provide a better understanding of the dynamic gradient damage model. We will highlight the properties of gradient damage models as a phase-field approach to fracture and mainly focus on the initiation and propagation phases of defect evolution. For that, an antiplane tearing numerical experiment is considered and an existing crack may initiate and then propagate along a predefined path in a mode-III condition. A comparison with the classical Griffith's theory will be conducted in the cases when the loading speed is of the same order, or small, with respect to the wave speed.

Specifically this paper is organized as follows. The variational formulation of the dynamic gradient damage model is firstly recalled. An open-source numerical implementation is proposed to solve the elastic-damage dynamic problem. Some theoretic concepts are introduced to bridge the link between gradient damage and fracture. We then present and discuss the simulation results corresponding to crack propagation in an antiplane tearing situation. In the last section we will summarize the findings and indicate possible further research directions.

General notation conventions adopted in this paper are summarized here. Scalar-valued quantities will be denoted by italic Roman or Greek letters like the damage field $\alpha_t$. Vectors, second-order tensors and their matrix representations will be represented by boldface letters such as the displacement field $\mathbf{u}_t$ and the stress tensor $\sigma_t$. Higher order tensors considered as linear operators will be indicated by sans-serif letters: the elasticity tensor $A$ for instance. Intrinsic notation is adopted and contraction on lower-order tensors will be written without dots $A\boldsymbol{\varepsilon}_t = A_{ijkl}\boldsymbol{\varepsilon}_{kl}$ (the summation convention is assumed). Inner products between two vectors or tensors of the same order will be denoted with a dot, such as $A\boldsymbol{\varepsilon}_t \cdot \boldsymbol{\varepsilon}_t = A_{ijkl}\boldsymbol{\varepsilon}_{kl}\boldsymbol{\varepsilon}_{ij}$ (the summation convention is assumed). Time dependence will be indicated at the subscripts of the involved quantities, like $\mathbf{u} : (t, \mathbf{x}) \mapsto \mathbf{u}_t(\mathbf{x})$.

# Methods

## Variational framework

The variational ingredients of the dynamic gradient damage model as well as the induced governing equations are firstly recalled in this section. Let us consider a two-dimensional isotropic body $\Omega$. We place ourselves under the small displacement condition. This is a plausible hypothesis for brittle materials when large rotations are also not expected. Damage and subsequent fracture occur more easily in tension than in compression, thus tension-compression asymmetry formulations in the sense of [10, 11] for instance are in general needed. Nevertheless the numerical experiments considered here do not require the use of such formulations. Hence we adhere to the original approach of [1, 19] where damage acts on the sound elastic energy density $\psi_0$ symmetrically under tension and compression. Assuming these hypotheses, the elastic energy of the domain $\Omega$ is given by

$$ 
\begin{aligned}
\mathcal{E}(\mathbf{u}_t, \alpha_t) &= \int_{\Omega} \psi(\boldsymbol{\varepsilon}(\mathbf{u}_t), \alpha_t) \, dx = \int_{\Omega} a(\alpha_t)\psi_0(\boldsymbol{\varepsilon}(\mathbf{u}_t)) \, dx \\
&= \int_{\Omega} \frac{1}{2} a(\alpha_t) A \boldsymbol{\varepsilon}(\mathbf{u}_t) \cdot \boldsymbol{\varepsilon}(\mathbf{u}_t) \, dx
\end{aligned} 
\quad (1) $$

where $A$ is the standard isotropic Hooke's elasticity tensor and $\boldsymbol{\varepsilon}(\mathbf{u}_t) = \frac{1}{2}(\nabla\mathbf{u}_t + \nabla^\mathrm{T}\mathbf{u}_t)$ the linearized strain. By definition, the stress tensor conjugate to the strain measure reads $\sigma_t = a(\alpha_t)A\boldsymbol{\varepsilon}(\mathbf{u}_t)$, with $\alpha \mapsto a(\alpha)$ a non-dimensional damage constitutive function describing stiffness degradation in the bulk. Concerning the kinetic energy, we admit total
mass conservation and no damage dependence of the material density, which leads to its classical definition

$$ \mathcal{K}(\dot{\mathbf{u}}_t) = \int_{\Omega} \kappa(\dot{\mathbf{u}}_t) \mathrm{d}\mathbf{x} = \int_{\Omega} \frac{1}{2} \rho \dot{\mathbf{u}}_t \cdot \dot{\mathbf{u}}_t \mathrm{d}\mathbf{x}. \qquad (2) $$

We now turn to the definition of the *dissipated energy* which quantifies the amount of energy consumed in the damage process. In the phase-field terminology, it is also called the regularized crack functional. It is due to the fact that this energy is closely related to the Griffith-like surface energy of the phase-field representation of cracks according to the Γ-convergence theory [2]. Dynamics should not influence the definition of such functional. We will hence use the same definition used in [7] for quasi-static calculations

$$ S(\alpha_t) = \int_{\Omega} \zeta(\alpha_t, \nabla\alpha_t) \mathrm{d}\mathbf{x} = \int_{\Omega} \frac{G_c}{c_w} \left( \frac{w(\alpha_t)}{\ell} + \ell \nabla\alpha_t \cdot \nabla\alpha_t \right) \mathrm{d}\mathbf{x}. \qquad (3) $$

Contrary to local strain-softening constitutive models, here the damage dissipation mechanism becomes non-local and localization is systematically accompanied by finite energy consumption, due to the presence of the gradient term. In (3), the function $\alpha \mapsto w(\alpha)$ is another non-dimensional damage constitutive law characterizing local damage dissipation. This function along with the former stiffness degradation function $a(\alpha)$ contribute to the damage constitutive behavior of the material and should also satisfy certain physical properties [6] which we do not reproduce here. In this paper we adopt the following constitutive functions

$$ a(\alpha) = (1 - \alpha)^2 \quad \text{and} \quad w(\alpha) = \alpha. \qquad (4) $$

Concerning local damage dissipation, a quadratic function $w(\alpha) = \alpha^2$ originally proposed in [19] is widely used among the phase-field community [8,13–15]. It can be regarded as the Ambrosio and Tortorelli elliptic regularization of the Griffith functional based on their work on image segmentation. In our brittle fracture modeling context however, the use of $w(\alpha) = \alpha$ should be preferred since it guarantees the existence of a purely elastic domain and provides a non-null threshold for damage evolution, see [6]. The link between the gradient damage description of cracks and the Griffith’s one lies in the definition of the normalization factor $c_w$. The fracture toughness $G_c$, i.e. the energy required to create a unit Griffith-like crack surface, can be identified in the gradient damage terminology as the energy dissipated during the optimal damage band creation. Using this identification along with a direct calculation in a 1-d setting, we obtain

$$ c_w = 4 \int_0^1 \sqrt{w(\beta)} \mathrm{d}\beta \qquad (5) $$

which gives $c_w = \frac{8}{3}$ for $w(\alpha) = \alpha$. We refer the readers again to [6,20] for a more detailed discussion on the relationship between these gradient damage and fracture parameters.

In this work, external loads are applied to the body $\Omega$ only through a prescribed displacement $t \mapsto \mathbf{U}_t$ on a subset $\partial\Omega_{\mathrm{U}}$ of the boundary. It will be defined in the admissible displacement space $C_t$ which reads

$$ C_t = \{ \mathbf{u}_t : \Omega \to \mathbb{R}^{\mathrm{dim}} | \mathbf{u}_t = \mathbf{U}_t \text{ on } \partial\Omega_{\mathrm{U}} \}. $$

Damage is modeled as an irreversible defect evolution. Its admissible space will be built from the current damage state $0 \le \alpha_t \le 1$ and it is defined by

$$ D(\alpha_t) = \{ \beta_t : \Omega \to [0, 1] | 0 \le \alpha_t \le \beta_t \le 1 \}. \qquad (6) $$
It can be seen that a damage field $\beta_t$ is admissible, if and only if it is accessible from the current damage state $\alpha_t$ verifying the irreversibility condition, i.e. the damage only grows. In order to formulate the temporal displacement-damage evolution as a boundary value problem (Hamilton's principle), we consider an arbitrary interval of time $I = [0, T]$ and fix the values of $(\mathbf{u}, \alpha)$ at both time ends denoted by $\mathbf{u}_{\partial I} = (\mathbf{u}_0, \mathbf{u}_T)$ and $\alpha_{\partial I} = (\alpha_0, \alpha_T)$. This leads to the following admissible evolution spaces

$$C(\mathbf{u}) = \{ \mathbf{v}: I \times \Omega \to \mathbb{R}^{\text{dim}} \mid \mathbf{v}_t \in C_t \text{ for all } t \in I \text{ and } \mathbf{v}_{\partial I} = \mathbf{u}_{\partial I} \}$$

and

$$D(\alpha) = \{ \beta: I \times \Omega \to [0, 1] \mid \beta_t \in D(\alpha_t) \text{ for all } t \in I \text{ and } \beta_{\partial I} = \alpha_{\partial I} \}.$$

With all the variational ingredients set, we are now in a position to form the space-time action integral given by

$$A(\mathbf{u}, \alpha) = \int_I \mathcal{L}(\mathbf{u}_t, \dot{\mathbf{u}}_t, \alpha_t) \, dt = \int_I (\mathcal{E}(\mathbf{u}_t, \alpha_t) + \mathcal{S}(\alpha_t) - \mathcal{K}(\dot{\mathbf{u}}_t)) \, dt \quad (7)$$

and announce the following three physical principles governing the coupled two-field $(\mathbf{u}, \alpha)$ time-continuous dynamic gradient damage problem:

1 Irreversibility the damage $t \mapsto \alpha_t$ is a non-decreasing function of time.

2 First-order stability the first-order action variation is non-negative with respect to arbitrary admissible displacement and damage evolutions

$$A'(\mathbf{u}, \alpha)(\mathbf{v} - \mathbf{u}, \beta - \alpha) \geq 0 \quad \text{for all } \mathbf{v} \in C(\mathbf{u}) \text{ and all } \beta \in D(\alpha). \quad (8)$$

3 Energy balance the only energy dissipation is due to damage

$$\mathcal{H}_t = \mathcal{H}_0 + \int_0^t \left( \int_{\Omega} (\boldsymbol{\sigma}_s \cdot \boldsymbol{\varepsilon}(\ddot{\mathbf{U}}_s) + \rho \ddot{\mathbf{u}}_s \cdot \ddot{\mathbf{U}}_s) \, dx \right) \, ds \quad (9)$$

where the total energy is defined by

$$\mathcal{H}_t = \mathcal{E}(\mathbf{u}_t, \alpha_t) + \mathcal{S}(\alpha_t) + \mathcal{K}(\dot{\mathbf{u}}_t).$$

**Remark** In the first-order stability condition (8), we evaluate the directional (Gâteaux) derivative of the action integral (7) at $(\mathbf{u}, \alpha)$, a possible solution to the dynamic evolution problem, in the direction $(\mathbf{v} - \mathbf{u}, \beta - \alpha)$ which corresponds to a perturbation. Formally, using the Lagrangian $\mathcal{L}$, we have

$$A'(\mathbf{u}, \alpha)(\mathbf{w}, \beta - \alpha) = \int_I \left( \frac{\partial \mathcal{L}}{\partial \mathbf{u}_t}(\mathbf{s}_t)(\mathbf{w}_t) + \frac{\partial \mathcal{L}}{\partial \dot{\mathbf{u}}_t}(\mathbf{s}_t)(\dot{\mathbf{w}}_t) + \frac{\partial \mathcal{L}}{\partial \alpha_t}(\mathbf{s}_t)(\beta_t - \alpha_t) \right) \, dt,$$

where $\mathbf{w} = \mathbf{v} - \mathbf{u}$ denotes a displacement evolution variation and $\mathbf{s}_t = (\mathbf{u}_t, \dot{\mathbf{u}}_t, \alpha_t)$ corresponds to a state of the dynamical system.

**Remark** Remark that (8) is written as a variational inequality to take into account the unilateral effects introduced by the irreversibility condition in the definition of the damage admissible space (6). It can be regarded as an extension of Hamilton's principle applied to systems with irreversible dissipation. The energy balance condition (9) complements
the first-order stability condition (8) which ensures that energy could only be dissipated through damage (or phase-field like fracture).

By developing the directional derivative of the action integral (7), further physical insights into the first-order stability condition (8) can be obtained if sufficient spatial and temporal regularities of the involved fields are assumed. Testing (8) with $\beta = \alpha$, we obtain after an integration by parts in the time domain

$$ \mathcal{A}'(\mathbf{u}, \alpha)(\mathbf{w}, 0) = \int_I \left( \int_{\Omega} (\sigma_t \cdot \boldsymbol{\varepsilon}(\mathbf{w}_t) + \rho \ddot{\mathbf{u}}_t \cdot \mathbf{w}_t) \, \mathrm{d}\mathbf{x} \right) \mathrm{d}t = 0 \quad \text{for all } \mathbf{w}_t \in \mathcal{C}_0, $$

where $\mathcal{C}_0$ is the associated linear space of $\mathcal{C}_t$, i.e. defined by

$$ \mathcal{C}_0 = \{ \mathbf{u}_t : \Omega \to \mathbb{R}^{\dim|\mathbf{u}_t = 0 \text{ on } \partial\Omega_{II}|} \} $$

This leads thus to the weak elastic-damage dynamic wave equation

$$ \int_{\Omega} (\sigma_t \cdot \boldsymbol{\varepsilon}(\mathbf{w}_t) + \rho \ddot{\mathbf{u}}_t \cdot \mathbf{w}_t) \, \mathrm{d}\mathbf{x} = 0 \quad \text{for all } \mathbf{w}_t \in \mathcal{C}_0. \quad (10) $$

Compared to the classical elastodynamic equation, we note that here the stress tensor is modulated by the stiffness degradation function $\sigma_t = a(\alpha_t)A\boldsymbol{\varepsilon}(\mathbf{u}_t)$.

We now turn to the governing equation for damage evolution induced from the first-order stability condition (8). We observe that the admissible damage space $D(\alpha_t)$ defined in (6) is convex. Due to the arbitrariness of the temporal variation of $\beta$, testing (8) now with $\mathbf{v} = \mathbf{u}$ gives the Euler's inequality condition stating the partial minimality of the total energy with respect to the damage variable under the irreversible constraint for every $t \in I$

$$ \mathcal{E}(\mathbf{u}_t, \alpha_t) + S(\alpha_t) \le \mathcal{E}(\mathbf{u}_t, \beta_t) + S(\beta_t) \quad \text{for all } \beta_t \in D(\alpha_t). \quad (11) $$

Although the same energy minimization principle (11) holds also for quasi-static gradient damage models [6], here the displacement field $\mathbf{u}_t$ is governed by the elastic-damage dynamic wave equation (10). As will be shown through subsequent numerical experiments, this has a direct impact on the apparent crack evolution when damage is propagating along a specific curve. In this work the energy minimization principle (11) will be numerically solved directly at the structural scale by a specific bound-constrained convex optimization algorithm to guarantee the irreversibility condition. The equivalent point-wise conditions of (11) and the energy balance condition (9) can be readily derived by evaluating the inequality $\mathcal{A}'(\mathbf{u}, \alpha)(0, \beta - \alpha) \ge 0$ and performing a temporal derivative of the total energy in (9). This yields a strong formulation in the form of the Kuhn-Tucker conditions which govern local damage evolution at a particular material point. Due to the presence of the damage gradient, the criterion is described by an elliptic type equation in space involving the Laplacian of the damage. Damage growth is not possible until a certain non-local threshold is reached and all energy dissipated in the body during such process until time $t$ is recorded by the crack functional $S(\alpha_t)$. For a detailed derivation of the pointwise conditions, interested readers are referred to [6].

## Numerical implementation

This section describes a numerical implementation of the above continuous two-field evolution problem. In practice it consists of solving numerically the elastic-damage dynamic wave equation (10) coupled with the total energy minimization (11). The irreversibility condition will be automatically enforced during the bound-constrained minimization
process. The time-discrete model which we describe below should converge to the continuous one when the time increment becomes small, see [21]. In particular, the energy balance condition (9) will be hence automatically satisfied.

A spatial discretization is performed based on a mesh $\Omega_h$ of the original domain $\Omega$. The ideal properties of this mesh are indicated in [6] in order to achieve a better modeling of fracture via phase-field approaches. The displacement $\mathbf{u}_t$ and the damage field $\alpha_t$ will be discretized by linear isoparametric finite elements. For two-dimensional plane problems, an arbitrary element possesses at every node 3 nodal degrees of freedom corresponding to 2 components of the displacement and 1 scalar value of the damage. The symbols $\mathbf{u}$ and $\boldsymbol{\alpha}$ are used to denote the current global displacement and damage nodal vectors. Inside a given element $\Omega_e \in \Omega_h$, their local nodal vectors $\mathbf{u}^e$ and $\boldsymbol{\alpha}^e$ are used to perform an interpolation of the displacement and damage fields as well as their derivatives

$$ \mathbf{u}_t(\mathbf{x}) = \mathbf{N}(\mathbf{x})\mathbf{u}^e \quad \text{and} \quad \boldsymbol{\varepsilon}(\mathbf{u}_t)(\mathbf{x}) = \mathbf{B}(\mathbf{x})\mathbf{u}^e, $$

$$ \alpha_t(\mathbf{x}) = \mathbf{N}_{\alpha}(\mathbf{x})\boldsymbol{\alpha}^e \quad \text{and} \quad \nabla\alpha_t(\mathbf{x}) = \mathbf{B}_{\alpha}(\mathbf{x})\boldsymbol{\alpha}^e $$

where $\mathbf{N}$'s and $\mathbf{B}$'s are respectively the interpolation and differentiation matrices.

After spatial discretization the elastic-damage dynamic wave equation (10) becomes

$$ M\ddot{\mathbf{u}} = -\mathbf{F}_{\text{int}}(\mathbf{u}, \boldsymbol{\alpha}) \qquad (12) $$

where $\mathbf{M}$ refers to the classical consistent mass matrix which will be lumped using the traditional row-sum technique described for example in [22]. The internal force vector $\mathbf{F}_{\text{int}}$ is assembled from the elementary vectors given by

$$ \mathbf{F}_{\text{int}}^{e} = \int_{\Omega_e} \mathbf{B}^{\top}\boldsymbol{\sigma}(\mathbf{B}\mathbf{u}^{e}, \mathbf{N}_{\alpha}\boldsymbol{\alpha}^{e}) \, \mathrm{d}\mathbf{x} = \int_{\Omega_e} \mathbf{B}^{\top}(\mathbf{a}(\mathbf{N}_{\alpha}\boldsymbol{\alpha}^{e})\mathbf{A}\mathbf{B}\mathbf{u}^{e}) \, \mathrm{d}\mathbf{x} \qquad (13) $$

where $\mathbf{A}$ is the Voigt representation of the elasticity tensor $\mathbf{A}$.

We now turn to the spatially-discretized damage problem. It can be observed that the use of the damage constitutive law (4) leads to a total damageable energy $\mathcal{E} + \mathcal{S}$ which is a second-order quadratic polynomial with respect to the damage vector $\boldsymbol{\alpha}$. The energy minimization principle (11) involves thus the following functional

$$ q_{\mathbf{u}}(\boldsymbol{\alpha}) = \frac{1}{2} \boldsymbol{\alpha}^{\top} \mathbf{H}(\mathbf{u}) \boldsymbol{\alpha} - \mathbf{b}(\mathbf{u})^{\top} \boldsymbol{\alpha}. \qquad (14) $$

where the Hessian matrix $\mathbf{H}$ and the second member vector $\mathbf{b}$ can be assembled from the elementary matrix and vector given by

$$ \mathbf{H}^{e} = \int_{\Omega_e} (2\psi_0(\mathbf{B}\mathbf{u}^{e})\mathbf{N}_{\alpha}^{\top}\mathbf{N}_{\alpha} + 2w_1\ell^2\mathbf{B}_{\alpha}^{\top}\mathbf{B}_{\alpha}) \, \mathrm{d}\mathbf{x}, $$

$$ \mathbf{b}^{e} = \int_{\Omega_e} (2\psi_0(\mathbf{B}\mathbf{u}^{e}) - w_1) \mathbf{N}_{\alpha} \, \mathrm{d}\mathbf{x} $$

with $w_1 = G_c/(c_w\ell)$. They remain constant during the solving process of the damage problem, since they depend solely on the current displacement state $\mathbf{u}$.

We now consider an arbitrary discretization ($t^n$) of the time interval of interest $I$ where the superscript $n$ denotes a quantity evaluated at the $n$-th time step. We will mainly focus on the time stepping procedures bringing the current known states ($\mathbf{u}^n, \dot{\mathbf{u}}^n, \ddot{\mathbf{u}}^n, \boldsymbol{\alpha}^n$) to the next time step ($\mathbf{u}^{n+1}, \dot{\mathbf{u}}^{n+1}, \ddot{\mathbf{u}}^{n+1}, \boldsymbol{\alpha}^{n+1}$). In the time-continuous model the elastic-damage dynamic wave equation (10) and the damage minimality condition (11) are coupled in the first-order stability principle (8). After temporal discretization $\mathbf{u}$ and $\boldsymbol{\alpha}$ evaluated at the last time step $t = t^n$ and the current time step $t = t^{n+1}$ are in general involved in an implicit
fashion. However, we observe that the energy minimization principle (11) for damage is not a genuine time evolution problem since time dependence is only introduced via the irreversibility condition. In the space-time discrete model at time $t = t^{n+1}$, it reads

$$q_{\underline{u}^{n+1}}(\underline{\alpha}^{n+1}) \le q_{\underline{u}^{n+1}}(\underline{\beta}) \quad \text{for all } \underline{\beta} \text{ that } 0 \le \underline{\alpha}^n \le \underline{\beta} \le 1 \quad (15)$$

where the Hessian matrix and the second member vector in (14) are evaluated at $\underline{u}^{n+1}$. The Eq. (15) can be interpreted as a numerical minimization problem of the quadratic functional $q$ under the irreversible constraint that the current sought damage state $\underline{\alpha}^{n+1}$ is pointwise within the bound $[\underline{\alpha}^n, 1]$

$$\underline{\alpha}^{n+1} = \underset{\underline{\beta}}{\operatorname{argmin}} q_{\underline{u}^{n+1}}(\underline{\beta}) \text{ subjected to the constraints } 0 \le \underline{\alpha}^n \le \underline{\alpha}^{n+1} \le 1. \quad (16)$$

The next damage state $\underline{\alpha}^{n+1}$ can thus be accurately calculated as long as the next deformation state $\underline{u}^{n+1}$ is known. The gradient projection conjugate gradient algorithm initially proposed in [23] is used to solve (16) in an iterative fashion. It is designed for quadratic bound-constrained minimization problems. Due to the bound constraint, approximate solutions $\underline{\beta}$ to (16) can be defined using the projected gradient $[g]$ of which the $i$-th component is given by

$$0 \stackrel{?}{\approx} [g]_i = \begin{cases} \partial_i q & \text{if } \underline{\beta}_i \in (\underline{\alpha}_i^n, 1), \\ \min(\partial_i q, 0) & \text{if } \underline{\beta}_i = \underline{\alpha}_i^n, \\ \max(\partial_i q, 0) & \text{if } \underline{\beta}_i = 1. \end{cases}$$

At each solving iteration, the method consists of several gradient projections to approximately identify the active nodes, i.e. those either $\underline{\alpha}_i^{n+1} = \underline{\alpha}_i^n$ or $\underline{\alpha}_i^{n+1} = 1$. Then it applies the preconditioned conjugate gradient method to minimize an unconstrained reduced problem of the free variables, i.e. those satisfying $\underline{\alpha}_i^n < \underline{\alpha}_i^{n+1} < 1$. The method proceeds to the next iteration until convergence. Interested readers are referred to [23] for a more detailed explanation of the algorithm. The GPCG method is implemented in the parallel linear algebra library PETSc [24]. We also use this library for manipulation of sparse matrices and vectors, similarly to the previous work of [5,6].

It remains then to specify the temporal discretization method used for the $\underline{u}$-problem. In this work we adopt the classical Newmark-$\beta$ integrator, which assumes the following time-stepping procedure

$$\underline{\dot{u}}^{n+1} = \underline{\dot{u}}^n + \frac{\Delta t}{2} (\underline{\ddot{u}}^n + \underline{\ddot{u}}^{n+1}), \quad (17)$$

$$\underline{u}^{n+1} = \underline{u}^n + \Delta t \underline{\dot{u}}^n + \frac{1-2\beta}{2} \Delta t^2 \underline{\ddot{u}}^n + \beta \Delta t^2 \underline{\dot{u}}^{n+1}. \quad (18)$$

The explicit method $\beta = 0$ should be preferred mainly in terms of computational efficiency for applications where the loading speed or the crack propagation speed is comparable to the material speed of sound. The implicit method $0 < \beta \le \frac{1}{2}$ may be suitable for intermediate situations between a quasi-static and an explicit dynamic calculation. From (12), the determination of the new acceleration $\underline{\ddot{u}}^{n+1}$ requires the knowledge of the new deformation state $\underline{u}^{n+1}$ which itself determines the new damage field at time $t = t^{n+1}$ via (16). For the implicit Newmark method $\beta \ne 0$, (18) can thus be regarded as a nonlinear equation in $\underline{u}^{n+1}$, where nonlinearity results from the irreversibility condition when minimizing the total energy (16). To decouple the $(\underline{u}^{n+1}, \underline{\alpha}^{n+1})$ problem, we use a staggered time-stepping procedure as used in [13, 14, 25] among others. The idea is to
update the acceleration $\ddot{\underline{u}}^{n+1}$ while fixing the damage state at its previous known value $\underline{\alpha}^n$. When a relatively small time-step is used, it is expected that the damage increment $\underline{\alpha}^{n+1} - \underline{\alpha}^n$ is bounded and the staggered time-discrete problem will converge to the continuous one, cf. [21]. Introducing the displacement prediction at time $t = t^{n+1}$

$$ \ddot{\underline{u}}^{n+1} = \underline{u}^n + \Delta t \dot{\underline{u}}^n + \frac{1-2\beta}{2} \Delta t^2 \ddot{\underline{u}}^n, $$

from (18) we obtain the linear system for $\underline{u}^{n+1}$

$$ M \frac{\underline{u}^{n+1} - \ddot{\underline{u}}^{n+1}}{\Delta t^2} = -\beta F_{\text{int}}(\underline{u}^{n+1}, \underline{\alpha}^n) = -\beta K(\underline{\alpha}^n) \underline{u}^{n+1} \quad (19) $$

where $K$ is the standard stiffness matrix corresponding to the previous damage state $\underline{\alpha}^n$. The time-stepping procedure for the dynamic gradient damage model based on the implicit Newmark-$\beta$ method in a prediction-correction form is summarized in Algorithm 1.

**Remark** After temporal discretization, the elastic-damage dynamic wave equation (10) and the damage minimality condition (11) can also be solved in a monolithic fashion as described for example in [8,14]. Due to the irreversible condition contained in (16), the GPCG method for instance should be included in the monolithic solver to ensure that the damage variable is subject to a bound constraint during solving. Future work could be devoted to a detailed analysis of these schemes in terms of computational efficiency.

In the explicit case when $\beta = 0$, it turns out that the time evolution system in $(\underline{u}, \underline{\alpha})$ is automatically decoupled and the two subproblems separately in $\underline{u}^{n+1}$ and in $\underline{\alpha}^{n+1}$ can be independently solved one from the other at every time step. Introducing the middle-step velocity

$$ \dot{\underline{u}}^{n+1/2} = \dot{\underline{u}}^n + \frac{\Delta t}{2} \ddot{\underline{u}}^n, $$

the explicit time-stepping procedure for the discretized dynamic gradient damage model is summarized in Algorithm 2.

The initialization phase for the above implicit and explicit time-stepping procedure is described in Algorithm 3. We observe that the initial damage is recomputed $\underline{\alpha}^{-1} \mapsto \underline{\alpha}^0$ in the step 2. The role of $\underline{\alpha}^{-1}$ is to bring some a priori knowledge of the damage field resulting from a previous calculation or more frequently to represent an initial crack $\underline{\alpha}^{-1} = 1$ on

**Algorithm 1** Staggered time-stepping procedure for the dynamic gradient damage model based on the implicit Newmark-$\beta$ method.

1: **for** every successive time step $n \ge 0$ **do**

2: Prediction $\ddot{\underline{u}}^{n+1} = \underline{u}^n + \Delta t \dot{\underline{u}}^n + \frac{1-2\beta}{2} \Delta t^2 \ddot{\underline{u}}^n$.

3: Prediction $\ddot{\underline{v}}^{n+1} = \dot{\underline{u}}^n + \frac{\Delta t}{2} \ddot{\underline{u}}^n$.

4: Solve $\underline{u}^{n+1}$ via (19).

5: Update $\ddot{\underline{u}}^{n+1} = (\underline{u}^{n+1} - \ddot{\underline{u}}^{n+1})/(\beta \Delta t^2)$.

6: Update $\dot{\underline{u}}^{n+1} = \ddot{\underline{v}}^{n+1} + \frac{\Delta t}{2} \ddot{\underline{u}}^{n+1}$.

7: Solve $\underline{\alpha}^{n+1}$ via (16).

8: **end for**
**Algorithm 2** Explicit time-stepping procedure for the dynamic gradient damage model.

1: **for** every successive time step $n \ge 0$ **do**

2: Update $\dot{\underline{u}}^{n+1/2} = \dot{\underline{u}}^n + \frac{\Delta t}{2} \ddot{\underline{u}}^n$.

3: Update $\underline{u}^{n+1} = \underline{u}^n + \Delta t \dot{\underline{u}}^{n+1/2}$.

4: Solve $\underline{\alpha}^{n+1}$ via (16).

5: Update $M\ddot{\underline{u}}^{n+1} = -F_{\text{int}}(\underline{u}^{n+1}, \underline{\alpha}^{n+1})$.

6: Update $\dot{\underline{u}}^{n+1} = \dot{\underline{u}}^{n+1/2} + \frac{\Delta t}{2} \ddot{\underline{u}}^{n+1}$.

7: **end for**

**Algorithm 3** Initialization for the implicit and explicit time-stepping procedure described in Algorithm 1 and 2.

1: Given initial conditions $\underline{u}^0, \dot{\underline{u}}^0$ and $\underline{\alpha}^{-1}$.

2: Reinitialize the damage $\underline{\alpha}^0 = \operatorname{argmin} q_{\underline{u}^0}(\cdot)$ subjected to constraints $0 \le \underline{\alpha}^{-1} \le \underline{\alpha}^0 \le 1$.

3: Initialize the acceleration $M\ddot{\underline{u}}^0 = F_{\text{ext}}^0 - F_{\text{int}}(\underline{u}^0, \underline{\alpha}^0)$.

$\Gamma_0$. The initial step 2 thus renders it compatible with the initial displacement condition and the energy minimization structure.

Both the implicit and explicit time-stepping Algorithms 1 and 2 are implemented as a Python package named "FEniCS Dynamic Gradient Damage". It is based on the FEniCS Project [26] for automated solution of partial differential equations.

## Energy release rate in dynamic gradient damage models

The definition of an energy release rate in gradient damage models which competes with the fracture toughness $G_c$ can be found in [27] under quasi-static conditions. It is found that the damage evolution, when seen as a propagating crack band concentrated along a certain path, is governed by Griffith's law in an asymptotic sense when the internal length $\ell$ is small compared to any other structural length. A theoretic derivation of these similar concepts in dynamics is presented in [28]. Here we will summarize our findings and introduce some useful quantities that establish the link from gradient damage to fracture.

We assume that damage is contained in a thin band described by a parametrized curve $\Gamma_t = \{x \in \Omega | \alpha_t(x) = 1\}$ representing the crack with its current tip $P_t$. Our definition of energy release rates in dynamics is based on shape derivative techniques, see for example [29] for an application of these methods in fracture mechanics. A virtual extension $\theta_t$ of the crack tip in the current propagation direction is introduced. This function $\theta_t$ should verify certain properties discussed in [29]. In particular, we have $\theta_t(P_t) = \tau_t$ where $\tau_t$ refers to the current crack propagation direction. Moreover, It does not alter the crack lip shape, that is $\theta_t \cdot n = 0$ on the crack lip $\Gamma_t$ with $n$ the unit normal vector. A widely used definition of the virtual perturbation is recalled as follows. Suppose that the crack $\Gamma_t$ lies on the $x$-axis and its current crack tip $P_t$ is propagating along the $e_1$ direction. The virtual perturbation $\theta_t$ which introduces a fictive crack advance admits the form $\theta_t = \theta_t e_1$. The construction of the continuous scalar field $0 \le \theta_t \le 1$ parametrized by two radii $r < R$ is given in Fig. 2.
<figure><img src="image_6.png" /><figcaption>**Fig. 2** A particular virtual perturbation $\theta_t = \theta_t e_1$ parametrized by two radii $r < R$. We have $\theta_t = 1$ inside the ball $B_r(\mathbf{P}_t)$, $\theta_t = 0$ outside the ball $B_R(\mathbf{P}_t)$, and a linear interpolation in between</figcaption></figure>

We then perform a reparametrization of the damage field $\alpha_t$ by the arc length $l_t$ of the crack $\Gamma_t$. When $\ell$ is small by comparison with the dimension of the body, the first-order stability condition (8), the energy balance (9) as well as an essential singularity analysis lead to a Griffith-like evolution law governing the phase-field crack length evolution

$$ \dot{l}_t \ge 0, \quad G_t^\alpha - \gamma_t \le 0 \quad \text{and} \quad (G_t^\alpha - \gamma_t) \dot{l}_t = 0. \qquad (20) $$

In (20), the conventional dynamic energy release rate $G_t^\alpha$, which plays the role of “G” in Griffith’s law, is defined by

$$ G_t^\alpha = \int_{\Omega} \left( \begin{aligned} & (\kappa(\dot{\mathbf{u}}_t) - \psi(\boldsymbol{\varepsilon}(\mathbf{u}_t), \alpha_t)) \operatorname{div} \theta_t + \boldsymbol{\sigma}_t \cdot (\nabla \mathbf{u}_t \nabla \theta_t) + \rho \dot{\mathbf{u}}_t \cdot \nabla \mathbf{u}_t \theta_t \\ & + \rho \dot{\mathbf{u}}_t \cdot \nabla \dot{\mathbf{u}}_t \theta_t \end{aligned} \right) \mathrm{d}\mathbf{x}. \qquad (21) $$

It is formally very similar to the classical dynamic energy release rate in linear elastic fracture mechanics (LEFM) written as a volume integral with the help of the virtual perturbation $\theta_t$, see for example [30]. Note however that in the gradient damage model the elastic energy density $\psi$ and the stress tensor $\boldsymbol{\sigma}_t$ are modulated by the stiffness degradation function $a(\alpha)$ and consequently stress singularity automatically disappears at the crack tip. Written in the form of (21), the conventional dynamic energy release rate involves an integral in the elements and hence is more convenient and accurate compared to a line integral (e.g. the J-integral) in a finite element calculation. The use of the traditional J-integral in the sense of [31] to calculate an effective energy release rate in a gradient damage modeling of fracture can be found in [17,20] for instance.

The quantity $\gamma_t$ in (20) which plays the role of “$G_c$” in the classical Griffith’s law is the damage dissipation rate

$$ \gamma_t = \int_{\Omega} (\zeta(\alpha_t, \nabla\alpha_t) \operatorname{div} \theta_t - 2w_1 \ell^2 \nabla\alpha_t \cdot \nabla\theta_t \nabla\alpha_t) \mathrm{d}\mathbf{x}, \qquad (22) $$

where $w_1 = G_c/(c_w \ell)$. Formally it is defined as the derivative of the dissipated energy (3) with respect to the crack length $l_t$, thus quantifying the energy dissipated due to damage per crack advance.

**Remark** The conventional dynamic energy release rate (21) and the damage dissipation rate (22) admit also a J-like path integral representation involving a generalized Eshelby tensor used frequently in configurational force approaches such as [32,33] and references therein. During crack propagation, we have
$$ \hat{J}_t = \lim_{\epsilon \to 0} \int_{C_\epsilon} \hat{J}_t \mathbf{n} \cdot \boldsymbol{\tau}_t \, ds = G_t^\alpha - \gamma_t, $$

where $\mathbf{n}$ denotes the normal pointing out of the ball of radius $\epsilon$ centered at the crack tip $B_\epsilon(\mathbf{P}_t)$ with $C_\epsilon = \partial B_\epsilon(\mathbf{P}_t)$ its boundary. The generalized dynamic $\hat{J}_t$ tensor is defined by

$$ \hat{J}_t = (\psi(\boldsymbol{\epsilon}(\mathbf{u}_t), \alpha_t) + \kappa(\dot{\mathbf{u}}_t) + \zeta(\alpha_t, \nabla\alpha_t)) \mathbb{I} - \nabla \mathbf{u}_t^\top \boldsymbol{\sigma}_t - 2w_1 \ell^2 \nabla\alpha_t \otimes \nabla\alpha_t, $$

where $w_1 = G_c/(c_w \ell)$. Interested readers are referred to [28] for a detailed discussion on this point.

With the help of a two-scale approach, the inner damage problem near the crack tip and the outer LEFM problem far from the crack band can be separated [27]. It can be shown that if the inner radius $r$ of the virtual perturbation defined in Fig. 2 is sufficiently big with respect to the internal length $\ell$, $G_t^\alpha$ defines an equivalent dynamic energy release rate $G_t$ corresponding to the outer mechanical fields. Similarly, the damage dissipation rate $\gamma_t$ will converge to the fracture toughness $G_c$ defined in the crack functional (3). The following asymptotic Griffith's law is obtained when the outer fields are considered as $r$ increases

$$ \dot{l}_t \ge 0, \quad G_t - G_c \le 0 \quad \text{and} \quad (G_t - G_c) \dot{l}_t = 0. \qquad (23) $$

## Results and discussion

In this section we will present and discuss a particular numerical experiment tailored to highlight the properties of the dynamic gradient damage model while focusing on the initiation and propagation phases of defect evolution. Specifically, we will investigate the fracture mechanics criterion for an existing phase-field crack to initiate, and then to propagate along a certain path.

We consider a mode-III antiplane tearing of a two dimensional plate $\Omega = (0, L) \times (-H, H)$ subject to a hard device $\mathbf{U}_t = \operatorname{sgn}(y)kt \mathbf{e}_3$ on its left border $x = 0$, see Fig. 3.

An initial damage field corresponding to an preexisting crack $\Gamma_0 = \{x \in \mathbb{R}^2 | \alpha_0(x) = 1\} = [0, l_0] \times \{0\}$ is present in the domain. For that we prescribe naturally $\underline{\alpha}^{-1} = 1$ on $\Gamma_0$ in Algorithm 3. The loading velocity $k$ will be varied and its effect on the crack propagation speed will be studied. With a modification of the damage dependence of the elastic energy $\mathcal{E}(\mathbf{u}_t, \alpha_t)$ proposed in [25], the crack tip $t \mapsto \mathbf{P}_t$ is enforced to propagate

<figure><img src="image_7.png" /><figcaption>**Fig. 3** Mode III antiplane tearing of a two dimensional plate $\Omega = (0, L) \times (-H, H)$ with a loading speed parametrized by $k$. An initial crack $[0, l_0] \times \{0\}$ is present in the domain. The crack is enforced to propagate along the constant direction $\mathbf{e}_1$</figcaption></figure>
along the constant direction $e_1$, which prohibits consequently crack kinking or branching. Denoting $u$ as the out-of-plane displacement, the elastic energy density in this particular situation reads

$$ \psi(\nabla u_t, \alpha_t) = \frac{1}{2}\mu \left( a(\alpha_t) \left( \frac{\partial u_t}{\partial x_2} \right)^2 + \left( \frac{\partial u_t}{\partial x_1} \right)^2 \right). \qquad (24) $$

where damage acts only on the $\frac{\partial u_t}{\partial x_2}$ component of the displacement gradient. We assume that this modification (24) can be regarded as a particular case of the original model (1) when the crack actually propagates along a specific path without kinking or branching. It allows us to focus on the crack propagation stage.

*Remark* If the original model is used, i.e. when the degradation function also acts on $\frac{\partial u_t}{\partial x_1}$, numerically it is observed that for low propagation speeds crack curving (including kinking and branching) does not take place and the modification (24) produces the same response as the original model. However for higher propagation speeds (for example due to a larger loading velocity $k$), crack curving is observed (see for example [25]) and these two models no longer predict the same crack evolution. Crack path prediction is exactly the *raison d'être* of phase-field models of fracture. A thorough investigation of crack kinking/branching phenomena (as a function of crack speed for instance) is a very important task to which future work will be devoted. Nevertheless, the current contribution focuses on the behavior of gradient damage models when these dynamic instabilities (kinking, branching) are somehow suppressed (see for example [34] for an experimental investigation on this point), which permits a direct comparison with the classical Griffith's theory of dynamic fracture.

This problem is initially raised in [25]. In their model the crack surface energy is approximated by the Ambrosio and Tortorelli elliptic regularization $w(\alpha) = \alpha^2$, whereas here the damage constitutive law (4) is used. The objective is to compare the crack evolution obtained in the dynamic gradient damage model with that predicted by Griffith's law which determines initiation and propagation of cracks. Two experiments will be considered:

1. In the first case, the fracture toughness $G_c$ is assumed to be homogeneous throughout the domain. The loading speed is of the same order of the material speed of sound $c = \sqrt{\mu/\rho}$ and we will use the explicit Newmark time-stepping method.

2. In the second case, $G_c$ may admit a spatial discontinuity in the propagation direction. We also prescribe a relatively small loading speed in order to investigate the quasi-static limit of the dynamic model. Depending on whether the crack propagation speed itself is smaller with respect to the speed of sound or not (the term *unstable* propagation often refers to this case), the implicit or the explicit Newmark method will be used.

A rescaling of the displacement and a normalization of the space/time scales are performed to obtain a non-dimensional problem. Specifically, a reference elastic constant $\bar{\mu}$, material density $\bar{\rho}$ and fracture toughness $\bar{G}_c$ have been chosen and the displacement is scaled by a factor of

$$ \bar{u} = \sqrt{\bar{G}_c H / \bar{\mu}}. $$
The height of the plate *H* and the corresponding time *H*/*c* for the reference elastic wave (with speed *c* = √μ/ρ) to travel such distance is used to normalize the spatial and temporal scales. We have thus

$$u(x, t) = \bar{u} \hat{u} \left( \frac{x}{H}, \frac{\bar{c}t}{H} \right) e_3,$$

$$\alpha(x, t) = \hat{\alpha} \left( \frac{x}{H}, \frac{\bar{c}t}{H} \right).$$

Rewriting Lagrangian defined in (7) using $\hat{u}$ and $\hat{\alpha}$ amounts to adopt the following non-dimensional quantities

$$\hat{\rho} = \frac{\rho}{\bar{\rho}}, \quad \hat{\mu} = \frac{\mu}{\bar{\mu}}, \quad \hat{G_c} = \frac{G_c}{\bar{G_c}}, \quad \text{and} \quad \hat{\ell} = \frac{\ell}{H}.$$

For notational simplicity, we drop the bar and use directly non-dimensional quantities in the sequel.

A structured crossed triangular mesh with a uniform discretization spacing $\Delta x = \Delta y = h$ is generated. For the explicit time-stepping method, the Courant–Friedrichs–Lewy (CFL) time-step is used

$$\Delta t_{\text{CFL}} = \frac{h}{c} = \frac{h}{\sqrt{\mu/\rho}}. \qquad (25)$$

The parameters adopted for all subsequent calculations are summarized in Table 1. A typical damage field obtained in this simulation is illustrated in Fig. 4, where the damage varies from 0 (blue zones) to 1 (red zones). Thanks to the Γ-convergence theory in the variational approach to fracture, the current crack length $l_t$ could be approximately derived from the damage dissipation energy

$$S(\alpha_t) \approx (G_c)_{\text{eff}} l_t \qquad (26)$$

with $(G_c)_{\text{eff}} = (1 + 3h/(8\ell))G_c$ the numerical amplified fracture toughness due to spatial discretization, see [2, 20]. This is the value that a phase-field crack actually dissipates per unit extension in the gradient damage model. However (26) does not immediately apply to the case where $G_c$ admits a spatial discontinuity. For consistency, the current crack tip $P_t = (l_t, 0)$ is located on the contour $\alpha = 0.5$. The crack speed can thus be obtained by a linear regression analysis during the steady propagation phase.

<table><caption>**Table 1** Geometric, material and numerical parameters for the antiplane tearing experiment</caption><thead><tr><th>L</th><th>H</th><th>l<sub>0</sub></th><th>μ</th><th>ρ</th><th>G<sub>c</sub></th><th>ℓ</th><th>h</th><th>Δt</th></tr></thead><tbody><tr><td>5</td><td>1</td><td>1</td><td>0.2</td><td>1</td><td>0.01</td><td>0.05</td><td>0.01</td><td>Δt<sub>CFL</sub></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 4** Typical damage field obtained in the antiplane tearing example. The damage varies from 0 (blue zones) to 1 (red zones)</figcaption></figure>
## Homogeneous fracture toughness case

In the first case a homogeneous plate will be considered. This antiplane tearing example is physically similar to the 1-d film peeling problem which can be studied using the classical Griffith's theory of dynamic fracture. According to [35] and [2], the crack speed, with respect to the loading displacement $U = kt$ or to the physical time $t$, as a function of the loading velocity $k$ is given by

$$ \frac{dl}{dU}(k) = \sqrt{\frac{\mu H}{G_c + \rho H k^2}} \quad \text{or} \quad \frac{dl}{dt}(k) = \sqrt{\frac{\mu H k^2}{G_c + \rho H k^2}} \qquad (27) $$

from which we retrieve the quasi-static limit $dl/dU(0) = \sqrt{\mu H/G_c}$ predicted in [2] and the dynamic limit as the shear wave speed $dl/dt(\infty) = \sqrt{\mu/\rho}$, which is a classical result in the Griffith's theory of dynamic fracture [36]. We also observe that for low loading speeds $k \approx 0$, the dynamic crack speed $dl/dt \approx k\sqrt{\mu H/G_c}$ scales linearly in $k$, which agrees with the remarks given in [25]. Comparisons between the numerical results using the dynamic gradient model and this theoretic result (27) with $G_c$ replaced by $(G_c)_{eff}$ are illustrated in Fig. 5.

Despite the transverse wave reflection present in the two-dimensional numerical model, a very good quantitative agreement is found between them. In particular, as it is also observed in [25], the numerically obtained crack speed indeed approaches the limiting shear wave speed when the loading speed increases. The explicit time-stepping Algorithm 2 as well as its implementation work fine even at supersonic loading speeds.

The conventional dynamic energy release rate (21) is numerically computed and the validity of the asymptotic Griffith's law is analyzed by varying the inner radius $r$ of virtual perturbations defined in Fig. 2. During the propagation phase $\dot{l}_t > 0$, three arbitrary time instants are taken when the crack length attains respectively $l_t \approx 1.6$, $l_t \approx 2$ and $l_t \approx 2.4$. An evident $r$-dependence of $G_t^\alpha$ is illustrated in Fig. 6, where the ratio $R/r = \frac{5}{2}$ is fixed.

In LEFM, it is known that the $J$-integral is directly related to the dynamic stress intensity factors at the crack tip, see [36]. In gradient damage models however, there is no more stress singularities. When $r$ is small, we go directly into the process zone dominated by damage-induced strain softening and $G_t^\alpha \to 0$ is expected as $r \to 0$. However, as $r$

<figure><img src="image_9.png" /><figcaption>**Fig. 5** Crack speeds as a function of the loading velocity. The crack speed $\frac{dl}{dt}$ with respect to $t$ is indicated in **a**, while in **b** the crack speed $\frac{dl}{dU}$ with respect to $U$ is shown. Comparison with the 1-d analytical solution (27) based on the Griffith's criterion</figcaption></figure>

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 6** Conventional dynamic energy release rate $G^\alpha$ as a function of the inner radius $r$ of the virtual perturbation $\theta_t$. Three arbitrary instants when the crack propagates $l_t > 0$ are chosen.</figcaption></figure>

increases, $G_t^\alpha$ captures well the outer mechanical fields. An equivalent energy release rate can thus be defined, and according to the asymptotic behaviors of $\gamma_t$, we have the desired result $G_t^\alpha = \gamma_t \rightarrow (G_c)_{\text{eff}}$.

We will then turn to the evolution of the conventional dynamic energy release rate when the existing crack initiates and further propagates. From the above $r$-dependence analysis, a fixed inner radius $r = 2\ell$ is used which should already correctly capture the far mechanical fields. The crack length $l_t$ given by (26) as well as the calculated $G_t^\alpha$ are given as a function of the loading displacement in Fig. 7, where three separate calculations corresponding to three loading speeds $k$ are reported.

Recall that an initial crack of length 1 is present in the body and we observe $G_t^\alpha = 0$ before the waves arrive at the initial crack tip. When the energy release rate $G_t^\alpha$ at the initial crack tip attains the fracture toughness $(G_c)_{\text{eff}}$, the existing crack initiates and then propagates with the equality $G_t^\alpha = (G_c)_{\text{eff}}$ if the spatial and temporal numerical discretization errors are ignored. Indeed this equality is not enforced algorithmically during the solving of the $(u, \alpha)$ evolution which is instead determined by Algorithm 2. We may conclude that

<figure><img src="image_12.png" /><figcaption>**Fig. 7** Conventional dynamic energy release rate $G^\alpha$ as a function of the loading displacement. Three loading speeds $k$ are used: $k = 0.1 \approx 0.2c$, $k = 0.2 \approx 0.4c$ and $k = 0.3 \approx 0.7c$</figcaption></figure>
the crack-tip evolution (initiation and propagation) is well governed by the asymptotic Griffith's law (23) in the dynamic gradient damage model, when outer fields are considered.

The internal length $\ell$ plays a rather subtle role during the propagation phase. The crack tip is governed by the asymptotic Griffith's law (23) if and only if a separation of scales between the inner damage problem and the outer LEFM is possible, i.e. only when the internal length is sufficiently small compared to any other structural length [27]. Although $\ell$ is indeed hidden in (23), the validity of the latter depends directly on it. Below we present the simulation results with a fixed loading speed $k = 0.2$ and three small enough internal lengths. As can be seen from Fig. 8, the crack evolution is globally conforming with Griffith's law, as long as the involved quantities are calculated with a virtual perturbation $\theta_t$ capturing correctly the far fields. Here according to Fig. 6, we use an inner radius adapted with the internal length $r = 2\ell$, which should produce an error less than 3%.

The stress distribution along a vertical slice $\{(x, y) \in \mathbb{R}^2 | x = l_t\}$ passing by the current crack tip $P_t$ should illustrate and highlight the separation of scales when $\ell$ is small. For the sake of simplicity, we consider a stationary crack $[0, 2] \times \{0\}$ and solve the static problem with the gradient damage model and the LEFM model (linear elastic body with a sharp crack embedded in the domain). We can verify from Fig. 9 that the LEFM develops a

<figure><img src="image_13.png" /><figcaption>**Fig. 8** Crack evolution as a function of the loading displacement. Three small enough internal lengths are used</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 9** Stress distribution along a vertical slice $\{(x, y) \in \mathbb{R}^2 | x = l_t\}$ passing by the current crack tip $P_t$. The $\sigma_{13}$ (a) and the $\sigma_{23}$ (b) components are indicated. The gradient damage model and the LEFM model are compared</figcaption></figure>

<figure><img src="image_15.png" /></figure>
well-known inverse square root singularity for the two stress components $\sigma_{13}$ and $\sigma_{23}$ and their near-tip fields are well approximated by the theoretic asymptotic solutions.

On the other hand, the gradient damage model provides a better modeling of the stress field near the crack tip as their values are bounded. A good matching can be observed far from the crack tip and the discrepancy with the outer LEFM model is concentrated within a process zone proportionally dependent on the internal length. When $\ell$ is very large, the process zone could cover the whole structural domain and a separation of scales is no longer possible. In this case the asymptotic Griffith's law (23) is not applicable since we are no longer dealing with a fracture mechanics problem.

## Discontinuous fracture toughness cases

As another illustration of the Griffith-conforming crack evolution obtained with the dynamic gradient damage model, we consider the quasi-static limits of the model in the presence of a possible fracture toughness discontinuity in the previous plate

$$G_c = \begin{cases} K_1 & x \le x_0, \\ K_2 & x > x_0. \end{cases}$$

A preexisting crack is always present and is introduced via an initial damage field. According to [37], when the loading speed $k$ is decreased the dynamic gradient damage model converges to the following first-order quasi-static gradient damage evolution model:

1. Irreversibility the damage $t \mapsto \alpha_t$ is a non-decreasing function of time.

2. First-order stability the first-order variation of the potential energy is non-negative with respect to arbitrary admissible displacement and damage fields

$$\mathcal{P}'(\mathbf{u}_t, \alpha_t)(\mathbf{v}_t - \mathbf{u}_t, \beta_t - \alpha_t) \geq 0 \quad \text{for all} \quad \mathbf{v}_t \in C_t \text{ and all} \quad \beta_t \in D(\alpha_t). \quad (28)$$

where in the absence of external forces the potential energy is given by

$$\mathcal{P}(\mathbf{u}_t, \alpha_t) = \mathcal{E}(\mathbf{u}_t, \alpha_t) + \mathcal{S}(\alpha_t)$$

3. Energy balance the only energy dissipation is due to damage

$$\mathcal{P}_t = \mathcal{P}_0 + \int_0^t \left( \int_{\Omega} \boldsymbol{\sigma}_s \cdot \boldsymbol{\varepsilon}(\dot{\mathbf{U}}_s) \, \mathrm{d}\mathbf{x} \right) \, \mathrm{d}s. \quad (29)$$

In addition, the proof of convergence is made under the hypothesis that the crack evolution $t \mapsto l_t$ is at least continuous in time (as in the classical Griffith theory). Numerically, it is the first-order stability condition (28) that is effectively implemented by the alternate minimization procedure [6] while the energy balance condition (29) can only be at best checked a posteriori.

The homogeneous antiplane tearing problem is firstly solved by the dynamic gradient damage model and the above first-order quasi-static gradient damage model. In the dynamic calculation a small loading speed $k = 0.001 \approx 0.2\%c$ is assumed and we use the unconditionally stable implicit Newmark scheme as described in Algorithm 1, with $\beta = \frac{1}{4}$. The time step is set to $\Delta t = 10\Delta t_{CFL}$, where $t_{CFL}$ is defined in (25). In Fig. 10 we plot the crack length evolution as well as the conventional energy release rate $G_t^\alpha$ both for the dynamic model and the first-order quasi-static model. It is recalled that the static $G_t^\alpha$ can be simply obtained by setting $\dot{\mathbf{u}}_t$ and $\ddot{\mathbf{u}}_t$ to zero in (21). We observe that these two solutions coincide, and both present a time-continuous crack evolution (initiation
<figure><img src="image_16.png" /><figcaption>**Fig. 10** Crack length and conventional energy release rate $G^\alpha$ for the homogeneous fracture toughness plate at a very slow loading speed. Comparison between the dynamic model and the first-order quasi-static model</figcaption></figure>

and propagation) conforming to the asymptotic Griffith's law (23). The numerically computed quasi-static crack speed (with respect to $U = kt$) is compared in Table 2 to the analytical value $\sqrt{\mu H/G_c}$ announced in [2]. A very good agreement can be found if the numerically amplified fracture toughness $(G_c)_{eff}$ is used in the formula.

We then turn to the case where the fracture toughness jumps suddenly from a lower value $K_1 = 0.01$ to a higher one $K_2 = 2K_1 = 0.02$ at $x = 2$. The unconditionally stable implicit Newmark scheme with $\beta = \frac{1}{4}$ is used again with a time increment $\Delta t = 10\Delta t_{CFL}$. As can be observed from Fig. 11 the convergence of the dynamic model toward the quasi-static one is verified and the crack initiates and propagates following Griffith's law. A temporary arrest phase is present shortly after the crack reaches the interface at $x = 2$. Due to continuous loading the energy release rate increases and the crack then restarts

<table><caption>**Table 2** Comparison of the numerically computed quasi-static crack speed in the homogeneous case with the theoretic one $\sqrt{\mu H/G_c}$ given in [2]</caption><thead><tr><td></td><td>Numerical</td><td>Theoretic</td><td>Error</td></tr></thead><tbody><tr><td>Quasi-static crack speed</td><td>4.326</td><td>4.391</td><td>1.5%</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>**Fig. 11** Crack length and conventional energy release rate $G^\alpha$ for the increasing fracture toughness plate at a very slow loading speed. Comparison between the dynamic model and the first-order quasi-static model. The numerically amplified fracture toughness $(G_c)_{eff}$ is calculated based on $K_1 = 0.01$</figcaption></figure>
and begins to propagate in the second material when the energy release rate $G_t^\alpha$ attains the higher fracture toughness $K_2$.

However, for the case where the fracture toughness $K_1 = 2K_2 = 0.02$ suddenly drops to a smaller value $K_2 = 0.01$ at $x = 1$ (exceptionally here the initial crack length is $1/4$), a relatively good matching can only be found before and after the jump phase produced at the discontinuity, both in terms of the crack length evolution and the energy release rate. It is exactly at the jump phase that these two models strongly disagree, cf. Fig. 12.

Here due to the unstable crack propagation during the jump, the explicit Newmark scheme is used for the dynamic calculation with $\Delta t = \Delta t_{CFL}$. When the crack arrives at the discontinuity, the *first-order* quasi-static numerical model underestimates the crack jump and predicts no further crack arrest, by relating directly the static energy release rate $G_t^\alpha$ to the fracture toughness $K_2$ just after the jump. For the dynamic model, the jump length is bigger and a subsequent temporary crack arrest is observed, as the dynamic energy release rate oscillates with a high frequency but remains smaller than the fracture toughness $K_2$ after the jump. We observe that in both cases the jump takes place at $x \approx 0.9$ somewhat prior to the fracture toughness discontinuity $x = 1$. We suspect that this is due to the damage regularization of cracks with a half-band $D = 2\ell = 0.1$ using the constitutive laws of (4). If this effect is ignored, the crack length after the jump is recorded in Table 3 for each case. From the static energy release rate evolution, we see that the crack length $l_m$ after the jump predicted in the *first-order* quasi-static numerical model is governed by $G(l_m) = G_c(l_m)$ from which authors of [35] find $l_m = \sqrt{K_1/K_2} = \sqrt{2}$. However their dynamic analysis shows that the crack length after the jump $l_c$ should instead be given by the total (quasi-static) energy conservation principle $P(1) = P(l_c)$, which results in $l_c = K_1/K_2 = 2$. We see from Table 3 that our dynamic gradient damage model indeed reproduces this correct value.

<figure><img src="image_18.png" /><figcaption>**Fig. 12** crack length and the conventional energy release rate $G^\alpha$ for the decreasing fracture toughness plate at a very slow loading speed. Comparison between the dynamic model and the first-order quasi-static model. The numerically amplified fracture toughness $(G_c)_{eff}$ is calculated based on $K_2 = 0.01$</figcaption></figure>

<table><caption>**Table 3** Comparison of the numerical crack lengths after the jump with the theoretic predictions</caption><thead><tr><td></td><td>Quasi-static</td><td>Dynamic</td></tr></thead><tbody><tr><td>Numerical</td><td>1.465</td><td>1.995</td></tr><tr><td>Theoretic</td><td>√2</td><td>2</td></tr><tr><td>Error</td><td>3.6%</td><td>0.25%</td></tr></tbody></table>
To better analyze the jump phase, energy evolutions are investigated against the crack length in Fig. 13.

In the quasi-static case we pick the total energy $P = \mathcal{E} + S$ while in the dynamic case we plot separately the static energy $P = \mathcal{E} + S$ and the kinetic one $K$. Before and sufficiently after the jump a good agreement between these two potential energies can be found. We observe that the (incorrect) quasi-static jump i.e., an unstable or brutal crack propagation) is accompanied by a slight loss of the total energy $\Delta P_{\text{stat}}$, contradicting the balance condition (29). This phenomenon has already been observed by several authors such as [2, 6, 10, 25]. On the one hand, it can be regarded as a numerical issue as the effective implementation of the quasi-static model is solely based on the first-order stability condition (28). For this particular problem based on quasi-static energy conservation we could predict a correct quasi-static crack evolution toward which the dynamic solution converges when the loading speed becomes small, see [35]. On the other hand, from a theoretic point of view, it is already known in [38] that there may not exist an energy-conserving evolution which also respects the stability criterion at every time. Moreover even equipped with the energy balance condition, the quasi-static model may still differ from the dynamic analysis [39]. A natural and physical remedy for all general unstable crack propagation cases is to introduce inertial effects. In Fig. 13 the dynamic jump process is continuous (the crack propagates at a finite speed bounded by the shear wave speed) compared to the quasi-static one where the jump occurs necessarily in a discontinuous fashion between two iterations. We verify the conclusions drawn in [35] that the kinetic energy $K$ plays only a transient role in this problem, as it attains a finite value during the jump and becomes again negligible after. The dynamic potential energy $P = \mathcal{E} + S$ after the jump is slightly bigger that its value before the jump, due to the fact that the loading speed $k = 0.001$ is small but not zero.

During the jump, the crack propagates at a speed comparable to the material speed of sound which, according to [35], is given by

$$v_{\text{jump}} = \frac{\left(\sqrt{\hat{K}_1 + \epsilon^2 + \epsilon}\right)^2 - \hat{K}_2}{\left(\sqrt{\hat{K}_1 + \epsilon^2 + \epsilon}\right)^2 + \hat{K}_2} \cdot c \quad (30)$$

<figure><img src="image_19.png" /><figcaption>Fig. 13 Energy variation as a function of the crack length for the decreasing fracture toughness plate at a very slow loading speed. Comparison between the dynamic model and the first-order quasi-static model</figcaption></figure>
with the non-dimensional fracture toughness $\hat{K}_i = K_i/(2\mu H)$ and the normalized loading speed $\epsilon = k/c$. The crack length evolution during the jump is illustrated in Fig. 14.

Due to transverse wave reflection in this 2-d problem, the crack propagates during this interval with a small fluctuation of period $T$ approximately corresponding to the first standing wave between the boundary and the crack $T \approx 2H/c \approx 4.5$. That is why we calculate from Fig. 14 only the initial crack speed at jump for comparison in Table 4. A good agreement can be found between the numerical and the theoretic ones.

## Conclusion

In this paper further physical insights into the dynamic gradient damage model are provided via a simple antiplane tearing experiment. As a phase-field approach to brittle fracture, it can indeed be regarded as a generalization or a superset of the LEFM theory, since the crack evolution is shown to be Griffith-conforming in several situations:

* In the dynamic tearing example of a homogeneous plate, it is verified that the crack evolution is governed by the asymptotic Griffith's law (23), as long as the material internal length is sufficiently small to establish a separation of scales between the inner damage problem and the outer LEFM problem. The conventional dynamic energy release rate is numerically computed and verified as a tool to translate gradient damage mechanics results in fracture mechanics terminology. We conduct a comparison with the 1-d peeling problem [35] analytically studied with the classical Griffith's theory of dynamic fracture. A good agreement between them can be found in terms of the crack speeds prediction as a function of the loading speed.

* We then investigate the quasi-static limits of the dynamic gradient damage model. In the absence of *brutal* or *unstable* crack propagation where the classical static Griffith's theory fails, the dynamic model converges to the first-order quasi-static gradient damage model, when the loading speed decreases. However, when the crack

<figure><img src="image_20.png" /><figcaption>Fig. 14 Zoom in time at the crack length jump due to sudden toughness decrease for the dynamic model</figcaption></figure>

<table><caption>**Table 4** Comparison of the numerically computed crack jump speed with the theoretic one (30) given in [35]</caption><thead><tr><td></td><td>Numerical</td><td>Theoretic</td><td>Error</td></tr></thead><tbody><tr><td>Relative jump speed v<sub>jump</sub>/c</td><td>0.3325</td><td>0.3396</td><td>2%</td></tr></tbody></table>
may propagate at a speed comparable to the material speed of sound, the dynamic model should be preferred to correctly account for inertial effects. The crack evolution in the dynamic gradient damage model is in quantitative accordance with the LEFM predictions on the 1-d peeling problem.

These numerical experiments provide hence a justification of the dynamic gradient damage model along with its current implementation, when it is used as a phase-field model for complex real-world dynamic fracture problems. Furthermore, the gradient damage model could also be used as a tool to explore and explain numerous dynamic fracture phenomena. Our future work will be devoted to this point.

**Abbreviations**

LEFM: linear elastic fracture mechanics; GPCG: gradient projection conjugate gradient algorithm.

**Authors' contributions**

TL and JJM have developed the theory part. TL has carried out the numerical implementation and has drafted the manuscript. JJM, DG and SP have supervised the different studies and the corrections of the draft. All authors read and approved the final manuscript.

**Author details**

¹Institute of Mechanical Sciences and Industrial Applications, UMR EDF-CNRS-CEA-ENSTA ParisTech 9219, Université Paris-Saclay, 7 boulevard Gaspard Monge, 91120 Palaiseau, France, ²Laboratoire de Mécanique des Solides, École Polytechnique, 91128 Palaiseau Cedex, France, ³DEN-Service d'Etudes Mécaniques et Thermiques (SEMT), CEA Saclay, 91191 Gif-sur-Yvette Cedex, France, ⁴EDF Lab Paris-Saclay, 7 boulevard Gaspard Monge, 91120 Palaiseau, France.

**Competing interests**

The authors declare that they have no competing interests.

**Availability of data and materials**

The datasets supporting the conclusions of this article are available in the FEniCS Dynamic Gradient Damage repository, https://bitbucket.org/litianyi/dynamic-gradient-damage.

Received: 11 May 2016 Accepted: 12 August 2016

Published online: 30 August 2016

**References**

1. Pham K, Marigo J-J. Approche variationnelle de l'endommagement : II. Les modèles à gradient. Comptes Rendus Mécanique. 2010;338(4):199-206.

2. Bourdin B, Francfort GA, Marigo J-J. The variational approach to fracture. J Elast. 2008;91(1-3):5-148.

3. Maurini C, Bourdin B, Gauthier G, Lazarus V. Crack patterns obtained by unidirectional drying of a colloidal suspension in a capillary tube: experiments and numerical simulations using a two-dimensional variational approach. Int J Fract. 2013;184(1-2):75-91. doi:10.1007/s10704-013-9824-5.

4. Baldelli AAL, Babadjian JF, Bourdin B, Henao D, Maurini C. A variational model for fracture and debonding of thin films under in-plane loadings. J Mech Phys Solids. 2014;70:320-48.

5. Mesgarnejad A, Bourdin B, Khonsari MM. Validation simulations for the variational approach to fracture. Comput Methods Appl Mech Eng. 2015;290:420-37. doi:10.1016/j.cma.2014.10.052.

6. Pham K, Amor H, Marigo J-J, Maurini C. Gradient damage models and their use to approximate brittle fracture. Int J Damage Mech. 2011;20(4):618-52.

7. Bourdin B, Marigo J-J, Maurini C, Sicsic P. Morphogenesis and propagation of complex cracks induced by thermal shocks. Phys Rev Lett. 2014;112(1):014301.

8. Schlüter A, Willenbücher A, Kuhn C, Müller R. Phase field approximation of dynamic brittle fracture. Comput Mech. 2014;54(5):1141-61. doi:10.1007/s00466-014-1045-x.

9. Kuhn C, Schlüter A, Müller R. On degradation functions in phase field fracture models. Comput Mater Sci. 2015;108:374-84.

10. Amor H, Marigo J-J, Maurini C. Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments. J Mech Phys Solids. 2009;57(8):1209-29. doi:10.1016/j.jmps.2009.04.011.

11. Li T, Marigo J-J, Guilbaud D, Potapov S. Gradient damage modeling of brittle fracture in an explicit dynamics context. Int J Num Methods Eng. 2016. doi:10.1002/nme.5262.

12. Li T, Marigo J-J, Guilbaud D, Potapov S. Variational approach to dynamic brittle fracture via gradient damage models. Appl Mech Mater. 2015;784:334-41.

13. Hofacker M, Miehe C. Continuum phase field modeling of dynamic fracture: variational principles and staggered FE implementation. Int J Fract. 2012;178(1-2):113-29.

14. Borden MJ, Verhoosel CV, Scott MA, Hughes TJR, Landis CM. A phase-field description of dynamic brittle fracture. Comput Methods Appl Mech Eng. 2012;217-220:77-95. doi:10.1016/j.cma.2012.01.008.

15. Dally T, Weinberg K. The phase-field approach as a tool for experimental validations in fracture mechanics. Contin Mech Thermodyn. 2015. doi:10.1007/s00161-015-0443-4.
16. Lorentz E, Cuvilliez S, Kazymyrenko K. Convergence of a gradient damage model toward a cohesive zone model. Comptes Rendus Mécanique. 2011;339(1):20–6. doi:10.1016/j.crme.2010.10.010.

17. Klinsmann M, Rosato D, Kamlah M, McMeeking RM. An assessment of the phase field formulation for crack growth. Comput Methods Appl Mech Eng. 2015;294:313–30.

18. Schlüter A, Kuhn C, Müller R, Gross D. An investigation of intersonic fracture using a phase field model. Archive Appl Mech. 2016;86(1):321–33. doi:10.1007/s00419-015-1114-4.

19. Bourdin B, Francfort GA, Marigo J-J. Numerical experiments in revisited brittle fracture. J Mech Phys Solids. 2000;48(4):797–826.

20. Hossain MZ, Hsueh C-J, Bourdin B, Bhattacharya K. Effective toughness of heterogeneous media. J Mech Phys Solids. 2014;71:15–32. doi:10.1016/j.jmps.2014.06.002.

21. Larsen CJ, Ortner C, Süli E. Existence of solutions to a regularized model of dynamic fracture. Math Models Methods Appl Sci. 2010;20(7):1021–48.

22. Hughes TJR. The finite element method: linear static and dynamic finite element analysis. Englewood Cliffs: Prentice-Hall; 1987.

23. Moré J, Toraldo G. On the solution of large quadratic programming problems with bound constraints. SIAM J Optim. 1991;1(1):93–113.

24. Balay S, Abhyankar S, Adams MF, Brown J, Brune P, Buschelman K, Dalcin L, Eijkhout V, Gropp WD, Kaushik D, Knepley MG, McInnes LC, Rupp K, Smith BF, Zampini S, Zhang H. PETSc users manual. In: Technical report ANL-95/11—Revision 3.6, Argonne National Laboratory; 2015. Accessed 1 May 2016. http://www.mcs.anl.gov/petsc

25. Bourdin B, Larsen CJ, Richardson CL. A time-discrete model for dynamic fracture based on crack regularization. Int J Fract. 2011;168(2):133–43.

26. Logg A, Mardal K-A, Wells G. Automated solution of differential equations by the finite element method: the FEniCS book. Berlin: Springer; 2012. doi:10.1007/978-3-642-23099-8.

27. Sicsic P, Marigo J-J. From gradient damage laws to Griffith's theory of crack propagation. J Elast. 2013;113(1):55–74.

28. Li T, Marigo J-J. Crack tip equation of motion in dynamic gradient damage models. J Elast. 2016. https://hal.archives-ouvertes.fr/hal-01344552/hal-01344552

29. Destuynder P, Djaoua M, Lescure S. Quelques remarques sur la mécanique de la rupture élastique. J de Mécanique Théorique et Appliquée. 1983;2(1):113–35.

30. Attigui M, Petit C. Numerical path independent integral in dynamic fracture mechanics. In: ECF 11—mechanisms and mechanics of damage and failure. 1996.

31. Rice JR. A path independent integral and the approximate analysis of strain concentration by notches and cracks. J Appl Mech. 1968;35(2):379–86.

32. Hakim V, Karma A. Laws of crack motion and phase-field models of fracture. J Mech Phys Solids. 2009;57(2):342–68. doi:10.1016/j.jmps.2008.10.012.

33. Kuhn C, Müller R. A continuum phase field model for fracture. Eng Fract Mech. 2010;77(18):3625–34.

34. Livne A, Ben-David O, Fineberg J. Oscillations in rapid fracture. Phys Rev Lett. 2007;98:124301. doi:10.1103/PhysRevLett.98.124301

35. Dumouchel P-E, Marigo J-J, Charlotte M. Dynamic fracture: an example of convergence towards a discontinuous quasistatic solution. Contin Mech Thermodyn. 2008;20(1):1–19. doi:10.1007/s00161-008-0071-3.

36. Freund LB. Dynamic fracture mechanics. Cambridge: Cambridge University Press; 1990. doi:10.1017/CBO9780511546761

37. Versieux H. A relation between a dynamic fracture model and quasi-static evolution. ESAIM: M2AN. 2016;50(1):77–91.

38. Pham K. Construction et analyse de modèles d'endommagement à gradient. PhD thesis, Université Pierre et Marie Curie (2010)

39. Lazzaroni G, Bargellini R, Dumouchel P-E, Marigo J-J. On the role of kinetic energy during unstable propagation in a heterogeneous peeling test. Int J Fract. 2012;175(2):127–50. doi:10.1007/s10704-012-9708-0

Submit your manuscript to a SpringerOpen journal and benefit from:

► Convenient online submission

► Rigorous peer review

► Immediate publication on acceptance

► Open access: articles freely available online

► High visibility within the field

► Retaining the copyright to your article

Submit your next manuscript at ► springeropen.com