# Enforcement of the irreversibility condition in the phase field damage model: Influence of the history variable approximation

Abdelmalek Barki, Jihed Zghal, Laurent Gallimard, Isabelle Bruant, Nicolas Chevaugeon, Luc Davenne

► To cite this version:

Abdelmalek Barki, Jihed Zghal, Laurent Gallimard, Isabelle Bruant, Nicolas Chevaugeon, et al. Enforcement of the irreversibility condition in the phase field damage model: Influence of the history variable approximation. Engineering Fracture Mechanics, 2025, 328, pp.111518. 10.1016/j.engfracmech.2025.111518. hal-05274120

HAL Id: hal-05274120

https://hal.science/hal-05274120v1

Submitted on 23 Sep 2025

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>

ELSEVIER

Contents lists available at ScienceDirect

Engineering Fracture Mechanics

journal homepage: www.elsevier.com/locate/engfracmech

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

# Enforcement of the irreversibility condition in the phase field damage model: Influence of the history variable approximation

Abdelmalek Barki a, Jihed Zghal a,b,*, Laurent Gallimard a, Isabelle Bruant a, Nicolas Chevaugeon b, Luc Davenne a

## ARTICLE INFO

**Keywords:**

Damage

Phase field

Irreversibility condition

History variable

## ABSTRACT

The phase field damage model has been widely used in recent decades, based on regularizing the local variational problem (Griffith criterion revisited). This original phase field damage model must satisfy the irreversibility condition of cracks, so a constrained variational inequality must be solved. Interest in this method within the mechanical community grew after the publication of Miehe et al., where the damage irreversibility condition is simplified through a history variable in the well-known PF-AT2 model. The damage is governed by this new variable, which prevents a decrease in the elastic potential. This simplification makes it possible to solve an unconstrained variational equality and consequently consider a simplified problem. For the PF-AT2 model, this simplification is consistent. Later, other phase field models (PF-AT1, PF-CZM, ...) were developed, and the use of the history variable was adapted by other researchers to these newer models. In simulating structural damage with the phase field method, many authors have employed the simplified formulation involving the history variable. Therefore, this paper presents a critical comparison of the original (consider damage irreversibility condition) and the simplified (use history variable to approximate the damage irreversibility condition) phase field damage model, through many supported by extensive numerical tests.

## 1. Introduction

Predicting the critical load a structure can withstand, along with the corresponding crack path, is crucial across many fields. Fracture mechanics models crack initiation and propagation in materials subjected to mechanical stress. Its origins lie in the Griffith criterion [1], which introduced two key energy-based assumptions: fracture energy is proportional to the crack surface, and crack propagation follows a minimizing energy principle. However, these assumptions fail to account for important phenomena such as crack nucleation and fatigue, crack growth. Moreover, Griffith's theory is limited to relatively simple cases. These limitations motivated the development of damage mechanics.

The variational approach to brittle fracture introduced by [2] is based on the principle of global energy minimization [1]. At each instant, the energy functional must be minimized over all admissible displacement fields and crack sets. This variational approach presents several limitations, both in theoretical and numerical implementation. Firstly, volume and surface loading cannot be applied: these boundary conditions result in an energy that is not bounded from below, making the problem mechanically and mathematically ill-posed. Only boundary conditions involving imposed displacements can be considered. Even in this case, the

<footnote>a Laboratoire Energetique Mecanique Electromagnetisme (LEME), Paris Nanterre University, 50 rue de Sevres, Ville d'Avray, 92410, France</footnote>

<footnote>b GeM, UMR 6183, 2 Chemin de la Houssiniere, Nantes, 44322, France</footnote>

<footnote>* Corresponding author.
E-mail addresses: barkia@parisnanterre.fr (A. Barki), jzghal@parisnanterre.fr (J. Zghal).</footnote>
Nomenclature

<table><tr><td>R<sup>n</sup></td><td>n-dimensional Euclidean space</td></tr><tr><td>Ω</td><td>Domain of the body</td></tr><tr><td>∂Ω</td><td>Boundary of the domain</td></tr><tr><td>∂Ω<sub>u</sub></td><td>Boundary where displacement is imposed</td></tr><tr><td>∂Ω<sub>t</sub></td><td>Boundary where surface traction is applied</td></tr><tr><td>Γ</td><td>Crack set</td></tr><tr><td>H<sup>n-1</sup>(Γ)</td><td>(n-1)-dimensional Hausdorff measure of the crack set Γ</td></tr><tr><td>X</td><td>Material point in Ω</td></tr><tr><td>t</td><td>Time</td></tr><tr><td>k</td><td>Time step index</td></tr><tr><td>u</td><td>Displacement field</td></tr><tr><td>U</td><td>Imposed displacement</td></tr><tr><td>b*</td><td>Body force</td></tr><tr><td>t*</td><td>Surface traction</td></tr><tr><td>F</td><td>Average reaction force</td></tr><tr><td>F<sub>c</sub></td><td>Average reaction force at failure state</td></tr><tr><td>U<sub>0</sub></td><td>Imposed displacement at which damage threshold is first exceeded</td></tr><tr><td>U<sub>c</sub></td><td>Imposed displacement required to reach failure</td></tr><tr><td>ε</td><td>Strain tensor</td></tr><tr><td>ε<sub>dev</sub></td><td>Deviatoric part of the strain tensor</td></tr><tr><td>σ</td><td>Cauchy stress tensor</td></tr><tr><td>σ<sub>0</sub></td><td>Stress tensor of the undamaged material</td></tr><tr><td>∇</td><td>Gradient operator</td></tr><tr><td>Δ</td><td>Laplacian operator</td></tr><tr><td>d</td><td>Damage field</td></tr><tr><td>d<sub>c</sub></td><td>Critical damage field</td></tr><tr><td>D<sub>c</sub></td><td>Half-width of the damage localization band</td></tr><tr><td>ℓ</td><td>Regularization length</td></tr><tr><td>γ(d, ∇d)</td><td>Crack surface density function</td></tr><tr><td>α(d)</td><td>Geometric crack function</td></tr><tr><td>g(d)</td><td>Degradation or softening function</td></tr><tr><td>H</td><td>History variable</td></tr><tr><td>H<sub>0</sub></td><td>Damage threshold</td></tr><tr><td>Y&#773;</td><td>Damage driving potential</td></tr><tr><td>ψ<sub>0</sub></td><td>Elastic potential of the undamaged material</td></tr><tr><td>ψ</td><td>Degraded elastic potential</td></tr><tr><td>ψ<sub>0</sub><sup>+</sup>, ψ<sub>0</sub><sup>-</sup></td><td>Tensile and compressive parts of ψ<sub>0</sub></td></tr><tr><td>ε</td><td>Total energy</td></tr><tr><td>ε<sub>e</sub></td><td>Elastic energy</td></tr><tr><td>ε<sub>d</sub></td><td>Fracture energy</td></tr><tr><td>ε<sub>ext</sub></td><td>External energy</td></tr><tr><td>dε</td><td>Weak form of the damage problem</td></tr><tr><td>ε<sub>d</sub><sup>c</sup></td><td>Critical fracture energy</td></tr><tr><td>C</td><td>Set of admissible displacement fields</td></tr><tr><td>D</td><td>Set of admissible damage fields</td></tr><tr><td>G<sub>c</sub></td><td>Critical energy release rate (fracture toughness)</td></tr><tr><td>c</td><td>Normalization constant for crack surface energy</td></tr><tr><td>Q(d)</td><td>Denominator term in the softening function g(d)</td></tr><tr><td>a<sub>1</sub>, a<sub>2</sub>, a<sub>3</sub></td><td>Parameters of softening type</td></tr><tr><td>ξ</td><td>Weighting parameter in α(d)</td></tr><tr><td>p</td><td>Exponent in the softening function g(d)</td></tr></table>
<table><tr><td>σ<sub>c</sub></td><td>Critical stress or failure strength</td></tr><tr><td>K<sub>n</sub></td><td>Bulk modulus</td></tr><tr><td>λ</td><td>First Lamé parameter</td></tr><tr><td>μ</td><td>Shear modulus (second Lamé parameter)</td></tr><tr><td>E</td><td>Young's modulus</td></tr><tr><td>ν</td><td>Poisson's ratio</td></tr><tr><td>Δ<sub>r</sub></td><td>Relative error</td></tr><tr><td>h</td><td>Finite element mesh size</td></tr><tr><td>CPU</td><td>Computational time</td></tr><tr><td>|·|</td><td>Absolute value</td></tr><tr><td>⟨a⟩<sup>±</sup></td><td>Positive/negative part of scalar a: ⟨a⟩<sup>±</sup> = <sup>1</sup>&frasl;<sub>2</sub>(a ± |a|)</td></tr><tr><td>PF-AT2</td><td>AT2 phase field model</td></tr><tr><td>PF-AT1</td><td>AT1 phase field model</td></tr><tr><td>PF-CZM</td><td>CZM phase field model</td></tr></table>

existence of solutions can only be ensured under restrictive assumptions. The main limitations are the non-convexity of the energy functional and the discontinuities in displacement fields, caused by one or multiple cracks. The variational approach to brittle fracture, with its limitations, has enabled the emergence of other approaches. Among these approaches is the regularization of the variational approach to brittle fracture in the sense of Γ-convergence, proposed by [3,4]. Their results demonstrate that the global minimum of the regularized functional converges to that of the original discontinuous formulation. This regularization introduces diffuse crack regions controlled by a scalar, the regularization length, thus smoothing the sharp crack surfaces in the solid. This framework has been refined and significantly extended by [5,6], evolving into the phase field damage model. This approach needed adaptations including the use of vector fields instead of scalar ones, taking into account irreversibility, boundary conditions, and more general loading. Due to the non-convexity of the regularized energy functional, resulting from the major difference between damage and displacement, finding a global minimum remains difficult, if not impossible. [5] introduced a new solution strategy in damage mechanics: alternate minimization. As the regularized energy functional is separately convex, this new solution strategy simplifies numerical implementation and tractability. The obtained minima are local minima of the energy functional rather than global minima. The Γ-convergence results do not necessarily apply to these local minima. Thus, this phase field damage model no longer approximates the variational approach to brittle fracture. Despite this, the numerical results from [5] are shown to be consistent with experimental observations. Moreover, [2,7] have explicitly noted that the assumption of global energy minimization does not always find justification from a physical perspective. This solution strategy allows considering boundary conditions of volume and surface loading. The development of phase field models involves choosing a damage parameter, its dependency on the stress-strain relationship, and an evolution law [8,9]. The stress equilibrium equation (quasi-static or dynamic) and a gradient-type evolution equation for the damage field that accounts for the irreversibility condition describe the cracking process.

In the decade following the initial works on the phase field damage model, interest in the mechanics community has increased significantly. In particular, the works of [10] have played a determining role. [10] simplify the damage irreversibility condition by introducing a history variable in the case of the well-known PF-AT2 model. This simplification is consistent because the PF-AT2 model does not have a threshold: damage is triggered by infinitesimal loading. The damage evolution is then governed by this new variable, which prevents a potential decrease in the elastic potential. This simplification allows for solving an unconstrained variational equality and consequently considering a simplified problem. Later, other phase field models (like PF-AT1, PF-CZM, etc.) have been developed, and the use of a history variable has been adapted for them by other researchers [11–13].

The ability of the phase field damage model to predict consistent damage nucleation has been discussed in recent publications. The first observation was made early by [6], where the material strength is not considered in the initial version of the phase field damage model (called PF-AT2). Then, the first cohesive zone phase field models appeared [14] and were further developed by [11]. In the cohesive zone model, the failure strength is taken into account in an attempt to represent the material strength. This addition improves the phase field damage model, but it remains insufficient. To overcome this limitation, various phase field models and elastic potential decompositions have been proposed. The works of [15,16] conclude that no phase field model or decomposition added to the initial version of the phase field damage model can capture the full range of phenomena from nucleation to propagation. To address this, [16] proposed revisiting the phase field damage model by introducing a driving force in the damage problem. This force should allow for consistent damage nucleation in the presence of large, small, or no pre-cracks [15–17].

In this paper, the focus is placed on the use of the history variable in the classical phase field damage model. A critical analysis between the original phase field damage model and the simplified phase field damage model (using the history variable) is provided through several numerical tests. The outline of this work is as follows. In Section 2, the theory behind the phase field damage model is presented, along with its implementation. Then, in Section 3, the influence of the history variable on the predictions of the phase field damage model is illustrated through several numerical tests. First, the differences between the original and simplified problems are recalled for comparison purposes. The numerical tests conducted are detailed in Section 3.1 to Section 3.5. Finally, the main conclusions are summarized in the last section.
## 2. Phase field damage model: review of the theoretical background

### 2.1. Approach variational to the brittle fracture

The phase field method relies on a revisitation of the Griffith criterion into a local variational formulation, as proposed by [2]. Consider a domain $\Omega \subset \mathbb{R}^n$ (with $n = 1, 2, 3$) containing a crack $\Gamma$ (Fig. 1(a)). The solid is characterized by an infinitesimal displacement field $u$, and the associated strain tensor is given by:

$$ \epsilon = \frac{1}{2} (\nabla u + \nabla u^T) \qquad (1) $$

The displacement is imposed on a part of the boundary denoted $\partial\Omega_u$, which is assumed to be fixed over time. The imposed displacement is denoted by $U$. The mechanical behavior of the material is described through the total energy $\mathcal{E}$ (Eq. (2)), which is expressed as the sum of the elastic energy $\mathcal{E}_e$ and the fracture energy $\mathcal{E}_d$:

$$ \mathcal{E}(u, \Gamma) = \underbrace{\int_{\Omega \setminus \Gamma} \psi_0(\epsilon(u)) dV}_{\mathcal{E}_e} + \underbrace{G_c \mathcal{H}^{n-1}(\Gamma)}_{\mathcal{E}_d} \qquad (2) $$

where:

* $\psi_0(\epsilon)$ is the elastic potential of sound material,

* $G_c$ is the fracture toughness of the material,

* $\mathcal{H}^{n-1}(\Gamma)$ is the $(n-1)$ dimensional Hausdorff measure of the cracks set $\Gamma$.

In accordance with the Griffith criterion, under quasi-static loading conditions, the total energy $\mathcal{E}$ must be minimized over all admissible displacement fields and crack sets, subject to an irreversibility condition on crack evolution. To model the quasi-static evolution of the system, an implicit time discretization is introduced. We consider a discrete sequence of time steps $\{t^k\}_{k=0}^N$, defined by $t^{k+1} = t^k + \Delta t$, with a time increment $\Delta t$. The displacement and crack set at time step $k$ are then denoted as $u^k := u(\mathbf{X}, t^k)$ and $\Gamma^k := \Gamma(t^k)$, with $(\mathbf{X}, t^k) \in \Omega \times [0, t^N]$. The minimization problem at time step $k+1$ can be written as:

$$ (u^{k+1}, \Gamma^{k+1}) = \arg\min_{u \in \mathcal{C}^{k+1}, \Gamma^k \subset \Gamma} \mathcal{E}(u, \Gamma) \qquad (3) $$

where the set of admissible displacements is defined by:

$$ C^{k+1} = \{u \in H^1(\Omega \setminus \Gamma) : u = U^{k+1} \text{ on } \partial\Omega_u \setminus \Gamma\} \qquad (4) $$

The implementation of the minimization of Eq. (2), as formulated in Eq. (3), is not feasible. To overcome these limitations, a regularization by a damage phase field and a strategy of alternate minimization are proposed.

### 2.2. Phase field damage model

To simplify the minimization of the total energy $\mathcal{E}$ (Eq. (2)), the crack set is represented by the damage field (Fig. 1(b)). The main idea is to approximate the fracture energy over the volume using a crack surface density function. The unified formulation proposed by [11] is used, which leads to Eq. (5):

$$ \mathcal{E}_d = \int_{\Gamma} G_c dS \approx \int_{\Omega} G_c \gamma(d, \nabla d) dV \qquad (5) $$

with

$$ \gamma(d, \nabla d) = \frac{1}{c\ell} (\alpha(d) + \ell^2 |\nabla d|^2) $$

where

* $d$ is the damage variable. $d \in [0, 1]$ with $d = 0$ for sound material and $d = 1$ for fully damaged material

* $\alpha(d)$ is the geometric crack function

* $c = 4 \int_0^1 \frac{1}{\sqrt{\alpha(\phi)}} d\phi$ is a scaling parameter to recover the sharp crack surface

* $\ell$ is the regularization length.

The total energy Eq. (2) becomes:

$$ \mathcal{E} = \underbrace{\int_{\Omega} g(d) \psi_0(\epsilon) dV}_{\mathcal{E}_e} + \underbrace{\int_{\Omega} \frac{G_c}{c\ell} (\alpha(d) + \ell^2 |\nabla d|^2) dV}_{\mathcal{E}_d} - \underbrace{\int_{\Omega} b^* u dV - \int_{\partial\Omega_t} t^* u dS}_{\mathcal{E}_{ext}} \qquad (6) $$

with
<figure><img src="image_4.png" /><figcaption>Fig. 1. Solid $\Omega$ with crack set $\Gamma$ and damage $d$.</figcaption></figure>

$g(d)$ is the softening or degradation function,

$b^*$ is the distributed body force and $t^*$ surface traction over $\partial\Omega_t$.

According to the literature, the phase field damage problem can be resolved in two main ways. The first approach is based on minimizing the total energy $\mathcal{E}$ under box constraints to ensure damage irreversibility. This is a classical method, in which the minimum damage bound is updated at each increment. The second approach was proposed by [10]. A history variable is used to approximate the damage irreversibility condition and to simplify the damage problem. This variable is defined as:

$$H(u, t) = \max_{0 \le \tau \le t} (\psi_0(\epsilon(u, \tau)), H_0) \quad (7)$$

where $H_0$ is the damage threshold, defined as a function of the chosen phase field model. In the following, we will refer to the model where the history variable is used as the **simplified problem**, and to the one where the history variable is not used as the **original problem**.

The alternate minimization proposed by [5] is used. It consists of separately solving for $u$ and $d$, which minimize the total energy $\mathcal{E}$.

### Displacement problem:

$$u^{k+1} = \arg \min_{u \in C^{k+1}} \mathcal{E}(u, d^k) \quad (8)$$

where the set of admissible displacements is defined by:

$$C^{k+1} = \{u \in H^1(\Omega) : u = U^{k+1} \text{ on } \partial\Omega_u\} \quad (9)$$

### Damage problem:

#### Original problem

$$d^{k+1} = \arg \min_{d \in D^{k+1}} \mathcal{E}(u^{k+1}, d) \quad (10)$$

where the set of admissible damage fields is defined by:

$$D^{k+1} = \{d \in H^1(\Omega) : d \ge d^k \text{ on } \Omega\} \quad (11)$$

#### Simplified problem

$$H^{k+1} = \max_{d \in H^1(\Omega)} (\psi_0^{k+1}, H^k, H_0) \\ d^{k+1} = \arg \min_{d \in H^1(\Omega)} \mathcal{E}(u^{k+1}, d; H^{k+1}) \quad (12)$$

The two solution approaches to the phase field problem can be summarized through the following strong form:

### Displacement problem:

$$\operatorname{div}(\sigma) + b^* = 0 \quad \text{in } \Omega \quad (13)$$

$$u = U \quad \text{on } \partial\Omega_u \quad (14)$$

$$\sigma \cdot n = t^* \quad \text{on } \partial\Omega_t \quad (15)$$
<table><caption>**Table 1**
Parameters of phase field models [11].</caption><thead><tr><th>Model</th><th>Softening type</th><th>c</th><th>α</th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>3</sub></th><th>p</th><th>H<sub>0</sub></th></tr></thead><tbody><tr><td>PF-AT2</td><td></td><td>2</td><td>d<sup>2</sup></td><td>2</td><td>-0.5</td><td>0</td><td>2</td><td>0</td></tr><tr><td>PF-AT1</td><td></td><td>8/3</td><td>d</td><td>2</td><td>-0.5</td><td>0</td><td>2</td><td>3G<sub>c</sub>/16ℓ</td></tr><tr><td rowspan="2">PF-CZM</td><td>Linear</td><td>π</td><td>2d - d<sup>2</sup></td><td>-</td><td>-0.5</td><td>0</td><td>2</td><td>σ<sub>c</sub><sup>2</sup>/2E</td></tr><tr><td>Cornelissen</td><td>π</td><td>2d - d<sup>2</sup></td><td>-</td><td>1.3868</td><td>0.6567</td><td>2</td><td>σ<sub>c</sub><sup>2</sup>/2E</td></tr></tbody></table>

## Damage problem:

### Original problem

$$ \dot{d} \geq 0 \quad \text{in } \Omega \qquad (16) $$

$$ g'(d)\psi_0(\epsilon) + \frac{G_c}{c\ell} (\alpha'(d) - 2\ell^2\Delta d) \geq 0 \quad \text{in } \Omega \qquad (17) $$

$$ \dot{d} \left[ g'(d) \psi_0(\epsilon) + \frac{G_c}{c\ell} (\alpha'(d) - 2\ell^2 \Delta d) \right] = 0 \quad \text{in } \Omega \qquad (18) $$

$$ \nabla d \cdot n \geq 0 \quad \text{on } \partial\Omega \qquad (19) $$

$$ \dot{d} \nabla d \cdot n = 0 \quad \text{on } \partial\Omega \qquad (20) $$

### Simplified problem

$$ H(u, t) = \max_{0 \le \tau \le t} (\psi_0(\epsilon(u, \tau)), H_0) \qquad (21) $$

$$ g'(d) H(\epsilon, t) + \frac{G_c}{c\ell} (\alpha'(d) - 2\ell^2 \Delta d) = 0 \quad \text{in } \Omega \qquad (22) $$

$$ \nabla d \cdot n = 0 \quad \text{on } \partial\Omega \qquad (23) $$

## 2.3. Implementation

Different phase field models were implemented using the unified formulation [11], by defining for each model:

### The degradation function:

$$ g(d) = \frac{(1-d)^p}{(1-d)^p + Q(d)} \qquad (24) $$

$$ Q(d) = a_1 d(1 + a_2 d + a_3 a_2 d^2), \quad Q(d) > 0, \quad \text{and } p \ge 2 $$

with $a_1 = \frac{4G_c E}{\pi \ell \sigma_c^2}$, and $a_2, a_3$ should be determined to fit the chosen cohesive law. $\sigma_c$ is the failure strength.

### The geometric crack function:

$$ \alpha(d) = \xi d + (1 - \xi)d^2 \qquad (25) $$

where Eq. (25) should verify the following conditions:

$$ \alpha(0) = 0, \quad \alpha(1) = 1, \quad \text{and } \alpha(d) \ge 0 \Rightarrow \xi \in [0, 2] $$

To switch from one model to another (PF-AT1, PF-AT2, PF-CZM...), we just have to choose the appropriate variables as described in Table 1.

For the PF-CZM model, $a_1$ should be calculated as shown below, and $a_2$ and $a_3$ should be identified from the traction-separation law [11].

The elastic potential decomposition aims to respect the unilateral contact condition. Damage initiates when the elastic potential exceeds a damage threshold. To ensure the unilateral contact condition, only the tensile stresses should be considered in the damage process. In the initial formulation without elastic potential decomposition (referred to here as the symmetric model), non-physical crack propagation can appear. [5] attributes this to the lack of consideration for material impenetrability. Following this rationale, the volumetric/deviatoric decomposition from [18] was developed. Simultaneously, the spectral decomposition from [10] was introduced to account only for the tensile stresses in the damage problem. Several decompositions exist, each with its own advantages and weaknesses. To our knowledge, no decomposition is optimal in a general context. The volumetric/deviatoric decomposition is used in this study.
<figure><img src="image_5.png" /><figcaption>Fig. 2. Implementation Strategy in FEniCS 2019 with PETSc Library Solvers [20].</figcaption></figure>

### Symmetric model [6]:

$$ \psi(u, d) = \frac{1}{2}\epsilon : \sigma = g(d)\frac{1}{2}\epsilon : \sigma_0 = g(d)\left[\frac{1}{2}\lambda\mathrm{tr}(\epsilon)^2 + \mu\epsilon : \epsilon\right] \quad (26) $$

$$ \sigma(u, d) = \frac{\partial\psi}{\partial\epsilon}(u, d) = g(d)\sigma_0 = g(d)[\lambda\mathrm{tr}(\epsilon)I + 2\mu\epsilon] \quad (27) $$

### Volumetric/deviatoric decomposition [18]: It consists in decompose $\psi_0$ to $\psi_0(\epsilon) = \psi_0(\epsilon)^+ + \psi_0(\epsilon)^-$ and to consider the deviatoric part of $\epsilon$:

$$ \epsilon^{\mathrm{dev}} = \epsilon - \frac{1}{3}\mathrm{tr}(\epsilon)I $$

The exponent + indicates the positive part and – the negative part. [18] have proposed:

$$ \begin{cases} \psi_0(\epsilon)^+ = \frac{1}{2}K_n\langle\mathrm{tr}(\epsilon)\rangle_+^2 + \mu(\epsilon^{\mathrm{dev}} : \epsilon^{\mathrm{dev}}) \\ \psi_0(\epsilon)^- = \frac{1}{2}K_n\langle\mathrm{tr}(\epsilon)\rangle_-^2 \end{cases} \Rightarrow \psi(\epsilon, d) = g(d)\psi_0(\epsilon)^+ + \psi_0(\epsilon)^- \quad (28) $$

with, $K_n = \lambda + \frac{2\mu}{n}$ and $\langle a \rangle_\pm = \frac{1}{2}(a \pm |a|)$. $K_n$ represent the bulk modulus and n the dimension of problem.

Thus the stress tensor $\sigma$ is also decomposed:

$$ \begin{cases} \sigma_0^+ = K_n \langle \mathrm{tr}(\epsilon) \rangle_+ I + 2\mu \epsilon^{\mathrm{dev}} \\ \sigma_0^- = K_n \langle \mathrm{tr}(\epsilon) \rangle_- I \end{cases} \Rightarrow \sigma = g(d)\sigma_0^+ + \sigma_0^- \quad (29) $$

Using this decomposition, the history variable is written as:

$$ H(u, t) = \max_{0 \le \tau \le t} (\psi_0(\epsilon(u, \tau)))^+, H_0 \quad (30) $$

The unified phase field damage model was implemented in FEniCS [19], following the strategy described in Fig. 2.

In this paper, the influence of the history variable on the damage prediction of different phase field models, such as PF-AT1, PF-AT2, and PF-CZM, is investigated.

### 3. Influence of history variable approximation on phase field damage model prediction

The main difference between the original and simplified problems lies in Eqs. (17) and (22). In the original problem, the elastic potential is represented by its classical expression, which requires solving a variational inequality constrained by the damage irreversibility condition. In the simplified problem, a numerical artefact, the history variable H, is introduced to approximate the damage irreversibility condition, thereby simplifying the problem into an unconstrained variational equality.

Finally, to observe the influence of the history variable H approximation on the phase field damage model, one should focus only on the damage problem, particularly on the part involving the elastic potential in both formulations. The weak form of the damage problem is defined as the sum of an elastic part, $d\mathcal{E}_e$, and a fracture part, $d\mathcal{E}_d$:

$$ d\mathcal{E} = \underbrace{\int_{\Omega} g'(d) \bar{Y} d^* dV}_{d\mathcal{E}_e} + \underbrace{\int_{\Omega} \frac{G_c}{c^\ell} (\alpha'(d) d^* + 2c^\ell^2 \nabla d \nabla d^*) dV}_{d\mathcal{E}_d} \quad (31) $$
<figure><img src="image_6.png" /><figcaption>Fig. 3. Uni-axial bar — geometry definition (dimensions are in millimeters).</figcaption></figure>

<table><caption>**Table 2**

Uni-axial bar — material properties.</caption><tr><td>Young's Modulus : E [MPa]</td><td>10,000</td></tr><tr><td>Poisson's ratio : v</td><td>0.3</td></tr><tr><td>Toughness : G<sub>c</sub> [kN/m]</td><td>1</td></tr><tr><td>Failure strength : σ<sub>c</sub> [MPa]</td><td>324.8</td></tr></table>

with,

$$ \bar{Y} = \begin{cases} \psi_0(\epsilon) & \text{for the original problem} \\ H(\epsilon, t) & \text{for the simplified problem} \end{cases} \quad \text{and } d^* \text{ is the test variable} $$

where the sets of admissible damage test is defined by :

$$ D^* = \begin{cases} \{d^* \in H^1(\Omega) : d^* \ge 0 \text{ on } \Omega\} & \text{for the original problem} \\ \{d^* \in H^1(\Omega)\} & \text{for the simplified problem} \end{cases} $$

In the original problem, $d\mathcal{E}_e$ and $d\mathcal{E}_d$ respectively represent the variations of the elastic energy $\mathcal{E}_e$ and the fracture energy $\mathcal{E}_d$ with respect to the damage variable $d$. Consequently, $d\mathcal{E}$ represents the variation of the total energy $\mathcal{E}$ with respect to $d$. In contrast, in the simplified problem, $d\mathcal{E}_e$ does not represent this variation due to the presence of the history variable $H$.

To compare the numerical results obtained using these two problems, the meshes are the same for each of the following simulations. To illustrate the influence of the history variable, two key aspects are studied: the temporal evolution of the damage and the spatial evolution of the damage. The fracture energy depends mainly on these aspects. Referring to Table 1, the damage threshold $H_0$ is clearly defined. In the PF-AT2 model, damage initiates from the very first loading increment, even when the applied load is infinitesimally small. As a result, the damage threshold $H_0$ is equal to zero, and in this study, the PF-AT2 model is considered part of the class of phase field models without a damage threshold. The PF-AT1 and PF-CZM models belong to the class of phase field models with a damage threshold. For each simulation, the convergence tolerance of the staggered algorithm is $1 \times 10^{-5}$. To simplify, the regularization length for each test is determined using the analytical expression computed for a uni-axial stress state in Appendix B.1.

In the following numerical tests, $U_0$ is the imposed displacement required to reach the damage threshold $H_0$ at least at one node. The critical values $U_c$ and $F_c$ represent the imposed displacement and the average reaction force required to reach the failure state, respectively.

## 3.1. Bar structures

In this section, the phase field damage predictions on two bars are analyzed. The first one is referred to as the uni-axial bar, while the second is a V-notch specimen.

### 3.1.1. Uni-axial bar

A bar subjected to uni-axial tension is considered, as illustrated in Fig. 3. The plane stress hypothesis is applied. The material properties are given in Table 2. The boundary conditions are described in Fig. 3. A displacement is applied in the x-direction, while the y-direction is left free. The increment of imposed displacement is $\Delta U = 1 \times 10^{-5}$ mm. The PF-AT2, PF-AT1, and PF-CZM (linear) models are used. Table 3 presents the different values of $\ell$ and the corresponding analytical critical values for a one-dimensional bar. The structure is discretized using a uniform mesh of 58 136 finite elements, resulting in a mesh size of $h = \ell_{AT2}/5 = 2 \times 10^{-3}$ mm.

### Results for the phase field models with threshold damage (PF-AT1 and PF-CZM)

For PF-AT1 and PF-CZM, the two approaches are compared. In Fig. 4, the force-displacement responses are the same for both approaches. Fig. 5 presents the evolution of fracture energy with respect to the loading. Before failure, the fracture energy is identical. However, at the failure point, the simplified approach overestimates the critical fracture energy. In Table 4, the relative errors $\Delta\ell(X)$
<table><caption>**Table 3**

Uni-axial bar — critical analytical values bar from Appendix B.1.</caption><thead><tr><th></th><th>ℓ [mm]</th><th>U<sub>c</sub> [mm]</th><th>F<sub>c</sub> [N]</th></tr></thead><tbody><tr><td>PF-AT1</td><td>0.036</td><td>0.032</td><td>324.8</td></tr><tr><td>PF-AT2</td><td>0.01</td><td>0.058</td><td>324.8</td></tr><tr><td>PF-CZM</td><td>0.036</td><td>0.0325</td><td>324.8</td></tr></tbody></table>

<table><caption>**Table 4**

Phase field model with threshold damage — relative errors Δ<sup>a</sup><sub>r</sub>, [%], compared with the analytical for the uni-axial bar.</caption><thead><tr><th>Model</th><th>Problem</th><th>Δ<sup>a</sup><sub>r</sub>(F<sub>c</sub>)</th><th>Δ<sup>a</sup><sub>r</sub>(U<sub>c</sub>)</th><th>Δ<sup>a</sup><sub>r</sub>(ε<sup>c</sup><sub>d</sub>)</th></tr></thead><tbody><tr><td rowspan="2">PF-AT1</td><td>Original</td><td>0.02</td><td>0.01</td><td>0.13</td></tr><tr><td>Simplified</td><td>0.02</td><td>0.01</td><td>31.9</td></tr><tr><td rowspan="2">PF-CZM</td><td>Original</td><td>0.02</td><td>0.01</td><td>0.05</td></tr><tr><td>Simplified</td><td>0.02</td><td>0.01</td><td>20.2</td></tr></tbody></table>

<table><caption>**Table 5**

Phase field model without threshold damage — relative errors Δ<sup>a</sup><sub>r</sub>, [%], compared with the analytical for the uni-axial bar.</caption><thead><tr><th>Model</th><th>Problem</th><th>Δ<sup>a</sup><sub>r</sub>(F<sub>c</sub>)</th><th>Δ<sup>a</sup><sub>r</sub>(U<sub>c</sub>)</th></tr></thead><tbody><tr><td rowspan="2">PF-AT2</td><td>Original</td><td>2.22</td><td>3.23</td></tr><tr><td>Simplified</td><td>2.22</td><td>3.23</td></tr></tbody></table>

for each problem are presented and compared with the analytical predictions from Appendix B.1 and Appendix C. The relative error is defined as:

$$ \Delta_r^a(X) = \left| \frac{X - X^a}{X^a} \right| \qquad (32) $$

where X<sup>a</sup> denotes the analytical value, and X corresponds to the result obtained from either the original or the simplified problem. The errors are the same for the critical force F<sub>c</sub> and displacement U<sub>c</sub>. However, the error in the critical fracture energy ε<sup>c</sup><sub>d</sub> is significantly larger for the simplified problem. The analytical (Appendix B.2) and numerical critical damage fields are shown in Fig. 6.

The original problem reproduces the same critical damage profile as the analytical prediction Appendix B.2. The critical damage field is much more diffuse in the simplified problem than in the original one. This increased diffusion of the critical damage field explains why the fracture energy is overestimated.

### Results for the phase field model without threshold damage (PF-AT2)

For the PF-AT2 model, the simplified and original problems yield similar force-displacement responses. In both cases, the fracture energy is overestimated because this phase-field model presents an infinite damage localization band, as indicated in Appendix B.2. The absence of a damage threshold explains this diffuse damage profile. The boundary condition enforcing zero damage at the edges (left and right) prevents damage from localizing at the boundaries. The shape of the numerical critical localization band can be explained by the lack of a threshold. Therefore, without a sharp angle causing a stress concentration, the comparison of fracture energy with the analytical prediction Appendix C is not straightforward. The weaknesses of the PF-AT2 model are discussed in [21]. The relative error in the critical fracture energy between the two problems is 1.53%. This difference is minimal because the difference in damage diffusion is negligible, as illustrated in Fig. 8. The relative errors in critical force and displacement are identical. In summary, the influence of the history variable H on the results is not significant. In this case, the simplified problem reproduces the same responses as the original one. (see Fig. 7 and see Table 5).

### 3.1.2. Bar with a V-notch

Consider a bar with a V-notch subjected to tensile loading, as illustrated in Fig. 9. The plane stress hypothesis is used. The boundary conditions are applied as shown in Fig. 9. A displacement is applied in the x-direction, while the y-direction is left free. The loading increments are ΔU = 1 × 10<sup>-5</sup> mm. The material properties are detailed in Table 6. The PF-AT2, PF-AT1, and PF-CZM (linear) models are used. The different values of ℓ are given in Table 7. The structure is discretized with a non-uniform mesh, refined in the notch zone, consisting of 101 605 finite elements, as illustrated in Fig. 10. The mesh size in the crack zone is equal to h = ℓ/AT2 = 2.8 × 10<sup>-4</sup> mm.
<figure><img src="image_7.png" /><figcaption>Fig. 4. F vs U for uni-axial bar.</figcaption></figure>

**Table 6**

Bar with a V-notch — material properties.

<table><tr><td>Young's Modulus : E [MPa]</td><td>10,000</td></tr><tr><td>Poisson's ratio : ν</td><td>0.3</td></tr><tr><td>Toughness : G<sub>c</sub> [kN/m]</td><td>1</td></tr><tr><td>Failure strength : σ<sub>c</sub> [MPa]</td><td>868.0</td></tr></table>

**Table 7**

Bar with a V-notch — regularization lengths.

<table><thead><tr><th></th><th>PF-AT1</th><th>PF-AT2</th><th>PF-CZM</th></tr></thead><tbody><tr><td>ℓ [mm]</td><td>0.005</td><td>0.0014</td><td>0.005</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** $\varepsilon_d$ vs $U$ for the uni-axial bar.</figcaption></figure>

Results for the phase field models with damage threshold (PF-AT1 and PF-CZM)

In Fig. 11(a), the force-displacement responses differ between the two approaches. The simplified problem reaches failure under a lower load and overestimates the critical fracture energy, as shown in Fig. 11(b). This overestimation is due to a more diffuse damage field in the simplified problem, as illustrated in Fig. 12.

Figs. 13 and 14 presents the quantities $d\varepsilon_e$, $d$, and $\bar{Y}$ in the damage zone for the PF-AT1 model. As a reminder, $\bar{Y}$ represents the elastic potential in the original problem and the history variable in the simplified problem. At the notch tip, the generalized driving force $d\varepsilon_e$ is the same for both the original and simplified problems (Fig. 13(a)). However, at the same location, the damage is more diffuse and of greater magnitude in the simplified problem, as shown in Fig. 13(b). This difference in damage leads, under the same loading, to a variation in the $\bar{Y}$ values between the two problems. Fig. 14 shows the evolution of $\bar{Y}$ at the notch tip before failure. After a certain loading point, $\bar{Y}$ increases more rapidly in the simplified problem. Unlike the case of the bar without a notch, the force-displacement responses here differ, because the damage intensity differs. The failure state is reached earlier in the simplified problem. This difference in critical fracture energy between the two problems is explained by the more diffuse damage in the simplified case. The relative errors, reported in Table 9, are slightly larger for the PF-AT1 model than for the PF-CZM model. There is no clear correlation between having more diffuse damage and reaching failure more quickly (see the uni-axial bar test for the PF-AT2 model, Fig. 3).
<figure><img src="image_9.png" /><figcaption>Fig. 6. d vs x for the uni-axial bar.</figcaption></figure>

A detailed comparison between computational cost and accuracy is presented in Table 8 for the PF-AT1 model applied to the notched bar problem. Relative errors in key quantities ($F_c$, $U_c$, and $\varepsilon_d^c$) are reported alongside total CPU times for each problem.

The relative errors are calculated with respect to the original problem, which is considered the reference. The relative error is defined as:

$$ \Delta_r(X) = \left| \frac{X^{\text{Simplified}} - X^{\text{Original}}}{X^{\text{Original}}} \right| \qquad (33) $$

where $X^{\text{Original}}$ denotes the value from the original problem, and $X^{\text{Simplified}}$ corresponds to the result obtained from the simplified problem. This definition is used for all subsequent simulations.

It is important to emphasize that the original problem can be solved using two different numerical strategies, as shown in Fig. 2. The first, Newton semi-smooth line search (ssls), solves the system over all degrees of freedom. The second, Newton reduced-space active-set line search (rsls), exploits the structure of the inequality constraint by identifying active and inactive sets, and solves the problem only on the active degrees of freedom. Both approaches solve the same mathematical problem, only the numerical methods differ, but the rsls method achieves a significant reduction in computational time without altering the solution.
<figure><img src="image_10.png" /><figcaption>Fig. 7. Results of PF-AT2 model for the uni-axial bar.</figcaption></figure>

In contrast, the simplified problem is based on an approximate formulation that is inherently less costly to solve. This reduction in computational time comes at the expense of accuracy, with moderate errors in force and displacement (around 7%) and a more significant discrepancy in the fracture energy (25.2%).

In conclusion, the simplified problem shows clear computational benefits when the original problem is solved using full-space methods like *ssls*. However, when an efficient reduced-space method such as *rsls* is employed, the original problem can be solved in a time comparable to, or even faster than, the simplified one, without any compromise on accuracy. In such cases, the simplified model becomes less relevant, as similar performance can be achieved with the full model using a more efficient numerical solver.

In addition to the previous results, an analysis of the computational time as a function of mesh refinement is included to highlight the potential advantages of the simplified problem.

The simplified problem significantly reduces computation time compared to the original problem solved with the *ssls* method, as illustrated in Fig. 15. However, solving the original problem using the Newton reduced-space active-set line search method (*rsls*) also leads to a considerable reduction in computational time.
<figure><img src="image_11.png" /><figcaption>Fig. 8. PF-AT2 model — *d* vs *x* for the uni-axial bar.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. Bar with a V-notch — geometry definition (dimensions are in millimeters).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Bar with a V-notch — mesh.</figcaption></figure>

<table><caption>**Table 8**
PF-AT1 model — Accuracy vs. computational cost for the bar with a V-notch.</caption><thead><tr><th>Model</th><th>Δ<sub>r</sub>(F<sub>c</sub>)</th><th>Δ<sub>r</sub>(U<sub>c</sub>)</th><th>Δ<sub>r</sub>(ε<sub>d</sub><sup>0</sup>)</th></tr></thead><tbody><tr><td>PF-AT1</td><td>7.17</td><td>6.91</td><td>25.2</td></tr></tbody></table>

Relative errors Δ<sub>r</sub>, [%], between the original (reference) and simplified problems.

<table><thead><tr><th></th><th>CPU time</th></tr></thead><tbody><tr><td>Original - sssl</td><td>1 d 17 m 59 s</td></tr><tr><td>Original - rssl</td><td>2 h 58 m 42 s</td></tr><tr><td>Simplified</td><td>3 h 26 m 46 s</td></tr><tr><td>Total CPU times.</td><td></td></tr></tbody></table>

As the mesh is refined, the gap in computational cost between the simplified and original problems (solved with sssl) becomes more pronounced. Nevertheless, the use of rssl narrows this gap substantially. As a result, the computational times of the original and simplified problems become comparable, and for certain mesh sizes, even slightly lower for the original problem when rssl is employed.

This observation is particularly important: both the sssl and rssl methods solve the same original problem, but the rssl method significantly reduces the computational cost. Consequently, the simplified problem loses its main advantage, faster
<figure><img src="image_14.png" /><figcaption>(a) F vs U</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) ε<sub>d</sub> vs U</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12. PF-AT1 model — Damage at failure for the bar with a V-notch.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 13. PF-AT1 model – when $U = 0.026$ mm – for the bar with a V-notch.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 14. Before failure – notch zone ($x = 0.5, y = 0.05$) – $\bar{Y}$ vs $U$ – for the bar with a V-notch.</figcaption></figure>

computation, when compared to the original problem solved efficiently with *rsls*. In this case, the simplified problem may no longer be justified, as similar performance can be achieved without altering the problem formulation.

As previously discussed, the behavior of the phase-field models PF-AT1 and PF-CZM with respect to the history variable $H$ is essentially identical. Therefore, only the results of the PF-AT1 model have been presented in detail. For the PF-CZM model, only the relative errors are reported in Table 9.
<figure><img src="image_19.png" /><figcaption>Fig. 15. PF-AT1 model – Computational time (CPU) vs Finite element number – for the bar with a V-notch.</figcaption></figure>

<caption>**Table 9**</caption>

PF-CZM model — Relative errors $\Delta_r$ [%] between the original (reference) and simplified problems for the bar with a V-notch.

<table><caption>**Table 10**</caption><thead><tr><th>Model</th><th>Δ<sub>r</sub>(F<sub>c</sub>) [%]</th><th>Δ<sub>r</sub>(U<sub>c</sub>) [%]</th><th>Δ<sub>r</sub>(ε<sup>c</sup><sub>d</sub>) [%]</th></tr></thead><tbody><tr><td>PF-CZM</td><td>4.16</td><td>4.02</td><td>17.3</td></tr></tbody></table>

Phase-field model without damage threshold — relative errors $\Delta_r$, [%], between the original (reference) and simplified problems for the bar with a V-notch.

<table><thead><tr><th>Model</th><th>Δ<sub>r</sub>(F<sub>c</sub>)</th><th>Δ<sub>r</sub>(U<sub>c</sub>)</th><th>Δ<sub>r</sub>(ε<sup>c</sup><sub>d</sub>)</th></tr></thead><tbody><tr><td>PF-AT2</td><td>0.00</td><td>0.00</td><td>0.65</td></tr></tbody></table>

## Results for the phase field model without threshold damage (PF-AT2)

The results for the PF-AT2 model are presented in Fig. 16 and Table 10. The responses are identical. Due to the presence of a notch, a stress concentration occurs, and damage localizes at the notch tip. For the PF-AT2 model, damage initiates under infinitesimal loading across the entire structure. Therefore, the original and simplified problems yield the same result at failure ($d = 1$), as the damage field is naturally constrained between 0 and 1. Finally, given that the simplified problem is solved more efficiently, the use of the history variable $H$ appears to be an acceptable compromise for the PF-AT2 model in this test case.

After the first two tests, two main conclusions can be drawn. First, the original and simplified problems are equivalent for the PF-AT2 model. Second, for phase-field models with a damage threshold, the damage is more diffuse and evolves more rapidly in the simplified problem than in the original one. This is because a variational equality is solved in the simplified problem. As a result, at the local scale, damage triggered at one node naturally influences adjacent nodes. Moreover, the overestimation of the critical fracture energy and the earlier onset of failure are direct consequences of this behavior. In the case of the uni-axial bar test, the only noticeable difference is in the critical fracture energy. For the bar without a notch, failure occurs at the first and only damage initiation point. The parameters $H_0$ and $\ell$ are appropriately calibrated for a one-dimensional bar, which explains why the force-displacement responses are similar for both problems. In contrast, for the notched bar, damage evolves more gradually. Before reaching the failure state, the stress and strain fields are non-homogeneous due to the presence of the notch. Thus, at the local scale, damage increases progressively with each loading increment. As a result, the history variable influences the global force-displacement response.

Damage evolution depends on material properties, geometry, and applied loads. In the following tests, we focus on the influence of material properties on the accuracy of the simplified model. Specifically, we vary $\sigma_c$ and study the ratio between $U_0$ and $U_c$.
<figure><img src="image_20.png" /><figcaption>Fig. 16. Results of PF-AT2 model for the bar with V-notch.</figcaption></figure>

## 3.2. Plate with a hole

The plate with a hole, subjected to tensile loading as illustrated in Fig. 17, is studied. The plane stress hypothesis is considered. The material properties are given in Table 11. The boundary conditions are described in Fig. 17. A tensile displacement is applied in the y-direction, while displacements in the x-direction are fixed. The imposed displacement increments are $\Delta U = 1 \times 10^{-5}$ mm. The PF-CZM linear model is used. The regularization length $\ell$ is fixed at 0.025 mm. The mesh consists of 82 500 finite elements, resulting in a mesh size of 0.005 mm in the damage zone. (see Fig. 18).

The values of the failure strength are presented in Table 11. For each value considered, the ratio $U_0/U_c$ is computed based on the original problem (Fig. 19). Then, the relative errors $\Delta_r(F_c)$, $\Delta_r(U_c)$, and $\Delta_r(\varepsilon_d^c)$ between the original and simplified problems are determined. Finally, these relative errors are plotted against the $U_0/U_c$ ratio in Fig. 20(b). From Figs. 19 and 20(b), it can be observed that as the imposed failure strength increases, the ratio $U_0/U_c$ also increases, along with the relative errors. The relative errors associated with the critical force and displacement exhibit similar trends in both direction and magnitude and remain below 10%. However, the relative error in the critical fracture energy increases significantly, reaching up to 371%. Therefore, it is essential
<figure><img src="image_21.png" /><figcaption>Fig. 17. Plate with a hole — geometry definition (dimensions are in millimeters).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 18. Plate with a hole — mesh.</figcaption></figure>
<table><caption>**Table 11**
Plate with a hole — material properties.</caption><tr><td>Young's Modulus : E [MPa]</td><td>10,000</td></tr><tr><td>Poisson's ratio : v</td><td>0.3</td></tr><tr><td>Toughness : G<sub>c</sub> [kN/m]</td><td>1</td></tr><tr><td>Failure strength : σ<sub>c</sub> [MPa]</td><td>{100, 200, 300, 400, 500}</td></tr></table>

<figure><img src="image_23.png" /><figcaption>**Fig. 19.** PF-CZM model — $\frac{U_0}{U_c}$ vs $\sigma_c$ for the plate with a hole.</figcaption></figure>

<table><caption>**Table 12**
Hat-shaped specimen — material properties.</caption><tr><td>Young's Modulus : E [MPa]</td><td>2000</td></tr><tr><td>Poisson's ratio : v</td><td>0.3</td></tr><tr><td>Toughness : G<sub>c</sub> [kN/m]</td><td>0.025</td></tr><tr><td>Failure strength : σ<sub>c</sub> [MPa]</td><td>{20, 40}</td></tr></table>

not to evaluate accuracy based solely on force-displacement differences. To provide an overall view of the influence of the failure strength, Fig. 20(a) shows that the greater the $U_0/U_c$ ratio, the greater the discrepancy between the original and simplified problems. These results are specific to this test case. For phase-field models with a damage threshold, the influence of the $U_0/U_c$ ratio on the accuracy of the simplified problem compared to the original one is clearly illustrated.

### 3.3. Hat-shaped specimen

Consider a hat-shaped specimen, studied in the thesis by [22], subjected to compressive loading as illustrated in Fig. 21. The plane strain assumption is used. The boundary conditions are applied as shown in Fig. 21: the specimen is compressed in the y-direction, while displacements in the x-direction are fixed. In addition, a symmetry condition is applied along the vertical centerline of the specimen, allowing the problem size to be reduced by modeling only half of the geometry. The imposed displacement increments are $\Delta U = 2 \times 10^{-5}$ mm. The PF-CZM linear model is used, and the regularization length $\ell$ is set to 0.01 mm. The structure is discretized with a non-uniform mesh (Fig. 22), refined in the damage zone with a mesh size of $h = \frac{\ell}{8} = 1.25 \times 10^{-3}$ mm. The mesh consists of 96 236 finite elements.

The material properties used are listed in Table 12. Two values of $\sigma_c$ are considered: $\sigma_c = 20$ MPa and $\sigma_c = 40$ MPa.

The results are presented in Fig. 23. When $\sigma_c = 20$ MPa, the force-displacement responses are nearly identical. However, when $\sigma_c = 40$ MPa, the failure state is reached earlier in the simplified problem compared to the original problem. As previously mentioned, as the ratio $\frac{U_0}{U_c}$ increases, the relative errors also increase — particularly for the critical fracture energy (Table 13). In this test, the error in the critical fracture energy becomes significant when the ratio reaches 33%, whereas in the plate with a hole test, the error becomes significant when the ratio reaches 83%. These differences arise because the evolution of the elastic energy potential $\psi_0$ depends on the material properties, geometry, and loading conditions.
<figure><img src="image_24.png" /><figcaption>Fig. 20. Results of PF-CZM model for the plate with a hole.</figcaption></figure>

<table><caption>**Table 13**
PF-CZM model — relative errors $\Delta_r$, [%], for the hat-shaped specimen.</caption><thead><tr><th>σ<sub>c</sub></th><th>U<sub>0</sub>/U<sub>c</sub> [%]</th><th>Δ<sub>r</sub>(F<sub>c</sub>)</th><th>Δ<sub>r</sub>(U<sub>c</sub>)</th><th>Δ<sub>r</sub>(ε<sup>c</sup><sub>d</sub>)</th></tr></thead><tbody><tr><td>20</td><td>26.1</td><td>0.52</td><td>0.49</td><td>2.26</td></tr><tr><td>40</td><td>33.4</td><td>5.02</td><td>4.41</td><td>34.3</td></tr></tbody></table>

An important aspect to highlight is that, despite numerical inaccuracies, the simplified problem remains a valid option in certain cases. This is exemplified by the three-point bending test performed on a mortar sample [23], in which numerical results are compared with experimental data.
<figure><img src="image_25.png" /><figcaption>Fig. 21. Hat-shaped specimen — geometry definition (dimensions are in millimeters).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 22. Hat-shaped specimen — mesh.</figcaption></figure>

### 3.4. Experimental test : 3 points bending

In this subsection, the validity of the simplified problem is assessed for the three-point bending configuration using experimental results from [23]. A three-point bending test on a mortar sample is considered, as illustrated in Fig. 24. This test is based on the experimental data from the thesis of [23] and the numerical results from the study by [11].

The plane strain assumption is adopted, with a specimen thickness of 100 mm. The material properties of the mortar are listed in Table 14. The boundary conditions applied are shown in Fig. 24. A compressive displacement is applied in the y-direction, while the x-direction is left free. The regularization length is set to $\ell = 1.5$ mm. The imposed displacement increments are $\Delta U = 2 \times 10^{-4}$ mm. The PF-CZM Cornelissen model is used. The structure is discretized using a non-uniform mesh, as shown in Fig. 25, with refinement
<figure><img src="image_27.png" /><figcaption>Fig. 23. PF-CZM model – F vs U – for the hat-shaped specimen.</figcaption></figure>

<caption>**Table 14**</caption>

3 points bending sample — material properties [23].

<table><tr><td>Young's Modulus : E [MPa]</td><td>20,000</td></tr><tr><td>Poisson's ratio : v</td><td>0.2</td></tr><tr><td>Toughness : G<sub>c</sub> [kN/m]</td><td>0.113</td></tr><tr><td>Failure strength : σ<sub>c</sub> [MPa]</td><td>2.4</td></tr></table>

in the damage zone. The mesh size in the refined region is $h = \frac{\ell}{6} = 0.25$ mm. The mesh consists of 104 978 finite elements. The numerical results obtained are compared to the experimental data from [23].
<figure><img src="image_28.png" /><figcaption>Fig. 24. 3 points bending sample — geometry definition (dimensions are in millimeter) [23].</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 25. 3 points bending sample — mesh.</figcaption></figure>

<caption>Table 15</caption>

PF-CZM model — relative errors $\Delta_r$ [%], for the 3 points bending.

<table><caption>Table 16</caption><thead><tr><th>U<sub>0</sub>/U<sub>c</sub> [%]</th><th>Δ<sub>r</sub>(F<sub>c</sub>)</th><th>Δ<sub>r</sub>(U<sub>c</sub>)</th><th>Δ<sub>r</sub>(E<sup>c</sup><sub>d</sub>)</th></tr></thead><tbody><tr><td>9.6</td><td>0.02</td><td>0.3</td><td>0.8</td></tr></tbody></table>

PF-CZM model — Computational time for the 3 points bending.

<table><thead><tr><th>Problem</th><th>Original - ssls</th><th>Original - rsls</th><th>Simplified</th></tr></thead><tbody><tr><td>CPU</td><td>7 d 6 h 45 m 46 s</td><td>5 d 3 h 39 m 24 s</td><td>5 d 3 h 26 m 23 s</td></tr></tbody></table>

As shown in Fig. 26, the curves are nearly identical for both the original and simplified problems. Furthermore, the experimental results are well reproduced. With a ratio $U_0/U_c$ of only 9.6%, the damage threshold is almost negligible. Consequently, the simplified problem yields the same response as the original one (see Table 15). In addition, the computational cost is reported in Table 16. Solving the original problem using the numerical method Newton reduced-space active-set line search (rsls) results in a computational time comparable to that of the simplified problem. However, when the original problem is solved using the Newton semi-smooth line search (ssls) method, the simplified problem shows a clear advantage in terms of computational time, without a significant loss of accuracy for the case shown. Therefore, in situations where a numerical method such as rsls is not available, solving the simplified problem may be a meaningful and efficient alternative for this test case.

Finally, all the simulations have illustrated the influence of the history variable H on both the force-displacement and fracture energy responses. The final test aims to highlight the influence of H on the predicted crack path.

## 3.5. Micro-structure

Consider a micro-structure with several cavities representing porosities (35%), as illustrated in Fig. 27. The plane strain assumption is used. The boundary conditions are applied as shown in Fig. 27. A tensile displacement is applied in the y-direction,
<figure><img src="image_30.png" /><figcaption>Fig. 26. PF-CZM model — F vs U for the 3 points bending.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Fig. 27. Micro-structure — geometry definition (1 mm × 1 mm).</figcaption></figure>

**Table 17**

Micro-structure — material properties.

<table><tr><td>Young's Modulus : E [MPa]</td><td>1000</td></tr><tr><td>Poisson's ratio : v</td><td>0.3</td></tr><tr><td>Toughness : G<sub>c</sub> [kN/m]</td><td>0.01</td></tr><tr><td>Failure strength : σ<sub>c</sub> [MPa]</td><td>7.9</td></tr></table>

while displacements in the x-direction are fixed. The imposed displacement increments are $ΔU = 1×10^{-5}$ mm. The material properties is detailed in Table 17. The PF-AT2, PF-AT1, and PF-CZM linear models are used. The values of $ℓ$ are listed in Table 18. A uniform mesh with a mesh size of $h = \frac{ℓ_{AT2}}{5} = 3.375×10^{-3}$ mm is used, resulting in a total of 136 259 finite elements.
<caption>**Table 18**</caption>

Micro-structure — regularization length for the different phase field models.

<table><thead><tr><th></th><th>PF-AT1</th><th>PF-AT2</th><th>PF-CZM</th></tr></thead><tbody><tr><th>ℓ [mm]</th><td>0.06</td><td>0.016875</td><td>0.06</td></tr></tbody></table>

<figure><img src="image_32.png" /><figcaption>**Fig. 28.** PF-AT1 model for the micro-structure.</figcaption></figure>

### Results for the phase field models with damage threshold (PF-AT1 and PF-CZM)

The two approaches yield different responses, as shown in Figs. 28(a), 28(b), 30(a) and 30(b). An earlier failure state and an overestimation of the critical fracture energy are observed for the simplified problem. In addition, the crack paths differ between the two approaches. The damage field is more diffuse in the simplified problem (Figs. 29(a) and 31(a)), leading to crack branches that are not present in the original problem (Figs. 29(b) and 31(b)).
<figure><img src="image_33.png" /><figcaption>Fig. 29. PF-AT1 model — Damage and crack path for the micro-structure.</figcaption></figure>

<caption>**Table 19**</caption>

Phase field models with damage threshold — relative errors $\Delta_r$, [%], for the micro-structure.

<table><thead><tr><th></th><th>U<sub>0</sub> [%]</th><th>Δ<sub>r</sub>(F<sub>c</sub>)</th><th>Δ<sub>r</sub>(U<sub>c</sub>)</th><th>Δ<sub>r</sub>(ε<sup>c</sup><sub>d</sub>)</th></tr></thead><tbody><tr><td>PF-AT1</td><td>37.9</td><td>10.2</td><td>9.13</td><td>38.0</td></tr><tr><td>PF-CZM</td><td>41.1</td><td>6.93</td><td>4.73</td><td>25.4</td></tr></tbody></table>

As observed previously, the error in the critical fracture energy is the most significant for the simplified problem, as shown in Table 19.

### Results for the phase field model without damage threshold (PF-AT2)

The global responses of both problems for the PF-AT2 model are shown in Fig. 32. The relative errors, reported in Table 19, are very small. The crack paths are also similar (Fig. 33(b)). Therefore, for this phase-field model, the simplified problem represents a suitable and efficient compromise (see Table 20).

The difference between the original and simplified problem results has been evaluated across various fracture scenarios. In the simplified problem, early failure and an overestimation of the critical fracture energy have been observed. Moreover, the quality of the simplified problem's results shows a clear dependence on the ratio $U_0/U_c$. Below a certain threshold for this ratio, the simplified
<figure><img src="image_34.png" /><figcaption>Fig. 30. PF-CZM model for the micro-structure.</figcaption></figure>

<caption>Table 20</caption>

Phase field model without damage threshold — relative errors $\Delta_r$, [%], for the micro-structure.

<table><thead><tr><th></th><th>$\frac{U_0}{U_c}$ [%]</th><th>$\Delta_r(F_c)$</th><th>$\Delta_r(U_c)$</th><th>$\Delta_r(\varepsilon_d^c)$</th></tr></thead><tbody><tr><td>PF-AT2</td><td>0</td><td>0.006</td><td>0</td><td>3.78</td></tr></tbody></table>

problem yields results comparable to those of the original problem. However, this threshold is not straightforward to determine for general context. As illustrated in the plate with a hole and the hat-shaped specimen, this ratio is test-case dependent, and as such, the simplified problem can lead to a loss of accuracy compared to the original problem.

In addition, the computational advantage of the simplified problem is not always evident. In this work, two numerical methods are used to solve the original problem. The first, the Newton semi-smooth line search (SSLS), solves the damage problem for each
<figure><img src="image_35.png" /><figcaption>Fig. 31. PF-CZM model — Damage and crack path for the micro-structure.</figcaption></figure>

degree of freedom. The second, the Newton reduced-space active-set line search (rsls), separates the degrees of freedom into active and inactive sets, solving the damage problem only for the active set. As a result, the simplified problem reduces computational cost compared to the original problem solved with the ssls method. However, when compared to the original problem solved with the rsls method, the computational cost is roughly equivalent.

Therefore, when a numerical method such as rsls is available, the use of the simplified problem is not recommended, as it offers no computational benefit and may introduce a loss of accuracy.

4. Conclusions

After all the tests, a global summary is provided to address several key aspects.

First, for the phase-field model without a damage threshold (PF-AT2), the history variable introduced by [10] can be considered consistent, as its influence on the results is minor.

In contrast, for phase-field models with a damage threshold, such as PF-AT1 and PF-CZM, the history variable significantly affects the global responses, including the force-displacement behavior, fracture energy evolution, and, in some cases, the crack path. This can be explained by the resolution of a variational equality at each increment. This promotes faster and more diffuse damage evolution. In the case of the bar with a notch, this accelerated damage evolution is clearly observed. For the bar without a notch, the damage evolution is identical between the original and simplified problem, because the damage threshold H₀ is properly calibrated, and the fields remain homogeneous prior to failure. This results in a failure state defined by U = U₀(= Uₑ). Ultimately,
<figure><img src="image_36.png" /><figcaption>Fig. 32. PF-AT2 model for the micro-structure.</figcaption></figure>

using the simplified problem may lead to an overestimation of the critical fracture energy and a premature prediction of failure, particularly for models with a damage threshold.

Following this test, the discussion shifts to identifying the conditions under which the history variable has little or no influence. To this end, several simulations are carried out by varying the material property $\sigma_c$, the failure strength. So, the influence of the history variable is found to depend on the ratio $\frac{U_0}{U_c}$. This phenomenon is further illustrated in the plate with a hole and hat-shaped specimen tests. For the PF-AT1 and PF-CZM models, when the ratio $\frac{U_0}{U_c}$ falls below a certain threshold, the fracture energy and force-displacement responses from both the original and simplified problems are nearly identical. In this case, the use of the history variable remains consistent and reliable. However, when the ratio exceeds this threshold, the simplified problem leads to faster damage evolution and a more diffuse damage field. As a consequence, the failure state is reached prematurely, the critical fracture energy is overestimated, and, in some cases, the crack paths are altered, as observed in the micro-structure test.

Although the simplified problem may lead to a reduction in computational time in certain cases, this benefit is not systematic and depends strongly on the numerical method used to solve the original problem. Furthermore, the simplified problem can result
<figure><img src="image_37.png" /><figcaption>**Fig. 33.** PF-AT2 model — Damage and crack path for the micro-structure.</figcaption></figure>

in a noticeable loss of accuracy. Therefore, its use should be limited to specific contexts where advanced numerical solvers are not available and precision requirements are relaxed.

In light of these findings, a re-examination of the definition and formulation of the history variable $H$ represents a necessary step for those intending to pursue the use of the simplified problem. As it stands, this approach entails a risk of accuracy loss without offering a significant advantage in terms of computational cost. A new calibration strategy for the damage threshold $H_0$ could be proposed to improve the force–displacement response. However, addressing the overestimation of fracture energy caused by the enforced equality condition when $d = 0$ remains a significant challenge.

Until now, the previous conclusions have been established within the framework of quasi-static loading. Extending this analysis to a dynamic context appears highly relevant. Given that the validity of the simplified problem depends on the ratio $U_0/U_e$, it is reasonable to assume that the loss of accuracy may be less significant under dynamic conditions. To confirm this hypothesis, a complementary study under dynamic loading naturally constitutes the next step of this work.

CRediT authorship contribution statement

**Abdelmalek Barki:** Writing – review & editing, Writing – original draft, Validation, Software, Project administration, Methodology, Investigation, Conceptualization.
**Jihed Zghal:** Writing – review & editing, Writing – original draft, Validation, Supervision, Software, Project administration, Methodology, Investigation, Conceptualization.
**Laurent Gallimard:** Writing – review & editing, Validation, Supervision, Methodology, Investigation.
**Isabelle Bruant:** Writing – review & editing, Writing – original draft, Supervision, Methodology.
**Nicolas Chevaugeon:** Writing – review & editing, Validation, Methodology, Investigation, Conceptualization.
**Luc Davenne:** Writing – review & editing, Investigation.
Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors thank N. Moës and A. Ammar for the interesting discussions and recommendations.

Appendix A. Staggered algorithm used: original and simplified problem

**Algorithm 1 Staggered Algorithm (original and simplified)**

1: **for** each time step *k* **do**

2: Set current state at time *k*: ($u^k$, $d^k$)

3: Initialize iteration counter: *i* ← 0

4: Initialize the first iteration: ($u_0$, $d_0$) ← ($u^k$, $d^k$)

5: **while** Not converged: $\max(\mathbf{err}_u, \mathbf{err}_d) > \text{tol}$ and $i < i_{\max}$ **do**

6: $u_{i+1} \leftarrow \underset{u \in C^{k+1}}{\arg \min} \mathcal{E}(u, d_i)$

7: **if** Original problem **then**

8: $d_{i+1} \leftarrow \underset{d \in D^{k+1}}{\arg \min} \mathcal{E}(u_{i+1}, d)$

9: **else if** Simplified problem **then**

10: $H_{i+1} \leftarrow \max(\psi_0(u_{i+1}), H^k, H_0)$

11: $d_{i+1} \leftarrow \underset{d \in H^1(\Omega)}{\arg \min} \mathcal{E}(u_{i+1}, d; H_{i+1})$

12: **end if**

13: Compute $\mathbf{err}_u \leftarrow \|u_{i+1} - u_i\|_{L^\infty}$

14: Compute $\mathbf{err}_d \leftarrow \|d_{i+1} - d_i\|_{L^\infty}$

15: Update: $(u_i, d_i) \leftarrow (u_{i+1}, d_{i+1})$

16: $i \leftarrow i + 1$

17: **end while**

18: Set converged solution: ($u^{k+1}$, $d^{k+1}$) ← ($u_i$, $d_i$)

19: **if** Original problem **then**

20: Apply irreversibility condition

21: $d_{\min} \leftarrow d^{k+1}$

22: **else if** Simplified problem **then**

23: Apply history variable update rule

24: $H^{k+1} \leftarrow \max(\psi_0(u^{k+1}), H^k, H_0)$

25: **end if**

26: **end for**

Appendix B. Uni-dimensional bar

The analytical results from [11,24,25] of a uni-dimensional bar to tensile loading are presented here.

B.1. Global responses for a homogeneous state

<table><caption>**Table B.1**

Analytical expressions for a homogeneous state.</caption><thead><tr><th></th><th>ℓ</th><th>H<sub>0</sub></th><th>U<sub>c</sub></th><th>F<sub>c</sub></th></tr></thead><tbody><tr><td>PF-AT2</td><td>27EG<sub>c</sub>/256σ<sub>c</sub><sup>2</sup></td><td>0</td><td>√(G<sub>c</sub>/3ℓE)</td><td>E(1/G<sub>c</sub> + ℓEε<sub>c</sub><sup>2</sup>)<sup>2</sup>ε<sub>c</sub></td></tr><tr><td>PF-AT1</td><td>3EG<sub>c</sub>/8σ<sub>c</sub><sup>2</sup></td><td>3G<sub>c</sub>/16ℓ</td><td>√(3G<sub>c</sub>/8ℓE)</td><td>√(3G<sub>c</sub>E/8ℓ)</td></tr><tr><td>PF-CZM</td><td>-</td><td>σ<sub>c</sub><sup>2</sup>/2E</td><td>σ<sub>c</sub>/E</td><td>σ<sub>c</sub></td></tr></tbody></table>
<caption>**Table B.2**</caption>

Critical damage field $d_c$ and half localization band $D_c$ (non homogeneous state).

<table><thead><tr><th></th><th>d<sub>c</sub></th><th>D<sub>c</sub></th></tr></thead><tbody><tr><td>PF-AT2</td><td>e<sup>-|x|/ℓ</sup></td><td>∞</td></tr><tr><td>PF-AT1</td><td>{ 0 if |x| > 2ℓ <br> (1 - |x|/ℓ)<sup>2</sup> else</td><td>2ℓ</td></tr><tr><td>PF-CZM</td><td>{ 0 if |x| > π/2ℓ <br> 1 - sin |x|/ℓ else</td><td>π/2ℓ</td></tr></tbody></table>

## B.2. Damage and localization

## Appendix C. Error in critical fracture energy

[6] show that the fracture energy is overestimated in simulations based on finite element discretization. This numerical error is given by:

$$G_c^{err} = \frac{G_c h}{c \ell} \qquad (C.1)$$

where $h$ is the mesh size. This error must be considered when determining the critical fracture energy $\varepsilon_d^c$, which leads to:

$$\varepsilon_d^c = (G_c^{err} + G_c)\Gamma \qquad (C.2)$$

## Data availability

Data will be made available on request.

## References

[1] Griffith AA. The phenomena of rupture and flow in solids. Philos Trans R. Soc Lond Ser A 1921;221:163-98.

[2] Francfort GA, Marigo J-J. Revisiting brittle fracture as an energy minimization problem. J Mech Phys Solids 1998;46(8):1319-42. http://dx.doi.org/10.1016/S0022-5096(98)00034-9, URL https://www.sciencedirect.com/science/article/pii/S0022509698000349.

[3] Ambrosio Luigi, Tortorelli Vincenzo Maria. Approximation of functional depending on jumps by elliptic functional via $\gamma$-convergence. Comm Pure Appl Math 1990;43(8):999-1036. http://dx.doi.org/10.1002/cpa.3160430805, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/cpa.3160430805.

[4] Ambrosio Luigi, Tortorelli VM. On the approximation of free discontinuity problems. Boll. Unione Mat. Ital. VII. Ser B 1992;6(1):105-23.

[5] Bourdin B, Francfort GA, Marigo J-J. Numerical experiments in revisited brittle fracture. J Mech Phys Solids 2000;48(4):797-826. http://dx.doi.org/10.1016/S0022-5096(99)00028-9, URL https://www.sciencedirect.com/science/article/pii/S0022509699000289.

[6] Bourdin Blaise, Francfort Gilles A, Marigo Jean-Jacques. The variational approach to fracture. J Elasticity 2008/04/01;91(1):5-148. http://dx.doi.org/10.1007/s10659-007-9107-3.

[7] Chambolle A, Francfort GA, Marigo J-J. When and how do cracks propagate? J Mech Phys Solids 2009;57(9):1614-22. http://dx.doi.org/10.1016/j.jmps.2009.05.009, URL https://www.sciencedirect.com/science/article/pii/S0022509609000726.

[8] Pham Kim, Marigo Jean-Jacques. Approche variationnelle de l'endommagement : I. Les concepts fondamentaux. Comptes Rendus Mécanique 2010;338(4):191-8. http://dx.doi.org/10.1016/j.crme.2010.03.009, URL https://www.sciencedirect.com/science/article/pii/S1631072110000410.

[9] Pham Kim, Marigo Jean-Jacques. Approche variationnelle de l'endommagement : II. Les modèles à gradient. Comptes Rendus Mécanique 2010;338(4):199-206. http://dx.doi.org/10.1016/j.crme.2010.03.012, URL https://www.sciencedirect.com/science/article/pii/S1631072110000446.

[10] Miehe Christian, Hofacker Martina, Welschinger Fabian. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Comput Methods Appl Mech Engrg 2010;199(45):2765-78. http://dx.doi.org/10.1016/j.cma.2010.04.011, URL https://www.sciencedirect.com/science/article/pii/S0045782510001283.

[11] Wu Jian-Ying. A unified phase-field theory for the mechanics of damage and quasi-brittle failure. J Mech Phys Solids 2017;103:72-99. http://dx.doi.org/10.1016/j.jmps.2017.03.015, URL https://www.sciencedirect.com/science/article/pii/S0022509616308341.

[12] Molnár Gergely, Gravouil Anthony. 2D and 3D abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture. Finite Elem Anal Des 2017;130:27-38. http://dx.doi.org/10.1016/j.finel.2017.03.002, URL https://www.sciencedirect.com/science/article/pii/S0168874X16304954.

[13] Martínez-Pañeda Emilio, Golahmar Alireza, Niordson Christian F. A phase field formulation for hydrogen assisted cracking. Comput Methods Appl Mech Engrg 2018;342:742-61. http://dx.doi.org/10.1016/j.cma.2018.07.021, URL https://www.sciencedirect.com/science/article/pii/S0045782518303529.

[14] Lorentz Eric, Cuvilliez S, Kazymyrenko K. Convergence of a gradient damage model toward a cohesive zone model. Comptes Rendus Mécanique 2011;339(1):20-6. http://dx.doi.org/10.1016/j.crme.2010.10.010, URL https://www.sciencedirect.com/science/article/pii/S1631072110001671.

[15] Lopez-Pamies Oscar, Dolbow John E, Francfort Gilles A, Larsen Christopher J. Classical variational phase-field models cannot predict fracture nucleation. Comput Methods Appl Mech Engrg 2025;433:117520. http://dx.doi.org/10.1016/j.cma.2024.117520, URL https://www.sciencedirect.com/science/article/pii/S0045782524007746.

[16] Kumar Aditya, Bourdin Blaise, Francfort Gilles A, Lopez-Pamies Oscar. Revisiting nucleation in the phase-field approach to brittle fracture. J Mech Phys Solids 2020;142:104027. http://dx.doi.org/10.1016/j.jmps.2020.104027, URL https://www.sciencedirect.com/science/article/pii/S0022509620302623.

[17] Kamarei Farhad, Kumar Aditya, Lopez-Pamies Oscar. The poker-chip experiments of synthetic elastomers explained. J Mech Phys Solids 2024;188:105683. http://dx.doi.org/10.1016/j.jmps.2024.105683, URL https://www.sciencedirect.com/science/article/pii/S0022509624001492.
[18] Amor Hanen, Marigo Jean-Jacques, Maurini Corrado. Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments. J Mech Phys Solids 2009;57(8):1209-29. http://dx.doi.org/10.1016/j.jmps.2009.04.011, URL https://www.sciencedirect.com/science/article/pii/S0022509609000659.

[19] Alnæs Martin, Blechta Jan, Hake Johan, Johansson August, Kehlet Benjamin, Logg Anders, et al. The FEniCS project version 1.5. Arch Numer Softw 2015;3(100). http://dx.doi.org/10.11588/ans.2015.100.20553.

[20] Balay Satish, Abhyankar Shrirang, Adams Mark F, Brown Jed, Brune Peter, Buschelman Kris, et al. PETSc users manual. Argonne Natl Lab 2019;(ANL-95/11 - Revision 3.12). URL https://petsc.org/release/docs/manual/.

[21] Tanné E, Li T, Bourdin B, Marigo J-J, Maurini C. Crack nucleation in variational phase-field models of brittle fracture. J Mech Phys Solids 2018;110:80-99. http://dx.doi.org/10.1016/j.jmps.2017.09.006, URL https://www.sciencedirect.com/science/article/pii/S0022509617306543.

[22] Gmati Hela. phase field modelling of fracture of elastic and elasto-viscoplastic solid materials (Ph.D. thesis), (2020HESAE010). HESAM Université; 2020, URL https://pastel.hal.science/tel-03026401.

[23] Rots Jan Gerrit. Computational modeling of concrete fracture (Ph.D. thesis), Technische Hogeschool Delft; 1988.

[24] Pham Kim, Amor Hanen, Marigo Jean-Jacques, Maurini Corrado. Gradient damage models and their use to approximate brittle fracture. Int J Damage Mech 2011;20(4):618-52. http://dx.doi.org/10.1177/1056789510386852.

[25] Marigo Jean-Jacques, Maurini Corrado, Pham Kim. An overview of the modelling of fracture by gradient damage models. Meccanica 2016/12/01;51(12):3107-28. http://dx.doi.org/10.1007/s11012-016-0538-4.