# Understanding regularized crack initiation through the lens of Finite Fracture Mechanics

Aurélien Doitrand, Gergely Molnár

► To cite this version:

Aurélien Doitrand, Gergely Molnár. Understanding regularized crack initiation through the lens of Finite Fracture Mechanics. International Journal of Fracture, 2025, 249 (1), pp.12. 10.1007/s10704-024-00837-9. hal-04890063

HAL Id: hal-04890063

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

Submitted on 16 Jan 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
# Understanding regularized crack initiation through the lens of Finite Fracture Mechanics

Aurelien Doitrand<fnref n="1" /> and Gergely Molnár<fnref n="2" />

As a remedy to pathological sharp crack configurations such as strong singularities or anti-plane shear cracks, where crack initiation is driven solely by energy, a regularized crack description can be adopted to study crack initiation. The nucleation of a regularized crack at a V-notch is studied using the coupled criterion through matched asymptotic expansions. The process zone around the crack is described by crack regularization usually employed in phase-field models. The effective crack length increases with increasing regularization length so that the incremental energy release rate decreases, which in turn increases the critical generalized stress intensity factor at initiation. Decreasing incremental energy release rate is also obtained with increasing Poisson's ratio. For a given material characteristic length, it is shown that the initiation crack length only depends on the V-notch angle and Poisson's ratio. For a given geometry and Poisson's ratio, the initiation length is proportional to the regularization length. The proposed description of regularized crack initiation shows good correspondence to the generalized stress intensity factor obtained by phase-field calculation, the only difference being in the description of the process zone prior to crack initiation.

Keywords Coupled criterion, regularized crack, matched asymptotic approach

## 1 Introduction

Phase field (PF) fracture models were originally developed to describe crack propagation by solving Griffith's approach (Griffith 1921; Griffith 1924) revisited as a global minimization problem (Francfort et al. 1998). This is done through a regularized description of the crack by a diffused damage zone which extent is controlled by a regularization length. This length is an additional input to the critical energy release rate (ERR) which is already used in linear elastic fracture mechanics (LEFM). PF models thus naturally retrieve Griffith's approach when LEFM assumptions are verified, i.e., for a sufficiently large pre-existing crack under pure opening mode (Molnár et al. 2022). In addition, these models also have the ability to deal with crack initiation, which is one of the main limitations of LEFM since in this case, not only energy but also a stress criterion must be considered (Leguillon 2002). An illustrative example of the latter assertion is crack initiation at a notch that has a V-shape with an angle $\beta$. When $\beta = 180$ deg. (i.e., the case of a crack-less flat edge), there is no singularity and a strength criterion prevails to study failure and crack initiation. On the contrary, a $\beta = 0$ deg. angle corresponds to a crack for which propagation may be assessed by an energy criterion. For intermediate $0 < \beta < 180$ deg. angles, crack initiation cannot be predicted employing solely a stress or an energy criterion. In this case, both criteria must be satisfied at a finite length to assess the V-notch crack initiation.

A particular case of crack initiation is when homogeneous stress fields are considered. For a given stress triaxiality, a value of the material strength can be defined. It results in the definition of the material strength surface, which thus represents the failure stresses for different triaxiality ratios in the principal stress or strain space. Under uniaxial tension, the PF regularization length can be determined so as to achieve nucleation at a stress magnitude corresponding to the material tensile strength, thus interpreting the length scale as a material parameter (Amor et al. 2009; Borden et al. 2012; Mesgarnejad et al. 2015; Pham et al. 2017). A reasonable comparison between PF predictions and experiments can then be achieved by setting the regularization length based on the material tensile strength only, for instance to assess crack nucleation at a notch (Tanné et al. 2018). Nevertheless, this description may only be valid for a specific uniaxial loading, i.e.,

<footnote n="1">¹ Univ Lyon, INSA Lyon, Université Claude Bernard Lyon 1, CNRS UMR5510, MATEIS, F-69621 Villeurbanne, France, aurelien.doitrand@insa-lyon.fr</footnote>

<footnote n="2">² CNRS, INSA Lyon, LaMCoS, UMR5259, 69621 Villeurbanne, France</footnote>
one particular point of the strength surface. In turn, such a choice of the regularization length may result in incorrect failure strength prediction for other triaxial loading conditions. As a matter of example, significantly different regularization lengths were identified for the same material loaded under different biaxiality ratios by Abaza et al. (2022).

It was previously shown that under homogeneous stress state, the regularization length can explicitly be written as the product between Irwin's length and a function depending on both Poisson's ratio and the principal stress ratios (Molnár et al. 2020; Molnár et al. 2022). Moreover, it was found that this relation remains verified even under opening, in-plane or out-of-plane shear mode loadings (Molnár et al. 2020; Molnár et al. 2024). Based on this finding, one way to define the strength surface in PF models consists in adapting the regularization length depending on the local stress state to follow the desired strength surface (Doitrand et al. 2023b). Kumar et al. (2020) explicitly defined the strength surface as an input parameter of the PF model, in addition to the critical ERR and the regularization length. This was done by adding an external driving force describing the strength surface. Reasonable comparison between the PF simulations and experiments was obtained. It required, for each confrontation to an experiment, the calibration of an additional parameter allowing to retrieve Griffith's approach in case when LEFM assumptions are fulfilled.

Another alternative to define the fracture strength surface in a PF model consists in adapting the energy decomposition (De Lorenzis et al. 2022; Vicentini et al. 2024) choosing for instance volumetric-deviatoric (Amor et al. 2009), spectral (Miehe et al. 2010b) or no-tension (Freddi et al. 2010) decompositions. It was shown that a given energy decomposition calibrated to a fixed strength surface, still conserving the relation between the regularization length and the strength magnitude for this strength surface shape. In particular, an energy decomposition was proposed to retrieve Drucker-Prager strength surface (De Lorenzis et al. 2022) by adding a parameter controlling the shear or compressive to tensile strength ratio.

All the approaches mentioned above mainly rely on the idea of a coupling between energy and strength aspects to assess crack nucleation. This idea was actually the cornerstone of the so-called coupled criterion (CC) developed by Leguillon (2002) to address the shortcomings of LEFM in evaluating crack initiation. So far, the coupled criterion appears as an efficient model for understanding the details of crack initiation (Weißgraeber et al. 2016; Doitrand et al. 2024). It is particularly emphasized when the CC is implemented through matched asymptotic expansions (Leguillon 2002; Yosibash et al. 2006; García et al. 2012; Leguillon et al. 2017; Felger et al. 2019; Doitrand et al. 2020; Jiménez-Alfaro et al. 2022; Doitrand et al. 2023a) as it yields direct access to the local loading in terms of Generalized Stress Intensity factors (GSIF) and provides a clear understanding of the stress and energy contribution to crack initiation. Previous works focused on comparing the CC to other approaches involving a regularized crack description, such as the Thick Level Set (TLS) Zghal et al. (2018). It was shown that TLS and CC provided close apparent initiation strengths for V-notches or cavities. The CC and PF approaches were also compared, for instance for crack initiation under Hertz indentation (Strobl et al. 2017; Strobl et al. 2020), in the case of thin ply laminates cracking (Reinoso et al. 2017) or bending failure of notched specimens (Abaza et al. 2022). Both models were also confronted under mode I, mode II and mode III crack initiation conditions (Molnár et al. 2020; Molnár et al. 2024). It was shown that the correlation between the regularization length and Irwin's length determined by the homogeneous PF solution could also be extended to configurations where the damage gradient is not negligible, evidencing that the PF approach also fulfills both energy and strength criteria. These confrontations highlighted the question whether initial phase-field boundary conditions should be prescribed or not to accurately represent crack initiation at a singularity. For instance, it was evidenced that initial phase field boundary conditions were required in order to retrieve Griffith's solution in the case of a sufficiently long crack (Molnár et al. 2020) or to retrieve the CC in the case of a V-notch (Doitrand et al. 2023b).

In some configurations, the initiation loading predicted by the CC is insensitive to the tensile strength, for instance in cases corresponding to strong singularities, i.e., when the ERR is a decreasing function of the crack length (Leguillon et al. 2000; Leguillon et al. 2012; Aranda et al. 2023). In these configurations, the IERR tends to infinity as the crack length vanishes and the stress criterion is always satisfied. Another example is a semi-infinite sharp crack under anti-plane
shear loading, for which the stress criterion does not influence the initiation loading (Doitrand et al. 2023a). Since the energy release rate is maximum for rectilinear propagation (Mittelman et al. 2014; Mittelman et al. 2015), it fails to reproduce the experimentally observed inclined facets (Pham et al. 2016). In these configurations, the CC reverts solely to an energy criterion and the effect of length disappears. In these kinds of configuration, an ingredient is missing to accurately capture the crack topology (Doitrand et al. 2023a; Molnár et al. 2024), which could possibly be overcome by regularization of the sharp crack. The objective of this work is thus to make a step forward in the understanding of regularized crack initiation through the scope of the coupled criterion. The initiation of a regularized crack ahead of a sharp V-notch tip is studied by means of matched asymptotic expansions. The configuration under investigation is presented in Section 2. Section 3 focuses on incremental energy release rate calculation based on matched asymptotic expansions. Section 4 presents the CC application to a regularized crack. The influence of the V-notch angle and Poisson's ratio is finally given in Section 5, together with a comparison to PF calculations.

# 2 Regularized crack

## 2.1 Crack initiation at a V-notch

We study the problem of 2D regularized crack initiation at a V-notch under plane strain conditions and pure opening mode in a homogeneous isotropic brittle material (Young's modulus $E$, Poisson's ratio $\nu$, critical ERR $G_c$, tensile strength $\sigma_c$), the V-notch angle is noted $\beta$ (Figure 1). The material is assumed to follow Rankine's strength surface (i.e., with a stress state independent critical tensile stress). Instead of a purely sharp description, the regularized crack is represented as a diffused

<figure><img src="image_2.png" /><figcaption>Figure 1: a) Representation of the regularized crack as a process zone (hatched zone described by the regularization length $l_c$) in which the damage variable $d$ varies between 1 (fully damaged) and 0 (pristine material) extending around a sharp crack (nominal length $l$) at a V-notch tip (angle $\beta$). b) Dimensionless regularization crack length $\xi = l_c/l$ when the coordinates are normalized by the sharp crack length (normalized length 1).</figcaption></figure>

damage zone around a sharp crack of nominal length $l$ Figure 1(a). The diffuse damage zone around the sharp crack is defined using the crack regularization usually employed in PF models (regularization length $l_c$), by prescribing $\hat{d} = 1$ PF boundary conditions along the sharp crack faces. The definition of the damage diffusion around the sharp crack ahead of the V-notch tip is done by solving the phase-field problem without the elastic energy. The diffuse damage topology results in a process zone around the crack where the damage variable $d$ varies between 1 (at the crack faces) and 0 (out of the process zone). The effective crack length $l_{eff}$ of the regularized crack is defined as:

$$l_{eff} = \int_{\Omega} \frac{1}{c_{\omega}l_{c}} (\omega(d) + l_{c}^{2}|\nabla d|^{2})d\Omega, \quad (1)$$

where $c_{\omega} = 8/3$ and $\omega(d) = d$. In this study we use the AT1 PF formulation. Implementation details can be found in the paper (Molnár et al. 2022). The dimensionless effective length
$\xi_{\text{eff}} = \ell_{\text{eff}}/\ell$ can be computed as a function of the dimensionless regularization length $\xi = \ell_c/\ell$. It can be noted that $\xi_{\text{eff}} \to 1$ when $\xi \to 0$, or equivalently $\ell_{\text{eff}} \to \ell$ when $\ell_c \to 0$, thus retrieving the sharp crack description for vanishing regularization length and process zone.

## 2.2 Regularization length

The proposed model enables implementing the coupled criterion to study the initiation of a regularized crack. It provides us with the option of using 3 parameters (i.e., critical energy release rate, tensile strength and regularization length). But the regularization length can also be chosen so as to retrieve a strength surface described in the principal stress space by $\sigma_f(\sigma_{II}/\sigma_I)$ (Molnár et al. 2020; Molnár et al. 2022) following the relation:

$$ \sigma_f \left( \frac{\sigma_{II}}{\sigma_I} \right) = \eta \left( v, \frac{\sigma_{II}}{\sigma_I} \right) \sqrt{\frac{E G_c}{\ell_c}}, \quad (2) $$

where $\eta$ is the normalized tensile strength which depends on the principal stress ratio and on the Poisson's ratio. Numerical implementations to compute $\eta$ are available in (Molnár et al. 2020; Molnár et al. 2022). In the case of a Rankine strength surface, $\sigma_f(\sigma_{II}/\sigma_I) = \sigma_c \forall \sigma_{II}/\sigma_I$, so that the regularization length can be written as a function of Irwin's length $\ell_{\text{mat}} = \frac{EG_c}{\sigma_c^2}$:

$$ \ell_c = \eta \left( v, \frac{\sigma_{II}}{\sigma_I} \right)^2 \ell_{\text{mat}}. \quad (3) $$

The function $\eta$ can be computed by determining the principal stress ratio given the V-notch angle. Since we focus on opening mode crack initiation, the principal stress ratio is computed along the V-notch bisector. It is thus independent of the distance to the V-notch tip and only depends on the V-notch angle and singularity exponent $\lambda$. The principal stress ratio writes (Yosibash et al. 2006):

$$ \frac{\sigma_{II}}{\sigma_I} = \frac{(1 - \lambda) \sin((2\pi - \beta)(1 - \lambda)/2) + (3 - \lambda) \sin((2\pi - \beta)(1 + \lambda)/2)}{-(1 - \lambda) \sin((2\pi - \beta)(1 - \lambda)/2) + (1 + \lambda) \sin((2\pi - \beta)(1 + \lambda)/2)}. \quad (4) $$

The characteristic exponent $\lambda$, solution of the equation $\sin((2\pi - \beta)\lambda) + \lambda \sin(2\pi - \beta) = 0$ (Leguillon et al. 1987), is given as a function of the V-notch angle in Figure 2(a).

<figure><img src="image_3.png" /><figcaption>Figure 2: a) Mode I characteristic exponent, b) Principal stress ratio ($\sigma_{II}/\sigma_I$) and c) function $\eta$ linking the regularization length and Irwin's length for different Poisson's ratios $v$, as a function of the V-notch angle $\beta$.</figcaption></figure>

Figure 2(b) shows the ratio $\sigma_{II}/\sigma_I$ as a function of the V-notch angle, which is a decreasing function varying from 1 for $\beta = 0$ deg. (equibiaxial stress state for a crack under tension) to 0 for $\beta = 180$ deg. (uniaxial stress state for a flat edge under tension). It determines the variation of $\eta$ as a function of the V-notch angle and Poisson's ratio values (Figure 2(c)). In particular, for $\beta = 180$ deg., the value of $\eta$ are the one corresponding to uniaxial tension. Given the material characteristic length $\ell_{\text{mat}}$, the V-notch angle (or equivalently the principal stress ratio) and the Poisson's ratio, the characteristic length is finally determined using Equation (3).
# 3 Incremental energy release rate

## 3.1 Matched asymptotic expansion

In the matched asymptotic approach, two configurations are studied. The first configuration is at the specimen scale, in which the regularized crack that initiates can be neglected since it is small with respect to the specimen dimensions. The second configuration focuses on a zone close to the initiation crack, disregarding the whole specimen geometry. The displacement fields are sought for both configurations and then matched in an intermediate zone to obtain the full solution considering both the whole specimen geometry and the small crack. We thus consider a two-scale problem to solve under plane strain conditions and linear elasticity. The first problem is written at the scale of the specimen. The displacement field $U^{\ell, \ell_c}$ (the superscripts $\ell$ and $\ell_c$ refer to the presence of a regularized crack, i.e., a process zone described by $\ell_c$ around a sharp crack of nominal length $\ell$) is the solution of the following set of equations:

$$ 
\left\{ 
\begin{array}{l}
-\nabla \cdot \underline{\sigma}(U^{\ell, \ell_c}) = 0, \\
\underline{\sigma}(U^{\ell, \ell_c}) = \underline{C} : \nabla U^{\ell, \ell_c}, \underline{C} \text{ is the stiffness tensor,} \\
\underline{\sigma}(U^{\ell, \ell_c}) \cdot \mathbf{n} = \mathbf{0} \text{ along the free edges, } \mathbf{n} \text{ is the normal vector to the free edges.}
\end{array}
\right. 
\qquad (5) 
$$

It is assumed that the initiation crack length $\ell$ and the process zone length $\ell_c$ are in the same order of magnitude, and relatively small compared to the specimen dimensions, this assumption has to be checked afterwards to ensure the validity of the matched asymptotic approach. Sufficiently far from the V-notch tip, the displacement and stress distributions are not much influenced by the presence of the small regularized crack. Therefore, the actual solution can be approximated by the solution without the regularized crack, which only requires a small correction to be brought close to the V-notch tip:

$$ U^{\ell, \ell_c}(x_1, x_2) = U^0(x_1, x_2) + \text{small correction}, \qquad (6) $$

where $U^0(x_1, x_2)$ is the solution of an idealized problem without crack. The small correction is actually determined by solving the second problem, close to the crack tip, detailed in the sequel. It decreases to 0 when $\ell$ tends to 0 and when moving away from the V-notch tip. This approximation constitutes the outer field, which is valid except near the V-notch tip. In the case of a V-notch, the mode I asymptotic displacement field depends on a characteristic exponent $\lambda$ and a corresponding mode $\mathbf{u}(\theta)$, that are obtained solving an eigenvalue problem (Leguillon et al. 1987), as $K_I r^\lambda \mathbf{u}(\theta) + o(r^\lambda)$. $K_I$ represents the V-notch Generalized Stress Intensity Factor (GSIF). The displacement field in the V-notch specimen can thus be written as:

$$ U^0(x_1, x_2) = U^0(0, 0) + K_I r^\lambda \mathbf{u}(\theta) + o(r^\lambda) \qquad (7) $$

To have a detailed form of the actual solution $U^{\ell, \ell_c}$, the initial domain is rescaled by $1/\ell$. The new dimensionless space variables are defined as $y_i = x_i/\ell$. As $\ell$ tends to 0, we obtain the inner domain, which is unbounded and in which the dimensionless crack length is 1 and the dimensionless regularization length is defined as $\xi = \ell_c/\ell$. The inner domain is depicted in Figure 1(b). The actual solution is assumed to expand in the following way:

$$ U^{\ell, \ell_c}(x_1, x_2) = U^{\ell, \ell_c}(\ell y_1, \ell y_2) = F_0(\ell) V^0(y_1, y_2, \xi) + F_1(\ell) V^1(y_1, y_2, \xi), \qquad (8) $$

with

$$ \lim_{\ell \to 0} \frac{F_1(\ell)}{F_0(\ell)} = 0. \qquad (9) $$

The $V^1$, which form the inner field, are solutions of the problems prescribed at infinity. They must match with the behaviour of the far field at infinity. As a consequence, there is an intermediate zone, close to the V-notch tip in the outer expansion and far from it the inner expansion where
both expansions given in Eq. (6) and Eq. (8) hold true. The matching of the terms in Eq. (7) and Eq. (8) results in :

$$ 
\left\{
\begin{aligned}
F_0(\ell) &= 1, \\
\mathbf{V}^0(y_1, y_2, \xi) &= \mathbf{U}^0(0, 0), \\
F_1(\ell) &= K_I \ell^\lambda, \\
\mathbf{V}^1(y_1, y_2, \xi) &\approx \gamma^\lambda \mathbf{u}(\theta),
\end{aligned}
\right.
\qquad (10) $$

with the normalized polar coordinate $\gamma = r/\ell$. The symbol $\approx$ means that $\mathbf{V}^1$ behaves like $\gamma^\lambda \mathbf{u}(\theta)$ at infinity, it can thus be written as:

$$ \mathbf{V}^1(y_1, y_2, \xi) = \gamma^\lambda \mathbf{u}(\theta) + \hat{\mathbf{V}}^1(y_1, y_2, \xi). \qquad (11) $$

It is necessary to prove that $\hat{\mathbf{V}}^1(y_1, y_2, \xi)$ exists and verifies the equilibrium equations. By combining Eqs. (8) and (10) into Eq. (5), and noting that the V-notch faces remain stress free in the inner domain, $\hat{\mathbf{V}}^1$ is solution to the following problem:

$$ 
\left\{
\begin{aligned}
-\nabla_y \cdot \hat{\underline{\sigma}} &= 0 \quad \text{where } \nabla_x = \frac{1}{\ell} \nabla_y, \\
\hat{\underline{\sigma}} &= \underline{C} : \nabla_y \hat{\mathbf{V}}^1, \\
\hat{\underline{\sigma}} \cdot \mathbf{n} &= 0 \quad \text{along the V-notch edges,} \\
\hat{\underline{\sigma}} \cdot \mathbf{n} &= -(\underline{C} : \nabla_y \gamma^\lambda \mathbf{u}(\theta)) \cdot \mathbf{n} \quad \text{along the crack edges,} \\
\hat{\mathbf{V}}^1 \quad &\text{decreases at infinity.}
\end{aligned}
\right.
\qquad (12) $$

There is a unique solution with finite energy to this system of equations (Laguillon et al. 1987) (extension of Lax-Milgram theorem to unbounded domains). The finite energy solution implies that the solution decreases to 0 at infinity. Finally, the expansion writes:

$$ \mathbf{U}^{\ell, \ell_c}(x_1, x_2) = \mathbf{U}^{\ell, \ell_c}(\ell y_1, \ell y_2) = \mathbf{U}^0(0, 0) + K_I \ell^\lambda [\gamma^\lambda \mathbf{u}(\theta) + \hat{\mathbf{V}}^1(y_1, y_2, \xi)]. \qquad (13) $$

The change in elastic strain energy due to regularized crack initiation is thus written as (Laguillon et al. 1987; Labossiere et al. 1999)

$$ -\delta W_{el} = \Psi(\mathbf{U}^{\ell, \ell_c}(x_1, x_2), \mathbf{U}^{0,0}(x_1, x_2)), \qquad (14) $$

where

$$ \Psi(f, g) = \frac{1}{2} \int_{\Gamma} [\underline{\sigma}(f) \cdot \mathbf{n} \cdot g - \underline{\sigma}(g) \cdot \mathbf{n} \cdot f] ds. \qquad (15) $$

The path independent integral $\Psi$ is defined on $\Gamma$, a closed contour surrounding the studied crack initiation location, starting and finishing on the V-notch free faces. The inwards normal to this contour is noted $\mathbf{n}$. The elastic strain energy variation thus writes:

$$ -\delta W_{el} = \frac{K_I^2}{E} \ell^{2\lambda} A(\xi, \nu), \qquad (16) $$

where $A(\xi, \nu)$ depends on the dimensionless regularization length, V-notch angle and Poisson's ratio.

The incremental energy release rate (IERR) for a regularized crack is defined as :

$$ \mathcal{G}_{\text{inc}} = \frac{-\delta W_{\text{el}}}{\ell_{\text{eff}}} = \frac{K_I^2}{E} \ell^{2\lambda} \frac{A(\xi, \nu)}{\ell_{\text{eff}}} = \frac{K_I^2}{E} \ell^{2\lambda-1} \frac{A(\xi, \nu)}{\xi_{\text{eff}}(\xi)} = \frac{K_I^2}{E} \ell^{2\lambda-1} A_{\text{eff}}(\xi, \nu). \qquad (17) $$

The dimensionless IERR $A_{\text{eff}} = A(\xi, \nu)/\xi_{\text{eff}}$ can be obtained by computing the elastic strain energy difference between uncracked and regularized cracked states in the inner domain for a given dimensionless regularization length, which is detailed in next section.
## 3.2 Numerical calculation of $A_{eff}$

The dimensionless IERR $A_{eff}$ is computed using the Finite Element (FE) method in the inner domain which is artificially bounded at a distance large with respect to $\xi$ and 1 (the dimensionless crack length in the inner domain). For a given V-notch angle, a regularization length and a Poisson's ratio, $A_{eff}$ is obtained by calculating the elastic strain energy difference between a configuration without crack and a configuration with a regularized crack. Note that this method is equivalent to directly computing the contour integral given in Equation (14). FE calculations are performed under small deformation assumptions. Dirichlet boundary conditions corresponding to mode I asymptotic displacement fields around a sharp V-notch are prescribed on the inner domain boundary. The mesh consists of linear 4-nodes elements. It is refined both near the initial sharp crack and in the process zone so that the mesh size is at least 1/10 of the dimensionless regularization length and at least 1/100 of the dimensionless sharp crack length. Differences in the IERR smaller than 1% are obtained for a finer mesh. The IERR requires both the calculation of the elastic strain energy difference and the dimensionless effective crack length. The latter is shown as a function of $\xi$ for several V-notch angles in Figure 3(a). The larger the dimensionless

<figure><img src="image_4.png" /><figcaption>Figure 3: a) Effective crack length $l_{eff}$ and b-c) normalized IERR obtained for b) different V-notch angles and c) different Poisson's ratio, as a function of the dimensionless regularization length $\xi$. The disk symbols in b) and c) represent the normalized IERR value for a sharp crack and the triangle symbols represent the initiation regularization length obtained when applying the CC (see Section 4).</figcaption></figure>

regularization length, the larger the effective crack length as it covers a larger process zone. For a given dimensionless regularization length, the effective length is also larger for smaller V-notch angle. It is worth noting that whatever the V-notch angle, as $\xi \to 0$, the effective crack length tends towards 1, i.e. the dimensionless sharp crack length in the inner domain. Thus, the proposed description of the regularized crack tends towards the description of a sharp crack when the regularization length vanishes. The dimensionless IERR $A_{eff}$ is shown as a function of $\xi$ for different V-notch angles ($v = 0.36$, Figure 3(b)) and for different Poisson's ratio ($\beta = 90$ deg., Figure 3(c)). Overall, it is a decreasing function of the dimensionless regularization length and of the Poisson's ratio. Nevertheless, it may exhibit an increasing ($\xi > 5$, Figure 3(b)) or non-monotonous ($\xi < 5$, Figure 3(b)) variation as a function of the V-notch angle. The FE calculation in the inner domain also provides the variation of the stress along the crack path before initiation. Choosing a normalization of the asymptotic displacement and stress fields so that $\sigma_{\theta\theta}(\theta = 0, \ell = 1) = 1$, the stress field writes, based on the expansion before crack initiation:

$$ \sigma(\ell) = K_I \ell^{\lambda-1}, \qquad (18) $$

The matched asymptotic approach thus enables the calculation of the stress variation along the crack path prior to crack initiation as well as the elastic strain energy change as a function of the dimensionless regularization length. These two functions are used in the sequel to implement the coupled criterion to assess crack initiation.

# 4 The coupled criterion

The CC states that crack initiation occurs if two conditions are fulfilled (Leguillon 2002):

* the stress must be larger than the material strength $\sigma_c$ along the crack path before initiation:
$$ \sigma(y) \ge \sigma_c \quad \forall y < \ell, $$
• the IERR must be larger than the critical ERR: $G_{\text{inc}}(\ell) \ge G_c$.

Applying the coupled criterion consists in coupling both above-mentioned conditions, which results in the following equation system to be solved:

$$ \left\{ \begin{array}{l} \frac{K_1^2}{E} \ell^{2\lambda-1} A_{\text{eff}}(\xi) \ge G_c, \\ K_1 \ell^{\lambda-1} \ge \sigma_c. \end{array} \right. \qquad (19) $$

The two unknowns are the crack length and the GSIF. Solving this system thus provides the initiation crack length $\ell_i$ and the apparent GSIF at initiation $K_{\text{lapp}}$. Combining both equations finally yields the equation that must be solved to determine the dimensionless regularization length $\xi_c$ and the corresponding initiation crack length $\ell_i$:

$$ \ell_i = \frac{\ell_c}{\xi_c} = \frac{E G_c}{A_{\text{eff}}(\xi_c) \sigma_c^2} = \frac{\ell_{\text{mat}}}{A_{\text{eff}}(\xi_c)}. \qquad (20) $$

As recalled in Section 2, the defined relation between $\ell_c$ and $\ell_{\text{mat}}$ through $\ell_c = \eta^2 \ell_{\text{mat}}$ enables determining the equation to be solved to determine the dimensionless regularization length $\xi_c$, and thus the initiation length $\ell_i$. The previous equation rewrites as:

$$ \frac{A_{\text{eff}}(\xi_c)}{\xi_c} = \frac{1}{\eta^2}. \qquad (21) $$

Eq. (21) is solved using Newton's method to determine the dimensionless regularization length and the corresponding initiation crack length. Finally, the apparent GSIF at initiation is obtained as:

$$ K_{\text{lapp}} = \left( \frac{E G_c}{A_{\text{eff}}(\xi_c)} \right)^{1-\lambda} \sigma_c^{2\lambda-1}. \qquad (22) $$

# 5 Results

The CC is applied to study regularized crack initiation at a V-notch. It first requires solving Equation (21) to determine the dimensionless initiation regularization length $\xi_c$. It can be noted that $\xi_c$ only depends on the V-notch angle and Poisson's ratio. It means that for fixed V-notch angle and Poisson's ratio, it is mathematically shown that the initiation length is proportional to the regularization length, which confirms a result previously observed through the comparison between CC and PF simulations (Molnár et al. 2020).

## 5.1 Influence of Poisson's ratio

Figure 4(a) shows $\xi_c$ as a function of the Poisson's ratio for $\beta = 90$ deg. V-notch angle, also represented in Figure 3(c). Increasing the Poisson's ratio induces an increase in the initiation dimensionless regularization length. As a matter of illustration, the process zone obtained for different normalized regularization length are shown in Figure 4(c). Considering or not the regularization, the dimensionless IERR is a decreasing function of the Poisson's ratio (Figure 4). Whatever the Poisson's ratio, the crack regularization induces a decrease in the dimensionless IERR compared to the case of a sharp crack. The decrease in the IERR is more important for large Poisson's ratios. A decrease in the IERR necessarily results in a larger initiation loading so that the IERR reaches the critical ERR. It is illustrated in Figure 5(a) which shows the ratio between the initiation GSIF obtained in the cases of a regularized ($K_{\text{lapp}}^{\text{PF}}$) or a sharp ($K_{\text{lapp}}^{0}$) crack. The increase in initiation GSIF due to the crack regularization is highlighted whatever the Poisson's ratio and it is larger for larger Poisson's ratio. The corresponding initiation length variation as a function of the Poisson's ratio is shown in Figure 5(b). Considering or not the regularization, the initiation length increases with increasing Poisson's ratio. Regularized crack initiation results in a larger initiation length than in the case of sharp crack initiation, the difference between both cases being higher for larger Poisson's ratio.
<figure><img src="image_5.png" /><figcaption>**Figure 4:** Dimensionless a) initiation regularization length $\xi_c$ and b) IERR as a function of the Poisson's ratio $\nu$ for fixed $\beta = 90$ deg. V-notch angle. c) Illustration of the process zone obtained for dimensionless initiation regularization lengths in the case of $\beta = 90$ deg. V-notch angle and (left) $\xi_c = 0.9$, (middle) $\xi_c = 2$ or (right) $\xi_c = 4$. The dimensionless sharp crack length is 1.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5:** Influence of the Poisson's ratio $\nu$ on a) the ratio between the initiation GSIF obtained for a regularized ($K_{\text{lapp}}^{\text{PF}}/K_{\text{lapp}}^{0}$) or a sharp ($K_{\text{lapp}}^{0}$) crack and b) the initiation length to material characteristic length ratio with (triangle symbols) or without (disk symbols) regularization.</figcaption></figure>

## 5.2 Influence of V-notch angle

Figure 6(a) shows $\xi_c$ as a function of the V-notch angle for $\nu = 0.36$ Poisson's ratio, also represented in Figure 3(b). The V-notch angle seems to have a moderate influence on the initiation dimensionless regularization length, which lies between 1.6 and 2.1 for the considered V-notch angles. The resulting process zone has nevertheless a significant influence on the normalized IERR which decreases compared to the case of a sharp crack (Figure 6(b)). The decrease in IERR is more important compared to the sharp crack case as the V-notch angle decreases. As a consequence, the ratio between the initiation GSIF obtained in the cases of a regularized or a sharp crack decreases with increasing V-notch angle and the presence of a process zone results in a larger initiation GSIF than for a sharp crack. Figure 7(b) shows the initiation length as a function of the V-notch angle. Again, the presence of the process zone results in a larger initiation
length compared to sharp crack initiation length, which increases with increasing V-notch angle.

## 5.3 Comparison with PF calculations

We now compare the proposed CC approach for regularized crack initiation to PF calculations, considering a three-point bending specimen studied in Abaza et al. (2022) as a matter of example. The specimen geometry is depicted in Figure 8(a). The specimen dimensions are $L = 8$ mm, $h = 2$ mm and $a_n = 0.4$ mm. The ceramic specimen is made of 3YSZ (Young's modulus: $E = 214$ GPa, Poisson's ratio $\nu = 0.31$, critical energy release rate $G_c = 110$ J/m$^2$ and tensile strength $\sigma_c = 583$ MPa). Four V-notch angles (0, 60, 90 and 120 deg.) are considered. The PF implementation proposed in (Molnár et al. 2022) is used with a spectral energy decomposition (Bernard et al. 2012) and AT1 model. In the PF calculation, the regularization length is adjusted depending on the V-notch angle so as to ensure a Rankine strength surface based on Equation (2). It results in the regularization lengths given in Table 1. The load-displacement curves corresponding to the studied cases obtained with the PF model are linear up to failure that corresponds to a sudden force drop to zero. The apparent GSIF at initiation is then calculated based on the maximum force before failure (Abaza et al. 2022). The initiation GSIF obtained by PF calculation are shown as a function of the V-notch angle in Figure 8(b), which also depicts the initiation GSIF obtained with the proposed approach using the CC, with and without crack regularization. Considering regularization or not, similar trends are obtained using both methods. In case regularization is disregarded, a lower initiation GSIF is obtained, the difference with the regularized crack initiation GSIF decreasing with increasing V-notch angle.

The proposed approach of the CC applied to regularized crack initiation provides less than 5%

<figure><img src="image_7.png" /><figcaption>**Figure 6:** Dimensionless a) initiation regularization length $\xi_c$ and b) IERR as a function of the V-notch angle $\beta$ for fixed $\nu = 0.36$ Poisson's ratio.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7:** Influence of the V-notch angle $\beta$ on a) the ratio between the initiation GSIF obtained for a regularized ($K_{Iapp}$) or a sharp ($K_{Iapp}^0$) crack and b) the initiation length to material characteristic length ratio with (triangle symbols) or without (disk symbols) regularization.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8:** a) Geometry of the three-point bending V-notched specimen and b) initiation GSIF obtained by (square) PF simulations or by application of the CC with (triangle symbols) or without (disk symbols) regularization.</figcaption></figure>

<table><caption>**Table 1:** Regularization lengths for different V-notch angles to ensure a Rankine strength surface in the PF calculations. Initiation length and regularization lengths obtained using the CC applied to regularized crack initiation.</caption><thead><tr><th>V-notch angle β (deg.)</th><th>0</th><th>60</th><th>90</th><th>120</th></tr></thead><tbody><tr><td>PF - Regularization length l<sub>c</sub> (μm)</td><td>26</td><td>28</td><td>34</td><td>36</td></tr><tr><td>CC - initiation length l<sub>i</sub> (μm)</td><td>21.6</td><td>21.4</td><td>22.3</td><td>24.1</td></tr><tr><td>CC - dimensionless regularization length ξ<sub>c</sub></td><td>1.42</td><td>1.63</td><td>1.77</td><td>1.75</td></tr><tr><td>CC - Regularization length ξ<sub>c</sub>l<sub>i</sub> (μm)</td><td>30.6</td><td>34.5</td><td>39.7</td><td>42.3</td></tr></tbody></table>

difference compared to GSIFs obtained using the PF approach (Figure 8(b)). Based on the initiation length and dimensionless regularization length obtained by solving the CC, it also yields an estimate of the initiation regularization length (Table 1). The regularization lengths obtained using the CC are around 20% larger than the one used as inputs in the PF calculations. The remaining differences between the proposed CC approach and PF calculations can be understood by observing the damage field corresponding to the maximum force before initiation in the PF calculations, which is depicted in Figure 9(a). As a complement, the damage field at the first time

<figure><img src="image_10.png" /><figcaption>**Figure 9:** a) Process zone (left) at the force maximum before sudden drop to zero and (right) when damage attains 1 for the first time for different V-notch angles. b) Corresponding variation of damage as a function of the distance to the sharp V-notch tip normalized by the regularization length obtained in PF simulations at the force maximum (dashed line) or when $d = 1$ for the first time (plain line).</figcaption></figure>

for which $d = 1$ is also depicted. It is observed that prior to crack initiation corresponding to the
sudden force drop, a damage zone is established so that the damage variable attains a maximum value around 0.3 to 0.4 depending on the V-notch angle. It can be seen from Figure 9(b) that the process zone extent is around $1.5\ell_c$ along the crack direction at this stage. This is around half the process zone size along the crack direction (i.e., around $3\ell_c$) when damage attains $d = 1$ for the first time. The initial process zone before crack initiation is not considered in the proposed CC approach, which explains the differences observed in terms of GSIF and regularization length. This small damage zone at a V-notch tip was also previously observed by Tanné et al. (2018). It was also discussed in several papers in the case of a crack (Klinsmann et al. 2015; Sargado et al. 2018; Molnár et al. 2020), where it was shown that prescribing phase-field Dirichlet boundary conditions or not along the crack faces resulted in different critical load for its propagation.

# 6 Conclusion

Regularized crack initiation can be described by the matched asymptotic approach of the CC in which a regularized crack is considered to represent a process zone around a sharp crack. While using AT1 PF regularization model in the proposed work, this approach could be generalized to any kind of regularization to describe the process zone around the crack. The larger this process zone, the larger the effective crack length and thus the smaller the normalized IERR, which depends on the Poisson's ratio, the V-notch angle and the regularization length. In case the regularization length is proportional to the material Irwin's length, such as in PF models as shown previously through the parameter $\eta$, it is mathematically shown that the initiation length is proportional to the regularization length. As a consequence, the normalized regularization length depends only on the Poisson's ratio and on the V-notch angle. It is found to be larger than the initiation length for the V-notch configurations under investigation. Considering a regularized crack results in larger initiation GSIF than in the case of a sharp crack without process zone. This difference is larger for decreasing V-notch angles. This result confirms previous findings that in PF calculations in the case of an initial crack, it is essential to provide Dirichlet PF conditions on the initial crack to retrieve Griffith's solution (Molnár et al. 2020; Molnár et al. 2022). The proposed approach of the CC for crack initiation of a regularized crack yields differences smaller than 5% in the initiation GSIF and regularization lengths slightly larger as compared to PF simulations. The remaining slight discrepancy between the proposed approach and PF simulations is likely due to the formation of a process zone prior to crack initiation, locally reducing the stiffness and strength, thus blunting the V-notch and requiring a slightly larger initiation loading. While linking the regularization length to the tensile strength similarly to PF models, these parameters could also be considered as independent so as to implement the CC application in presence of a process zone through a 3-parameter CC formulation. Such a formulation will enable making the regularization parameter independent of the material strength surface and critical ERR. Future works will also cover the formation of the process zone before initiation, which could also be incorporated in the present approach and may be of primary importance for energy driven configurations such as strong singularities or mode III cracking.
## References

Abaza, A., J. Laurencin, A. Nakajo, S. Meille, J. Debayle, and D. Leguillon (2022). "Prediction of crack nucleation and propagation in porous ceramics using the phase-field approach". *Theoretical and Applied Fracture Mechanics* 119, p. 103349. DOI: 10.1016/j.tafmec.2022.103349

Amor, H., J. Marigo, and C. Maurini (2009). "Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments". *Journal of the Mechanics and Physics of Solids* 57.8, pp. 1209-1229. DOI: 10.1016/j.jmps.2009.04.011

Aranda, M. and D. Leguillon (2023). "Prediction of failure of hybrid composites with ultra-thin carbon/epoxy layers using the Coupled Criterion". *Engineering Fracture Mechanics* 281, p. 109053. DOI: 10.1016/j.engfracmech.2023.109053

Bernard, P., N Moes, and N. Chevaugeon (2012). "Damage growth modeling using the Thick Level Set (TLS) approach: Efficient discretization for quasi-static loadings". *Computer Methods in Applied Mechanics and Engineering* 233-236, pp. 11-27. DOI: https://doi.org/10.1016/j.cma.2012.02.020

Borden, M., C. Verhoosel, M. Scott, T. Hughes, and C. Landis (2012). "A phase-field description of dynamic brittle fracture". *Computer Methods in Applied Mechanics and Engineering* 217-220, pp. 77-95. DOI: 10.1016/j.cma.2012.01.008

De Lorenzis, L. and C. Maurini (2022). "Nucleation under multi-axial loading in variational phase-field models of brittle fracture". *International Journal of Fracture* 237.1, pp. 61-81. DOI: 10.1007/s10704-021-00555-6

Doitrand, A., T. Duminy, H. Girard, and X. Chen (2024). "A review of the coupled criterion". DOI: https://hal.science/hal-04023438

Doitrand, A., D. Leguillon, G. Molnár, and V. Lazarus (2023a). "Revisiting facet nucleation under mixed mode I+III loading with T-stress and mode-dependent fracture properties". *International Journal of Fracture* 242, pp. 85-106. DOI: https://doi.org/10.1007/s10704-023-00703-0

Doitrand, A., E. Martin, and D. Leguillon (2020). "Numerical implementation of the coupled criterion: Matched asymptotic and full finite element approaches". *Finite Element in Analysis and Design* 168, p. 103344. DOI: 10.1016/j.finel.2019.103344

Doitrand, A., G. Molnár, R. Estevez, and A. Gravouil (2023b). "Strength-based regularization length in phase field fracture". *Theoretical and Applied Fracture Mechanics* 124, p. 103728. DOI: https://doi.org/10.1016/j.tafmec.2022.103728

Felger, J., P. Rosendahl, D. Leguillon, and W. Becker (2019). "Predicting crack patterns at bi-material junctions: A coupled stress and energy approach". *International Journal of Solids and Structures* 164, pp. 191-201. DOI: https://doi.org/10.1016/j.ijsolstr.2019.01.015

Francfort, G. and J. Marigo (1998). "Revisiting brittle fracture as an energy minimization problem". *Journal of the Mechanics and Physics of Solids* 46.8, pp. 1319-1342. DOI: 10.1016/S0022-5096(98)00034-9

Freddi, F. and G. Royer-Carfagni (2010). "Regularized variational theories of fracture: A unified approach". *Journal of the Mechanics and Physics of Solids* 58.8, pp. 1154-1174. DOI: 10.1016/j.jmps.2010.02.010

García, I. and D. Leguillon (2012). "Mixed-mode crack initiation at a v-notch in presence of an adhesive joint". *International Journal of Solids and Structures* 49.15, pp. 2138-2149. DOI: https://doi.org/10.1016/j.ijsolstr.2012.04.018

Griffith, A. (1921). "The Phenomena of Rupture and Flow in Solids". *Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences* 221.582-593, pp. 163-198

Griffith, A. (1924). "The theory of rupture". *First Int. Cong. Appl. Mech*, pp. 55-63

Jiménez-Alfaro, S. and D. Leguillon (2022). "Modelling of glass matrix composites by the Coupled Criterion and the Matched Asymptotics approach. The role of a single platelet". *Theoretical and Applied Fracture Mechanics* 122, p. 103650. DOI: https://doi.org/10.1016/j.tafmec.2022.103650

Klinsmann, M., D. Rosato, M. Kamlah, and M. RM (2015). "An assessment of the phase field formulation for crack growth". *Computer Methods in Applied Mechanics and Engineering* 294, pp. 313-330. DOI: https://doi.org/10.1016/j.cma.2015.06.009
Kumar, A., B. Bourdin, G. 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, p. 104027. DOI: 10.1016/j.jmps.2020.104027

Labossiere, P. E. and M. L. Dunn (1999). "Stress intensities at interface corners in anisotropic bimaterials". *Engineering Fracture Mechanics* 62.6, pp. 555-576. DOI: https://doi.org/10.1016/S0013-7944(99)00005-3

Leguillon, D. (2002). "Strength or toughness? A criterion for crack onset at a notch". *European Journal of Mechanics - A/Solids* 21(1), pp. 61-72. DOI: 10.1016/S0997-7538(01)01184-6

Leguillon, D., C. Lacroix, and E. Martin (2000). "Interface debonding ahead of a primary crack". *Journal of the Mechanics and Physics of Solids* 48.10, pp. 2137-2161. DOI: 10.1016/S0022-5096(99)00101-5

Leguillon, D. and S. Murer (2012). "Fatigue crack nucleation at a stress concentration point". *CP 2012 Conference Proceedings*. Vol. 46. https://www.gruppofrattura.it/ocs/index.php/esis/CP2012/paper/viewFile/9235/5996

Leguillon, D. and E. Sanchez-Palencia (1987). *Computation of Singular Solutions in Elliptic Problems and Elasticity*. Wiley, USA

Leguillon, D. and Z. Yosibash (2017). "Failure initiation at V-notch tips in quasi-brittle materials". *International Journal of Solids and Structures* 122-123, 1-13. DOI: https://doi.org/10.1016/j.ijsolstr.2017.05.036

Mesgarnejad, A., B. Bourdin, and M. Khonsari (2015). "Validation simulations for the variational approach to fracture". *Computer Methods in Applied Mechanics and Engineering* 290, pp. 420-437. DOI: https://doi.org/10.1016/j.cma.2014.10.052

Miehe, C., M. Hofacker, and F. Welschinger (2010b). "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, pp. 2765-2778. DOI: 10.1016/j.cma.2010.04.011

Mittelman, B. and Z. Yosibash (2014). "Asymptotic analysis of the potential energy difference because of a crack at a V-notch edge in a 3D domain". *Engineering Fracture Mechanics* 131, 232-256. DOI: http://dx.doi.org/10.1016/j.engfracmech.2014.07.031

Mittelman, B. and Z. Yosibash (2015). "Energy release rate cannot predict crack initiation orientation in domains with a sharp V-notch under mode III loading". *Engineering Fracture Mechanics* 141, pp. 230-241. DOI: 10.1016/j.engfracmech.2015.05.008

Molnár, G., A. Doitrand, R. Estevez, and A. Gravouil (2020). "Toughness or strength? Regularization in phase-field fracture explained by the coupled criterion". *Theoretical and Applied Fracture Mechanics* 109, p. 102736. DOI: 10.1016/j.tafmec.2020.102736

Molnár, G., A. Doitrand, A. Jaccon, B. Prabel, and A. Gravouil (2022). "Thermodynamically consistent linear-gradient damage model in Abaqus". *Engineering Fracture Mechanics* 266, p. 108390. DOI: 10.1016/j.engfracmech.2022.108390

Molnár, G., A. Doitrand, and V. Lazarus (2024). "Phase-field simulation and coupled criterion link echelon cracks to internal length in antiplane shear". *Journal of the Mechanics and Physics of Solids*, p. 105675. DOI: https://doi.org/10.1016/j.jmps.2024.105675

Pham, K. and K. Ravi-Chandar (2016). "On the growth of cracks under mixed-mode I + III loading". *International Journal of Fracture* 199, 105-134. DOI: 10.1007/s10704-016-0098-6

Pham, K., K. Ravi-Chandar, and C. Landis (2017). "Experimental validation of a phase-field model for fracture". *International Journal of Fracture* 205.1, pp. 83-101. DOI: 10.1007/s10704-017-0185-3

Reinoso, J., A. Arteiro, M. Paggi, and P. Camanho (2017). "Strength prediction of notched thin ply laminates using finite fracture mechanics and the phase field approach". *Composites Science and Technology* 150, pp. 205-216. DOI: 10.1016/j.compscitech.2017.07.020

Sargado, J., E. Keilegavlen, I. Berre, and J. Nordbotten (2018). "High-accuracy phase-field models for brittle fracture based on a new family of degradation functions". *Journal of the Mechanics and Physics of Solids* 111, pp. 458-489. DOI: https://doi.org/10.1016/j.jmps.2017.10.015

Strobl, M., P. Dowgiallo, and T. Seelig (2017). "Analysis of Hertzian indentation fracture in the framework of finite fracture mechanics". *International Journal of Fracture* 206, 67-79. DOI: 10.1007/s10704-017-0201-7

Strobl, M. and T. Seelig (2020). "Phase field modeling of Hertzian indentation fracture". *Journal of the Mechanics and Physics of Solids* 143, p. 104026. DOI: https://doi.org/10.1016/j.jmps.2020.104026
Tanné, E., T. Li, B. Bourdin, 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**, pp. 80–99. doi: 10.1016/j.jmps.2017.09.006

Vicentini, F., C. Zolesi, P. Carrara, C. Maurini, and L. De Lorenzis (2024). "On the energy decomposition in variational phase-field models for brittle fracture under multi-axial stress states". *International Journal of Fracture*. doi: https://doi.org/10.1007/s10704-024-00763-w

Weißgraeber, P., D. Leguillon, and W. Becker (2016). "A review of Finite Fracture Mechanics: crack initiation at singular and non-singular stress raisers". *Archive Appl. Mech.* **86**(1-2), 375–401. doi: 10.1007/s00419-015-1091-7

Yosibash, Z., E. Priel, and D Leguillon (2006). "A failure criterion for brittle elastic materials under mixed-mode loading". *International Journal of Fracture* **141**, 291–312. doi: https://doi.org/10.1007/s10704-006-0083-6

Zghal, J., K. Moreau, N. Moës, D. Leguillon, and C. Stolz (2018). "Analysis of the failure at notches and cavities in quasi-brittle media using the Thick Level Set damage model and comparison with the coupled criterion". *International Journal of Fracture* **211**, pp. 253–280. doi: https://doi.org/10.1007/s10704-018-0287-6

**Authors' contributions** AD performed numerical simulations and drafted the manuscript. GM and AD developed the methodology, conceived of the study, and participated in its design, coordination, and critical review of the manuscript. All authors read and approved the final manuscript.