# Size effects in the toughening of brittle materials by heterogeneities: a
non-linear analysis of front deformations

Mathias Lebihain  Manish Vasoya  Véronique Lazarus  Navier, Ecole des Ponts,
Univ Gustave Eiffel, CNRS, Marne-la-Vallée, France  Department of Aerospace
Engineering, Texas A&M University, College Station, TX 77843, USA  IMSIA,
ENSTA Paris, CNRS, EDF, CEA, Institut Polytechnique de Paris, Palaiseau,
France

###### Abstract

Traditional computational approaches in simulating crack propagation in
perfectly brittle materials rely on the estimate of stress intensity factors
along the rupture front. This proves highly challenging in 3D when the crack
geometry departs from very specific cases for which analytical solutions are
available, like e.g. the penny-shaped crack geometry. Here, we extend the
first-order theory of Gao and Rice (1987b), and predict the distribution of
the mode I stress intensity factor $K_{\mathrm{I}}$ along the front of a
tensile coplanar crack that is slightly perturbed from a reference penny-
shaped configuration, up to second order in the perturbation amplitude. Our
theory is validated against analytical solutions available for embedded
elliptical cracks, and its range of validity is further assessed using
numerical simulations performed on cosine front perturbations of varying mode
and amplitude. It is then used to develop a homogenization framework for the
toughness of weakly disordered media. The effective toughness and its
fluctuations are bridged quantitatively to the intensity of the toughness
fluctuations and their spatial structure. Our theoretical predictions are
compared to the results of $\sim 1$ million simulations of crack propagation
building on our second-order theory and Fast Fourier Transforms. We show that
the impact of toughness heterogeneities is size-dependent, as they generally
weaken the material when the crack size is lower or comparable to the typical
heterogeneity size, but reinforces it otherwise. It results in an apparent
R-curve behavior of the brittle composite at the macroscale.

###### keywords:

Brittle failure , three-dimensional fracture , stress intensity factor ,
coplanar perturbation , circular crack , effective toughness

††journal: International Journal of Solids and Structures

Second-order variation of the SIF for small perturbations of a circular crack
is derived.

It is used to get the effective toughness of heterogeneous materials in the
weak pinning regime.

This differs from the average of the toughness distribution, predicted by the
linear theory.

The effective toughness depends on the crack versus heterogeneity size.

From weaker than the average for small cracks, it becomes stronger for large
ones, suggesting a R-curve behavior.

##  1 Introduction

Despite the rise of powerful numerical methods, like e.g. cohesive zone models
(CZMs) (Barenblatt, 1962; Xu and Needleman, 1994; Camacho and Ortiz, 1996) or
phase-field models (PFMs) (Francfort and Marigo, 1998; Hakim and Karma, 2005;
Bourdin et al., 2008) that can compute spontaneous three-dimensional crack
initiation and propagation in complex structures and multiphysics
environments, more “traditional” simulation techniques based on the Linear
Elastic Fracture Mechanics (LEFM) like e.g. the eXtended Finite Element Method
(XFEM) (Moës et al., 1999), have still some value in modeling crack
propagation in elaborate settings (Paul et al., 2018; Lebihain et al., 2021).
This is particularly true when dealing with sharp discontinuities of material
properties, as the characteristic length scale associated with the failure
process in CZMs and PFMs may interact with that of the fluctuations. In that
case, numerical simulations may fail in reproducing the predictions of LEFM,
even in simple situations of crack deflection at a bi-material interface
(Henry, 2019), and ah-hoc numerical compensation methods are required (Hansen-
Dörr et al., 2020). Meanwhile, LEFM-based methods, which build on the _small-
scale yielding_ assumption, do not suffer these limitations due to the innate
localized nature of the dissipation (Lebihain et al., 2020). Yet, in these
models, crack extension is not spontaneous, but must be predicted from the
combination of ad-hoc propagation (Griffith, 1921; Irwin, 1958) and direction
(Erdogan and Sih, 1963; Hussain et al., 1974; Gol’dstein and Salganik, 1974)
criteria. They often require the a priori knowledge of the stress intensity
factors (SIFs) along the whole crack front. Stress intensity factors are known
analytically for a limited number of simple crack geometries and loading
conditions (Tada et al., 2000), and they can be computed numerically from e.g.
the J-integral (Rice, 1968), or the virtual crack extension (also called
G-theta) method (Destuynder et al., 1981; deLorenzi, 1982). However, when a
crack propagates in a 3D composite material, its front is locally pinned by
material heterogeneities, and adopts a tortuous shape that departs from
“‘standard” geometries (Gao and Rice, 1989). As a result, the SIF distribution
along the front cannot be inferred from exact analytical solutions, and the
computational cost of SIF evaluation often proves too high to model 3D
propagation from the scale of the smallest heterogeneity to that of the
structure.

A good compromise can be found by resorting to perturbative approaches of LEFM
(Lazarus, 2011), in which the stress intensity factors along a front that is
slightly distorted within its plane from a reference configuration can be
estimated from the knowledge of that in the unperturbed configuration.
Building on Bueckner-Rice’s weight function theory (Bueckner, 1987; Rice,
1989), Rice (1985) derived the mode I stress intensity factor distribution
along a semi-infinite planar crack whose front is slightly perturbed from its
reference straight configuration, up to first order in the perturbation
amplitude. The approach was then extended to (i) mixed mode conditions by Gao
and Rice (1986), (ii) other crack geometries, like e.g. the penny-shaped crack
(Gao and Rice, 1987b; Borodachev, 2007), the circular connection (Gao and
Rice, 1987a), tunnel crack (Leblond et al., 1996; Lazarus and Leblond, 2002a,
b), interaction of two tunnel cracks (Pindra et al., 2010; Legrand and
Leblond, 2010), the crack lying between two plates (Legrand et al., 2011),
(iii) out-of-plane perturbations (Movchan et al., 1998; Leblond et al., 2011),
(iv) dynamic ruptures (Perrin and Rice, 1994; Ramanathan and Fisher, 1997),
and (v) cohesive materials (Lebihain et al., 2022). The unparalleled
computational performances offered by the perturbative methods, through the
use of the Fast Fourier Transform (FFT), allowed to investigate the
propagation of a crack front in heterogeneous media (Gao and Rice, 1989), its
roughening under the action of material disorder (Schmittbuhl et al., 1995;
Rosso and Krauth, 2002), and the intermittent statistics that emerge from its
interaction with material heterogeneities (Ponson and Bonamy, 2010; Laurson et
al., 2010). It has given quantitative means to rationalize the crack-inclusion
interactions observed during failure experiments in patterned (Bower and
Ortiz, 1991; Dalmas et al., 2009; Chopin et al., 2011; Patinet et al., 2013a)
and disordered materials (Delaplace et al., 1999; Barés et al., 2018). It has
also provided basic ingredients to formulate a homogenization framework for
coplanar (Roux et al., 2003; Patinet et al., 2013b; Démery et al., 2014a, b)
and non-coplanar (Lebihain et al., 2021) brittle failure. Note that
perturbative approaches have also been used in modeling contact and adhesion
along homogeneous and heterogeneous surfaces (Adda-Bedia and Mahadevan, 2006;
Xia et al., 2012; Budzik and Jensen, 2014; Xia et al., 2015; Argatov, 2021;
Sanner and Pastewka, 2022).

However, the experimental study of Vasoya et al. (2016b) suggests that
_higher-order_ theories may help to reconcile theoretical predictions with
experimental observations of crack front deformation. Moreover, material
toughening by heterogeneities has been shown to be of _second-order_ in the
amplitude of the toughness fluctuations (Patinet et al., 2013b; Démery et al.,
2014b), which is of the same order as the terms omitted by these authors in
their asymptotic expansion of the stress intensity factor. As a result, their
models only draw conclusions on how material toughening scales with the
amplitude of the fluctuations, but they may fail in predicting it
quantitatively. Higher-order theories for the crack front perturbations are
thus crucial to improve our understanding of the interactions between cracks
and heterogeneities. Yet, up to this day, second-order perturbative approaches
have only been derived for the semi-infinite crack (Adda-Bedia et al., 2006;
Leblond et al., 2012; Vasoya et al., 2013, 2016b) loaded in tensile mode I.
Their extension to other geometries (like e.g. the penny-shaped crack or the
circular connection) may bring valuable insights on a wide variety of physical
problems, like e.g. the onset of the fingering instability occurring during
crack propagation (Vasoya et al., 2016a) or during the contact between an
indenter and soft elastic films (Yu and Jiang, 2021), the adhesion hysteresis
observed along chemically heterogeneous interfaces (Sanner and Pastewka,
2022), or even the shape of fluid-driven shear ruptures along frictional
interfaces governed by Coulomb’s friction (Sáez et al., 2022).

In this work, we extend the linear model of Gao and Rice (1987b) for the
penny-shaped crack loaded in tensile mode I to second-order, following the
reasoning of Leblond et al. (2012). We build next on our non-linear theory to
investigate material reinforcement arising from the presence of randomly
arranged heterogeneities of toughness, and focus in particular on the size
effects emerging from the finiteness of the penny-shaped crack geometry.

The paper is organized as follows: in Section 2, we derive a non-linear theory
that predicts the mode I stress intensity factor variations arising from
coplanar perturbations of the crack front from its reference circular
configuration, up to second-order in the perturbation amplitude. This is
performed by deriving first the linear expansion of the fundamental kernel
associated to the penny-shaped crack geometry. In Section 3, we validate our
theory by comparing its predictions with analytical solutions obtained for
elliptical cracks loaded by a uniform tensile stress. The validity range of
the second-order model is then assessed, building on numerical SIF estimates
along wavy crack fronts obtained by the method of Lazarus (2003) and David and
Lazarus (2022) that computes the SIF variations at all orders. Section 4 is
devoted to the analysis of quasi-static crack propagation in random fields of
toughness, with a particular focus on the toughening arising from these
fluctuations. Namely, we discuss how (i) accounting for non-linear terms in
the SIF expansion and (ii) the finiteness of the penny-shaped crack geometry
may change our results with respect to those available in the literature for
the semi-infinite crack (Patinet et al., 2013b; Démery et al., 2014a).

##  2 Derivation of the second-order theory

###  2.1 Generalities

We consider a planar penny-shaped crack embedded in an infinite medium made of
some isotropically linearly elastic material of Young’s modulus $E$ and
Poisson’s ratio $\nu$. The crack front $\mathcal{F}_{0}$ describes a circle of
radius $a_{0}$ centered on a point $O$. The crack is contained within the
$xOz$ plane. The pair $(r,\theta)$ denotes the coordinates of a point of
observation $M$ in the polar coordinate system centered in $O$, and
$(\mathbf{e_{r}},\mathbf{e_{\theta}})$ the polar coordinate basis vectors at
this point. The crack is loaded in pure mode I through some system of forces,
giving rise to a distribution $K_{\mathrm{I}}^{0}(a_{0},\theta)$ of the
unperturbed stress intensity factor along the reference circular crack front
$\mathcal{F}_{0}$.

The crack front is now perturbed within its plane in the direction
$\mathbf{e_{r}}$ normal to its reference configuration by a small amount:

$\delta a(\theta)=A\phi(\theta),$ (1)

where $A$ is an infinitesimally small parameter quantifying the non-
circularity of the perturbed front $\mathcal{F}$, and $\phi$ is a shape
function (see Fig. 1). $\delta K_{\mathrm{I}}$ denotes the variations of the
mode I SIF arising from the perturbation $\delta a$. Evaluated at a point
$M_{1}$ of $\mathcal{F}$ indexed by $\theta_{1}$, $\delta K_{\mathrm{I}}$
reads (Panasyuk, 1962; Gao and Rice, 1987b):

$\delta K_{\mathrm{I}}(\theta_{1})=\dfrac{\partial
K_{\mathrm{I}}^{0}}{\partial a}(a_{0},\theta_{1})\delta
a(\theta_{1})+\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)K_{\mathrm{I}}^{0}(a_{0},\theta)\left[\delta
a(\theta)-\delta a(\theta_{1})\right]a_{0}d\theta+\mathcal{O}(A^{2}),$ (2)

at first order in the perturbation $\delta a$. Here, the symbol $\mathrm{PV}$
denotes a Cauchy principal value, and $Z^{0}(a_{0};\theta_{1},\theta_{2})$ is
the _kernel_ associated with the reference circular crack of radius $a_{0}$,
evaluated at any points $M^{0}_{1}$ and $M^{0}_{2}$ of $\mathcal{F}^{0}$,
indexed by $\theta_{1}$ and $\theta_{2}$ respectively (see Fig. 1). It takes
the general form (Rice, 1989):

$Z^{0}(a_{0};\theta_{1},\theta_{2})=\dfrac{1}{2\pi}\dfrac{W^{0}(a_{0};\theta_{1},\theta_{2})}{\mathcal{D}(M^{0}_{1},M^{0}_{2})^{2}},$
(3)

where $\mathcal{D}(M^{0}_{1},M^{0}_{2})$ is the Euclidean distance between the
points $M^{0}_{1}$ and $M^{0}_{2}$, and $W^{0}$ is the _fundamental kernel_
111Note the subtle difference of vocabulary: _fundamental kernel_ standing for
the dimensionless quantity $W$ and _kernel_ for $Z$ appearing as the kernel of
the integral in Eq. 2. associated with the reference crack front geometry. For
a circular crack front, $W^{0}(a_{0};\theta_{1},\theta_{2})=1$, so that
$Z^{0}$ reads:

$Z^{0}(a_{0};\theta_{1},\theta_{2})=\dfrac{1}{8\pi
a_{0}^{2}\sin^{2}\left[(\theta_{2}-\theta_{1})/2\right]}.$ (4)

Our goal here is to go beyond the linear expansion (2) of Gao and Rice
(1987b), and derive the expression of the perturbed SIF $K_{\mathrm{I}}$ at
_second order_ in the perturbation $\delta a$. We follow here the approach of
Leblond et al. (2012), and calculate first the variations of the kernels $W$
and $Z$ at _first order_ in $\delta a$.

Figure 1: The reference penny-shaped crack with a circular front
$\mathcal{F}_{0}$ of radius $a_{0}$ (in dashed black line) is perturbed within
its plane $(xOz)$ by a small amount $\delta a(\theta)$ in the direction
$\mathbf{e_{r}}$, giving rise to a perturbed crack front $\mathcal{F}$
parametrized by $a(\theta)=a_{0}+\delta a(\theta)$ (in solid red line). The
points $M_{1}$, $M_{2}$, and $M$, indexed by $\theta_{1}$, $\theta_{2}$ and
$\theta$ respectively, are located on the perturbed front $\mathcal{F}$.
$M^{0}_{1}$, $M^{0}_{2}$ and $M^{0}$ are their orthogonal projections on
$\mathcal{F}_{0}$.

###  2.2 First-order expansion of the fundamental kernel

For tensile cracks, the first-order expansion of the fundamental kernel $W$
has been derived by Rice (1989). Provided that the perturbation $\delta a$
vanishes at $M_{1}$ and $M_{2}$ indexed by $\theta_{1}$ and $\theta_{2}$, the
variations $\delta W$ of the fundamental kernel write as:

$\dfrac{\delta
W(\theta_{1},\theta_{2})}{\mathcal{D}(M_{1},M_{2})^{2}}=\dfrac{1}{2\pi}\mathrm{PV}\int_{0}^{2\pi}\dfrac{W^{0}(a_{0};\theta_{1},\theta)}{\mathcal{D}(M_{1},M^{0})^{2}}\dfrac{W^{0}(a_{0};\theta_{2},\theta)}{\mathcal{D}(M_{2},M^{0})^{2}}\delta
a(\theta)\,a_{0}d\theta+\mathcal{O}(A^{2}),$ (5)

where $M^{0}$ is the point indexed by $\theta$ on $\mathcal{F}_{0}$ .

In the general case where $\delta a$ does not vanish at $M_{1}$ and $M_{2}$,
one must find a perturbation $\delta a_{*}$ (i) for which the variations of
the fundamental kernel $\delta W^{*}$ are known a priori, and (ii) that takes
the same values as $\delta a$ at two points $M_{1}$ and $M_{2}$ (see Fig. 2):

$\delta a_{*}(\theta_{1})=\delta a(\theta_{1})\text{, and }\delta
a_{*}(\theta_{2})=\delta a(\theta_{2}).$ (6)

One then has (Favier et al., 2006a):

$\dfrac{\delta
W(\theta_{1},\theta_{2})}{\mathcal{D}(M^{0}_{1},M^{0}_{2})^{2}}=\dfrac{\delta
W^{*}(\theta_{1},\theta_{2})}{\mathcal{D}(M^{0}_{1},M^{0}_{2})^{2}}+\dfrac{1}{2\pi}\mathrm{PV}\int_{0}^{2\pi}\dfrac{W^{0}(a_{0};\theta_{1},\theta)}{\mathcal{D}(M^{0}_{1},M^{0})^{2}}\dfrac{W^{0}(a_{0};\theta_{2},\theta)}{\mathcal{D}(M^{0}_{2},M^{0})^{2}}\left[\delta
a(\theta)-\delta a_{*}(\theta)\right]a_{0}d\theta+\mathcal{O}(A^{2}),$ (7)

where $M^{0}_{1}$ and $M^{0}_{2}$ are the locations defined by the orthogonal
projection on $\mathcal{F}_{0}$ along the direction $\mathbf{e_{r}}$ of the
points $M_{1}$ and $M_{2}$, i.e. their position before the application of the
transformation $\delta a$ (see Fig. 1). The variation of the kernel $Z$ then
reads:

$\displaystyle\delta Z(\theta_{1},\theta_{2})$
$\displaystyle=\dfrac{W^{0}(a_{0};\theta_{1},\theta_{2})+\delta
W(\theta_{1},\theta_{2})}{2\pi\mathcal{D}(M_{1},M_{2})^{2}}-\dfrac{W^{0}(a_{0};\theta_{1},\theta_{2})}{2\pi\mathcal{D}(M^{0}_{1},M^{0}_{2})^{2}}$
(8)
$\displaystyle=Z^{0}(a_{0};\theta_{1},\theta_{2})\left[\left(\dfrac{\mathcal{D}(M_{1},M_{2})}{\mathcal{D}(M^{0}_{1},M^{0}_{2})}\right)^{-2}-1\right]+\dfrac{\delta
W^{*}(\theta_{1},\theta_{2})}{2\pi\mathcal{D}(M^{0}_{1},M^{0}_{2})^{2}}$
$\displaystyle\hskip
8.53581pt+\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)Z^{0}(a_{0};\theta_{2},\theta)\left[\delta
a(\theta)-\delta a_{*}(\theta)\right]a_{0}d\theta+\mathcal{O}(A^{2}),$

at first order in $A$. In the case of the semi-infinite crack, the first term
of $\delta Z$ is of second order in $A$, and reduces thus to zero in the
first-order calculations of Leblond et al. (2012). This is not the case for a
reference circular crack front, for which it writes as :

$Z^{0}(a_{0};\theta_{1},\theta_{2})\left[\left(\dfrac{\mathcal{D}(M_{1},M_{2})}{\mathcal{D}(M^{0}_{1},M^{0}_{2})}\right)^{-2}-1\right]=-\dfrac{\delta
a(\theta_{1})+\delta a(\theta_{2})}{8\pi
a_{0}^{3}\sin^{2}\left[(\theta_{2}-\theta_{1})/2\right]}+\mathcal{O}(A^{2}).$
(9)

The second and third terms of $\delta Z$ in Eq. (8) depend on the actual
choice of the perturbation $\delta a_{*}$. To compute the second term
analytically, the fundamental kernel of the crack front parametrized by
$a^{*}(\theta)=a_{0}+\delta a^{*}(\theta)$ must be known explicitly. For
unbounded solids, Rice (1989) proposed to express $\delta a^{*}$ as the
combination of translations, rotations or scaling, leaving the fundamental
kernel unchanged, i.e. $\delta W^{*}=0$. Such a transformation has been used
by Leblond et al. (2012) to derive the first-order variation of the
fundamental kernel $W$ for the semi-infinite coplanar crack loaded in tensile
mode I. However, one can show that, in the case investigated here, this
transformation contains a rotation so that the perturbation $\delta a_{*}$
cannot be directly expressed as a normal extension in the direction
$\mathbf{e_{r}}$.

Here, we follow another route and take $\delta a_{*}$ as the orthogonal
projection in the direction $\mathbf{e_{r}}$, on the reference crack front
$\mathcal{F}_{0}$, of another circular crack front $\mathcal{F}^{*}$ that goes
by $M_{1}$ and $M_{2}$ (see Fig. 2). $\mathcal{F}^{*}$ corresponds to
$\mathcal{F}_{0}$ dilated by a factor $(1+\Delta a_{*}/a_{0})$, and translated
by an amount $\Delta r_{*}$ in the direction $\pi/2+\theta_{*}$. The perturbed
crack front $\mathcal{F}^{*}$ is thus a circle of radius $a_{0}+\Delta a_{*}$,
centered on the point $O^{*}$ of Cartesian coordinates $(-\Delta
r_{*}\sin\theta_{*},\Delta r_{*}\cos\theta_{*})$. The associated perturbation
$\delta a_{*}$ can be expressed as:

$\delta a_{*}(\theta)=\Delta r_{*}\sin(\theta-\theta_{*})+\sqrt{(a_{0}+\Delta
a_{*})^{2}-{\Delta r_{*}}^{2}\cos^{2}(\theta-\theta_{*})}-a_{0}.$ (10)

Since the crack front $\mathcal{F}^{*}$ is circular, its fundamental kernel
$W^{*}$ is also equal to $1$. As a result, $\delta W^{*}=0$, and the second
term of $\delta Z$ in Eq. (8) vanishes.

Figure 2: The first-order variations of the kernel
$Z^{0}(a_{0};\theta_{1},\theta_{2})$, associated with the reference circular
crack front (in dashed black line), are calculated from a shifted circular
crack front $\mathcal{F}_{*}$ (in dashed blue line) of radius $a_{0}+\Delta
a_{*}$ and centered in $O^{*}$, that passes by $M_{1}$ and $M_{2}$ located on
the perturbed crack front (in red solid line). The transformation from the
reference crack front $\mathcal{F}_{0}$ and shifted crack front
$\mathcal{F}^{*}$ is parametrized by the perturbation $\delta
a_{*}\left(\theta\right)$ defined in Eq. (10).

It leaves us with the calculation of the principal value integral in Eq. (8).
There is an infinite number of circles that go by the points $M_{1}$ and
$M_{2}$. Consequently, $\delta a_{*}$ is a 3-parameters transformation
$(\Delta a_{*},\Delta r_{*},\theta_{*})$ that must satisfy the sole 2
conditions of Eq. (6). One may then carefully choose an appropriate
$\theta_{*}$ so that (i) $\delta a_{*}$ of Eq. (10) is well-defined, and (ii)
one can compute the third term of $\delta Z$ analytically. Here, we take:

$\theta_{*}=\dfrac{\theta_{1}+\theta_{2}}{2}.$ (11)

Note that, in equation (8), one only needs an expression of $\delta a_{*}$ at
first order in $A$ to compute the variations of $Z$. The conditions (6) then
yield:

$\Delta r_{*}=\dfrac{\delta a(\theta_{2})-\delta
a(\theta_{1})}{2\sin\left[(\theta_{2}-\theta_{1})/2\right]}+\mathcal{O}(A^{2}),\text{and
}\Delta a_{*}=\dfrac{\delta a(\theta_{1})+\delta
a(\theta_{2})}{2}+\mathcal{O}(A^{2}).$ (12)

The perturbation $\delta a_{*}$ can finally be expressed as:

$\displaystyle\delta a_{*}(\theta)=\dfrac{\delta a(\theta_{1})+\delta
a(\theta_{2})}{2}+\dfrac{\delta a(\theta_{2})-\delta
a(\theta_{1})}{2}\cdot\dfrac{\sin\left[\theta-(\theta_{1}+\theta_{2})/2\right]}{\sin\left[(\theta_{2}-\theta_{1})/2\right]}+\mathcal{O}(A^{2})$
(13) $\displaystyle\hskip 2.84526pt=\delta a(\theta_{1})+\dfrac{\delta
a(\theta_{2})-\delta
a(\theta_{1})}{2}\left[\dfrac{\cos\left[(\theta_{2}-\theta_{1})/2\right]}{\sin\left[(\theta_{2}-\theta_{1})/2\right]}\sin(\theta-\theta_{1})+2\sin^{2}[(\theta-\theta_{1})/2]\right]+\mathcal{O}(A^{2})$
$\displaystyle\hskip 2.84526pt=\delta a(\theta_{2})+\dfrac{\delta
a(\theta_{2})-\delta
a(\theta_{1})}{2}\left[\dfrac{\cos\left[(\theta_{2}-\theta_{1})/2\right]}{\sin\left[(\theta_{2}-\theta_{1})/2\right]}\sin(\theta-\theta_{2})-2\sin^{2}[(\theta-\theta_{2})/2]\right]+\mathcal{O}(A^{2}).$

The second and third expressions of $\delta a_{*}(\theta)$ in Eq. (13) only
involves $(\theta-\theta_{1})$ and $(\theta-\theta_{2})$ respectively, thanks
to the careful choice of $\theta_{*}$. They are analogous to Leblond et al.
(2012)’s Eq. (7), and prove key in the evaluation of the principal value
integral of our Eq. (8) (see A for details).

Next, we use another trick, analogous to the partial fraction decomposition of
Leblond et al. (2012)’s Eq. (9), and propose the following decomposition:

$\displaystyle\dfrac{1}{\sin^{2}[(\theta-\theta_{1})/2]\sin^{2}[(\theta-\theta_{2})/2]}=$
$\displaystyle\dfrac{1}{\sin^{2}\left[(\theta_{2}-\theta_{1})/2\right]}\left[\dfrac{\cos^{2}[(\theta-\theta_{1})/2]}{\sin^{2}[(\theta-\theta_{1})/2]}+\dfrac{\cos^{2}[(\theta-\theta_{2})/2]}{\sin^{2}[(\theta-\theta_{2})/2]}\right]$
(14)
$\displaystyle+\dfrac{2\cos[(\theta_{2}-\theta_{1})/2]}{\sin^{3}[(\theta_{2}-\theta_{1})/2]}\left[\dfrac{\cos^{3}[(\theta-\theta_{1})/2]}{\sin[(\theta-\theta_{1})/2]}-\dfrac{\cos^{3}[(\theta-\theta_{2})/2]}{\sin[(\theta-\theta_{2})/2]}\right]$
$\displaystyle+\dfrac{2}{\sin^{2}\left[(\theta_{2}-\theta_{1})/2\right]}\left[\cos^{2}[(\theta-\theta_{1})/2]+\cos^{2}[(\theta-\theta_{2})/2]\right]$

Combining Eqs. (8), (13-14), and after some simplifications, one gets:

$\displaystyle\begin{aligned} \delta
Z~{}(\theta_{1},\theta_{2})&=-\dfrac{5}{32\pi a_{0}^{3}}\dfrac{\delta
a(\theta_{1})+\delta
a(\theta_{2})}{\sin^{2}\left[(\theta_{2}-\theta_{1})/2\right]}+\dfrac{1}{(8\pi)^{2}a_{0}^{3}\sin^{2}\left[(\theta_{2}-\theta_{1})/2\right]}\int_{0}^{2\pi}2\left[\cos^{2}\left(\dfrac{\theta-\theta_{1}}{2}\right)+\cos^{2}\left(\dfrac{\theta-\theta_{2}}{2}\right)\right]\delta
a(\theta)d\theta\\\
&+\dfrac{1}{(8\pi)^{2}a_{0}^{3}\sin^{2}\left[(\theta_{2}-\theta_{1})/2\right]}\times\Bigg{\\{}\\\
&\hskip
28.45274pt\left.\mathrm{PV}\int_{0}^{2\pi}\dfrac{\cos^{2}[(\theta-\theta_{1})/2]}{\sin^{2}[(\theta-\theta_{1})/2]}\left[\delta
a(\theta)-\delta
a(\theta_{1})\right]d\theta+\mathrm{PV}\int_{0}^{2\pi}\dfrac{\cos^{2}[(\theta-\theta_{2})/2]}{\sin^{2}[(\theta-\theta_{2})/2]}\left[\delta
a(\theta)-\delta a(\theta_{2})\right]d\theta\right.\\\ &\hskip
28.45274pt\left.+\dfrac{2\cos[(\theta_{2}-\theta_{1})/2]}{\sin^{3}[(\theta_{2}-\theta_{1})/2]}\cdot\mathrm{PV}\int_{0}^{2\pi}\left[\dfrac{\cos^{3}[(\theta-\theta_{1})/2]}{\sin[(\theta-\theta_{1})/2]}-\dfrac{\cos^{3}[(\theta-\theta_{2})/2]}{\sin[(\theta-\theta_{2})/2]}\right]\delta
a(\theta)d\theta\right\\}+\mathcal{O}(A^{2})\end{aligned}$

(15)

Details on the derivation of Eq. (15) are given in A. Note that the two first
terms of Eq. (15) typically emerge from the finiteness of the crack geometry,
and go to zero in the limit $a_{0}\rightarrow+\infty$. One can show that
Leblond et al. (2012)’s Eq. (10) for the semi-infinite crack is retrieved in
the same limit.

We derived here the first-order variation of the kernel $Z$ for a penny-shaped
crack loaded in mode I. However, the determination of $\delta a_{*}$ provided
here enables extending our results to mixed mode I+II+III (Favier et al.,
2006a), to other geometries like e.g. the circular connection (Gao and Rice,
1987a), or to circular cracks propagating in cohesive materials (Lebihain et
al., 2022).

###  2.3 Second-order expansion of the SIF

We consider again a crack front $\mathcal{F}$ that has been perturbed from a
circular configuration of radius $a_{0}$ by an amount $\delta
a(\theta)=A\phi(\theta)$, where $A$ is now a _small_ albeit _not
infinitesimal_ parameter. The stress intensity factor $K_{\mathrm{I}}$ and the
kernel $Z$ along the perturbed crack front $\mathcal{F}$ can be expanded at
several orders in $A$. Their expansion writes:

${\color[rgb]{0,0,0}\begin{cases}\begin{aligned}
K_{\mathrm{I}}(a_{0},\left[A\phi\right];\theta_{1})=K_{\mathrm{I}}^{0}(a_{0},\theta_{1})&+A\,K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})\\\
&+A^{2}K_{\mathrm{I}}^{2}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})+\mathcal{O}(A^{3})\end{aligned}\\\
Z(a_{0},\left[A\phi\right];\theta_{1},\theta_{2})=Z^{0}(a_{0};\theta_{1},\theta_{2})+AZ^{1}(a_{0},\left[\phi\right];\theta_{1},\theta_{2})+\mathcal{O}(A^{2}),\end{cases}}$
(16)

where $K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})$ is
the first-order variation of the stress intensity factor induced at
$\theta_{1}$ by the perturbation $\phi$ applied to a reference circular
configuration of radius $a_{0}$ and SIF $K_{\mathrm{I}}^{0}$, and
$K_{\mathrm{I}}^{2}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})$ its second-
order variation. Similarly,
$Z^{1}(a_{0},\left[\phi\right];\theta_{1},\theta_{2})$ is the first-order
variation of the kernel $Z^{0}$ induced at $(\theta_{1},\theta_{2})$ by the
perturbation $\phi$ from the reference configuration of radius $a_{0}$.
$K_{\mathrm{I}}^{1}$ and $Z^{0}$ have been derived by Gao and Rice (1987b),
and are recalled in Eqs. (2-3). $Z^{1}$ has been derived here in Eq. (15). Our
goal here is to provide an analytical expression for $K_{\mathrm{I}}^{2}$.

To this end, we follow the steps of Leblond et al. (2012), and superimpose to
the primary perturbation $\delta a$ a _proportional_ yet _infinitesimal_
secondary perturbation $\delta A\phi$, where $\delta A$ is an infinitesimally
small parameter. By definition, one has:

$\displaystyle K_{\mathrm{I}}(a_{0},\left[(A+\delta
A)\phi\right];\theta_{1})-K_{\mathrm{I}}(a_{0},\left[A\phi\right];\theta_{1})$
$\displaystyle=\dfrac{\partial K_{\mathrm{I}}}{\partial
A}(a_{0},\left[A\phi\right];\theta_{1})\delta A+\mathcal{O}(\delta A^{2})$
(17)
$\displaystyle=\left\\{K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})+2AK_{\mathrm{I}}^{2}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})+\mathcal{O}(A^{2})\right\\}\delta
A+\mathcal{O}(\delta A^{2})$

Alternatively, one may consider the front $\mathcal{F^{\prime}}$ of reference
radius $a^{\prime}_{0}=a_{0}+\delta A\phi(\theta_{1})$ perturbed by the
quantity $\delta a(\theta)=A\phi(\theta)$, and perturbed it further by a
infinitesimal quantity $\delta A\phi(\theta)-\delta A\phi(\theta_{1})$. Yet,
this secondary perturbation is not performed in the direction
$\mathbf{e_{r}}$, but rather in the direction normal to the perturbed front
$\mathcal{F^{\prime}}$. Given that the error in the resulting front position
is of second order in both in $\delta A$ (see Leblond et al. (2012))222Note
that one must ensure that the secondary perturbation is zero at
$\theta=\theta_{1}$, otherwise the transformation induces a phase-shift that
is first-order in $\delta A$. This was overlooked by Leblond et al. (2012) and
Vasoya et al. (2013), but their results remain valid as these first-order
corrections are zero when $K_{\mathrm{I}}^{0}$ is independent of the position
along the front., one has (Rice, 1989):

$\displaystyle K_{\mathrm{I}}(a_{0},\left[(A+\delta
A)\phi\right];\theta_{1})-K_{\mathrm{I}}(a_{0},\left[A\phi\right];\theta_{1})=K_{\mathrm{I}}(a_{0}+\delta
A\phi(\theta_{1}),\left[A\phi\right];a_{0};\theta_{1})$ (18)
$\displaystyle\hskip 56.9055pt+\mathrm{PV}\int_{0}^{2\pi}Z(a_{0}+\delta
A\phi(\theta_{1}),\left[A\phi\right];\theta_{1},\theta)K_{\mathrm{I}}(a_{0}+\delta
A\phi(\theta_{1}),\left[A\phi\right];\theta)\left[\delta A\phi(\theta)-\delta
A\phi(\theta_{1})\right]ds(\theta)$ $\displaystyle\hskip 56.9055pt-
K_{\mathrm{I}}(a_{0},\left[A\phi\right];\theta_{1})+\mathcal{O}(\delta
A^{2}),$

where $s$ is the curvilinear abscissa along the perturbed front
$\mathcal{F^{\prime}}$. Equating Eqs. (17) and (18), one may identify the
terms of zero order in $A$ in the expression of $\partial
K_{\mathrm{I}}/\partial A$:

$K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})=\dfrac{\partial
K_{\mathrm{I}}^{0}}{\partial
a}(a_{0},\theta_{1})\phi(\theta_{1})+\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)K_{\mathrm{I}}^{0}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]a_{0}d\theta,$
(19)

which is Eq. (2). Identifying the terms of first order in $A$ in the
expression of $\partial K_{\mathrm{I}}/\partial A$, one has after some
simplifications:

$\displaystyle
K_{\mathrm{I}}^{2}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})=$
$\displaystyle\frac{1}{2}\dfrac{\partial^{2}K_{\mathrm{I}}^{0}}{\partial
a}(a_{0},\theta_{1})\phi(\theta_{1})^{2}+\frac{1}{2}\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)\dfrac{\partial
K_{\mathrm{I}}^{0}}{\partial
a}(a_{0},\theta)\left[\phi(\theta)^{2}-\phi(\theta_{1})^{2}\right]a_{0}d\theta$
(20)
$\displaystyle+\frac{1}{2}\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)K_{\mathrm{I}}^{0}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]^{2}a_{0}d\theta$
$\displaystyle+\frac{1}{2}\mathrm{PV}\int_{0}^{2\pi}Z^{1}(a_{0},\phi;\theta_{1},\theta)K_{\mathrm{I}}^{0}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]a_{0}d\theta$
$\displaystyle+\frac{1}{2}\mathrm{PV}\int_{0}^{2\pi}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)Z^{0}(a_{0};\theta,\theta^{\prime})K_{\mathrm{I}}^{0}(a_{0},\theta^{\prime})\left[\phi(\theta)-\phi(\theta_{1})\right]\left[\phi(\theta^{\prime})-\phi(\theta)\right]a_{0}^{2}d\theta
d\theta^{\prime}.$

Details on the derivation of Eq. (20) can be found in B. Equations (19) and
(20) are general equations that prove unfit for analytical calculations, and
one has to resort to Fourier series to evaluate them. We shall do it here in
the simplified case where the crack is loaded through an axisymmetric system
of forces, i.e. when $K_{\mathrm{I}}^{0}$ is independent of $\theta$. Note
however that this is enough to derive second-order analytical solutions for
crack propagation in heterogeneous media, even when $K_{\mathrm{I}}^{0}$
actually depends on $\theta$, as these fluctuations are already of second
order in the first-order theory of Gao and Rice (1987b). However, for the sake
of completeness, the derivation for a non-axisymmetric $K_{\mathrm{I}}^{0}$ is
given in B.

Next, we assume that $K_{\mathrm{I}}^{0}$ only depends on the radius $a_{0}$
of the reference penny-shaped crack. One can decompose $\phi$ in Fourier
series. It reads:

$\phi(\theta)=\sum_{k=-\infty}^{+\infty}\hat{\phi}_{k}e^{ik\theta}\text{,
where
}\hat{\phi}_{k}=\dfrac{1}{2\pi}\int_{0}^{2\pi}\phi(\theta)e^{-ik\theta}d\theta.$
(21)

Then, the first-order contribution of Eq. (19) writes as:

$K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})=\sum_{k=-\infty}^{+\infty}\left[\dfrac{dK_{\mathrm{I}}^{0}}{da}(a_{0})-\dfrac{|{k}|}{2a_{0}}K_{\mathrm{I}}^{0}(a_{0})\right]\hat{\phi}_{k}e^{ik\theta_{1}},$
(22)

and the second-order contribution of Eq. (20) reads:

$\displaystyle
K_{\mathrm{I}}^{2}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})=\sum_{k=-\infty}^{+\infty}\sum_{k^{\prime}=-\infty}^{+\infty}\left[\dfrac{1}{2}\dfrac{d^{2}K_{\mathrm{I}}^{0}}{da^{2}}(a_{0})-\dfrac{1}{4a_{0}}\dfrac{dK_{\mathrm{I}}^{0}}{da}(a_{0})|{k+k^{\prime}}|+\dfrac{R(k,k^{\prime})}{16a_{0}^{2}}K_{\mathrm{I}}^{0}(a_{0})+\dfrac{|{k}|+|{k^{\prime}}|}{4a_{0}^{2}}K_{\mathrm{I}}^{0}(a_{0})\right]\hat{\phi}_{k}\hat{\phi}_{k^{\prime}}e^{i(k+k^{\prime})\theta_{1}},$

(23)

where
$R(k,k^{\prime})=2|{k+k^{\prime}}|(|{k}|+|{k^{\prime}}|)-(k+k^{\prime})^{2}-k^{2}-{k^{\prime}}^{2}$
is twice the function $F$ of Leblond et al. (2012), written in the form
simplified by Vasoya et al. (2016b). Note that Eq. (20) differs from Eq. (8)
of Vasoya et al. (2013) only by the term
$(|{k}|+|{k^{\prime}}|)K_{\mathrm{I}}^{0}(a_{0})/4a_{0}^{2}$ that typically
emerges from the finiteness of the penny-shaped crack geometry. As expected,
this term goes to $0$ when $a_{0}\rightarrow+\infty$, and one retrieves, in
this limit, the results of Leblond et al. (2012) and Vasoya et al. (2013) for
the semi-infinite crack.

One can alternatively express Eq. (16a) in terms of the Fourier coefficient
$\widehat{\delta a}_{k}$ of the perturbation $\delta a$:

$K_{\mathrm{I}}(a_{0},\left[\delta
a\right];\theta_{1})=K_{\mathrm{I}}^{0}(a_{0})\left[1+\sum_{k}L(k)\dfrac{\widehat{\delta
a}_{k}}{a_{0}}e^{ik\theta_{1}}+\sum_{k,k^{\prime}}H(k,k^{\prime})\dfrac{\widehat{\delta
a}_{k}}{a_{0}}\dfrac{\widehat{\delta
a}_{k^{\prime}}}{a_{0}}e^{i(k+k^{\prime})\theta_{1}}\right]+\mathcal{O}(A^{3}),$
(24)

where the first- and second-order convolution kernels $L$ and $H$ read:

$\begin{cases}L(k)=\dfrac{a_{0}}{K_{\mathrm{I}}^{0}(a_{0})}\dfrac{dK_{\mathrm{I}}^{0}}{da}(a_{0})-\dfrac{|{k}|}{2}\\\
H(k,k^{\prime})=\dfrac{a_{0}^{2}}{2K_{\mathrm{I}}^{0}(a_{0})}\dfrac{d^{2}K_{\mathrm{I}}^{0}}{da^{2}}(a_{0})-\dfrac{a_{0}}{K_{\mathrm{I}}^{0}(a_{0})}\dfrac{dK_{\mathrm{I}}^{0}}{da}(a_{0})\dfrac{|{k+k^{\prime}}|}{4}+\dfrac{R(k,k^{\prime})}{16}+\dfrac{|{k}|+|{k^{\prime}}|}{4}.\end{cases}$
(25)

##  3 Validation of the proposed theory

We aim now at validating the second-order expansion of the stress intensity
factor given in Eq. (24). One may consider that retrieving Leblond et al.
(2012)’s and Vasoya et al. (2013)’s results in the limit of very large crack
radius $a_{0}$ proves enough to validate our calculations. Yet, Leblond et al.
(2012)’s second-order formula has been found conflicting with those previously
derived by Adda-Bedia et al. (2006) and Katzav et al. (2007), and its validity
has been checked through _numerical simulations_ only. As a consequence, one
should seek further validation cases to assess the correctness of Eq. (24).
Here, we follow two different routes: first, (i) we consider elliptical cracks
as perturbed circular cracks and test our formula against _analytical results_
obtained for this geometry (Irwin, 1958). It provides a critical test for our
second-order theory, and yields a _first analytical validation_ of Leblond et
al. (2012)’s formula for the semi-infinite crack. Next, (ii) we use Lazarus
(2003)’s model that computes the mode I SIF variations at all orders to
compare our second-order expansions to numerical results obtained for circular
cracks deformed by cosinusoidal perturbations of varying modes and amplitudes.
This ultimately allows us to estimate the validity range of our second-order
theory.

###  3.1 SIF distribution along the front of an elliptical crack

In contrast with models based on the semi-infinite crack geometry (Leblond et
al., 2012; Vasoya et al., 2013), our asymptotic theory derived for penny-
shaped cracks can be tested against analytical solutions obtained for planar
cracks with elliptical front shapes. Here, we consider a crack whose front
describes an ellipsis centered in $O$, of short semi-axis $c$ aligned with
$\mathbf{e_{z}}$, and long semi-axis $b$ aligned with $\mathbf{e_{x}}$. In the
polar coordinate system of Fig. 1, the front position is parametrized by the
function:

$a(\theta)=\dfrac{bc}{\sqrt{b^{2}\sin^{2}(\theta)+c^{2}\cos^{2}(\theta)}}.$
(26)

The crack is loaded by a uniform remote tensile stress $\sigma$. The
distribution of mode I stress intensity factor along the elliptical crack
front is known _exactly_ from the work of Irwin (1958). It reads:

$K_{\mathrm{I}}^{e}(\theta)=\dfrac{\sigma\sqrt{\pi
c}}{\mathrm{E}(\kappa)}\left(\dfrac{c^{4}/b^{4}\cos^{2}(\theta)+\sin^{2}(\theta)}{c^{2}/b^{2}\cos^{2}(\theta)+\sin^{2}(\theta)}\right)^{1/4},$
(27)

where $\mathrm{E}(\kappa)$ is the elliptical integral of the second kind
evaluated at $\kappa=\left(1-c^{2}/b^{2}\right)^{1/2}$. When $b=c=a_{0}$, one
retrieves the case of a circular crack front of radius $a_{0}$:

$K_{\mathrm{I}}^{0}(a_{0})=\dfrac{2}{\pi}\sigma\sqrt{\pi a_{0}},$ (28)

so that the convolution kernels of Eq. (25) read:

$\begin{cases}L(k)=-\dfrac{|{k}|}{2}+\dfrac{1}{2}\\\
\vspace{2mm}H(k,k^{\prime})=\dfrac{R(k,k^{\prime})}{16}+\dfrac{|{k}|+|{k^{\prime}}|}{4}-\dfrac{|{k+k^{\prime}}|}{8}-\dfrac{1}{8}.\end{cases}$
(29)

To derive asymptotic values of $K_{\mathrm{I}}^{e}$ from our second-order
theory, one needs to choose (i) a small parameter $\epsilon$ that quantifies
the deviation from circularity of the elliptical crack front, and (ii) the
radius $a_{0}$ of the reference circular configuration. We follow here the
reasoning of Gao and Rice (1987b), and take $\epsilon=1-c/b$. However, they
considered a reference radius $a_{0}$ equal to the crack front position at the
point of SIF evaluation, i.e. $a_{0}=a(\theta_{1})$. Consequently, the
perturbation $\delta a(\theta)=a(\theta)-a_{0}$ is zero at
$\theta=\theta_{1}$, and the double sum of Eq. (23) containing terms in
$|{k}|+|{k^{\prime}}|$:

$\sum_{k,k^{\prime}}\dfrac{|{k}|+|{k^{\prime}}|}{4}\widehat{\delta
a}_{k}\widehat{\delta
a}_{k^{\prime}}e^{i(k+k^{\prime})\theta_{1}}=\left[\sum_{k}\dfrac{|{k}|}{2}\widehat{\delta
a}_{k}e^{i(k+k^{\prime})\theta_{1}}\right]\times\delta a(\theta_{1})$ (30)

is also null. Yet, these terms are precisely those we need to check first, as
they vanish in the limit of very large crack radii, and cannot be tested
against Vasoya et al. (2013)’s formula for the semi-infinite crack. Here, we
rather take the length $c$ of the short semi-axis as reference radius $a_{0}$.
Following Eq. (26), the perturbation $\delta a=a-a_{0}$ can be expanded as:

$\displaystyle\delta a(\theta)/a_{0}=$ $\displaystyle\epsilon
A_{1}(\theta)+\epsilon^{2}A_{2}(\theta)+\mathcal{O}(\epsilon^{3})$ (31)
$\displaystyle=$
$\displaystyle\epsilon\left[\frac{1}{2}+\frac{1}{2}\cos(2\theta)\right]+\epsilon^{2}\left[\frac{5}{16}+\frac{1}{2}\cos(2\theta)+\frac{3}{16}\cos(4\theta)\right]+\mathcal{O}(\epsilon^{3}),$

where $A_{1}$ and $A_{2}$ are the first- and second-order expansions of
$\delta a/a_{0}$. Similarly, the exact SIF distribution $K_{\mathrm{I}}^{e}$
given in Eq. (27) can be expanded to second-order in $\epsilon$:

$K_{\mathrm{I}}^{e}(\theta)/K_{\mathrm{I}}^{0}\left(a_{0}\right)=1+\epsilon\left[\frac{1}{4}-\frac{1}{4}\cos(2\theta)\right]+\epsilon^{2}\left[\frac{11}{64}-\frac{3}{16}\cos(2\theta)-\frac{11}{64}\cos(4\theta)\right]+\mathcal{O}(\epsilon^{3}).$
(32)

One seeks now to find back Eq. (32) from our non-linear theory. Using Eq. (32)
in Eq. (24), the second-order estimate $K_{\mathrm{I}}^{p}$ of the mode I SIF
reads:

$K_{\mathrm{I}}^{p}(\theta)/K_{\mathrm{I}}^{0}(a_{0})=1+\epsilon\left[\sum_{k}L(k)\hat{A}^{1}_{k}e^{ik\theta}\right]+\epsilon^{2}\left[\sum_{k}L(k)\hat{A}^{2}_{k}e^{ik\theta}+\sum_{k,k^{\prime}}H(k,k^{\prime})\hat{A}^{1}_{k}\hat{A}^{1}_{k^{\prime}}e^{i(k+k^{\prime})\theta}\right]+\mathcal{O}(\epsilon^{3}),$
(33)

where $\hat{A}^{1}_{k}$ and $\hat{A}^{2}_{k}$ are the Fourier coefficients of
$A_{1}$ and $A_{2}$. Using Eqs. (29) and (31), one finds:

$\displaystyle\sum_{k}L(k)\hat{A}^{1}_{k}e^{ik\theta}=\dfrac{1}{4}-\dfrac{1}{4}\cos(2\theta),$
(34)
$\displaystyle\sum_{k}L(k)\hat{A}^{2}_{k}e^{ik\theta}=\dfrac{5}{32}-\dfrac{1}{4}\cos(2\theta)-\dfrac{9}{32}\cos(4\theta),$
and
$\displaystyle\sum_{k,k^{\prime}}H(k,k^{\prime})\hat{A}^{1}_{k}\hat{A}^{1}_{k^{\prime}}e^{i(k+k^{\prime})\theta}=\dfrac{1}{64}+\dfrac{1}{16}\cos(2\theta)+\dfrac{7}{64}\cos(4\theta),$

which yields the result expected from Eq. (32). Note that we also checked that
$K_{\mathrm{I}}^{e}$ and our expansion $K_{\mathrm{I}}^{p}$ are equal up to
second-order when one takes $a_{0}=a(\theta_{1})$ as Gao and Rice (1987b) did.
Further validations of our second-order theory are provided in C, even for the
case of a non-axisymmetric $K_{\mathrm{I}}^{0}$.

The correctness of our Eq. (24) in turn provides a first analytical validation
of Leblond et al. (2012)’s formula, which is obtained in the limit
$a_{0}\rightarrow+\infty$. Our non-linear perturbative approach yields
estimates of the stress intensity factor along elliptical crack fronts for
$b/c\simeq 1$, even when it arises from discontinuous stress distributions
(e.g. a point-force loading). It is complementary to other methods available
in the literature that provide $K_{\mathrm{I}}$ estimates along elliptical
crack fronts for arbitrary $b/c$ ratio, but considering polynomial stress
distributions (Shah and Kobayashi, 1971; Kassir, 1975; Atroshchenko et al.,
2009). However, as we will see next, its scope is much wider as it can be used
on perturbed fronts of arbitrary shapes.

###  3.2 SIF distribution along a sinusoidal crack front

Now that our second-order theory has been validated considering elliptical
crack fronts, one needs to assess its performances in estimating the stress
intensity factor distribution along randomly perturbed crack fronts. As
elliptical crack fronts are associated with second-order perturbations
consisting of the superposition of several deformation modes (see Eq. (31)),
they do not allow us to quantify the errors made by our second-order model for
a given perturbation mode of fixed amplitude. We rather investigate a wavy
crack front described by:

$a(\theta)=a_{0}\left[1+\epsilon\cos(k\theta)\right],$ (35)

where the integer $k\geq 2$ is the mode of the perturbation, and $\epsilon$ is
a small parameter that quantifies its amplitude. We assume next that the crack
parametrized by Eq. (35) is loaded by a uniform remote tensile stress
$\sigma$. As there is no analytical formula available for the mode I stress
intensity factor along the front of such wavy cracks, one needs to compute
$K_{\mathrm{I}}$ numerically. Here, we use the method of Lazarus (2003) and
David and Lazarus (2022) that computes $K_{\mathrm{I}}$ along an arbitrary
crack front using the perturbative approach of Rice (1989) iteratively.
Namely, (i) we start from a penny-shaped configuration for which the mode I
SIF $K_{\mathrm{I}}$ and the fundamental kernel $W$ are known. Then, (ii) we
discretize the crack front in $N$ equidistant points, and sequentially deform
it by a small amount $\delta a(\theta)=\Delta\epsilon\cos(k\theta)$ with
$\Delta\epsilon=0.0025$, and (iii) update both $K_{\mathrm{I}}$ and $W$ using
formulae similar to Eqs. (2) and (5). As $\Delta\epsilon\ll 1$, the first-
order assumption of Rice (1989) applies, and $K_{\mathrm{I}}$ is computed at
_all orders_ in $\epsilon$. The update of the stress intensity factor and the
fundamental kernel are costly procedures, but this method still demands less
computational resources than e.g. simulations based on the finite element
method (FEM) or the boundary element method (BEM), as only the crack front
needs to be meshed. This method is also powerful in modeling far-field loading
applied at infinity or hypotheses of infinite body.

Using our second-order model of Eq. (24), together with the convolution
kernels of Eq. (29), the stress intensity factor variations $\delta
K_{\mathrm{I}}^{p}$ from the uniform $K_{\mathrm{I}}^{0}$ reference reads:

$\delta
K_{\mathrm{I}}^{p}(\theta)/K_{\mathrm{I}}^{0}(a_{0})=-\dfrac{k-1}{2}\epsilon\cos(k\theta)-\dfrac{k^{2}-4k+1}{16}\epsilon^{2}+\dfrac{k^{2}+2k-1}{16}\epsilon^{2}\cos(2k\theta)+\mathcal{O}(\epsilon^{3}),$
(36)

where $K_{\mathrm{I}}^{0}(a_{0})$ is defined in Eq. (28).

Figure 3: Normalized stress intensity factor variations $\delta
K_{\mathrm{I}}/\epsilon K_{\mathrm{I}}^{0}$ along a crack front perturbed by
$\delta a(\theta)/a_{0}=\epsilon\cos(k\theta)$ for (a) $k=2$ and (b) $k=5$,
and $\epsilon=0.01$ (left panel), $\epsilon=0.1$ (middle panel), and
$\epsilon=0.25$ (right panel). The SIF variations
$K_{\mathrm{I}}^{\mathrm{num}}$ computed from Lazarus (2003)’s numerical model
(in black solid line) are compared to the first-order estimate of Gao and Rice
(1987b) (in dotted orange line) and our second-order prediction (in dashed red
line). The errors on the SIF variations between the numerical simulations (in
solid line) and the first-order (in dotted line) and second-order (in dashed
line) models are compared using (c) $L^{0}$ (maximum) or (d) $L^{2}$
(Euclidean) norms. They increase with both the perturbation mode $k$ and
amplitude $\epsilon$. Deviations of the numerical results from the first- and
second-order predictions can be directly observed on (e) the average
$\left<K_{\mathrm{I}}\right>$ of the stress-intensity factor along the crack
front, or (f) its amplitude
$K_{\mathrm{I}}^{\mathrm{max}}-K_{\mathrm{I}}^{\mathrm{min}}$, for which first
and second-order estimates coincide.

The numerical results using the method of Lazarus (2003) (in black solid line)
are compared to the second-order prediction of Eq. (36) (in dashed red line)
and to its first-order truncation (in dotted orange line) in Fig. 3. For $k=2$
and $\epsilon=0.01$ (see Fig. 3a), the three quantities are close. When
$\epsilon=0.1$, the first-order assumption is no longer valid, and the linear
estimate departs from the numerical results, while the second-order prediction
remains accurate. As $\epsilon$ grows even larger ($\epsilon=0.25$), the
second-order prediction starts overestimating the stress intensity factor on
the regions where the crack is less elongated. Note however that the second-
order model do predict an upward shift in the average of the SIF distribution,
while the first-order prediction remains symmetric with respect to $\delta
K_{\mathrm{I}}=0$. For $k=5$ (see Fig. 3b), the validity range of our theory
shrinks, and departure of the second-order prediction from the numerical
results is already noticeable when $\epsilon=0.1$. When $\epsilon=0.25$, the
discrepancies are fairly large, and our model fails in predicting the strong
decrease of SIF amplitude at $\theta=\pi/5$. We observe that the errors made
by our second-order model around this point are even larger than those made by
the first-order theory. These preliminary observations are further
investigated in Fig. 3c-d where we show the evolution with $\epsilon$ of the
(maximal) error using the $L^{0}$ norm:

$||\,\mathrm{error}\,||_{0}=\max_{\theta\in\left[0,2\pi\right]}|{\delta
K_{\mathrm{I}}^{p}(\theta)-\delta
K_{\mathrm{I}}^{\mathrm{num}}(\theta)}|/K_{\mathrm{I}}^{0}(a_{0}),$

and the (average) error using the $L^{2}$ norm:

$||\,\mathrm{error}\,||_{2}=\dfrac{1}{2\pi/K_{\mathrm{I}}^{0}(a_{0})}\int_{0}^{2\pi}\left(\delta
K_{\mathrm{I}}^{p}(\theta)-\delta
K_{\mathrm{I}}^{\mathrm{num}}(\theta)\right)^{2}d\theta.$

Here, $\delta K_{\mathrm{I}}^{\mathrm{num}}$ is the mode I SIF variation
predicted using Lazarus (2003)’s method. Values of $\epsilon$ for which the
errors made by the first- and second-order models exceeds a threshold value
are listed in Table 1. Note that we report errors made by our second-order
prediction against a numerical estimate of $K_{\mathrm{I}}$, since its exact
value is not known. Nonetheless, for the discretization adopted here
($N=\left\\{252,378,680,850\right\\}$ for $k=\left\\{2,3,4,5\right\\}$, and
$\delta a=0.0025$), the numerical error made in evaluating $\delta
K_{\mathrm{I}}^{\mathrm{num}}$ can be estimated to be less than $0.2\%$. This
estimation builds on the error analysis made by Lazarus (2003) on
$K_{\mathrm{I}}$ evaluation along elliptical crack fronts with $b/c=3$ (i.e.
$\epsilon=0.66$), that necessitates more small perturbations iterations than
here where $\epsilon\leq 0.25$.

$||\cdot||_{0}$ $||\cdot||_{2}$ 1% 5% 20% 1% 5% 20% $k=2$ 0.2275 0.3375 0.48
0.3725 $>$0.5 $>$0.5 $k=3$ 0.1425 0.2225 0.34 0.2425 0.3825 $>$0.5 $k=4$ 0.105
0.1675 0.265 0.18 0.295 0.4875 $k=5$ 0.085 0.135 0.215 0.1425 0.2375 0.405

Table 1: Values of perturbation amplitude $\epsilon$ ($\pm 0.0025$) above
which the errors in the SIF variation $\delta K_{\mathrm{I}}$ made by the
second-order theory exceeds a certain threshold value
$\in\left\\{1\%,5\%,25\%\right\\}$. Values of $\epsilon$ above $0.5$ have not
been investigated in this study.

We observe nonetheless that our second-order formula generally yields more
accurate predictions than the linear theory of Gao and Rice (1987b), and that
quantitative agreement between the numerical results and the second-order
expansion of Eq. (36) is found for $\epsilon\lesssim 0.1$ and
$k\in\left[0,5\right]$. This is further confirmed by the evolution of the
average SIF $\left<K_{\mathrm{I}}\right>$ with $\epsilon$ (see Fig. 3e):

$\left<K_{\mathrm{I}}\right>=\dfrac{1}{2\pi}\int_{0}^{2\pi}K_{\mathrm{I}}(\theta)d\theta=K_{\mathrm{I}}^{0}(a_{0})\times\left[1-\dfrac{k^{2}-4k+1}{16}\epsilon^{2}\right]+\mathcal{O}(\epsilon^{3}),$
(37)

and the SIF amplitude
$K_{\mathrm{I}}^{\mathrm{max}}-K_{\mathrm{I}}^{\mathrm{min}}$ (see Fig. 3f):

$K_{\mathrm{I}}^{\mathrm{max}}-K_{\mathrm{I}}^{\mathrm{min}}=\max_{\theta\in\left[0,2\pi\right]}K_{\mathrm{I}}(\theta)-\min_{\theta\in\left[0,2\pi\right]}K_{\mathrm{I}}(\theta)=K_{\mathrm{I}}^{0}(a_{0})\times(k-1)\epsilon+\mathcal{O}(\epsilon^{3}).$
(38)

The change in average SIF $\left<K_{\mathrm{I}}\right>$ appears indeed of
second order in $\epsilon$, and is thus overlooked by the first-order theory.
Moreover, we observe a change in the sign of the second-order terms in
$\left<K_{\mathrm{I}}\right>$ between $k=3$ and $4$. On the contrary, there is
no second-order contribution to the SIF amplitude
$K_{\mathrm{I}}^{\mathrm{max}}-K_{\mathrm{I}}^{\mathrm{min}}$, and higher-
order theories would be required to predict the strong decrease of
$K_{\mathrm{I}}$ near $\theta=\pi/k$ observed in Fig. 3b. We observed that the
abrupt change of slope in the SIF amplitude computed using Lazarus (2003)’s
method is associated with a change in the point where $K_{\mathrm{I}}$ is
minimal, from $\theta=\pi/4k$ (see the third panel of Fig. 3b) to
$\theta=\pi/k$.

We conclude that (i) our second-order theory generally yields better estimates
of the SIF variations along perturbed crack fronts than the first-order model
of Gao and Rice (1987b), but (ii) its accuracy decreases with increasing
perturbation amplitude $\epsilon$ and mode $k$. As we will see next, the
latter may be rather inconsequential, as the amplitude of the front
perturbations is usually a fast decreasing function ($\propto 1/k$) of their
mode $k$.

##  4 Coplanar crack propagation in weakly disordered materials

The influence of microscopic fluctuations of toughness on the critical load at
which a material may fail has been studied quite extensively for the semi-
infinite crack geometry since the seminal study of Roux et al. (2003). Various
methods have been used in the literature to address this problem, like e.g.
phase-field models (Hossain et al., 2014; Brach et al., 2019) or minimal
surfaces (Ernesti and Schneider, 2021; Michel and Suquet, 2022). But
perturbative approaches have played an important role in homogenizing fracture
properties. Efforts have been made in both proving the existence of a finite
effective toughness for random media (Dirr et al., 2011; Dondl and Jesenko,
2020), and predicting its value from the statistical (Roux and Hild, 2008;
Patinet et al., 2013b; Démery et al., 2014b) and spatial (Démery et al.,
2014a; Lebihain, 2021) distribution of toughness $K_{\mathrm{Ic}}(z,x)$. In
particular, it has been shown that the material is toughened by the presence
of disorder, i.e. the effective toughness $K_{\mathrm{Ic}}^{\mathrm{eff}}$ is
always larger than the average value $K_{\mathrm{Ic}}^{0}$ of the local
toughness field. A major pitfall of these studies is that they predict a
toughening that scales as the _second-order_ of the disorder intensity
(measured by e.g. the standard deviation $\sigma$ of the toughness field) from
a _first-order_ expansion of the stress-intensity factor. Hence, they neglect
potential contributions arising from the non-linearities in $K_{\mathrm{I}}$
and may be inaccurate in predicting $K_{\mathrm{Ic}}^{\mathrm{eff}}$.
Moreover, the semi-infinite crack geometry considered in these works also
prevents investigating size-effects emerging from the finite length of the
crack.

Here, we investigate crack propagation when $K_{\mathrm{Ic}}$ is fluctuating
randomly within the failure plane $y=0$, and we focus on the effective
toughness $K_{\mathrm{Ic}}^{\mathrm{eff}}$ of such composite materials. We
propose a rigorous homogenization framework for weakly disordered materials
that correctly predicts the non-linear toughening arising from local
fluctuations of toughness from a second-order expansion of the stress
intensity factor. We also examine the influence of the finiteness of the crack
geometry on material toughening, as previous studies mainly focused on the
semi-infinite crack geometry. Our analytical predictions are finally compared
to the results of $\sim 1$ millions of efficient numerical simulations of
crack propagation in heterogeneous brittle media.

###  4.1 Problem statement and field generation

We consider a penny-shape crack of radius $a_{\mathrm{ini}}$ centered in $O$.
The crack propagates along the failure plane located at $y=0$, and interacts
with the non-uniform toughness field:

$K_{\mathrm{Ic}}(r,\theta)=K_{\mathrm{Ic}}^{0}+\sigma f(r,\theta),$ (39)

where $K_{\mathrm{Ic}}^{0}$ corresponds to the spatial average of the
toughness field, $\sigma$ is its standard deviation, and $f$ is a fluctuation
term of zero average and unit variance. Fluctuations are either described by
the function $f$ in polar coordinates or by $f^{*}$ in Cartesian coordinates.
The two functions are linked through the equation:

$f(r,\theta)=f^{*}(z,x)\text{, where }z=r\cos\theta\text{ and }x=r\sin\theta.$
(40)

An example of toughness field is given in Fig. 4a. It is characterized its
probability density function $\mathcal{P}$ that measures the statistical
distribution of the local toughness. Here, we consider that the values of
$K_{\mathrm{Ic}}$ are uniformly distributed between two extremal values
$K_{\mathrm{Ic}}^{\mathrm{min}}$ and $K_{\mathrm{Ic}}^{\mathrm{max}}$ (see
Fig. 4b), so that:

$\mathcal{P}(K_{\mathrm{Ic}})=\dfrac{1}{K_{\mathrm{Ic}}^{\mathrm{max}}-K_{\mathrm{Ic}}^{\mathrm{min}}}.$
(41)

Consequently, its mean value $K_{\mathrm{Ic}}^{0}$ reads:

$K_{\mathrm{Ic}}^{0}=\dfrac{K_{\mathrm{Ic}}^{\mathrm{max}}+K_{\mathrm{Ic}}^{\mathrm{min}}}{2},$
(42)

and the intensity $\bar{\sigma}=\sigma/K_{\mathrm{Ic}}^{0}$ of the material
disorder is equal to:

$\bar{\sigma}=\dfrac{1}{\sqrt{3}}\dfrac{K_{\mathrm{Ic}}^{\mathrm{max}}-K_{\mathrm{Ic}}^{\mathrm{min}}}{K_{\mathrm{Ic}}^{\mathrm{max}}+K_{\mathrm{Ic}}^{\mathrm{min}}}.$
(43)

The toughness field is further characterized by its spatial structure,
described e.g. by the two-points correlation function of $f$. In the
following, $\mathrm{E}\left[\cdot\right]$ denotes the expectation of a random
variable over multiple disorder realizations. We consider (i) _statistically
homogeneous_ and (ii) _ergodic_ materials, so that ensemble averages (i) do
not depend on the observation point $M$, and (ii) are equivalent to spatial
average over large enough surfaces. Using these notations, the two-points
correlation function of the fluctuations reads:

$\mathrm{E}\left[{f^{*}(z,x)f^{*}(z^{\prime},x^{\prime})}\right]=\mathcal{F}_{z}(|{z-z^{\prime}}|)\mathcal{F}_{x}(|{x-x^{\prime}}|),$
(44)

where $\mathcal{F}_{z}$ and $\mathcal{F}_{x}$ are the correlation functions in
the $z$ and $x$ direction, respectively. Here we assume that the material is
isotropic and that the correlations follows a bell-curve, so that:

$\mathcal{F}_{z}(u)=\mathcal{F}_{x}(u)=\mathcal{F}(u)=e^{-u^{2}/d^{2}},$ (45)

where $d$ is the characteristic size of the disorder (see Fig. 4c).

Figure 4: (a) Front positions (in solid lines) of a crack propagating in a
disordered toughness field $K_{\mathrm{Ic}}$ under the action of a point force
$P$ applied at the crack center. The fronts are plotted for a constant load
increment. (b) The values of $K_{\mathrm{Ic}}$ are uniformly distributed
around a mean value $K_{\mathrm{Ic}}^{0}$ with a probability density function
$\mathcal{P}(K_{\mathrm{Ic}})=1/(K_{\mathrm{Ic}}^{\mathrm{max}}-K_{\mathrm{Ic}}^{\mathrm{min}})$,
where $K_{\mathrm{Ic}}^{\mathrm{max}}$ and $K_{\mathrm{Ic}}^{\mathrm{min}}$
are the maximum and minimum value of the distribution. (c) The field
$K_{\mathrm{Ic}}$ is spatially correlated at a typical heterogeneity scale $d$
correlations $\mathcal{F}(\Delta r/d)$ (in solid orange line) following
$e^{-(\Delta r/d)^{2}}$ (in dashed black line). (d) The load $P$ at which the
crack propagates (in solid red line), normalized by
$P^{0}_{d}=K_{\mathrm{Ic}}^{0}(\pi d)^{3/2}$, fluctuates around its reference
value obtained when the toughness field is uniform to $K_{\mathrm{Ic}}^{0}$.
(e) As a result, one may estimate an effective toughness
$K_{\mathrm{Ic}}^{\mathrm{eff}}$ from the critical loading required to
propagate a crack of average radius $\left<a\right>$.

Realizations of $K_{\mathrm{Ic}}$ fields are produced following the procedure
of Albertini et al. (2021), suitably extended to two-dimensional surfaces.
First, a Gaussian random field of zero average, unit variance, and controlled
power-spectrum is constructed using the Python package FyeldGenerator (Cadiou,
2022). We then apply a non-linear mapping from the Gaussian cumulative
distribution function to that of the uniform distribution. This transformation
is expected not to affect much the correlation shape $\mathcal{F}$ and the
correlation length $d$, as $\mathcal{F}$ is positive (Grigoriu, 2002) (see
Fig. 4c).

###  4.2 Computation of the crack front position

The crack is loaded by a pair of symmetric tensile forces of magnitude $P$
applied at the crack center $O$. It gives rise to a stress intensity
$K_{\mathrm{I}}^{0}$ along the front of a penny-shaped crack of radius
$a_{0}$, which reads (Tada et al., 2000):

$K_{\mathrm{I}}^{0}(a_{0})=\dfrac{P}{(\pi a_{0})^{3/2}}.$ (46)

The magnitude $P$ of the tensile forces slowly increases in time, so that the
crack extends quasi-statically due to the stabilizing nature of the loading.
The successive crack front positions are computed from Irwin’s criterion:

$K_{\mathrm{I}}(\theta)=K_{\mathrm{Ic}}(a(\theta),\theta).$ (47)

$K_{\mathrm{I}}(\theta)$ is usually computed from the front deformation
$\delta a$ from its reference circular configuration of radius $a_{0}$ using
the first-order expansion of Rice (1985), while
$K_{\mathrm{Ic}}(a(\theta),\theta)$ is computed exactly i.e. at all orders in
the front deformation. The non-linear Eq. (47) is then solved by building on a
viscous regularization of the evolution equation (Patinet et al., 2013b;
Lebihain, 2021), or by resorting to a comparison theorem to compute
efficiently the successive state equilibria (Rosso and Krauth, 2002; Démery et
al., 2014a). Here, we follow a different approach, and expand the front
position to second-order in the disorder intensity
$\bar{\sigma}=\sigma/K_{\mathrm{Ic}}^{0}$:

$a(\theta)=a_{0}\left[1+\bar{\sigma}A_{1}(a_{0},\theta)+\bar{\sigma}^{2}A_{2}(a_{0},\theta)\right]+\mathcal{O}(\bar{\sigma}^{3}),$
(48)

where $a_{0}$ is the radius of a crack propagating on a homogeneous medium of
toughness $K_{\mathrm{Ic}}^{0}$ under the action of the loading $P$. $A_{1}$
and $A_{2}$ are fluctuation terms that emerge from the local variations $f$ of
the toughness field. One can then expand the local SIF
$K_{\mathrm{I}}(\theta)$ at second-order using our Eq. (16):

$K_{\mathrm{I}}(\theta)=K_{\mathrm{I}}^{0}(a_{0})+\bar{\sigma}K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[a_{0}A_{1}];\theta)+\bar{\sigma}^{2}\left(K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[a_{0}A_{2}];\theta)+K_{\mathrm{I}}^{2}(a_{0},[K_{\mathrm{I}}^{0}],[a_{0}A_{1}];\theta)\right)+\mathcal{O}(\bar{\sigma}^{3}),$
(49)

where $K_{\mathrm{I}}^{1}$ and $K_{\mathrm{I}}^{2}$ are given by Eq. (24),
with the following convolution kernels:

$\begin{cases}\begin{aligned} &L(k)=-\dfrac{|{k}|}{2}-\dfrac{3}{2}\\\
\vspace{2mm}&H(k,k^{\prime})=\dfrac{R(k,k^{\prime})}{16}+\dfrac{|{k}|+|{k^{\prime}}|}{4}+\dfrac{3|{k+k^{\prime}}|}{8}+\dfrac{15}{8}.\end{aligned}\end{cases}$
(50)

Similarly, the local toughness $K_{\mathrm{Ic}}(a(\theta),\theta)$ can be
expanded at second-order following:

$K_{\mathrm{Ic}}(a(\theta),\theta)=K_{\mathrm{Ic}}^{0}\left(1+\bar{\sigma}f(a_{0},\theta)+\bar{\sigma}^{2}\dfrac{\partial
f}{\partial r}(a_{0},\theta)a_{0}A_{1}(a_{0},\theta)\right).$ (51)

The crack front shape is then computed by solving Irwin’s criterion of Eq.
(47) in cascade, i.e. at successive but increasing orders in $\bar{\sigma}$.
Starting with zero order, one gets from Eq. (46):

$K_{\mathrm{I}}^{0}(a_{0})=K_{\mathrm{Ic}}^{0}\Rightarrow
a_{0}=\dfrac{1}{\pi}\left(\dfrac{P}{K_{\mathrm{Ic}}^{0}}\right)^{2/3},$ (52)

which sets the instantaneous value of the loading $P$. One could have
alternatively set the magnitude $P$ of the opening forces, and solve Eq. (52)
to find the reference radius $a_{0}$. At first order, one has:

$K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[a_{0}A_{1}];\theta)=K_{\mathrm{Ic}}^{0}f(a_{0},\theta)\Rightarrow\sum_{k}L(k)\hat{A}^{1}_{k}(a_{0})e^{ik\theta}=\sum_{k}\hat{f}_{k}(a_{0})e^{ik\theta},$
(53)

so that the $k$-th Fourier coefficient of $A_{1}$ reads:

$\hat{A}^{1}_{k}(a_{0})=\dfrac{\hat{f}_{k}(a_{0})}{L(k)},$ (54)

where $L(k)$ is given in Eq. (50). We observe here that the perturbation
amplitude is a fast decreasing function of its mode $k$ ($\propto 1/k$), and
that crack fronts are generally stiffer to small wavelength perturbations.
Collecting the second-order terms in Eq. (47), one finds:

$\displaystyle
K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[a_{0}A_{2}];\theta)+K_{\mathrm{I}}^{2}(a_{0},[K_{\mathrm{I}}^{0}],[a_{0}A_{1}];\theta)=K_{\mathrm{Ic}}^{0}a_{0}\dfrac{\partial
f}{\partial r}(a_{0},\theta)A_{1}(\theta)$ (55) $\displaystyle\Rightarrow$
$\displaystyle\sum_{k}L(k)\hat{A}^{2}_{k}(a_{0})e^{ik\theta}=\sum_{k,k^{\prime}}a_{0}\widehat{\left(\dfrac{\partial
f}{\partial
r}\right)}_{k^{\prime}}\,(a_{0})\hat{A}^{1}_{k-k^{\prime}}(a_{0})e^{ik\theta}-\sum_{k,k^{\prime}}H(k^{\prime},k-k^{\prime})\hat{A}^{1}_{k^{\prime}}(a_{0})\hat{A}^{1}_{k-k^{\prime}}(a_{0})e^{ik\theta},$

where $H(k,k^{\prime})$ is given in Eq. (50). We observe that each mode
$\hat{A}^{2}_{k}$ of the second-order fluctuations $A_{2}$ are linked to all
perturbation modes of $A_{1}$. Injecting Eq. (54), the $k$-th Fourier
coefficient of $A_{2}$ reads:

$\hat{A}^{2}_{k}(a_{0})=\dfrac{1}{L(k)}\sum_{k^{\prime}}\left(\dfrac{a_{0}}{L(k-k^{\prime})}\widehat{\left(\dfrac{\partial
f}{\partial
r}\right)}_{k^{\prime}}\,(a_{0})\hat{f}_{k-k^{\prime}}(a_{0})-\dfrac{H(k^{\prime},k-k^{\prime})}{L(k^{\prime})L(k-k^{\prime})}\hat{f}_{k^{\prime}}(a_{0})\hat{f}_{k-k^{\prime}}(a_{0})\right).$
(56)

In practice, the crack front shape $a(\theta)$ and the critical loading $P$
are computed using Eqs. (48), (52), (54), and (56). For a given reference
crack radius $a_{0}$, the magnitude of the applied force $P$ is computed from
Eq. (52). The crack front is then discretized in $n=1024$ intervals of equal
length, and Eqs. (54), and (56) are solved using fast Fourier Transforms. The
computation time for one reference front position $a_{0}$ scales as
$\mathcal{O}(n^{2}\ln n)$. Examples of crack front shapes and loading
evolution are given in Fig. 4a and 4d.

Our numerical approach differs from those available in the literature, as only
second-order contributions of the toughness fluctuations are taken into
account in coherence with the second-order expansion of $K_{\mathrm{I}}$. We
expect it to break as soon as the front deformation $\delta
a(\theta)=a(\theta)-a_{0}$ gets larger than the characteristic length scale
$d$ of the toughness variations. Indeed, for larger front distortions, the
second-order expansion of the local toughness in Eq. (51) may not be valid
anymore. This strongly limits (i) the intensity $\bar{\sigma}$ of the material
disorder, and (ii) the ratio $a_{0}/d$ between the crack radius to
heterogeneity size considered in the simulations. Namely, we measured from our
numerical computations that we could not model crack radii larger than
$a_{0}=20d$ for $\bar{\sigma}=0.1$.

###  4.3 Effective toughness of weakly disordered materials

Due to the strong decrease of $K_{\mathrm{I}}^{0}$ with crack advance (see Eq.
(46)), we do not observe any snapback instability in the evolution of $P$ with
the average crack radius $\left<a\right>$ (see Fig. 4d), which reads:

$\left<a\right>/a_{0}=1+\bar{\sigma}\hat{A}^{1}_{0}(a_{0})+\bar{\sigma}^{2}\hat{A}^{2}_{0}(a_{0})+\mathcal{O}(\bar{\sigma}^{3}).$
(57)

As a result, all states computed from Eqs. (52), (54), and (56) are state
equilibria. This suggests that the instability nuclei, often referred to as
Larkin domains, are always larger than the crack perimeter $\mathcal{P}\simeq
2\pi\left<a\right>$. Accordingly, our simulations should fall down in the
_weak pinning_ regime identified by Roux et al. (2003) for small disorder
intensity $\bar{\sigma}$. In this regime, the material toughening:

$T=K_{\mathrm{Ic}}^{\mathrm{eff}}/K_{\mathrm{Ic}}^{0}-1$ (58)

is expected to be zero (Roux et al., 2003). Our goals here are (i) to show
that this is wrong if one considers a second-order expansion of both the
stress-intensity factor and the toughness field, and (ii) to assess the
dependence of the observed toughening on the crack size $a_{0}$.

At the macro-scale, one may overlook the undulations of the crack front and
consider the perturbed crack as a circular crack of radius $\left<a\right>$
(see Fig. 4a), which departs from its reference value $a_{0}$ under the action
of the material disorder (see Eq. (57)). As the crack propagates under the
action of the tensile forces, one may then define an apparent, or _effective_
, toughness of the heterogeneous media from the magnitude $P$ of the tensile
forces required to make a circular crack of radius $\left<a\right>$ propagate:

$K_{\mathrm{Ic}}^{\mathrm{eff}}=\dfrac{P}{(\pi\left<a\right>)^{3/2}},$ (59)

so that:

$K_{\mathrm{Ic}}^{\mathrm{eff}}=K_{\mathrm{I}}^{0}(a_{0})\left(\dfrac{a_{0}}{\left<a\right>}\right)^{3/2}=K_{\mathrm{Ic}}^{0}\left(\dfrac{a_{0}}{\left<a\right>}\right)^{3/2}.$
(60)

We observe that the toughening at the macroscale emerges from the difference
between the average front position $\left<a\right>$ and its reference value
$a_{0}$. If the front is pinned by the disorder ($\left<a\right><a_{0}$), the
composite appears tougher at the macroscale than the reference homogeneous
material of toughness $K_{\mathrm{Ic}}^{0}$. On the contrary, if the disorder
facilitates crack propagation ($\left<a\right>>a_{0}$), the material is
weakened by the heterogeneities. Using Eq. (57), one finds:

$K_{\mathrm{Ic}}^{\mathrm{eff}}/K_{\mathrm{Ic}}^{0}=K_{\mathrm{I}}^{0}(\left<a\right>)/K_{\mathrm{I}}^{0}(a_{0})=1+\bar{\sigma}L(0)\hat{A}^{1}_{0}(a_{0})+\bar{\sigma}^{2}\left(L(0)\hat{A}^{2}_{0}(a_{0})+H(0,0)\hat{A}^{1}_{0}(a_{0})^{2}\right),$
(61)

where $L(0)=a_{0}(\partial K_{\mathrm{I}}^{0}/\partial
a)(a_{0})/K_{\mathrm{I}}^{0}(a_{0})$ and
$H(0,0)=a_{0}^{2}(\partial^{2}K_{\mathrm{I}}^{0}/\partial
a^{2})(a_{0})/2K_{\mathrm{I}}^{0}(a_{0})$. An example of the evolution of
$K_{\mathrm{Ic}}^{\mathrm{eff}}$ with the average crack radius
$\left<a\right>$ is given in Fig. 4e.

We observe in Fig. 4d that the effective toughness is a random variable that
seems to fluctuate around $K_{\mathrm{Ic}}^{0}$ and between the two extremal
values $K_{\mathrm{Ic}}^{\mathrm{min}}$ and $K_{\mathrm{Ic}}^{\mathrm{max}}$
of the toughness distribution. To have a better perspective on the statistical
distribution of the random variable $K_{\mathrm{Ic}}^{\mathrm{eff}}$, we
perform more than $1$ million Monte-Carlo simulations of crack propagation in
heterogeneous media, following the procedure described in Section 4.2. Namely,
the disorder intensity is varied between 11 different values
$\bar{\sigma}\in\left[0.001,0.01\right]$. $N=100^{\prime}000$ independent and
identically distributed realizations of the toughness field
$K_{\mathrm{Ic}}(z,x)$ are generated for each value of $\bar{\sigma}$. For
each disorder realization, the crack front position is computed from 16
reference radius $a_{0}\in\left[0.02\,d,20\,d\right]$. It results on $17.6$
millions of crack fronts from which one can estimate the average front
position $\left<a\right>$ using Eq. (57), and the effective toughness
$K_{\mathrm{Ic}}^{\mathrm{eff}}$ following Eq. (61).

Figure 5: (a) At constant disorder intensity
$\sigma/K_{\mathrm{Ic}}^{0}=10^{-1}$, the simulated values of effective
toughness $K_{\mathrm{Ic}}^{\mathrm{eff}}$ are distributed between the two
extremal values $K_{\mathrm{Ic}}^{\mathrm{min}}$ and
$K_{\mathrm{Ic}}^{\mathrm{max}}$ of the local toughness distribution, and its
distribution depends on the ratio between the reference crack radius $a_{0}$
and the heterogeneity size $d$. (b) Similarly, at constant size ratio
$a_{0}/d=0.1$, the distribution of $K_{\mathrm{Ic}}^{\mathrm{eff}}$ depends on
the disorder intensity $\sigma/K_{\mathrm{Ic}}^{0}$ (the distribution plotted
in bold is the same to that shown in 6a). As a result, (c) the average
$\mathrm{E}\left[T\right]$ and (d) the standard deviation
$\mathrm{std}\left[T\right]$ of the toughening
$T=K_{\mathrm{Ic}}^{\mathrm{eff}}/K_{\mathrm{Ic}}^{0}-1$ depends on both
$\sigma/K_{\mathrm{Ic}}^{0}$ and $a_{0}/d$. In particular,
$\mathrm{E}\left[T\right]$ scales as $(\sigma/K_{\mathrm{Ic}}^{0})^{2}$ and
its sign depends on $a_{0}/d$, while $\mathrm{std}\left(T\right)$ is linear in
$\sigma/K_{\mathrm{Ic}}^{0}$.

Figure 6: (a) When normalized by the square of the disorder intensity
$(\sigma/K_{\mathrm{Ic}}^{0})^{2}$, the values of the average
$\mathrm{E}\left[T\right]$ of the toughening showed in Figure 6c all fall down
onto a master curve predicted by Eq. (62). (b) Similarly, when normalized by
$(\sigma/K_{\mathrm{Ic}}^{0})$, the values of the standard deviation
$\mathrm{std}\left[T\right]$ presented in Figure 6d follow the master curve of
Eq. (63) controlled by $a_{0}/d$.

The results of our statistical analysis are summarized in Fig. 6. For
$\bar{\sigma}=0.1$, we observe that the effective toughness is indeed bounded
by the two extremal values $K_{\mathrm{Ic}}^{\mathrm{min}}$ and
$K_{\mathrm{Ic}}^{\mathrm{max}}$ of the toughness distribution (see Fig. 6a),
as expected from Irwin’s criterion. For a small crack radius $a_{0}=0.1\,d$,
the effective toughness follows the probability distribution of the local
toughness $K_{\mathrm{Ic}}$ (see Fig. 4b), as the crack front is often
embedded into one single heterogeneity. As the crack grows larger, the
distribution of the effective toughness departs from that of the local
toughness and strongly resemble a Gaussian distribution. Interestingly, the
variance of the distribution looks to decrease with increasing crack to
heterogeneity size ratio $a_{0}/d$. One can also notice that the average value
of the toughness distribution looks slightly below $K_{\mathrm{Ic}}^{0}$ for
$a_{0}/d=1$, and closer to it for $a_{0}/d=10$. Departure of this average
value from $K_{\mathrm{Ic}}^{0}$ seems in general very small with respect to
the fluctuations of $K_{\mathrm{Ic}}^{\mathrm{eff}}$, so that both quantities
may not be of the same order in $\bar{\sigma}$. We also observe in Fig. 6b
that the fluctuations of the effective toughness increase with the disorder
intensity $\bar{\sigma}$.

In the following, we do not try to predict the exact probability distribution
of $K_{\mathrm{Ic}}^{\mathrm{eff}}$ in Fig. 6a-b, but rather focus on its
statistical descriptors, namely its average value and its standard deviation.
As the values of $K_{\mathrm{Ic}}^{\mathrm{eff}}$ seems to gather around
$K_{\mathrm{Ic}}^{0}$, we concentrate instead on the statistics of the
toughening $T=K_{\mathrm{Ic}}^{\mathrm{eff}}/K_{\mathrm{Ic}}^{0}-1$. By
performing ensemble averages on Eq. (61), one can show that its average reads:

$\mathrm{E}\left[{T}\right]=\left(\dfrac{\sigma}{K_{\mathrm{Ic}}^{0}}\right)^{2}\left[\sum_{k}\left(\dfrac{1}{L(k)}\hat{G}_{k}(a_{0}/d)-\dfrac{H(k,-k)}{L(k)^{2}}\hat{F}_{k}(a_{0}/d)\right)+\dfrac{H(0,0)}{L(0)^{2}}\hat{F}_{0}(a_{0}/d)\right]+\mathcal{O}\left(\left(\dfrac{\sigma}{K_{\mathrm{Ic}}^{0}}\right)^{3}\right),$
(62)

while its standard deviation writes as:

$\mathrm{std}\left[{T}\right]=\left(\dfrac{\sigma}{K_{\mathrm{Ic}}^{0}}\right)\hat{F}_{0}(a_{0}/d)^{1/2}+\mathcal{O}\left(\left(\dfrac{\sigma}{K_{\mathrm{Ic}}^{0}}\right)^{2}\right),$
(63)

where $F$ and $G$ are the functions:

$F(a_{0}/d,\theta)=e^{-4(a_{0}/d)^{2}\sin^{2}(\theta/2)}\text{, and
}G(a_{0}/d,\theta)=-4(a_{0}/d)^{2}\sin^{2}(\theta/2)e^{-4(a_{0}/d)^{2}\sin^{2}(\theta/2)},$
(64)

which both relate to the correlation shape $\mathcal{F}$ of the fluctuations
$f$. Details on the derivation of Eqs. (62) and (63) are given in D. Several
comments are in order:

  * 1.

We observe that the average toughening $\mathrm{E}\left[{T}\right]$ scales as
$\bar{\sigma}^{2}$, while $\mathrm{std}\left[{T}\right]$ increases linearly
with $\sigma$. As a consequence, it is extremely difficult to measure
$\mathrm{E}\left[{T}\right]$ from our “naive” Monte-Carlo simulations, as the
central limit theorem states that the fluctuations of
$\mathrm{E}\left[{T}\right]$ are proportional to $\bar{\sigma}/\sqrt{N}$,
where $N$ is the number of realizations. To circumvent this issue, we use a
variance reduction technique, and resort to antithetic variates. Let us call
$\Omega_{1}$ the ensemble containing all $N=100^{\prime}000$ realizations
generated here for fixed values of $\bar{\sigma}$ and $a_{0}/d$. We attribute
to each realization an integer $i$, and note $T_{1}^{i}$ the value of $T$
measured at $a_{0}$ for a crack propagating in the toughness field
$K_{\mathrm{Ic}}^{i}=K_{\mathrm{Ic}}^{0}(1+\bar{\sigma}f_{i})$. We now
consider the ensemble $\Omega_{2}$ containing the 100’000 realizations
corresponding to the toughness field
$\bar{K_{\mathrm{Ic}}}^{i}=K_{\mathrm{Ic}}^{0}(1-\bar{\sigma}f_{i})$ (opposite
fluctuations), and note $T_{2}^{i}$ the value of $T$ measured at $a_{0}$ in
this case. $T_{2}^{i}$ can be computed at almost no cost, as it corresponds to
a change of sign in front of $A_{1}$ only. As both $f_{i}$ and $-f_{i}$ are
identically distributed (zero average, unit variance, uniform probability
distribution, same two-points correlation function), one has:

$\mathrm{E}\left[{T}\right]=\dfrac{\mathrm{E}\left[{T_{1}}\right]+\mathrm{E}\left[{T_{2}}\right]}{2}\simeq\dfrac{\sum_{i\in\Omega_{2}}T_{1}^{i}+\sum_{i\in\Omega_{2}}T_{2}^{i}}{2N},$
(65)

but one expects that the covariance of the random variables $T_{1}$ and
$T_{2}$ is negative, so that the variance of $\tilde{T}=(T_{1}+T_{2})/2$ is
greatly reduced. Using antithetic variates, one observe from Fig. 6c that the
average toughening $\mathrm{E}\left[{T}\right]$ do scales as
$\bar{\sigma}^{2}$, while its standard deviation
$\mathrm{std}\left[{T}\right]$ is only linear in $\bar{\sigma}$ (see Fig. 6d).
Note that numerical validations of these scalings can only be obtained through
the use of variance reduction methods, as no convergence was found even for
$N=100^{\prime}000$. Using our antithetic variates, one observes satisfactory
scaling for $N$ as low as $N\simeq 250$. Our trick may find some use if one
wanted to compare the results of perturbative calculations to that of e.g. BEM
simulations.

  * 2.

In contrast with the predictions of Roux et al. (2003), Patinet et al.
(2013b), and Démery et al. (2014b), the effective toughness
$K_{\mathrm{Ic}}^{\mathrm{eff}}$ is not equal to the spatial average
$K_{\mathrm{Ic}}^{0}$ of the toughness field while being in the weak pinning
regime. In the semi-infinite crack limit ($a_{0}/d\rightarrow+\infty$), the
average toughening $\mathrm{E}\left[{T}\right]$ is equal to
$\mathrm{E}\left[{T}\right]=0.5(\sigma/K_{\mathrm{Ic}}^{0})^{2}$ (see Fig.
6a). This toughening contribution is dictated by non-linearities in
$K_{\mathrm{I}}$, as only the terms arising from the second-order expansion of
$K_{\mathrm{I}}$ actually play a role in determining
$\mathrm{E}\left[{T}\right]$ (through the terms involving the convolution
kernel $H$ in Eq. (62)). Moreover, the toughening predicted here in the _weak
pinning_ regime is comparable to that measured by Démery et al. (2014b) and
Démery et al. (2014a) in the _strong pinning_ regime, where crack propagation
is no more smooth and occurs as the succession of depinning instabilities
(Roux et al., 2003). Indeed, the predictions $\mathrm{E}\left[{T}\right]$ of
Eq. (62) for a semi-infinite crack amounts to one fourth of the toughening
measured by Démery et al. (2014b) (see their Eq. (20)), and half that
predicted by Démery et al. (2014a) (see their Eqs. (34) and (35)). We conclude
that one cannot neglect non-linearities in the stress intensity factor
variations when estimating the effective toughness of brittle composites.

This calls for a new derivation of the theoretical predictions of Patinet et
al. (2013b) and Démery et al. (2014b) for a semi-infinite crack in the strong
pinning regime, involving now the second-order terms in the expansion of
$K_{\mathrm{I}}$ or $G$ derived by Vasoya et al. (2013). To do so, one can
follow the route proposed by Démery et al. (2014a), which is based on a
second-order expansion of the front position similar to ours, but taking into
account a finite driving velocity. These predictions should then be compared
to the results of simulations performed in the strong pinning regime, in which
the crack jumps from one equilibrium position to the other. One could build on
the approach recently proposed by Sanner and Pastewka (2022), and derive a
variational formulation of the problem using a third-order expansion of the
elastic energy. The successive state equilibria can then be computed from a
trust-region Newton conjugate-gradient algorithm. Another promising route
would be to extend the numerical method of Rosso and Krauth (2002) to second-
order, but one would need first to assess the validity of Middleton’s
comparison theorem for the non-linear convolution integrals of Vasoya et al.
(2013).

  * 3.

We observe in Fig. 6c and Fig. 6a that the presence of material disorder does
not always result in an overall toughening of the material, in contrast with
the results obtained for the semi-infinite crack (Patinet et al., 2013b;
Démery et al., 2014b; Lebihain, 2021). Indeed, for very small crack radii
$a_{0}\ll d$, the average toughening $\mathrm{E}\left[{T}\right]$ is zero, as
the crack propagates in a single heterogeneity of uniform toughness, whose
probability distribution is given by $\mathcal{P}$ of Eq. (41) and average
value is $K_{\mathrm{Ic}}^{0}$. For larger but still small crack radii
$a_{0}\lesssim 6d$, the average toughening $\mathrm{E}\left[{T}\right]$ is
negative. The presence of material disorder leads to an overall weakening of
the composite from its average toughness $K_{\mathrm{Ic}}^{0}$. On the
contrary, the material is reinforced by it for larger crack radii $a_{0}$, in
agreement with the results obtained for the semi-infinite crack. These
variations are accurately captured by Eq. (62). We show in Fig. 6a the
normalized average $\mathrm{E}\left[{T}\right]/\bar{\sigma}^{2}$ of the
toughening as measured in our numerical simulations. All data points collapse
to the master curve described by Eq. (62). We further observe from Fig. 6a
that the average toughening $\mathrm{E}\left[{T}\right]$ saturates in the
limit of very large crack radii $a_{0}\gg d$. This increase of the macroscopic
toughness with crack length up to a saturation value bears striking
similarities with the R-curve behavior observed in a wide variety of materials
like e.g. concrete (Bažant and Jirásek, 1993), ceramics (Ebrahimi et al.,
2000), rocks (Funatsu et al., 2004), bone (Nalla et al., 2005), compacted
granular materials (Girardot et al., 2023), etc.

A natural consequence of this R-curve behavior is that material disorder may
promote micro-cracking below $\simeq d$. Indeed, the extension of pre-existing
microscopic defects is facilitated by the overall decrease of the effective
toughness up to $a_{0}=d$, while their growth is further stabilized by the
apparent increase of the fracture resistance. Heterogeneous materials may thus
display an increased ductility at the macroscale.

  * 4.

We see in Fig. 6b that the fluctuations of the toughening $T$ (or equivalently
that of the effective toughness $K_{\mathrm{Ic}}^{\mathrm{eff}}$) decreases
with increasing crack radius $a_{0}$. The dependence of
$\mathrm{std}\left[{T}\right]$ in $a_{0}/d$ is well captured by our theory, as
one can see from Fig. 6b, where numerical estimates of the normalized
fluctuations $\mathrm{std}\left[{T}\right]/\bar{\sigma}$ all collapse to the
master curve described by Eq. (63). Again, one notices that
$\mathrm{std}\left[{K_{\mathrm{Ic}}^{\mathrm{eff}}}\right]=\sigma=\mathrm{std}\left[{K_{\mathrm{Ic}}}\right]$
for small crack radii $a_{0}/d\ll 1$, as the crack propagates in a single
heterogeneity only. Moreover, for large crack radii $a_{0}/d\gg 1$,
$\mathrm{std}\left[{T}\right]\propto\sigma(a_{0}/d)^{-1/2}$. This finite size
effect may be reminiscent of having Larkin domains smaller than the front
perimeter (Démery et al., 2014b). As
$\mathrm{std}\left[{K_{\mathrm{Ic}}^{\mathrm{eff}}}\right]\underset{a_{0}/d\rightarrow+\infty}{\rightarrow}0$,
it is tempting to think that our theory proves the existence of an _intrinsic_
(deterministic) effective toughness for weakly disordered brittle composites
in the limit $a_{0}/d\rightarrow+\infty$. However, one must keep in mind that
Eq. (63) only provides a first-order expansion of the standard deviation
$\mathrm{std}\left[{T}\right]$ of the toughening, so that the second-order
fluctuations of the effective toughness are not necessarily zero in that
limit. To access these second-order fluctuations, one should derive a third-
order expansion of the variance $\mathrm{Var}\left[{T}\right]$ of the
toughening. Note that this is feasible here, even within a second-order
expansion of $K_{\mathrm{I}}$, but it would involve the three-points
correlation function of $f$, a quantity that is often poorly constrained and
difficult to compute (Jiao et al., 2009). Those calculations are out-of-scope
of the present study, and are left for further works.

  * 5.

The expressions of the average of the toughening $\mathrm{E}\left[{T}\right]$
and its standard deviation $\mathrm{std}\left[{T}\right]$ involve the
functions $F$ and $G$ that relate to the two-points correlation function
$\mathcal{F}$ of the toughness fluctuations $f$. This was already noticed by
Démery et al. (2014a) for the semi-infinite crack, and it suggests that
playing on the spatial structure of the toughness field could potentially lead
to a change in material toughening. An interesting route would be to
investigate cases where the correlation length $d_{z}$ of the toughness field
in the direction $z$ is different from that $d_{x}$ in the direction $x$. One
could then expand the front shape $a(\theta)$ with the disorder intensity
$\bar{\sigma}$, but also with a small parameter $\epsilon=d_{x}/d_{z}-1$ that
quantifies the toughness anisotropy. One could then observe how it modifies
both the average front shape, and the overall toughening of the material. Note
that such an expansion is only possible within our second-order theory, as
terms of first-order in $\epsilon$ are necessarily zero (absence of disorder).

  * 6.

$\mathrm{E}\left[{T}\right]$ involves sums over the Fourier convolution
kernels $L(k)$ and $H(k,k^{\prime})$. As the latter both depend the reference
SIF $K_{\mathrm{I}}^{0}$ generated by the loading, and its derivatives
$\partial K_{\mathrm{I}}^{0}/\partial a$ and
$\partial^{2}K_{\mathrm{I}}^{0}/\partial a^{2}$ with crack advance (see Eq.
(25)), one expects that the effective toughness
$K_{\mathrm{Ic}}^{\mathrm{eff}}$ generally depends on the loading conditions.
This was already noted by Démery et al. (2014a) and Lebihain et al. (2021) in
the case of the semi-infinite crack. This dependency breaks down as soon as
the structural length scales set by both the crack size $a_{0}$ and the
loading conditions:

$\mathcal{L}_{1}=K_{\mathrm{I}}^{0}/(dK_{\mathrm{I}}^{0}/da)\text{ and
}\mathcal{L}_{2}=\sqrt{K_{\mathrm{I}}^{0}/(d^{2}K_{\mathrm{I}}^{0}/da^{2})}.$
(66)

get much larger than the heterogeneity size $d$. Note however that the
standard deviation of the toughening $\mathrm{std}\left[{T}\right]$ given in
Eq. (63) is independent of the loading conditions.

Overall, our findings shed new lights on the toughening of brittle composites
by heterogeneities, as well as on the size effects emerging from the
finiteness of the crack geometry. However, one must keep in mind that they are
only valid as long as a second-order expansion of $K_{\mathrm{Ic}}$ remains
valid, i.e. when the front perturbation $\delta a$ is smaller than the
characteristic length scale $d$ of the toughness field ($\delta a\lesssim d$).
As a result, one should expect departure of the numerical results and our
theory as soon as the crack enters the so-called _strong pinning_ regime, in
which crack propagation articulates as the succession of depinning
instabilities (Démery et al., 2014b). But one should also be aware that taking
into account size effects and non-linearities in the strong pinning regime may
change the conclusions obtained by Patinet et al. (2013b) and Démery et al.
(2014b).

##  5 Conclusion

This paper address the second-order variation of the local mode I stress
intensity factor (SIF) arising from a coplanar perturbation of the front of a
penny-shaped crack embedded in some isotropically linearly elastic infinite
body. We build next on our model to determine material reinforcement arising
from local fluctuations of toughness, and account for both non-linear and
finite-size effects.

Section 2 was devoted to the derivation of our theoretical model, following
the approach proposed by Rice (1989) and used by Leblond et al. (2012) for the
semi-infinite crack. First, we derived the first-order variations of the
fundamental kernel $Z$ arising from the perturbation of the crack front from
its reference circular configuration. Next, we determined the first-order
variations of the derivative $\partial K_{\mathrm{I}}/\partial A$ of the mode
I SIF $K_{\mathrm{I}}$ with respect to the amplitude $A$ of the perturbation,
which ultimately yields our second-order formula of Eqs. (24-25) by
integration. This was performed in the case where the unperturbed crack is
loaded by some axisymmetric system of forces, but the general derivation was
proposed in B.

In Section 3, we validated our second-order theory by comparing its
predictions to the second-order expansion of an analytical formula derived by
Irwin (1958) for the SIF distribution along the front of an elliptical crack
loaded by a uniform tensile stress. The validity range of our model was then
assessed by comparing its output to the local SIF distributions along wavy
cracks of varying mode and amplitude, which were computed using the numerical
method of Lazarus (2003). We showed that the predicting abilities of our
second-order theory are improved from the linear model of Gao and Rice
(1987b), but that errors are increasing with both the perturbation mode and
amplitude.

In Section 4, we use our non-linear model to investigate crack propagation in
heterogeneous media, in which the local toughness was fluctuating at random.
We only considered moderate fluctuations of material toughness, so that crack
propagation was stable and occurred as the continuous succession of
equilibrium positions (weak pinning regime). The toughening arising from
toughness fluctuations has been widely investigated in the literature for the
semi-infinite coplanar crack (Roux et al., 2003), even beyond the case of weak
pinning, (Patinet et al., 2013b; Démery et al., 2014b; Lebihain, 2021).
However, these authors predicted a zero toughening for weak pinning, building
a _first-order_ expansion of the local stress intensity factor
$K_{\mathrm{I}}$. Here, we showed that this was no more true as soon as one
considers non-linearities in the expression of $K_{\mathrm{I}}$. More
precisely, we showed that the effective toughness of disordered brittle media
is a random variable, whose statistical descriptors, like e.g. its average and
its variance, can be predicted rigorously at second order in $\bar{\sigma}$
using our non-linear theory. We validated our theoretical predictions on the
results of $\sim 1$ million numerical simulations of crack propagation in
weakly disordered materials. One of our key findings is that material
toughening by heterogeneities of typical size $d$ is _size-dependent_. Namely,
cracks of average radius $\left<a\right>$ smaller than $\simeq 6\,d$ are
experiencing an effective toughness $K_{\mathrm{Ic}}^{\mathrm{eff}}$ that is
lower than the average $K_{\mathrm{Ic}}^{0}$ of the toughness field. On the
contrary, larger cracks propagates in a seemingly tougher material
($K_{\mathrm{Ic}}^{\mathrm{eff}}>K_{\mathrm{Ic}}^{0}$), in accordance with
previous results obtained for the semi-infinite crack. It leads to an apparent
R-curve behavior of the brittle composite at the macroscale. A natural
consequence of this observation is that material disorder may promote the
development of micro-cracks at the heterogeneity scale, while preventing them
to extend further, leading to an overall increase in material ductility. We
also proved that material toughening scales, on average, as the square
$\bar{\sigma}^{2}$ of the fluctuations amplitude, and that non-linear terms in
the expansion of $K_{\mathrm{I}}$ do play a role in the expression of the
associated prefactor. Moreover, this contribution in toughening is comparable
to that emerging from the intermittent dynamics of crack front at moderate
disorder intensity (Démery et al., 2014a). This drive the needs to revisit the
results of Patinet et al. (2013b) and Démery et al. (2014b) in a non-linear
setting, but appropriate numerical methods must be developed to validate
theoretical predictions.

Our work also raises interesting perspectives to model crack propagation in a
multi-physics environment. First, we provided here a second-order expansion of
the stress intensity factor along quasi-circular crack fronts. However, the
perturbative approach of Rice (1989) also yields the first-order variations of
local aperture, a quantity that has rarely been used in the past. One could
use a reasoning similar to that of Leblond et al. (2012), and derive its
second-order expansion. Combined with that of the crack face weight function,
it should provide means to investigate crack propagation in heterogeneous yet
cohesive materials, whose strength $\sigma_{\mathrm{c}}$ depends on the local
crack opening (Lebihain et al., 2022). The knowledge of the aperture of
perturbed cracks also gives access to the crack volume, which may be relevant
to describe 3D hydraulic fracturing (Savitski and Detournay, 2002) and its
interplay with material heterogeneities. Second, one could extend our mode I
calculations to mixed mode I+II+III using the perturbative framework of Favier
et al. (2006a). Preliminary calculations show that it should allow
rationalizing the quasi-elliptical front shapes of fluid-driven shear ruptures
propagating along frictional interfaces governed by Coulomb’s friction (Sáez
et al., 2022). Third, the knowledge of the perturbation $\delta a_{*}$, which
has been derived in Section 2.1 to determine the first-order variation of the
fundamental kernel, permits us to extend our second-order model to nearly
circular connections (Gao and Rice, 1987a) at relatively small cost. One may
then apply this theory to investigate adhesion along chemically heterogeneous
interfaces (Sanner and Pastewka, 2022) or contact between a soft material and
a rough indenter (Argatov, 2021).

## CRediT authorship contribution statement

Mathias Lebihain: Conceptualization, Investigation, Formal analysis, Numerical
simulations, Software, Visualization, Writing - Original Draft. Manish Vasoya:
Investigation, Numerical simulations, Writing - Review & Editing. Véronique
Lazarus: Investigation, Software, Writing - Review & Editing.

## Online access to data

The Python script used to generate the figures and the associated data are
available on the French open data plateform recherche.data.gouv at
http://dx.doi.org/10.57745/FJFLYB.

## Acknowledgments

The authors gratefully thank Jean-Baptiste Leblond for stimulating discussions
on the second-order derivation, Antoine Sanner and Lars Pastewka for fruitful
discussions on the extension of our approach to other physical systems and on
variational formulations, Sébastien Brisard for suggesting variance reduction
methods. We also thank the anonymous reviewer for their constructive and
insightful comments. ML thanks Antoine Sanner for a critical reading of the
manuscript.

For the purpose of Open Access, a CC-BY public copyright licence has been
applied by the authors to the present document and will be applied to all
subsequent versions up to the Author Accepted Manuscript arising from this
submission.

Distributed under a Creative Commons Attribution — 4.0 International licence

##  Appendix A Calculation of the first-order variations of the kernel $\delta
Z$

In this Appendix, our aim is to provide useful indications on the derivation
of Eq. (15) that gives the first-order variations $\delta Z$ of the kernel,
starting from the general formula (8). Following Eqs. (4) and (9), the formula
(8) writes as:

$\displaystyle\begin{aligned} \delta
Z(\theta_{1},\theta_{2})&=Z^{0}(a_{0};\theta_{1},\theta_{2})\left[\left(\dfrac{\mathcal{D}(M_{1},M_{2})}{\mathcal{D}(M^{0}_{1},M^{0}_{2})}\right)^{-2}-1\right]+\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)Z^{0}(a_{0};\theta_{2},\theta)\left[\delta
a(\theta)-\delta a_{*}(\theta)\right]a_{0}d\theta+\mathcal{O}(A^{2})\\\
&=-\dfrac{\delta a(\theta_{1})+\delta a(\theta_{2})}{8\pi
a_{0}^{3}\sin^{2}\left[(\theta_{2}-\theta_{1})/2\right]}+\dfrac{1}{64\pi^{2}a_{0}^{3}}\mathrm{PV}\int_{0}^{2\pi}\dfrac{1}{\sin^{2}[(\theta-\theta_{1})/2]\sin^{2}[(\theta-\theta_{2})/2]}\left[\delta
a(\theta)-\delta
a_{*}(\theta)\right]a_{0}d\theta+\mathcal{O}(A^{2}).\end{aligned}$

(67)

Next, we use the first-order expansion of $\delta a_{*}$ of Eq. (13), and the
decomposition of the sinus fraction given in Eq. (14):

$\displaystyle\begin{aligned} \delta Z(\theta_{1},\theta_{2})=&-\dfrac{\delta
a(\theta_{1})+\delta a(\theta_{2})}{8\pi
a_{0}^{3}\sin^{2}\left[(\theta_{2}-\theta_{1})/2\right]}+\dfrac{1}{(8\pi)^{2}a_{0}^{3}\sin^{2}\left[(\theta_{2}-\theta_{1})/2\right]}\times\left\\{\rule{0.0pt}{28.45274pt}\right.\\\
&\mathrm{PV}\int_{0}^{2\pi}\dfrac{\cos^{2}[(\theta-\theta_{1})/2]}{\sin^{2}[(\theta-\theta_{1})/2]}\left[\delta
a(\theta)-\delta
a(\theta_{1})\right]d\theta+\mathrm{PV}\int_{0}^{2\pi}\dfrac{\cos^{2}[(\theta-\theta_{2})/2]}{\sin^{2}[(\theta-\theta_{2})/2]}\left[\delta
a(\theta)-\delta a(\theta_{2})\right]d\theta\\\ &-\dfrac{\delta
a(\theta_{2})-\delta
a(\theta_{1})}{2}\cdot\dfrac{\cos[(\theta_{2}-\theta_{1})/2]}{\sin[(\theta_{2}-\theta_{1})/2]}\cdot\left(\mathrm{PV}\int_{0}^{2\pi}\underbrace{\dfrac{\cos^{2}[(\theta-\theta_{1})/2]}{\sin^{2}[(\theta-\theta_{1})/2]}\sin(\theta-\theta_{1})}_{\text{odd
integrand}}d\theta+\mathrm{PV}\int_{0}^{2\pi}\underbrace{\dfrac{\cos^{2}[(\theta-\theta_{2})/2]}{\sin^{2}[(\theta-\theta_{2})/2]}\sin(\theta-\theta_{2})}_{\text{odd
integrand}}d\theta\right)\\\ &-\dfrac{\delta a(\theta_{2})-\delta
a(\theta_{1})}{2}\cdot\left(\underbrace{\int_{0}^{2\pi}\cos^{2}[(\theta-\theta_{1})/2]d\theta-\int_{0}^{2\pi}\cos^{2}[(\theta-\theta_{2})/2]d\theta}_{\text{cancel
out}}\right)\\\
&+\dfrac{2\cos[(\theta_{2}-\theta_{1})/2]}{\sin[(\theta_{2}-\theta_{1})/2]}\cdot\mathrm{PV}\int_{0}^{2\pi}\left[\dfrac{\cos^{3}[(\theta-\theta_{1})/2]}{\sin[(\theta-\theta_{1})/2]}-\dfrac{\cos^{3}[(\theta-\theta_{2})/2]}{\sin[(\theta-\theta_{2})/2]}\right]\delta
a(\theta)d\theta\\\
&-\left(\mathrm{PV}\int_{0}^{2\pi}\underbrace{\dfrac{\cos^{3}[(\theta-\theta_{1})/2]}{\sin[(\theta-\theta_{1})/2]}}_{\text{odd
integrand}}d\theta\right)\delta
a(\theta_{1})-\left(\mathrm{PV}\int_{0}^{2\pi}\underbrace{\dfrac{\cos^{3}[(\theta-\theta_{2})/2]}{\sin[(\theta-\theta_{2})/2]}}_{\text{odd
integrand}}d\theta\right)\delta a(\theta_{2})\\\ &-\dfrac{\delta
a(\theta_{2})-\delta
a(\theta_{1})}{2}\cdot\dfrac{2\cos^{2}[(\theta_{2}-\theta_{1})/2]}{\sin^{2}[(\theta_{2}-\theta_{1})/2]}\cdot\left(\underbrace{\int_{0}^{2\pi}\dfrac{\cos^{3}[(\theta-\theta_{1})/2]}{\sin[(\theta-\theta_{1})/2]}\sin(\theta-\theta_{1})d\theta-\int_{0}^{2\pi}\dfrac{\cos^{3}[(\theta-\theta_{2})/2]}{\sin[(\theta-\theta_{2})/2]}\sin(\theta-\theta_{2})d\theta}_{\text{cancel
out}}\right)\\\ &-\dfrac{\delta a(\theta_{2})-\delta
a(\theta_{1})}{2}\cdot\dfrac{2\cos[(\theta_{2}-\theta_{1})/2]}{\sin[(\theta_{2}-\theta_{1})/2]}\cdot\left(\int_{0}^{2\pi}\underbrace{\dfrac{\cos^{3}[(\theta-\theta_{1})/2]}{\sin[(\theta-\theta_{1})/2]}2\sin^{2}[(\theta-\theta_{1})/2]}_{\text{odd
integrand}}d\theta+\int_{0}^{2\pi}\underbrace{\dfrac{\cos^{3}[(\theta-\theta_{2})/2]}{\sin[(\theta-\theta_{2})/2]}2\sin^{2}[(\theta-\theta_{2})/2]}_{\text{odd
integrand}}d\theta\right)\\\
&+\int_{0}^{2\pi}2\cos^{2}[(\theta-\theta_{1})/2]\left[\delta a(\theta)-\delta
a(\theta_{1})\right]d\theta+\int_{0}^{2\pi}2\cos^{2}[(\theta-\theta_{2})/2]\left[\delta
a(\theta)-\delta a(\theta_{2})\right]d\theta\\\ &-\dfrac{\delta
a(\theta_{2})-\delta
a(\theta_{1})}{2}\cdot\dfrac{\cos[(\theta_{2}-\theta_{1})/2]}{\sin[(\theta_{2}-\theta_{1})/2]}\left(\int_{0}^{2\pi}2\underbrace{\cos^{2}[(\theta-\theta_{1})/2]\sin(\theta-\theta_{1})}_{\text{odd
integrand}}d\theta+\int_{0}^{2\pi}2\underbrace{\cos^{2}[(\theta-\theta_{2})/2]\sin(\theta-\theta_{2})}_{\text{odd
integrand}}d\theta\right)\\\ &\left.-\dfrac{\delta a(\theta_{2})-\delta
a(\theta_{1})}{2}\cdot\left(\underbrace{\int_{0}^{2\pi}4\cos^{2}[(\theta-\theta_{1})/2]\sin^{2}[(\theta-\theta_{1})/2]d\theta-\int_{0}^{2\pi}4\cos^{2}[(\theta-\theta_{2})/2]\sin^{2}[(\theta-\theta_{2})/2]d\theta}_{\text{cancel
out}}\right)\right\\}+\mathcal{O}(A^{2}).\end{aligned}$

(68)

After simplifications of the integrals that are zero due to the oddness of the
integrand, and those that cancel out, one has:

$\displaystyle\begin{aligned} \delta Z(\theta_{1},\theta_{2})=&-\dfrac{\delta
a(\theta_{1})+\delta a(\theta_{2})}{8\pi
a_{0}^{3}\sin^{2}\left[(\theta_{2}-\theta_{1})/2\right]}+\dfrac{1}{(8\pi)^{2}a_{0}^{3}\sin^{2}\left[(\theta_{2}-\theta_{1})/2\right]}\times\Bigg{\\{}\\\
&\hskip
42.67912pt\mathrm{PV}\int_{0}^{2\pi}\dfrac{\cos^{2}[(\theta-\theta_{1})/2]}{\sin^{2}[(\theta-\theta_{1})/2]}\left[\delta
a(\theta)-\delta
a(\theta_{1})\right]d\theta+\mathrm{PV}\int_{0}^{2\pi}\dfrac{\cos^{2}[(\theta-\theta_{2})/2]}{\sin^{2}[(\theta-\theta_{2})/2]}\left[\delta
a(\theta)-\delta a(\theta_{2})\right]d\theta\\\ &\hskip
42.67912pt+\dfrac{2\cos[(\theta_{2}-\theta_{1})/2]}{\sin[(\theta_{2}-\theta_{1})/2]}\cdot\mathrm{PV}\int_{0}^{2\pi}\left[\dfrac{\cos^{3}[(\theta-\theta_{1})/2]}{\sin[(\theta-\theta_{1})/2]}-\dfrac{\cos^{3}[(\theta-\theta_{2})/2]}{\sin[(\theta-\theta_{2})/2]}\right]\delta
a(\theta)d\theta\\\ &\hskip
42.67912pt\left.+\int_{0}^{2\pi}2\left[\cos^{2}\left(\dfrac{\theta-\theta_{1}}{2}\right)+\cos^{2}\left(\dfrac{\theta-\theta_{2}}{2}\right)\right]\delta
a(\theta)d\theta-2\pi\left[\delta a(\theta_{1})+\delta
a(\theta_{2})\right]\right\\}+\mathcal{O}(A^{2}),\end{aligned}$

(69)

from which one retrieves Eq. (15). Note that the simplifications in Eq. (68)
arise from the careful choice of $\delta a_{*}$ in Eq. (10) and that of
$\theta_{*}$ in Eq. (11), which allows for the alternate expressions of
$\delta a_{*}$ in Eq. (13) that pair well with the decomposition of Eq. (14).

##  Appendix B Calculation of the second-order variations of the SIF $\delta
K$

In this Appendix, we provide useful indications on how to calculate the
second-order $K_{\mathrm{I}}^{2}$ perturbations of the mode I stress intensity
factor arising from perturbations $\delta a=A\phi$ of the circular front of
radius $a_{0}$. For the sake of completeness, we do it here in the general
case where the reference SIF $K_{\mathrm{I}}^{0}$, applied to the unperturbed
circular crack, depends on $\theta$. The calculations presented in the main
text only address cases where $K_{\mathrm{I}}^{0}$ is axisymmetric.

As established in Eq. (18), the SIF variations arising from an additional
infinitesimal perturbation $\delta A\phi$ reads:

$\displaystyle K_{\mathrm{I}}(a_{0},\left[(A+\delta
A)\phi\right];\theta_{1})-K_{\mathrm{I}}(a_{0},\left[A\phi\right];\theta_{1})=K_{\mathrm{I}}(a_{0}+\delta
A\phi(\theta_{1}),\left[A\phi\right];a_{0};\theta_{1})-K_{\mathrm{I}}(a_{0},\left[A\phi\right];\theta_{1})$
(70) $\displaystyle\hskip 56.9055pt+\mathrm{PV}\int_{0}^{2\pi}Z(a_{0}+\delta
A\phi(\theta_{1}),\left[A\phi\right];\theta_{1},\theta)K_{\mathrm{I}}(a_{0}+\delta
A\phi(\theta_{1}),\left[A\phi\right];\theta)\left[\delta A\phi(\theta)-\delta
A\phi(\theta_{1})\right]ds(\theta)+\mathcal{O}(\delta A^{2}),$

where $s$ is the curvilinear abscissa along the perturbed front
$\mathcal{F^{\prime}}$, so that $ds$ reduces here to:

${\color[rgb]{0,0,0}ds(\theta)=a_{0}+\delta
A\phi(\theta_{1})+A\phi(\theta)+\mathcal{O}(A^{2}),}$ (71)

Using Eq. (2), one has:

$\displaystyle K_{\mathrm{I}}(a_{0}+\delta
A\phi(\theta_{1}),\left[A\phi\right];a_{0};\theta_{1})-K_{\mathrm{I}}(a_{0},\left[A\phi\right];\theta_{1})$
(72) $\displaystyle\hskip 56.9055pt=\delta A\times\left\\{\dfrac{\partial
K_{\mathrm{I}}^{0}}{\partial
a}(a_{0},\theta_{1})\phi(\theta_{1})+A\dfrac{\partial^{2}K_{\mathrm{I}}^{0}}{\partial
a^{2}}(a_{0},\theta_{1})\phi(\theta_{1})^{2}\right.$ $\displaystyle\hskip
85.35826pt+A\,\mathrm{PV}\int_{0}^{2\pi}\dfrac{\partial Z^{0}}{\partial
a}(a_{0};\theta_{1},\theta)K_{\mathrm{I}}^{0}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]\phi(\theta_{1})a_{0}d\theta$
$\displaystyle\hskip
85.35826pt+A\,\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)\dfrac{\partial
K_{\mathrm{I}}^{0}}{\partial
a}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]\phi(\theta_{1})a_{0}d\theta$
$\displaystyle\hskip
85.35826pt\left.+A\,\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)K_{\mathrm{I}}^{0}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]\phi(\theta_{1})d\theta+\mathcal{O}(A^{2})\right\\}+\mathcal{O}(\delta
A^{2}).$

The calculation of the second term in the right-hand side of Eq. (70) is also
straightforward. From Eqs. (16) and (71), one has:

$\displaystyle\mathrm{PV}\int_{0}^{2\pi}$ $\displaystyle Z(a_{0}+\delta
A\phi(\theta_{1}),\left[A\phi\right];\theta_{1},\theta)K_{\mathrm{I}}(a_{0},\left[A\phi\right];\theta)\left[\delta
A\phi(\theta)-\delta A\phi(\theta_{1})\right]ds(\theta)$ (73)
$\displaystyle=\mathrm{PV}\int_{0}^{2\pi}$ $\displaystyle
Z(a_{0},\left[A\phi\right];\theta_{1},\theta)K_{\mathrm{I}}(A;a_{0};\theta)\left[\delta
A\phi(\theta)-\delta A\phi(\theta_{1})\right]ds(\theta)+\mathcal{O}(\delta
A^{2})$ $\displaystyle=\delta
A\times\left\\{\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)K_{\mathrm{I}}^{0}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]a_{0}d\theta\right.$
$\displaystyle\hskip
17.07164pt+A\,\mathrm{PV}\int_{0}^{2\pi}Z^{1}(a_{0},\phi;\theta_{1},\theta)K_{\mathrm{I}}^{0}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]a_{0}d\theta$
$\displaystyle\hskip
17.07164pt+A\,\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)\dfrac{\partial
K_{\mathrm{I}}^{0}}{\partial
a}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]\phi(\theta)a_{0}d\theta$
$\displaystyle\hskip
17.07164pt+A\,\mathrm{PV}\int_{0}^{2\pi}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)Z^{0}(a_{0};\theta,\theta^{\prime})K_{\mathrm{I}}^{0}(a_{0},\theta^{\prime})\left[\phi(\theta)-\phi(\theta_{1})\right]\left[\phi(\theta^{\prime})-\phi(\theta)\right]a_{0}^{2}d\theta
d\theta^{\prime}$ $\displaystyle\hskip
17.07164pt\left.+A\,\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)K_{\mathrm{I}}^{0}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]\phi(\theta)d\theta+\mathcal{O}(A^{2})\right\\}+\mathcal{O}(\delta
A^{2}).$

Identifying the terms in $\delta A$ in the expansions of Eqs. (17) and (70),
one finds:

$K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})=\dfrac{\partial
K_{\mathrm{I}}^{0}}{\partial
a}(a_{0},\theta_{1})\phi(\theta_{1})+\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)K_{\mathrm{I}}^{0}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]a_{0}d\theta,$
(74)

which is Eq. (2) of Gao and Rice (1987b). The terms in $\delta A\times A$
yield:

$\displaystyle
2K_{\mathrm{I}}^{2}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})=$
$\displaystyle\dfrac{\partial^{2}K_{\mathrm{I}}^{0}}{\partial
a^{2}}(a_{0},\theta_{1})\phi(\theta_{1})^{2}+\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)\dfrac{\partial
K_{\mathrm{I}}^{0}}{\partial
a}(a_{0},\theta)\left[\phi(\theta)^{2}-\phi(\theta_{1})^{2}\right]a_{0}d\theta$
(75)
$\displaystyle+\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)K_{\mathrm{I}}^{0}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]^{2}d\theta$
$\displaystyle+\mathrm{PV}\int_{0}^{2\pi}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)Z^{0}(a_{0};\theta,\theta^{\prime})K_{\mathrm{I}}^{0}(a_{0},\theta^{\prime})\left[\phi(\theta)-\phi(\theta_{1})\right]\left[\phi(\theta)-\phi(\theta^{\prime})\right]a_{0}^{2}d\theta
d\theta^{\prime}$
$\displaystyle+\mathrm{PV}\int_{0}^{2\pi}Z^{1}(a_{0},\left[\phi\right];\theta_{1},\theta)K_{\mathrm{I}}^{0}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]a_{0}d\theta,$

which is Eq. (20). Here we used $\partial Z^{0}/\partial a=-2Z^{0}/a_{0}$ from
Eq. (4).

We need now to calculate Eqs. (74) and (75) by resorting to Fourier series.
The calculations are heavy, and we will show them here only for
$K_{\mathrm{I}}^{1}$. Eq. (74) yields:

$\displaystyle\begin{aligned}
K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})&=\dfrac{\partial
K_{\mathrm{I}}^{0}}{\partial
a}(a_{0},\theta_{1})\phi(\theta_{1})+\mathrm{PV}\int_{0}^{2\pi}Z^{0}(a_{0};\theta_{1},\theta)K_{\mathrm{I}}^{0}(a_{0},\theta)\left[\phi(\theta)-\phi(\theta_{1})\right]a_{0}d\theta\\\
&=\sum_{k,k^{\prime}}\widehat{\left(\dfrac{\partial
K_{\mathrm{I}}^{0}}{\partial
a}\right)}_{k}\\!(a_{0})\,\hat{\phi}_{k^{\prime}}e^{i(k+k^{\prime})\theta_{1}}+\sum_{k,k^{\prime}}\widehat{K_{\mathrm{I}k}^{0}}(a_{0})\hat{\phi}_{k^{\prime}}\left[\dfrac{1}{8a_{0}}\mathrm{PV}\int_{0}^{2\pi}\dfrac{e^{ik^{\prime}(\theta-\theta_{1})}-1}{\sin^{2}[(\theta-\theta_{1})/2]}e^{ik(\theta-\theta_{1})}d\theta\right]e^{i(k+k^{\prime})\theta_{1}}\\\
&=\sum_{k,k^{\prime}}\left[\widehat{\left(\dfrac{\partial
K_{\mathrm{I}}^{0}}{\partial
a}\right)}_{k}\\!(a_{0})-\dfrac{|{k+k^{\prime}}|-|{k}|}{2a_{0}}\widehat{K_{\mathrm{I}k}^{0}}(a_{0})\right]\hat{\phi}_{k^{\prime}}e^{i(k+k^{\prime})\theta_{1}},\end{aligned}$

(76)

where $\hat{\phi}_{k}$, $\widehat{K_{\mathrm{I}k}^{0}}(a_{0})$, and
$\widehat{\left(\dfrac{\partial K_{\mathrm{I}}^{0}}{\partial
a}\right)}_{k}\\!(a_{0})$ are the $k$-th Fourier coefficient of $\phi$,
$K_{\mathrm{I}}^{0}(a_{0},\cdot)$, and $\dfrac{\partial
K_{\mathrm{I}}}{\partial a}(a_{0},\cdot)$ respectively. The terms into brack
in Eq. (76) reduces to $L(k)$ as given in Eq. (25a) when $K_{\mathrm{I}}^{0}$
is independent of $\theta$. In deriving (76), we used the following integral:

$\mathrm{PV}\int_{0}^{2\pi}\dfrac{e^{iku}-1}{\sin^{2}(u/2)}du=-4\pi|{k}|,$
(77)

which is calculated using the residue theorem. The calculations for
$K_{\mathrm{I}}^{2}$ are similar, although more complicated as they involve
triple sums. They can be simplified by building on the following expression of
the first-order kernel variations $Z^{1}$:

$\displaystyle\begin{aligned}
Z^{1}(a_{0},\left[\phi\right];\theta_{1},\theta_{2})=-\dfrac{1}{16\pi
a_{0}^{2}\sin^{2}\left[(\theta_{2}-\theta_{1})/2\right]}\times\left[\rule{0.0pt}{17.07182pt}\right.&\sum_{k}\hat{\phi}_{k}\left(e^{ik\theta_{2}}+e^{ik\theta_{1}}\right)\left(|{k}|+2-\dfrac{\Delta_{|{k}|1}}{4}\right)\\\
&+\left.i\,\dfrac{\cos\left[(\theta_{2}-\theta_{1})/2\right]}{\sin\left[(\theta_{2}-\theta_{1})/2\right]}\sum_{k}\hat{\phi}_{k}\left(e^{ik\theta_{2}}-e^{ik\theta_{1}}\right)\left(1-\Delta_{k0}-\dfrac{\Delta_{|{k}|1}}{4}\right)\mathrm{sign}(k)\right],\end{aligned}$

(78)

as well as the following integrals:

$\displaystyle\mathrm{PV}\int_{0}^{2\pi}\dfrac{2\cos^{3}(u/2)}{\sin(u/2)}e^{iku}du=4i\pi\,\mathrm{sign}(k)\left(1-\dfrac{1}{4}\Delta_{|{k}|1}-\Delta_{k0}\right)$
(79)
$\displaystyle\mathrm{PV}\int_{0}^{2\pi}\dfrac{\cos(u/2)}{\sin^{3}(u/2)}(e^{iku}-1)(e^{ik^{\prime}u}-1)du=-4i\pi\left[k(|{k+k^{\prime}}|-|{k}|)+k^{\prime}(|{k+k^{\prime}}|-|{k^{\prime}}|)\right],$

where $\Delta_{ij}$ is the Kronecker symbol that is equal to $1$ when $i=j$
and $0$ otherwise, and $\mathrm{sign}(k)$ is the sign of the integer $k$.
Using Eqs. (77), (78), and (79), Eq. (75) yields:

$K_{\mathrm{I}}^{2}(a_{0},[K_{\mathrm{I}}^{0}],[\phi];\theta_{1})=\sum_{k,k^{\prime},k^{\prime\prime}}\mathcal{H}\left(\left[K_{\mathrm{I}}^{0}\right];k,k^{\prime},k^{\prime\prime}\right)\hat{\phi}_{k^{\prime}}\hat{\phi}_{k^{\prime\prime}}e^{i(k+k^{\prime}+k^{\prime\prime})\theta_{1}},$
(80)

where the convolution kernel $\mathcal{K}$ reads:

$\displaystyle\begin{aligned}
\mathcal{H}\left(\left[K_{\mathrm{I}}^{0}\right];k,k^{\prime},k^{\prime\prime}\right)=&\dfrac{1}{2}\widehat{\left(\dfrac{\partial^{2}K_{\mathrm{I}}^{0}}{\partial
a^{2}}\right)}_{k}\\!(a_{0})-\dfrac{1}{4a_{0}}\widehat{\left(\dfrac{\partial
K_{\mathrm{I}}^{0}}{\partial
a}\right)}_{k}\\!(a_{0})\left(|{k+k^{\prime}+k^{\prime\prime}}|-|{k}|\right)-\dfrac{1}{4a_{0}^{2}}\widehat{K_{\mathrm{I}k}^{0}}(a_{0})\left(|{k+k^{\prime}+k^{\prime\prime}}|-|{k+k^{\prime}}|-|{k+k^{\prime\prime}}|+|{k}|\right)\\\
&-\dfrac{1}{16a_{0}^{2}}\widehat{K_{\mathrm{I}k}^{0}}(a_{0})\left((k+k^{\prime})^{2}+(k+k^{\prime\prime})^{2}\right)+\dfrac{1}{16a_{0}^{2}}\widehat{K_{\mathrm{I}k}^{0}}(a_{0})\left(|{k+k^{\prime}}|+|{k+k^{\prime\prime}}|\right)\left(|{k+k^{\prime}+k^{\prime\prime}}|+|{k}|\right)\\\
&-\dfrac{1}{8a_{0}^{2}}\widehat{K_{\mathrm{I}k}^{0}}(a_{0})|{k}||{k+k^{\prime}+k^{\prime\prime}}|+\dfrac{1}{4a_{0}^{2}}\widehat{K_{\mathrm{I}k}^{0}}(a_{0})\left(|{k+k^{\prime}+k^{\prime\prime}}|-|{k}|\right)+\dfrac{1}{16a_{0}^{2}}\widehat{K_{\mathrm{I}k}^{0}}(a_{0})\left(|{k^{\prime}}||{k+k^{\prime\prime}}|+|{k^{\prime\prime}}||{k+k^{\prime}}|\right)\\\
&-\dfrac{1}{16a_{0}^{2}}\widehat{K_{\mathrm{I}k}^{0}}(a_{0})|{k}|\left(|{k^{\prime}}|+|{k^{\prime\prime}}|\right)-\dfrac{1}{16a_{0}^{2}}\widehat{K_{\mathrm{I}k}^{0}}(a_{0})k\left(|{k+k^{\prime}+k^{\prime\prime}}|-|{k+k^{\prime}}|-|{k+k^{\prime\prime}}|+|{k}|\right)\left(\mathrm{sign}(k^{\prime})+\mathrm{sign}(k^{\prime\prime})\right)\\\
&-\dfrac{1}{16a_{0}^{2}}\widehat{K_{\mathrm{I}k}^{0}}(a_{0})\left(|{k+k^{\prime}+k^{\prime\prime}}|(|{k^{\prime}}|+|{k^{\prime\prime}}|)-|{k^{\prime}}||{k+k^{\prime}}|-|{k^{\prime\prime}}||{k+k^{\prime\prime}}|\right)\mathrm{sign}(k^{\prime}k^{\prime\prime}).\end{aligned}$

(81)

The set of equations (B.11-12) is validated next in C on elliptical cracks
loaded by a non-axisymmetric loading. Considering the limiting case of a
reference SIF $K_{\mathrm{I}}^{0}$ independent of $\theta$, one has:

$\displaystyle\mathcal{H}\left(\left[K_{\mathrm{I}}^{0}\right];k,k^{\prime},k^{\prime\prime}\right)$
$\displaystyle=\dfrac{1}{2}\dfrac{d^{2}K_{\mathrm{I}}}{da^{2}}\\!(a_{0})-\dfrac{1}{a_{0}}\dfrac{dK_{\mathrm{I}}}{da}(a_{0})\dfrac{|{k+k^{\prime}}|}{4}+\dfrac{1}{a_{0}^{2}}K_{\mathrm{I}}^{0}(a_{0})\dfrac{|{k}|+|{k^{\prime}}|}{4}$
(82)
$\displaystyle+\dfrac{1}{8a_{0}^{2}}K_{\mathrm{I}}^{0}(a_{0})|{k}||{k^{\prime}}|+\dfrac{1}{16a_{0}^{2}}K_{\mathrm{I}}^{0}(a_{0})\left(|{k+k^{\prime}}|(|{k}|+|{k^{\prime}}|)-k^{2}-{k^{\prime}}^{2}\right)\left(1-\mathrm{sign}(kk^{\prime})\right)$
$\displaystyle=\dfrac{1}{2}\dfrac{d^{2}K_{\mathrm{I}}}{da^{2}}\\!(a_{0})-\dfrac{1}{a_{0}}\dfrac{dK_{\mathrm{I}}}{da}(a_{0})\dfrac{|{k+k^{\prime}}|}{4}+\dfrac{1}{a_{0}^{2}}K_{\mathrm{I}}^{0}(a_{0})\dfrac{|{k}|+|{k^{\prime}}|}{4}+\dfrac{1}{16a_{0}^{2}}K_{\mathrm{I}}^{0}(a_{0})R(k,k^{\prime}),$

from which one finds back Eq. (20).

The general expression (B.11-12) of the second-order variations of stress
intensity factor for a perturbed crack front given prove rather complex, and
its scope of application seems limited. Indeed, from a theoretical point of
view, the influence of a non-uniformity of macroscopic loading on mode I crack
propagation is already of second-order within the linear expansion of Gao and
Rice (1987b) (see Eq. (2)). Moreover, the triple sum over the Fourier
coefficients related to $K_{\mathrm{I}}^{0}$, its derivatives with crack
advance $(\partial^{n}K_{\mathrm{I}}^{0}/\partial a^{n})$, and the
perturbation shape $\phi$ makes it unfit for numerical calculations, as
estimating second-order contributions would take a computation time in
$\mathcal{O}(n^{3}\ln n)$, where $n$ is the number of discretization points
along the crack front. This is already larger than more conventional numerical
approaches, like e.g. the boundary element method, that can simulate crack
propagation _at all orders_ in “only” $\mathcal{O}(n^{2}\ln n)$.

However, we did not perform here the rigorous calculation for a general
$K_{\mathrm{I}}^{0}$ for the sake of completeness only. Indeed, for shear
ruptures, the mode II and mode III stress intensity factors applied to a
circular crack depends on $\theta$, even when the applied stress is
axisymmetric (Gao, 1988). Propagating internal shear cracks often adopt quasi-
elliptical front shapes (Sáez et al., 2022) that may be described with a
perturbative approach to LEFM. We provide thus in this Appendix the basic
steps needed to derive an equivalent second-order theory for shear cracks.
Future work will be devoted to the construction of this model and its
application to the propagation of e.g. fluid-driven shear ruptures.

##  Appendix C Validation of the general formula (80) on elliptical cracks
loaded by non-axisymmetric loadings

This Appendix aims at validating the general expression of the second-order
variations of stress intensity factor along a perturbed crack front obtained
in Eqs. (B.11-12). We provided in Section 3.1 a validation of Eq. (20) or
equivalently (82) by comparing our second-order predictions to the second-
order expansions of the exact stress intensity factor distribution along
elliptical crack fronts loaded by a uniform remote stress $\sigma$.
Fortunately, analytical predictions for the SIF distribution along elliptical
crack fronts are available in the literature for more complex loading
situations (Shah and Kobayashi, 1971; Kassir, 1975; Atroshchenko et al.,
2009). Here, a comparison between our model’s predictions and second-order
expansions of exact analytical solutions is provided for three non-
axisymmetric loading cases.

As in Section 3.1, we consider a crack whose front $\mathcal{F}$ describes an
ellipsis centered in $O$, of short semi-axis $c$ aligned with
$\mathbf{e_{z}}$, and long semi-axis $b$ aligned with $\mathbf{e_{x}}$. The
crack front $\mathcal{F}$ follows the contour parametrized by the variable
$\varphi$: $\left(z=b\cos\varphi,x=c\sin\varphi\right)$. Its position
$a(\theta)$ in polar coordinates is given by Eq. (26), so that:

${\color[rgb]{0,0,0}a(\theta)\cos\theta=b\cos\varphi\text{, and
}a(\theta)\sin\theta=c\sin\varphi.}$ (83)

We follow here the reasoning of Gao and Rice (1987b), and take as small
expansion parameter $\epsilon=1-c/b$, and a reference radius $a_{0}$ equal to
the front extension at the point $\theta_{1}$ of SIF evaluation (i.e.
$a_{0}=a(\theta_{1})$). In that case, the perturbation $\delta a$ reads:

$\delta a(\theta)=a_{0}\left(1+\epsilon
A_{1}(\theta)+\epsilon^{2}A_{2}(\theta)\right),$ (84)

where

$\displaystyle
A_{1}(\theta)=-\dfrac{1}{2}\cos(2\theta_{1})+\dfrac{1}{2}\cos(2\theta),$ (85)
and $\displaystyle
A_{2}(\theta)=\dfrac{1}{8}-\dfrac{1}{4}\cos(2\theta_{1})-\dfrac{1}{16}\cos(4\theta_{1})+\left(\dfrac{1}{4}-\dfrac{1}{4}\cos(2\theta_{1})\right)\cos(2\theta)+\dfrac{3}{16}\cos(4\theta).$

Next, we introduce $k_{\mathrm{I}}$, the crack face weight function of a
circular crack. $k_{\mathrm{I}}(a,\theta;r,\theta^{\prime})$ corresponds to
the mode I SIF generated at the point $\theta$ along the circular front of
radius $a$ by a pair of unitary symmetric tensile forces acting on the point
interior to the crack indexed by $(r,\theta^{\prime})$. It reads (Tada et al.,
2000):

$k_{\mathrm{I}}(a,\theta;r,\theta^{\prime})=\dfrac{\sqrt{a^{2}-r^{2}}}{\pi^{3/2}a^{1/2}\left[a^{2}+r^{2}-2ar\cos(\theta-\theta^{\prime})\right]}.$
(86)

Considering that the crack is loaded by a distribution of tensile stress

$\sigma(r,\theta)=\sum_{k}\hat{\sigma}_{k}(r)e^{ik\theta},$

acting on the failure plane $y=0$ of the intact body, the mode I SIF generated
along the front of a circular crack of radius $a$ reads:

$\displaystyle K_{\mathrm{I}}^{0}(a,\theta)$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{a}\sigma(r,\theta^{\prime})k_{\mathrm{I}}(a,\theta;r,\theta^{\prime})rdrd\theta^{\prime}$
(87) $\displaystyle=\sum_{k}\left[\dfrac{2}{\pi}\sqrt{\pi
a}\int_{0}^{1}\dfrac{u^{|{k}|+1}}{\sqrt{1-u^{2}}}\hat{\sigma}_{k}(au)du\right]e^{ik\theta}.$

We have at our disposal all the basic ingredients to use our second-order
model of Eqs. (76) and (80).

#### Linear stress field dependent on $z$

We first consider the case where the elliptical crack is loaded by a linear
stress field dependent on $z$:

$\sigma(z,x)=\sigma_{0}\dfrac{z}{b}\Rightarrow\sigma(r,\theta)=\sigma_{0}\dfrac{r}{b}\cos(\theta).$
(88)

Using Eq. (87) the associated SIF distribution along the front of a circular
crack of radius $a_{0}$ reads:

$K_{\mathrm{I}}^{0}(a,\theta)=\dfrac{4}{3\sqrt{\pi}}\dfrac{a^{3/2}}{b}\sigma_{0}\cos(\theta).$
(89)

The exact SIF distribution along the elliptical crack front has been
calculated by (Shah and Kobayashi, 1971) and (Atroshchenko et al., 2009). It
reads:

$\displaystyle K_{\mathrm{I}}^{e}(\varphi)=-\sqrt{\pi
c}\sigma_{0}\cos(\varphi)\dfrac{\kappa^{2}}{(1-2\kappa^{2})E(\kappa)-{\kappa^{\prime}}^{2}K(\kappa)}\left((c/b)^{4}\cos^{2}(\varphi)+\sin^{2}(\varphi)\right)^{1/4}$
(90) $\displaystyle\Rightarrow$ $\displaystyle
K_{\mathrm{I}}^{e}(\theta)=-\dfrac{\sqrt{\pi}c^{3/2}}{b}\sigma_{0}\cos(\theta)\dfrac{\kappa^{2}}{(1-2\kappa^{2})E(\kappa)-{\kappa^{\prime}}^{2}K(\kappa)}\dfrac{\left((c/b)^{4}\cos^{2}(\theta)+\sin^{2}(\theta)\right)^{1/4}}{\left((c/b)^{2}\cos^{2}(\theta)+\sin^{2}(\theta)\right)^{3/4}}.$

where $\kappa=\left(1-c^{2}/b^{2}\right)^{1/2}$ and $\kappa^{\prime}=c/b$.

One may then use Eqs. (89) and (90) to expand
$K_{\mathrm{I}}^{e}/K_{\mathrm{I}}^{0}$ to second order in $\epsilon$ at
$\epsilon=0$ and $\theta=\theta_{1}$. One finds:

$\displaystyle K_{\mathrm{I}}^{e}(A_{1})/K_{\mathrm{I}}^{0}(a_{0},\theta_{1})$
$\displaystyle=-\dfrac{3\pi}{4}\dfrac{\kappa^{2}}{(1-2\kappa^{2})E(\kappa)-{\kappa^{\prime}}^{2}K(\kappa)}\left((c/b)^{4}\cos^{2}(\theta_{1})+\sin^{2}(\theta_{1})\right)^{1/4}$
(91)
$\displaystyle=1+\epsilon\left[\dfrac{1}{4}-\dfrac{1}{2}\cos(2\theta_{1})\right]+\epsilon^{2}\left[\dfrac{5}{16}-\dfrac{3}{8}\cos(2\theta_{1})-\dfrac{3}{16}\cos(4\theta_{1})\right]+\mathcal{O}(\epsilon^{3}).$

Our aim here is to retrieve Eq. (91) from our general second-order theory of
Eqs. (76) and (80). The second-order estimate $K_{\mathrm{I}}^{p}$ writes as:

$K_{\mathrm{I}}^{p}(\theta_{1})=K_{\mathrm{I}}^{0}(a_{0},\theta_{1})+\epsilon
K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[A_{1}];\theta_{1})+\epsilon^{2}\left[K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[A_{2}];\theta_{1})+K_{\mathrm{I}}^{2}(a_{0},[K_{\mathrm{I}}^{0}],[A_{1}];\theta_{1})\right].$
(92)

Combining Eqs. (76) and (85), one has

$\displaystyle
K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[A_{1}];\theta_{1})=K_{\mathrm{I}}^{0}(a_{0},\theta_{1})\times\left[\dfrac{1}{4}-\dfrac{1}{2}\cos(2\theta_{1})\right],$
(93) and $\displaystyle
K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[A_{2}];\theta_{1})=K_{\mathrm{I}}^{0}(a_{0},\theta_{1})\times\left[\dfrac{5}{32}-\dfrac{3}{16}\cos(2\theta_{1})-\dfrac{1}{4}\cos(4\theta_{1})\right].$

Similarly, the combination of Eqs. (80) and (85) yields:

$\displaystyle
K_{\mathrm{I}}^{2}(a_{0},[K_{\mathrm{I}}^{0}],[A_{1}];\theta_{1})=K_{\mathrm{I}}^{0}(a_{0},\theta_{1})\times\left[\dfrac{5}{32}-\dfrac{3}{16}\cos(2\theta_{1})+\dfrac{1}{16}\cos(4\theta_{1})\right].$
(94)

We observe that the two expansions coincide at $\theta=\theta_{1}$.

#### Linear stress field dependent on $x$

We consider next the case where the elliptical crack is loaded by a linear
stress field dependent on $x$:

$\sigma(z,x)=\sigma_{0}\dfrac{x}{c}\Rightarrow\sigma(r,\theta)=\sigma_{0}\dfrac{r}{c}\sin(\theta).$
(95)

As the procedure is similar to the previous case, we should be brief here. One
has:

$K_{\mathrm{I}}^{0}(a,\theta)=\dfrac{4}{3\sqrt{\pi}}\dfrac{a^{3/2}}{c}\sigma_{0}\sin(\theta).$
(96)

The exact SIF distribution along the elliptical crack front has been obtained
by (Shah and Kobayashi, 1971) and (Atroshchenko et al., 2009). It reads:

$\displaystyle K_{\mathrm{I}}^{e}(\varphi)=\sqrt{\pi
c}\sigma_{0}\sin(\varphi)\dfrac{\kappa^{2}}{(1+\kappa^{2})E(\kappa)-{\kappa^{\prime}}^{2}K(\kappa)}\left((c/b)^{4}\cos^{2}(\varphi)+\sin^{2}(\varphi)\right)^{1/4}$
(97) $\displaystyle\Rightarrow$ $\displaystyle
K_{\mathrm{I}}^{e}(\theta)=\sqrt{\pi
c}\sigma_{0}\sin(\theta)\dfrac{\kappa^{2}}{(1+\kappa^{2})E(\kappa)-{\kappa^{\prime}}^{2}K(\kappa)}\dfrac{\left((c/b)^{4}\cos^{2}(\theta)+\sin^{2}(\theta)\right)^{1/4}}{\left((c/b)^{2}\cos^{2}(\theta)+\sin^{2}(\theta)\right)^{3/4}}.$

The expansion of $K_{\mathrm{I}}^{e}/K_{\mathrm{I}}^{0}$ at
$\theta=\theta_{1}$ to second order in $\epsilon$ at $\epsilon=0$ and
$\theta=\theta_{1}$ yields:

$\displaystyle K_{\mathrm{I}}^{e}(A_{1})/K_{\mathrm{I}}^{0}(a_{0},\theta_{1})$
$\displaystyle=\dfrac{3\pi}{4}\dfrac{\kappa^{2}}{(1+\kappa^{2})E(\kappa)-{\kappa^{\prime}}^{2}K(\kappa)}\left((c/b)^{4}\cos^{2}(\theta_{1})+\sin^{2}(\theta_{1})\right)^{1/4}$
(98)
$\displaystyle=1+\epsilon\left[-\dfrac{1}{4}-\dfrac{1}{2}\cos(2\theta_{1})\right]+\epsilon^{2}\left[\dfrac{1}{16}-\dfrac{1}{8}\cos(2\theta_{1})-\dfrac{3}{16}\cos(4\theta_{1})\right]+\mathcal{O}(\epsilon^{3}).$

Our goal is to retrieve again Eq. (98) from Eq. (92). One has:

$\displaystyle
K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[A_{1}];\theta_{1})=K_{\mathrm{I}}^{0}(a_{0},\theta_{1})\times\left[-\dfrac{1}{4}-\dfrac{1}{2}\cos(2\theta_{1})\right],$
(99) and $\displaystyle
K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[A_{2}];\theta_{1})=K_{\mathrm{I}}^{0}(a_{0},\theta_{1})\times\left[-\dfrac{3}{32}-\dfrac{5}{16}\cos(2\theta_{1})-\dfrac{1}{4}\cos(4\theta_{1})\right],$

and:

$\displaystyle
K_{\mathrm{I}}^{2}(a_{0},[K_{\mathrm{I}}^{0}],[A_{1}];\theta_{1})=K_{\mathrm{I}}^{0}(a_{0},\theta_{1})\times\left[\dfrac{5}{32}+\dfrac{3}{16}\cos(2\theta_{1})+\dfrac{1}{16}\cos(4\theta_{1})\right].$
(100)

The two expansions are found equal once again.

#### Quadratic stress field odd in $z$ and $x$

We conclude this Appendix with a last validation case, in which the stress
field is quadratic and odd in both $z$ and $x$:

$\sigma(z,x)=\sigma_{0}\dfrac{xy}{bc}\Rightarrow\sigma(r,\theta)=\sigma_{0}\dfrac{r^{2}}{2bc}\sin(2\theta).$
(101)

For which one has:

$K_{\mathrm{I}}^{0}(a,\theta)=\dfrac{8}{15\sqrt{\pi}}\dfrac{a^{5/2}}{bc}\sigma_{0}\sin(2\theta).$
(102)

The exact SIF distribution along the elliptical crack front has been obtained
by (Kassir, 1975) and (Atroshchenko et al., 2009). It reads:

$\displaystyle K_{\mathrm{I}}^{e}(\varphi)=\sqrt{\pi
c}\sigma_{0}\sin(2\varphi)\dfrac{\kappa^{4}}{2\left[2({\kappa^{\prime}}^{4}-{\kappa^{\prime}}^{2}+1)E(\kappa)-({\kappa^{\prime}}^{4}+{\kappa^{\prime}}^{2})K(\kappa)\right]}\left((c/b)^{4}\cos^{2}(\varphi)+\sin^{2}(\varphi)\right)^{1/4}$
(103) $\displaystyle\Rightarrow$ $\displaystyle
K_{\mathrm{I}}^{e}(\varphi)=\dfrac{\sqrt{\pi}c^{3/2}}{2b}\sigma_{0}\sin(2\theta)\dfrac{\kappa^{4}}{2({\kappa^{\prime}}^{4}-{\kappa^{\prime}}^{2}+1)E(\kappa)-({\kappa^{\prime}}^{4}+{\kappa^{\prime}}^{2})K(\kappa)}\dfrac{\left((c/b)^{4}\cos^{2}(\theta)+\sin^{2}(\theta)\right)^{1/4}}{\left((c/b)^{2}\cos^{2}(\theta)+\sin^{2}(\theta)\right)^{5/4}}.$

The expansion of $K_{\mathrm{I}}^{e}/K_{\mathrm{I}}^{0}$ at
$\theta=\theta_{1}$ to second order in $\epsilon$ at $\epsilon=0$ and
$\theta=\theta_{1}$ then yields:

$\displaystyle K_{\mathrm{I}}^{e}(A_{1})/K_{\mathrm{I}}^{0}(a_{0},\theta_{1})$
$\displaystyle=\dfrac{15\pi}{16}\dfrac{\kappa^{4}}{2({\kappa^{\prime}}^{4}-{\kappa^{\prime}}^{2}+1)E(\kappa)-({\kappa^{\prime}}^{4}+{\kappa^{\prime}}^{2})K(\kappa)}\left((c/b)^{4}\cos^{2}(\theta_{1})+\sin^{2}(\theta_{1})\right)^{1/4}$
(104)
$\displaystyle=1+\epsilon\left[-\dfrac{1}{2}\cos(2\theta_{1})\right]+\epsilon^{2}\left[\dfrac{3}{32}-\dfrac{1}{4}\cos(2\theta_{1})-\dfrac{3}{16}\cos(4\theta_{1})\right]+\mathcal{O}(\epsilon^{3}).$

We need now to compute the three different terms of our second-order equation
(92). The terms arising from the linearization of $K_{\mathrm{I}}$ yields:

$\displaystyle
K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[A_{1}];\theta_{1})=K_{\mathrm{I}}^{0}(a_{0},\theta_{1})\times\left[-\dfrac{1}{2}\cos(2\theta_{1})\right],$
(105) and $\displaystyle
K_{\mathrm{I}}^{1}(a_{0},[K_{\mathrm{I}}^{0}],[A_{2}];\theta_{1})=K_{\mathrm{I}}^{0}(a_{0},\theta_{1})\times\left[-\dfrac{1}{16}-\dfrac{1}{4}\cos(2\theta_{1})-\dfrac{1}{4}\cos(4\theta_{1})\right],$

while the second-order contribution read:

$\displaystyle
K_{\mathrm{I}}^{2}(a_{0},[K_{\mathrm{I}}^{0}],[A_{1}];\theta_{1})=K_{\mathrm{I}}^{0}(a_{0},\theta_{1})\times\left[\dfrac{5}{32}+\dfrac{1}{16}\cos(4\theta_{1})\right].$
(106)

The two expansions are found equal once more, although the shape of the stress
field is more complex.

##  Appendix D Ensemble averages of crack front propagating in weakly
disordered materials