# How to introduce an initial crack in phase field simulations to accurately predict the linear elastic fracture propagation threshold?

Flavien Loiseau, Veronique Lazarus

► To cite this version:

Flavien Loiseau, Veronique Lazarus. How to introduce an initial crack in phase field simulations to accurately predict the linear elastic fracture propagation threshold?. Journal of Theoretical, Computational and Applied Mechanics, 2025, pp.1-16. 10.46298/jtcam.15198. hal-04931758v2

HAL Id: hal-04931758

https://hal.science/hal-04931758v2

Submitted on 28 Jun 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
**Identifiers**

DOI 10.46298/jtcam.15198

HAL hal-04931758v2

**History**

Received Feb 6, 2025

Accepted Apr 23, 2025

Published Jun 30, 2025

**Associate Editor**

Phu NGUYEN

**Reviewers**

Tushar K. MANDAL

Md Masiur RAHAMAN

**Open Review**

HAL hal-05080459

**Supplementary Material**

Data

DOI 10.5281/zenodo.14749906

**Licence**

CC BY 4.0

©The Authors

# How to introduce an initial crack in phase field simulations to accurately predict the linear elastic fracture propagation threshold?

Flavien LOISEAU and Véronique LAZARUS

IMSIA, CNRS, EDF, ENSTA Paris, Institut Polytechnique de Paris, 91120 Palaiseau, France

Variational phase field fracture models are now widely used to simulate crack propagation in structures. A critical aspect of these simulations is the correct determination of the propagation threshold of pre-existing cracks, as it highly relies on how the initial cracks are implemented. While prior studies briefly discuss initial crack implementation techniques, we present here a systematic investigation. Various techniques to introduce initial cracks in phase field fracture simulations are tested, from the crack explicit meshing to the replacement by a fully damaged phase field, including different variants for the boundary conditions. Our focus here is on phase field models aiming to approximate, in the Γ-convergence limit, Griffith quasi-static propagation in the framework of Linear Elastic Fracture Mechanics. Therefore, a sharp crack model from classic linear elastic fracture mechanics based on Griffith criterion is the reference in this work. To assess the different techniques to introduce initial cracks, we rely on path-following methods to compute the sharp crack and the phase field smeared crack solutions. The underlying idea is that path-following ensures staying at equilibrium at each instant so that any difference between phase field and sharp crack models can be attributed to numerical artifacts. Thus, by comparing the results from both models, we can provide practical recommendations for reliably incorporating initial cracks in phase field fracture simulations. The comparison shows that an improper initial crack implementation often requires the smeared crack to transition to a one-element-wide phase band to adequately represent a displacement jump along a crack. This transition increases the energy required to propagate the crack, leading to a significant overshoot in the force-displacement response. The take-home message is that to predict the propagation threshold accurately and avoid artificial toughening, the crack must be initialized either setting the phase field to its damage state over a one-element-wide band or meshing the crack explicitly as a one-element-wide slit and imposing the fully cracked state on the crack surface.

**Keywords:** Variational phase-field approach, Initial crack implementation, Linear Elastic Fracture Mechanics, Griffith theory.

## 1 Introduction

Phase field fracture models based on the variational approach to fracture, proposed by Francfort and Marigo (1998) and implemented by Bourdin et al. (2000), have become increasingly popular in the fracture mechanics community. These models are based on the idea of Griffith (1920) that the mechanical state of a structure is governed by the balance between the elastic energy and the fracture surface energy. This idea can also be formalized as the minimization of the potential energy among any possible crack extension, that is any sharp discontinuity of the displacement field (Francfort and Marigo 1998). As the crack surface is not known a priori, performing minimization on this set of sharp cracks is generally not possible in practice. The phase field fracture models replace the sharp crack with a smeared crack represented by a phase field and use a regularization of the elastic and fracture energies based on the work of Ambrosio and Tortorelli (1990); Ambrosio and Tortorelli (1992). For variational phase field models, the regularized functional has been shown to Γ-converge towards its sharp counterpart (Bellettini and Coscia 1994; Bourdin 1999; Giacomini and Ponsiglione 2003; Chambolle 2004). This property means that
the results obtained using phase field fracture models, with a small-enough regularization length and sufficient refinement, converges towards the sharp crack model using Griffith theory.

Recent studies have highlighted the limitations of classical variational phase field models in describing crack nucleation from ab-nihilo (Kumar et al. 2020; Bažant et al. 2022; Lopez-Pamies et al. 2024). However, these models effectively handle crack propagation problems. To deal with the failure of a structure using a variational phase field model, it is therefore necessary to introduce an initial crack, as in standard Linear Elastic Fracture Mechanics (LEFM). But, how initial cracks are introduced in phase field simulations can significantly bias the observed mechanical response, in particular the propagation threshold. For instance, the work of Singh et al. (2016) shows that an overshoot of the peak force occurs in the force-displacement response of the structure without a specific treatment applied to the crack tip. Similarly, Kristensen et al. (2021) also observed an artificial increase of the critical energy release rate when introducing crack geometrically in the simulation domain. It underscores the importance of having a clear and reliable method to introduce initial cracks as it affects when and how initial cracks propagates.

The problem of implementing initial cracks in phase field fracture simulation is discussed in the literature, but in a scattered manner and often as a side study. A first technique, used by Klinsmann et al. (2015); Singh et al. (2016); Kristensen et al. (2021), among others, consists in explicitly embedding the initial cracks into the mesh by duplicating nodes, resulting in an infinitely sharp crack. This technique is largely employed in Finite Element fracture simulations (Anderson and Anderson 2017, Section 12.4.). However, its application to phase field fracture simulations raises several questions. Indeed, one challenge is determining the appropriate boundary condition to impose on the initial crack. For the displacement field, stress-free boundary conditions has to be used for unloaded crack. However, different boundary conditions in crack phase coexist in the literature. The phase field variable can be: (a) not set explicitly (zero Neumann boundary condition), (b) set to a fully cracked state at the crack tip, or (c) set to a fully cracked state on the whole crack. Similarly, Tanné et al. (2018) proposed to replace the initial crack with a slightly open sharp V-notch. A second technique consists in including the initial crack in the initial phase field. With this implementation technique, the initial crack is not meshed explicitly, but the crack is introduced instead by setting the initial phase field. More specifically, the crack phase is set to a fully cracked state along the initial crack (Nguyen et al. 2015; Makvandi et al. 2019; Lo et al. 2019; Yoshioka et al. 2020; Kristensen et al. 2021). A third technique can be applied in phase field models where a history variable governs the crack phase evolution (Miehe et al. 2010). The initial crack phase field is imposed indirectly through the history variable. This method has been proposed by Borden et al. (2012) and employed by Klinsmann et al. (2015); K. H. Pham et al. (2017); Liu et al. (2020), among others. As this method only applies to the phase field models using a history variable, this approach will only be briefly discussed in this work. For all those methods, a crucial factor is the initial crack thickness. This aspect is often unspecified in the literature (or even overlooked in some cases). However, it can significantly impact the numerical predictions of phase field fracture simulations as it will be shown here.

In this article, we explore systematically the impact of various initial crack implementation techniques on the mechanical response of structures. The sharp crack model will act as our reference as the phase field smeared crack model Γ-converges towards it. Our aim is to determine whether certain crack implementation technique introduce a bias into the mechanical results, in particular the critical load, and why. The comparison of the phase field fracture with the reference sharp crack model enables us to identify the most reliable technique to introduce initial cracks in phase field simulations. The findings will help users of the phase field method to choose the most accurate approach. Focus will be placed on the quasi-static propagation that satisfies Griffith criterion at each moment using a path-following method. For each crack growth increment, the corresponding critical load will be determined explicitly in sharp crack models or through a path-following constraint in phase field simulations. Following the quasi-static equilibrium path during crack propagation enables us to highlight the errors induced by an improper crack initialization irrespective of any dynamical effect.

This paper is organized as follows. Section 2 introduces the sharp crack model through the Griffith theory, and its formalization in the variational approach to fracture. The smeared crack model (i.e., the phase field fracture model) is then presented as a regularization of the
variational approach to fracture. Section 3 defines the quasi-static equilibrium path and explains how it is obtained, analytically and numerically, for the sharp and smeared crack. Section 4 presents various techniques to introduce initial cracks in phase field fracture simulation. Section 5 compares eight different techniques to implement initial cracks in phase field simulations on a simple problem (Single Edge Notched Tension test). The results obtained in this section enable us to discuss the different techniques and determine which ones are unbiased. The whole study has been done for straight cracks and phase field models without history variables. Section 6 briefly discusses how the results can be extended beyond. Finally, Section 7 concludes the study and provides practical recommendations to introduce initial cracks in phase field fracture simulations.

## 2 From Griffith theory to variational phase field models

Starting from the original idea of Griffith (1920), we recall the variational approach to fracture proposed by Francfort and Marigo (1998). This section covers the theoretical grounds of this work with an emphasis on the link between Griffith theory and variational phase field models for fracture. Bourdin et al. (2008) provides a more detailed discussion on this subject.

### 2.1 Griffith theory

For a static elastic problem, the displacement field of an elastic body submitted to external loads minimizes its potential energy $P = \mathcal{E} - W_{\text{ext}}$, where $\mathcal{E}$ is the body elastic energy and $W_{\text{ext}}$ is the external work. In a 2D setting, these quantities are defined by unit thickness. Griffith (1920) proposed to extend this minimization principle to study the rupture of a solid with a crack $\Gamma$ of length $a$. To this aim, Griffith (1920) introduced the surface energy $G_c$, also called the fracture energy, defined as the energy required to generate new free surfaces<fnref n="1" />. From this statement, the state of the solid at equilibrium must minimize the energy $P + G_c a$. Deriving this energy with respect to the crack length $a$ leads to the definition of the energy release rate

$$G = -\frac{\partial P}{\partial a} \quad \text{at constant load.} \qquad (1)$$

To ensure the stability of a given crack equilibrium state, $G \le G_c$ must be imposed, as $G > G_c$ implies the existence of an excess of energy, hence dynamical propagation. Moreover, introducing an irreversibility constraint and deriving this energy with respect to the crack surface give the following quasi-static crack propagation conditions:

$$(G - G_c)\dot{a} = 0, \quad G \le G_c, \quad \dot{a} \ge 0. \qquad (2)$$

It means that two distinct cases can occur. If $G < G_c$, then $\dot{a} = 0$, which means that the crack does not propagate. Conversely, if $G = G_c$, then $\dot{a} \ge 0$, indicating that the crack propagation becomes possible.

### 2.2 Variational approach to fracture

Francfort and Marigo (1998) extended this work by formalizing this idea and writing Griffith theory in a variational framework as follows.

Let us consider an elastic body $\Omega$ with an initial crack $\Gamma_0$. Body forces $f$ are prescribed in the bulk. The body boundary, denoted $\partial\Omega$, includes two distinct subsets: $\partial\Omega_u$ and $\partial\Omega_t$, where displacement $u$, respectively, surface forces are prescribed

$$u = u_{\text{imp}} \text{ on } \partial\Omega_u \text{ and } \sigma \cdot n = t \text{ on } \partial\Omega_t. \qquad (3)$$

Crack surfaces are assumed to be free of any external load. Denote $\mathcal{U}_\Gamma$ the space of admissible displacements defined as

$$\mathcal{U}_\Gamma = \{u \in H^1(\Omega \setminus \Gamma) \mid \forall x \in \partial\Omega_u, u(x) = u_{\text{imp}}(x)\}, \qquad (4)$$

<footnote n="1">¹ While considered as a reversible energy in (Griffith 1920), we consider here that $G_c$ includes all the energy dissipated during propagation.</footnote>
where H¹(Ω \ Γ) is a classic Sobolev space. In this case, the elastic energy and the external work are

$$ 
\begin{aligned}
\mathcal{E}(\mathbf{u}, \Gamma) &= \int_{\Omega \setminus \Gamma} \frac{1}{2} \varepsilon(\mathbf{u}) : \mathbf{E} : \varepsilon(\mathbf{u}) \, d\mathbf{x}, \\
\mathcal{W}_{\text{ext}}(\mathbf{u}) &= \int_{\Omega} f \cdot \mathbf{u} \, d\mathbf{x} + \int_{\partial\Omega_t} t \cdot \mathbf{u} \, d\mathbf{x},
\end{aligned}
\quad (5) $$

where $\varepsilon(\mathbf{u}) = \frac{1}{2}(\nabla\mathbf{u} + \nabla^T\mathbf{u})$ is the strain, $\mathbf{E}$ is the elasticity tensor. Griffith theory, in Equation (2), can be rewritten as the minimization problem

$$ (\mathbf{u}, \Gamma) = \arg \min_{\substack{\mathbf{u}' \in \mathcal{U}_f \\ \Gamma' \supseteq \Gamma_0}} \mathcal{E}(\mathbf{u}', \Gamma') - \mathcal{W}_{\text{ext}}(\mathbf{u}') + \int_{\Gamma'} G_c \, d\mathbf{x}. \quad (6) $$

The constraint $\Gamma' \supseteq \Gamma_0$ enforces the crack growth (i.e., the cracking irreversibility preventing cracks from healing). Note that in this framework, a crack is represented by a stress-free surface that can propagate in the domain $\Omega$. It induces a sharp discontinuity in the displacement field over which a displacement jump can occur.

## 2.3 Regularization as variational phase field model

As the crack corresponds to a discontinuity in the elastic body, the minimization problem in Equation (6) is a free-discontinuity problem, which is a well-studied problem (Ambrosio et al. 2000). In the context of image segmentation, free-discontinuity problems have been regularized so that sharp discontinuities are replaced by smeared band of non-vanishing width (Mumford and Shah 1989; Ambrosio and Tortorelli 1992). Francfort and Marigo (1998) proposed to apply the same regularization to the minimization problem in Equation (6). Applying the regularization replace the explicit description of the discontinuity surface $\Gamma$ with a new field: the (crack) phase field, denoted $\alpha(\mathbf{x})$. By convention, we consider that the material at coordinate $\mathbf{x}$ is uncracked if $\alpha(\mathbf{x}) = 0$ and fully crack if $\alpha(\mathbf{x}) = 1$. The minimization problem in Equation (6) governing the state of the domain becomes

$$ (\mathbf{u}, \alpha) = \arg \min_{\substack{\mathbf{u} \in \mathcal{U} \\ \alpha \in \mathcal{A}}} \mathcal{E}(\mathbf{u}, \alpha) + \mathcal{D}(\alpha) - \mathcal{W}_{\text{ext}}(\mathbf{u}) \quad (7) $$

where $\mathcal{U}$ and $\mathcal{A}$ are respectively the admissibility space of the displacement field $\mathbf{u}$ and of the crack phase field $\alpha$. The irreversibility of cracking is imposed through an irreversibility constraint on the phase field. Given an initial phase field $\alpha_0(\mathbf{x})$, the phase field can stay constant or only grow, leading to the admissibility space

$$ \mathcal{A} = \{ \alpha \in H^1(\Omega) \mid \forall \mathbf{x} \in \Omega, \alpha(\mathbf{x}) \ge \alpha_0(\mathbf{x}) \}. \quad (8) $$

The term $\mathcal{D}(\alpha)$ is the regularization of the dissipation integral over the (unknown) crack surface into a domain integral. The different energies are defined as

$$ \mathcal{E}(\mathbf{u}, \alpha) = \int_{\Omega} \frac{1}{2} a(\alpha) \varepsilon(\mathbf{u}) : \mathbf{E} : \varepsilon(\mathbf{u}) \, d\mathbf{x} \quad \text{and} \quad \mathcal{D}(\alpha) = \frac{G_c}{c_w} \int_{\Omega} \frac{w(\alpha)}{\ell} + \ell \nabla \alpha \cdot \nabla \alpha \, d\mathbf{x}. \quad (9) $$

In addition to the classic elastic parameters and the critical energy release rate $G_c$, the regularization introduces a new parameter: the regularization length $\ell$, to which the thickness of the crack phase profile is proportional. The function $a(\alpha)$ describes how the crack phase affects the initial elastic properties. The function $w(\alpha)$ describes the local dissipation due to the crack phase growth. The scaling constant $c_w$ is deriving from the choice of $w$. For the choice of these functions, we retain the AT1 model, proposed by K. Pham et al. (2011) and shown to be pertinent against experiments by Tanné et al. (2018). In the AT1 model, the functions are

$$ a(\alpha) = (1 - \alpha)^2, \quad w(\alpha) = \alpha, \quad c_w = 8/3. \quad (10) $$
The energy functional involved in the regularized minimization, in Equation (7), i.e. the smeared crack model, has been shown to Γ-converge towards its sharp crack counterpart, in Equation (6), when $\ell \to 0$ (Braides 1998). In practice, it means that the solution of Equation (7) converges to the solution of Griffith theory, in Equation (2), when $\ell \to 0$. Therefore, the convergence should also be recovered on the numerical solution when the numerical errors vanish.

### 3 Quasi-static equilibrium path

At this point, we presented two models for Griffith theory: the sharp crack model and the smeared crack model using phase field. As specified in the introduction, this work focuses on the study of crack propagation in quasi-static conditions. To track the whole crack propagation, the crack must propagate steadily: it must verify the Griffith criterion ($G = G_c$) throughout the crack propagation. The series of mechanical states verifying $G = G_c$ during the whole crack propagation is called the quasi-static equilibrium path. Following the equilibrium path provides a mean to simulate the whole crack propagation without instability. In this work, it enables capturing what happens at the beginning of crack propagation, where potential instabilities can occur. In practice, following this path consists in adapting the load during the crack propagation to stay at the crack propagation threshold. In this section, we briefly show how the quasi-static equilibrium path of the structure can be recovered, with both the sharp crack and the phase field models. More complete description and discussion of path-following methods in the context of phase field models for fracture will be proposed in another work dedicated to this subject (Loiseau and Lazarus 2025).

#### 3.1 Sharp crack model

To follow the quasi-static equilibrium path with the sharp crack model, the idea consists in adapting the load to verify the Griffith criterion $G = G_c$ throughout the crack propagation. We employ the same method as Triclot et al. (2024) to obtain the displacement field at equilibrium for different crack length. It is recalled in the following paragraph. Note that the crack path must be known a priori to apply this method.

Let us consider an elastic body $\Omega$ with an imposed displacement $u_{imp} = \lambda \bar{u}_{imp}$ on its boundary $\partial\Omega_u$, where $\lambda$ is the load factor (or amplitude). Linear elasticity ensures that the whole displacement field in the body is proportional to the load factor: $u(x) = \lambda \bar{u}(x)$, where $\bar{u}(x)$ is the displacement field induced by the unitary prescribed displacement $\bar{u}_{imp}$. The domain contains an initial crack of length $a(t = 0) = a_0$, which evolves during the (proportional) loading. The crack propagates when the energy release rate $G$ reaches the critical energy release rate $G_c$. As the energy release rate $G$ is quadratic in the displacement $u$, we can calculate it from the energy release rate $\bar{G}$ corresponding to $\lambda = 1$ writing $G = \lambda^2 \bar{G}$. Therefore, the critical load factor $\lambda$ can be recovered for each crack length $a$ by:

1. Solving the elastic problem to obtain the displacement field $\bar{u}$ for a unitary load ($\lambda = 1$),

2. Calculating the associated energy release rate $\bar{G}(a)$,

3. Calculating the critical load factor using Griffith criterion

$$ \lambda^2(a) \bar{G}(a) = G_c \implies \lambda(a) = \sqrt{\bar{G}(a)/G_c}. \qquad (11) $$

Knowing the critical loads for any crack length $a$ enables to recover the quasi-static equilibrium path $u(x, a) = \lambda(a) \bar{u}(x, a)$.

In practice, we employ Finite Element simulations to compute the displacement field $\bar{u}(x, a)$ for a unitary load ($\lambda = 1$) and a given crack length $a$. The energy release rate $\bar{G}$ is calculated using the $G(\Theta)$ method, originally proposed by Destuynder et al. (1981) and later discussed by Moran and Shih (1987); Suo and Combescure (1992). Algorithm 1 summarizes the methods and its numerical implementation.
**Require** Mesh containing a crack of length *a*, initial crack length *a*<sub>0</sub>, crack length increment Δ*a* and the maximum crack length *a*<sub>max</sub>.

1. Initialize the crack length *a* = *a*<sub>0</sub>

2. **while** *a* < *a*<sub>max</sub> **do**

3. Solve the elastic problem for a unitary loading (λ = 1) → **u**(*x*)

4. Compute the (unitary) energy release rate → **G**

5. Compute the critical load factor → λ = (G/G<sub>c</sub>)<sup>1/2</sup>

6. Compute the displacement field at equilibrium → **u**(*x*) = λ**u**(*x*)

7. Increment the crack length → *a* = *a* + Δ*a*

8. **end while**

**Algorithm 1** Computation of the quasi-static equilibrium path for the sharp crack model using Finite Element elastic simulations.

3.2 Smeared crack model

Due to the regularization of the crack discontinuity in phase field models, the method proposed
in the previous section can not be directly applied. To obtain the quasi-static equilibrium path
in phase field models, the idea consists in applying an indirect load control method (Rastiello
et al. 2022). Those methods, initially employed for geometric non-linearities (Wempner 1971;
Riks 1972; Riks 1979; Crisfield 1981), consist in introducing the load factor λ as an unknown in
the minimization problem, and its associated control equation. For crack propagation problem,
the control equation is arbitrarily chosen to limit crack propagation and enforces a steady
propagation. For instance, Gutiérrez (2004) limited the elastic energy release between two load
steps, whereas Singh et al. (2016) directly controlled the fracture dissipation (which is equivalent
to controlling the crack length). Various control equations have been reviewed by Rastiello et al.
(2022). In this work, we chose to control the load factor λ by imposing the maximum strain
increment (over the domain Ω) between two steps as proposed by Chen and Schreyer (1991).
It restricts the maximum strain increment in the domain between two load steps to limit the
crack growth. Given the state (u₀, α₀) at the previous load step, the minimization problem of
Equation (7) becomes

$$
\begin{equation}
\begin{aligned}
(\bar{u}, \alpha) &= \operatorname*{arg\,min}_{\substack{u \in \mathcal{U} \\ \alpha \in \mathcal{A}}} \mathcal{E}(\lambda \bar{u}, \alpha) + \mathcal{D}(\alpha) - \mathcal{W}_{\text{ext}}(\lambda \bar{u}) \\
\text{subject to} \quad & \max_{x \in \Omega} (\varepsilon(\lambda \bar{u}(x)) - \varepsilon(u_0)) = \Delta \varepsilon_{\text{imp}},
\end{aligned}
\tag{12}
\end{equation}
$$

where Δε<sub>imp</sub> is an arbitrary maximum strain increment which (indirectly) controls the load step.

To solve the problem in Equation (12), we extend the classic solution method of phase
field fracture models. Let us start by describing how the problem without path-following in
Equation (7) is solved. As the functional is poly-convex, the minimization problem can be split
into two sub-problems on the displacement field and the crack phase (Bourdin et al. 2000). The
solution method consists of alternatively solving those two sub-problems, hence the alternate
minimization name. The equation governing the stationary points of each sub-problem is
obtained by calculating the directional derivatives of the energy functional and equating them
to zero for any direction. It gives the variational formulation for the two sub-problems. From
a practical point of view, the discretization of the problem and its resolution are performed
via the Finite Element Method using FEniCSx (Alnæs et al. 2014; Baratta et al. 2023). The
irreversibility constraint on the phase field is imposed through variational inequality using the
solver vinnewonrsls from the PETSc library.

When introducing the path-following constraint in Equation (12), the problem can still
be solved alternatively. Instead of solving the displacement sub-problem for given boundary
conditions, it is solved with unitary boundary conditions (λ = 1 = ⇒ **u**<sub>imp</sub> = **u**<sub>imp</sub>) giving the
displacement field **u**(*x*). Then, the control equation is solved to get the critical load factor λ,
and the displacement field is rescaled: **u**(*x*) = λ**u**(*x*). The phase field sub-problem remains
unchanged. The algorithm to solve one load step of this problem is summarized in Algorithm 2.
**Require** Domain mesh, initial state ($u_0(x), \alpha_0(x)$), load increment $\Delta\varepsilon_{imp}$.

1 **while** not converged **do**

2 Solve the elastic sub-problem for a unitary loading ($\lambda = 1) \rightarrow \bar{u}(x)$

3 Solve the control equation to obtain the critical load factor $\lambda \rightarrow \lambda$

4 Compute the displacement field at equilibrium: $\rightarrow u(x) = \lambda \bar{u}(x)$

5 Solve the phase field sub-problem to get the phase field $\rightarrow \alpha(x)$

6 **end while**

**Algorithm 2** Computation of a load step of the quasi-static equilibrium path for the smeared crack model.

## 4 Initial crack implementation techniques in phase field fracture

This section presents the different crack initialization techniques for phase field fracture simulations that will be compared and analyzed. The techniques are classified into two categories: the geometric initial cracks (explicitly represented in the mesh) and the phase initial cracks (embedded in the initial phase field). For each category, different variants of the technique will be proposed. An identifier is associated with each technique variant for clarity purposes. Figure 1 provides a summary and illustrates each of the variants considered in this work.

### 4.1 Geometric initial crack (GEO)

The geometric initial cracks are specified in the geometry (i.e., in the mesh) of the domain. This technique corresponds to incorporating initial cracks as sharp (non-regularized) discontinuities with traction-free boundary conditions, as it is usually done in fracture mechanics.

Different variants of the geometric technique will be tested, starting with changes of the crack thickness. Two thicknesses are compared: an infinitely thin crack (To) and a one-element-wide crack (T1). Note that infinitely thin crack (To) corresponds to duplicating nodes along the crack path, as often employed in fracture mechanics Finite Element simulations. For all geometric initial cracks, the displacement boundary condition is set as zero Neumann boundary condition, corresponding to traction-free crack surfaces. Then, the effect of boundary conditions on the initial crack $\alpha$ is also investigated. Three different boundary conditions in crack phase $\alpha$ are considered: zero Neumann boundary conditions (NEU), crack phase imposed to 1 at the crack tip (TIP), and crack-phase imposed to 1 along whole crack (WHL).

### 4.2 Phase field initial cracks (PHA)

The phase field initial cracks are obtained by initializing the crack phase field. It corresponds to adding the initial cracks to the domain as regularized discontinuities.

In this technique, the width of initial crack will also be investigated for the phase field initial crack. Two variants are considered: the case of the infinitely thin initial crack (To), and the case of a one-element-wide initial crack (T1).

Let us describe how the phase field initial crack are implemented in practice. For clarity purposes, we consider the case of a single straight initial crack. Curved initial cracks are discussed in Section 6. To implement phase field initial cracks, two steps are required. First, the nodal values of the phase field are set to one along the initial crack. The case of an infinitely thin (To) and one-element-wide (T1) initial crack must be distinguished. For the (To) initial crack, the phase field $\alpha(x) = 1$ is prescribed at the nodes along the initial crack. It requires the mesh nodes to be rigorously aligned along the initial crack. For a (T1) initial crack, the phase field $\alpha(x) = 1$ is prescribed on the nodes of elements crossed by the initial crack. In this case, the elements must be aligned to represent a straight crack accurately so that the initial crack goes through their centers. The second step consists in regularizing the sharp crack to make it smeared. This can be done by solving the phase field sub-problem (see step 5 of Algorithm 2) with zero loads.

Figure 2 illustrates how the initial phase field is defined for both variants. This figure shows the phase field for infinitely thin and one-element-wide initial cracks. The values before and after regularization are shown. The nodal values of the phase field are represented as square markers in the figure to highlight the prescribed field values before the regularization.
<figure><img src="image_2.png" /><figcaption>**Figure 1** Illustration of the different crack initialization techniques. The green line on the GEO-To illustrates corresponds to the line where the node are duplicated. The color corresponds to the phase field with blue being uncracked ($\alpha = 0$) and red being fully cracked ($\alpha = 1$). The abbreviation BC corresponds to boundary condition.</figcaption></figure>

# 5 Comparison of initial crack implementation techniques

## 5.1 Presentation of the numerical benchmark

The Single Edge Notched Tensile (SENT) test, which is often used a numerical benchmark for phase field simulations, is employed to compare the different crack initialization techniques. The geometry and boundary conditions, illustrated in Figure 3, correspond to those used by Singh
<figure><img src="image_3.png" /><figcaption>**Figure 2** Illustration of the definition of phase field initial crack (PHA). The red color corresponds to the crack ($\alpha$) whereas the blue color corresponds to the uncracked parts ($\alpha = 0$).</figcaption></figure>

et al. (2016). The SENT specimen is a square with a side length $L = 1$ mm and an initial crack of length $a_0 = L/2 = 0.5$ mm. The material is linear elastic with a Young modulus $E = 230.77$ GPa, a Poisson ratio $\nu = 0.43$, a critical energy release rate $G_c = 2700$ J m$^{-2}$, and a regularization length $\ell = 0.015$ mm. In this section, we use two types of mesh for the phase field simulations. First, an ideal structured mesh with elements aligned with the initial crack and the expected crack path is used to concentrate fully on the influence of the initial crack implementation technique. The techniques satisfyingly representing initial cracks with this mesh are then confronted with an unstructured mesh (still verifying the conditions required to apply the technique). For the sake of brevity, the meshes are represented in the results. The mesh size is set to $h = \ell/6$ in the region $y \in [L/2 - 2\ell, L/2 + 2\ell]$, where $y$ is the coordinate along the vertical direction. Larger elements with a size $h_{far} = 2\ell$ are used outside this region.

<figure><img src="image_4.png" /><figcaption>**Figure 3** Geometry of the SENT specimen and boundary conditions.</figcaption></figure>

Following the recommendations of Bourdin et al. (2008), the critical energy release rate in the phase field simulations is adjusted to obtain the correct effective energy release rate

$$G_c^{\text{eff}} = G_c^{\text{PF}} \left(1 + \frac{h}{c_w \ell}\right). \qquad (13)$$

Hence, in the phase field simulations, we choose the critical energy release rate $G_c^{PF}$ so that the effective critical energy release rate $G_c^{\text{eff}} = G_c \implies G_c^{PF} = 2541$ J m$^{-2}$.

To compare the techniques, we will consider the force-displacement curves from the sharp and smeared crack models in Figures 4(a), 5(a) and 6(a). In these curves, the displacement $u$ is the vertical component of the imposed displacement $u_{imp}$. The force $F$ is the vertical component of the reaction force $F$ calculated by integrating the traction vector (obtained by post-processing the displacement field $u$) over the upper face.
## 5.2 Infinitely thin initial cracks (To)

We start by presenting the results with infinitely thin initial cracks (variants with To). The comparison will be based on two types of observations both represented in Figure 4. The first type of observations, in Figure 4(a), is the force-displacement curves. The force-displacement obtained with the different initial crack implementation techniques are represented along with the reference curve of the sharp crack model. The second type of observation is the phase field around the initial crack tip after the propagation. They are provided in Figures 4(b) to 4(e).

<figure><img src="image_5.png" /><figcaption>(a) Force-displacement curves for an infinitely thin initial crack (To). The inset illustrates thickening of the crack phase band occurring for PHA-To during the snap-back.</figcaption></figure>

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

**Figure 4** Comparison of the phase field fracture solutions with the different crack initialization techniques for an infinitely thin initial crack (To). The top part shows the force-displacement curves along with the reference solution from the sharp crack model. The lower part shows the crack phase field after crack propagation with blue being uncracked ($\alpha = 0$) and red being fully cracked ($\alpha = 1$).

Before commenting the peak and post-peak part of Figure 4(a), let us note that the elastic phase is biased by the smeared crack model. Indeed, the diffusion of the crack induces a loss of macroscopic rigidity. Thus, it induces a small decrease of the slope in the elastic part. This loss decreases when reducing the regularization length $\ell$. This applies for all the following results (except for the PHA-To technique, which is discussed later). Moreover, our simulations are stopped when the elastic energy reaches 0.2% of its maximum value. Hence, the load does not drop to zero after failure because the crack has not fully propagated; a small ligament remains. As we focus on the peak load, this choice does not impact our analysis.

The first observation, in Figure 4(a), is that all the geometric methods (GEO) present a peak on the force-displacement curve, which does not appear on the reference solution. As the path-following method employed here enables to properly capture the snap-back, we can conclude that this peak is a numerical artifact. This type of overshoot has already been observed directly in the work of Singh et al. (2016) and indirectly in the work of Klinsmann et al. (2015).
The peak reduces when the crack phase is set to 1 either on the crack tip or on the whole crack boundary. Nevertheless, it persists with all the geometric methods, even when the crack phase is imposed to one on the whole crack (WHL). In Figures 4(b) to 4(d), it can be observed that the crack phase is initially equal to one on a (duplicated) line of nodes. However, after the propagation, the crack phase equals one on a line of elements (i.e., two lines of nodes). Indeed, for the displacement “discontinuity” to occur in an element, the whole element needs to be cracked. Thus, a slight bifurcation of the crack is necessary to transition from an infinitely thin initial crack to a one-element-wide crack that enables the displacement jump. This transition requires an artificial excess of energy to propagate the infinitely thin initial crack compared to a proper initial crack.

The case of the PHA initial crack is different. The associated force-displacement curve in Figure 4(a) shows a Young modulus that is initially too high and a snap-back is observed in what is supposed to be the elastic phase. As shown in the inset images of Figure 4(a), the observed snap-back corresponds to the thickening of the initial crack. The reason is that a whole line of elements must be cracked for the displacement “discontinuity” to occur. Thus, the initial crack field needs to thicken to allow the displacement jump across the crack. Once this thickening is achieved, the remainder of the curves matches the reference results, and no force overshoot is observed.

From those first results, we can conclude that infinitely thin initial cracks should not be used in phase field fracture. It means that a duplication of the nodes on the crack faces is not sufficient to adequately represent an initial crack (even when setting the crack phase to one along the crack). Note that the proposition of Tanné et al. (2018) to introduce the initial crack as a sharp V-notch will show similar issues if the notch tip is one node. Indeed, the problem of having to transition from one node to a line of cracked elements persists.

## 5.3 One-element-wide initial cracks (T1)

Let us now present the results for one-element-wide initial cracks. The results are shown in Figure 5. The Figure 5(a) provides a quantitative comparison of the initial crack implementation techniques through the force-displacement curves. Reference results obtained via the sharp crack model are also reported in this figure. The phase fields around the initial crack tip after propagation are also provided in Figures 5(b) to 5(e).

The technique GEO-T1-NEU shows a significant overshoot on the force-displacement curve compared to the reference solution. This artificial overshoot is due to the use of the Neumann boundary conditions. Indeed, it imposes that the derivative of the phase field normal to the crack is null, forcing multiple elements to crack at once to start the crack propagation. Once again, it generates an artificial increase in the energy required for the crack to propagate.

The GEO-T1-TIP also shows a (smaller) overshoot on the force-displacement curve. In Figure 5(c), we observe a cone shape at the beginning of the crack propagation. It can be compared to other methods where the overshoot is not observed, e.g. GEO-T1-WHL in Figure 5(d). One possible explanation is that the regularization region of the phase is not developed as it would for a propagated crack. What happens before the crack tip likely affects the energy required to propagate the initial crack, hence the crack propagation threshold.

The crack initialization techniques GEO-T1-WHL and PHA-T1 can not be distinguished from each other in terms of mechanical response and phase field. Their force-displacement curves have the same shape as the reference solution. Moreover, we can observe, in Figures 5(d) and 5(e), that the initial crack and the propagated part of the crack are indistinguishable for both implementation techniques. To find which method is the best choice, we investigated the numerical performances of both methods. The number of iterations required for the alternate minimization to converge gives an indication of the overall numerical performances of both methods. The geometric initialization technique cumulates 11 055 iterations, whereas the phase initialization technique has 11 734. While the geometric method has a lower total number of iterations, no method clearly outperforms the other.
<figure><img src="image_7.png" /><figcaption>(a) Force-displacement curves for a one-element-wide initial crack (T1).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) GEO-T1-NEU.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) GEO-T1-TIP.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) GEO-T1-WHL.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(e) PHA-T1.</figcaption></figure>

**Figure 5** Comparison of the phase field fracture solutions with the different crack initialization techniques for a one-element-wide initial crack (T1). The top part shows the force-displacement curves along with the reference solution from the sharp crack model. The lower part shows the crack phase field after crack propagation with blue being uncracked ($α = 0$) and red being fully cracked ($α = 1$).

## 5.4 Extension to unstructured mesh

The crack initialization techniques GEO-T1-WHL and PHA-T1 properly represent an initial crack on a structured mesh. Now, we want to check whether this result extends to unstructured meshes. Figure 6 shows the results obtained with unstructured mesh. We observe that the force-displacement curves in Figure 6(a) perfectly superimpose. While more oscillations in the post-peak can be observed, the peak force remains mostly unchanged. Note that the oscillations are due to the irregular mesh preventing the crack propagation in a perfect straight line, as shown in Figures 6(c) and 6(e). Those oscillations decrease with mesh refinement (Loiseau and Lazarus 2025).

# 6 Discussion on the extensions to more complex cases

This study only tackles the case of rectilinear initial cracks and phase field models not relying on a history variable. Extensions to curved cracks and models governed by a history variable can be envisaged and are discussed in the following paragraphs.

The extension of GEO-T1-WHL and PHA-T1 to curved initial cracks is primarily a technical difficulty. For the geometric technique, the crack boundaries (either internal or external) must align with the crack path, and the crack itself must remain one-element wide. For the phase field initial crack PHA-T1, the nodal values of the phase field must be set to $α(\mathbf{x}) = 1$ on all the nodes of elements intersecting with the crack path. Ideally, the mesh nodes should be aligned with the initial crack so that the crack path passes through the centers of these elements.
<figure><img src="image_12.png" /><figcaption>(a) Force-displacement curves for a one-element-wide initial crack.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) GEO-T1-WHL (before propagation).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) GEO-T1-WHL (after propagation).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) PHA-T1 (before propagation).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(e) PHA-T1 (after propagation).</figcaption></figure>

**Figure 6** Comparison of the phase field fracture solutions on unstructured meshes with the crack initialization techniques GEO-T1-WHL and PHA-T1. The top part shows the force-displacement curves along with the reference solution from the sharp crack model. The fields on the lower part show the phase field before and after propagation for both techniques with blue being uncracked ($α = 0$) and red being fully cracked ($α = 1$).

For a phase field model governed by a history variable, such as the model proposed by Miehe et al. (2010), all the proposed techniques can be applied. However, the initial phase field can also be indirectly imposed through the history variable as it drives the crack phase growth. As this method does not apply to every phase field model, it has not been studied in this study. However, it corresponds to indirectly setting the initial crack phase through the history variable. The outcome should be similar to that achieved using the initial phase field approach (PHA).

# 7 Conclusion

This study assesses various techniques for incorporating initial cracks in variational phase field fracture simulations so that it provides a good approximation of Griffith theory. Different crack initialization techniques for phase field fracture simulations are compared to a sharp crack model, towards which the phase field model should theoretically Γ-converges. This study shows that representing initial cracks as infinitely thin, for instance by duplicating nodes, leads to invalid propagation thresholds in phase field simulations. Indeed, this approach induces peaks in the force-displacement curves and the critical load is significantly overestimated. This conclusion is expected as, for phase field fracture simulations using continuous finite elements, cracks must be at least one element wide to properly represent the displacement jump across the crack. Therefore, any initial crack must also be at least one element wide to maintain consistency and to avoid the artificial excess of energy required to transition to a one-element-wide zone. Moreover, this study shows that the crack phase field must be set to 1 (fully cracked state) along the entire
initial crack boundary to avoid artificial toughening effects.

This study enables us to identify two unbiased techniques to incorporate initial cracks in phase field fracture simulations:

* Geometric technique (GEO-T1-WHL): This approach consists of embedding a one-element-wide crack in the geometry (i.e., in the mesh) and imposing the crack phase field to the fully cracked state ($\alpha = 1$) on it,

* Initial phase field technique (PHA-T1): This method initialized the phase field to the fully cracked state on the elements that the initial crack crosses.

Both techniques give equivalent results for a similar numerical cost, so both methods can be employed indifferently.

## References

Alnæs, M. S., A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells (2014). Unified form language: A domain-specific language for weak formulations of partial differential equations. *ACM Transactions on Mathematical Software* 40(2):1–37. [DOI], [ARXIV].

Ambrosio, L., N. Fusco, and D. Pallara (2000). *Functions of Bounded Variation and Free Discontinuity Problems*. Oxford University Press. [DOI].

Ambrosio, L. and V. M. Tortorelli (1990). Approximation of functional depending on jumps by elliptic functional via t-convergence. *Communications on Pure and Applied Mathematics* 43(8):999–1036. [DOI].

Ambrosio, L. and V. M. Tortorelli (1992). On the Approximation of Free Discontinuity Problems. *Bollettino Della Unione Matematica Italiana* 6B(1):105–123.

Anderson, T. L. and T. L. Anderson (2017). *Fracture mechanics: fundamentals and applications*. 4th ed. CRC Press. [DOI].

Baratta, I. A., J. P. Dean, J. S. Dokken, M. Habera, J. S. Hale, C. N. Richardson, M. E. Rognes, M. W. Scroggs, N. Sime, and G. N. Wells (2023). *DOLFINx: The next generation FEniCS problem solving environment*. [DOI], [OA].

Bažant, Z. P., H. T. Nguyen, and A. Abdullah Dönmez (2022). Critical Comparison of Phase-Field, Peridynamics, and Crack Band Model M7 in Light of Gap Test and Classical Fracture Tests. *Journal of Applied Mechanics* 89(6). [DOI], [OA].

Bellettini, G. and A. Coscia (1994). Discrete approximation of a free discontinuity problem. *Numerical Functional Analysis and Optimization* 15(3-4):201–224. [DOI].

Borden, M. J., C. V. Verhoosel, M. A. Scott, T. J. Hughes, and C. M. Landis (2012). A phase-field description of dynamic brittle fracture. *Computer Methods in Applied Mechanics and Engineering* 217-220:77–95. [DOI].

Bourdin, B., G. Francfort, and J.-J. Marigo (2000). Numerical experiments in revisited brittle fracture. *Journal of the Mechanics and Physics of Solids* 48(4):797–826. [DOI], [HAL].

Bourdin, B. (1999). Image segmentation with a finite element method. *ESAIM: Mathematical Modelling and Numerical Analysis* 33(2):229–244. [DOI], [OA].

Bourdin, B., G. A. Francfort, and J.-J. Marigo (2008). The Variational Approach to Fracture. *Journal of Elasticity* 91(1-3):5–148. [DOI], [HAL].

Braides, A. (1998). *Approximation of Free-Discontinuity Problems*. Vol. 1694. Lecture Notes in Mathematics. Berlin, Heidelberg: Springer Berlin Heidelberg. [DOI].

Chambolle, A. (2004). An approximation result for special functions with bounded deformation. *Journal de Mathématiques Pures et Appliquées* 83(7):929–954. [DOI], [OA].

Chen, Z. and H. Schreyer (1991). Secant structural solution strategies under element constraint for incremental damage. *Computer Methods in Applied Mechanics and Engineering* 90(1-3):869–884. [DOI].

Crisfield, M. (1981). A fast incremental/iterative solution procedure that handles “snap-through”. *Computers & Structures* 13(1-3):55–62. [DOI].

Destuynder, P., P. E. M. Djaoua, L. Chesnay, and J. C. Nedelec (1981). Sur une interprétation mathématique de l’intégrale de Rice en théorie de la rupture fragile. *Mathematical Methods in the Applied Sciences* 3(1):70–87. [DOI].

Francfort, G. and J.-J. Marigo (1998). Revisiting brittle fracture as an energy minimization problem.
*Journal of the Mechanics and Physics of Solids* 46(8):1319–1342. [DOI], [HAL].

Giacomini, A. and M. Ponsiglione (2003). A Discontinuous Finite Element Approximation of Quasi-static Growth of Brittle Fractures. *Numerical Functional Analysis and Optimization* 24(7-8):813–850. [DOI], [ARXIV].

Griffith, A. A. (1920). The phenomena of rupture and flow in solids. *Philosophical Transactions of the Royal Society of London. Series A* 221(582-593):163–198. [DOI], [HAL].

Gutiérrez, M. A. (2004). Energy release control for numerical simulations of failure in quasi-brittle solids. *Communications in Numerical Methods in Engineering* 20(1):19–29. [DOI].

Klinsmann, M., D. Rosato, M. Kamlah, and R. M. McMeeking (2015). An assessment of the phase field formulation for crack growth. *Computer Methods in Applied Mechanics and Engineering* 294:313–330. [DOI].

Kristensen, P. K., C. F. Niordson, and E. Martínez-Pañeda (2021). An assessment of phase field fracture: crack initiation and growth. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* 379(2203):20210021. [DOI], [OA].

Kumar, A., B. Bourdin, G. A. Francfort, and O. Lopez-Pamies (2020). Revisiting nucleation in the phase-field approach to brittle fracture. *Journal of the Mechanics and Physics of Solids* 142:104027. [DOI], [OA].

Liu, J., X. Liang, Y. Xue, Y. Fu, K. Yao, and F. Dou (2020). Investigation on crack initiation and propagation in hydraulic fracturing of bedded shale by hybrid phase-field modeling. *Theoretical and Applied Fracture Mechanics* 108:102651. [DOI].

Lo, Y.-S., M. J. Borden, K. Ravi-Chandar, and C. M. Landis (2019). A phase-field model for fatigue crack growth. *Journal of the Mechanics and Physics of Solids* 132:103684. [DOI], [OA].

Loiseau, F. and V. Lazarus (2025). *Path-following methods for quasi-static crack propagation simulated by variational phase field approach*. Preprint. [HAL].

Lopez-Pamies, O., J. E. Dolbow, G. A. Francfort, and C. J. Larsen (2024). *Classical variational phase-field models cannot predict fracture nucleation*. [DOI], [OA].

Makvandi, R., S. Duczek, and D. Juhre (2019). A phase-field fracture model based on strain gradient elasticity. *Engineering Fracture Mechanics* 220:106648. [DOI], [HAL].

Miehe, C., M. Hofacker, and F. Welschinger (2010). A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. *Computer Methods in Applied Mechanics and Engineering* 199(45-48):2765–2778. [DOI].

Moran, B. and C. Shih (1987). Crack tip and associated domain integrals from momentum and energy balance. *Engineering Fracture Mechanics* 27(6):615–642. [DOI].

Mumford, D. and J. Shah (1989). Optimal approximations by piecewise smooth functions and associated variational problems. *Communications on Pure and Applied Mathematics* 42(5):577–685. [DOI], [OA].

Nguyen, T., J. Yvonnet, Q.-Z. Zhu, M. Bornert, and C. Chateau (2015). A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure. *Engineering Fracture Mechanics* 139:18–39. [DOI], [HAL].

Pham, K. H., K. Ravi-Chandar, and C. M. Landis (2017). Experimental validation of a phase-field model for fracture. *International Journal of Fracture* 205(1):83–101. [DOI].

Pham, K., H. Amor, J.-J. Marigo, and C. Maurini (2011). Gradient Damage Models and Their Use to Approximate Brittle Fracture. *International Journal of Damage Mechanics* 20(4):618–652. [DOI], [HAL].

Rastiello, G., H. L. Oliveira, and A. Millard (2022). Path-following methods for unstable structural responses induced by strain softening: a critical review. *Comptes Rendus. Mécanique* 350(G2):205–236. [DOI], [OA].

Riks, E. (1972). The Application of Newton's Method to the Problem of Elastic Stability. *Journal of Applied Mechanics* 39(4):1060–1065. [DOI].

Riks, E. (1979). An incremental approach to the solution of snapping and buckling problems. *International Journal of Solids and Structures* 15(7):529–551. [DOI].

Singh, N., C. Verhoosel, R. de Borst, and E. van Brummelen (2016). A fracture-controlled path-following technique for phase-field modeling of brittle fracture. *Finite Elements in Analysis and Design* 113:14–29. [DOI], [OA].

Suo, X.-Z. and A. Combescure (1992). On the application of G(Θ) method and its comparison with
De Lorenzi's approach. *Nuclear Engineering and Design* 135(2):207-224. [DOI], [HAL].

Tanné, E., T. Li, B. Bourdin, J.-J. Marigo, and C. Maurini (2018). Crack nucleation in variational phase-field models of brittle fracture. *Journal of the Mechanics and Physics of Solids* 110:80-99. [DOI], [OA].

Triclot, J., T. Corre, A. Gravouil, and V. Lazarus (2024). Toughening effects of out-of-crack-path architected zones. *International Journal of Fracture*. [DOI], [HAL].

Wempner, G. A. (1971). Discrete approximations related to nonlinear theories of solids. *International Journal of Solids and Structures* 7(11):1581-1599. [DOI].

Yoshioka, K., D. Naumov, and O. Kolditz (2020). On crack opening computation in variational phase-field models for fracture. *Computer Methods in Applied Mechanics and Engineering* 369:113210. [DOI].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

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

**Authors' contributions** F. Loiseau and V. Lazarus contributed equally to the conceptualization and methodology of the study. F. Loiseau was responsible for data collection, initial data analysis, and drafting the original manuscript. Both authors contributed to data interpretation and manuscript revision. V. Lazarus supervised the research and secured funding. Both authors reviewed and approved the final version of the manuscript.

**Supplementary Material** The associated dataset is available at the permalink 10.5281/zenodo.14749906. This directory contains the images and CSV files used to generate the figures in the accompanying article. These processed data files allow for reproducing the plots and analyses presented in the paper.

**Funding** The work was supported by Agence de l'Innovation de Défense - AID - via Centre Interdisciplinaire d'Etudes pour la Défense et la Sécurité - CIEDS - (projects 2022 - FracAddi), and by Agence Nationale de la Recherche (ANR-23-CE51-0054 3FAM).

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.