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

Mathias Lebihain; Manish Vasoya; V\'eronique Lazarus

arXiv:2302.12083·cond-mat.mtrl-sci·June 27, 2023

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

Mathias Lebihain, Manish Vasoya, V\'eronique Lazarus

PDF

Open Access

TL;DR

This paper extends a theoretical model to predict stress intensity factors along slightly perturbed crack fronts in brittle materials, validating it with simulations and applying it to understand size-dependent toughening effects due to heterogeneities.

Contribution

It develops a second-order theory for crack front perturbations, validated against analytical solutions and simulations, and introduces a homogenization framework for toughness in disordered media.

Findings

01

Toughness heterogeneities weaken small or comparable cracks.

02

Heterogeneities reinforce larger cracks, leading to size-dependent effects.

03

The theory predicts an R-curve behavior at the macroscale.

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 (1987), and predict the distribution of the mode I stress intensity factor $K_{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…

Equations321

δ a (θ) = A ϕ (θ),

δ a (θ) = A ϕ (θ),

δ K_{I} (θ_{1}) = \frac{\partial K _{I}^{0}}{\partial a} (a_{0}, θ_{1}) δ a (θ_{1}) + PV \int_{0}^{2 π} Z^{0} (a_{0}; θ_{1}, θ) K_{I}^{0} (a_{0}, θ) [δ a (θ) - δ a (θ_{1})] a_{0} d θ + O (A^{2}),

δ K_{I} (θ_{1}) = \frac{\partial K _{I}^{0}}{\partial a} (a_{0}, θ_{1}) δ a (θ_{1}) + PV \int_{0}^{2 π} Z^{0} (a_{0}; θ_{1}, θ) K_{I}^{0} (a_{0}, θ) [δ a (θ) - δ a (θ_{1})] a_{0} d θ + O (A^{2}),

Z^{0} (a_{0}; θ_{1}, θ_{2}) = \frac{1}{2 π} \frac{W ^{0} ( a _{0} ; θ _{1} , θ _{2} )}{D ( M _{1}^{0} , M _{2}^{0} ) ^{2}},

Z^{0} (a_{0}; θ_{1}, θ_{2}) = \frac{1}{2 π} \frac{W ^{0} ( a _{0} ; θ _{1} , θ _{2} )}{D ( M _{1}^{0} , M _{2}^{0} ) ^{2}},

Z^{0} (a_{0}; θ_{1}, θ_{2}) = \frac{1}{8 π a _{0}^{2} sin ^{2} [ ( θ _{2} - θ _{1} ) /2 ]} .

Z^{0} (a_{0}; θ_{1}, θ_{2}) = \frac{1}{8 π a _{0}^{2} sin ^{2} [ ( θ _{2} - θ _{1} ) /2 ]} .

\frac{δ W ( θ _{1} , θ _{2} )}{D ( M _{1} , M _{2} ) ^{2}} = \frac{1}{2 π} PV \int_{0}^{2 π} \frac{W ^{0} ( a _{0} ; θ _{1} , θ )}{D ( M _{1} , M ^{0} ) ^{2}} \frac{W ^{0} ( a _{0} ; θ _{2} , θ )}{D ( M _{2} , M ^{0} ) ^{2}} δ a (θ) a_{0} d θ + O (A^{2}),

\frac{δ W ( θ _{1} , θ _{2} )}{D ( M _{1} , M _{2} ) ^{2}} = \frac{1}{2 π} PV \int_{0}^{2 π} \frac{W ^{0} ( a _{0} ; θ _{1} , θ )}{D ( M _{1} , M ^{0} ) ^{2}} \frac{W ^{0} ( a _{0} ; θ _{2} , θ )}{D ( M _{2} , M ^{0} ) ^{2}} δ a (θ) a_{0} d θ + O (A^{2}),

δ a_{*} (θ_{1}) = δ a (θ_{1}), and δ a_{*} (θ_{2}) = δ a (θ_{2}) .

δ a_{*} (θ_{1}) = δ a (θ_{1}), and δ a_{*} (θ_{2}) = δ a (θ_{2}) .

\frac{δ W ( θ _{1} , θ _{2} )}{D ( M _{1}^{0} , M _{2}^{0} ) ^{2}} = \frac{δ W ^{*} ( θ _{1} , θ _{2} )}{D ( M _{1}^{0} , M _{2}^{0} ) ^{2}} + \frac{1}{2 π} PV \int_{0}^{2 π} \frac{W ^{0} ( a _{0} ; θ _{1} , θ )}{D ( M _{1}^{0} , M ^{0} ) ^{2}} \frac{W ^{0} ( a _{0} ; θ _{2} , θ )}{D ( M _{2}^{0} , M ^{0} ) ^{2}} [δ a (θ) - δ a_{*} (θ)] a_{0} d θ + O (A^{2}),

\frac{δ W ( θ _{1} , θ _{2} )}{D ( M _{1}^{0} , M _{2}^{0} ) ^{2}} = \frac{δ W ^{*} ( θ _{1} , θ _{2} )}{D ( M _{1}^{0} , M _{2}^{0} ) ^{2}} + \frac{1}{2 π} PV \int_{0}^{2 π} \frac{W ^{0} ( a _{0} ; θ _{1} , θ )}{D ( M _{1}^{0} , M ^{0} ) ^{2}} \frac{W ^{0} ( a _{0} ; θ _{2} , θ )}{D ( M _{2}^{0} , M ^{0} ) ^{2}} [δ a (θ) - δ a_{*} (θ)] a_{0} d θ + O (A^{2}),

δ Z (θ_{1}, θ_{2})

δ Z (θ_{1}, θ_{2})

= Z^{0} (a_{0}; θ_{1}, θ_{2}) [(\frac{D ( M _{1} , M _{2} )}{D ( M _{1}^{0} , M _{2}^{0} )})^{- 2} - 1] + \frac{δ W ^{*} ( θ _{1} , θ _{2} )}{2 π D ( M _{1}^{0} , M _{2}^{0} ) ^{2}}

+ PV \int_{0}^{2 π} Z^{0} (a_{0}; θ_{1}, θ) Z^{0} (a_{0}; θ_{2}, θ) [δ a (θ) - δ a_{*} (θ)] a_{0} d θ + O (A^{2}),

Z^{0} (a_{0}; θ_{1}, θ_{2}) [(\frac{D ( M _{1} , M _{2} )}{D ( M _{1}^{0} , M _{2}^{0} )})^{- 2} - 1] = - \frac{δ a ( θ _{1} ) + δ a ( θ _{2} )}{8 π a _{0}^{3} sin ^{2} [ ( θ _{2} - θ _{1} ) /2 ]} + O (A^{2}) .

Z^{0} (a_{0}; θ_{1}, θ_{2}) [(\frac{D ( M _{1} , M _{2} )}{D ( M _{1}^{0} , M _{2}^{0} )})^{- 2} - 1] = - \frac{δ a ( θ _{1} ) + δ a ( θ _{2} )}{8 π a _{0}^{3} sin ^{2} [ ( θ _{2} - θ _{1} ) /2 ]} + O (A^{2}) .

δ a_{*} (θ) = Δ r_{*} sin (θ - θ_{*}) + (a_{0} + Δ a_{*})^{2} - Δ r_{*}^{2} cos^{2} (θ - θ_{*}) - a_{0} .

δ a_{*} (θ) = Δ r_{*} sin (θ - θ_{*}) + (a_{0} + Δ a_{*})^{2} - Δ r_{*}^{2} cos^{2} (θ - θ_{*}) - a_{0} .

θ_{*} = \frac{θ _{1} + θ _{2}}{2} .

θ_{*} = \frac{θ _{1} + θ _{2}}{2} .

Δ r_{*} = \frac{δ a ( θ _{2} ) - δ a ( θ _{1} )}{2 sin [ ( θ _{2} - θ _{1} ) /2 ]} + O (A^{2}), and Δ a_{*} = \frac{δ a ( θ _{1} ) + δ a ( θ _{2} )}{2} + O (A^{2}) .

Δ r_{*} = \frac{δ a ( θ _{2} ) - δ a ( θ _{1} )}{2 sin [ ( θ _{2} - θ _{1} ) /2 ]} + O (A^{2}), and Δ a_{*} = \frac{δ a ( θ _{1} ) + δ a ( θ _{2} )}{2} + O (A^{2}) .

δ a_{*} (θ) = \frac{δ a ( θ _{1} ) + δ a ( θ _{2} )}{2} + \frac{δ a ( θ _{2} ) - δ a ( θ _{1} )}{2} \cdot \frac{sin [ θ - ( θ _{1} + θ _{2} ) /2 ]}{sin [ ( θ _{2} - θ _{1} ) /2 ]} + O (A^{2})

δ a_{*} (θ) = \frac{δ a ( θ _{1} ) + δ a ( θ _{2} )}{2} + \frac{δ a ( θ _{2} ) - δ a ( θ _{1} )}{2} \cdot \frac{sin [ θ - ( θ _{1} + θ _{2} ) /2 ]}{sin [ ( θ _{2} - θ _{1} ) /2 ]} + O (A^{2})

= δ a (θ_{1}) + \frac{δ a ( θ _{2} ) - δ a ( θ _{1} )}{2} [\frac{cos [ ( θ _{2} - θ _{1} ) /2 ]}{sin [ ( θ _{2} - θ _{1} ) /2 ]} sin (θ - θ_{1}) + 2 sin^{2} [(θ - θ_{1}) /2]] + O (A^{2})

= δ a (θ_{2}) + \frac{δ a ( θ _{2} ) - δ a ( θ _{1} )}{2} [\frac{cos [ ( θ _{2} - θ _{1} ) /2 ]}{sin [ ( θ _{2} - θ _{1} ) /2 ]} sin (θ - θ_{2}) - 2 sin^{2} [(θ - θ_{2}) /2]] + O (A^{2}) .

\frac{1}{sin ^{2} [( θ - θ _{1} ) /2 ] sin ^{2} [( θ - θ _{2} ) /2 ]} =

\frac{1}{sin ^{2} [( θ - θ _{1} ) /2 ] sin ^{2} [( θ - θ _{2} ) /2 ]} =

+ \frac{2 cos [( θ _{2} - θ _{1} ) /2 ]}{sin ^{3} [( θ _{2} - θ _{1} ) /2 ]} [\frac{cos ^{3} [( θ - θ _{1} ) /2 ]}{sin [( θ - θ _{1} ) /2 ]} - \frac{cos ^{3} [( θ - θ _{2} ) /2 ]}{sin [( θ - θ _{2} ) /2 ]}]

+ \frac{2}{sin ^{2} [ ( θ _{2} - θ _{1} ) /2 ]} [cos^{2} [(θ - θ_{1}) /2] + cos^{2} [(θ - θ_{2}) /2]]

\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}

\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}

{\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}}

{\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}}

K_{I} (a_{0}, [(A + δ A) ϕ]; θ_{1}) - K_{I} (a_{0}, [A ϕ]; θ_{1})

K_{I} (a_{0}, [(A + δ A) ϕ]; θ_{1}) - K_{I} (a_{0}, [A ϕ]; θ_{1})

= {K_{I}^{1} (a_{0}, [K_{I}^{0}], [ϕ]; θ_{1}) + 2 A K_{I}^{2} (a_{0}, [K_{I}^{0}], [ϕ]; θ_{1}) + O (A^{2})} δ A + O (δ A^{2})

K_{I} (a_{0}, [(A + δ A) ϕ]; θ_{1}) - K_{I} (a_{0}, [A ϕ]; θ_{1}) = K_{I} (a_{0} + δ A ϕ (θ_{1}), [A ϕ]; a_{0}; θ_{1})

K_{I} (a_{0}, [(A + δ A) ϕ]; θ_{1}) - K_{I} (a_{0}, [A ϕ]; θ_{1}) = K_{I} (a_{0} + δ A ϕ (θ_{1}), [A ϕ]; a_{0}; θ_{1})

+ PV \int_{0}^{2 π} Z (a_{0} + δ A ϕ (θ_{1}), [A ϕ]; θ_{1}, θ) K_{I} (a_{0} + δ A ϕ (θ_{1}), [A ϕ]; θ) [δ A ϕ (θ) - δ A ϕ (θ_{1})] d s (θ)

- K_{I} (a_{0}, [A ϕ]; θ_{1}) + O (δ A^{2}),

K_{I}^{1} (a_{0}, [K_{I}^{0}], [ϕ]; θ_{1}) = \frac{\partial K _{I}^{0}}{\partial a} (a_{0}, θ_{1}) ϕ (θ_{1}) + PV \int_{0}^{2 π} Z^{0} (a_{0}; θ_{1}, θ) K_{I}^{0} (a_{0}, θ) [ϕ (θ) - ϕ (θ_{1})] a_{0} d θ,

K_{I}^{1} (a_{0}, [K_{I}^{0}], [ϕ]; θ_{1}) = \frac{\partial K _{I}^{0}}{\partial a} (a_{0}, θ_{1}) ϕ (θ_{1}) + PV \int_{0}^{2 π} Z^{0} (a_{0}; θ_{1}, θ) K_{I}^{0} (a_{0}, θ) [ϕ (θ) - ϕ (θ_{1})] a_{0} d θ,

K_{I}^{2} (a_{0}, [K_{I}^{0}], [ϕ]; θ_{1}) =

K_{I}^{2} (a_{0}, [K_{I}^{0}], [ϕ]; θ_{1}) =

+ \frac{1}{2} PV \int_{0}^{2 π} Z^{0} (a_{0}; θ_{1}, θ) K_{I}^{0} (a_{0}, θ) [ϕ (θ) - ϕ (θ_{1})]^{2} a_{0} d θ

+ \frac{1}{2} PV \int_{0}^{2 π} Z^{1} (a_{0}, ϕ; θ_{1}, θ) K_{I}^{0} (a_{0}, θ) [ϕ (θ) - ϕ (θ_{1})] a_{0} d θ

+ \frac{1}{2} PV \int_{0}^{2 π} \int_{0}^{2 π} Z^{0} (a_{0}; θ_{1}, θ) Z^{0} (a_{0}; θ, θ^{'}) K_{I}^{0} (a_{0}, θ^{'}) [ϕ (θ) - ϕ (θ_{1})] [ϕ (θ^{'}) - ϕ (θ)] a_{0}^{2} d θ d θ^{'} .

ϕ (θ) = k = - \infty \sum + \infty \hat{ϕ}_{k} e^{ik θ}, where \hat{ϕ}_{k} = \frac{1}{2 π} \int_{0}^{2 π} ϕ (θ) e^{- ik θ} d θ .

ϕ (θ) = k = - \infty \sum + \infty \hat{ϕ}_{k} e^{ik θ}, where \hat{ϕ}_{k} = \frac{1}{2 π} \int_{0}^{2 π} ϕ (θ) e^{- ik θ} d θ .

K_{I}^{1} (a_{0}, [K_{I}^{0}], [ϕ]; θ_{1}) = k = - \infty \sum + \infty [\frac{d K _{I}^{0}}{d a} (a_{0}) - \frac{∣ k ∣}{2 a _{0}} K_{I}^{0} (a_{0})] \hat{ϕ}_{k} e^{ik θ_{1}},

K_{I}^{1} (a_{0}, [K_{I}^{0}], [ϕ]; θ_{1}) = k = - \infty \sum + \infty [\frac{d K _{I}^{0}}{d a} (a_{0}) - \frac{∣ k ∣}{2 a _{0}} K_{I}^{0} (a_{0})] \hat{ϕ}_{k} e^{ik θ_{1}},

K_{I}^{2} (a_{0}, [K_{I}^{0}], [ϕ]; θ_{1}) = k = - \infty \sum + \infty k^{'} = - \infty \sum + \infty [\frac{1}{2} \frac{d ^{2} K _{I}^{0}}{d a ^{2}} (a_{0}) - \frac{1}{4 a _{0}} \frac{d K _{I}^{0}}{d a} (a_{0}) ∣ k + k^{'} ∣ + \frac{R ( k , k ^{'} )}{16 a _{0}^{2}} K_{I}^{0} (a_{0}) + \frac{∣ k ∣ + ∣ k ^{'} ∣}{4 a _{0}^{2}} K_{I}^{0} (a_{0})] \hat{ϕ}_{k} \hat{ϕ}_{k^{'}} e^{i (k + k^{'}) θ_{1}},

K_{I}^{2} (a_{0}, [K_{I}^{0}], [ϕ]; θ_{1}) = k = - \infty \sum + \infty k^{'} = - \infty \sum + \infty [\frac{1}{2} \frac{d ^{2} K _{I}^{0}}{d a ^{2}} (a_{0}) - \frac{1}{4 a _{0}} \frac{d K _{I}^{0}}{d a} (a_{0}) ∣ k + k^{'} ∣ + \frac{R ( k , k ^{'} )}{16 a _{0}^{2}} K_{I}^{0} (a_{0}) + \frac{∣ k ∣ + ∣ k ^{'} ∣}{4 a _{0}^{2}} K_{I}^{0} (a_{0})] \hat{ϕ}_{k} \hat{ϕ}_{k^{'}} e^{i (k + k^{'}) θ_{1}},

K_{I} (a_{0}, [δ a]; θ_{1}) = K_{I}^{0} (a_{0}) 1 + k \sum L (k) \frac{δ a _{k}}{a _{0}} e^{ik θ_{1}} + k, k^{'} \sum H (k, k^{'}) \frac{δ a _{k}}{a _{0}} \frac{δ a _{k^{'}}}{a _{0}} e^{i (k + k^{'}) θ_{1}} + O (A^{3}),

K_{I} (a_{0}, [δ a]; θ_{1}) = K_{I}^{0} (a_{0}) 1 + k \sum L (k) \frac{δ a _{k}}{a _{0}} e^{ik θ_{1}} + k, k^{'} \sum H (k, k^{'}) \frac{δ a _{k}}{a _{0}} \frac{δ a _{k^{'}}}{a _{0}} e^{i (k + k^{'}) θ_{1}} + O (A^{3}),

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsComposite Material Mechanics · Rock Mechanics and Modeling · Seismic Imaging and Inversion Techniques

Full text

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

{highlights}

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:

[TABLE]

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):

[TABLE]

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):

[TABLE]

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 kernel111Note 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:

[TABLE]

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$ .

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:

[TABLE]

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):

[TABLE]

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

[TABLE]

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:

[TABLE]

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 :

[TABLE]

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:

[TABLE]

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.

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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):

[TABLE]

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$ :

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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 [math] 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$ :

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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}$ :

[TABLE]

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

[TABLE]

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}}|$ :

[TABLE]

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:

[TABLE]

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$ :

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

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:

[TABLE]

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

[TABLE]

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$ .

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):

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

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):

[TABLE]

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:

[TABLE]

$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}$ :

[TABLE]

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):

[TABLE]

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

[TABLE]

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

[TABLE]

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):

[TABLE]

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:

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

so that:

[TABLE]

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:

[TABLE]

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).

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:

[TABLE]

while its standard deviation writes as:

[TABLE]

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

[TABLE]

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:

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

Using Eq. (2), one has:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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}$ :

[TABLE]

as well as the following integrals:

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

where

[TABLE]

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):

[TABLE]

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

[TABLE]

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:

[TABLE]

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$ :

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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$ :

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

and:

[TABLE]

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$ :

[TABLE]

For which one has:

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

while the second-order contribution read:

[TABLE]

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

This Appendix is dedicated to the statistical analysis of equation (61) defining the effective toughness of a random medium. In the following, we consider a statistical ensemble $\Omega$ containing all the possible realizations of the heterogeneous medium considered in Section 4. To each realization, we associate a specific real number $\omega$ . We note $p~{}:\omega\mapsto p\left(\omega\right)$ the probability function associated to the ensemble $\Omega$ . Thus, the probability that the variable $\omega^{\prime}$ lies in some neighborhood of $\omega$ of measure $d\omega$ is $p(\omega)d\omega$ . The mathematical expectation $E\left[u(x_{1},x_{2})\right]$ of any spatial observable $u:(x_{1},x_{2})\mapsto u(x_{1},x_{2})$ is defined as:

[TABLE]

Using these notations, the ensemble averaging of Eq. (61) yields:

[TABLE]

so that its determination requires the prior calculation of $\mathrm{E}\left[{\hat{A}^{1}_{0}(a_{0})}\right]$ , $\mathrm{E}\left[{\hat{A}^{2}_{0}(a_{0})}\right]$ , and $\mathrm{E}\left[{\hat{A}^{1}_{0}(a_{0})^{2}}\right]$ . We start here with $\hat{A}^{1}_{0}(a_{0})$ , for which the calculation is the simplest and illustrates well the reasoning behind the ulterior ones. For an arbitrary mode $k$ , $\mathrm{E}\left[{\hat{A}^{1}_{k}(a_{0})}\right]$ reads:

[TABLE]

where we used the expression of $\hat{A}^{1}_{k}(a_{0})$ given in Eq. (54). Even if only the value of $\mathrm{E}\left[{\hat{A}^{1}_{0}(a_{0})}\right]$ to estimate the average value of the effective toughness, our calculations show that:

[TABLE]

meaning that the first-order contribution to the crack front deformations are on average zero.

We proceed similarly for $A_{2}$ , and start by averaging Eq. (56):

[TABLE]

Let us first concentrate on:

[TABLE]

To calculate Eq. (112), it is convenient to switch to Cartesian coordinates. The pairs $(r,\theta)$ and $(r^{\prime},\theta^{\prime})$ are then equivalent to $(z=a_{0}\cos\theta,x=a_{0}\sin\theta)$ and $(z^{\prime}=a_{0}\cos\theta^{\prime},x^{\prime}=a_{0}\sin\theta^{\prime})$ . With these notations, one has:

[TABLE]

where we used the expressions of $\mathcal{F}_{z}$ and $\mathcal{F}_{x}$ given in Eq. (45), as well as the following identities:

[TABLE]

Injecting Eq. (113) into (112), one finds:

[TABLE]

where $\hat{F}_{k}$ is the $k$ -th Fourier coefficient of the $2\pi$ -periodic function:

[TABLE]

which relates to the correlation shape $\mathcal{F}$ of Eq. (45). We must now calculate the second average term in the sum of Eq.(111). It reads:

[TABLE]

The calculation of equation (116) proves more difficult than that of Eq. (113). In the following, use will be made of Fourier transforms. The definition adopted here for the single Fourier transform $\tilde{g}(p)$ of an arbitrary function $g(x)$ is:

[TABLE]

and the double $(z,x)$ Fourier transform $\tilde{\raisebox{2.0pt}[0.85pt]{$ \tilde{h} $}}(m,n)$ of an arbitrary function $h(z,x)$ reads:

[TABLE]

One then has:

[TABLE]

As shown by Favier et al. (2006b) in their Appendix B, one has:

[TABLE]

where $\delta$ is the Dirac function, and $\widetilde{\mathcal{F}}_{z}$ and $\widetilde{\mathcal{F}}_{x}$ are the Fourier transforms of the correlation functions $\mathcal{F}_{z}$ and $\mathcal{F}_{x}$ . Injecting Eq. (120) into (119), one finds:

[TABLE]

where we used the following identities:

[TABLE]

and:

[TABLE]

The combination of Eqs. (116) and (121) finally yields:

[TABLE]

where $\hat{G}_{k}$ is the $k$ -th Fourier coefficient of the $2\pi$ -periodic function:

[TABLE]

As a result, all $\mathrm{E}\left[{\hat{A}^{2}_{k}(a_{0})}\right]$ are zero, except for $k=0$ , for which one gets:

[TABLE]

This means that the average second-order contribution to the front deformations is a simple dilation of the reference circular front of radius $a_{0}$ by a quantity:

[TABLE]

We only need to calculate the third and last missing term of Eq. (108). Using Eq. (114), one has:

[TABLE]

so that the average of the effective toughness finally reads:

[TABLE]

One can additionally show that:

[TABLE]

Bibliography96

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Adda-Bedia et al. (2006) Adda-Bedia, M., Katzav, E., Vandembroucq, D., 2006. Second-order variation in elastic fields of a tensile planar crack with a curved front. Physical Review E 73, 035106. URL: https://link.aps.org/doi/10.1103/Phys Rev E.73.035106 , doi: 10.1103/Phys Rev E.73.035106 . · doi ↗
2Adda-Bedia and Mahadevan (2006) Adda-Bedia, M., Mahadevan, L., 2006. Crack-front instability in a confined elastic film. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 462, 3233–3251. URL: https://royalsocietypublishing.org/doi/full/10.1098/rspa.2006.1708 , doi: 10.1098/rspa.2006.1708 . · doi ↗
3Albertini et al. (2021) Albertini, G., Karrer, S., Grigoriu, M.D., Kammer, D.S., 2021. Stochastic properties of static friction. Journal of the Mechanics and Physics of Solids 147, 104242. URL: http://www.sciencedirect.com/science/article/pii/S 0022509620304555 , doi: 10.1016/j.jmps.2020.104242 . · doi ↗
4Argatov (2021) Argatov, I.I., 2021. Controlling the adhesive pull-off force via the change of contact geometry. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 379, 20200392. URL: https://royalsocietypublishing.org/doi/10.1098/rsta.2020.0392 , doi: 10.1098/rsta.2020.0392 . · doi ↗
5Atroshchenko et al. (2009) Atroshchenko, E., Potapenko, S., Glinka, G., 2009. Stress intensity factor for an embedded elliptical crack under arbitrary normal loading. International Journal of Fatigue 31, 1907–1910. URL: https://www.sciencedirect.com/science/article/pii/S 014211230800282 X , doi: 10.1016/j.ijfatigue.2008.12.004 . · doi ↗
6Barenblatt (1962) Barenblatt, G.I., 1962. The mathematical theory of equilibrium cracks in brittle fracture, in: Dryden, H.L., von Kármán, T., Kuerti, G., van den Dungen, F.H., Howarth, L. (Eds.), Advances in Applied Mechanics. Elsevier. volume 7, pp. 55–129.
7Barés et al. (2018) Barés, J., Dubois, A., Hattali, L., Dalmas, D., Bonamy, D., 2018. Aftershock sequences and seismic-like organization of acoustic events produced by a single propagating crack. Nature Communications 9, 1253. URL: https://www.nature.com/articles/s 41467-018-03559-4 , doi: 10.1038/s 41467-018-03559-4 . · doi ↗
8Bažant and Jirásek (1993) Bažant, Z.P., Jirásek, M., 1993. R-curve modeling of rate and size effects in quasibrittle fracture. International Journal of Fracture 62, 355–373. URL: https://doi.org/10.1007/BF 00017241 , doi: 10.1007/BF 00017241 . · doi ↗

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Taxonomy

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

Abstract

keywords:

1 Introduction

2 Derivation of the second-order theory

2.1 Generalities

2.2 First-order expansion of the fundamental kernel

2.3 Second-order expansion of the SIF

3 Validation of the proposed theory

3.1 SIF distribution along the front of an elliptical crack

3.2 SIF distribution along a sinusoidal crack front

4 Coplanar crack propagation in weakly disordered materials

4.1 Problem statement and field generation

4.2 Computation of the crack front position

4.3 Effective toughness of weakly disordered materials

5 Conclusion

CRediT authorship contribution statement

Online access to data

Acknowledgments

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

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

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

Linear stress field dependent on zzz

Linear stress field dependent on xxx

Quadratic stress field odd in zzz and xxx

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

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

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

Linear stress field dependent on $z$

Linear stress field dependent on $x$

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