# Review of the matched asymptotic approach of the coupled criterion

Sara Jiménez-Alfaro, Israel García García, Aurélien Doitrand

► To cite this version:

Sara Jiménez-Alfaro, Israel García García, Aurélien Doitrand. Review of the matched asymptotic approach of the coupled criterion. Comptes Rendus. Mécanique, 2025, 353 (G1), pp.339-357. 10.5802/crmeca.285. hal-04935607

HAL Id: hal-04935607

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

Submitted on 12 Feb 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Mécanique

Sara Jiménez-Alfaro, Israel García García and Aurélien Doitrand

**Review of the matched asymptotic approach of the coupled criterion**

Volume 353 (2025), p. 339-357

Online since: 6 February 2025

https://doi.org/10.5802/crmeca.285

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

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/
Review article / *Article de synthèse*

# Review of the matched asymptotic approach of the coupled criterion

## Revue de l'approche asymptotique du critère couplé

Sara Jiménez-Alfaro<sup>●,*,a,b</sup>, Israel García García<sup>●,c</sup> and Aurélien Doitrand<sup>●,d</sup>

E-mails: sara.jimenezalfaro@eng.ox.ac.uk, s.jimenez-alfaro@imperial.ac.uk (S. Jiménez-Alfaro), israelgarcia@us.es (I. G. García), aurelien.doitrand@insa-lyon.fr (A. Doitrand)

**Abstract.** Matched Asymptotics is a powerful mathematical technique with broad applicability in various engineering fields. One of its key uses is in Fracture Mechanics, where it provides accurate approximations in the vicinity of the crack tip with low computational complexity. This method can be seamlessly integrated with the Coupled Criterion (CC), which enables the prediction of crack nucleation and propagation in brittle materials. Hence, this paper deeply explains how the MA technique can be applied together with the CC in the context of Fracture Mechanics, providing a detailed literature review of the advances made in the last decade.

**Résumé.** Les développements asymptotiques raccordés constituent une technique mathématique puissante, largement applicable dans divers domaines de l'ingénierie. L'une de leurs principales utilisations se situe en mécanique de la rupture, où ils permettent d'obtenir des approximations précises à proximité de la pointe des fissures tout en maintenant une faible complexité de calcul. Cette méthode peut être intégrée de manière fluide au critère couplé (CC), qui permet de prédire l'amorçage et la propagation des fissures dans les matériaux fragiles. Cet article explique comment la technique des développements asymptotiques raccordés peut être utilisée conjointement avec le critère couplé dans le cadre de la mécanique de la rupture, tout en offrant une revue détaillée de la littérature sur les avancées réalisées au cours de la dernière décennie.

**Keywords.** Matched asymptotic expansion, Coupled criterion, Fracture mechanics.

**Mots-clés.** Développement asymptotique raccordé, Critère couplé, Mécanique de la rupture.

**Funding.** European Union's Horizon 2020 research and innovation programme under Marie Sklodowska-Curie grant agreement No. 861061-NEWFRAC, Iberdrola Foundation under the Marie Sklodowska-Curie Grant Agreement No 101034297, Ministerio de Ciencia e Innovación de España (Project PID2020-117001GB-100/AEI/10.13039/5011000111033).

Manuscript received 29 October 2024, revised and accepted 14 January 2025.

<footnote><sup>a</sup> Department of Civil and Environmental Engineering, Imperial College London, Exhibition Road, London, SW7 2AZ, UK</footnote>

<footnote><sup>b</sup> Department of Engineering Science, University of Oxford, Parks Road, Oxford, OX1 3PJ, UK</footnote>

<footnote><sup>c</sup> Departamento de Mecánica de Medios Continuos y Teoría de Estructuras, Escuela Técnica Superior de Ingeniería, Escuela Politécnica Superior, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Sevilla, Spain</footnote>

<footnote><sup>d</sup> Université Lyon, INSA-Lyon, UCBL, CNRS, MATEIS, UMR5510, F-69621 Villeurbanne, France</footnote>

<footnote>*Corresponding author</footnote>
# 1. Introduction

The matched asymptotic (MA) expansion method is an approach that enables solving an equation or a system of equations [1–5]. It is well adapted to solve singularly perturbed differential equations, for which different approximate solutions are determined, each of which being accurate for a given part of the domain under investigation. These solutions are then combined to give a single approximate solution that is accurate for the whole domain under investigation. The domain may generally be divided into two subdomains. In the first one, the solution (called the outer solution) is accurately approximated by an asymptotic series representing a regular perturbation (i.e. by setting to zero a small parameter representing, e.g., a singular perturbation). The second one consists of a region in which this first approximation is inaccurate, due to perturbation terms that are not negligible. This constitutes the inner solution. The outer and inner solutions are then combined through a process called “matching” in such a way that an approximate solution for the whole domain is finally obtained.

Asymptotic expansions were used to define the elastic constitutive law of the homogeneous equivalent material of a composite when a tangential slip is allowed on the fiber/matrix interface [6]. It was shown that a limit slip coefficient exists beyond which the stiffness of the material rapidly decreased. They were also used in the framework of homogenization as an alternative to the multiple scalings approach [7, 8]. MA expansions were used in combination to the singularity theory to determine the elastic displacements and stress fields corresponding to a class of junctions between rods and bulk bodies modeled as a flexible clamping in the framework of two-dimensional elasticity [9]. Leguillon analysed the problem of crack branching in a homogeneous elastic but non isotropic material. Based on asymptotic expansions, the energy release rate was computed and a revisited Griffith's criterion including anisotropic fracture properties was suggested [10]. Sicsic and Marigo studied the propagation of a crack band and derived the conditions for which it behaves like a Griffith's crack [11]. MA expansions were also used to study the behavior of interface cracks, for instance to further analyse the “Cook and Gordon” [12] interface debonding effect ahead of a primary crack [13], edge debonding in laminates [14] or to analyze the role of residual thermal stresses on the crack deflection or penetration at a bimaterial interface. The 2D and 3D singularities at a bimaterial interface were derived [15], also considering contact and friction between two anisotropic materials [16]. The mode III asymptotic expansions for a crack in or along a joint enabled defining an apparent toughness of the interface to be used for crack propagation [17]. It was also used to derive the stress intensity factors near an angular point on the front of an interface crack [18].

Moreover, the character of the stress singularity at the tip of a classical crack in a homogeneous material was approximated by an asymptotic series for cracks in Mode I, Mode II and Mode III. The first two modes were studied in the work of Williams [19], which is well known by the scientific community, since asymptotic solutions for free-free, clamped-clamped and free-clamped boundary conditions are given therein.

MA expansions are particularly relevant when studying fracture and especially crack initiation in a structure. Indeed, the latter can be studied in the framework of MA expansions as the unbroken problem corrected by the crack that initiates (provided its smallness with respect to the structure characteristic dimensions). This idea was actually made effective by Leguillon [20, 21] who proposed to study crack initiation by coupling a stress criterion and an energy criterion. This approach has spread and is now a common way to study crack initiation, as evidenced by numerous applications summarized in the two review papers [22, 23]. The CC can be implemented through several ways, for instance by solving an implicit equation if analytical solutions can be provided for the stress fields and the energy release rate variation as a function of the crack surface [22, 24–26]. A second way to implement the CC is through finite element
<figure><img src="image_3.png" /><figcaption>**Figure 1.** V-notched three-point bending specimen subjected to a force *F* with a crack (length *ℓ*) initiating at the V-notch tip. Notice that *e* represents the notch depth.</figcaption></figure>

(FE) simulations of the full structure under investigation including the crack that initiates [27–30]. In cases where analytical solutions are not available and so as to achieve a numerically more efficient approach than full FE implementation, an effective way to implement the CC is to use MA expansions. The objective of this paper is to give an overview of the matched asymptotic approach of the CC. We first recall the general idea (Section 2) and the formulation (Section 3) of the MA approach. Then, we describe its numerical implementation (Section 4) and provide some application examples (Section 5).

## 2. The idea behind matched asymptotics

Before presenting the mathematical formalism of the matched asymptotic (MA) approach of the Coupled Criterion (CC), this section is dedicated to provide the philosophy behind it to further understand how it can be set up and define the main required ingredients. In the sequel, the CC is formulated under linear elasticity and small deformation assumptions. Both inertial effects [31, 32] or dissipation mechanisms other than cracking that may occur during initiation, such as, e.g., plasticity [33], diffuse damage [34] or viscous effects, are disregarded. The MA approach of the CC is useful to efficiently study the problem of a small crack initiating in a complex structure subjected to a mechanical or thermal loading. The objective is to determine the loading level at which the crack is likely to initiate as well as the initiation crack length. As a matter of example, we consider the problem of a crack of length *ℓ* that initiates at the tip of a V-notch (angle *β*) in a specimen loaded under three-point bending (Figure 1). Notice that this technique is only valid provided *ℓ* is smaller than a characteristic dimensions of specimen (*ℓ* ≪ *e* in Figure 1), an initial assumption that should be checked after the implementation, once the actual value of *ℓ* is obtained using the coupled criterion.

### 2.1. *The coupled criterion*

The main idea behind the CC arises from the following observations:

* Considering an energy criterion only, it enables assessing the propagation of a crack based on the material critical energy release rate $G_c$ [35–37] but generally fails to study its initiation.

* Considering a stress criterion only, it enables assessing crack initiation based on the material tensile strength $\sigma_c$ except in the presence of a singular point.

Stress and energy criteria thus appear as complementary and their combination enables assess-
ing crack initiation in many configurations. The stress criterion of the CC is a condition estab-
lished in the initial domain before crack initiation (thus without crack). It states that the stress
normal to the future crack path must be larger than the material strength attained under a simi-
lar principal stress state. For instance, it reverts to comparing the opening stress to material ten-
sile strength under uniaxial tensile loading. For the sake of simplicity, we will consider a brittle
material that exhibits a Rankine strength surface in the sequel, which enables defining the material strength surface based on a single parameter, its tensile strength. The stress criterion thus requires the calculation of the stress field before crack initiation. In the vicinity of a V-notch, the stress tensor actually writes as an expansion in powers of $r$, William's expansion in this case since it is a singularity, Taylor's expansion for a smooth stress field:

$$ \underline{\underline{\sigma}} = Kr^{\lambda-1}\underline{s}(\theta) + o(r^{\lambda-1}), \qquad (1) $$

where $r$ and $\theta$ are polar coordinates, $\lambda$ is the characteristic exponent of the singularity and $\underline{s}$ is the stress function derived from the characteristic mode of the singularity $\underline{u}$. The characteristic exponent and mode of the singularity are obtained by solving an eigen-value problem [2]. The parameter $K$ is the Generalized Stress Intensity Factor (GSIF) of the singularity, it represents the magnitude of the loading around the V-notch. Notice that in (1) only the dominant term has been represented, assuming that it is real and has multiplicity one, due to the symmetry of the problem represented in Figure 1. However, this is not the case of mixed mode loadings, for example, where two or more singular terms should be considered, with the associated GSIFs.

The energy criterion of the CC is obtained from the energy equilibrium between the states prior to and after crack initiation. The crack surface creation energy $\mathcal{G}_c S$, where $S$ is the crack surface, must be balanced by the variation in external force work ($W_{ext}$) and in elastic strain energy ($W_{el}$) so that:

$$ \Delta W_{el} + \mathcal{G}_c S = \Delta W_{ext} \qquad (2) $$

When solving the CC, the objective is to determine the initiation crack surface $S_c$ and initiation imposed loading (for instance the initiation force $F_c$ based on the example provided in Figure 1) by simultaneously fulfilling both stress and energy criteria. We thus need (i) one calculation on the structure without the crack to compute the stress fields and (ii) several calculations with different crack surfaces to establish the energy equilibrium. If we are considering small cracks in a large structure, this may be computationally costly as fine meshes are required in the area close the crack location. The MA approach provides an alternative and efficient method to apply the CC, which is described in the sequel. Notice that in a bidimensional problem (2) can be expressed as

$$ \Delta W_{el} + \mathcal{G}_c \ell = \Delta W_{ext} \qquad (3) $$

where $\ell$ is the newly created crack length (a priori unknown). At the initiation imposed loading $\ell = \ell_c$, the initiation crack length. Moreover, it is important to highlight that in problems where there are notches or pre-existing cracks, the crack nucleation is frequently determined by the initiation GSIFs of the singularity, denoted as $K_i$. These parameters depend on the initiation imposed loading and the geometry of the problem. In the problem represented in Figure 1, there is only one leading term, see (1), and therefore only one initiation GSIF.

## 2.2. The matched asymptotic approach

The MA approach of the CC is based on the fact that the crack can be considered as a small perturbation to the elasticity problem where the structure is subjected to a given loading. It consists in successively considering two problems to be solved at two scales. The first problem, solved in the so-called outer domain, is obtained by considering the full structure and neglecting the crack that initiates. In complement, the second problem focuses only on the inner domain around the crack initiation point, independently of the whole structure under investigation. The final solution is then obtained by matching both problems to obtain the stress and energy balance required for the CC application.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Outer domain where the crack is disregarded, the contour Γ can be used to calculate the Generalized Stress Intensity Factor acting at the V-notch for a given force F. (b) Inner domain where the whole structure is disregarded, the normalized crack length is 1 and the arrows represent the imposed asymptotic displacement fields prescribed at a fictitious boundary sufficiently far from the crack.</figcaption></figure>

**Outer domain.** In the outer problem, the perturbation (i.e. the crack) is neglected and a solution of the problem without perturbation is provided in the outer domain (i.e. the structure without crack). This solution is valid everywhere except in a zone near the crack initiation location, for which a correction to this solution must be brought. The outer domain corresponding to the example given in Figure 1, is shown in Figure 2a. In the outer domain, the loading is described in terms of prescribed displacement or force. Then, the GSIF of the singular point (here, the V-notch tip) can be calculated for a given applied force or displacement. Under the assumptions of small deformation and linear elasticity, the GSIF is proportional to the imposed force. For a given imposed force, the GSIF can be computed using a contour integral [2] on a closed path surrounding the singular point (e.g., Γ in Figure 2a). The GSIF calculation based on the contour integral can be implemented in 3D [38] or in 2D for isotropic [21, 39] or anisotropic [40, 41] materials, for multi-material configurations [40, 42–44], or even based on displacement fields measured experimentally by digital image correlation [45]. Other approaches also exist to compute the GSIF, such as the quasidual function method [46, 47], least square fitting [48] or an extraction from the strain energy density [49]. The solution obtained in the outer domain is valid except near the crack initiation location, which requires a correction representative of the initial problem (Figure 1).

**Inner domain.** The correction to the solution obtained in the outer domain without the perturbation is obtained through the second problem which is solved in the inner domain. It consists in focusing only in a zone near the crack initiation location, providing a detailed description of the crack around the singular point, regardless of the entire structure itself. The inner domain thus corresponds to the singular point that would lie in an infinite medium and would be subjected to remote asymptotic displacement or stress fields. The prescribed loading is thus described in terms of GSIF. An example of inner domain corresponding to the problem depicted in Figure 1 is shown in Figure 2b. In the inner domain, the space variables are normalized with respect to the crack length so that the normalized initiation crack length is 1. Since the whole structure geometry and boundary conditions are disregarded in the inner domain, the asymptotic displacement fields are prescribed as boundary conditions in order to obtain the stress and energy balance required to solve the CC. The solution derived in the inner domain is thus accurate in a zone near the crack initiation location.

**Matching inner and outer problem solutions.** Solving the problems in the outer and inner domains yields two solutions (displacement fields) that accurately represent the initial problem
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Examples of configurations: (a) Inclusion, (b) crack ahead of a V-notch or (c) cavity close to a free edge, that can be studied applying the MA approach of the CC. For display purposes in the representation $\ell$ is purposely not small compared to any dimensions of the structure.</figcaption></figure>

of the structure containing a small crack respectively far from and close to the crack initiation location. Matching both solutions also requires a common description of the applied loading. Since it is only described by the GSIF in the inner domain, it justifies the need of calculating the relation between the GSIF and the applied force or displacement in the outer domain. The next step in the MA approach consists in combining both solutions to finally solve the initial problem. This is done by matching both displacement fields in a zone that is (i) sufficiently far from the singular point in the inner domain and (ii) sufficiently close to it in the outer domain. The matching conditions finally enable obtaining the stress and energy balance corresponding to the initial problem under investigation (Figure 1) and further apply the CC for studying crack initiation.

**Solving the CC.** The matching of the inner and outer problem solutions provide a general solution that is accurate over the whole domain under investigation. It yields the displacement fields in the whole structure in presence of a crack. It thus enables calculating the stress fields before crack initiation (Equation (1)) as well as the elastic strain energy variation due to crack initiation (Equation (2)) for a given loading. It finally yields all the ingredients required to solve the CC. The remaining step consists in determining the minimum imposed loading and the corresponding crack length for which both stress and energy conditions are fulfilled.

### 3. Formulation of the approach

The matched asymptotic expansion is used in mechanical engineering to predict the solution, i.e. the displacement field $U_\ell(x_1, x_2)$ (where $(x_1, x_2)$ represents the Cartesian coordinates) in the vicinity of an element that can be an inclusion, a crack or a cavity, see Figure 3. This element is frequently called perturbation, since it is assumed that its size $\ell$ is small compared to any dimensions of the structure.

As an example to illustrate the formulation of the problem, a small cavity located at the tip of a V-notched is considered, see Figure 4, where the notation of the problem that is approximated is represented. The domain $\Omega_\ell$ has an outer contour $\Gamma = \Gamma_V \cup \Gamma_N \cup \Gamma_D \cup \Gamma_\ell$. The contour $\Gamma_D$ is characterized by an imposed displacement $\bar{U}$, whereas the contours $\Gamma_N$, $\Gamma_V$ and $\Gamma_\ell$ have a stress-free boundary conditions. The notation $\Gamma_V$ is referred to the contour of the V-notch and $\Gamma_\ell$ to the one of the small perturbation.

Hence, the set of equations that defines the actual problem is

$$-\nabla_x \cdot \underline{\underline{\sigma}}_\ell = 0 \quad \text{in } \Omega_\ell, \qquad (4)$$

$$\underline{\underline{\sigma}}_\ell = \mathbf{C} : \nabla_x \underline{\underline{U}}_\ell, \qquad (5)$$

$$\underline{\underline{\sigma}}_\ell \cdot \underline{\underline{n}} = \bar{h} \quad \text{on } \Gamma_N, \qquad (6)$$
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Representation of the notations in (a) the inner domain and (b) the outer domain for the example of a cavity ahead of a V-notch.</figcaption></figure>

$$ \underline{\sigma_{\ell} \cdot n} = 0 \quad \text{on } \Gamma_V \cup \Gamma_{\ell}, \qquad (7) $$

$$ \underline{U_{\ell}} = \underline{\bar{U}} \quad \text{on } \Gamma_D, \qquad (8) $$

where $\nabla_x$ is referred to the coordinates system of the actual domain $x_1, x_2$. In the MA approach, a twofold representation of $U_{\ell}(x_1, x_2)$ is proposed in the form of an outer and an inner expansion. Notice that it is assumed that the specimen is in the absence of body forces.

**Outer expansion.** In this approximation, the actual solution is represented as

$$ \underline{U}_{\ell}(x_1, x_2) = \underline{U}_{0}(x_1, x_2) + \cdots \qquad (9) $$

where $U_0(x_1, x_2)$ is the solution of the same elasticity problem considering that the perturbation is not observable in the domain, i.e., solved in an unperturbed domain $\Omega_0$. As an example, Figure 4b represents $\Omega_0$ associated with the actual domain of Figure 4a. The second term in (9) denoted with an ellipsis is a “small correction” that decreases to 0 as $\ell \to 0$. The solution $U_0(x_1, x_2)$ is a good approximation of $U_{\ell}(x_1, x_2)$ far away from the perturbation. For this reason, it is called the outer field. The set of equations that defines $U_0(x_1, x_2)$ is

$$ -\nabla_x \cdot \underline{\sigma_0} = 0 \quad \text{in } \Omega_0, \qquad (10) $$

$$ \underline{\sigma_0} = \underline{C} : \nabla_x \underline{U_0}, \qquad (11) $$

$$ \underline{\sigma_0 \cdot n} = \underline{h} \quad \text{on } \Gamma_N, \qquad (12) $$

$$ \underline{\sigma_0 \cdot n} = 0 \quad \text{on } \Gamma_V, \qquad (13) $$

$$ \underline{U_0} = \underline{\bar{U}} \quad \text{on } \Gamma_D. \qquad (14) $$

Notice that a better approximation of the outer expansion can be achieved by considering higher order terms. Particularly, Leguillon et al. considered the second outer term in [50].

**Inner expansion.** A second expansion can be used to approximate the actual solution by introducing the change of variables $y_i = x_i/\ell$ and $\rho = r/\ell$. In the limit when $\ell \to 0$ we obtain an unbounded domain $\Omega_{in}$ in which the dimensionless characteristic length of the perturbation is now equal to 1, see Figure 5 as an example, where the chosen characteristic length is the diameter of the cavity.

The inner expansion is therefore expressed as

$$ \underline{U_{\ell}}(x_1, x_2) = \underline{U_{\ell}}(\ell y_1, \ell y_2) = F_0(\ell) \underline{V_0}(y_1, y_2) + F_1(\ell) \underline{V_1}(y_1, y_2) + \cdots \qquad (15) $$
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Scheme of the inner problem.</figcaption></figure>

The set of equations related to the two terms $\underline{V_0}(y_1, y_2)$ and $\underline{V_1}(y_1, y_2)$ are

$$ 
\begin{array}{l@{\qquad}l}
-\nabla_y \cdot \underline{\tilde{\sigma}_0} = 0 & \text{in } \Omega_{\text{in}}, & -\nabla_y \cdot \underline{\tilde{\sigma}_1} = 0 \text{ in } \Omega_{\text{in}}, \\
\underline{\tilde{\sigma}_0} = \underline{\mathbf{C}} : \nabla_y \underline{V_0} & & \underline{\tilde{\sigma}_1} = \underline{\mathbf{C}} : \nabla_y \underline{V_1} \\
\underline{\tilde{\sigma}_0} \cdot \underline{n} = 0 & \text{on } \Gamma_V \cup \Gamma_\ell & \underline{\tilde{\sigma}_1} \cdot \underline{n} = 0 \text{ on } \Gamma_V \cup \Gamma_\ell
\end{array}
$$

The problems indicated above are well-posed when the so-called matching conditions are added to these sets of equations. As a results, it is obtained an inner expansion that is a good approximation of the actual solution $\underline{U}_\ell(x_1, x_2)$ in the neighbourhood of the perturbation.

**Matching conditions.** Since the outer expansion is a good approximation of the actual solution far away from the location of the perturbation and the inner expansion is a good approximation in the vicinity of the perturbation, there must exist an intermediate region where both expansions are valid. In that region the matching conditions are defined. The behaviour of the far field near the origin can be described by an expansion in powers of the distance to the singular point $r$, that can be the Taylor's expansion in the case of a smooth stress field or the Williams' expansion in case of a singularity. Assuming the example of the cavity in a V-notch highlighted in Figure 4, the Williams' expansion can be applied, normally expressed in polar coordinates as

$$ \underline{U}(r, \theta) = \underline{U}(0, 0) + K r^{\lambda} \underline{u}(\theta) + o(r^{\lambda}), \quad (16) $$

assuming that the dominant term is real and have multiplicity one. The matching conditions can be expressed as

$$ F_0(\ell) \underline{V}_0(y_1, y_2) \approx \underline{U}(0, 0), \quad \text{when } \rho \to \infty \qquad (17) $$

$$ F_1(\ell) \underline{V}_1(y_1, y_2) \approx K \ell^{\lambda} \rho^{\lambda} \underline{u}(\theta), \quad \text{when } \rho \to \infty \qquad (18) $$

where the term $\approx$ means “behaves like” and $\rho = r/\ell = \sqrt{y_1^2 + y_2^2}$. It can thus be set:

$$ F_0(\ell) = 1 \quad \text{and} \quad \underline{V}_0(y_1, y_2) \approx \underline{U}(0, 0), \quad \text{when } \rho \to \infty \qquad (19) $$

$$ F_1(\ell) = K \ell^{\lambda} \quad \text{and} \quad \underline{V}_1(y_1, y_2) \approx \rho^{\lambda} \underline{u}(\theta), \quad \text{when } \rho \to \infty \qquad (20) $$

However, it can be shown that the matching condition over $\underline{V_1}(y_1, y_2)$ does not fulfill the Lax-Milgram theorem, since it has an infinite energy in the unbounded domain $\Omega_{\text{in}}$, while it should decrease to 0 at infinity to have a finite energy. For this reason, the superposition principle is applied,

$$ \underline{V}_1(y_1, y_2) = \rho^{\lambda} \underline{u}(\theta) + \hat{\underline{V}}_1(y_1, y_2) \qquad (21) $$

where $\hat{\underline{V}}_1(y_1, y_2)$ is the solution to a well-posed problem. The set of equations that defines the new term $\hat{\underline{V}}_1(y_1, y_2)$ is:

$$ 
\begin{array}{l}
-\nabla_y \cdot \underline{\hat{\sigma}_1} = 0 \quad \text{in } \Omega_{\text{in}}, \\
\underline{\hat{\sigma}_1} = \underline{\mathbf{C}} : \nabla_y \underline{\hat{V}_1},
\end{array}
$$
<figure><img src="image_8.png" /><figcaption>Figure 6. FEM model of the outer domain problem.</figcaption></figure>

$$ 
\begin{array}{l}
\hat{\sigma}_1 \cdot \underline{n} = 0 \quad \text{on } \Gamma_\ell, \\
\hat{V}_1 \approx 0 \quad \text{as } \rho \to \infty, \\
\hat{\sigma}_1 \cdot \underline{n} = -\mathbf{C} : \nabla_y (\rho \underline{u}(\theta)) \cdot \underline{n} \quad \text{on } \Gamma_V
\end{array}
$$

Note that the terms $\hat{V}_1(y_1, y_2)$ and $V_1(y_1, y_2)$ are independent of the global geometry and the applied load, that are included in the GSIF $K$. Finally, it yields a new expression for the inner expansion, where $\hat{V}_1(y_1, y_2)$ has finite energy at infinity,

$$ \underline{U}_\ell(x_1, x_2) = \underline{U}_\ell(\ell y_1, \ell y_2) = F_0(\ell) \underline{U}(0, 0) + K \ell^\lambda [\rho^\lambda \underline{u}(\theta) + \hat{V}_1(y_1, y_2)] + o(\ell^\lambda). \quad (22) $$

The MA approach thus enables determining the displacement field accounting for the perturbation by the initiation crack length. It yields all the ingredients required to compute the initiation loading and crack length by further implementation of the CC, i.e. the stress field before crack initiation and the elastic strain energy release due to the crack.

## 4. Numerical implementation of the matched asymptotic approach

The strategy outlined and described in previous sections can be implemented using the typical computational tools employed in solid mechanics, such as the Finite Element Method and the Boundary Element Method. In this work, the applications will focus on the Finite Element Method, but the idea directly applies to other computational methods. The numerical implementation will be applied initially to the problem described in Figure 1 as a simple case and after it will be extended to more complex cases, where some other aspects have to be taken into account.

As described in Section 2.1, the implementation of the coupled criterion requires the evaluation of the stress and energy criteria separately. Thus, the objective of the numerical implementation of the MA will be the evaluation of these two criteria. Following the strategy of the MA approach, two domains, inner and outer, are used for this objective. The numerical implementation will consist on using the Finite Element Methods to compute the necessary elastic solutions in the two domains. Some of them will be used to impose matching between the two domains and others to obtain the stresses or the change in elastic strain energy necessary to evaluate the stress and energy criteria.

In this sense, the steps of a CC analysis assisted by MA and FEM for the problem described in Section 2.1 are the following:

(1) Generation of a FEM model of the outer domain without crack. The mesh should be fine enough around the V-notch in order to approximate accurately the displacement
solution around this point, see Figure 6. The minimum required mesh size can be determined based on a mesh convergence analysis.

(2) Computation of the outer domain model for a certain load. If the problem can be considered linear, i.e. material behavior and boundary conditions are linear and it is possible to assume small deformations, the solution will be linear with the load. In this case, a unit load can be applied and the results obtained can be multiplied by the load to get any result.

(3) The displacements $\underline{u}^{FE}(\rho, \theta)$ and traction vector $\underline{t}^{FE}(\rho, \theta)$ in the vicinity of the V-notch are extracted from the outer-domain model in a closed contour around (noted $\Gamma$), typically a circle. Since this first problem contains a V-notch, these results can be directly used to compute the value of the GSIF by using the property of orthogonality of modes through the next contour integral:

$$K \approx \frac{\int_{\Gamma} \underline{t}^{FE}(\rho, \theta) \cdot \rho^{\lambda} \underline{u}^{-}(\theta) - \underline{t}^{-}(\rho, \theta) \cdot \underline{u}^{FE}(\rho, \theta) d\Gamma}{\int_{\Gamma} \underline{t}(\rho, \theta) \cdot \rho^{\lambda} \underline{u}^{-}(\theta) - \underline{t}^{-}(\rho, \theta) \cdot \rho^{\lambda} \underline{u}(\theta) d\Gamma}, \quad (23)$$

where $\underline{u}(\theta)$ is the displacement function of the singular mode corresponding to the exponent $\lambda$ and the GSIF $K$ and $\underline{u}^{-}(\theta)$ is the corresponding dual one. The traction vectors $\underline{t}$ and $\underline{t}^{-}$ respectively correspond to the displacement functions of the singular and dual modes. Once $K$ is known, the stresses before the crack initiation (necessary for the evaluation of the stress criterion) can be extracted directly from the singular expansion, assuming the crack length at initiation is sufficiently small compared with the size of the region governed by the first term of the William's expansion.

(4) Generation of a FE model of the inner domain with crack, see Figure 7. The external radius should be much larger than the crack length. Typically, the dimensionless crack length is set to 1, and the external radius should be at least 200. The displacements given by the singular mode are prescribed at the external boundaries.

(5) Two versions of this model are computed:

(a) Submodel 0: This submodel corresponds to the state just before the crack initiation. The two crack faces (and then their corresponding nodes) are tied to each other, in order to model the situation without crack.

(b) Submodel 1: This submodel corresponds to the state just after the crack initiation. Crack faces are stress-free.

From these two submodels, the change in potential elastic energy $\Delta W_{el}$ necessary for the evaluation of the energy criterion (Equation (2)). Several strategies can be used:

(a) Crack closure technique: The change in potential elastic energy can be obtained by the work of the virtual forces necessary to close the cracks for problems with displacement control loading. Assuming linearity in the process, the change in potential elastic energy can be computed with the following expression:

$$\Delta W_{\text{ext}} - \Delta W_{\text{el}} = \int_{\Gamma_{\text{crack}}} \frac{1}{2} \underline{t}^{0+} \cdot (\underline{u}^{1+} - \underline{u}^{1-}) d\Gamma_{\text{crack}}. \quad (24)$$

where $\underline{t}^{0+}$ is the traction vector at the submodel 0 at one of the crack faces, named positive. The terms $\underline{t}^{1+}$ and $\underline{t}^{1-}$ refer to the displacements at the positive and the opposite faces respectively in the submodel 1. In case of the existence of force control loads, the below expression should include the work done by the external forces during crack closure. The fact of having to take into account the external forces takes away the advantages of this technique, so the next technique would be recommended in case of force control loading.
(b) Technique based on the change in the work of the external forces applied at the external boundary. The change in potential energy can be computed by:

$$ \Delta W_{\text{ext}} - \Delta W_{\text{el}} = \int_{\Gamma_{\text{ext}}} \frac{1}{2}(\underline{t}^0 \underline{u}^1 - \underline{t}^1 \underline{u}^0)d\Gamma_{\text{ext}}, \quad (25) $$

where $\underline{t}^0$ and $\underline{t}^1$ are the traction vectors at the external boundary in the submodels 0 and 1 respectively. In both strategies, it is more convenient to extract the nodal forces directly from the FEM results and multiply them by the nodal displacement, instead of computing the contour integral.

(c) Technique based on the change in elastic strain energy: This technique is based on computing directly the change in elastic strain energy. This is especially useful when the problem is displacement controlled, otherwise it will be necessary to add the work done by external forces during crack initiation. Thus, the next general expression is:

$$ \Delta W_{\text{ext}} - \Delta W_{\text{el}} = \int_{\Gamma_N} \frac{1}{2} \underline{t}(\underline{u}^1 - \underline{u}^0) d\Gamma + \int_{\Omega} \frac{1}{2} (\underline{\sigma}^0 : \underline{\varepsilon}^0 - \underline{\sigma}^1 : \underline{\varepsilon}^1) d\Omega \quad (26) $$

where $\Gamma_N$ corresponds to the boundary with Neumann boundary condition. The first term would vanish in case of displacement control.

(6) Since in this problem $\ell$ is the only characteristic length and the space variables are made dimensionless with respect to $\ell$ in the inner domain, the dimensionless crack length in the inner domain is 1. Therefore, it is only needed to compute $\Delta W_{el}$ for a dimension crack length of 1. In fact, according to the Dimensional Analysis of the inner domain, the change in potential elastic energy in this problem should follow the next expression:

$$ \Delta W_{\text{ext}} - \Delta W_{\text{el}} = \frac{K^2}{E} \ell^{2\lambda} Ah \quad (27) $$

where $h$ is the thickness, and $A$ is a dimensionless parameter that can be obtained by the steps described before for a unit imposed GSIF $K$. In case other characteristic lengths are present in the inner domain, such as a blunt notch radius $r_b$, or other unknown about the crack initiation geometry, such as deviation angle $\alpha$, the term $A$ becomes a dimensionless function that contains this dependence in terms of dimensionless parameters, such as $r_b/\ell$ or $\alpha$. In next sections more complex problems will be presented in this sense.

(7) Once the value of $K$ is estimated, the coupled criterion of the finite fracture mechanics can be applied in a quite straightforward manner:

* Stress criterion: Combining the first term of the expression in (1) and the condition outlined for the stress criterion, this criterion can be expressed as:

$$ K r^{\lambda-1} s(\theta) \geq \sigma_c \quad \forall r, \quad 0 \leq r \leq \ell \quad (28) $$

Since for singular cases ($\lambda < 1$) this function is decreasing with $r$, this condition can also be also written as:

$$ K \ell^{\lambda-1} s(\theta) \geq \sigma_c. \quad (29) $$

* Energy criterion: Combining the expression in (27) with the condition in (2), the energy criterion can be expressed as:

$$ \frac{K^2}{E} \ell^{2\lambda} Ah \geq \mathcal{G}_c \ell h \quad (30) $$

The two criteria have to be fulfilled simultaneously. In this case, since the left term in Equation (29) is a decreasing function of $\ell$ and the left term in Equation (30) is an increasing function of $\ell$, the minimum value of $K$ for which the two conditions are
<figure><img src="image_9.png" /><figcaption>**Figure 7.** FEM model of the inner domain problem.</figcaption></figure>

fulfilled is given solving the conditions as two equations with two unknowns $\ell$ and $K$. The solution, for a critical length $\ell_i$, is:

$$ \ell_i = \frac{\mathcal{G}_c E s(\theta)^2}{A \sigma_c^2} \qquad (31) $$

and for the critical GSIF $K_i$,

$$ K_i = \left( \frac{\mathcal{G}_c E}{A} \right)^{\lambda-1} \left( \frac{\sigma_c}{s(\theta)} \right)^{2\lambda-1} \qquad (32) $$

Since $K_i$ is directly proportional to the external load, the critical load for crack initiation can be calculated directly from the critical value for $K_i$.

# 5. Implementation in complex cases and examples

Previous sections were dedicated to present the MA approach of the CC in a relatively simple 2D configuration involving a single length parameter (the crack length) under opening mode, such as in [51]. In this section, we provide a detailed overview of more complex configurations and examples to thoroughly demonstrate the full potential of the MA approach.

## 5.1. *Cases with two lengths involved*

The problem under investigation may involve more than one characteristic length parameter. An emblematic example is the case of crack initiation at the tip of a blunted V-notch [39, 52, 53], including as length parameters both the V-notch radius $r_V$ and the crack length $\ell$. In such kind of problems, one must choose which parameter will be used for the inner expansion. The correction of the outer expansion will thus be obtained either in the inner domain where the initiation crack length is 1 and the dimensionless V-notch radius is $r_V/\ell$, or in the inner domain where the V-notch radius is 1 and the dimensionless initiation crack length is $\ell/r_V$. Either approach is strictly equivalent, so that the final solution is independent of the choice of the normalization parameter. The normalization parameter can thus be chosen for practicality reasons, such as the use of a single (expansion with respect to the V-notch radius) instead of several FE meshes and calculations (expansion with respect to the crack length). Other examples of the MA approach with at least two length parameters include crack initiation at a shallow notch [54], a blunted U-notch [55, 56], or a pore crack initiation [57, 58].
## 5.2. Mixed mode loading

Considering mixed mode loading crack initiation means that the outer and inner expansions involve more than one term in the form $K_j r^{\lambda_j} \underline{u}_j(\theta)$, such as for instance:

$$ \underline{U}(r, \theta) = \underline{U}(0, 0) + K_a r^{\lambda_a} \underline{u}_a(\theta) + K_b r^{\lambda_b} \underline{u}_b(\theta) + \dots \quad (33) $$

This is classically encountered when studying mixed opening and shear mode crack initiation, for instance at a V-notch in 2D [51, 53, 59–61] or in 3D [56, 62]. Higher order terms, such as the T-stress, may also be involved for instance when studying 2D crack deviation [63] or 3D crack front segmentation into facets [64]. In such configurations, the stress field is composed of the sum of the terms corresponding to each mode. The elastic strain energy (and thus the Incremental Energy Release Rate) is the sum of terms corresponding to each mode (in the form $\alpha_j K_j r^{2\lambda_j}$) as well as coupling terms (in the form $\alpha_{jk} K_j K_k r^{\lambda_j + \lambda_k}$). These terms can be computed either by directly prescribing mixed-mode loading boundary conditions in the inner domain, or successively prescribing single-mode loading boundary conditions then using the superposition principle.

## 5.3. Interface cracking

The MA of the CC can be applied to study crack initiation at an interface between two materials, possibly in presence of a singular point [59, 65]. The characteristic exponent and corresponding mode depend on the elastic property contrast and singularity geometry [66, 67]. It may even result in complex characteristic exponents and modes [68], which does not prevent the implementation of the MA approach as the displacement field remains real but may requires the manipulation of complex numbers. Applications of the MA of the CC to interface cracking include bimaterial joint failure [69–72] with the consideration of residual stresses [73], crack initiation in microelectronics structures [74, 75], bond failure of a SiC–SiC brazed assembly [76] and crack deflection in layered ceramics [77] or ceramic matrix composites [78].

## 5.4. Comparison with other models

The MA approach of the CC was compared to other fracture models such as Cohesive Zone Model (CZM), Phase Field (PF) approach for fracture, Thick Level Set (TLS) and Strain Energy Density (SED) approach. It was shown that the MA approach of the CC and Dugdale CZM yielded similar initiation GSIF when applied to predict crack initiation at a V-notch [79, 80], which was a particular case of the more general result that the CZM traction-separation profile corresponding to the CC actually depends on the geometry, the type of loading, the cracking mechanism and the adopted stress criterion [23]. The comparison between the TLS and the matched asymptotic approach of the CC revealed similar apparent strengths for all cases provided the assumptions of the MA approach are satisfied [81], as well as a dependence of the TLS results to the choice of the stress decrease as a function of the crack opening. This dependence can be put in parallel to the traction-separation profile of CZM so that one stress decrease function could be identified in the TLS to retrieve the CC results. Abaza et al. [82] showed that the PF regularization length could be calibrated so that the apparent SIF at crack nucleation were similar to those obtained with the CC. Similar variations of the apparent SIF at crack nucleation in notched ceramic specimens were then obtained for different notch geometries using CC and PF. The SED and CC comparison in the case of a V-notch loaded under in-plane shear revealed that the analytically computed initiation GSIF were proportional to powers of $K_{Ic}$ and $\sigma_c$ for both methods. The proportionality factor was a function of the notch angle for the CC whereas it was a function of the Poisson’s ratio for the
SED approach. Basing the control volume over which the strain energy density was averaged, both approaches predicted similar apparent SIF at crack initiation [83].

### 5.5. 3D

There is no restriction regarding the 3D implementation of the MA approach [84]. The displacement field is based on 3D William's expansions [85, 86] which takes the same general form as the 2D one. The 3D singular exponent and associated eigenmode can be computed in a similar way to the 2D one [2, 87]. The main difference compared to the 2D case is related to the dual mode. If the primal mode is $r^{\lambda}u(\theta, \phi)$, then the dual mode is $r^{-\lambda-1}u_{3D}^{-}(\theta, \phi)$ in 3D instead of $r^{-\lambda}u_{2D}^{-}(\theta)$ in 2D. Otherwise, the space coordinates expansion is also performed with respect to a characteristic length parameter, which can be chosen as the crack extension in a given direction [56, 64, 84] or with respect to the square root of the crack surface [62]. Overall, the CC has takes the same form as in 2D and the same contour integral can be used to compute the GSIF [38], except that the integration domain is a 2D surface in the volume encompassing the singular point instead of a 1D curve in a surface. The main difficulty for the 3D CC application is the crack shape determination. A 2D crack can generally be described by its length and angle, but there is an infinite number of possible 3D crack surface. Even assuming a planar crack path, this plane must be determined, which can be done for instance by maximizing the stress criterion [56] or based on both stress and energy requirements by minimizing the initiation loading [64]. Then, in the crack plane there is still an infinite number of 2D curves to describe the crack front. An option to overcome this difficulty is to adopt a parameterized description of the crack shape [62, 64] or derive the crack shapes from the stress isocontours [56].

### 5.6. Other applications

Some works including the MA approach of the CC also cover studying the presence of a process zone or damage zone ahead of a V-notch or a crack [88–90]. It also provided a robust method for small crack detection based on displacement fields measured by digital image correlation [50, 91]. Several works about ceramic failure were also proposed, such as cracking in layered ceramics [92], platelet-based ceramics [41, 93, 94] or surface defects in polycrystalline ceramics [95]. Other works concerned brittle fracture size effect [96], strength anisotropy of 3D printed materials [97], multicracking of a stiff inclusion in a soft matrix [98] or elliptical hole-induced crack initiation [99].

## 6. Conclusions and further developments

Crack initiation is an unsolved problem in Linear Elastic Fracture Mechanics, that does manage correctly to predict crack propagation in diverse material systems. However, crack initiation requires the employment of more complex, computational costly and controversial tools. This review shows how Matched Asymptotics and the Coupled Criterion of the Finite Fracture Mechanics have been combined in the literature to predict crack initiation near stress singularities, such as V-notches or multimaterial corners, or related to, such as U-notches.

The main advantages of this combination are: (i) it allows to obtain predictions with a very low computational cost, where typically only linear models are involved and quasianalytical expressions can be obtained for the load prediction, (ii) The results are easily generalized for variations in material properties and even geometry, (iii) it is physically based, thus results can be interpreted, explained and tailored following a physical reasoning and finally (iv) according to the literature review it presents a good agreement with experiments.
The main disadvantages are (i) the solutions are still limited to crack length at onset that has to be very small compared with the characteristic lengths of the problem and (ii) it is necessary to introduce special assumptions when nonlinearities are involved.

Other methods have been used to predict crack initiation in the stress singularities or in the problems target of the combination of matched asymptotics and finite fracture mechanics: Cohesive Zone Models (CZM) prescribe a cohesive law between a pair of surfaces, relating force and separation, see e.g. [80]. This method is very versatile and presents good agreement with experiments, but requires setting crack geometry before initiation and typically involve nonlinear computational models. In the last decades gradient-damage-based models such as that named Phase Field have been extensively developed. These models are based on the regularization of the crack through a regularization length. It has been proven that when this length vanishes, the result of LEFM is recovered. Then, for V-notches and related problems, these models present the same problems as LEFM. To overcome this issue, several strategies have been proposed, such as assuming that the regularization length is a material parameter [100], defining a CZM [101], or understanding this regularization length in the context of the CC [90].

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

## Dedication

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

## Acknowledgments

This paper is entirely dedicated to our friend and colleague, Pr. Dominique Leguillon.

SJ-A and IGG acknowledge the funding received from the European Union's Horizon 2020 research and innovation programme under Marie Sklodowska-Curie grant agreement No. 861061-NEWFRAC. SJ-A also acknowledges the Iberdrola Foundation under the Marie Sklodowska-Curie Grant Agreement No 101034297. IGG also acwnoledges the support of Ministerio de Ciencia e Innovación de España (Project PID2020-117001GB-I00/AEI/10.13039/501100011033).

## References

[1] M. Van Dyke, *Perturbation Methods in Fluid Mechanics*, Academic Press: New York, 1964.

[2] D. Leguillon and E. Sanchez-Palencia, *Computation of Singular Solutions in Elliptic Problems and Elasticity*, Wiley: USA, 1987.

[3] P. Lagerstrom, *Matched Asymptotic Expansions: Ideas and Techniques*, Applied Mathematical Sciences, Springer: Berlin, 1988.

[4] J. Kevorkian and J. Cole, *Multiple scale and singular perturbation methods*, Springer Science & Business Media: New York, 2012.

[5] R. O'Malley, *The Method of Matched Asymptotic Expansions and Its Generalizations*, Historical Developments in Singular Perturbations, Springer: Cham, 2014.

[6] F. Lene and D. Leguillon, "Homogenized constitutive law for a partially cohesive composite material", *Int. J. Solids Struct.* **18** (1982), no. 5, pp. 443-458.
[7] A. Bensoussan, L. Lions and G. Papanicolaou, *Asymptotic Analysis for Periodic Structures. Studies in Mathematics and its Application*, North-Holland: Amsterdam, 1978.

[8] D. Leguillon, "Comparison of matched asymptotics, multiple scalings and averages in homogenization of periodic structures", *Math. Models Methods Appl. Sci.* **07** (1997), no. 05, pp. 663–680.

[9] D. Leguillon and E. Sanchez-Palencia, "Approximation of a two dimensional problem of junctions", *Comput. Mech.* **6** (1990), pp. 435–455.

[10] D. Leguillon, "Asymptotic and numerical analysis of a crack branching in non-isotropic materials", *Eur. J. Mech. Solids* **12** (1993), pp. 33–51.

[11] P. Sicsic and J. Marigo, "From gradient damage laws to Griffith's theory of crack propagation", *J. Elast.* **113** (2013), no. 1, pp. 55–74.

[12] J. Cook and J. Gordon, "A mechanism for the control of crack propagation in all-brittle systems", *Proc. R. Soc. A* **282** (1964), pp. 508–520.

[13] D. Leguillon, C. Lacroix and E. Martin, "Interface debonding ahead of a primary crack", *J. Mech. Phys. Solids* **48** (2000), no. 10, pp. 2137–2161.

[14] D. Leguillon, "A method based on singularity theory to predict edge delamination of laminates", *Int. J. Fract.* **100** (1999), no. 1, pp. 105–120.

[15] D. Leguillon and E. Sanchez-Palencia, "On 3D cracks intersecting a free surface in laminated composites", *Int. J. Fract.* **99** (1999), no. 1, pp. 25–40.

[16] D. Leguillon, "Interface crack tip singularity with contact and friction", *C. R. Acad. Sci. IIB* **327** (1999), no. 5, pp. 437–442.

[17] D. Leguillon and R. Abdelmoula, "Mode III near and far fields for a crack lying in or along a joint", *Intl J. Solids Struct.* **37** (2000), no. 19, pp. 2651–2672.

[18] J. Leblond and D. Leguillon, "The stress intensity factors near an angular point on the front of an interface crack", *Eur. J. Mech. A Solids* **18** (1999), no. 5, pp. 837–857.

[19] M. Williams, "Stress singularities resulting from various boundary condition sinangular corners of plates inex- tension", *J. Appl. Mech.* **19** (1952), pp. 526–528.

[20] D. Leguillon, "A criterion for crack nucleation at a notch in homogeneous materials", *C. R. Acad. Sci. IIB* **329** (2001), no. 2, pp. 97–102.

[21] D. Leguillon, "Strength or toughness? A criterion for crack onset at a notch", *Eur. J. Mech. A Solids* **21** (2002), no. 1, pp. 61–72.

[22] P. Weißgräber, D. Leguillon and W. Becker, "A review of finite fracture mechanics: crack initiation at singular and non-singular stress raisers", *Arch. Appl. Mech.* **86** (2016), no. 1–2, pp. 375–401.

[23] A. Doitrand, T. Duminy, H. Girard and X. Chen, *A review of the coupled criterion*, preprint, 2024. Online at https://hal.science/hal-04023438.

[24] A. Sapora, A. Torabi, S. Etesam and P. Cornetti, "Finite Fracture Mechanics crack initiation from a circular hole", *Fatigue Fract. Eng. Mater. Struct.* **41** (2018), no. 7, pp. 1627–1636.

[25] P. Cornetti and A. Sapora, "Penny-shaped cracks by Finite Fracture Mechanics", *Int. J. Fract.* **219** (2019), pp. 153–159.

[26] T. Methfessel, C. El Yaakoubi-Mesbah and W. Becker, "Failure analysis of crack-prone joints with Finite Fracture Mechanics using an advanced modeling approach for the adhesive layer", *Int. J. Adhesion Adhesives* **130** (2024), article no. 103608.

[27] M. Kashtalyan, I. García and V. Mantič, "Evolution of crack density in cross-ply laminates—application of a coupled stress and energy criterion", *Key Eng. Mater.* **713** (2016), pp. 262–265.

[28] A. Doitrand, C. Fagiano, N. Carrére, V. Chiaruttini and M. Hirsekorn, "Damage onset modeling in woven composites based on a coupled stress and energy criterion", *Eng. Fract. Mech.* **169** (2017), pp. 189–200.

[29] J. Li, E. Martin, D. Leguillon and C. Dupin, "A finite fracture model for the analysis of multi-cracking in woven ceramic matrix composites", *Composites B* **139** (2018), pp. 75–83.

[30] J. Vereecke, C. Bois, J. Wahl, T. Briand, L. Ballère and F. Lavelle, "Explicit modelling of meso-scale damage in laminated composites—comparison between finite fracture mechanics and cohesive zone model", *Compos. Sci. Technol.* **253** (2024), article no. 110640.

[31] A. Doitrand, G. Molnár, D. Leguillon, E. Martin and N. Carrère, "Dynamic crack initiation assessment with the coupled criterion", *Eur. J. Mech. A Solids* **93** (2022), article no. 104483.

[32] X. Chen, A. Doitrand, N. Godin and C. Fusco, "Dynamic crack initiation in circular hole PMMA plates considering nonlinear elastic material behavior", *Theor. Appl. Fract. Mech.* **124** (2023), article no. 103783.

[33] A. Doitrand, S. Rubeck, S. Meille, J. Chevalier, P. Steyer and S. Gallois-Garreignot, "Influence of debonding and substrate plasticity on thin film multicracking", *Theor. Appl. Fract. Mech.* **131** (2024), article no. 104375.

[34] J. Li, D. Leguillon, E. Martin and X. Zhang, "Numerical implementation of the coupled criterion for damaged materials", *Int. J. Solids Struct.* **165** (2019), pp. 93–103.
[35] A. Griffith, "The theory of rupture", in *First International Congress on Applied Mechanics, Farnborough*, 1924, pp. 55-63.

[36] G. Irwin, "Analysis of stresses and strains near the end of a crack traversing a plate", *J. Appl. Mech.* **24** (1957), no. 3, pp. 361-364.

[37] G. Irwin, *Fracture*, Springer: Berlin, Heidelberg, 1958, pp. 551-590.

[38] A. Doitrand, D. Leguillon and E. Martin, "Computation of generalized stress intensity factors of 3D singularities", *Int. J. Solids Struct.* **190** (2020), pp. 271-280.

[39] D. Leguillon, D. Quesada, C. Putot and E. Martin, "Prediction of crack initiation at blunt notches and cavities—size effects", *Eng. Fract. Mech.* **74** (2007), no. 15, pp. 2420-2436.

[40] P. E. Labossiere and M. L. Dunn, "Stress intensities at interface corners in anisotropic bimaterials", *Eng. Fract. Mech.* **62** (1999), no. 6, pp. 555-576.

[41] T. Duminy, V. Ayyalasomayajula, A. Doitrand and S. Meille, "Influence of elastic and toughness anisotropy on crack initiation", *Int. J. Solids Struct.* **302** (2024), article no. 112950.

[42] M. Herrera-Garrido, V. Mantič and A. Barroso, "Computational semi-analytic code for stress singularity analysis", *Procedia Struct. Integr.* **42** (2022), pp. 958-966. 23 European Conference on Fracture.

[43] M. Herrera-Garrido, V. Mantič and A. Barroso, "A powerful matrix formalism for stress singularities in anisotropic multi-material corners. Homogeneous (orthogonal) boundary and interface conditions", *Theor. Appl. Fract. Mech.* **119** (2022), article no. 103271.

[44] M. Herrera-Garrido, V. Mantič and A. Barroso, "A semi-analytical matrix formalism for stress singularities in anisotropic multi-material corners with frictional boundary and interface conditions", *Theor. Appl. Fract. Mech.* **129** (2024), article no. 104160.

[45] A. Doitrand, D. Leguillon and R. Estevez, "Experimental determination of generalized stress intensity factors from full-field measurements", *Eng. Fract. Mech.* **230** (2020), article no. 106980.

[46] Z. Yosibash and Y. Schapira, "Edge stress intensity functions in polyhedral domains and their extraction by a quasidual function method", *Int. J. Fract.* **136** (2005), pp. 37-73.

[47] Z. Yosibash and Y. Schapira, "Edge stress intensity functions along elliptic and part-elliptic 3D cracks", *Eng. Fract. Mech.* **245** (2021), article no. 107477.

[48] M. Dunn, W. Suwito and S. Cunningham, "Fracture initiation at sharp notches: correlation using critical stress intensities", *Int. J. Solids Struct.* **34** (1997), no. 29, pp. 3873-3883.

[49] P. Lazzarin, F. Berto and M. Zappalorto, "Rapid calculations of notch stress intensity factors based on averaged strain energy density from coarse meshes: theoretical bases and applications", *Int. J. Fatigue* **32** (2010), no. 10, pp. 1559-1567.

[50] D. Leguillon, "Determination of the length of a short crack at a v-notch from a full field measurement", *Int. J. Solids Struct.* **48** (2011), no. 6, pp. 884-892.

[51] Z. Yosibash, E. Priel and D. Leguillon, "A failure criterion for brittle elastic materials under mixed-mode loading", *Int. J. Fract.* **141** (2006), pp. 291-312.

[52] D. Leguillon and Z. Yosibash, "Crack onset at a v-notch. Influence of the notch tip radius", *Int. J. Fract.* **122** (2003), pp. 1-21.

[53] E. Priel, Z. Yosibash and D. Leguillon, "Failure initiation at a blunt V-notch tip under mixed mode loading", *Int. J. Fract.* **149** (2008), pp. 143-173.

[54] P. Cornetti, D. Taylor and A. Carpinteri, "An asymptotic matching approach to shallow-notched structural elements", *Eng. Fract. Mech.* **77** (2010), no. 2, pp. 348-358.

[55] D. Picard, D. Leguillon and C. Putot, "A method to estimate the influence of the notch-root radius on the fracture toughness measurement of ceramics", *J. Eur. Ceram. Soc.* **26** (2006), no. 8, pp. 1421-1427.

[56] A. Doitrand and D. Leguillon, "Numerical modeling of the nucleation of facets ahead of a primary crack under mode I+III", *Int. J. Fract.* **123** (2018), no. 1, pp. 37-50.

[57] D. Leguillon and R. Piat, "Fracture of porous materials—influence of the pore size", *Eng. Fract. Mech.* **75** (2008), no. 7, pp. 1840-1853. Critical Distance Theories of Fracture.

[58] A. Doitrand and D. Leguillon, "Asymptotic analysis of pore crack initiation near a free edge", *Theor. Appl. Fract. Mech.* **116** (2021), article no. 103125.

[59] I. García and D. Leguillon, "Mixed-mode crack initiation at a v-notch in presence of an adhesive joint", *Int. J. Solids Struct.* **49** (2012), no. 15, pp. 2138-2149.

[60] A. Doitrand, E. Martin and D. Leguillon, "Numerical implementation of the coupled criterion: matched asymptotic and full finite element approaches", *Finite Elem. Anal. Des.* **168** (2020), article no. 103344.

[61] P. Cornetti, A. Sapora and A. Carpinteri, "Mode mixity and size effect in V-notched structures", *Int. J. Solids Struct.* **50** (2013), no. 10, pp. 1562-1582.

[62] Z. Yosibash and B. Mittelman, "A 3-D failure initiation criterion from a sharp V-notch edge in elastic brittle structures", *Eur. J. Mech. A Solids* **60** (2016), pp. 70-94.

[63] D. Leguillon and S. Murer, "Crack deflection in a biaxial stress state", *Int. J. Fract.* **150** (2008), pp. 75-90.
[64] A. Doitrand, D. Leguillon, G. Molnár and V. Lazarus, "Revisiting facet nucleation under mixed mode I+III loading with T-stress and mode-dependent fracture properties", *Int. J. Fract.* **242** (2023), pp. 85–106.

[65] V. Tran, D. Leguillon and A. Krishnan, "Interface crack initiation at V-notches along adhesive bonding in weakly bonded polymers subjected to mixed-mode loading", *Int. J. Fract.* **176** (2012), pp. 65–79.

[66] D. Leguillon and E. Martin, "The strengthening effect caused by an elastic contrast—part I: the bimaterial case", *Int. J. Fract.* **179** (2013), pp. 157–167.

[67] D. Leguillon and E. Martin, "The strengthening effect caused by an elastic contrast—part II: stratification by a thin stiff layer", *Int. J. Fract.* **179**(1–2) (2013), pp. 169–178.

[68] D. Leguillon and S. Murer, "A criterion for crack kinking out of an interface", *Key Eng. Mater.* **385-387** (2008), pp. 9–12.

[69] D. Leguillon, J. Laurencin and M. Dupeux, "Failure initiation in an epoxy joint between two steel plates", *Eur. J. Mech. A Solids* **22** (2003), no. 4, pp. 509–524.

[70] A. Moradi, D. Leguillon and N. Carrére, "Influence of the adhesive thickness on a debonding—an asymptotic model", *Eng. Fract. Mech.* **114** (2013), pp. 55–68.

[71] J. Felger, P. Rosendahl, D. Leguillon and W. Becker, "Predicting crack patterns at bi-material junctions: a coupled stress and energy approach", *Int. J. Solids Struct.* **164** (2019), pp. 191–201.

[72] J. Felger, N. Stein, C. Frey and W. Becker, "Scaling laws for the adhesive composite butt joint strength derived by finite fracture mechanics", *Compos. Struct.* **208** (2019), pp. 546–556.

[73] C. Henninger and D. Leguillon, "Adhesive fracture of an epoxy joint under thermal and mechanical loadings", *J. Therm. Stresses* **31** (2007), no. 1, pp. 59–76.

[74] H. Brillet-Rouxel, E. Arfan, D. Leguillon, M. Dupeux, M. Braccini and S. Orain, "Crack initiation in Cu-interconnect structures, Microelectronic Engineering", *Phys. Scr.* **83** (2007), pp. 2297–2302.

[75] J. Schlosser, E. Martin, C. Henninger, et al., "CFC/Cu bond damage in actively cooled plasma facing components", *Phys. Scr.* (2007), pp. 204–209.

[76] L. Nguyen, D. Leguillon, O. Gillia and E. Riviere, "Bond failure of a SiC/SiC brazed assembly", *Mech. Mater.* **50** (2012), pp. 1–8.

[77] D. Leguillon, S. Tariolle, E. Martin, T. Chartier and J. Besson, "Prediction of crack deflection in porous/dense ceramic laminates", *J. Eur. Ceram. Soc.* **26** (2006), no. 3, pp. 343–349.

[78] C. Lacroix, D. Leguillon and E. Martin, "The influence of an interphase on the deflection of a matrix crack in a ceramic-matrix composite", *Compos. Sci. Technol.* **62** (2002), no. 4, pp. 519–523.

[79] C. Henninger, D. Leguillon and E. Martin, "Crack initiation at a V-notch—comparison between a brittle fracture criterion and the Dugdale cohesive model", *C. R. Méc.* **335** (2007), no. 7, pp. 388–393.

[80] S. Murer and D. Leguillon, "Static and fatigue failure of quasi-brittle materials at a V-notch using a Dugdale model", *Eur. J. Mech. A Solids* **29** (2010), no. 2, pp. 109–118.

[81] J. Zghal, K. Moreau, N. Moës, D. Leguillon and C. Stolz, "Analysis of the failure at notches and cavities in quasi-brittle media using the Thick Level Set damage model and comparison with the coupled criterion", *Int. J. Fract.* **211** (2018), pp. 253–280.

[82] A. Abaza, J. Laurencin, A. Nakajo, S. Meille, J. Debayle and D. Leguillon, "Prediction of crack nucleation and propagation in porous ceramics using the phase-field approach", *Theor. Appl. Fract. Mech.* **119** (2022), article no. 103349.

[83] A. Campagnolo, F. Berto and D. Leguillon, "Fracture assessment of sharp V-notched components under Mode II loading: a comparison among some recent criteria", *Theor. Appl. Fract. Mech.* **85** (2016), pp. 217–226.

[84] D. Leguillon, "An attempt to extend the 2D coupled criterion for crack nucleation in brittle materials to the 3D case", *Theor. Appl. Fract. Mech.* **74** (2014), pp. 7–17.

[85] M. Williams, "On the stress distribution at the base of a stationary crack", *ASME J. App. Mech.* **24** (1957), pp. 109–114.

[86] T. Apel, D. Leguillon, C. Pester and Z. Yosibash, "Edge singularities and structure of the 3-D Williams expansion", *C. R. Méc.* **336** (2008), no. 8, pp. 629–635.

[87] D. Leguillon, "Computation of 3d-singularities in elasticity", in *Boundary Value Problems and Integral Equations in Nonsmooth Domains*, Lecture Notes in Pure and Applied Mathematics, Marcel Dekker: Paris, 1995, pp. 161–170.

[88] D. Leguillon, "A damage model based on singular elastic fields", *C. R. Méc.* **336** (2008), no. 3, pp. 283–288.

[89] D. Leguillon and Z. Yosibash, "Failure initiation at V-notch tips in quasi-brittle materials", *Int. J. Solids Struct.* **122-123** (2017), pp. 1–13.

[90] A. Doitrand and G. Molnár, *Understanding regularized crack initiation through the lens of Finite Fracture Mechanics*, preprint, 2024. Online at https://doi.org/10.21203/rs.3.rs-4583166/v1.

[91] R. Romani, M. Bornert, D. Leguillon, R. Le Roy and K. Sab, "Detection of crack onset in double cleavage drilled specimens of plaster under compression by digital image correlation—theoretical predictions based on a coupled criterion", *Eur. J. Mech. A Solids* **51** (2015), pp. 172–182.
[92] D. Leguillon, E. Martin, O. Ševeček and R. Bermejo, "Application of the coupled stress-energy criterion to predict the fracture behaviour of layered ceramics designed with internal compressive stresses", *Eur. J. Mech. A Solids* **54** (2015), pp. 94–104.

[93] S. Jiménez-Alfaro and D. Leguillon, "Modelling of glass matrix composites by the Coupled Criterion and the Matched Asymptotics approach. The role of a single platelet", *Theor. Appl. Fract. Mech.* **122** (2022), article no. 103650.

[94] S. Jiménez-Alfaro and D. Leguillon, "Modelling of glass matrix composites by the Coupled Criterion and the Matched Asymptotic Approach. The effect of residual stresses and volume fraction", *Theor. Appl. Fract. Mech.* **128** (2023), article no. 104112.

[95] D. Leguillon and E. Martin, "Prediction of multi-cracking in sub-micron films using the coupled criterion", *Int. J. Fract.* **209** (2018), pp. 187–202.

[96] S. Jiménez-Alfaro and D. Leguillon, "Finite fracture Mechanics at the micro-scale. Application to bending tests of micro cantilever beams", *Eng. Fract. Mech.* **258** (2021), article no. 108012.

[97] L. Xu and D. Leguillon, "Dual-notch void model to explain the anisotropic strengths of 3D printed polymers", *J. Eng. Mater. Technol.* (2019), article no. 014501.

[98] D. Quesada, D. Leguillon and C. Putot, "Multiple failures in or around a stiff inclusion embedded in a soft matrix under a compressive loading", *Eur. J. Mech. A Solids* **28** (2009), no. 4, pp. 668–679.

[99] J. Felger, N. Stein and W. Becker, "Asymptotic finite fracture mechanics solution for crack onset at elliptical holes in composite plates of finite-width", *Eng. Fract. Mech.* **182** (2017), pp. 621–634.

[100] E. Tanné, T. Li, B. Bourdin, J.-J. Marigo and C. Maurini, "Crack nucleation in variational phase-field models of brittle fracture", *J. Mech. Phys. Solids* **110** (2018), pp. 80–99.

[101] J.-Y. Wu and V. P. Nguyen, "A length scale insensitive phase-field damage model for brittle fracture", *J. Mech. Phys. Solids* **119** (2018), pp. 20–42.