# A review of the coupled criterion

Aurélien Doitrand, Thomas Duminy, Hugo Girard, Xi Chen

► To cite this version:

Aurélien Doitrand, Thomas Duminy, Hugo Girard, Xi Chen. A review of the coupled criterion. 2023. [hal-04023438](https://hal.science/hal-04023438)

HAL Id: hal-04023438

https://hal.science/hal-04023438

Preprint submitted on 10 Mar 2023

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
# A review of the coupled criterion

Aurélien Doitrand, Thomas Duminy, Hugo Girard, Xi Chen

► To cite this version:

Aurélien Doitrand, Thomas Duminy, Hugo Girard, Xi Chen. A review of the coupled criterion. Journal of Theoretical, Computational and Applied Mechanics, 2024, 10.46298/jtcam.11072. hal-04023438v2

HAL Id: hal-04023438

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

Submitted on 21 Nov 2024

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
# A review of the coupled criterion

Aurélien Doitrand¹, Thomas Duminy¹, Hugo Girard¹,², and Xi Chen¹

¹ Univ Lyon, INSA Lyon, Université Claude Bernard Lyon 1, CNRS UMR5510, MATEIS, F-69621 Villeurbanne, France, aurelien.doitrand@insa-lyon.fr

² Axel'One (Solvay), 87 Avenue des Frères Perret, CS 70061, 69192 Saint Fons, France

We provide a review of recent works about the Coupled Criterion (CC), originally developed to assess the initiation of a crack and thus overcome a main limitation of linear elastic fracture mechanics. This paper includes theoretical and implementation aspects of the CC, new fundamental developments such as extensions to 3D, dynamic aspects, nonlinearities or fatigue. It also focuses on confrontations between the CC and other fracture modeling approaches such as cohesive zone models or phase field approach for fracture. Originally developed in 2D, for linear elastic materials under the assumption of small deformations, the CC has now expanded and can be applied to a wide range of applications from macroscopic to atomic scale.

Keywords Coupled criterion ;

## 1 Introduction

Crack growth can generally be studied based on Linear Elastic Fracture Mechanics (LEFM) using Griffith's theory *Griffith* (1921); *Griffith* (1924) or Irwin's approach *Irwin* (1958). However, both theories are limited to study the propagation of an existing crack and cannot predict crack initiation.

The framework of Finite Fracture Mechanics (FFM) was thus developed to fill this gap in fracture mechanics by addressing crack initiation based on finite size cracks (Hashin 1996; Zhou et al. 1999; Nairn 2000). Since then, several approaches emerged in the FFM framework such as the Theory of Critical Distances (TCD) (Taylor et al. 2005). The FFM framework was enriched by Leguillon (2002) who proposed the coupled criterion (CC) for crack initiation in brittle materials. This criterion states that crack initiation can occur provided an energy-based criterion and a stress-based criterion are simultaneously fulfilled. This criterion thus basically requires as inputs the elastic stiffness tensor and two fracture parameters: the material critical energy release rate and its tensile strength.

So far, the CC has proved its efficiency to predict crack initiation in many situations including at weak stress singularities such as V-notches (Leguillon 2002; Leguillon et al. 2003b; Picard et al. 2006; Cornetti et al. 2006; Yosibash et al. 2006; Hebel et al. 2008), at adhesive joints (Leguillon et al. 2003a; Müller et al. 2006; Moradi et al. 2013; Ernesto Mendoza-Navarro et al. 2013; Hell et al. 2014) or at interfaces (Martin et al. 2004; Leguillon et al. 2013a; Leguillon et al. 2013b), in composites (Hebel et al. 2008; Martin et al. 2010; García et al. 2014) or considering thermal residual stresses (Henninger et al. 2007a; Leguillon et al. 2015b).

All these applications are described in the review paper by Weißgraeber et al. (2016). In this paper, several possible extensions of the CC were pointed out, such as considering fatigue loadings, dynamic loadings, or extension to the 3D case. The aim of the present paper is to give an overview of the recent developments and applications related to the CC. It starts with recalling the CC theory in Section 2, then describes the new fundamental developments of the CC in Section 3, provides comparisons with other fracture approaches in Section 4 and focuses on some CC applications in Section 5.

## 2 Theory

The objective of the CC is to assess the initiation of a crack, i.e. to determine at which loading magnitude the crack initiates and what is the initiation crack shape. The CC is founded on
**Identifiers**

DOI 10.46298/jtcam.11072

HAL hal-04023438v2

**History**

Received Mar 10, 2023

Accepted Aug 23, 2024

Published Nov 24, 2024

**Associate Editor**

Julien RÉTHORÉ

**Reviewers**

Alberto SAPORA

Anonymous

**Open Review**

HAL hal-04770360

**Licence**

CC BY 4.0

©The Authors

# A review of the coupled criterion

Aurélien DOITRAND<fnref n="1" />, Thomas DUMINY<fnref n="1" />, Hugo GIRARD<fnref n="1" />,<fnref n="2" />, and Xi CHEN<fnref n="1" />

We provide a review of recent works about the Coupled Criterion (CC), originally developed to assess the initiation of a crack and thus overcome a main limitation of linear elastic fracture mechanics. This paper includes theoretical and implementation aspects of the CC, new fundamental developments such as extensions to 3D, dynamic aspects, nonlinearities or fatigue. It also focuses on confrontations between the CC and other fracture modeling approaches such as cohesive zone models or phase field approach for fracture. Originally developed in 2D, for linear elastic materials under the assumption of small deformations, the CC has now expanded and can be applied to a wide range of applications from macroscopic to atomic scales.

**Keywords:** Coupled criterion; 3D formulation; nonlinearities; dynamic scenarios; fatigue; comparison with other fracture models

## 1 Introduction

Crack growth can generally be studied based on Linear Elastic Fracture Mechanics (LEFM) using Griffith's theory (1921), or Irwin's approach (1958). However, both theories are limited to study the propagation of an existing crack and cannot predict crack initiation.

The framework of Finite Fracture Mechanics (FFM) was thus developed to fill this gap in fracture mechanics by addressing crack initiation based on finite size cracks (Hashin 1996; Zhou et al. 1999; Nairn 2000). Since then, several approaches emerged in the FFM framework such as the Theory of Critical Distances (TCD) in which the input parameters are a characteristic length and a strength (Taylor et al. 2005). The FFM framework was enriched by Leguillon (2002) who proposed the coupled criterion (CC) for crack initiation in brittle materials. This criterion states that crack initiation can occur provided an energy-based criterion and a stress-based criterion are simultaneously fulfilled. This criterion thus basically requires as inputs the elastic stiffness tensor and two fracture parameters: the material critical energy release rate and its tensile strength.

So far, the CC has proved its efficiency to predict crack initiation in many situations including at weak stress singularities such as V-notches (Leguillon 2002; Leguillon and Yosibash 2003; Picard et al. 2006; Cornetti et al. 2006; Yosibash et al. 2006; Hebel and Becker 2008), at adhesive joints (Legaillon et al. 2003; Müller et al. 2006; Moradi et al. 2013; Ernesto Mendoza-Navarro et al. 2013; Hell et al. 2014) or at interfaces (Martin and Leguillon 2004; Leguillon and Martin 2013a; Leguillon and Martin 2013b), in composites (Hebel and Becker 2008; Martin et al. 2010; García et al. 2014) or considering thermal residual stresses (Henninger and Leguillon 2007; Leguillon, Martin, Ševeček, et al. 2015).

All these applications are described in the review paper by Weißgraeber, Leguillon, et al. (2016). In this paper, several possible extensions of the CC were pointed out, such as considering fatigue loadings, dynamic loadings, or extension to the 3D case. The aim of the present paper is to give an overview of the recent developments and applications related to the CC. It starts with recalling the CC theory in Section 2, then describes the new fundamental developments of the CC in Section 3, provides comparisons with other fracture approaches in Section 4 and focuses on some CC applications in Section 5.

<footnote n="1">¹ Univ Lyon, INSA Lyon, Université Claude Bernard Lyon 1, CNRS UMR5510, MATEIS, F-69621 Villeurbanne, France</footnote>

<footnote n="2">² Axel'One (Solvay), 87 Avenue des Frères Perret, CS 70061, 69192 Saint Fons, France</footnote>
the main idea that crack initiation occurs if two conditions are fulfilled. These conditions are described in the following.

## 2.1 Stress condition

The stress condition of the CC, as originally proposed by Leguillon (2002), states that before initiation, the area corresponding to a newly formed crack (described by a surface in 3D or curve in 2D, denoted $\Gamma$, whose size is noted $S$ (area, in 3D) or $l$ (length, in 2D)) must undergo a sufficiently large stress level. From a general perspective, it writes:

$$f(\underline{\sigma}(\mathbf{x}, U)) \geq 0 \forall \mathbf{x} \in \Gamma, \qquad (1)$$

where $f$ describes the material strength surface, $\mathbf{x}$ the position and $U$ the applied displacement or loading. The material strength surface is defined as the surface in the principal stress space determined by the critical failure stress under monotonically increasing uniform stress state. For crack initiation in homogeneous isotropic materials following a Rankine-like strength surface (Figure 1), the stress condition reverts to comparing the stress normal to the crack plane ($\sigma_{nn}$) to the material tensile strength ($\sigma_c$) (Laguillon 2002; Rosendahl et al. 2017; Doitrand et al. 2020f):

$$\sigma_{nn}(\mathbf{x}, U) \geq \sigma_c \forall \mathbf{x} \in \Gamma. \qquad (2)$$

It may also be written using a principal stress criterion (Weissgraeber et al. 2013; Stein et al. 2015) involving the shear stress ($\sigma_{nt}$):

$$\frac{\sigma_{nn}(\mathbf{x}, U)}{2} + \sqrt{\left(\frac{\sigma_{nn}(\mathbf{x}, U)}{2}\right)^2 + \sigma_{nt}(\mathbf{x}, U)^2} \geq \sigma_c \forall \mathbf{x} \in \Gamma. \qquad (3)$$

Other criteria involving the material shear strength ($\tau_c$) can be used, such as power ellipse

<figure><img src="image_2.png" /><figcaption>Figure 1: a) Rankine-like, b) power-ellipse ($q = 2, 3, 5$ or $10$) or c) maximum stress based strength surfaces.</figcaption></figure>

criterion (Figure 1(b), generally used with $q=2$) (Tran et al. 2012; Modniks et al. 2015; Carrère et al. 2015; Muñoz-Reja et al. 2016; Felger et al. 2017b; Felger et al. 2019a; Doitrand et al. 2020b; Le Pavic et al. 2020):

$$\left(\frac{\sigma_{nn}(\mathbf{x}, U)^q}{\sigma_c} + \frac{\sigma_{nt}(\mathbf{x}, U)^q}{\tau_c}\right)^{1/q} \geq 1 \forall \mathbf{x} \in \Gamma, \qquad (4)$$

or maximum stress criterion (Figure 1(c)) (Doitrand et al. 2018b; Doitrand et al. 2020b):

$$\max \left( \frac{\sigma_{nn}(\mathbf{x}, U)}{\sigma_c}, \frac{|\sigma_{nt}(\mathbf{x}, U)|}{\tau_c} \right) \geq 1 \forall \mathbf{x} \in \Gamma. \qquad (5)$$

These criteria are usually employed for configurations in which the stress state over the presupposed crack path includes both tensile and shear components, such as for interface cracking or crack initiation under mixed mode loading (Laguillon et al. 2008b; Doitrand et al. 2023a).

The two unknowns in the stress criterion are the prescribed displacement or loading and the crack surface. It allows determining, for a given loading, the set of cracks $\Gamma_{\sigma}(U)$ for which the stress criterion is satisfied.

$$\Gamma_{\sigma}(U) = \{\Gamma | f(\underline{\sigma}(\mathbf{x}, U)) \geq 0 \forall \mathbf{x} \in \Gamma\}. \qquad (6)$$
## 2 Coupled criterion theory

The objective of the CC is to assess the initiation of a crack, i.e. to determine at which loading magnitude the crack initiates and what is the initiation crack shape. The CC is founded on the main idea that crack initiation occurs if two conditions are fulfilled. The stress and the energy conditions are described in the following.

### 2.1 Stress condition

The stress condition of the CC, as originally proposed by Leguillon (2002), states that before initiation, the area corresponding to a newly formed crack (described by a surface in 3D or curve in 2D, denoted $\Gamma$, whose size is noted $S$ (area, in 3D) or $l$ (length, in 2D)) must undergo a sufficiently large stress level. From a general perspective, it reads

$$f(\underline{\sigma}(\mathbf{x}, U)) \geq 0, \quad \forall \mathbf{x} \in \Gamma, \qquad (1)$$

where $f$ is a functional so that $f = 0$ describes the material strength surface, $\mathbf{x}$ the position and $U$ the applied displacement or loading. The material strength surface is defined as the surface in the principal stress space determined by the critical failure stress under monotonically increasing uniform stress state. For crack initiation in homogeneous isotropic materials following a Rankine-like strength surface, see Figure 1, the stress condition reverts to comparing the stress normal to the crack plane $\sigma_{nn}$ to the material tensile strength $\sigma_c$ (Laguillon 2002; Rosendahl et al. 2017; Doitrand, Martin, et al. 2020):

$$\sigma_{nn}(\mathbf{x}, U) \geq \sigma_c, \quad \forall \mathbf{x} \in \Gamma. \qquad (2)$$

It may also be written using a principal stress criterion (Weißgraeb er and Becker 2013; Stein et al. 2015) involving the shear stress $\sigma_{nt}$:

$$\frac{\sigma_{nn}(\mathbf{x}, U)}{2} + \sqrt{\left(\frac{\sigma_{nn}(\mathbf{x}, U)}{2}\right)^2 + \sigma_{nt}(\mathbf{x}, U)^2} \geq \sigma_c, \quad \forall \mathbf{x} \in \Gamma. \qquad (3)$$

Other criteria involving the material shear strength $\tau_c$ can be used, such as power ellipse criterion

<figure><img src="image_2.png" /><figcaption>**Figure 1** Examples of material strength surfaces obtained under an homogeneous stress state: (a) Rankine-like, (b) power-ellipse ($q = 2, 3, 5$ or 10) or (c) maximum stress based strength surfaces.</figcaption></figure>

as in Figure 1(b), generally used with $q = 2$ (Tran et al. 2012; Modniks et al. 2015; Carrère et al. 2015; Muñoz-Reja et al. 2016; Felger et al. 2017b; Felger, Rosendahl, et al. 2019; Doitrand, Henry, Saad, et al. 2020; Le Pavic et al. 2020):

$$\left(\frac{\sigma_{nn}(\mathbf{x}, U)^q}{\sigma_c} + \frac{\sigma_{nt}(\mathbf{x}, U)^q}{\tau_c}\right)^{1/q} \geq 1, \quad \forall \mathbf{x} \in \Gamma, \qquad (4)$$

or maximum stress criterion as in Figure 1(c) (Doitrand and Leguillon 2018b; Doitrand, Henry, Saad, et al. 2020):

$$\max\left(\frac{\sigma_{nn}(\mathbf{x}, U)}{\sigma_c}, \frac{|\sigma_{nt}(\mathbf{x}, U)|}{\tau_c}\right) \geq 1, \quad \forall \mathbf{x} \in \Gamma. \qquad (5)$$
If crack initiation occurs near a stress singularity or a stress concentrator, the stress criterion is generally a decreasing function of the distance to the singular point or stress concentrator and $\Gamma_{\sigma}(U)$ corresponds to an upper bound of the admissible crack surfaces (Figure 2(a)). Some exceptions may be encountered in which the stress variation is not a monotonically decreasing function of the distance (Figure 2(b)). This situation may for instance be encountered for scarf joints (Doitrand et al. 2018b; Felger et al. 2019a). As a consequence, $\Gamma_{\sigma}(U)$ is the union of all areas for which the stress criterion is fulfilled (filled zone in Figure 2(b)). Other configurations where the stress variation exhibits a maximum can also be encountered, for instance under bending of a specimen having compressive residual stresses at the surface (Figure 2(c)). The initiation crack

<figure><img src="image_3.png" /><figcaption>Figure 2: Stress variation as a function of the crack surface, a) monotonically decreasing variation or b-c) including a local minimum or a local maximum. The filled zones correspond to the crack surfaces for which the stress criterion is fulfilled.</figcaption></figure>

surface is finally determined by the energy criterion (see Section 2.2) as the one minimizing the initiation imposed loading. In practice, the crack surface may emanate from the stress singularity or concentrator (Nguyen et al. 2012) but it was also shown that it is possible that crack initiates only over the hatched zone in Figure 2(b) for instance in scarf joint configurations (Doitrand et al. 2018b). It is worth noting that the stress condition as originally proposed by Leguillon (2002) is not a point stress criterion but it is a non-local condition which is more restrictive since it has to be fulfilled everywhere over $\Gamma$ before initiation. An alternative approach consists in replacing the stress by the average stress over the crack surface (Cornetti et al. 2006). Contrary to the stress criterion (Equation (1)), which is a non-local condition, the average stress criterion reverts to a point-stress condition. As a consequence, the criterion becomes less restrictive than the one originally proposed by Leguillon since it allows a stress level smaller than the critical one over a part of the area where the crack nucleates. For instance, it does not require any particular treatment in the above-mentioned configurations of non-monotonically decreasing stress. Several works compared the use of original stress or average stress criterion, generally resulting in failure load differences smaller than 10% (Sapora et al. 2018b; Cornetti et al. 2019b; Doitrand et al. 2021a). Larger differences are obtained when decreasing the specimen size with respect to the material characteristic length (Baldassari et al. 2023). Similar adhesive joint failure loads were also obtained with both approaches in case the stress criterion shows a monotonic behavior, whereas larger differences were obtained for non-monotonic stress variations (Talmon l'Armée et al. 2020). An alternative to the stress criterion was proposed by Rosendahl et al. (2019) who used a stretch criterion in the case of silicone and Li et al. (2019) who used a strain-based criterion to determine the potential crack initiation zones and further apply the energy criterion in these zones. Further investigations could be established to compare stress or strain-based criterion.

## 2.2 Energy condition

Crack propagation can be studied in the framework of LEFM based on Griffith's approach (Griffith 1921) using the crack propagation condition:

$$ \mathcal{G} \geq \mathcal{G}_c, \qquad (7) $$

where $\mathcal{G}$ is the Energy Release Rate (ERR) and $\mathcal{G}_c$ the critical ERR. However, this criterion is inappropriate to assess the nucleation of a crack since the ERR at a singular point or stress
These criteria are usually employed for configurations in which the stress state over the presupposed crack path includes both tensile and shear components, such as for interface cracking or crack initiation under mixed mode loading (Leguillon and Murer 2008; Doitrand, Leguillon, et al. 2023).

The two unknowns in the stress criterion are the prescribed displacement or loading at crack initiation and initiation the crack surface. It allows determining, for a given loading, the set of cracks $\Gamma_{\sigma}(U)$ for which the stress criterion is satisfied:

$$ \Gamma_{\sigma}(U) = \{\Gamma \mid f(\underline{\sigma}(x, U)) \ge 0, \forall x \in \Gamma\}. \qquad (6) $$

If crack initiation occurs near a stress singularity or a stress concentrator, the stress criterion is generally a decreasing function of the distance to the singular point or stress concentrator. In case there is a unique crack configuration corresponding to a given crack surface, $\Gamma_{\sigma}(U)$ corresponds to an upper bound of the admissible crack surfaces, Figure 2(a). Some exceptions may be encountered in which the stress variation is not a monotonically decreasing function of the distance, Figure 2(b). This situation may for instance be encountered for scarf joints (Doitrand and Leguillon 2018b; Felger, Rosendahl, et al. 2019). As a consequence, $\Gamma_{\sigma}(U)$ is the union of all areas for which the stress criterion is fulfilled (filled zone in Figure 2(b)). Other configurations where the stress variation exhibits a maximum can also be encountered, for instance under bending of a specimen having compressive residual stresses at the surface, see Figure 2(c).

<figure><img src="image_3.png" /><figcaption>**Figure 2** Stress variation as a function of the crack surface: (a) monotonically decreasing variation or (b)-(c) including a local minimum or a local maximum. The filled zones correspond to the crack surfaces for which the stress criterion is fulfilled.</figcaption></figure>

The initiation crack surface is finally determined by the energy criterion (see Section 2.2) as the one minimizing the initiation imposed loading. In practice, the crack surface may emanate from the stress singularity or concentrator (L.-M. Nguyen et al. 2012) but it was also shown that it is possible that crack initiates only over the hatched zone in Figure 2(b) for instance in scarf joint configurations (Doitrand and Leguillon 2018b). It is worth noting that the stress condition as originally proposed by Leguillon (2002) is not a point stress criterion but it is a non-local condition which is more restrictive since it has to be fulfilled everywhere over $\Gamma$ before initiation. An alternative approach consists in replacing the stress by the average stress over the crack surface (Cornetti et al. 2006). In contrast to the stress criterion in Equation (1), which is a non-local condition, the average stress criterion reverts to a point-stress condition. As a consequence, the criterion becomes less restrictive than the one originally proposed by Leguillon (2002) since it allows a stress level smaller than the critical one over a part of the area where the crack nucleates. For instance, it does not require any particular treatment in the above-mentioned configurations of non-monotonically decreasing stress. Several works compared the use of original stress or average stress criterion, generally resulting in failure load differences smaller than 10% (Sapora et al. 2018; Cornetti and Sapora 2019; Doitrand, Cornetti, et al. 2021). Larger differences are obtained when decreasing the specimen size with respect to the material characteristic length (Baldassari et al. 2023). Similar adhesive joint failure loads were also obtained with both approaches in case the stress criterion shows a monotonic behavior, whereas larger differences were obtained for non-monotonic stress variations (Talmon l'Armée and Becker 2020). An alternative to the stress criterion was proposed by Rosendahl et al. (2019) who used a stretch criterion in the case of silicone and Li et al. (2019) who used a strain-based criterion to determine the potential crack
concentrator generally tends towards 0 when the crack surface tends towards 0. This result implies that the crack must initiate over a finite rather than an infinitesimal surface. The energy condition of the coupled criterion derives from the energy conservation principle before and after crack nucleation over a finite surface. It consists in a balance of the variation of external work forces ($\Delta W_{ext}$), kinetic energy ($\Delta W_k$), elastic strain energy ($\Delta W_{el}$), crack surface creation energy ($G_cS$, where S the crack surface) and energy dissipated into other mechanisms ($\Delta W_d$) such as plasticity, diffused damage or friction for instance:

$$ \Delta W_k + \Delta W_{el} + \Delta W_d + G_cS = \Delta W_{ext}. \qquad (8) $$

This energy balance rewrites in a form similar to Griffith's criterion in which the ERR is replaced by the Incremental Energy Release Rate (IERR) $G_{inc}$:

$$ G_{inc}(S, U) = \frac{\Delta W_{ext} - \Delta W_k - \Delta W_{el} - \Delta W_d}{S} = G_c. \qquad (9) $$

It allows determining, for a given loading, the set of cracks $\Gamma_G$ for which the criterion is satisfied.

$$ \Gamma_G(U) = \{\Gamma | G_{inc}(S, U) \geq G_c\}. \qquad (10) $$

The IERR and ERR can be obtained from each other based on the following expressions:

$$ \left\{ \begin{array}{l} G_{inc}(S) = \frac{1}{S} \int_{\Gamma} G(s) ds, \\ G(S) = G_{inc}(S) + S \frac{dG_{inc}}{dS}. \end{array} \right. \qquad (11) $$

It is worth highlighting that even if the CC was originally developed to assess crack nucleation (Leguillon 2002), it can also be used to study crack propagation. Indeed, as the crack surface tends towards 0, the IERR tends towards the ERR as $\lim_{S \to 0} G_{inc}(S) = G(S)$. The energy condition given in Equation (9) actually reverts to Griffith's criterion when the crack surface tends towards 0, therefore LEFM is actually included in the CC. As a consequence, in the case of an initial crack, its propagation can be studied using either Equation (9) or Equation (7), which induces differences on the critical propagation loading smaller than 10% whatever the material properties (Molnár et al. 2020a). For singularities other than a crack or stress concentrators, Equation (9) must be used since the ERR tends towards 0 when the crack length tends towards 0. Unlike for the stress condition, no alternative energy criterion are proposed since Equation (9) results from an indisputable energy conservation principle (Equation (8)).

In some works, $G_c$ is not considered as a constant but may vary. For instance, Catalanotti et al. (2013) proposed to consider increasing $G_c$ as a function of the crack length during crack propagation in composite laminates that follows a R-curve obtained experimentally. Configurations where $G_c$ depends on the crack velocity may be also be encountered (Kanninen et al. 1985; Dally et al. 1985; Molnár et al. 2020), which should be possibly implemented if dynamic aspects are considered (Doitrand et al. 2022). Another example consists in considering mode-dependent incremental and critical ERR in the energy criterion. Two approaches are classically used. The first one consists in separately considering opening and shear incremental to critical ERR contributions, which requires the calculation of both contributions and introduces the opening ($G_{Ic}$) and shear ($G_{IIc}$) critical ERR (Tran et al. 2012; Modniks et al. 2015):

$$ \frac{G_{inc}^{I}(S)}{G_{Ic}} + \frac{G_{inc}^{II}(S)}{G_{IIc}} = 1. \qquad (12) $$

The second approach consists in introducing the dependency of the critical ERR to the mode mixity $\psi$ and comparing the IERR to the average critical ERR over the crack path $\overline{G}_c(S)$ (Mantič 2009; García et al. 2012; Carrère et al. 2015):

$$ G_{inc}(S) \geq \overline{G}_c(S), \qquad (13) $$

where $\overline{G}_c$ is defined as:

$$ \overline{G}_c(S) = \frac{1}{S} \int_{\Gamma} G_c(\psi(x)) dS. \qquad (14) $$
initiation zones and further apply the energy criterion in these zones. Further investigations could be established to compare stress or strain-based criterion.

## 2.2 Energy condition

Crack propagation can be studied in the framework of LEFM based on Griffith's approach (Griffith 1921; Griffith 1924) using the crack propagation condition

$$ \mathcal{G} \geq \mathcal{G}_c \qquad (7) $$

where $\mathcal{G}$ is the Energy Release Rate (ERR) and $\mathcal{G}_c$ the critical ERR. However, this criterion is inappropriate to assess the nucleation of a crack since the ERR at a weak singular point or stress concentrator generally tends towards 0 when the crack surface tends towards 0. This result implies that the crack must initiate over a finite rather than an infinitesimal surface. The energy condition of the coupled criterion derives from the energy conservation principle between the states before and after crack nucleation over a finite surface. It consists in a balance of the variation of external work forces ($\Delta W_{ext}$), kinetic energy ($\Delta W_k$), elastic strain energy ($\Delta W_{el}$), crack surface creation energy ($\Delta W_{crack} = \int_{\Gamma} \mathcal{G}_c(\mathbf{x}) \, dS$) and energy dissipated into other mechanisms ($\Delta W_d$) such as plasticity, diffused damage or friction for instance

$$ \Delta W_k + \Delta W_{el} + \Delta W_d + \Delta W_{crack} = \Delta W_{ext}. \qquad (8) $$

This energy balance can be rewritten by introducing the Incremental Energy Release Rate (IERR) $\mathcal{G}_{inc}$ which depends on the surface $S$ of the crack $\Gamma$:

$$ \mathcal{G}_{inc}(\Gamma, U) = \frac{\Delta W_{ext} - \Delta W_k - \Delta W_{el} - \Delta W_d}{S} = \bar{\mathcal{G}}_c(\Gamma) \quad \text{where} \quad \bar{\mathcal{G}}_c(\Gamma) = \frac{1}{S} \int_{\Gamma} \mathcal{G}_c(\mathbf{x}) \, dS. \qquad (9) $$

It allows determining, for a given loading, the set of cracks $\Gamma_G$ for which the criterion is satisfied:

$$ \Gamma_G(U) = \{\Gamma | \mathcal{G}_{inc}(\Gamma, U)) \geq \bar{\mathcal{G}}_c(\Gamma)\}. \qquad (10) $$

For the sake of simplicity, the dependency of $\mathcal{G}_{inc}$ to the crack surface $S$ will be adopted in the sequel. Note that this notation is valid as long as there is a unique crack configuration corresponding to a given crack surface. Otherwise, the notation $\mathcal{G}_{inc}(\Gamma)$ should be preferred as it allows comparing two different crack configurations corresponding to the same crack surface $S$. The IERR and ERR can be obtained from each other based on the expressions

$$ \mathcal{G}_{inc}(S) = \frac{1}{S} \int_{\Gamma} \mathcal{G}(s) \, ds \quad \text{and} \quad \mathcal{G}(S) = \mathcal{G}_{inc}(S) + S \frac{d\mathcal{G}_{inc}}{dS}. \qquad (11) $$

It is worth highlighting that even if the CC was originally developed to assess crack nucleation (Leguillon 2002), it can also be used to study crack propagation. Indeed, as the crack surface tends towards 0, the IERR tends towards the ERR as

$$ \lim_{S \to 0} \mathcal{G}_{inc}(S) = \mathcal{G}(S). \qquad (12) $$

The energy condition given in Equation (9) thus reverts to Griffith's criterion when the crack surface tends towards 0, therefore LEFM is actually included in the CC. As a consequence, in the case of an initial crack, its propagation can be studied using either Equation (9) or Equation (7), which induces differences on the critical propagation loading smaller than 10% whatever the material properties (Molnár, Doitrand, et al. 2020). For singularities other than a crack or stress concentrators, Equation (9) must be used since the ERR tends towards 0 when the crack length tends towards 0. Unlike for the stress condition, no alternative energy criterion are proposed since Equation (9) results from an indisputable energy conservation principle, see Equation (8).

In some works, $\mathcal{G}_c$ is not considered as a constant but may vary. For instance, Catalanotti and Camanho (2013) proposed to consider increasing $\mathcal{G}_c$ as a function of the crack length during crack propagation in composite laminates that follows a R-curve obtained experimentally. Configurations where $\mathcal{G}_c$ depends on the crack velocity may be also encountered (Kanninen
The mode mixity $\psi$ can be defined by either shear to opening IERR or stress ratio. Carrère et al. (2015) proposed to define the mode mixity as the ratio between the opening and total IERR:

$$ \psi(S) = \frac{\mathcal{G}_{\text{inc}}^{I}(S)}{\mathcal{G}_{\text{inc}}(S)}. \qquad (15) $$

While basing the mode mixity definition on IERR contributions is congruent with the definition of the dependency of $\mathcal{G}_c$ to the mode mixity, it requires the separate calculation of opening and shear IERR contributions as functions of the crack surface.

A stress-based definition of the mode mixity was also adopted by several authors (Mantič 2009; García et al. 2012; García et al. 2014) so that:

$$ \psi(\mathbf{x}) = \arctan\left(\frac{\sigma_{nt}(\mathbf{x})}{\sigma_{nn}(\mathbf{x})}\right). \qquad (16) $$

The advantage of a stress-based mode mixity is the possibility to calculate $\psi$, and thus $\mathcal{G}_c$ and $\overline{\mathcal{G}}_c$ along the presupposed crack path directly based on the stress fields before initiation, which is consistent with the hypothesis of abrupt crack initiation (García et al. 2012). Nevertheless, Mantič (2009) also proposed to calculate a stress-based mode mixity $\psi(\ell)$ by considering the stress fields ahead of a crack tip of length $\ell$.

Once $\psi$ is determined, the critical ERR can be expressed as a function of $\psi$ following several models such as a power-law model (exponent $\varphi$) (Carrère et al. 2015):

$$ \mathcal{G}_c(S) = \left[ \left( \frac{\psi(S)}{\mathcal{G}_{\text{Ic}}} \right)^{\varphi} + \left( \frac{1 - \psi(S)}{\mathcal{G}_{\text{IIc}}} \right)^{\varphi} \right]^{\frac{-1}{\varphi}}, \qquad (17) $$

or Hutchinson and Suo (HS) model (Hutchinson et al. 1992):

$$ \mathcal{G}_c(\psi) = \mathcal{G}_{\text{Ic}}(1 + \tan^2((1 - \lambda)\psi)), \qquad (18) $$

where

$$ \lambda = 1 - \frac{2}{\pi} \tan^{-1}\left(\sqrt{\frac{\mathcal{G}_{\text{IIc}}}{\mathcal{G}_{\text{Ic}}}} - 1\right). \qquad (19) $$

When $\lambda = 1$, $\mathcal{G}_{\text{IIc}} = \mathcal{G}_{\text{Ic}}$ whereas $\mathcal{G}_{\text{IIc}}$ tends towards infinity when $\lambda \to 0$. So far, the choice of a quadratic or HS model seems to have a limited influence on the results as illustrated by Muñoz-Reja et al. (2016) in case of fiber-matrix debonding.

In configurations for which $\mathcal{G}_{\text{inc}}$ is monotonically increasing, the energy condition yields a lower bound of the admissible initiation crack surfaces (Figure 3(a)). Other configurations

<figure><img src="image_4.png" /><figcaption>Figure 3: Normalized IERR (thick line) and ERR (thin line) variations as a function of the crack surface for a) monotonic or b,c) non-monotonic variation of the IERR exhibiting b) a global or c) a local maximum. The filled zone in a) represents the crack surfaces for which the energy criterion is fulfilled. The circle, square and pentagon in b) respectively correspond to the crack surface at initiation, the lower and the upper bound of the crack arrest surface. The hatched zone in c) represents the prohibited crack initiation surfaces.</figcaption></figure>

where $\mathcal{G}_{\text{inc}}$ exhibits a global or a local maximum can also be encountered (Figures 3(b) and 3(c))
and Popelar 1987; Dally et al. 1985; Molnár, Gravouil, et al. 2020), which should be possibly implemented if dynamic aspects are considered (Doitrand et al. 2022). Another example consists in considering mode-dependent incremental and critical ERR in the energy criterion. Two approaches are classically used. The first one consists in separately considering opening and shear incremental to critical ERR contributions, which requires the calculation of both contributions and introduces the opening $G_{Ic}$ and shear $G_{IIc}$ critical ERR (Tran et al. 2012; Modniks et al. 2015):

$$ \frac{G_{\text{inc}}^{\text{I}}(S)}{G_{\text{Ic}}} + \frac{G_{\text{inc}}^{\text{II}}(S)}{G_{\text{IIc}}} = 1. \qquad (13) $$

The second approach consists in introducing the dependency of the critical ERR to the mode mixity $\psi$ and comparing the IERR to the average critical ERR over the crack path $\overline{G}_c(S)$ (Mantič 2009; García and Leguillon 2012; Carrère et al. 2015):

$$ G_{\text{inc}}(S) \geq \overline{G}_c(S), \quad \text{where} \quad \overline{G}_c = \overline{G}_c(S) = \frac{1}{S} \int_{\Gamma} G_c(\psi(\mathbf{x})) \, dS, \qquad (14) $$

and accounts for the local mode mixity $\psi$. The mode mixity can be defined by either shear to opening IERR or stress ratio. Carrère et al. (2015) proposed to define the mode mixity as the ratio between the opening and total IERR:

$$ \psi(S) = \frac{G_{\text{inc}}^{\text{I}}(S)}{G_{\text{inc}}(S)}. \qquad (15) $$

While basing the mode mixity definition on IERR contributions is congruent with the definition of the dependency of $G_c$ to the mode mixity, it requires the separate calculation of opening and shear IERR contributions as functions of the crack surface.

A stress-based definition of the mode mixity was also adopted by various authors (Mantič 2009; García and Leguillon 2012; García et al. 2014) so that

$$ \psi(\mathbf{x}) = \arctan\left(\frac{\sigma_{\text{nt}}(\mathbf{x})}{\sigma_{\text{nn}}(\mathbf{x})}\right). \qquad (16) $$

The advantage of a stress-based mode mixity is the possibility to calculate $\psi$, and thus $G_c$ and $\overline{G}_c$ along the presupposed crack path directly based on the stress fields before initiation, which is consistent with the hypothesis of abrupt crack initiation (García and Leguillon 2012). Nevertheless, Mantič (2009) also proposed to calculate a stress-based mode mixity $\psi(\ell)$ by considering the stress fields ahead of a crack tip of length $\ell$.

Once $\psi$ is determined, the critical ERR can be expressed as a function of $\psi$ following various models such as a power-law model (exponent $\varphi$) (Carrère et al. 2015)

$$ G_c(S) = \left[ \left( \frac{\psi(S)}{G_{Ic}} \right)^{\varphi} + \left( \frac{1 - \psi(S)}{G_{IIc}} \right)^{\varphi} \right]^{\frac{-1}{\varphi}} \qquad (17) $$

or the Hutchinson and Suo (HS) model (Hutchinson and Suo 1992)

$$ G_c(\psi) = G_{Ic}(1 + \tan^2((1 - \lambda)\psi)) \qquad (18) $$

where

$$ \lambda = 1 - \frac{2}{\pi} \tan^{-1}\left(\sqrt{\frac{G_{IIc}}{G_{Ic}} - 1}\right). \qquad (19) $$

The parameter $\lambda$ varies between 0 and 1. When $\lambda = 1$, $G_{IIc} = G_{Ic}$ whereas $G_{IIc}$ tends towards infinity when $\lambda \to 0$. So far, the choice of a quadratic or HS model seems to have a limited influence on the results as illustrated by Muñoz-Reja et al. (2016) in case of fiber-matrix debonding. Nevertheless, considering mode dependency of $G_c$ allows identifying both the opening and shear critical ERR, for instance of fiber-matrix interface (Girard et al. 2023; Girard, Doitrand, Koohbor, Rinaldi, Godin, and Bikard 2024).
(Leguillon et al. 2015a; Carrère et al. 2015; Leguillon et al. 2015c; Leguillon et al. 2016; Weißgräeber et al. 2016; Doitrand et al. 2017a; Wei et al. 2022). Basically, the initiation crack surface cannot correspond to a value for which $G_{inc}$ is decreasing, or more generally for which $G_{inc}$ is smaller than the maximum value attained for smaller crack surfaces. Indeed, such a situation would imply that the energy criterion is actually fulfilled for a smaller applied loading and a smaller crack surface. An example of prohibited initiation crack surfaces is shown in Figure 3(c). As a consequence, if $G_{inc}$ exhibits a global maximum (Figure 3(b)), the admissible initiation surfaces are crack surfaces smaller than the crack surface maximizing $G_{inc}$. If $G_{inc}$ exhibits a local maximum (Figure 3(c)), the admissible initiation surfaces are crack surfaces smaller than the crack surface maximizing $G_{inc}$ or corresponding to $G_{inc}$ values larger than the local maximum. The initiation crack surface is actually determined by coupling the stress and the energy criterion, which is detailed in next section. Once the initiation surface $S_c$ is determined, subsequent crack propagation can be assessed using Griffith's criterion (Equation (7)). Three situations may be encountered:

* In case of monotonically increasing $G_{inc}$ (Figure 3(a)), it can be deduced from Equation (11) that $dG/dS(S_c) > 0$ and $G(S_c) \ge G_{inc}(S_c)$. Since $G_{inc}(S_c) = G_c$, therefore $G(S_c) \ge G_c$ which means that unstable crack propagation occurs just after initiation and results in the final specimen failure.

* In case where $G_{inc}$ exhibits a maximum and that the initiation crack surface maximizes $G_{inc}$, then $G(S_c) = G_{inc}(S_c) = G_c$. Moreover, $G$ is decreasing for crack surface larger than the initiation crack surface. As a consequence stable crack propagation occurs after initiation when increasing the loading, which means that the initiation crack surface is also an arrest surface.

* In case where $G_{inc}$ exhibits a maximum and that the initiation crack surface is smaller than the one maximizing $G_{inc}$ (Figure 3(b)), $G_{inc}(S_c) = G_c$, $dG/dS(S_c) > 0$ and $G(S_c) \ge G_c$. Therefore, unstable crack propagation occurs at least until an arrest surface lower bound $S_{min}$ verifying $G(S_{min}) < G_c$. If all the available energy is used for crack propagation, the crack may possibly propagate up to an arrest surface upper bound $S_{max}$ veryfing $\int_{S_c}^{S_{min}} (G - G_c)dS = \int_{S_{min}}^{S_{max}} (G_c - G)dS$ (Leguillon et al. 2016; Weißgräeber et al. 2016; Doitrand et al. 2017b), see Figure 3(b).

A last category of configurations which can sometimes be encountered is the case of decreasing $G_{inc}$, for instance for small crack lengths under remote shear stress (Molnár et al. 2020a), for a crack approaching an interface in the case of a weak singularity (Leguillon et al. 2013a) or in the case of a strong singularity (Leguillon et al. 2000; Leguillon et al. 2012a; Aranda et al. 2023) such as a crack in a stiff material that impinges an interface with a more compliant material Aranda et al. 2023. The energy criterion thus provides an upper bound for admissible initiation crack surfaces, the admissible initiation surfaces are finally determined as the surfaces for which both criteria are fulfilled.

## 2.3 The coupled criterion

The CC consists in combining the above-mentioned stress (Equation (1)) and energy (Equation (9)) conditions in order to determine the initiation loading level $U_c$ and crack topology $\Gamma_c$ (corresponding crack surface $S_c$). In a general way, it can be formulated as an optimization problem to determine the minimum loading for which both conditions are fulfilled:

$$ \left\{ \begin{array}{l} U_c = \min\{U, \Gamma_{\sigma}(U) \cap \Gamma_G(U) \neq \emptyset\} \\ \Gamma_c = \Gamma_{\sigma}(U_c) \cap \Gamma_G(U_c) \end{array} \right. \qquad (20) $$

In the often-encountered particular case of linear elasticity and small deformation assumption, the stress and IERR are proportional respectively to the loading and the square loading. The stress and energy conditions thus write:

$$ \left\{ \begin{array}{l} G_{inc}(S, U) = A(S)U^2 \ge G_c, \\ \sigma(S, U) = k(S)U \ge \sigma_c, \end{array} \right. \qquad (21) $$
<figure><img src="image_4.png" /><figcaption>**Figure 3** Normalized IERR (thick line) and ERR (thin line) variations as a function of the crack surface for (a) monotonic or (b,c) non-monotonic variation of the IERR exhibiting (b) a global or (c) a local maximum. The filled zone in (a) represents the crack surfaces for which the energy criterion is fulfilled. The circle, square and pentagon in (b) respectively correspond to the crack surface at initiation, the lower and the upper bound of the crack arrest surface. The filled zone in (c) represents the prohibited crack initiation surfaces.</figcaption></figure>

In configurations for which $G_{inc}$ is monotonically increasing, the energy condition yields a lower bound of the admissible initiation crack surfaces, see Figure 3(a). Other configurations where $G_{inc}$ exhibits a global or a local maximum can also be encountered, see Figures 3(b) and 3(c) (Leguillon, Martin, and Lafarie-Frenot 2015; Carrère et al. 2015; Leguillon, Sevecek, et al. 2015; Leguillon et al. 2016; Weißgraeber, Hell, et al. 2016; Doitrand, Fagiano, Carrère, et al. 2017; Wei et al. 2022). Basically, the initiation crack surface cannot correspond to a value for which $G_{inc}$ is decreasing, or more generally for which $G_{inc}$ is smaller than the maximum value attained for smaller crack surfaces. Indeed, such a situation would imply that the energy criterion is actually fulfilled for a smaller applied loading and a smaller crack surface. An example of prohibited initiation crack surfaces is shown in Figure 3(c). As a consequence, if $G_{inc}$ exhibits a global maximum, in Figure 3(b), the admissible initiation surfaces are crack surfaces smaller than the crack surface maximizing $G_{inc}$. If $G_{inc}$ exhibits a local maximum, in Figure 3(c), the admissible initiation surfaces are crack surfaces smaller than the crack surface maximizing $G_{inc}$ or corresponding to $G_{inc}$ values larger than the local maximum. The initiation crack surface is actually determined by coupling the stress and the energy criterion, which is detailed in next section. Once the initiation surface $S_c$ is determined, subsequent crack propagation can be assessed using Griffith's criterion, see Equation (7). Three situations may be encountered:

* In case of monotonically increasing $G_{inc}$, Figure 3(a), it can be deduced from Equation (11) that $dG/dS(S_c) > 0$ and $G(S_c) \ge G_{inc}(S_c)$. Since $G_{inc}(S_c) = G_c$, then $G(S_c) \ge G_c$ which means that unstable crack propagation occurs just after initiation and results in the final specimen failure.

* In case where $G_{inc}$ exhibits a maximum and that the initiation crack surface maximizes $G_{inc}$, then $G(S_c) = G_{inc}(S_c) = G_c$. Moreover, $G$ is decreasing for crack surface larger than the initiation crack surface. As a consequence stable crack propagation occurs after initiation when increasing the loading, which means that the initiation crack surface is also an arrest surface.

* In case where $G_{inc}$ exhibits a maximum and that the initiation crack surface is smaller than the one maximizing $G_{inc}$, Figure 3(b), $G_{inc}(S_c) = G_c$, $dG/dS(S_c) > 0$ and $G(S_c) \ge G_c$. Therefore, unstable crack propagation occurs at least until an arrest surface lower bound $S_{min}$ verifying $G(S_{min}) < G_c$. As illustrated in Figure 3(b), if all the available energy is used for crack propagation, the crack may possibly propagate up to an arrest surface upper bound $S_{max}$ verifying (Leguillon et al. 2016; Weißgraeber, Hell, et al. 2016; Doitrand, Fagiano, Hild, et al. 2017)

$$ \int_{S_c}^{S_{min}} (G - G_c) \, dS = \int_{S_{min}}^{S_{max}} (G_c - G) \, dS. \qquad (20) $$

A last category of configurations which can sometimes be encountered is the case of decreasing $G_{inc}$, for instance for small crack lengths under remote shear stress (Molnár, Doitrand, et al. 2020), for a crack approaching an interface in the case of a weak singularity (Leguillon and Martin 2013a) or in the case of a strong singularity (Leguillon et al. 2000; Leguillon and Martin 2012; Aranda and Leguillon 2023) such as a crack in a stiff material that impinges an interface with a more compliant material (Aranda and Leguillon 2023). The energy criterion thus provides
where A and k are functions depending on the geometry and material properties. As a consequence, the loading that must be applied to fulfill the stress (U<sub>σ</sub>) or the energy (U<sub>G</sub>) criterion for a crack surface S can be derived:

$$
\left\{
\begin{array}{l}
U_G(S) = \sqrt{\frac{G_c}{A(S)}}, \\
U_\sigma(S) = \frac{\sigma_c}{k(S)}.
\end{array}
\right.
\qquad (22)
$$

The problem given in Equation (20) thus reverts to

$$
\left\{
\begin{array}{l}
U_c = \min_S \{ \max\{U_G(S), U_\sigma(S)\} \}, \\
S_c = \{S, \max\{U_G(S), U_\sigma(S)\} = U_c\}.
\end{array}
\right.
\qquad (23)
$$

For monotonic variations of the stress and IERR, it can be shown that the problem given in Equation (23) admits at most one solution S<sub>c</sub>, which verifies:

$$
\frac{A(S_c)}{k^2(S_c)} = \frac{G_c}{\sigma_c^2} \qquad (24)
$$

There are two equivalent ways to represent the CC solution, either through stress and energy
criterion variations as a function of the crack surface (Martin et al. 2010; Martin et al. 2012;
Leguillon et al. 2015b) (Figures 4(a) and 4(c)), or through variation of the required loadings

<figure><img src="image_5.png" /><figcaption>**Figure 4:** a,c) Normalized stress (dashed line) and IERR (solid lines) variations as a function of the crack surface obtained for a loading either smaller than (thin lines) or equal to (thick lines) the initiation loading for either a) monotonic or c) non-monotonic variations of the IERR. b,d) Required loading to fulfill the stress (dashed line) or the energy (solid line) criteria as a function of the crack surface evidencing the configurations for which the CC is fulfilled (hatched zone) for either b) monotonic or d) non-monotonic variations of the IERR.</figcaption></figure>
an upper bound for admissible initiation crack surfaces, the admissible initiation surfaces are finally determined as the surfaces for which both criteria are fulfilled.

## 2.3 Coupling stress and energy criteria

The CC consists in combining the above-mentioned stress (Equation (1)) and energy (Equation (9)) conditions in order to determine the initiation loading level $U_c$ and crack topology $\Gamma_c$ (corresponding crack surface $S_c$). In a general way, it can be formulated as an optimization problem to determine the minimum loading for which both conditions are fulfilled:

$$ \begin{cases} U_c = \min\{U, \Gamma_\sigma(U) \cap \Gamma_G(U) \neq \emptyset\}, \\ \Gamma_c = \Gamma_\sigma(U_c) \cap \Gamma_G(U_c). \end{cases} \qquad (21) $$

In the often-encountered particular case of linear elasticity and small deformation assumption, the stress and IERR are proportional respectively to the loading and the square loading. The stress and energy conditions thus read

$$ \begin{cases} \mathcal{G}_{\text{inc}}(S, U) = A(S)U^2 \geq \mathcal{G}_c, \\ \sigma(S, U) = k(S)U \geq \sigma_c \end{cases} \qquad (22) $$

where $A$ and $k$ are functions depending on the geometry and material properties. As a consequence, the loading that must be applied to fulfill the stress $U_\sigma$ or the energy $U_G$ criterion for a crack surface $S$ can be derived:

$$ U_G(S) = \sqrt{\frac{\mathcal{G}_c}{A(S)}}, \quad U_\sigma(S) = \frac{\sigma_c}{k(S)}. \qquad (23) $$

The problem given in Equation (21) thus reverts to

$$ \begin{cases} U_c = \min_S \{\max\{U_G(S), U_\sigma(S)\}\}, \\ S_c = \{S, \max\{U_G(S), U_\sigma(S)\} = U_c\}. \end{cases} \qquad (24) $$

For monotonic variations of the stress and IERR, it can be shown that the problem given in Equation (24) admits at most one solution $S_c$, which verifies

$$ \frac{A(S_c)}{k^2(S_c)} = \frac{\mathcal{G}_c}{\sigma_c^2}. \qquad (25) $$

There are two equivalent ways to represent the CC solution, either through stress and energy criterion variations as a function of the crack surface (Martin et al. 2010; Martin et al. 2012; Leguillon, Martin, Ševeček, et al. 2015), see Figures 4(a) and 4(c), or through variation of the required loadings to fulfill the stress and the energy criteria (Leguillon and Piat 2008; Leguillon 2011; Doitrand, Fagiano, Carrère, et al. 2017), see Figures 4(b) and 4(d). They are represented schematically in Figure 4 for either monotonic, Figures 4(a) and 4(b), or non-monotonic, Figures 4(c) and 4(d), variations of the IERR.

Figure 4(a) shows an example of monotonic normalized stress and IERR as a function of the crack surface obtained for a loading either smaller than or equal to the initiation loading. For a loading smaller than the initiation loading, there is no crack surfaces for which the stress and the energy criterion are fulfilled since the admissible initiation crack surface upper bound given by the stress criterion is smaller than the lower bound given by the energy criterion. Increasing the loading results in an increase of the upper bound given by the stress criterion and a decrease of the lower bound given by the energy criterion until both criteria are fulfilled for a crack surface $S_c$ that corresponds to the initiation crack surface. An equivalent way to determine the initiation loading and crack surface consists in studying the variation of the loading that must be applied to fulfill the stress or the energy criterion, Figure 4(b). For a given crack surface $S$, the loading required to fulfill the CC corresponds to the maximum loading between $U_G(S)$ and $U_\sigma(S)$, see Equation (24). As a consequence, all the configurations for which the CC is fulfilled are depicted
to fulfill the stress and the energy criteria (Leguillon et al. 2008a; Leguillon 2011; Doitrand et al. 2017a) (Figures 4(b) and 4(d)). They are represented schematically in Figure 4 for either monotonic (Figures 4(a) and 4(b)) or non-monotonic (Figures 4(c) and 4(d)) variations of the IERR.

Figure 4(a) shows an example of monotonic normalized stress and IERR variations as a function of the crack surface obtained for a loading either smaller than or equal to the initiation loading. For a loading smaller than the initiation loading, there is no crack surfaces for which the stress and the energy criterion are fulfilled since the admissible initiation crack surface upper bound given by the stress criterion is smaller than the lower bound given by the energy criterion. Increasing the loading results in an increase of the upper bound given by the stress criterion and a decrease of the lower bound given by the energy criterion until both criteria are fulfilled for a crack surface $S_c$ that corresponds to the initiation crack surface. An equivalent way to determine the initiation loading and crack surface consists in studying the variation of the loading that must be applied to fulfill the stress or the energy criterion (Figure 4(b)). For a given crack surface $S$, the loading required to fulfill the CC corresponds to the maximum loading between $U_G(S)$ and $U_\sigma(S)$ (Equation (23)). As a consequence, all the configurations for which the CC is fulfilled are depicted as the hatched zone. Among all these configurations, the initiation crack surface and loading are determined as the ones corresponding to the minimum applied loading for which the CC is fulfilled.

Both CC solution representations are shown in the case of a non-monotonic IERR (Figures 4(c) and 4(d)). If the initiation surface is smaller than the surface maximizing $G_{inc}$, the same representation as in Figures 4(a) and 4(b) can be used. Therefore, we focus on an initiation configuration for which the stress criterion is fulfilled on a crack surface larger than the one maximizing $G_{inc}$. For a loading smaller than the initiation loading, the stress criterion still provides an upper bound of the admissible crack surfaces, however there is no crack surface for which the energy criterion is fulfilled. Increasing the applied loading results in determining that the CC is fulfilled for the crack surface maximizing $G_{inc}$ (Figure 4(c)). This result is also obtained by studying the variation of the loading that must be applied to fulfill the stress or the energy criterion (Figure 4(d)). The hatched zone corresponds to all the configurations for which the CC is fulfilled, among which the crack surface maximizing $G_{inc}$ corresponds to the configuration for which the CC is fulfilled for the minimum applied loading.

## 2.4 Practical implementation of the CC

This section focuses on the different options to practically implement the CC.

### 2.4.1 Analytical application

An efficient way to implement the CC can be achieved provided the stress and IERR can be obtained from analytical formulae, which are generally available for simple geometries and loadings such as, e.g., circular holes (Sapora et al. 2018b; Sapora et al. 2018a; Cornetti et al. 2019b) or spherical cavities (Chao Correas et al. 2021; Ferrian et al. 2022). For instance, Stress Intensity Factor (SIF) variation as a function of the crack length can be obtained from handbook. Then, the ERR is obtained based on the relation between the SIF and the ERR $G = K_I^2/E'$ ($E' = E$ under plane stress assumption and $E' = E/(1 - v^2)$ under plane strain assumption or in 3D), which integration over the crack path gives the IERR (Equation (11)). Then the two criteria are combined following the approach presented in Section 2.3 to determine the initiation length and the corresponding remote loading. The analytical implementation of the CC thus only requires a minimization algorithm to solve Equation (23).

### 2.4.2 Matched Asymptotic (MA) approach

The MA approach was originally proposed for the CC implementation in the founding paper (Leguillon 2002). It consists in solving a two-scale problem. The first one is called the outer problem, in which the small crack that initiates is disregarded. This problem enables determining the relation between the displacement or loading applied to the structure under investigation, and the Generalized Stress Intensity Factors (GSIF) at the singular point or stress concentrator at which the crack initiates. The GSIF $K_j$ thus corresponds to the magnitude of the displacement
<figure><img src="image_5.png" /></figure>

**Figure 4** Representation of the coupled criterion solution: (a,c) Normalized stress (dashed line) and IERR (solid lines) as a function of the crack surface obtained for a loading either smaller than (thin lines) or equal to (thick lines) the initiation loading for either (a) monotonic or (c) non-monotonic variations of the IERR. (b,d) Required loading to fulfill the stress (dashed line) or the energy (solid line) criteria as a function of the crack surface evidencing the configurations for which the CC is fulfilled (hatched zone) for either (b) monotonic or (d) non-monotonic variations of the IERR.

as the hatched zone. Among these configurations, the initiation crack surface and loading are the ones corresponding to the minimum applied loading for which the CC is fulfilled.

Both CC solution representations are shown in the case of a non-monotonic IERR, Figures 4(c) and 4(d). If the initiation surface is smaller than the surface maximizing $G_{\text{inc}}$, the same representation as in Figures 4(a) and 4(b) can be used. Therefore, we focus on an initiation configuration for which the stress criterion is fulfilled on a crack surface larger than the one maximizing $G_{\text{inc}}$. For a loading smaller than the initiation loading, the stress criterion still provides an upper bound of the admissible crack surfaces, however there is no crack surface for which the energy criterion is fulfilled. Increasing the applied loading results in determining that the CC is fulfilled for the crack surface maximizing $G_{\text{inc}}$, Figure 4(c). This result is also obtained by studying the variation of the loading that must be applied to fulfill the stress or the energy criterion, Figure 4(d). The hatched zone corresponds to all the configurations for which the CC is fulfilled, among which the crack surface maximizing $G_{\text{inc}}$ corresponds to the configuration for which the CC is fulfilled for the minimum applied loading.

## 2.4 Practical implementation of the CC

This section focuses on the different options to practically implement the CC.

### 2.4.1 Analytical application

An efficient way to implement the CC can be achieved provided the stress and IERR can be obtained from analytical formulae, which are generally available for simple geometries and loadings such as, e.g., circular holes (Sapora et al. 2018; Sapora and Cornetti 2018; Cornetti and Sapora 2019) or spherical cavities (Chao Correas et al. 2021; Ferrian et al. 2022), or determined by
field around a singular point corresponding to a given mode, which writes:

$$u(r, \theta) = K_j r^{\lambda_j} u_j(\theta), \qquad (25)$$

where $\lambda_j$ is the characteristic exponent of the mode. It is obtained by solving an eigenvalue problem depending on the singular point geometry and material properties (Leguillon et al. 1987). The vector $u_j$ is the mode eigenvector corresponding to the eigenvalue $\lambda_j$. For instance, for a V-notch in 2D, the two first modes correspond to opening and shear modes, the characteristic exponents are given in Figure 5. The GSIF in the outer problem can be obtained based on a

<figure><img src="image_6.png" /><figcaption>Figure 5: (a) V-notch configuration and corresponding opening and shear mode (b) eigenvalues as a function of the V-notch angle and (c) eigenvectors as a function of the polar angle for $\beta = 90$ deg.</figcaption></figure>

contour integral calculation (Leguillon et al. 1987; Labossiere et al. 1999; Leguillon et al. 2007; Doitrand et al. 2020f) in which the dual mode intervenes. The exponent of the dual mode is $-\lambda_j$ in 2D and $-\lambda_j - 1$ in 3D. It can be determined based on the same eigenvalue problem as for the primal mode Leguillon et al. (1987). Even if the 3D exponents differ from the 2D exponents, the contour integral can be calculated in 3D (Doitrand et al. 2020e) and the same CC formulation can be adopted (Leguillon et al. 2014). It is also possible to evaluate the GSIF experimentally using the contour integral applied to displacement fields measured by digital image correlation (Doitrand et al. 2020d). Other approaches also exist to determine the GSIF, such as, e.g. the quasidual function method (Yosibash et al. 2005; Yosibash et al. 2021), boundary collocation method (Gross et al. 1972) or least square fitting method (Yu et al. 2012; Yi et al. 2017).

The second problem to solve in the MA approach is defined in the inner domain, a region close to the crack where the geometry of the structure far from the crack is disregarded. The boundary conditions applied to the inner problem correspond to the asymptotic displacement fields, which magnitude is determined by the GSIF. The stress variation and the potential energy difference are thus calculated in the inner domain including or not the crack. Both problems are finally matched by considering an area, close to the singular point in the outer problem, and far enough from it in the inner problem, where both solutions hold true. It yields the initiation crack length and GSIF from which the initiation applied loading or displacement is obtained. More details about the MA implementation of the CC for crack initiation at a V-notch are given in (Doitrand et al. 2020f) together with implementation examples.

It is worth noting that the MA approach is valid provided the initiation crack length is small with respect to the characteristic dimensions of the structure under consideration (Martin et al. 2016; Doitrand et al. 2020f). The MA approach of the CC was applied for crack initiation assessment at a hole (Felger et al. 2017b; Doitrand et al. 2021d), inclusion (Martin et al. 2020; Jiménez-Alfaro et al. 2022a), or bimaterial junction (Felger et al. 2019a), at a V-notch (Leguillon et al. 2017b; Zghal et al. 2018; Doitrand et al. 2020f), to handle a periodic array of cracks in thin layers (Leguillon et al. 2018a) or under mode I+III loading (Doitrand et al. 2018c). It also enabled rationalizing the meaning of the tensile stress in polycrystalline ceramics (Leguillon et al. 2018b), study their thermal shock induced cracking (Ricardo et al. 2020) and their small scale failure and size effect (Jiménez-Alfaro et al. 2021). The MA approach is also convenient to consider the influence of a specific mode such as, for instance, the T-stress (Leguillon et al. 2001; Leguillon et al. 2008a; Sapora et al. 2016a; Sapora et al. 2016b).
solving the governing equations (Methfessel et al. 2024). For instance, Stress Intensity Factor (SIF) variation as a function of the crack length can be obtained from handbook. Then, the ERR is obtained based on the relation between the SIF and the ERR $G = K_I^2/E'$ ($E' = E$ under plane stress assumption and $E' = E/(1 - v^2)$ under plane strain assumption or in 3D), which integration over the crack path gives the IERR, see Equation (11). Then the two criteria are combined following the approach presented in Section 2.3 to determine the initiation length and the corresponding remote loading. The analytical implementation of the CC thus only requires a minimization algorithm to solve Equation (24).

### 2.4.2 Matched Asymptotic approach

The Matched Asymptotic (MA) approach was originally proposed for the CC implementation in the founding paper (Leguillon 2002). It consists in solving a two-scale problem. The first one is called the outer problem, in which the small crack that initiates is disregarded. This problem enables determining the relation between the displacement or loading applied to the structure under investigation, and the Generalized Stress Intensity Factor (GSIF) at the singular point or stress concentrator at which the crack initiates. The GSIF $K_j$ thus corresponds to the magnitude of the displacement field around a singular point corresponding to a given mode, which reads

$$ \mathbf{u}(r, \theta) = K_j r^{\lambda_j} \mathbf{u}_j(\theta) \qquad (26) $$

where $\lambda_j$ is the characteristic exponent of the mode. It is obtained by solving an eigenvalue problem depending on the singular point geometry and material properties (Leguillon and Sanchez-Palencia 1987). The vector $\mathbf{u}_j$ is the mode eigenvector corresponding to the eigenvalue $\lambda_j$. For instance, for a V-notch in 2D, the two first modes correspond to opening and shear modes, the characteristic exponents are given in Figure 5. The GSIF in the outer problem can be

<figure><img src="image_6.png" /><figcaption>**Figure 5** Matched asymptotic approach of the coupled criterion: (a) V-notch configuration and corresponding opening and shear mode (b) eigenvalues as a function of the V-notch angle and (c) eigenvectors as a function of the polar angle for $\beta = 90^\circ$.</figcaption></figure>

obtained based on a contour integral calculation (Leguillon and Sanchez-Palencia 1987; Labossiere and Dunn 1999; Leguillon et al. 2007; Doitrand, Martin, et al. 2020) in which the dual mode intervenes. The exponent of the dual mode is $-\lambda_j$ in 2D and $-\lambda_j - 1$ in 3D. It can be determined based on the same eigenvalue problem as for the primal mode Leguillon and Sanchez-Palencia (1987). Even if the 3D exponents differ from the 2D exponents, the contour integral can be calculated in 3D (Doitrand, Leguillon, and Martin 2020) and the same CC formulation can be adopted (Leguillon et al. 2014). It is also possible to evaluate the GSIF experimentally using the contour integral applied to displacement fields measured by digital image correlation (Doitrand, Leguillon, and Estevez 2020). Other approaches also exist to determine the GSIF, such as, e.g. the quasidual function method (Yosibash et al. 2005; Yosibash and Schapira 2021), boundary collocation method (Gross and Mendelson 1972) or least square fitting method (Yu and Shi 2012; Yi et al. 2017).

The second problem to solve in the MA approach is defined in the inner domain, a region close to the crack tip where the geometry of the structure far from the crack is disregarded. The
### 2.4.3 Full finite element (FFE) approach

In configurations for which analytical solutions are not available and the validity conditions of the CC MA approach are not met, an alternative consists in performing a full finite element (FFE) calculation of the structure for several crack lengths. It consists in considering the full structure in which the crack is described instead of focusing only on the crack initiation location, as in the MA approach inner domain. Details and implementation example of the FFE approach are given in Doitrand et al. (2020). Some authors compared the MA and FFE implementation of the CC. For instance, Yosibash et al. (2006) and Priel et al. (2008) obtained differences smaller than 1% for the potential energy variations computed by both approaches for a sharp or a blunted V-notch for small enough cracks, these differences were larger for larger cracks. Martin et al. (2016) showed that the range of validity of the asymptotic approach can be determined provided that the material characteristic length $\ell_{mat} = E'G_c/\sigma_c^2$ (Irwin's length) is small enough compared to the characteristic dimension of the structure ($\ell_{mat} < 0.05h_s$ for the studied case of a bimaterial specimen submitted to a four point bending test, where $h_s$ is the substrate thickness). Doitrand et al. (2020) also compared FFE and MA approaches for crack initiation assessment at a V-notch. While similar results were obtained for small enough material characteristic length, significant differences were obtained with increasing $\ell_{mat}$, which was related to the increase in crack initiation length with increasing $\ell_{mat}$ and V-notch angle. They showed that the initiation length must be at least one order of magnitude smaller than the notch length to remain within the framework of the asymptotic analysis. A too large relative $\ell_{mat}$ may even lead to non-physical variations of the initiation imposed loading as a function of the V-notch angle, as also observed by Carpinteri et al. (2011).

### 2.4.4 Other implementations

Some alternative implementations to the above-mentioned approaches of the CC were established. Mantič (2014) proposed the so-called Principle of Minimum Total Energy subjected to a Stress Constraint (PMTE-SC). This approach relies on a reformulation of the CC as a constrained optimization problem that can be seen as a generalization of LEFM. It consists in minimizing the total energy (defined as the sum of the potential energy for a given crack surface and the energy dissipated up to the formation of the new crack) while constraining the space of admissible crack surfaces based on the stress criterion. Muñoz-Reja et al. (2022) illustrated, based on double cantilever beam test simulations, that the PMTE-SC implementation is equivalent to the classical CC implementation. They highlighted that the initiation finite crack jump can be associated to a tunneling effect across the total energy barrier.

An implementation of the CC dedicated to study interface crack initiation and propagation was proposed by combining Linear Elastic Brittle Interface Model (LEBIM) together with FFM, which consisted in describing the interface as a bed of linear elastic springs. It can be applied to cracked or uncracked geometries, for instance to the case of pull-push shear (Cornetti et al. 2012) or double cantilever beam (Muñoz-Reja et al. 2022) tests. It was extended to consider mixed mode loading using different criteria (Muñoz-Reja et al. 2016; Muñoz-Reja et al. 2020a) and implemented in a Boundary Element Method (BEM) code (Muñoz-Reja et al. 2020b), focusing on the required mesh refinement in the zone close to the crack tip in case of stiff interfaces (Muñoz-Reja et al. 2018).

Li et al. (2018) set up a numerical implementation of the CC using an iterative algorithm that first determines the potential crack locations based on a watershed flooding technique as the areas where the stress criterion is fulfilled. Then, the potential crack paths are obtained within these areas following a principal stress trajectory. If $N$ potential cracks are identified, then $2^N$ calculations are needed to find the cracked configuration maximizing the IERR. Then, the elements in interception with the predicted crack path are removed from the mesh. It thus requires a sufficiently fine mesh to well describe the crack path through element removal. The procedure is then repeated until either the energy or the stress criterion is not fulfilled, in which case the loading is increased. The originality of this implementation lies in the ability to describe both the initiation and the growth of one or more cracks simultaneously. It was successfully applied to crack initiation assessment at a V-notch, interaction between two cracks and laminate or woven composite multicracking (Li et al. 2018a). It also enabled considering nonlinear material
boundary conditions applied to the inner problem correspond to the asymptotic displacement fields, which magnitude is determined by the GSIF. The stress variation and the potential energy difference are thus calculated in the inner domain including or not the crack. Both problems are finally matched by considering an area, close to the singular point in the outer problem, and far enough from it in the inner problem, where both solutions hold true. It yields the initiation crack length and GSIF from which the initiation applied loading or displacement is obtained. More details about the MA implementation of the CC for crack initiation at a V-notch are given in (Doitrand, Martin, et al. 2020) together with implementation examples.

It is worth noting that the MA approach is valid provided the initiation crack length is small with respect to the characteristic dimensions of the structure under consideration (Martin et al. 2016; Doitrand, Martin, et al. 2020). The MA approach of the CC was applied for crack initiation assessment at a hole (Felger et al. 2017b; Doitrand and Leguillon 2021), inclusion (Martin et al. 2020; Jiménez-Alfaro and Leguillon 2022; Jiménez-Alfaro and Leguillon 2023), or bimaterial junction (Felger, Rosendahl, et al. 2019), at a V-notch (Leguillon and Yosibash 2017; Zghal et al. 2018; Doitrand, Martin, et al. 2020), to handle a periodic array of cracks in thin layers (Leguillon and Martin 2018) or under mode I+III loading (Doitrand and Leguillon 2018c). It also enabled rationalizing the meaning of the tensile stress in polycrystalline ceramics (Leguillon et al. 2018), study their thermal shock induced cracking (Faria Ricardo et al. 2020) and their small scale failure and size effect (Jiménez-Alfaro and Leguillon 2021). The MA approach is also convenient to consider the influence of a specific mode such as, for instance, the T-stress (Leguillon et al. 2001; Leguillon and Piat 2008; Sapora et al. 2016; Sapora and Mantič 2016).

### 2.4.3 Full finite element (FFE) approach

In configurations for which analytical solutions are not available and the validity conditions of the MA approach of the CC are not met, an alternative consists in performing a full finite element (FFE) calculation of the structure for several crack lengths. It consists in considering the full structure in which the crack is described instead of focusing only on the crack initiation location, as in the MA approach inner domain. Details and implementation example of the FFE approach are given in Doitrand, Martin, et al. (2020). Some authors compared the MA and FFE implementation of the CC. For instance, Yosibash et al. (2006) and Priel et al. (2008) obtained differences smaller than 1% for the potential energy variations computed by both approaches for a sharp or a blunted V-notch for small enough cracks, these differences were larger for larger cracks. Martin et al. (2016) showed that the range of validity of the asymptotic approach can be determined provided that the material characteristic length $\ell_{\text{mat}} = E'G_c/\sigma_c^2$ (Irwin's length) is small enough compared to the characteristic dimension of the structure ($\ell_{\text{mat}} < 0.05h_s$ for the studied case of a bimaterial specimen submitted to a four point bending test, where $h_s$ is the substrate thickness). Doitrand, Martin, et al. (2020) also compared FFE and MA approaches for crack initiation assessment at a V-notch. While similar results were obtained for small enough material characteristic length, significant differences were obtained with increasing $\ell_{\text{mat}}$, which was related to the increase in crack initiation length with increasing $\ell_{\text{mat}}$ and V-notch angle. They showed that the initiation length must be at least one order of magnitude smaller than the notch length to remain within the framework of the asymptotic analysis. A too large relative $\ell_{\text{mat}}$ may even lead to non-physical variations of the initiation imposed loading as a function of the V-notch angle, as also observed by Carpinteri et al. (2011).

### 2.4.4 Other implementations

Some alternative implementations to the above-mentioned approaches of the CC were established. Mantič (2014) proposed the so-called Principle of Minimum Total Energy subjected to a Stress Constraint (PMTE-SC). This approach relies on a reformulation of the CC as a constrained optimization problem that can be seen as a generalization of LEFM. It consists in minimizing the total energy (defined as the sum of the potential energy for a given crack surface and the energy dissipated up to the formation of the new crack) while constraining the space of admissible crack surfaces based on the stress criterion. Muñoz-Reja et al. (2022) illustrated, based on double cantilever beam test simulations, that the PTME-SC implementation is equivalent to the classical
behavior (Li et al. 2019).

Zhang et al. (2018) implemented the CC under a peridynamic framework to study crack initiation at circular holes. They showed that except near the hole edge, similar strain and stress values were obtained using either peridynamics or FE modeling. As a consequence, they showed that peridynamics yielded similar failure stresses as the one obtained using FFM (Zhang et al. 2018; Zhang et al. 2019).

Le Pavic et al. (2018) combined the CC implementation with a surrogate model to predict bonded joint failure. While the CC was actually calculated on a limited number of geometry configurations, the model response was obtained through spatial interpolation using a Kriging model, which allowed interpolating failure loading for configurations for which the CC was not directly calculated, and thus reduced the computational costs when several parameters varied in a given range.

Rettl et al. (2022) set-up a two scale methodology to predict crack initiation under opening mode from 2D notches that contain a sharp edge or crack. Crack models were first pre-computed for five deformation modes and then scaled to an actual notch using dimensional analysis, linear superposition and meta-modeling. The IERR was calculated using a relaxation method for the boundary conditions. On one hand, this method is computationally efficient and yields failure load differences smaller than 15% compared to FFE approach. On the other hand, it is limited to study small initiation crack corresponding to Irwin lengths smaller than a certain threshold.

## 2.5 Prescribed displacement or force

The CC FFE implementation requires boundary conditions that are usually implemented as imposed displacement or force, see Figure 6(a). There is not much difference in the CC

<figure><img src="image_7.png" /><figcaption>**Figure 6:** a) Three-point bending configuration under imposed displacement or force. b) IERR to critical ERR ratio as a function of dimensionless crack length and c) failure stress as a function of dimensionless specimen size obtained using imposed displacement or force.</figcaption></figure>

implementation except that, for imposed displacement, $\Delta W_{ext}$ vanishes in Equation (9). Before crack initiation, both approaches yield the same stress fields and thus the same stress criterion. Nevertheless, using one or the other approach results in different IERR (Figure 6(b)). Indeed, if the imposed displacement is kept constant while the crack length increases, it may lead to a decrease in the corresponding force. On the contrary, keeping the same force level while increasing the crack length results in an increase of the application point displacement. Figure 6 shows an example of IERR obtained under bending using either imposed displacement or force as a function of the crack length to specimen width ratio. Differences smaller than 1% are obtained provided the crack length remains sufficiently small ($\ell < 0.1h$). For larger crack lengths, major differences are observed since the IERR is monotonically increasing for imposed force whereas it exhibits a maximum for imposed displacement. It means that whatever the material properties, unstable crack propagation necessarily occurs under imposed force, whereas either unstable (for $\ell_c$ smaller than the length which maximizes $G_{inc}$) or stable (if $\ell_c$ maximizes $G_{inc}$) crack propagation occurs under imposed displacement. The initiation crack length depends on the material properties through $\ell_{mat}$. Figure 6(c) shows the normalized failure stress variation as a function of the normalized specimen width obtained using either imposed displacement or force (Cornetti et al. 2006; Doitrand et al. 2021c). No major differences are observed provided the
CC implementation. They highlighted that the initiation finite crack jump can be associated to a tunneling effect across the total energy barrier.

An implementation of the CC dedicated to study interface crack initiation and propagation was proposed by combining Linear Elastic Brittle Interface Model together with FFM, which consisted in describing the interface as a bed of linear elastic springs. It can be applied to cracked or uncracked geometries, for instance to the case of pull-push shear (Cornetti et al. 2012) or double cantilever beam (Muñoz-Reja et al. 2022) tests. It was extended to consider mixed mode loading using different criteria (Muñoz-Reja et al. 2016; Muñoz-Reja, Cornetti, et al. 2020) and implemented in a Boundary Element Method code (Muñoz-Reja, Távara, Mantič, and Cornetti 2020), focusing on the required mesh refinement in the zone close to the crack tip in case of stiff interfaces (Muñoz-Reja et al. 2018).

Li et al. (2018) set up a numerical implementation of the CC using an iterative algorithm that first determines the potential crack locations based on a watershed flooding technique as the areas where the stress criterion is fulfilled. Then, the potential crack paths are obtained within these areas following a principal stress trajectory. If $N$ potential cracks are identified, then $2^N$ calculations are needed to find the cracked configuration maximizing the IERR. Then, the elements in interception with the predicted crack path are removed from the mesh. It thus requires a sufficiently fine mesh to well describe the crack path through element removal. The procedure is then repeated until either the energy or the stress criterion is not fulfilled, in which case the loading is increased. The originality of this implementation lies in the ability to describe both the initiation and the growth of one or more cracks simultaneously. It was successfully applied to crack initiation assessment at a V-notch, interaction between two cracks and laminate or woven composite multicracking (Li and Leguillon 2018). It also enabled considering nonlinear material behavior (Li et al. 2019).

H. Zhang and Qiao (2018) implemented the CC under a peridynamic framework to study crack initiation at circular holes. They showed that except near the hole edge, similar strain and stress values were obtained using either peridynamics or FE modeling. As a consequence, they showed that peridynamics yielded similar failure stresses as the one obtained using FFM (H. Zhang and Qiao 2018; H. Zhang et al. 2019).

Le Pavic et al. (2018) combined the CC implementation with a surrogate model to predict bonded joint failure. While the CC was actually calculated on a limited number of geometry configurations, the model response was obtained through spatial interpolation using a Kriging model, which allowed interpolating failure loading for configurations for which the CC was not directly calculated, and thus reduced the computational costs when several parameters varied in a given range.

Rettl et al. (2022) set-up a two scale methodology to predict crack initiation under opening mode from 2D notches that contain a sharp edge or crack. Crack models were first pre-computed for five deformation modes and then scaled to an actual notch using dimensional analysis, linear superposition and meta-modeling. The IERR was calculated using a relaxation method for the boundary conditions. On one hand, this method is computationally efficient and yields failure load differences smaller than 15% compared to FFE approach. On the other hand, it is limited to study small initiation crack corresponding to Irwin lengths smaller than a certain threshold.

Martin and Carrère (2023) proposed an analytical description of the size effect law for notched structures, which gives the notch strength as a function of a characteristic structure length. This formulation only depends on three parameters that must be identified, requiring only three finite element calculations. The proposed approach is thus an efficient procedure for applying the CC.

## 2.5 Prescribed displacement or force

The FFE implementation of the CC requires boundary conditions that are usually implemented as imposed displacement or force, see Figure 6(a). There is not much difference in the CC implementation except that, for imposed displacement, $\Delta W_{ext}$ vanishes in Figure 6(a). Before crack initiation, both approaches yield the same stress fields and thus the same stress criterion. Nevertheless, using one or the other approach results in different IERR (Figure 6(b)). Indeed, if the imposed displacement is kept constant while the crack length increases, it may lead to a decrease
material characteristic length remains small enough compared to the specimen dimensions. These
differences are highlighted either for materials exhibiting a large characteristic length, either
when decreasing the specimen size, as discussed by Jiménez-Alfaro et al. (2021). The boundary
conditions may thus carefully be chosen for these configurations. An experimental validation to
evaluate the failure force differences between imposed force or displacement is still missing. Such
a validation is not straightforward since usually, force-controlled experiments actually rely on
using a retroaction loop to adapt the imposed displacement so as to follow a prescribed force rate.
As a consequence, if the characteristic crack initiation time is smaller than the retroaction loop
duration, it is not trivial to really ensure a force control during crack initiation.

## 2.6 Unknown crack path

The CC implementation becomes challenging when the crack path is a priori unknown, which
means that both the initiation location and crack geometry must be determined. In 2D, it consists
in determining the crack length and orientation, which can be achieved by applying the CC
for different crack angles. The initiation angle will be determined as the one minimizing the
imposed loading obtained with the CC, which also provides the initiation crack length. This
situation is for instance encountered when studying mixed mode crack initiation in homogeneous
material (Leguillon et al. 2007; Priel et al. 2008; Doitrand et al. 2020f; Molnár et al. 2020a). If
there are several potential crack locations, the above-mentioned procedure must be repeated for
each location, such as, e.g. to study the competition between interface or bulk cracking in a
bimaterial configuration (Felger et al. 2019a) or to assess symmetric or asymmetric crack initiation
configurations (García et al. 2015; Doitrand et al. 2019b). The potential crack paths may be
determined based on the locations where the stress is maximum (Li et al. 2018b), since it is likely
that crack initiation may be more favorable at locations where the stress is larger. Nevertheless,
configurations in which crack initiation is driven by energy only can be encountered (Leguillon
2002; Martin et al. 2008; Leguillon et al. 2013a; Weißgraeber et al. 2016; Doitrand et al. 2017a), i.e.
in which the stress criterion is actually fulfilled on a larger zone than the one corresponding to
the initiation crack. In such configurations, the crack path may be determined based on the
energy criterion (Doitrand et al. 2017a), which is however computationally less efficient. Indeed,
it consists in determining, among all possible crack shapes, the one that mimizes the imposed
loading required to fulfill both criteria.

# 3 Fundamental developments of the CC

In the previous review paper by Weißgraeber et al. (2016), several challenges and possible
extensions of the CC were mentioned, including fatigue loadings, 3D configurations or dynamic
loadings, to which nonlinearities can also be added. Some works tried to address these challenges,
they are summarized in this section.

## 3.1 3D

The 3D CC formulation mainly differs from the 2D CC formulation as crack surfaces are calculated
instead of crack lengths, otherwise it remains similar. The major challenge to face in 3D is
that the crack path may be *a priori* unknown. In 2D, the configurations maximizing the energy
dissipation generally correspond to the ones maximizing the stress prior to crack initiation
(Erdogan et al. 1963; Goldstein et al. 1974; Pham et al. 2017). This is not always the case in 3D: the
initiation crack location minimizing the imposed loading based on the CC was found to lie near,
but not at the exact same location than the one corresponding to the stress criterion maximum in
woven composites (Doitrand et al. 2017a). Besides, it was shown that the ERR alone was not
able to provide the crack direction observed experimentally under mode III loading (Mittelman
et al. 2015). Moreover, a 2D crack is generally described by two parameters: the crack length
and the crack angle. The 3D crack surface definition is way more complicated since the crack
could theoretically be described by an infinite number of 3D surfaces. Even assuming a planar
crack path, it is necessary to first determine the crack plane, then the crack front shape in this
plane, which can also be described by an infinite number of 2D curves. Mittelman et al. (2014)
showed that the IERR corresponding to a V-notch elliptical shaped crack may vary significantly
in the corresponding force. On the contrary, keeping the same force level while increasing the
crack length results in an increase of the application point displacement. Figure 6 shows an
example of IERR obtained under bending using either imposed displacement or force as a function
of the crack length to specimen width ratio. Differences smaller than 1% are obtained provided
the crack length remains sufficiently small ($\ell < 0.1h$). For larger crack lengths, major differences
are observed since the IERR is monotonically increasing for imposed force whereas it exhibits a
maximum for imposed displacement. It means that whatever the material properties, unstable
crack propagation necessarily occurs under imposed force, whereas either unstable (for $\ell_c$ smaller
than the length which maximizes $G_{inc}$) or stable (if $\ell_c$ maximizes $G_{inc}$) crack propagation occurs
under imposed displacement. The initiation crack length depends on the material properties
through $\ell_{mat}$. Figure 6(c) shows the normalized failure stress variation as a function of the
normalized specimen width obtained using either imposed displacement or force (Cornetti
et al. 2006; Doitrand, Henry, and Meille 2021). No major differences are observed provided the
material characteristic length remains small enough compared to the specimen dimensions.
These differences are highlighted either for materials exhibiting a large characteristic length,
either when decreasing the specimen size, as discussed by Jiménez-Alfaro and Leguillon (2021).
The boundary conditions may thus carefully be chosen for these configurations. An experimental
validation to evaluate the failure force differences between imposed force or displacement is still
missing. Such a validation is not straightforward since usually, force-controlled experiments
actually rely on using a retroaction loop to adapt the imposed displacement so as to follow a
prescribed force rate. As a consequence, if the characteristic crack initiation time is smaller
than the retroaction loop duration, it is not trivial to really ensure a force control during crack
initiation.

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

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

<figure><img src="image_9.png" /><figcaption>**Figure 6** Fracture size effect explained by the coupled criterion: (a) Three-point bending configuration under imposed displacement or force, (b) IERR to critical ERR ratio as a function of dimensionless crack length and (c) failure stress as a function of dimensionless specimen size obtained using imposed displacement or force.</figcaption></figure>

## 2.6 Unknown crack path

The CC implementation becomes challenging when the crack path is a priori unknown, which means that both the initiation location and crack geometry must be determined. In 2D, it consists in determining the crack length and orientation, which can be achieved by applying the CC for different crack angles. The initiation angle will be determined as the one minimizing the imposed loading obtained with the CC, which also provides the initiation crack length. This situation is for instance encountered when studying mixed mode crack initiation in homogeneous material (Leguillon et al. 2007; Priel et al. 2008; Doitrand, Martin, et al. 2020; Molnár, Doitrand, et al. 2020). If there are several potential crack locations, the above-mentioned procedure must be repeated for each location, such as, e.g. to study the competition between interface or bulk cracking in a bimaterial configuration (Felger, Rosendahl, et al. 2019) or to assess symmetric or asymmetric crack initiation configurations (García et al. 2015; Doitrand et al. 2019b). The potential crack paths may be determined based on the locations where the stress is maximum (Li et al. 2018), since it is likely that crack initiation may be more favorable at locations where the
depending on the ellipse aspect ratio and that the choice of the crack shape is thus crucial in a 3D
analysis. Thus, how to proceed in 3D to define the crack shape and what are the differences
brought by the 3D approach compared to the 2D case?

**3.1.1 Assuming the 3D crack path**

A first way to implement the CC in 3D consists in assuming a crack path and applying stress and energy criteria corresponding to increasingly large cracks following this path. This strategy was for instance adopted to simulate 3D transverse crack initiation in 0/90/0 composite laminates (García et al. 2016; García et al. 2018; Hamam et al. 2022). The crack geometry was set to be in the form of rectangles with slightly rounded corners, thus being described by two scalar variables corresponding to crack extensions along the ply length and thickness. García et al. (2016) showed, based on the 3D application of the CC, that for thick laminates, crack initiation at the free edge is preferred whereas for thin laminates, no significant preference is found. This result may be explained by the stress concentration in the ply at the specimen free faces since crack initiation is rather driven by the stress criterion for thick laminates and by the energy criterion for thin laminates.

Transverse crack initiation was also studied in woven composites by Doitrand et al. (2017); Doitrand et al. (2017). A tentative to base the possible crack surfaces on stress isocontours was not successful since crack initiation was found to be controlled only by the energy criterion. As a consequence, idealized crack surfaces spanning the whole ply thickness and having straight crack fronts were assumed. Debonding ahead of the crack tip was also modeled assuming straight debonding fronts. The stress concentration at specimen free faces was disregarded due to the use of a periodic representative unit cell. The proposed approach enabled modeling the initiation and propagation of multiple transverse cracks and debondings considering interactions between cracks and determining crack and debonding density variations as a function of the applied loading.

Yosibash et al. (2016) applied the CC in 3D to assess crack initiation at a V-notch under mode I+II+III loading. They first determined the crack plane based on the stress criterion only and noticed that crack shape determination in such a plane based on the stress criterion would require different crack shapes depending on the crack plane orientation. They thus assumed circular sector shaped crack whatever the crack orientation, which means that the stress criterion was not homogeneously met on the whole crack surface prior to initiation. As a consequence, an average stress criterion over the circular sector shaped crack was used to implement the CC. Note that the classical stress criterion could also be applied in such a situation: it means that the minimum stress over the considered crack shape before initiation must be larger or equal than the material tensile strength.

Papšík et al. (2023) assumed elliptical crack shape to simulate tunneling crack in layered ceramics containing thermal residual stresses. Exploiting the configuration symmetry, the crack was modeled by releasing Dirichlet boundary conditions on the symmetry plane. They showed that the thickness size effect is governed by both stress and energy criteria. Moreover, they highlighted that crack initiation may be prevented if the energy criterion is not fulfilled, even if the stress reaches a 99% probability of failure according to Weibull’s analysis. They thus determined the conditions in terms of layer thickness and residual stresses for which edge crack or tunneling cracks are likely to initiate.

**3.1.2 Stress isocontour-based 3D crack path**

A natural way to determine the possible crack shapes consists in ensuring that they strictly fulfill the stress criterion, i.e. using shapes based on the stress criterion isocontours. It was first proposed by Leguillon (2014) who extended the CC to 3D interface corner crack prediction at a bimaterial junction. The 3D extension was established based the 3D singularity theory (Leguillon et al. 1987; Leguillon 1995; Yosibash 1997b; Yosibash 1997a) and matched asymptotic expansions. Despite some differences in its establishment compared to the 2D formulation, the 3D CC formulation was found to finally be similar to the 2D case. The crack plane was a priori known since it corresponded to the bimaterial interface. An additional hypothesis on the crack geometry was needed compared to the 2D case, consisting in defining the possible crack surfaces based on the tensile stress isocontours which were slightly modified to consider bounded crack
stress is larger. Nevertheless, configurations in which crack initiation is driven by energy only
can be encountered (Leguillon 2002; Martin et al. 2008; Leguillon and Martin 2013a; Weißgraeber,
Hell, et al. 2016; Doitrand, Fagiano, Carrère, et al. 2017), i.e. in which the stress criterion is
actually fulfilled on a larger zone than the one corresponding to the initiation crack. In such
configurations, the crack path may be determined based on the energy criterion (Doitrand,
Fagiano, Carrère, et al. 2017), which is however computationally less efficient. Indeed, it consists
in determining, among all possible crack shapes, the one that mimizes the imposed loading
required to fulfill both criteria. Girard, Doitrand, Koohbor, Rinaldi, Godin, Long, and Bikard (2024)
tested all possible fiber-matrix interface single debonding configurations so as to evaluate the
optimal crack shapes. It was found that for small interface brittleness number, the stress-based
debonding shapes are the optimal one, whereas for large interface brittleness number, the optimal
crack shapes are those that maximize the energy criterion. Finally, for intermediate interface
brittleness number, the optimal initiation debonding shapes are neither those predicted by
the stress criterion nor those predicted by the energy criterion, but lie in between these two
configurations.

3 Fundamental developments of the CC

In the previous review paper by Weißgraeber, Leguillon, et al. (2016), several challenges and
possible extensions of the CC were mentioned, including fatigue loadings, 3D configurations or
dynamic loadings, to which nonlinearities can also be added. Some works tried to address these
challenges, they are summarized in this section.

3.1 3D formulation

The 3D CC formulation mainly differs from the 2D CC formulation as crack surfaces are calculated
instead of crack lengths, otherwise it remains similar. The major challenge to face in 3D is that the
crack path may be a *priori* unknown. In 2D, the configurations maximizing the energy dissipation
generally correspond to the ones maximizing the stress prior to crack initiation (Erdogan and Sih
1963; Gol'dstein and Salganik 1974; K. H. Pham et al. 2017). This is not always the case in 3D: the
initiation crack location minimizing the imposed loading based on the CC was found to lie near,
but not at the exact same location than the one corresponding to the stress criterion maximum
in woven composites (Doitrand, Fagiano, Carrère, et al. 2017). Besides, it was shown that the
ERR alone was not able to provide the crack direction observed experimentally under mode III
loading (Mittelman and Yosibash 2015). Moreover, a 2D crack is generally described by two
parameters: the crack length and the crack angle. The 3D crack surface definition is way more
complicated since the crack could theoretically be described by an infinite number of 3D surfaces.
Even assuming a planar crack path, it is necessary to first determine the crack plane, then the
crack front shape in this plane, which can also be described by an infinite number of 2D curves.
Mittelman and Yosibash (2014) showed that the IERR corresponding to a V-notch elliptical shaped
crack may vary significantly depending on the ellipse aspect ratio and that the choice of the crack
shape is thus crucial in a 3D analysis. Thus, how to proceed in 3D to define the crack shape and
what are the differences brought by the 3D approach compared to the 2D case?

3.1.1 Assuming the 3D crack path

A first way to implement the CC in 3D consists in assuming a crack path and applying stress and
energy criteria corresponding to increasingly large cracks following this path. This strategy was
for instance adopted to simulate 3D transverse crack initiation in 0/90/0 composite laminates
(García et al. 2016; García et al. 2018; Hamam et al. 2022). The crack geometry was set to be in the
form of rectangles with slightly rounded corners, thus being described by two scalar variables
corresponding to crack extensions along the ply length and thickness. García et al. (2016) showed,
based on the 3D application of the CC, that for thick laminates, crack initiation at the free edge is
preferred whereas for thin laminates, no significant preference is found. This result may be
explained by the stress concentration in the ply at the specimen free faces since crack initiation is
surfaces. Note that the MA approach of the CC yields the critical GSIF at initiation. The relation between the GSIF at a singular point in a real specimen and the applied loading can be obtained by implementing a contour integral around the singular point. This approach was implemented by Doitrand et al. (2020), which provided the critical GSIF variation at the bimaterial corner as a function of the specimen width. The bimaterial interface corner crack configuration was further studied by Doitrand et al. (2018) based on FFE rather than MA approach. The possible initiation crack surfaces were also determined based on the tensile stress isocontours which were directly integrated in the specimen geometry and mesh at the bimaterial interface (Figure 7). The same

<figure><img src="image_8.png" /><figcaption>**Figure 7:** a) Normal stress distribution at a bimaterial interface (Labossiere et al. 2001) and mesh of the bimaterial interface including crack fronts corresponding to stress isocontours b) based on the stress isocontour topology or c) following the random mesh element edges.</figcaption></figure>

approach of defining the possible crack surfaces based on the stress isocontours and integrate the corresponding lines in the specimen geometry was further used to assess 3D crack initiation in several configurations such as interface scarf joint failure (Doitrand et al. 2018b), mode I+III crack front segmentation (Doitrand et al. 2018c), failure in rhombus hole specimens under compression (Doitrand et al. 2019b), mixed-mode interface failure in ceramic composites (Doitrand et al. 2020b), failure occurring from a pore in bonded joints (Carrère et al. 2021), small-scale specimen failure (Doitrand et al. 2020c), or strut failure in open ceramic foams (Sevecek et al. 2019). Another approach for crack front definition consists in unbuttoning the nodes contained in the area defined by the stress isocontours. It is thus easier to set-up since it does not require that the mesh elements follow the isocontour topology. An example of such crack front definition is shown (Figure 7), which evidences that it results in non-smooth crack fronts following the element edges. As a consequence, a relatively small mesh size may be required to provide similar results as the isoline geometry definition. A comparison between both approaches is still awaited in order to define the conditions of crack shape definitions in this case.

In the above-mentioned works, the crack plane is either assumed or determined based on the stress criterion. An improvement was recently proposed in order to determine the crack plane based on both stress and energy requirements (Doitrand et al. 2023a).

It is not straightforward to validate experimentally the stress isocontours-based 3D crack shapes. Actually, in many configurations, either crack initiation results in the final failure of the specimen or the initiated crack instantaneously propagates to a larger surface (See Section 2.2). The arrest crack surfaces were compared to crack shapes observed experimentally in a few works. Similar curved crack front (Doitrand et al. 2019b) or elliptical crack shape (Doitrand et al. 2018c) as observed experimentally were obtained based on the stress isocontours. Further validation for other configurations is awaited.

### 3.1.3 2D-3D confrontation

In some configurations, it is clearly not possible to simplify the geometry and loading to 2D so that the 3D CC implementation is essential (Yosibash et al. 2016; Doitrand et al. 2018c; Doitrand et al. 2020b). Otherwise, one may question the necessity to implement the CC in 3D compared to the 2D case since it requires extra computational efforts.

Several authors compared 2D and 3D applications of the CC in cases where the 3D crack
rather driven by the stress criterion for thick laminates and by the energy criterion for thin laminates.

Transverse crack initiation was also studied in woven composites by Doitrand, Fagiano, Carrère, et al. (2017) and Doitrand, Fagiano, Hild, et al. (2017). A tentative to base the possible crack surfaces on stress isocontours was not successful since crack initiation was found to be controlled only by the energy criterion. As a consequence, idealized crack surfaces spanning the whole ply thickness and having straight crack fronts were assumed. Debonding ahead of the crack tip was also modeled assuming straight debonding fronts. The stress concentration at specimen free faces was disregarded due to the use of a periodic representative unit cell. The proposed approach enabled modeling the initiation and propagation of multiple transverse cracks and debondings considering interactions between cracks and determining crack and debonding density variations as a function of the applied loading.

Yosibash and Mittelman (2016) applied the CC in 3D to assess crack initiation at a V-notch under mode I+II+III loading. They first determined the crack plane based on the stress criterion only and noticed that crack shape determination in such a plane based on the stress criterion would require different crack shapes depending on the crack plane orientation. They thus assumed circular sector shaped crack whatever the crack orientation, which means that the stress criterion was not homogeneously met on the whole crack surface prior to initiation. As a consequence, an average stress criterion over the circular sector shaped crack was used to implement the CC. Note that the classical stress criterion could also be applied in such a situation: it means that the minimum stress over the considered crack shape before initiation must be larger or equal than the material tensile strength.

Papšík, Ševeček, Hofer, et al. (2023) assumed elliptical crack shape to simulate tunneling crack in layered ceramics containing thermal residual stresses. Exploiting the configuration symmetry, the crack was modeled by releasing Dirichlet boundary conditions on the symmetry plane. They showed that the thickness size effect is governed by both stress and energy criteria. Moreover, they highlighted that crack initiation may be prevented if the energy criterion is not fulfilled, even if the stress reaches a 99% probability of failure according to Weibull's analysis. They thus determined the conditions in terms of layer thickness and residual stresses for which edge crack or tunneling cracks are likely to initiate.

### 3.1.2 Stress isocontour-based 3D crack path

A natural way to determine the possible crack shapes consists in ensuring that they strictly fulfill the stress criterion, i.e. using shapes based on the stress criterion isocontours. It was first proposed by Leguillon (2014) who extended the CC to 3D interface corner crack prediction at a bimaterial junction. The 3D extension was established based on the 3D singularity theory (Leguillon and Sanchez-Palencia 1987; Leguillon 1995; Yosibash 1997b; Yosibash 1997a) and matched asymptotic expansions. Despite some differences in its establishment compared to the 2D formulation, the 3D CC formulation was found to finally be similar to the 2D case. The crack plane was a priori known since it corresponded to the bimaterial interface. An additional hypothesis on the crack geometry was needed compared to the 2D case, consisting in defining the possible crack surfaces based on the tensile stress isocontours which were slightly modified to consider bounded crack surfaces. Note that the MA approach of the CC yields the critical GSIF at initiation. The relation between the GSIF at a singular point in a real specimen and the applied loading can be obtained by implementing a contour integral around the singular point. This approach was implemented by Doitrand, Leguillon, and Martin (2020), which provided the critical GSIF variation at the bimaterial corner as a function of the specimen width. The bimaterial interface corner crack configuration was further studied by Doitrand and Leguillon (2018a) based on FFE rather than MA approach. The possible initiation crack surfaces were also determined based on the tensile stress isocontours which were directly integrated in the specimen geometry and mesh at the bimaterial interface (Figure 7). The same approach of defining the possible crack surfaces based on the stress isocontours and integrate the corresponding lines in the specimen geometry was further used to assess 3D crack initiation in several configurations such as interface scarf joint failure (Doitrand and Leguillon 2018b), mode I+III crack front segmentation (Doitrand
surfaces can be obtained by extrusion of the 2D configuration (García et al. 2016; Doitrand et al. 2018b; Carrère et al. 2021). García et al. (2016) showed that the ply thickness influence on transverse crack initiation loading in composite laminates obtained in 3D was quantitatively similar as in 2D. The maximum difference between the critical stresses at crack initiation obtained in 2D and 3D was 10%. Doitrand et al. (2018) obtained similar variations of the tensile stress or bending moment at crack initiation in scarf joint specimens, the 3D predictions being slightly closer to experiments. They also estimated the influence of using either 2D or 3D CC in order to determine the fracture parameters by inverse confrontation to experimental measurements. They obtained similar tensile strengths, differences smaller than 10% on the critical ERR and a larger difference on the identified shear strength. Carrère et al. (2021) obtained relatively similar 2D and 3D normalized IERR and stress variations as a function of the crack extension near an adhesive single pore and concluded that in this configuration, similar 2D or 3D initiation crack extensions were obtained.

(Doitrand et al. 2019b) compared 2D and 3D crack initiation at rhombus hole in specimens loaded under compression. In such a configuration, the tensile stress isocontours were curved so that the crack extension is larger in the specimen middle plane than at the specimen free edges. Different forces at crack initiation were obtained, they were smaller by up to 15% in 3D than in 2D. Nevertheless, the 2D CC application under plane strain assumption provided a good estimate of the 3D initiation and arrest crack extension in the specimen middle plane.

## 3.2 Dynamic extensions of the CC

The CC (and more generally FFM (Nairn et al. 1993; Hashin 1996)) was originally conceived (Leguillon 2002), and afterwards mainly used in a quasi-static framework, assuming that crack nucleation occurs instantaneously thus disregarding the dynamic formation of the newly created crack. Indeed, Hashin (1996) explicitly dismissed the dynamic aspects stating that “New cracks appear in a very short time and it is not possible or of interest to follow the history of their development” when developing the FFM framework. Under quasi-static loading, the kinetic energy is zero before crack initiation, which implies a positive variation of the kinetic energy in the energy balance (Equation (8)), $\Delta W_k \ge 0$. As a consequence, in a quasi-static framework, the energy condition no longer writes as an equality (Equation (9)) but rather as an inequality: $G_{inc} \ge G_c$. A quasi-static framework makes the CC implementation easier since it avoids calculating the kinetic energy variation during crack initiation. Despite this simplification, the CC is able to predict that crack nucleation occurs instantaneously over a given length at a given loading level similarly to experimental observations (Doitrand et al. 2019b). Extending the CC to dynamic crack nucleation was mentioned as a challenge in the review paper by Weißgraeber et al. (2016), which was addressed by a few works.

A first idea to include loading rate effects in the CC was proposed by Le Pavic et al. (2020) who modified the stress criterion by considering loading rate dependent tensile and shear strengths in order to study tubular bonded structure failure at several loading rates. The advantage of this method is that there is no modification in the IERR and the only difference compared to the classical approach consists in identifying the loading rate dependence of the shear and tensile strengths, which was performed based on Arcan tests. A reasonable agreement of the CC predictions with experimental measurement of tubular bonded sample failure load was obtained based on this approach.

A second study about dynamic crack nucleation was brought by Laschuetza et al. (2021) who investigated dynamic cohesive fracture of a plate with hole under static pre-stress compared to quasi-static FFM results. They pointed out that the quasi-static CC approach predicts an excess of ERR to critical ERR ratio at crack initiation: $G_c > G_{inc}$, possibly resulting in underestimating the failure load. By assuming that this energy excess actually corresponds to the neglected inertial effects in the quasi-static energy balance and based on Freund’s analysis of dynamic crack propagation (Freund 1998), they estimated that the dynamic ERR ($G_{dyn} \approx (1 - \dot{l}/c_R)G_c$, where $\dot{l}$ is the crack velocity and $c_R$ the Rayleigh wave speed) should actually be equal to $G_c$ at initiation. They were thus able to estimate that the crack tip velocity during the finite initiation crack increment increased with increasing hole size, which was consistent with crack velocities measured using CZM.
<figure><img src="image_10.png" /><figcaption>**Figure 7** Meshing strategies to include the crack topology: (a) Normal stress distribution at a bimaterial interface (Labossiere and L. Dunn 2001) and mesh of the bimaterial interface including crack fronts corresponding to stress isocontours (b) based on the stress isocontour topology or (c) following the random mesh element edges.</figcaption></figure>

and Leguillon 2018c), failure in rhombus hole specimens under compression (Doitrand et al. 2019b), mixed-mode interface failure in ceramic composites (Doitrand, Henry, Saad, et al. 2020), failure occurring from a pore in bonded joints (Carrère et al. 2021), small-scale specimen failure (Doitrand, Henry, Zacharie-Aubrun, et al. 2020), or strut failure in open ceramic foams (Ševeček et al. 2019). Another approach for crack front definition consists in unbuttoning the nodes contained in the area defined by the stress isocontours. It is thus easier to set-up since it does not require that the mesh elements follow the isocontour topology. An example of such crack front definition is shown (Figure 7), which evidences that it results in non-smooth crack fronts following the element edges. As a consequence, a relatively small mesh size may be required to provide similar results as the isoline geometry definition. A comparison between both approaches is still awaited in order to define the conditions of crack shape definitions in this case.

In the above-mentioned works, the crack plane is either assumed or determined based on the stress criterion. An improvement was recently proposed in order to determine the crack plane based on both stress and energy requirements (Doitrand, Leguillon, et al. 2023).

It is not straightforward to validate experimentally the stress isocontours-based 3D crack shapes. Actually, in many configurations, either crack initiation results in the final failure of the specimen or the initiated crack instantaneously propagates to a larger surface (See Section 2.2). The arrest crack surfaces were compared to crack shapes observed experimentally in a few works. Similar curved crack front (Doitrand et al. 2019b) or elliptical crack shape (Doitrand and Leguillon 2018c) as observed experimentally were obtained based on the stress isocontours. Further validation for other configurations is awaited.

### 3.1.3 2D-3D confrontation

In some configurations, it is clearly not possible to simplify the geometry and loading to 2D so that the 3D CC implementation is essential (Yosibash and Mittelman 2016; Doitrand and Leguillon 2018c; Doitrand, Henry, Saad, et al. 2020). Otherwise, one may question the necessity to implement the CC in 3D compared to the 2D case since it requires extra computational efforts.

Several authors compared 2D and 3D applications of the CC in cases where the 3D crack surfaces can be obtained by extrusion of the 2D configuration (García et al. 2016; Doitrand and Leguillon 2018b; Carrère et al. 2021). García et al. (2016) showed that the ply thickness influence on transverse crack initiation loading in composite laminates obtained in 3D was quantitatively similar as in 2D. The maximum difference between the critical stresses at crack initiation obtained in 2D and 3D was 10%. Doitrand and Leguillon (2018b) obtained similar variations of the tensile stress or bending moment at crack initiation in scarf joint specimens, the 3D predictions being slightly closer to experiments. They also estimated the influence of using either 2D or 3D CC in order to determine the fracture parameters by inverse confrontation to experimental measurements. They obtained similar tensile strengths, differences smaller than 10% on the critical ERR and a larger difference on the identified shear strength. Carrère et al.
Chao Correas et al. (2022) developed a CC formulation able to account for failure load dependency on the loading rate under dynamic loading regimes. It consists in introducing a characteristic time, τ, which is an empirical, material dependent, input parameter interpreted by the authors as the minimum time required for microcracks to coalesce into a single macrocrack. The proposed CC formulation consists in determining the minimum time for which the stress and the energy criterion are fulfilled for an imposed loading level corresponding to the average imposed loading during the period τ. The authors successfully applied the proposed formulation to rock failure in specimens exhibiting either stress concentration or singularity under dynamic loading. After fitting the characteristic time, they were able to predict failure load dependence on the loading rate for both time to fracture larger than or smaller than the identified characteristic time.

A last extension of the CC considering dynamic crack nucleation was addressed by Doitrand et al. (2022). Instead of assuming an instantaneous crack initiation, it considers the kinetic energy creation due to crack nucleation in the energy balance. Compared to the classical quasi-static approach, it thus requires as extra input the material density as well as the initiation crack velocity profile, i.e. the crack length (or velocity) variation as a function of time during initiation. Dynamic calculations were performed to calculate the potential and kinetic energy variations as a function of the crack length and thus the dynamic IERR. Compared to the quasi-static approach, the IERR decreased with increasing mean crack velocity. The velocity profile also had a significant influence on the predicted failure stress (Doitrand et al. 2022; Chen et al. 2023). The dynamic extension of the CC enabled explaining the differences in holed plate failure stress variation as a function of hole size obtained between experimental measurements and numerical predictions (Li et al. 2006; Leite et al. 2021).

## 3.3 Nonlinearities

The CC was originally developed under linear elasticity assumption and small deformation framework. Even if these assumptions make its implementation easier, there is no restriction to consider nonlinearities in the CC. From a general perspective, considering geometrical or material nonlinearities results in a loss of proportionality between the stress and the applied loading on one hand, and between the IERR and the square loading on the other hand (Equation (21)). As a consequence, the nonlinear CC implementation requires several calculations at different imposed loading magnitudes instead of one given magnitude if nonlinearities are disregarded. For each loading magnitude, stress and energy criteria can be calculated to determine whether the CC is fulfilled or not. If material nonlinearities such as, e.g., plasticity or diffuse damage are considered, the energy change due to plasticity or damage increase must be considered in the energy balance (Equation (9)).

### 3.3.1 Material nonlinearities

Leguillon et al. (2017) applied the CC considering a small plastic or damaged zone ahead of a V-notch in quasi-brittle materials based on the MA approach of the CC. A damage model was weakly coupled to the CC in the sense that the damage locally modified the material Young's modulus, strength and toughness and thus influenced crack initiation. The application of the CC to damaged materials was also implemented by Li et al. (2019) who combined damage field determination based on continuum damage mechanics and discontinuous crack initiation based on the CC. In the stress criterion, on one hand, the stress is replaced by a larger effective stress considering the local damage state, which actually reverts to considering a linear decrease of the tensile strength with increasing damage variable. On the other hand, in the energy criterion, the strain energy variation depends not only on the strain state but also on the damage state. Based on the assumption of a crack path perpendicular to the principal stress and using element erosion to represent the discontinuous crack determined by the CC, crack initiation and propagation in damaged materials was thus assessed.

Crack initiation at a circular hole in PMMA specimens considering material nonlinearity was studied experimentally (Leite et al. 2021) and using the CC (Leite et al. 2021; Chen et al. 2023). It was shown that the fracture stresses decreased with increasing hole diameters. Quasi-static linear elastic (LE) and nonlinear elastic (NLE) applications of the CC were implemented. Minor
(2021) obtained relatively similar 2D and 3D normalized IERR and stress variations as a function of the crack extension near an adhesive single pore and concluded that in this configuration, similar 2D or 3D initiation crack extensions were obtained.

Doitrand et al. (2019b) compared 2D and 3D crack initiation at rhombus hole in specimens loaded under compression. In such a configuration, the tensile stress isocontours were curved so that the crack extension is larger in the specimen middle plane than at the specimen free edges. Different forces at crack initiation were obtained, they were smaller by up to 15% in 3D than in 2D. Nevertheless, the 2D CC application under plane strain assumption provided a good estimate of the 3D initiation and arrest crack extension in the specimen middle plane.

## 3.2 Dynamic extensions

The CC (and more generally FFM (Nairn et al. 1993; Hashin 1996)) was originally conceived (Laguillon 2002), and afterwards mainly used in a quasi-static framework, assuming that crack nucleation occurs instantaneously thus disregarding the dynamic formation of the newly created crack. Indeed, Hashin (1996) explicitly dismissed the dynamic aspects stating that “New cracks appear in a very short time and it is not possible or of interest to follow the history of their development” when developing the FFM framework. Under quasi-static loading, the kinetic energy is zero before crack initiation, which implies a positive variation of the kinetic energy in the energy balance (Equation (8)), $\Delta W_k \ge 0$. As a consequence, in a quasi-static framework, the energy condition no longer writes as an equality (Equation (9)) but rather as an inequality: $G_{inc} \ge G_c$. A quasi-static framework makes the CC implementation easier since it avoids calculating the kinetic energy variation during crack initiation. Despite this simplification, the CC is able to predict that crack nucleation occurs instantaneously over a given length at a given loading level similarly to experimental observations (Doitrand et al. 2019b). Extending the CC to dynamic crack nucleation was mentioned as a challenge in the review paper by Weißgraeber, Leguillon, et al. (2016), which was addressed by a few works.

A first idea to include loading rate effects in the CC was proposed by Le Pavic et al. (2020) who modified the stress criterion by considering loading rate dependent tensile and shear strengths in order to study tubular bonded structure failure at several loading rates. The advantage of this method is that there is no modification in the IERR and the only difference compared to the classical approach consists in identifying the loading rate dependence of the shear and tensile strengths, which was performed based on Arcan tests. A reasonable agreement of the CC predictions with experimental measurement of tubular bonded sample failure load was obtained based on this approach.

A second study about dynamic crack nucleation was brought by Laschuetza and Seelig (2021) who investigated dynamic cohesive fracture of a plate with hole under static pre-stress compared to quasi-static FFM results. They pointed out that the quasi-static CC approach predicts an excess of ERR to critical ERR ratio at crack initiation: $G_c > G_{inc}$, possibly resulting in underestimating the failure load. By assuming that this energy excess actually corresponds to the neglected inertial effects in the quasi-static energy balance and based on Freund’s analysis of dynamic crack propagation (Freund 1998), they estimated that the dynamic ERR ($G_{dyn} \approx (1 - \dot{l}/c_R)G$, where $\dot{l}$ is the crack velocity and $c_R$ the Rayleigh wave speed) should actually be equal to $G_c$ at initiation. They were thus able to estimate that the crack tip velocity during the finite initiation crack increment increased with increasing hole size, which was consistent with crack velocities measured using cohesive zone models (CZM).

Chao Correas et al. (2022) developed a CC formulation able to account for failure load dependency on the loading rate under dynamic loading regimes. It consists in introducing a characteristic time, $\tau$, which is an empirical, material dependent, input parameter interpreted by the authors as the minimum time required for microcracks to coalesce into a single macrocrack. The proposed CC formulation consists in determining the minimum time for which the stress and the energy criterion are fulfilled for an imposed loading level corresponding to the average imposed loading during the period $\tau$. The authors successfully applied the proposed formulation to rock failure in specimens exhibiting either stress concentration or singularity under dynamic loading. After fitting the characteristic time, they were able to predict failure load dependence on
differences on the failure stress were obtained using either NLE or LE material behavior, both
resulting in predicted failure stresses smaller than those measured experimentally. These
differences were further reduced by considering the kinetic energy in the energy balance (Chen
et al. 2023). A similar conclusion that either LE or NLE material behavior results in small
differences in terms of force at crack initiation was drawn by Doitrand et al. (2020) in the case of
crack initiation in Brazilian disk specimens.

3.3.2 Geometrical nonlinearities

The required additional effort to consider geometrical nonlinearities in the CC was highlighted by Talmon l'Armée et al. (2017) who studied the IERR of composite single-lap joints using a nonlinear crack opening integral. They showed the necessity to integrate the forces over the displacement to calculate the IERR in case the load reduction shows nonlinear behavior. Neglecting the nonlinear load reduction otherwise results in either underestimating or overestimating the IERR. The influence of the geometrical nonlinearity on the fracture load in the same configuration was studied by Wei et al. (2022). They highlighted that the consideration of the nonlinearity improved the fracture load prediction since IERR obtained from a linear analysis overestimates the IERR obtained considering geometrical nonlinearities, especially regarding the IERR shear mode contribution.

Rosendahl et al. (2019) coupled equivalent average mixed-mode strain and energy criteria to study crack initiation of a structural silicone adhesive. The use of a strain criterion instead of the classical stress criterion proved more robust to address crack initiation under large deformations. They showed that the IERR exhibited a maximum and thus that crack arrest after initiation was well reproduced by the CC as observed experimentally. The maximum in the IERR was due to mode I IERR contribution, whereas mode II IERR contribution was found to be monotonically increasing. The obtained arrest crack lengths were found to be only slightly dependent on the joint geometry and dimensions.

3.3.3 Other approaches

Torabi et al. (2019) applied the CC to predict crack initiation under moderate or large scale yielding regimes. Instead of considering the actual ductile material behavior, they defined an equivalent brittle material having the same fracture toughness, thus artificially increasing its tensile strength. The proposed approach does not change the CC formulation but introduces the material nonlinearities due to the ductile material behavior through an adjustment of the material tensile strength. Yosibash et al. (2022) applied the classical CC approach disregarding nonlinearities induced by the plastic zone ahead of a V-notch in high strength steel alloys. Similar failure forces were predicted using as critical stress parameter either the yield strength or the ultimate tensile strength. They obtained underestimated failure forces compared to experimental measurements, the differences being larger for larger V-notch angles which was explained by a larger plastic zone size for larger V-notch angles. The plastic strain energy variation due to crack initiation was disregarded in this analysis.

## 3.4 Fatigue loadings

A first model to address crack initiation at a V-notch under fatigue cyclic loading in the framework of the CC was established by Murer et al. (2010); Leguillon et al. (2012). It consisted in considering a gradual linear degradation of the toughness along the presupposed crack path at each fatigue cycle (the cycle SIF amplitude being K<sub>Im</sub>) so that the critical ERR varies between G<sub>IC</sub> (the material critical ERR before any loading) and G<sub>Im</sub>, defined as the critical ERR corresponding to failure at K<sub>Im</sub> SIF. The adjustable parameter in the toughness degradation function was identified so that the rate of advance coincides with that of a Paris law in case of a pre-existing long crack. This model predicted an intermittent crack growth which provided an explanation for the striations observed in experiments.

Wang et al. (2016) assessed fatigue cracking at a welded T-joint by considering a coupled stress and energy criterion involving as input parameters the fatigue threshold and the fatigue limit. The stress criterion stated that fatigue failure cannot occur if cyclic stress range along the potential crack path lies below the fatigue limit of the base material. The energy criterion was written in a similar manner as the classical energy criterion except that G_c was replaced
the loading rate for both time to fracture larger than or smaller than the identified characteristic time.

A last extension of the CC considering dynamic crack nucleation was addressed by Doitrand et al. (2022). Instead of assuming an instantaneous crack initiation, it considers the kinetic energy creation due to crack nucleation in the energy balance. Compared to the classical quasi-static approach, it thus requires as extra input the material density as well as the initiation crack velocity profile, i.e. the crack length (or velocity) variation as a function of time during initiation. Dynamic calculations were performed to calculate the potential and kinetic energy variations as a function of the crack length and thus the dynamic IERR. Compared to the quasi-static approach, the IERR decreased with increasing mean crack velocity. The velocity profile also had a significant influence on the predicted failure stress (Doitrand et al. 2022; Chen et al. 2023). The dynamic extension of the CC enabled explaining the differences in holed plate failure stress variation as a function of hole size obtained between experimental measurements and numerical predictions (Li and X.-B. Zhang 2006; Leite et al. 2021).

## 3.3 Nonlinearities

The CC was originally developed under linear elasticity assumption and small deformation framework. Even if these assumptions make its implementation easier, there is no restriction to consider nonlinearities in the CC. From a general perspective, considering geometrical or material nonlinearities results in a loss of proportionality between the stress and the applied loading on one hand, and between the IERR and the square loading on the other hand (Equation (22)). Nonlinearities influence both the stress and the energy conditions. As a consequence, the nonlinear CC implementation requires several calculations at different imposed loading magnitudes instead of one given magnitude if nonlinearities are disregarded. For each loading magnitude, stress and energy criteria can be calculated to determine whether the CC is fulfilled or not. If material nonlinearities such as, e.g., plasticity or diffuse damage are considered, the energy change due to plasticity or damage increase must be considered in the energy balance (Equation (9)).

### 3.3.1 Material nonlinearities

Leguillon and Yosibash (2017) applied the CC considering a small plastic or damaged zone ahead of a V-notch in quasi-brittle materials based on the MA approach of the CC. A damage model was weakly coupled to the CC in the sense that the damage locally modified the material Young's modulus, strength and toughness and thus influenced crack initiation. The application of the CC to damaged materials was also implemented by Li et al. (2019) who combined damage field determination based on continuum damage mechanics and discontinuous crack initiation based on the CC. In the stress criterion, on one hand, the stress is replaced by a larger effective stress considering the local damage state, which actually reverts to considering a linear decrease of the tensile strength with increasing damage variable. On the other hand, in the energy criterion, the strain energy variation depends not only on the strain state but also on the damage state. Based on the assumption of a crack path perpendicular to the principal stress and using element erosion to represent the discontinuous crack determined by the CC, crack initiation and propagation in damaged materials was thus assessed.

Pletz and Arbeiter (2022) considered an elastic-plastic material behavior with isotropic hardening to simulate fracture in a multi-layered polymer. They were able to consider both crack initiation, arrest and re-initiation as well as permanent deformation of the layers. Crack initiation at a circular hole in PMMA specimens considering material nonlinearity was studied experimentally (Leite et al. 2021) and using the CC (Leite et al. 2021; Chen et al. 2023). It was shown that the nonlinear stress concentration factor is lower than the linear one, so that the fracture stresses decreased with increasing hole diameters but less than under linear elastic conditions. Quasi-static linear elastic (LE) and nonlinear elastic (NLE) applications of the CC were implemented. Minor differences on the failure stress were obtained using either NLE or LE material behavior, both resulting in predicted failure stresses smaller than those measured experimentally. These differences were further reduced by considering the kinetic energy in the energy balance (Chen et al. 2023). A similar conclusion that either LE or NLE material behavior
by a fatigue threshold $\Delta G_{\text{th}}$. Based on this approach, they concluded that residual stresses on fatigue behavior of welded T-joints residual stresses have noticeable effect on the crack deflection and on fatigue strength only for a certain range of stress ratios. The same formulation as in (Wang et al. 2016) was adopted by Liu et al. (2020) to predict the fatigue limit of plates with a central circular hole under mode I loading in comparison with fatigue limits predicted using TCD. They determined unified correlation between the critical distance of TCD and the finite crack increment of the CC.

Sapora et al. (2019); Sapora et al. (2020) studied the fatigue limit of specimens with a center-through thickness sharp crack and a circular notch. Similarly to Wang et al. (2016); Liu et al. (2020), the proposed approach involves two parameters: the plain specimen fatigue limit and the threshold value of the SIF range for fatigue crack growth. Based on this approach, the authors determined that the plain fatigue limit does not depend on the shape of a feature provided its size is sufficiently small, and even becomes feature insensitive for a smaller feature size. The same approach was further extended to predict the fatigue limit of specimens with V- or U-notches (Sapora et al. 2021a) and then extended to fatigue failure under mode III loading (Campagnolo et al. 2021).

# 4 Comparison with other fracture approaches

## 4.1 Cohesive zone models

The main difference between CZM and CC concerns the description of the cracking process. The CC mainly relies on a binary description of fracture considering two possible states, namely undamaged or cracked material. The CZM defines an intermediate state: the process zone, through the description of a traction-separation profile. No separation of the surfaces occurs while the traction remains smaller than a critical stress (the tensile strength if only opening is considered), then the traction decreases with increasing separation between the two surfaces until a critical separation is attained, corresponding to the nucleation of a crack locally. The area under the traction-separation curve corresponds to the material critical ERR. Apart from the traction-separation profile, CZM and the CC share similar input parameters, namely tensile and shear strengths, and opening and shear critical ERR. Another main difference between both approaches is that the stress can locally be larger than the material tensile strength in the CC, not in CZM for which the stress variation at the crack tip is bounded by the material tensile or shear strength. Several authors compared the CC and CZM ability to assess crack initiation. Overall, both methods generally yield similar trends in terms of remote loading at crack initiation variation.

### 4.1.1 Comparison between CC and CZM

Dimitri et al. (2017) obtained differences smaller than a few unit percents on the maximum loads in double cantilever beam specimens predicted by CZM compared to the CC. They highlighted that both approaches provided similar trends and were able to capture the transition between the strength-governed and the fracture energy-ruled debonding process, which was also evidenced for instance by Gentieu et al. (2018); Chao Correas et al. (2021). García et al. (2014); Cornetti et al. (2015); Stein et al. (2015); Távara et al. (2016); Muñoz-Reja et al. (2020) also obtained relatively similar maximum loads using the CC compared to CZM to assess interface crack initiation. Cornetti et al. (2016) obtained a relative error smaller than 5% between the failure stress estimated by the CC and CZM for short crack propagation or crack initiation at a V-notch. Good agreement between V-notch initiation GSIF obtained using either Dudgale's CZM and the CC was also obtained by Murer et al. (2010). Martin et al. (2016) assessed edge debonding initiation at the interface of a bimaterial specimen loaded under four-point bending using either the CC or CZM with bilinear traction-separation profile. They studied the conditions for which CC and CZM provide similar force-displacement curves depending on the interface characteristic fracture length $E_s G_c / \sigma_c^2$, where $E_s$ is the substrate Young's modulus. A good agreement between both methods was only obtained for fracture lengths sufficiently small compared to the specimen characteristic thickness, which ensured a sufficiently small process zone. CZM produced a linear force displacement curve before initiation for small values of the fracture length and shifted to a
results in small differences in terms of force at crack initiation was drawn by Doitrand and Sapora (2020) in the case of crack initiation in Brazilian disk specimens.

### 3.3.2 Geometrical nonlinearities

The required additional effort to consider geometrical nonlinearities in the CC was highlighted by Talmon l'Armée et al. (2017) who studied the IERR of composite single-lap joints using a nonlinear crack opening integral. They showed the necessity to integrate the forces over the displacement to calculate the IERR in case the load reduction shows nonlinear behavior. Neglecting the nonlinear load reduction otherwise results in either underestimating or overestimating the IERR. The influence of the geometrical nonlinearity on the fracture load in the same configuration was studied by Wei et al. (2022). They highlighted that the consideration of the nonlinearity improved the fracture load prediction since IERR obtained from a linear analysis overestimates the IERR obtained considering geometrical nonlinearities, especially regarding the IERR shear mode contribution.

Rosendahl et al. (2019) coupled equivalent average mixed-mode strain and energy criteria to study crack initiation of a structural silicone adhesive. The use of a strain criterion instead of the classical stress criterion proved more robust to address crack initiation under large deformations. They showed that the IERR exhibited a maximum and thus that crack arrest after initiation was well reproduced by the CC as observed experimentally. The maximum in the IERR was due to mode I IERR contribution, whereas mode II IERR contribution was found to be monotonically increasing. The obtained arrest crack lengths were found to be only slightly dependent on the joint geometry and dimensions. Cavitation and crack nucleation in thin hyperelastic adhesives was assessed by Rheinschmidt et al. (2024), who showed that the onset of cavitation cannot be considered as a material failure criterion, otherwise greatly reducing the regime of applicability of this adhesive.

### 3.3.3 Other approaches

Torabi et al. (2019) applied the CC to predict crack initiation under moderate or large scale yielding regimes. Instead of considering the actual ductile material behavior, they defined an equivalent brittle material having the same fracture toughness, thus artificially increasing its tensile strength. The proposed approach does not change the CC formulation but introduces the material nonlinearities due to the ductile material behavior through an adjustment of the material tensile strength. Yosibash et al. (2022) applied the classical CC approach disregarding nonlinearities induced by the plastic zone ahead of a V-notch in high strength steel alloys. Similar failure forces were predicted using as critical stress parameter either the yield strength or the ultimate tensile strength. They obtained underestimated failure forces compared to experimental measurements, the differences being larger for larger V-notch angles which was explained by a larger plastic zone size for larger V-notch angles. The plastic strain energy variation due to crack initiation was disregarded in this analysis.

## 3.4 Fatigue loadings

A first model to address crack initiation at a V-notch under fatigue cyclic loading in the framework of the CC was established by Murer and Leguillon (2010) and Leguillon and Murer (2012). It consisted in considering a gradual linear degradation of the toughness along the presupposed crack path at each fatigue cycle (the cycle SIF amplitude being $K_{Im}$) so that the critical ERR varies between $G_{Ic}$ (the material critical ERR before any loading) and $G_{Im}$, defined as the critical ERR corresponding to failure at $K_{Im}$ SIF. The adjustable parameter in the toughness degradation function was identified so that the rate of advance coincides with that of a Paris law in case of a pre-existing long crack. This model predicted an intermittent crack growth which provided an explanation for the striations observed in experiments.

Wang et al. (2016) assessed fatigue cracking at a welded T-joint by considering a coupled stress and energy criterion involving as input parameters the fatigue threshold and the fatigue limit. The stress criterion stated that fatigue failure cannot occur if cyclic stress range along the potential crack path lies below the fatigue limit of the base material. The energy criterion
non linear response for larger values of the fracture length. These larger fracture lengths resulted
in larger process zones which introduced softening and postponed crack initiation compared to
the CC. Similar conclusions were drawn by Gentieu et al. (2018) who studied the particle size
effect on particle-matrix debonding. They showed that CC and CZM provided similar debonding
critical load for large particles, i.e. when failure is rather governed by the stress criterion. However,
they highlighted different critical stress variations for FFM (1/√R asymptote, R is the particle
radius) and CZM (1/R asymptote) for decreasing particle size. The difference between both
asymptotic behaviors was explained by the presence of a process zone prior to debonding, which
enabled energy dissipation before crack initiation. Chao Correas et al. (2021) studied the size
effect for brittle materials with spherical cavities by modeling annular crack initiation around a
spherical void using both CC and CZM with Dugdale traction-separation profile. Both models
allowed obtaining the gradual transition from stress-driven extreme solutions, the transition
being obtained within the same void radius range. For small voids, failure was driven by a pure
stress criterion and both models matched but larger differences were obtained for larger voids.
Since a good agreement between the CC and Dudgale CZM was obtained for other cracked
geometries (Cornetti et al. 2016; Cornetti et al. 2019a), this result highlighted that different CZM
traction separation profiles corresponding to different geometries may be used to retrieve the CC
failure load predictions.

4.1.2 Influence and identification of the traction-separation profile

The CZM traction-separation profile has no influence on stable rectilinear crack propagation under Griffith's conditions (Acary et al. 2006). However, differences are observed when studying phenomenon such as crack branching or crack initiation. For instance, Dimitri et al. (2017) showed that the difference between the failure load predicted by CC and CZM differed as the critical separation increased, which was also evidenced by Martin et al. (2016). As a consequence, an equivalence may possibly be determined between the CC and a given traction-separation profile, which was studied by several authors.

For instance, similar results were obtained on one hand for bilinear traction-separation profile
and average stress criterion (Dimitri et al. 2017) and on the other hand for Dugdale profile and
classical stress criterion (Henninger et al. 2007b; Cornetti et al. 2016). Rosendahl et al. (2017)
obtained initiation loading and crack lengths in the same order of magnitude at open-holes under
tensile and in-plane bending loading using CC and CZM. They showed that the CC with classical
stress criterion resembles either bilinear (for tensile dominated loading) or trilinear (for bending
dominated loading) traction-separation profile depending on the applied loading and cracking
mechanism. A better agreement with bilinear model was obtained using the average stress
criterion. Doitrand et al. (2019) deepened the CZM traction-separation profile identification
with respect to CC by comparing not only the initiation loading but also the crack arrest length
after initiation at a rhombus hole. The range of initiation forces and crack lengths obtained
with various traction-separation profiles respectively comprised the initiation force and lower
bound for crack arrest obtained using the CC. Similar initiation forces and crack arrest lengths as
those predicted using the CC were obtained using a bilinear cohesive zone model. Cornetti et al.
(2016) observed that for an initial crack in infinite medium or at a V-notch, the CZM process
zone size and initiation length obtained using the CC were significantly different but followed
almost identical trends with respect to the normalized initial crack length. They proposed to
redefine the CZM crack as a zone where the crack opening is larger than a fraction of the critical
opening δc (0.35δc for the studied configuration). This approach enabled matching the initiation
crack length predicted using the CC, except for vanishing initial cracks. This analysis was latter
refined by the same authors (Cornetti et al. 2019a) who introduced a weight function in the stress
condition of the CC which can be tuned or identified in order to match the CZM. The CZM
exhibiting cohesive laws with earlier softening were found to satisfactorily match with CC stress
conditions modified by weight functions higher close to the crack tip and vice-versa. In summary,
it appears that there is not a unique CZM traction-separation profile that enables retrieving the
failure loading and crack length predicted using the CC. The CZM traction-separation profile
corresponding to the CC actually depends on the geometry, the type of loading, the cracking
mechanism and thus has to be identified for a given configuration.
was written in a similar manner as the classical energy criterion except that $G_c$ was replaced by a fatigue threshold $\Delta G_{th}$. Based on this approach, they concluded that residual stresses on fatigue behavior of welded T-joints residual stresses have noticeable effect on the crack deflection and on fatigue strength only for a certain range of stress ratios. The same formulation as in (Wang et al. 2016) was adopted by Liu et al. (2020) to predict the fatigue limit of plates with a central circular hole under mode I loading in comparison with fatigue limits predicted using TCD. They determined unified correlation between the critical distance of TCD and the finite crack increment of the CC.

Sapora et al. (2019) and Sapora et al. (2020) studied the fatigue limit of specimens with a center-through thickness sharp crack and a circular notch. Similarly to Wang et al. (2016) and Liu et al. (2020), the proposed approach involves two parameters: the plain specimen fatigue limit and the threshold value of the SIF range for fatigue crack growth. Based on this approach, the authors determined that the plain fatigue limit does not depend on the shape of a feature provided its size is sufficiently small, and even becomes feature insensitive for a smaller feature size. The same approach was further extended to predict the fatigue limit of specimens with V- or U-notches (Sapora, Cornetti, Campagnolo, et al. 2021) and then extended to fatigue failure under mode III loading (Campagnolo and Sapora 2021). Going beyond the fatigue limit prediction, Mirzaei et al. (2023) proposed a criterion for lifetime estimation of notched components under uniaxial cycling loading, with the capability of considering the influence of both the notch radius and the loading conditions. This criterion was further extended in the case of notched laminated composites (Mirzaei et al. 2024).

# 4 Comparison with other fracture approaches

## 4.1 Cohesive zone models

The main difference between CZM and the CC concerns the description of the cracking process. The CC mainly relies on a binary description of fracture considering two possible states, namely undamaged or cracked material. The CZM defines an intermediate state: the process zone, through the description of a traction-separation profile. No significant separation of the surfaces occurs while the traction remains smaller than a critical stress (the tensile strength if only opening is considered), then the traction decreases with increasing separation between the two surfaces until a critical separation is attained, corresponding to the nucleation of a crack locally. The area under the traction-separation curve corresponds to the material critical ERR. Apart from the traction-separation profile, CZM and the CC share similar input parameters, namely tensile and shear strengths, and opening and shear critical ERR. Another main difference between both approaches is that the stress can locally be larger than the material tensile strength in the CC, not in CZM for which the stress variation at the crack tip is bounded by the material tensile or shear strength. Several authors compared the CC and CZM ability to assess crack initiation. Overall, both methods generally yield similar trends in terms of remote loading at crack initiation variation.

### 4.1.1 Comparison between CC and CZM

Dimitri et al. (2017) obtained differences smaller than a few unit percents on the maximum loads in double cantilever beam specimens predicted by CZM compared to the CC. They highlighted that both approaches provided similar trends and were able to capture the transition between the strength-governed and the fracture energy-ruled debonding process, which was also evidenced for instance by Gentieu et al. (2018) and Chao Correas et al. (2021). García et al. (2014), Cornetti et al. (2015), Stein et al. (2015), Távara et al. (2016), and Muñoz-Reja, Cornetti, et al. (2020) also obtained relatively similar maximum loads using the CC compared to CZM to assess interface crack initiation. Cornetti et al. (2016) obtained a relative error smaller than 5% between the failure stress estimated by the CC and CZM for short crack propagation or crack initiation at a V-notch. Good agreement between V-notch initiation GSIF obtained using either Dudgale's CZM and the CC was also obtained by Murer and Leguillon (2010). Martin et al. (2016) assessed edge debonding initiation at the interface of a bimaterial specimen loaded under four-point
### 4.1.3 Combined used of CC and CZM

Since it is possible to determine a CZM traction-separation profile that yields similar results as the CC in terms of crack initiation remote loading, CC and CZM can be used in a complementary manner to take advantage of both methods. The CC only requires solving one or some (generally linear elastic) calculations whereas CZM leads to nonlinear problems that may increase the overall computational time. Moreover, the fracture parameters $G_c$ and $\sigma_c$ are only implemented during the post-processing phase for the CC, contrary to CZM for which they are a priori defined. It means that once the calculations are performed for the CC, the CC solution can be directly obtained for several values of $G_c$ and $\sigma_c$ without any supplementary FE calculations. The CC is thus computationally more efficient than other methods for which the fracture parameters are a priori defined, e.g. in the context of parameter inverse identification. Doitrand et al. (2020) combined the CC and CZM for $UO_2$ micron scale specimen fracture assessment. Inverse identification of fracture properties was performed using the CC for fresh or irradiated $UO_2$ along three crystal orientations, taking advantage of the fact that $G_c$ and $\sigma_c$ are only implemented during the post-processing phase. Then, CZM was used to assess fracture in similar crystal orientations considering the presence of either one large pore or a distribution of small pores. Ricardo et al. (2020) predicted surface cracking in ceramics due to quenching using CC or CZM. They obtained similar crack initiation loads and spacings, however the CC appeared to be computationally more efficient than CZM, whereas the latter enabled a precise prediction of the crack network density doubling, which was not straightforward to obtain using the CC.

## 4.2 Phase-field models for fracture

Phase-field (PF) fracture approach is a common method to assess fracture problems. Originally developed as an implementation of Griffith's problem revisited as a global minimization problem (Francfort et al. 1998; Bourdin et al. 2000b), it has improved to account for plasticity (Ambati et al. 2015), dynamic effects (Molnár et al. 2020b), fatigue (Lo et al. 2019), interfacial damages (Nguyen et al. 2019) or phase transformation-induced fracture (Djeumen et al. 2022). The main idea of PF approach consists in approximating the sharp crack discontinuity by a smeared damage field description through the use of an internal length ($l_{PF}$). This length controls the magnitude of the damage diffusion. Initially, $l_{PF}$ was used to facilitate the numerical solution and avoid mesh dependence of the crack path, with the idea of reducing $l_{PF}$ to 0 to retrieve the original Griffith theory, thus simply considering $l_{PF}$ as a numerical parameter without physical meaning. Some authors considered this parameter as a material internal length that must be identified, for instance based on the failure loading measured experimentally (Pham et al. 2010a; Freddi et al. 2010) or based on fracture surface observations (Nguyen et al. 2016a).

Similarly to CZM, the CC and PF approaches can model the same fracture problems. Based on the main idea proposed by Leguillon (2002), i.e. considering not only the material critical ERR but also its strength, Kumar et al. (2020) proposed to revisit nucleation in the PF approach by including a stress criterion. Reinoso et al. (2017) simulated notched thin ply laminate fracture using CC and PF approaches. They highlighted the ability of the CC to capture the size effect and provide accurate failure stress compared to experimental data. PF approach provided failure stresses only slightly underestimating experimental measurements. Strobl et al. (2017); Strobl et al. (2020) used CC and PF to simulate Hertzian indentation induced fracture (also studied by Papšík et al. (2023)), resulting in similar trends in terms of crack position with respect to the indenter and critical indenter displacement for fracture initiation. Abaza et al. (2022) obtained similar variations of the apparent SIF at crack nucleation in notched ceramic specimens using CC and PF for different notch geometries. Jiménez-Alfaro et al. (2022) showed that, in the case of small-scale specimen, the critical displacement or force are mainly dependent on the critical ERR. They used the CC as a first step to determine the load range to use in the PF model. They also highlighted that the PF approach could be used as a first step of the CC in case the prescribed crack path is not known a priori.

CC and PF approaches share, as input parameter, the critical ERR, the difference being that the tensile strength in the CC is replaced by the regularization length in PF models. Reinoso et al. (2017) determined the PF regularization length in such a way that the fully degraded states around the stress concentrator were obtained for the final imposed displacement measured
bending using either the CC or CZM with bilinear traction-separation profile. They studied the conditions for which the CC and CZM provide similar force-displacement curves depending on the interface characteristic fracture length $E_sG_c/\sigma_c^2$, where $E_s$ is the substrate Young's modulus. An agreement between both methods was only obtained for fracture lengths sufficiently small compared to the specimen characteristic thickness, which ensured a sufficiently small process zone. CZM produced a linear force displacement curve before initiation for small values of the fracture length and shifted to a non linear response for larger values of the fracture length. These larger fracture lengths resulted in larger process zones which introduced softening and postponed crack initiation compared to the CC. Similar conclusions were drawn by Gentieu et al. (2018) who studied the particle size effect on particle-matrix debonding. They showed that CC and CZM provided similar debonding critical load for large particles, i.e. when failure is rather governed by the stress criterion. However, they highlighted different critical stress variations for FFM ($1/\sqrt{R}$ asymptote, R is the particle radius) and CZM ($1/R$ asymptote) for decreasing particle size. The difference between both asymptotic behaviors was explained by the presence of a process zone prior to debonding, which enabled energy dissipation before crack initiation. Chao Correas et al. (2021) studied the size effect for brittle materials with spherical cavities by modeling annular crack initiation around a spherical void using both CC and CZM with Dugdale traction-separation profile. Both models allowed obtaining the gradual transition from stress-driven extreme solutions, the transition being obtained within the same void radius range. For small voids, failure was driven by a pure stress criterion and both models matched but larger differences were obtained for larger voids. Since a good agreement between the CC and Dugdale CZM was obtained for other cracked geometries (Cornetti et al. 2016; Cornetti et al. 2019), this result highlighted that different CZM traction separation profiles corresponding to different geometries may be used to retrieve the CC failure load predictions.

### 4.1.2 Influence and identification of the traction-separation profile

The CZM traction-separation profile has no influence on stable rectilinear crack propagation under Griffith's conditions (Acary and Monerie 2006). However, differences are observed when studying phenomenon such as crack branching or crack initiation. For instance, Dimitri et al. (2017) showed that the difference between the failure load predicted by CC and CZM differed as the critical separation increased, which was also evidenced by Martin et al. (2016). As a consequence, an equivalence may possibly be determined between the CC and a given traction-separation profile, which was studied by several authors.

For instance, similar results were obtained on one hand for bilinear traction-separation profile and average stress criterion (Dimitri et al. 2017) and on the other hand for Dugdale profile and classical stress criterion (Henninger et al. 2007; Cornetti et al. 2016). Rosendahl et al. (2017) obtained initiation loading and crack lengths in the same order of magnitude at open-holes under tensile and in-plane bending loading using CC and CZM. They showed that the CC with classical stress criterion resembles either bilinear (for tensile dominated loading) or trilinear (for bending dominated loading) traction-separation profile depending on the applied loading and cracking mechanism. A better agreement with bilinear model was obtained using the average stress criterion. Doitrand et al. (2019a) deepened the CZM traction-separation profile identification with respect to CC by comparing not only the initiation loading but also the crack arrest length after initiation at a rhombus hole. The range of initiation forces and crack lengths obtained with various traction-separation profiles respectively comprised the initiation force and lower bound for crack arrest obtained using the CC. Similar initiation forces and crack arrest lengths as those predicted using the CC were obtained using a bilinear cohesive zone model. Cornetti et al. (2016) observed that for an initial crack in infinite medium or at a V-notch, the CZM process zone size and initiation length obtained using the CC were significantly different but followed almost identical trends with respect to the normalized initial crack length. They proposed to redefine the CZM crack as a zone where the crack opening is larger than a fraction of the critical opening $\delta_c$ (0.35$\delta_c$ for the studied configuration). This approach enabled matching the initiation crack length predicted using the CC, except for vanishing initial cracks. This analysis was latter refined by the same authors (Cornetti et al. 2019) who introduced a weight function in the stress
experimentally. Abaza et al. (2022) calibrated the regularization length so that the apparent SIF at crack nucleation were similar to those obtained with the CC. Strobl et al. (2020) used of the PF homogeneous solution under uniaxial tension to determine the regularization length corresponding to a given tensile strength in the PF model.

The PF homogeneous solution consists in neglecting the damage gradient, which enables deriving an analytical relation between the internal length and the material tensile strength (Nguyen et al. 2016a; Strobl et al. 2020; Kumar et al. 2020). It is generally expressed as

$$ \ell_{PF} = \zeta \frac{E \mathcal{G}_c}{\sigma_c^2}. \qquad (26) $$

where $\zeta$ depends on the chosen degradation function (Strobl et al. 2020). Molnár et al. (2020) showed that this relation can actually be generalized to multiaxial loadings, resulting in a surface describing the correlation between $\sigma_c$ and $\ell_{PF}$. This correlation, that depends on the local principal stress ratio and the Poisson's ratio (Molnár et al. 2020a), writes:

$$ \ell_{PF} = \eta(v, \frac{\sigma_{II}}{\sigma_I})^2 \frac{E \mathcal{G}_c}{\sigma_c^2}, \qquad (27) $$

where $\eta$ is a function that can be understood as a normalized tensile strength. It allows determining, for given material properties and local stress state, the regularization length that yields the same fracture stress as in the CC (Figure 8).

<figure><img src="image_9.png" /><figcaption>Figure 8: a) Variation of the normalized tensile strength $\eta$ as a function of the principal stress ratio and b) correlation between $\ell_{PF}$ and $\sigma_c$ for different principal stress ratios obtained for $v = 0.3$.</figcaption></figure>

The established correlation between $\ell_{PF}$ and $\sigma_c$ was supported by confrontations of the CC and PF in other configurations than the homogeneous solution such as mode I crack propagation, shear fracture and crack arrest configuration (Molnár et al. 2020a). Contrary to CZM in which the stress is bounded by $\sigma_c$, it is worth noting that for a given regularization length, stress levels larger than tensile strength determined using Equation (27) are locally attained, similarly to the CC. It was thus highlighted that the PF actually fulfills both stress and energy criteria as in the CC.

Inspired by the correlation between $\ell_{PF}$ and $\sigma_c$, a length-free (LF) implementation of PF approach to fracture was proposed by Doitrand et al. (2023). The input parameters of the LF-PF model are, similarly to the CC, $\mathcal{G}_c$ and $\sigma_c$. This formulation consists in locally adapting the regularization length based on previously determined correlation between $\ell_{PF}$ and $\sigma_c$. Thus, the LF-PF model overcomes the problem of identifying different internal lengths for different testing configurations for the same material, since there is no need to choose or identify an internal length in this formulation. The LF-PF formulation yields crack initiation results close to that obtained with the CC, still conserving the possibility to deal with subsequent crack propagation. The LF-PF can thus be considered as a PF implementation of the CC.
condition of the CC which can be tuned or identified in order to match the CZM. The CZM exhibiting cohesive laws with earlier softening were found to satisfactorily match with CC stress conditions modified by weight functions higher close to the crack tip and vice-versa. In summary, it appears that there is not a unique CZM traction-separation profile that enables retrieving the failure loading and crack length predicted using the CC. The CZM traction-separation profile corresponding to the CC actually depends on the geometry, the type of loading, the cracking mechanism and thus has to be identified for a given configuration.

### 4.1.3 Combined use of CC and CZM

Since it is possible to determine a CZM traction-separation profile that yields similar results as the CC in terms of crack initiation remote loading, CC and CZM can be used in a complementary manner to take advantage of both methods. The CC only requires solving one or some (generally linear elastic) calculations whereas CZM leads to nonlinear problems that may increase the overall computational time. Moreover, the fracture parameters $G_c$ and $\sigma_c$ are only implemented during the post-processing phase for the CC, contrary to CZM for which they are a priori defined. It means that once the calculations are performed for the CC, the CC solution can be directly obtained for several values of $G_c$ and $\sigma_c$ without any supplementary FE calculations. The CC is thus computationally more efficient than other methods for which the fracture parameters are a priori defined, e.g. in the context of parameter inverse identification. Doitrand, Henry, Zacharie-Aubrun, et al. (2020) combined the CC and CZM for UO$_2$ micron scale specimen fracture assessment. Inverse identification of fracture properties was performed using the CC for fresh or irradiated UO$_2$ along three crystal orientations, taking advantage of the fact that $G_c$ and $\sigma_c$ are only implemented during the post-processing phase. Then, CZM was used to assess fracture in similar crystal orientations considering the presence of either one large pore or a distribution of small pores. Faria Ricardo et al. (2020) predicted surface cracking in ceramics due to quenching using CC or CZM. They obtained similar crack initiation loads and spacings, however the CC appeared to be computationally more efficient than CZM, whereas the latter enabled a precise prediction of the crack network density doubling, which was not straightforward to obtain using the CC.

## 4.2 Phase-field models for fracture

Phase-field (PF) fracture approach is a common method to assess fracture problems. Originally developed as an implementation of Griffith's problem revisited as a global minimization problem (Francfort and Marigo 1998; Bourdin et al. 2000), it has improved to account for plasticity (Ambati et al. 2015), dynamic effects (Molnár, Gravouil, et al. 2020), fatigue (Lo et al. 2019), interfacial damages (T.-T. Nguyen et al. 2019) or phase transformation-induced fracture (Djeumen et al. 2022). The main idea of PF approach consists in approximating the sharp crack discontinuity by a smeared damage field description through the use of an internal length ($\ell_{PF}$). This length controls the magnitude of the damage diffusion. Initially, $\ell_{PF}$ was used to facilitate the numerical solution and avoid mesh dependence of the crack path, with the idea of reducing $\ell_{PF}$ to 0 to retrieve the original Griffith theory, thus simply considering $\ell_{PF}$ as a numerical parameter without physical meaning. Some authors considered this parameter as a material internal length that must be identified, for instance based on the failure loading measured experimentally (K. Pham and Marigo 2010; Freddi and Royer-Carfagni 2010) or based on fracture surface observations (T.-T. Nguyen et al. 2016).

Similarly to CZM, the CC and PF approaches can model the same fracture problems. Based on the main idea proposed by Leguillon (2002), i.e. considering not only the material critical ERR but also its strength, Kumar et al. (2020) proposed to revisit nucleation in the PF approach by including a stress criterion. Reinoso et al. (2017) simulated notched thin ply laminate fracture using CC and PF approaches. They highlighted the ability of the CC to capture the size effect and provide accurate failure stress compared to experimental data. PF approach provided failure stresses only slightly underestimating experimental measurements. Strobl et al. (2017) and Strobl and Seelig (2020) used CC and PF to simulate Hertzian indentation induced fracture (also studied by Papšík, Ševeček, Martin, et al. (2023)), resulting in similar trends in terms of crack position
## 4.3 Other models

In the framework of FFM, approaches based on the TCD are usually employed for engineering failure predictions (Taylor et al. 2005; Taylor 2007). Based on the fact that a maximum stress criterion is not deemed suitable to predict experimentally observed size effect even for non-singular stress fields (Chao Correas et al. 2021), the TCD consists in comparing the stress at a certain distance from a stress concentration or singular point to the material tensile strength. The major difference between the TCD and the CC is that the intrinsic length is an input parameter of the TCD, whereas in the CC it is an output that is obtained by combining stress and energy conditions. Liu et al. (2020) predicted the fatigue limit of V-notch specimens using the TCD and the CC. They derived a relation between the critical crack advance for both approaches. They found a relation that is independent of the geometry and material parameters, which is consistent with previous CC studies highlighting that the initiation crack length in the CC is proportional to Irwin's length (Martin et al. 2018; Doitrand et al. 2020f) and provided that the TCD critical distance is also based on Irwin's length. Whatever the notch radius, the critical distance in the TCD is constant whereas it decreases with increasing notch radius using the CC. The TCD Line method (i.e. evaluating the average stress over a characteristic distance instead of the stress at a certain distance from the stress concentration or singularity) and CC average stress criterion were found approximately equivalent in predicting the fatigue limit of specimens with various notch sizes. Differences up to 10% are obtained between both methods depending on the notch radius. Chao Correas et al. (2021) showed that both the TCD and the CC allowed obtaining the gradual transition from the stress-driven extreme solutions for crack initiation at a spherical void. Both approaches placed the transition within the same range of void radii compared to the material characteristic length. They also showed that the CC provided similar brittle ceramic failure stress variation as a function of the pore size compared to atomistic simulation predictions.

Campagnolo et al. (2016) compared the Strain Energy Density (SED) approach to the CC. SED model consists in considering as critical parameter the strain energy density measured over a control volume around the crack initiation location. The comparison was based on crack initiation at a V-notch under in-plane shear loading. The apparent SIF at crack initiation, computed analytically, were shown to be proportional to powers of $K_{Ic}$ and $\sigma_c$ for both methods. They only differed in the proportionality factor, which is a function of the notch angle in the CC whereas it is a function of the Poisson's ratio in the SED approach. Both approaches predicted similar apparent SIF at crack initiation in this particular configuration. The radius of the control volume, over which the strain energy density is averaged, is determined based on Irwin's length (Yosibash et al. 2004; Campagnolo et al. 2016), which raises the question about the applicability of SED when considering small scale specimens.

Zghal et al. (2018) compared the Thick Level Set (TLS) approach to the matched asymptotic approach of the CC. The TLS model is based on the evolution of the damage field with a dissipation control algorithm, introducing a characteristic length which represents the smallest possible distance between a point where there is no damage and a point fully damaged. Sharp or blunted notches and cavities were considered. TLS results were slightly influenced by the choice of the material model, i.e. the stress decrease as a function of the crack opening. TLS and CC resulted in close apparent strengths for all cases provided the assumptions of the MA approach are satisfied.

Carrère et al. (2015) confronted the CC to Continuum Damage Modeling (CDM) to assess adhesively bonded joint failure. Despite different definitions of the final failure, similar failure loads were obtained using both models provided small displacement assumption is valid since crack initiation occurs just before specimen final failure. Otherwise, the CC was found to provide conservative failure loads compared to CDM.

Sapora et al. (2021) confronted Gradient Elasticity (GE) model and CC to predict borehole crack initiation under uniform biaxial loading and pressure. The CC can be considered as local in the constitutive law and non-local in the failure criterion, since crack initiation occurs when stress and energy conditions are simultaneously met at a certain distance of the stress concentration or singular point. On the opposite, GE model is nonlocal in the constitutive law and local in the failure criterion, since it considers the stress concentration factor as the governing failure parameter. It is based on a characteristic internal length which quantifies the distance over
with respect to the indenter and critical indenter displacement for fracture initiation. Abaza et al. (2022) obtained similar variations of the apparent SIF at crack nucleation in notched ceramic specimens using CC and PF for different notch geometries. Jiménez-Alfaro et al. (2022) showed that, in the case of small-scale specimen, the critical displacement or force are mainly dependent on the critical ERR. They used the CC as a first step to determine the load range to use in the PF model. They also highlighted that the PF approach could be used as a first step of the CC in case the prescribed crack path is not known a priori.

CC and PF approaches share, as input parameter, the critical ERR, the difference being that the tensile strength in the CC is replaced by the regularization length in PF models. Reinoso et al. (2017) determined the PF regularization length in such a way that the fully degraded states around the stress concentrator were obtained for the final imposed displacement measured experimentally. Abaza et al. (2022) calibrated the regularization length so that the apparent SIF at crack nucleation were similar to those obtained with the CC. Strobl and Seelig (2020) used of the PF homogeneous solution under uniaxial tension to determine the regularization length corresponding to a given tensile strength in the PF model.

The PF homogeneous solution consists in neglecting the damage gradient, which enables deriving an analytical relation between the internal length and the material tensile strength (T.-T. Nguyen et al. 2016; Strobl and Seelig 2020; Kumar et al. 2020). It is generally expressed as

$$ \ell_{PF} = \zeta \frac{E \mathcal{G}_c}{\sigma_c^2} \qquad (27) $$

where $\zeta$ depends on the chosen degradation function (Strobl and Seelig 2020). Molnár, Doitrand, et al. (2020) showed that this relation can actually be generalized to multiaxial loadings, resulting in a surface describing the correlation between $\sigma_c$ and $\ell_{PF}$. This correlation, which depends on the local principal stress ratio and the Poisson's ratio (Molnár, Doitrand, et al. 2020), reads

$$ \ell_{PF} = \eta \left( \nu, \frac{\sigma_{II}}{\sigma_I} \right)^2 \frac{E \mathcal{G}_c}{\sigma_c^2} \qquad (28) $$

where $\eta$ is a function that can be understood as a normalized tensile strength. It allows determining, for given material properties and local stress state, the regularization length that yields the same fracture stress as in the CC (Figure 8).

<figure><img src="image_11.png" /><figcaption>**Figure 8** Correlation between the PF regularization length and the material tensile strength: (a) Variation of the normalized tensile strength $\eta$ as a function of the principal stress ratio and (b) correlation between $\ell_{PF}$ and $\sigma_c$ for different principal stress ratios obtained for $\nu = 0.3$.</figcaption></figure>

The established correlation between $\ell_{PF}$ and $\sigma_c$ was supported by confrontations of the CC and PF in other configurations than the homogeneous solution such as mode I crack propagation, shear fracture and crack arrest configuration (Molnár, Doitrand, et al. 2020). Contrary to CZM in which the stress is bounded by $\sigma_c$, it is worth noting that for a given regularization length, stress levels larger than tensile strength determined using Equation (28) are locally attained, similarly to the CC. It was highlighted that the PF actually fulfills both stress and energy criteria as in the CC.
which non-locality acts, thus smoothing high variations in the elastic stress field. It appears
to be not applicable below a threshold size as the internal length becomes comparable to the
hole radius. Sapora et al. (2021) showed that almost identical failure stress predictions could
be obtained provided the internal length in the GE model is properly chosen compared to the
material characteristic length.

Doitrand et al. (2021) assessed fracture size effect based on the CC or Weibull approaches. Whereas both approaches well reproduced the size effect in gypsum specimens, only the CC succeeded in correctly predicting the failure stress variation for ZnO specimens. Fracture surface observations actually revealed that the critical pore size range lied around 50–250 microns whatever the specimen size, not retrieving increasing critical pore size with increasing specimen size (Uhl et al. 2022). It thus questioned the basis assumption of increasing critical flaw size with increasing specimen size associated to Weibull’s approach. The CC and Weibull approach were combined by Leguillon et al. 2015a to predict ceramic bending failure and rationalize the differences between the measured tensile and bending strength, concluding that only the tensile strength can be considered as a material parameter.

5 CC applications

The CC is particularly adapted to assess crack initiation in many configurations already described
by Weißgraeber et al. (2016), such as in presence of pores, holes or cavities (Rosendahl et al. 2016;
Felger et al. 2017a; Sapora et al. 2018b; Sapora et al. 2018a; Torabi et al. 2018; Torabi et al. 2017;
Carrère et al. 2021; Uhl et al. 2022; Chao Correas et al. 2021; Doitrand et al. 2021d), U-notches
(Picard et al. 2006; Leguillon et al. 2007; Doitrand et al. 2020b; Baldassari et al. 2023) or V-notches
(Felger et al. 2017c; Martin et al. 2018; Xu et al. 2019), bonded joints or interfaces (Stein et al.
2016a; Stein et al. 2016b; Muñoz-Reja et al. 2017; Dolling et al. 2019; Felger et al. 2019b; Aranda
et al. 2020; Birro et al. 2020; Birro et al. 2021; Frey et al. 2021; Dolling et al. 2021) or less classical
applications such as snow slab avalanches prediction (Rosendahl et al. 2020). We herein focus on
two particular applications: multi-cracking and small scale failure.

**5.1 Multi-cracking**

The CC may be applied to predict multiple crack initiations, such as multiple transverse cracks in laminates (Kashtalyan et al. 2016) or in a stiff inclusion embedded in a soft matrix under a compressive loading (Quesada et al. 2009). This configuration results in multiple crack initiations in the inclusion due to an excess of available energy at first crack initiation. A main difficulty that arises when handling multiple cracks is related to the number of possible initiation configurations ($2^N$ for $N$ potential crack locations). Actually, the number of possible initiation configurations can be reduced by only considering potential interactions between two cracks that are close to each other. For instance, Doitrand et al. (2017) showed that it is necessary to include the initiation and propagation of both cracks in the energy balance to consider the interaction between two cracks that are close to each other. Nevertheless, they highlighted that the first crack has a negligible influence on the second crack initiation if both cracks are separated by more than 1 mm in a woven composite. Based on this observation, it was possible to reduce the number of possible configurations to assess the initiation and propagation of 14 cracks in a woven composite representative volume element by considering the interactions between only 6 pairs of cracks.

Another strategy was employed by Leguillon et al. (2016) to predict multiple surface cracks
in an oxidized polymer under bending. It consists in modeling a periodic array of cracks by
a representative cell (RC) containing a single crack. As a consequence, the spacing between
two cracks is determined by the RC length and the interaction between two neighbor cracks is
considered through periodic boundary conditions. Crack density increase is then obtained by
crack network subdivision through adding a second crack in the RC. The same approach was also
used to assess multi-cracking in coatings or thin layers (Leguillon et al. 2017a; Leguillon et al.
2018a; Coleman et al. 2020), inclined facet array initiation ahead of a parent crack front under
mode I+III loading (Doitrand et al. 2018c) or thermal shock-induced multi-cracking (Ricardo et al.
2020).
Inspired by the correlation between $\ell_{PF}$ and $\sigma_c$, a length-free (LF) implementation of PF approach to fracture was proposed by Doitrand, Molnár, et al. (2023). The input parameters of the LF-PF model are, similarly to the CC, $G_c$ and $\sigma_c$. This formulation consists in locally adapting the regularization length based on previously determined correlation between $\ell_{PF}$ and $\sigma_c$. Thus, the LF-PF model overcomes the problem of identifying different internal lengths for different testing configurations for the same material, since there is no need to choose or identify an internal length in this formulation. The LF-PF formulation yields crack initiation results close to that obtained with the CC, still conserving the possibility to deal with subsequent crack propagation. The LF-PF can thus be considered as a PF implementation of the CC.

## 4.3 Other models

In the framework of FFM, approaches based on the TCD are usually employed for engineering failure predictions (Taylor et al. 2005; Taylor 2007). Based on the fact that a maximum stress criterion is not deemed suitable to predict experimentally observed size effect even for non-singular stress fields (Chao Correas et al. 2021), the TCD consists in comparing the stress at a certain distance from a stress concentration or singular point to the material tensile strength. The major difference between the TCD and the CC is that the intrinsic length is an input parameter of the TCD, whereas in the CC it is an output that is obtained by combining stress and energy conditions. Liu et al. (2020) predicted the fatigue limit of V-notch specimens using the TCD and the CC. They derived a relation between the critical crack advance for both approaches. They found a relation that is independent of the geometry and material parameters, which is consistent with previous CC studies highlighting that the initiation crack length in the CC is proportional to Irwin's length (Martin et al. 2018; Doitrand, Martin, et al. 2020) and provided that the TCD critical distance is also based on Irwin's length. Whatever the notch radius, the critical distance in the TCD is constant whereas it decreases with increasing notch radius using the CC. The TCD Line method (i.e. evaluating the average stress over a characteristic distance instead of the stress at a certain distance from the stress concentration or singularity) and CC average stress criterion were found approximately equivalent in predicting the fatigue limit of specimens with various notch sizes. Differences up to 10% are obtained between both methods depending on the notch radius. Chao Correas et al. (2021) showed that both the TCD and the CC allowed obtaining the gradual transition from the stress-driven extreme solutions for crack initiation at a spherical void. Both approaches placed the transition within the same range of void radii compared to the material characteristic length. They also showed that the CC provided similar brittle ceramic failure stress variation as a function of the pore size compared to atomistic simulation predictions.

Campagnolo et al. (2016) compared the Strain Energy Density (SED) approach to the CC. SED model consists in considering as critical parameter the strain energy density measured over a control volume around the crack initiation location. The comparison was based on crack initiation at a V-notch under in-plane shear loading. The apparent SIF at crack initiation, computed analytically, were shown to be proportional to powers of $K_{Ic}$ and $\sigma_c$ for both methods. They only differed in the proportionality factor, which is a function of the notch angle in the CC whereas it is a function of the Poisson's ratio in the SED approach. Both approaches predicted similar apparent SIF at crack initiation in this particular configuration. The radius of the control volume, over which the strain energy density is averaged, is determined based on Irwin's length (Yosibash et al. 2004; Campagnolo et al. 2016), which raises the question about the applicability of SED when considering small scale specimens.

Zghal et al. (2018) compared the Thick Level Set (TLS) approach to the matched asymptotic approach of the CC. The TLS model is based on the evolution of the damage field with a dissipation control algorithm, introducing a characteristic length which represents the smallest possible distance between a point where there is no damage and a point fully damaged. Sharp or blunted notches and cavities were considered. TLS results were slightly influenced by the choice of the material model, i.e. the stress decrease as a function of the crack opening. TLS and CC resulted in close apparent strengths for all cases provided the assumptions of the MA approach are satisfied.

Carrère et al. (2015) confronted the CC to Continuum Damage Modeling (CDM) to assess adhesively bonded joint failure. Despite different definitions of the final failure, similar failure
## 5.2 Small scale fracture

There is no limitation to address small scale specimen fracture using the CC. In practice, fracture of small scale specimens may rather be driven by the energy criterion (Leguillon 2002) since the available stored energy in the system depends on the specimen volume. For instance, micron-scale alumina platelet apparent bending strength much larger than its tensile strength was measured experimentally (Feilden et al. 2017), due to failure mainly being driven by the energy criterion, and thus the material critical ERR (Doitrand et al. 2020a). Microscopic scale interface crack initiation in nacre-like composites made of similar alumina platelets embedded in a secondary glassy phase acting as an interface between the platelets was also addressed using the CC (Doitrand et al. 2020b; Duminy et al. 2023). Duminy et al. (2023) highlighted interface crack containment phenomenon that resulted in interface crack initiation being driven by the energy criterion.

Situations where the failure stress exceeds the critical stress were also encountered at the atomic scale by means of molecular simulations of graphene failure initiation (Brochard et al. 2016), the obtained failure behavior was consistent with the one predicted by the CC.

Gallo et al. (2020) implemented the CC to predict crack initiation of single crystal silicon notched cantilever beam tested in a transmission electron microscope. Whereas the CC correctly provided the experimentally measured failure loads, they highlighted a strong deviation of failure loads predicted by LEFM compared to experimentally measured loads due to the small studied scale. Crack initiation in micro-scale cantilever beam was also addressed by Doitrand et al. (2020); Jiménez-Alfaro et al. (2021). Jiménez-Alfaro et al. (2021) highlighted the limitation of the MA approach at small scales when the assumption for crack length smallness with respect to the studied structure is no longer valid. They also showed that small scale specimen fracture is actually weakly sensitive to the material tensile strength.

# 6 Conclusion and perspectives

Among the challenges mentioned in a previous review paper by Weißgraeber et al. (2016), several major advances have been proposed such as the CC application in 3D, extension to non linearities and fatigue, or considering dynamic loading or kinetic energy variation during crack initiation. The CC has also been successfully confronted to other fracture approaches, thus enabling using them in a complementary manner. Some questions are nevertheless still awaiting answers or further analysis.

If the energy balance, from which the energy criterion derives, seems difficult to refute, there are more possibilities for the stress criterion. Beyond the classical discussion about averaging or not the stress, the use of a strain criterion instead of a stress criterion could also be discussed. Although similar results can be expected with respect to stress criterion for linear elastic materials, it could bring significant differences in case of nonlinear material behavior, for instance metals exhibiting plasticity or 3D printed materials.

How to define the crack shapes in configurations where failure is controlled by the energy criterion? Are the stress isocontour-based crack shapes still valid in such configurations? The definition of crack shapes remains a crucial point to further develop 3D CC applications. Further experimental confrontation of the 3D crack shapes based on stress isovalues would be necessary, which requires determining stable crack initiation configurations with no unstable crack propagation after initiation.

Experimental evidences of sudden crack jumps at initiation have been given, which supports the CC results that a crack initiates over a finite length. Nevertheless, in many situations this jump includes both crack initiation and subsequent unstable crack propagation. Even with a rapid camera, it seems difficult to separate both phases experimentally. Therefore, experimental results in stable crack initiation configurations, driven by the energy criterion, are awaited to definitely confirm the existence of a crack jump solely corresponding to crack initiation.

Regarding the crack jump, despite the recent advances to consider the dynamic crack formation resulting in a kinetic energy increase, the dynamic extension of the CC requires the knowledge of the crack velocity profile. A challenge will consist in either determining it experimentally or predicting it numerically, which seems unavoidable to make this approach applicable without adding another parameter to identify.