\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRMECA,Unicode,biblatex,published]{cedram}
\addbibresource{crmeca20240881.bib}
\newenvironment{noXML}{}{}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\sfrac#1#2{{#1}/{#2}}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\usepackage{etoolbox}
\usepackage{upgreek}
\usepackage[scr=boondoxo, scrscaled=1]{mathalfa}
\csdef{Seqnsplit}{\\}
\newcounter{runlevel}
\let\MakeYrStrItalic\relax
\def\jobid{crmeca20240881}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}} 
\def\xsection#1{}
\def\0{\phantom{0}}
\def\mn{\phantom{$-$}}
\def\refinput#1{}
\def\back#1{}
\def\botline{\\\hline}
\def\tabnote#1{\vskip3pt\parbox{33pc}{#1}}
\newenvironment{inftab}{}{}
\usepackage{multirow} 
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\def\xmorerows#1#2{#2}
\DOI{10.5802/crmeca.285}
\datereceived{2024-10-29}
\daterevised{2025-01-14}
\dateaccepted{2025-01-14}
\ItHasTeXPublished
\def\bPhi{\boldsymbol{\Phi}}
\allowdisplaybreaks
}
{\documentclass[crmeca]{article}
\usepackage{upgreek}
\usepackage{stmaryrd}
\newtheorem{rema}{Remark}
\def\CDRdoi{10.5802/crmeca.285}
\def\newline{\unskip\break}
\def\xmorerows#1#2{\morerows{#1}{#2}}
\def\hyperlink#1#2{#2}
\def\hypertarget#1#2{#2}
\def\sfrac#1#2{{#1}/{#2}}
}
\makeatother

\dateposted{2025-02-06}
\begin{document}

\newcommand{\G}{\mathcal{G}}
\newcommand{\Ginc}{\mathcal{G}_{\mathrm{inc}}}
\newcommand{\Gc}{\mathcal{G}_{\mathrm{c}}}
\newcommand{\KIc}{K_{\mathrm{Ic}}}
\newcommand{\lc}{\ell_{\mathrm{c}}}
\newcommand{\lmat}{\ell_{\mathrm{mat}}}
\newcommand{\sigmac}{\sigma_{\mathrm{c}}}
\newcommand{\rV}{r_{\mathrm{V}}}
\def\uuline#1{\underline{\underline{#1}}}

\begin{noXML}

\CDRsetmeta{articletype}{review}

\title{Review of the matched asymptotic approach of the coupled
criterion}

\alttitle{Revue de l'approche asymptotique du crit\`ere coupl\'e}

\author{\firstname{Sara} \lastname{Jim\'enez-Alfaro}\CDRorcid{0000-0001-7657-2388}\IsCorresp}
\address{Department of Civil and Environmental Engineering, Imperial College
London, Exhibition Road, London, SW7 2AZ, UK}
\address{Department of Engineering Science, University of Oxford, Parks Road,
Oxford, OX1 3PJ, UK}
\email{sara.jimenezalfaro@eng.ox.ac.uk}
\email[S. Jim\'enez-Alfaro]{s.jimenez-alfaro@imperial.ac.uk}

\author{\firstname{Israel} \middlename{Garc\'ia} \lastname{Garc\'ia}\CDRorcid{0000-0002-3894-5304}}
\address{Departamento de Mec\'anica de Medios Continuos y Teor\'ia de
Estructuras, Escuela T\'ecnica Superior de Ingenier\'ia, Escuela
Polit\'ecnica Superior, Universidad de Sevilla, Camino de los
Descubrimienos s/n, 41092 Sevilla, Spain}
\email[I. G. Garc\'ia]{israelgarcia@us.es}

\author{\firstname{Aur\'elien} \lastname{Doitrand}\CDRorcid{0000-0003-1000-3507}}
\address{Universit\'e Lyon, INSA-Lyon, UCBL, CNRS, MATEIS, UMR5510,
F-69621 Villeurbanne, France}
\email[A. Doitrand]{aurelien.doitrand@insa-lyon.fr}

\begin{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. 
\end{abstract}

\begin{altabstract}
Les d\'eveloppements asymptotiques raccord\'es constituent une technique
math\'ematique puissante, largement applicable dans divers domaines de
l'ing\'enierie. L'une de leurs principales utilisations se situe en
m\'ecanique de la rupture, o\`u ils permettent d'obtenir des
approximations pr\'ecises \`a proximit\'e de la pointe des fissures tout
en maintenant une faible complexit\'e de calcul. Cette m\'ethode peut
\^etre int\'egr\'ee de mani\`ere fluide au crit\`ere coupl\'e (CC), qui
permet de pr\'edire l'amor\c{c}age et la propagation des fissures dans les
mat\'eriaux fragiles. Cet article explique comment la technique des
d\'eveloppements asymptotiques raccord\'es peut \^etre utilis\'ee
conjointement avec le crit\`ere coupl\'e dans le cadre de la m\'ecanique
de la rupture, tout en offrant une revue d\'etaill\'ee de la litt\'erature
sur les avanc\'ees r\'ealis\'ees au cours de la derni\`ere d\'ecennie.
\end{altabstract}

\shortrunauthors

\keywords{\kwd{Matched asymptotic expansion}
          \kwd{Coupled criterion}
          \kwd{Fracture mechanics}}

\altkeywords{\kwd{D\'eveloppement asymptotique raccord\'e}
          \kwd{Crit\`ere coupl\'e}
          \kwd{M\'ecanique de la rupture}}

\thanks{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\'on de Espa\~na (Project
PID2020-117001GB-I00/AEI/10.13039/501100011033)}

\maketitle

\end{noXML}

\section{Introduction}\label{sec1}
The matched asymptotic (MA) expansion method is an approach that
enables solving an equation or a system of
equations~\cite{VanDyke1964,Leguillon1987,Lagerstrom1988,Kevorkian2012,OMalley2014}. 
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~\cite{Lene1982}. 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~\cite{Bensoussan1978,Leguillon1997Homogenization}.  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~\cite{Leguillon1990}.  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~\cite{Leguillon1993}. Sicsic and Marigo
studied the propagation of a crack band and derived the conditions for
which it behaves like a Griffith's crack~\cite{Sicsic2013}.  MA
expansions were also used to study the behavior of interface cracks,
for instance to further analyse the ``Cook and Gordon''~\cite{Cook1964}
interface debonding effect ahead of a primary
crack~\cite{Leguillon2000}, edge debonding in
laminates~\cite{Leguillon1999} 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~\cite{Leguillon1999_3D_edgecrack}, also considering
contact and friction between two anisotropic
materials~\cite{Leguillon1999contact}. 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
~\cite{Leguillon2000modeIII}. It was also used to derive the stress
intensity factors near an angular point on the front of an interface
crack~\cite{Leblond1999}. 

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~\cite{Williams1952}, 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~\cite{Leguillon2001CRAS,Leguillon2002} 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~\cite{Weissgraeber2016review,Doitrand2024JTCAM}.
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~\cite{Weissgraeber2016review,Sapora2018,Cornetti2019Penny,Methfessel2024}. 
A second way to implement the CC is through finite element (FE)
simulations of the full structure under investigation including the
crack that
initiates~\cite{Kashtalyan2016,Doitrand2017EFM,Li2018,Vereecke2024}. 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~\ref{sec:Idea_MA}) and the formulation
(Section~\ref{sec:Formulation}) of the MA approach. Then, we describe
its numerical implementation (Section~\ref{sec:Implementation}) and
provide some application examples (Section~\ref{sec:Application}). 

\section{The idea behind matched asymptotics}
\label{sec:Idea_MA}
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~\cite{Doitrand2022EJMAS,Chen2023} or
dissipation mechanisms other than cracking that may occur during
initiation, such as, e.g., plasticity~\cite{Doitrand2024TAFMEC},
diffuse damage~\cite{Li2019} 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 $\ell$ that initiates at the tip of a V-notch (angle $\beta$) in
a specimen loaded under three-point bending
(Figure~\ref{fig:general_problem}). 
Notice that this technique is only valid provided $\ell$ is smaller
than a characteristic dimensions of specimen ($\ell\ll e$ in
Figure~\ref{fig:general_problem}), an initial assumption that should be
checked after the implementation, once the actual value of $\ell$ is
obtained using the coupled criterion.

\begin{figure}
\includegraphics{fig01}
\caption{V-notched three-point bending specimen subjected to a force
$F$ with a crack (length $\ell$) initiating at the V-notch tip. Notice
that $e$ represents the notch depth.\label{fig:general_problem}}
\end{figure}

\subsection{The coupled criterion}\label{ssec:CC}

The main idea behind the CC arises from the following observations: 
\begin{itemize}
\item Considering an energy criterion only, it enables assessing the
propagation of a crack based on the material critical energy release
rate $\Gc$~\cite{Griffith1924,Irwin1957,Irwin1958} but generally fails
to study its initiation.
\item Considering a stress criterion only, it enables assessing
crack initiation based on the material tensile strength $\sigmac$
except in the presence of a singular point.
\end{itemize}
Stress and energy criteria thus appear as complementary and their
combination enables assessing crack initiation in many configurations.
The stress criterion of the CC is a condition established 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 similar principal stress state. For
instance, it reverts to comparing the opening stress to material
tensile 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:
{\begin{equation}
\underline{\underline{\sigma}} = K r^{\lambda-1} \underline{\underline{s}}(\theta) +
o(r^{\lambda-1}),\label{eq:stress}
\end{equation}}\unskip
where $r$ and $\theta$ are polar coordinates, $\lambda$ is the
characteristic exponent of the singularity and
$\underline{\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~\cite{Leguillon1987}  . 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 (\ref{eq:stress}) 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~\ref{fig:general_problem}. 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 $\Gc S$, where $S$ is the crack surface, must
be balanced by the variation in external force work ($W_{\mathrm{ext}}$)
and in elastic strain energy  ($W_{\mathrm{el}}$) so that: 
{\begin{equation}
\Delta W_{\mathrm{el}} + \Gc S = \Delta W_{\mathrm{ext}}
\label{MA:energy_criterion_first}
\end{equation}}\unskip

When solving the CC, the objective is to determine the initiation crack
surface $S_{\mathrm{c}}$ and initiation imposed loading (for instance the
initiation force $F_{\mathrm{c}}$ based on the example provided in
Figure~\ref{fig:general_problem}) 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
(\ref{MA:energy_criterion_first}) can be expressed as  
{\begin{equation}
\Delta W_{\mathrm{el}} + \Gc \ell = \Delta W_{\mathrm{ext}}
\label{MA:energy_criterion_firstBi}
\end{equation}}\unskip
where $\ell$ is the newly created crack length (a priori unknown). At
the initiation imposed loading $\ell = \ell_{\rm{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_{\mathrm{i}}$. These parameters depend on the initiation
imposed loading and the geometry of the problem. In the problem
represented in Figure~\ref{fig:general_problem}, there is only one
leading term, see (\ref{eq:stress}), and therefore only one initiation
GSIF.  

\subsection{The matched asymptotic approach}\label{ssec:outer_domain}

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. 

\paragraph*{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~\ref{fig:general_problem}, is shown in
Figure~\ref{fig:outer_domain}a.
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~\cite{Leguillon1987} on a closed path surrounding the singular
point (e.g.,~$\Gamma$ in Figure~\ref{fig:outer_domain}a). The GSIF
calculation based on the contour integral can be implemented in
3D~\cite{Doitrand2020IJSS} or in 2D for
isotropic~\cite{Leguillon2002,Leguillon2007} or
anisotropic~\cite{Labossiere1999,Duminy2024} materials, for
multi-material
configurations~\cite{Labossiere1999,Herrera2022,Herrera2022b,Herrera2024},  
or even based on displacement fields measured experimentally by digital
image correlation~\cite{Doitrand2020EFM}. Other approaches also exist
to compute the GSIF, such as the quasidual function
method~\cite{Yosibash2005, Yosibash2021}, least square
fitting~\cite{Dunn1997} or an extraction from the strain energy
density~\cite{Lazzarin2010}. 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~\ref{fig:general_problem}). 

\begin{figure}
\includegraphics{fig02}
\caption{(a) Outer domain where the crack is disregarded, the contour
$\Gamma$ 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.  \label{fig:outer_domain}}
\end{figure}

\paragraph*{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~\ref{fig:general_problem}  is shown in
Figure~\ref{fig:outer_domain}b. 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.

\paragraph*{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 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~\ref{fig:general_problem}) and further apply the CC for
studying crack initiation.

\paragraph*{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~(\ref{eq:stress})) as well as
the elastic strain energy variation due to crack initiation
(Equation~(\ref{MA:energy_criterion_first})) 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. 

\section{Formulation of the approach} \label{sec:Formulation}
The matched asymptotic expansion is used in mechanical engineering to
predict the solution, i.e.~the displacement field
$\underline{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~\ref{fig:examples_matched_crack}. This element is frequently
called perturbation, since it is assumed that its size $\ell$ is small
compared to any dimensions of the structure. 

\begin{figure}
\includegraphics{fig03}
\caption{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.\label{fig:examples_matched_crack}}
\end{figure}

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~\ref{fig:examples_actualdomain}, where the notation of the
problem that is approximated is represented. The domain
$\Omega_{{\ell}}$ has an outer contour $\Gamma =
\Gamma_{\mathrm{V}}\cup \Gamma_{{N}}\cup
\Gamma_{\mathrm{D}}\cup\Gamma_{{\ell}}$. The contour
$\Gamma_{\mathrm{D}}$ is characterized by an imposed displacement
$\bar{U}$, whereas the contours $ \Gamma_{{N}}$,
$\Gamma_{\mathrm{V}}$ and $\Gamma_{{\ell}}$ have a stress-free
boundary conditions. The notation $\Gamma_{\mathrm{V}}$ is referred to
the contour of the V-notch and $\Gamma_{{\ell}}$ to the one of the
small perturbation.

\begin{figure}
\includegraphics{fig04}
\caption{Representation of the notations in (a) the inner domaine and (b)
the outer domain for the example of a cavity ahead of a V-notch. 
\label{fig:examples_actualdomain}}
\end{figure}

Hence, the set of equations that defines the actual problem is 
{\begin{eqnarray}
- \nabla_{{x}} \cdot \, \uuline{\sigma_{{\ell}}} &=& 0 \quad \text{in }\Omega_{{\ell}}, \Seqnsplit
\uuline{\sigma_{{\ell}}} &=& \mathbf{C} : \nabla_{{x}} \underline{U_{{\ell}}} , \Seqnsplit
\uuline{\sigma_{{\ell}}} \cdot \underline{n} &=& \bar{h} \quad \text{on } \Gamma_{{N}},  \Seqnsplit
\uuline{\sigma_{{\ell}}} \cdot \underline{n} &=& 0 \quad \text{on }\Gamma_{\mathrm{V}} \cup \Gamma_{{\ell}},  \Seqnsplit
\underline{U_{{\ell}}} &=& \underline{\bar{U}} \quad \text{on } \Gamma_{\mathrm{D}},
\end{eqnarray}}\unskip
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 $\underline{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. 

\paragraph*{Outer expansion} 
In this approximation, the actual solution is represented as
{\begin{equation}
\underline{U_l}(x_1, x_2) = \underline{U_0}(x_1, x_2) 
+ \cdots \label{eq:outerexpansion}
\end{equation}}\unskip
where $\underline{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~\ref{fig:examples_actualdomain}b represents
$\Omega_0$ associated with the actual domain of
Figure~\ref{fig:examples_actualdomain}a. The second term in
(\ref{eq:outerexpansion}) denoted with an ellipsis is a ``small
correction'' that decreases to $0$ as $\ell \rightarrow 0$. The
solution $\underline{U_0}(x_1, x_2)$ is a good approximation of
$\underline{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 $\underline{U_0}(x_1, x_2)$ is
{\begin{eqnarray}
- \nabla_{{x}} \cdot \, \uuline{\sigma_0} &=& 0 \quad \text{in }\Omega_{0}, \label{eq:outer_equilibrium}\Seqnsplit
\uuline{\sigma_0} &=& \mathbf{C} : \nabla_{{x}} \, \underline{U_0}, \label{eq:outer_constitutive} \Seqnsplit
\uuline{\sigma_0} \cdot \underline{n} &=& \underline{h} \quad \text{on }\Gamma_{{N}}, \label{eq:outer_neumannbcs} \Seqnsplit
\uuline{\sigma_0} \cdot \underline{n} &=& 0 \quad \text{on }\Gamma_{\mathrm{V}},   \Seqnsplit
\underline{U_0} &=& \underline{\bar{U}} \quad \text{on }\Gamma_{\mathrm{D}}. \label{eq:outer_dirichletbcs} 
\end{eqnarray}}\unskip

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~\cite{Leguillon2011}.

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

\begin{figure}
\includegraphics{fig05}
\caption{Scheme of the inner problem.     
\label{fig:examples_innerdomain}}
\end{figure}

The inner expansion is therefore expressed as 
{\begin{equation}
\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 
\label{chap1 eq:innerexpansion}
\end{equation}}\unskip
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{eqnarray*} 
\begin{array}{rcl}
- \nabla_{{y}} \cdot \, \uuline{\tilde{\sigma}_0} &=& 0 \quad \text{in }\Omega_{\mathrm{in}},\\
\uuline{\tilde{\sigma}_0} &=& \mathbf{C} : \nabla_{{y}} \, \underline{V_0} \\ 
\uuline{\tilde{\sigma}_0} \cdot \underline{n} &=& 0 \quad \text{on }\Gamma_{\mathrm{V}}\cup\Gamma_{{\ell}}  \\
\end{array}\qquad
\begin{array}{rcl}
- \nabla_{{y}} \cdot \, \uuline{\tilde{\sigma}_1} &=& 0 \quad \text{in }\Omega_{\mathrm{in}}, \\
\uuline{\tilde{\sigma}_1} &=& \mathbf{C} : \nabla_{{y}} \, \underline{V_1} \\
\uuline{\tilde{\sigma}_1} \cdot \underline{n} &=& 0 \quad \text{on } \Gamma_{\mathrm{V}}\cup\Gamma_{{\ell}} 
\end{array}
\end{eqnarray*}}\unskip
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. 

\paragraph*{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~\ref{fig:examples_actualdomain}, the William's expansion can be
applied, normally expressed in polar coordinates as
{\begin{equation} \label{chap1 eq:expansion}
\underline{U}(r, \theta) = \underline{U}(0,0) + K r^\lambda \underline{u}(\theta) + o(r^\lambda),
\end{equation}}\unskip
assuming that the dominant term is real and have multiplicity one. The
matching conditions can be expressed as
{\begin{eqnarray}
F_0 (\ell) \underline{V_0} (y_1, y_2) &\approx& \underline{U}(0,0), 
\quad \text{when }\rho \rightarrow \infty \label{chap1 eq:matching1} \Seqnsplit
F_1 (\ell) \underline{V_1} (y_1, y_2) &\approx& 
K \ell^\lambda \rho^\lambda \underline{u}(\theta), \quad \text{when }\rho \rightarrow \infty
\label{chap1 eq:matching2} 
\end{eqnarray}}\unskip
where the term $\approx$ means ``behaves like'' and
$\rho=r/\ell=\sqrt{y_1^2+y_2^2}$. It can thus be set:
{\begin{eqnarray}
F_0 (\ell) &=& 1 \quad \text{and} \quad \underline{V_0} (y_1, y_2) \approx \underline{U}(0,0), 
\quad \text{when }\rho \rightarrow \infty \label{chap1 eq:matching1_set} \Seqnsplit
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 \rightarrow \infty \label{chap1 eq:matching2_set} 
\end{eqnarray}}\unskip
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_{\mathrm{in}}$, while it
should decrease to $0$ at infinity to have a finite energy. For this
reason, the superposition principle is applied, 
{\begin{equation}
\underline{V_1} (y_1, y_2) = \rho^\lambda \underline{u}(\theta) + \underline{\hat{V}_1} (y_1, y_2) 
\end{equation}}\unskip
where $\underline{\hat{V}_1} (y_1, y_2)$ is the solution to a well-posed
problem. The set of equations that defines the new term
$\underline{\hat{V}_1} (y_1, y_2)$ is:
{\begin{eqnarray*} \label{chap1 eq:inner_system}
- \nabla_{{y}} \cdot \, \uuline{\hat{\sigma}_1} &=& 0 \quad \text{in } \Omega_{\mathrm{in}}, \\
\uuline{\hat{\sigma}_1} &=& \mathbf{C} : \nabla_{{y}} \, \underline{\hat{V}_1} ,  \\
\uuline{\hat{\sigma}_1} \cdot \underline{n} &=& 0 \quad \text{on }\Gamma_{{\ell}} , \\
\underline{\hat{V}_1} &\approx& 0 \quad \text{as }\rho\rightarrow\infty, \\
\uuline{\hat{\sigma}_1} \cdot \underline{n} &=& -\mathbf{C} : \nabla_{{y}}
\left(\rho\underline{u}(\theta)\right) \cdot \underline{n} \quad \text{on }\Gamma_{\mathrm{V}}
\end{eqnarray*}}\unskip
Note that the terms $\underline{\hat{V}_1} (y_1, y_2)$ and $\underline{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 $\underline{\hat{V}_1} (y_1,
y_2)$ has finite energy at infinity,
{\begin{equation}
\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 \left[ \rho^\lambda \underline{u}(\theta) + \underline{\hat{V}_1} 
(y_1, y_2) \right] + o(\ell^\lambda). \label{chap1 eq:innerexpansion2}
\end{equation}}\unskip


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. 

\section{Numerical implementation of the matched asymptotic approach}
\label{sec:Implementation}
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~\ref{fig:general_problem} 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~\ref{ssec:CC}, 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.

\begin{figure}
\includegraphics{fig06}
\caption{FEM model of the outer domain problem.
\label{fig:outer_domain_mesh}}
\end{figure}

In this sense, the steps of a CC analysis assisted by MA and FEM for
the problem described in Section~\ref{ssec:CC} are the following:
\begin{enumerate}[(1)]
\item 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~\ref{fig:outer_domain_mesh}. The minimum required mesh size can
be determined based on a mesh convergence analysis.
\item 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.
\item The displacements $\underline{u}^{\mathrm{FE}}(\rho,\theta)$ and traction
vector $\underline{t}^{\mathrm{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:
{\begin{equation}
K \approx \frac{\int_\Gamma \underline{t}^{\mathrm{FE}} (\rho,\theta) \cdot 
\rho^\lambda \underline{u}^{-}(\theta) - \underline{t}^{-} (\rho,\theta) 
\cdot \underline{u}^{\mathrm{FE}} (\rho,\theta) {\mathrm{d}} 
\Gamma}{\int_\Gamma \underline{t} (\rho,\theta) \cdot \rho^\lambda \underline{u}^{-} 
(\theta) - \underline{t}^{-} (\rho,\theta) \cdot \rho^\lambda 
\underline{u}(\theta) {\mathrm{d}} \Gamma},
\end{equation}}\unskip
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.

\item Generation of a FE model of the inner domain with crack, see
Figure~\ref{fig:inner_domain_mesh}. 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. 
\item Two versions of this model are computed: 
\begin{enumerate}[(a)]
\item 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.
\item Submodel 1: This submodel corresponds to the state just after the
crack initiation. Crack faces are stress-free.
\end{enumerate}
From these two submodels, the change in potential elastic energy
$\Delta W_{\rm{el}}$ necessary for the evaluation of the energy
criterion (Equation~(\ref{MA:energy_criterion_first})). Several strategies
can be used:
\begin{enumerate}[(a)]
\item 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:
{\begin{equation}
\Delta W_{\mathrm{ext}} - \Delta W_{\rm{el}} = \int_{\Gamma_{\rm{crack}}} \frac{1}{2} 
\underline{t}^{0+}  \cdot (\underline{u}^{1+} - \underline{u}^{1-}) {\mathrm{d}} \Gamma_{\rm{crack}}.
\end{equation}}\unskip
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.
\item 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:
{\begin{equation}
\Delta W_{\mathrm{ext}} - \Delta W_{\rm{el}} = \int_{\Gamma_{\mathrm{ext}}} 
\frac{1}{2} (\underline{t}^{0} \underline{u}^{1} -  \underline{t}^{1} 
\underline{u}^{0} ) {\mathrm{d}} \Gamma_{\mathrm{ext}},
\end{equation}}\unskip
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.
\item 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:
{\begin{equation}
\Delta W_{\mathrm{ext}} - \Delta W_{\rm{el}} = 
\int_{\Gamma_{\rm{N}}} \frac{1}{2}  \underline{t} \left(\underline{u}^{1} - 
\underline{u}^{0} \right) {\mathrm{d}} \Gamma + 
\int_\Omega \frac{1}{2} \left( \uuline{\sigma}^0:\uuline{\varepsilon}^0 - 
\uuline{\sigma}^1:\uuline{\varepsilon}^1  \right){\mathrm{d}}\Omega
\end{equation}}\unskip
where $\Gamma_{\rm{N}}$ corresponds to the boundary with Neumann
boundary condition. The first term would vanish in case of displacement
control.  
\end{enumerate}
\item 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_{\rm{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:
{\begin{equation}
\Delta W_{\mathrm{ext}} - \Delta W_{\rm{el}} = \frac{K^2}{E} \ell^{2 \lambda} A h
\label{eq-Kdirect}
\end{equation}}\unskip
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_{\rm{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_{\rm{b}}/\ell$ or $\alpha$. In next sections more complex problems
will be presented in this sense.

\item Once the value of $K$ is estimated, the coupled criterion of the
finite fracture mechanics can be applied in a quite straightforward
manner:
\begin{itemize}
\item Stress criterion: Combining the first term of the expression in
(\ref{eq:stress}) and the condition outlined for the stress criterion,
this criterion can be expressed as:
{\begin{equation}
K r^{\lambda - 1} s(\theta) \geq \sigma_{\rm{c}} \, \forall~r,
\quad 0 \leqslant r \leqslant \ell
\end{equation}}\unskip
Since for singular cases ($\lambda < 1$) this function is decreasing
with $r$, this condition can also be also written as:
{\begin{equation}
K \ell^{\lambda - 1} s(\theta) \geq \sigma_{\rm{c}}.
\label{MA:SCdef}
\end{equation}}\unskip
\item Energy criterion: Combining the expression in (\ref{eq-Kdirect})
with the condition in (\ref{MA:energy_criterion_first}), the energy
criterion can be expressed as:
{\begin{equation}
\frac{K^2}{E} \ell^{2 \lambda} A h \geq \Gc \ell h
\label{MA:ECdef}
\end{equation}}\unskip
\end{itemize}
The two criteria have to be fulfilled simultaneously. In this case,
since the left term in Equation~(\ref{MA:SCdef}) is a decreasing function of
$\ell$ and the left term in Equation~(\ref{MA:ECdef}) is an increasing
function of $\ell$, the minimum value of $K$ for which the two
conditions are fulfilled is given solving the conditions as two
equations with two unknowns $\ell$ and $K$. The solution, for a
critical length $\ell_{\mathrm{i}}$, is:
{\begin{equation}
\ell_{\mathrm{i}} = \frac{\Gc E s(\theta)^2}{A \sigma_{\rm{c}}^2}
\end{equation}}\unskip
\end{enumerate}
and for the critical GSIF $K_{\mathrm{i}}$,
{\begin{equation}
    K_{\mathrm{i}} = \left(\frac{ \Gc E}{A} \right)^{\lambda - 1} 
    \left(\frac{\sigma_{\rm{c}}}{s(\theta)} \right)^{2 \lambda - 1}
\end{equation}}\unskip
Since $K_{\mathrm{i}}$ is directly proportional to the external load,
the critical load for crack initiation can be calculated directly from
the critical value for $K_{\mathrm{i}}$.

\begin{figure}
\includegraphics{fig07}
\caption{FEM model of the inner domain problem.
\label{fig:inner_domain_mesh}}
\end{figure}

\section{Implementation in complex cases and examples}
\label{sec:Application}
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~\cite{Yosibash2006}. In this section, we provide a detailed overview
of more complex configurations and examples to thoroughly demonstrate
the full potential of the MA approach.

\subsection{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~\cite{Leguillon2003,Leguillon2007,Priel2008}, including as
length parameters both the V-notch radius $\rV$ 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 $\rV/\ell$,
or in the inner domain where the V-notch radius is $1$ and the
dimensionless initiation crack length is $\ell/\rV$. 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~\cite{Cornetti2010}, a
blunted U-notch~\cite{Picard2006,Doitrand2018IJF}, or a pore crack
initiation~\cite{leguillon2008,Doitrand2021TAFMEC}.

\subsection{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: 
{\begin{equation} \label{chap1 eq:mixedmode}
\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)+\cdots.
\end{equation}}\unskip
This is classically encountered when studying mixed opening and shear
mode crack initiation, for instance at a V-notch in
2D~\cite{Yosibash2006,Priel2008,Garcia2012,Doitrand2020FINEL,Cornetti2013} 
or in 3D~\cite{Yosibash2016,Doitrand2018IJF}. Higher order terms,
such as the T-stress, may also be involved for instance when studying
2D crack deviation~\cite{leguillon2008b} or 3D crack front segmentation
into facets~\cite{Doitrand2023IJF}. 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. 

\subsection{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~\cite{Garcia2012,Tran2012}. The characteristic exponent and
corresponding mode depend on the elastic property contrast and
singularity geometry~\cite{Leguillon2013,Leguillon2013b}. It may even
result in complex characteristic exponents and modes
\cite{leguillon2008Murer}, 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~\cite{Leguillon2003b,Moradi2013,Felger2019,Felger2019b} with
the consideration of residual stresses~\cite{Henninger2008}, crack
initiation in microelectronics
structures~\cite{Brillet2007,Schlosser2007}, bond failure of a SiC--SiC
brazed assembly~\cite{Nguyen2012} and crack deflection in layered
ceramics~\cite{Leguillon2006} or ceramic matrix
composites~\cite{Lacroix2002}.

\subsection{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~\cite{Henninger2007, Murer2010}, 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~\cite{Doitrand2024JTCAM}. 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~\cite{Zghal2018}, 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.~\cite{Abaza2022} 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
$\KIc$ and $\sigmac$ 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~\cite{Campagnolo2016}.

\subsection{3D}
There is no restriction regarding the 3D implementation of the MA
approach~\cite{Leguillon20143D}. The displacement field is based on 3D
William's expansions~\cite{Williams1957,Apel2008} 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~\cite{Leguillon1987,Leguillon1995}. The main difference compared to
the 2D case is related to the dual mode. If the primal mode is
$r^{\lambda} \underline{u}(\theta,\phi)$, then the dual mode is
$r^{-\lambda-1} \underline{u}_{3\mathrm{D}}^{-}(\theta,\varphi)$ in
3D instead of $r^{-\lambda} \underline{u}^{-}_{2\mathrm{D}}(\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~\cite{Leguillon20143D,Doitrand2018IJF,Doitrand2023IJF} or
with respect to the square root of the crack
surface~\cite{Yosibash2016}. Overall, the CC has takes the same form as
in 2D and the same contour integral can be used to compute the
GSIF~\cite{Doitrand2020IJSS}, except that the integration domain is a
2D surface in the volume encompassing the singular point instead of a
$1$D 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~\cite{Doitrand2018IJF} or based on both stress and energy
requirements by minimizing the initiation
loading~\cite{Doitrand2023IJF}. 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~\cite{Yosibash2016,Doitrand2023IJF} or derive the crack
shapes from the stress isocontours~\cite{Doitrand2018IJF}.

\subsection{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~\cite{Leguillon2008damage,Leguillon2017,Doitrand2024IJF}. It also
provided a robust method for small crack detection based on
displacement fields measured by digital image
correlation~\cite{Leguillon2011,Romani2015}. Several works about
ceramic failure were also proposed, such as cracking in layered
ceramics~\cite{Leguillon2015Layer}, platelet-based
ceramics~\cite{Jimenez2023,Jimenez2023residual,Duminy2024} or surface
defects in polycristalline ceramics~\cite{Leguillon2018}. Other works
concerned brittle fracture size effect~\cite{jimenez2021}, strength
anisotropy of 3D printed materials~\cite{RoyXu2019}, multicracking of a
stiff inclusion in a soft matrix~\cite{Quesada2009} or elliptical
hole-induced crack initiation~\cite{FELGER2017621}. 

\section{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.~\cite{Murer2010}. 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~\cite{TANNE201880}, defining a CZM~\cite{WU201820}, or
understanding this regularization length in the context of the
CC~\cite{Doitrand2024IJF}. 

\section*{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.

\section*{Dedication}
The manuscript was written through contributions of all authors. All
authors have given approval to the final version of the manuscript.

\section*{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\'on de Espa\~na (Project
PID2020-117001GB-I00/AEI/10.13039/501100011033).

\CDRGrant[EU]{861061-NEWFRAC}
\CDRGrant[Iberdrola Foundation]{101034297}
\CDRGrant[MICIU]{PID2020-117001GB-I00/AEI/10.13039/501100011033}

\back{}

\printbibliography
\refinput{crmeca20240881-reference.tex}

\end{document}
