Configurational stability of a crack propagating in a material with mode-dependent fracture energy - Part I: Mixed-mode I+III

TL;DR

This paper investigates how the stability of a propagating crack in a material depends on the relationship between fracture energy and mode mixity, proposing a new model that aligns better with experimental observations.

Contribution

It introduces a constitutive relationship between fracture energy and mode mixity ratio, improving the prediction of crack instability thresholds in mixed-mode I+III fracture.

Findings

Lowered the predicted instability threshold to match experimental data.

Derived formulas for stress intensity factors and energy release rate distributions.

Showed that mode-dependent fracture energy significantly influences crack stability.

Abstract

In a previous paper (Leblond et al., 2011), we proposed a theoretical interpretation of the experimentally well known instability of coplanar crack propagation in mode I+III. The interpretation relied on a stability analysis based on analytical expressions of the stress intensity factors for a crack slightly perturbed both within and out of its original plane, due to Gao and Rice (1986) and Movchan et al. (1998), coupled with a double propagation criterion combining Griffith's energetic condition and principle of local symmetry. Under such assumptions instability modes were indeed evidenced for values of the mode mixity ratio of the mode III to mode I stress intensity factors applied remotely larger than some threshold depending only on Poisson's ratio. Unfortunately, the predicted thresholds were much larger than those generally observed for typical values of this material parameter. While the subcritical character of the nonlinear bifurcation from coplanar to fragmented fronts has been proposed as a possible explanation for this discrepancy (Chen et al., 2015), we propose here an alternative explanation based on the introduction of a constitutive relationship between the fracture energy and the mode mixity ratio, which is motivated by experimental observations. By reexamining the linear stability analysis of a planar propagating front, we show that such a relationship suffices, provided that it is strong enough, to lower significantly the threshold value of the mode mixity ratio for instability so as to bring it in a range more consistent with experiments. Interesting formulae are also derived for the distributions of the perturbed stress intensity factors and energy release rate, in the special case of perturbations of the crack surface and front obeying the principle of local symmetry and having reached a stationary state.

Equations

$${\rho }^{\mathrm{cr}}={\left[\frac{(1-\nu )(2-3\nu )}{3(2-\nu )-4\sqrt{2}(1-2\nu )}\right]}^{1/2}$$

$$\delta K_p(x,z)={\delta }_x{K}_p(x,z)+{\delta }_y{K}_p(x,z)(p=I,II,III).$$

$$\left\{\begin{array}[]{lll}\displaystyle\delta_{x}K_{I}(x,z)&=&\displaystyle\frac{K_{I}^{0}}{2\pi}\,PV\int_{-\infty}^{+\infty}\frac{\phi_{x}(x,z^{\prime})-\phi_{x}(x,z)}{(z^{\prime}-z)^{2}}\,dz^{\prime}\\[8.53581pt] \displaystyle\delta_{x}K_{II}(x,z)&=&\displaystyle-\frac{2}{2-\nu}K_{III}^{0}\,\frac{\partial\phi_{x}}{\partial z}(x,z)+\frac{2-3\nu}{2-\nu}\,\frac{K_{II}^{0}}{2\pi}\,PV\int_{-\infty}^{+\infty}\frac{\phi_{x}(x,z^{\prime})-\phi_{x}(x,z)}{(z^{\prime}-z)^{2}}\,dz^{\prime}\\[8.53581pt] \displaystyle\delta_{x}K_{III}(x,z)&=&\displaystyle\frac{2(1-\nu)}{2-\nu}K_{II}^{0}\frac{\partial\phi_{x}}{\partial z}(x,z)+\frac{2+\nu}{2-\nu}\,\frac{K_{III}^{0}}{2\pi}\,PV\int_{-\infty}^{+\infty}\frac{\phi_{x}(x,z^{\prime})-\phi_{x}(x,z)}{(z^{\prime}-z)^{2}}\,dz^{\prime}\\ \end{array}\right.$$

$$\left\{\begin{array}[]{lll}\displaystyle\delta_{y}K_{I}(x,z)&=&\displaystyle-\frac{3}{2}K_{II}^{0}\frac{\partial\phi_{y}}{\partial x}(x,z)-2K_{III}^{0}\,\frac{\partial\phi_{y}}{\partial z}(x,z)\\[8.53581pt] {}&{}&\displaystyle-\frac{K_{II}^{0}}{2\pi}\,PV\int_{-\infty}^{+\infty}\frac{\phi_{y}(x,z^{\prime})-\phi_{y}(x,z)}{(z^{\prime}-z)^{2}}\,dz^{\prime}+\delta_{y}K_{I}^{\rm skew}(x,z)\\[8.53581pt] \displaystyle\delta_{y}K_{II}(x,z)&=&\displaystyle\frac{K_{I}^{0}}{2}\,\frac{\partial\phi_{y}}{\partial x}(x,z)-\frac{2-3\nu}{2-\nu}\frac{K_{I}^{0}}{2\pi}\,PV\int_{-\infty}^{+\infty}\frac{\phi_{y}(x,z^{\prime})-\phi_{y}(x,z)}{(z^{\prime}-z)^{2}}\,dz^{\prime}\\[8.53581pt] \displaystyle\delta_{y}K_{III}(x,z)&=&\displaystyle\frac{2(1-\nu)^{2}}{2-\nu}K_{I}^{0}\,\frac{\partial\phi_{y}}{\partial z}(x,z).\\ \end{array}\right.$$

$${\delta }_y{K}_I^{\mathrm{skew}}(x,z)=\frac{\sqrt{2}}{4\pi }\frac{1-2\nu }{1-\nu }\mathrm{Re}\left\{{\int }_{-\infty }^{x}d{x}^{\prime }{\int }_{-\infty }^{+\infty }\frac{[K_{III}^0-i(1-\nu )K_{II}^0](\partial {\varphi }_y/\partial z)(x^{\prime },z^{\prime })}{(x-x^{\prime }{)}^{1/2}{\left[x-x^{\prime }+i(z-z^{\prime })\right]}^{3/2}}d{z}^{\prime }\right\}$$

$$G_c(x,z)\equiv G_c[\rho (x,z)],\rho (x,z)\equiv \frac{K_{III}(x,z)}{K_I(x,z)}(>0).$$

$$G_c(\rho )\equiv G_{cI}(1+\gamma {\rho }^{\kappa })$$

$$\left\{\begin{array}[]{lll}\phi_{x}(x,z)&=&e^{\lambda x}\,\psi_{x}(z)\\ \phi_{y}(x,z)&=&e^{\lambda x}\,\psi_{y}(z)\end{array}\right.$$

$$\chi (z)={\int }_{-\infty }^{+\infty }\hat{\chi }(k)e^{ikz}dk\iff \hat{\chi }(k)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }\chi (z)e^{-ikz}dz.$$

$$\xi \equiv \frac{\lambda }{|k|}(>0);$$

$$\widehat{{\delta }_y{K}_I^{\mathrm{skew}}}(x,k)=\frac{i}{\sqrt{2}}\frac{1-2\nu }{1-\nu }{K}_{III}^0{e}^{\lambda x}F(\xi )k\widehat{{\psi }_y}(k)\mathrm{where}F(\xi )\equiv \frac{1}{\sqrt{1+\xi }}.$$

$$\left\{\begin{array}[]{lll}\widehat{\delta K_{I}}(x,k)&=&\displaystyle e^{\lambda x}\left\{-K_{I}^{0}\frac{|k|}{2}\widehat{\psi_{x}}(k)+iK_{III}^{0}\,\left[-2+\frac{1-2\nu}{\sqrt{2}(1-\nu)}F(\xi)\right]k\widehat{\psi_{y}}(k)\right\}\\[8.53581pt] \widehat{\delta K_{II}}(x,k)&=&\displaystyle e^{\lambda x}\left\{-iK_{III}^{0}\frac{2}{2-\nu}k\widehat{\psi_{x}}(k)+K_{I}^{0}\left[\frac{2-3\nu}{2(2-\nu)}|k|+\frac{\lambda}{2}\right]\widehat{\psi_{y}}(k)\right\}\\[8.53581pt] \widehat{\delta K_{III}}(x,k)&=&\displaystyle e^{\lambda x}\left\{-K_{III}^{0}\,\frac{2+\nu}{2(2-\nu)}|k|\widehat{\psi_{x}}(k)+iK_{I}^{0}\frac{2(1-\nu)^{2}}{2-\nu}k\widehat{\psi_{y}}(k)\right\}.\end{array}\right.$$

$$\widehat{{\psi }_y}(k)=\frac{4i}{2-3\nu +(2-\nu )\xi }{\rho }^0\mathrm{sgn}(k)\widehat{{\psi }_x}(k),$$

$${\rho }^0\equiv \frac{K_{III}^0}{K_I^0}$$

$$\left\{\begin{array}[]{lll}\displaystyle\frac{\widehat{\delta G}(k)}{G^{0}}&=&-e^{\lambda x}f(\rho^{0};\xi)\,|k|\widehat{\psi_{x}}(k)\\ \widehat{\delta\rho}(k)&=&-e^{\lambda x}g(\rho^{0};\xi)\,|k|\widehat{\psi_{x}}(k)\end{array}\right.$$

$$G^0\equiv \frac{1-{\nu }^2}{E}(K_I^0{)}^2+\frac{1+\nu }{E}(K_{III}^0{)}^2$$

$$\left\{\begin{array}[]{lll}f(\rho;\xi)&\equiv&\displaystyle\frac{1-\nu}{1-\nu+\rho^{2}}\left\{1-\frac{3(2-\nu)-4\sqrt{2}(1-2\nu)F(\xi)-(2+\nu)\xi}{(1-\nu)\left[2-3\nu+(2-\nu)\xi\right]}\,\rho^{2}\right\}\\[14.22636pt] g(\rho;\xi)&\equiv&\displaystyle\frac{\rho}{2-3\nu+(2-\nu)\xi}\left\{4-5\nu+\nu\xi+2\left[4-\sqrt{2}\,\frac{1-2\nu}{1-\nu}F(\xi)\right]\rho^{2}\right\}\end{array}\right.\quad(\xi>0).$$

$$\left\{\begin{array}[]{lll}G(x,z)&=&G^{0}+\delta G(x,z)\\ G_{c}(x,z)&=&\displaystyle G_{c}(\rho^{0})+\frac{dG_{c}}{d\rho}(\rho^{0})\,\delta\rho(x,z)\end{array}\right.$$

$$G^0=G_c({\rho }^0)(=G_{cI}[1+\gamma ({\rho }^0{)}^{\kappa }]).$$

$$\delta G(x,z)=\frac{d{G}_c}{d\rho }({\rho }^0)\delta \rho (x,z)\Rightarrow$$

$$G^0{e}^{\lambda x}f({\rho }^0;\xi )|k|\widehat{{\psi }_x}(k)=\frac{d{G}_c}{d\rho }({\rho }^0)e^{\lambda x}g({\rho }^0;\xi )|k|\widehat{{\psi }_x}(k)(\forall k).$$

$$f({\rho }^0;\xi )-\frac{d(\ln G_c)}{d\rho }({\rho }^0)g({\rho }^0;\xi )=0$$

$$\xi=\frac{N(\rho^{0};\xi)}{D(\rho^{0})}\ \ \mbox{where}\left\{\begin{array}[]{lll}N(\rho;\xi)&\equiv&\displaystyle-(1-\nu)(2-3\nu)+\left[3(2-\nu)-4\sqrt{2}(1-2\nu)F(\xi)\right]\rho^{2}\\[8.53581pt] {}&{}&\displaystyle+X(\rho)\left(1-\nu+\rho^{2}\right)\left\{4-5\nu+2\left[4-\sqrt{2}\frac{1-2\nu}{1-\nu}F(\xi)\right]\rho^{2}\right\}\\[8.53581pt] D(\rho)&\equiv&\displaystyle(1-\nu)(2-\nu)+(2+\nu)\rho^{2}-X(\rho)\nu\left(1-\nu+\rho^{2}\right)\\[8.53581pt] X(\rho)&\equiv&\displaystyle\frac{\kappa\gamma\rho^{\kappa}}{1+\gamma\rho^{\kappa}}\,.\end{array}\right.$$

$$\widehat{{\psi }_x}(k^{\prime })\equiv \frac{A_x}{2}\left[e^{i\theta }\delta (k^{\prime }-k)+e^{-i\theta }\delta (k^{\prime }+k)\right]$$

$${\psi }_x(z)={\int }_{-\infty }^{+\infty }\widehat{{\psi }_x}(k^{\prime })e^{i{k}^{\prime }z}dk=\frac{A_x}{2}\left[e^{i(kz+\theta )}+e^{-i(kz+\theta )}\right]=A_x\cos (kz+\theta ).$$

$$\widehat{{\psi }_y}(k^{\prime })=\frac{4i}{2-3\nu +(2-\nu )\xi }{\rho }^0\frac{A_x}{2}\left[e^{i\theta }\delta (k^{\prime }-k)-e^{-i\theta }\delta (k^{\prime }+k)\right].$$

$$\begin{array}[]{lll}\psi_{y}(z)&=&\displaystyle\int_{-\infty}^{+\infty}\widehat{\psi_{y}}(k^{\prime})\,e^{ik^{\prime}z}dk=\frac{4i}{2-3\nu+(2-\nu)\xi}\,\rho^{0}\,\frac{A_{x}}{2}\left[e^{i(kz+\theta)}-e^{-i(kz+\theta)}\right]\\ {}&=&\displaystyle-\frac{4}{2-3\nu+(2-\nu)\xi}\,\rho^{0}A_{x}\sin(kz+\theta).\end{array}$$

$$\left\{\begin{array}[]{lll}\phi_{x}(x,z)&=&A_{x}e^{\xi kx}\cos(kz+\theta)\\ \phi_{y}(x,z)&=&A_{y}e^{\xi kx}\sin(kz+\theta)\end{array}\right.\quad{\rm where}\quad\frac{A_{y}}{A_{x}}\equiv-\frac{4}{2-3\nu+(2-\nu)\lambda/{k}}\,\rho^{0}.$$

$$N(\rho ;0)\simeq -(1-\nu )(2-3\nu )+\frac{\kappa \gamma {\rho }^{\kappa }}{1+\gamma {\rho }^{\kappa }}(1-\nu )(4-5\nu ).$$

$${\rho }^{\mathrm{cr}}\simeq {\left[\frac{2-3\nu }{(4-5\nu )\kappa -(2-3\nu )}\frac{1}{\gamma }\right]}^{1/\kappa }(\gamma \gg 1).$$

Full text

Keywords : Configurational stability; mode I+III; mode-dependent fracture energy; Griffith’s criterion; principle of local symmetry

, ,

,

1 Introduction

Exploring the configurational stability of cracks loaded under mixed-mode conditions is of major interest, as it controls not only the patterns produced by fracture but also the resistance of the material to crack propagation. It amounts to studying the stability of the geometrical configuration of a crack as it propagates under some arbitrary combination of modes. (The problem of configurational stability thus defined should not be confused with that of propagation stability, consisting of the alternative between arrest and acceleration of a crack propagating under constant loading).

The instability of the coplanar configuration of a crack propagating under mixed mode I+III conditions is well-documented experimentally; see the works of Sommer (1969); Knauss (1970); Palaniswamy and Knauss (1975); Hourlier and Pineau (1979); Pollard et al. (1982); Suresh and Tschegg (1987); Pollard and Aydin (1988); Yates and Miller (1989); Hull (1993, 1995); Cooke and Pollard (1996); Lazarus (1997); Lazarus et al. (2001b, 2008); Lin et al. (2010); Goldstein and Osipenko (2012); Pham and Ravi-Chandar (2014), to quote just a few. Cracks loaded in mode I+III are known to propagate in the form of small facets inclined over the average plane of propagation. Materials as diverse as glass (Sommer, 1969; Pham and Ravi-Chandar, 2014), steels (Hourlier and Pineau, 1979; Yates and Miller, 1989; Lazarus, 1997), rocks (Pollard et al., 1982; Pollard and Aydin, 1988; Cooke and Pollard, 1996), alumina (Suresh and Tschegg, 1987), PMMA (Lazarus et al., 2008), gypsum and cheese (!) (Goldstein and Osipenko, 2012) have been considered in experiments; this brings strength to the idea that the causes of the fragmentation of the crack front are not rooted in microscopic, material-dependent mechanisms, so that the standard macroscopic theory of Linear Elastic Fracture Mechanics (LEFM) should be able to handle the problem. But the breaking of translational invariance in the direction of the crack front resulting from the formation of facets, and the very complex, fully 3D associated crack geometry, greatly complicate the analysis.

Following a long period during which theoretical analyses of the problem remained scarce and incomplete (the works of Lazarus (1997) and Lazarus et al. (2001a, b) are typical examples), a more satisfactory - though still imperfect - approach was proposed by Leblond et al. (2011), drawing inspiration from the results of numerical simulations of Pons and Karma (2010) based on a phase-field model of Karma et al. (2001). The work consisted of a rigorous linear stability analysis of coplanar propagation of the crack under mixed mode I+III conditions, which relied on two basic elements:

• •

rigorous first-order expressions of the stress intensity factors (SIFs) along the front of a crack slightly perturbed both within and out of its original plane, due to Gao and Rice (1986) for the in-plane perturbation of the crack front, and Movchan et al. (1998) for the out-of-plane perturbation of the crack surface;

• •

use of a double propagation criterion combining Griffith (1920)’s energetic condition and Goldstein and Salganik (1974)’s Principle of Local Symmetry (PLS).

Instability modes - perturbed configurations of the crack surface and front growing without bound in time while continuously satisfying the double criterion at all instants and positions along the crack front - were found for values of the mode mixity ratio $\rho^{0}$ - ratio of the unperturbed mode III SIF $K_{III}^{0}$ to the unperturbed mode I SIF $K_{I}^{0}$ - larger than a threshold value $\rho^{\rm cr}$ given by

[TABLE]

where $\nu$ denotes Poisson’s ratio. This quantity is in the range of $0.4$ to $0.5$ for typical values of $\nu$ for glass and polymers. Unfortunately it is much larger than most threshold values of the mode mixity ratio actually observed. Although an exceptionally high value of $0.57$ has been reported by Eberlein et al. (2017) in some aluminium alloy, most of the time small or very small thresholds have been observed: for instance Sommer (1969)’s estimate for glass amounted to a few percent, and Pham and Ravi-Chandar (2014) even claimed that there was no threshold at all in both Homalite 100 and glass.

Several tentative explanations of this discrepancy have been proposed. A natural idea is that imperfections of the crack front geometry and/or loading could promote deviations from the fundamental coplanar solution below the threshold. This idea was first put to test by Leblond and Lazarus (2015) using Cotterell and Rice (1980)’s heuristic 2D directional stability criterion, suitably extended to the 3D case; it was shown that small accidental undulations of the crack front within the crack plane, due for instance to small heterogeneities of toughness, may generate a positive non-singular stress in the direction of propagation, thus causing the crack to deviate from coplanarity. A more complete investigation of the idea was provided by Chen et al. (2015) in the form of an extensive numerical study of non-coplanar solutions using Karma et al. (2001)’s phase-field model, which has been shown by Hakim and Karma (2009) to be equivalent to combining the Griffith condition and PLS in isotropic media. This study demonstrated that the bifurcation accompanying the transition from coplanar to fragmented front is strongly subcritical, suggesting that large enough perturbations may allow the system to jump from the stable branch to the unstable one well below the theoretical threshold.

The essential aim of this paper is to explore still another possibility, namely the fact that the critical-energy-release rate $G_{c}$ might not be a constant but depend upon the mode mixity ratio $\rho\equiv K_{III}/K_{I}$. The dependence of the fracture energy upon the ratio $K_{II}/K_{I}$ is well documented for interfacial fracture between two distinct materials (Freund et al., 2003); this effect plays an important role on the type of fracture patterns emerging from the delamination of thin films from a rigid substrate (Faou et al., 2017). For bulk fracture under mixed-mode I+III, such an assumption is suggested by the general observation that crack propagation is more difficult under such conditions than in pure mode I. The effect was for instance quantified by Lin et al. (2010); Fig. 1 shows these authors’ experimental measurement of the critical stress intensity factors at initiation under mode I + III conditions. One clearly sees that the fracture energy required to propagate a crack is higher in mode I + III than in pure mode I. The assumption of a mode-dependent $G_{c}$ is also supported indirectly by the fact that as mentioned above, widely different thresholds have been observed in different materials; the modest variations of Poisson’s ratio from one material to another do not seem capable of explaining such large variations, which suggests that the threshold may depend on additional material parameters.

The dependence of the critical energy-release-rate $G_{c}$ upon $\rho=K_{III}/K_{I}$, just like its dependence upon $K_{II}/K_{I}$ in interfacial fracture, might arise from the influence of shear on the damage mechanisms occurring within the process zone, at scales disregarded by LEFM. It could also originate from shear-dependent mechanisms taking place at a scale larger than that of the process zone, though still smaller than the “local” scale of the stability analysis. For example, possible formation under mode I+III conditions of tilted facets at a sub-continuum-mechanics scale could be invoked, a view suggested by the observation of self-similar fragmentation patterns consisting of facets at various length scales along the crack front reported by Pham and Ravi-Chandar (2014) and also apparent in Chen et al. (2015)’s experimental observations. The aim of this paper is not to study the physical origin of the dependence of $G_{c}$ on $\rho$, but to assume its existence as a purely heuristic hypothesis, and explore its consequences on linear stability. The physical intuitive idea is that a shear-dependent fracture energy may promote fragmentation, since fragmentation may result in smaller local values of $K_{III}$ and consequently less energy required for fracture.

Furthermore, to link our theoretical predictions to experiments, we use the experimental data of Lin et al. (2010) to estimate the strength of the dependence of $G_{c}$ on $\rho$ at the local scale as shown in Fig. 1. This procedure is only approximate given that those measurements are made on a crack front that is already fragmented, at least on some small initial scale, before actual propagation of the crack. However, a possible justification of it lies in the aforementioned possibility that the $G_{c}(\rho)$ relationship may originate from some shear-dependent mechanism operating at scales larger than that of the process zone but smaller than that of the instability wavelength, assumed to set a larger facet scale for propagating fragmented fracture fronts. This hypothesized separation of scales between the mechanism underlying the $G_{c}(\rho)$ relationship and the initial front fragmentation scale induced by mode I+III crack propagation would presumably continue to hold during the subsequent facet coarsening process that further increases the fragmentation scale (Chen et al., 2015; Leblond et al., 2015). However facet coarsening is not investigated in the present study, limited to a linear stability analysis applicable only to the beginning of the development of unstability modes.

The paper is organized as follows:

• •

Section 2 presents general hypotheses and notations, and formulae - serving as a basis for the analyses to follow - for the SIFs along the front of a slightly perturbed semi-infinite crack, due to Gao and Rice (1986) for the in-plane perturbation of the crack front, and Movchan et al. (1998) for the out-of-plane perturbation of the crack surface.

• •

Then Section 3 extends the stability analysis of Leblond et al. (2011) of coplanar crack propagation in mixed-mode I+III by including a dependence of the fracture energy upon mode mixity, and examines the impact of this dependence upon the threshold value of the mode mixity ratio.

• •

As a complement, Section 4 investigates some remarkable properties of the perturbed SIFs and energy-release-rate along the crack front in the special case of stationary perturbations obeying the PLS, corresponding to instability modes in near-threshold conditions.

2 First-order perturbation of a semi-infinite crack in an infinite body

We consider an infinite body made of some isotropic elastic material and containing a semi-infinite crack. In the initial, unperturbed configuration (Fig. 2), this crack is planar and its front is straight. The usual convention for axes is used: the origin $O$ is chosen arbitrarily within the crack plane, the axis $Ox$ oriented along the direction of propagation, the axis $Oy$ perpendicular to the crack plane and the axis $Oz$ parallel to the crack front. The crack is loaded in general mixed-mode I+II+III conditions through some system of forces independent of the coordinate $z$, so that the unperturbed SIFs $K_{I}^{0}$, $K_{II}^{0}$, $K_{III}^{0}$ are independent of the position of their point of observation along the crack front.

In the perturbed configuration (Fig. 3), the crack front is perturbed within the original crack plane by a small amount $\phi_{x}(x,z)$ in the direction $Ox$ (Fig. 3(a)), and the crack surface is perturbed out of the original crack plane by a small amount $\phi_{y}(x,z)$ in the direction $Oy$ (Fig. 3(b)); the argument $x$ in the functions here represents the average position of the perturbed crack front in the direction of propagation. (Note that in Fig. 3(b) the perturbation is limited to the immediate vicinity of the crack front for legibility but can in fact extend over the entire crack surface). The perturbation of the $p$-th SIF ($p=I,II,III$) resulting from the in-plane perturbation $\phi_{x}(x,z)$ of the crack front is denoted $\delta_{x}K_{p}(x,z)$, and that resulting from the out-of-plane perturbation $\phi_{y}(x,z)$ of the crack surface, $\delta_{y}K_{p}(x,z)$. To first order in the pair $(\phi_{x},\,\phi_{y})$, these perturbations are additive so that the total perturbation of the $p$-th SIF is simply

[TABLE]

We shall apply the results of Gao and Rice (1986) on the in-plane perturbations of the SIFs, $\delta_{x}K_{p}(x,z)$, and those of Movchan et al. (1998) on their out-of-plane perturbations, $\delta_{y}K_{p}(x,z)$. To do so, we introduce the assumption that the characteristic length defined by the loading, in the absence of any length scale defined by the infinite geometry itself, is much larger than the typical distances of variation of the perturbations $\phi_{x}(x,z)$ and $\phi_{y}(x,z)$ of the crack front and surface. It is then possible to retain only these terms, in Gao and Rice (1986)’s and Movchan et al. (1998)’s formulae, which involve the unperturbed SIFs, discarding those involving the unperturbed non-singular stresses and higher-order constants characterizing the initial near-front stress field. Gao and Rice (1986)’s formulae then read

[TABLE]

where the symbol $PV$ denotes a Cauchy Principal Value; and Movchan et al. (1998)’s formulae read

[TABLE]

In the expression of $\delta_{y}K_{I}(x,z)$ here the quantity $\delta_{y}K_{I}^{\rm skew}(x,z)$ is connected to Bueckner (1987)’s skew-symmetric crack-face weight functions - whence the notation. Movchan et al. (1998)’s calculation of this term for a perturbation $\phi_{y}(x,z)$ independent of $x$ has been extended by Leblond et al. (2011) to an arbitrary perturbation, with the following result:

[TABLE]

where the cut of the complex power function is along the half-line of negative real numbers.

3 Linear stability analysis

3.1 Generalities

From now on the unperturbed SIF $K_{II}^{0}$ of mode II will be assumed to be zero so that the crack will be loaded in mixed-mode I+III. Without any loss of generality, we may assume the unperturbed SIF $K_{III}^{0}$ of mode III to be positive like that of mode I, $K_{I}^{0}$.

Like in our previous work (Leblond et al., 2011), the prediction of the crack path will be based on a double propagation criterion enforced at all points of the crack front and all instants, consisting of

• •

Griffith (1920)’s condition $G(x,z)=G_{c}(x,z)$ where $G(x,z)$ denotes the local energy-release-rate, and $G_{c}(x,z)$ the local critical value of this quantity;

• •

Goldstein and Salganik (1974)’s PLS stipulating that the local SIF $K_{II}(x,z)$ of mode II must be zero.

However, whereas the critical energy-release-rate $G_{c}$ was supposed to be a constant in the previous work, it will be assumed here to be a function of the mode mixity ratio $\rho$ (ratio of the local mode III to mode I SIFs):

[TABLE]

Special attention will be paid to the case where the function $G_{c}(\rho)$ is of the form

[TABLE]

where $G_{cI}$ denotes the value of $G_{c}$ in pure mode I ($\rho=0$), and $\gamma$ and $\kappa$ dimensionless material parameters.111Note that formula (7) is given only for positive values of the mode mixity ratio $\rho$; for arbitrary values it should be applied with $|\rho|$ instead of $\rho$, the function $G_{c}(\rho)$ being even for obvious symmetry reasons. These parameters are assumed to be positive, meaning that the presence of mode III is assumed to increase the value of the critical energy-release-rate, in line with Freund et al. (2003)’s observation that the presence of shear systematically results in an increase of the cracking resistance of interfaces between dissimilar materials. A quadratic variation of the fracture energy with mode mixity - $\kappa=2$ - would be the most natural assumption, as it corresponds to the leading term in a development of the even function $G_{c}(\rho)$ for small values of $\rho$; such a value is generally chosen to describe the variations of $G_{c}$ with the ratio $K_{II}/K_{I}$ in interfacial fracture, see e.g. Faou et al. (2017). But we choose here to keep $\kappa$ as a free parameter so as not to restrict the generality of the stability analysis.

The analysis will be based on consideration of instability modes consisting of in-plane perturbations of the crack front and out-of-plane perturbations of the crack surface of the form

[TABLE]

where $\lambda$ is a positive parameter homogeneous to the inverse of a length, characterizing the growth rate of the mode in the direction of propagation, and $\psi_{x}(z)$, $\psi_{y}(z)$ functions to be determined. (The case $\lambda<0$ could be studied as well, but is of little interest since it corresponds to banal stability of the coplanar configuration of the propagating crack).

Use will be made of Fourier transforms in the direction $z$ of the crack front; the definition adopted here for the Fourier transform $\widehat{\chi}(k)$ of an arbitrary function $\chi(z)$ is

[TABLE]

We then define a dimensionless “normalized growth rate” $\xi$ of the instability mode by the formula

[TABLE]

this parameter compares the growth rate of the instability mode, $\lambda$, in the direction $x$, to the wavenumber of the Fourier component considered, $|k|$, in the direction $z$.

3.2 Fourier transforms of the perturbations of the energy-release-rate and the mode mixity ratio

To calculate the Fourier transforms $\widehat{\delta G}$, $\widehat{\delta\rho}$ of the first-order variations $\delta G$, $\delta\rho$ of the energy-release-rate and mode mixity ratio, the first task is to calculate the Fourier transforms $\widehat{\delta K_{p}}$ of the variations of the SIFs $\delta K_{p}$ ($p=I,II,III$). The simplest way of doing so is to insert the expressions (8) of the perturbations $\phi_{x}(x,z)$, $\phi_{y}(x,z)$ into equations (2 - 5), use the expressions (9)1 of the functions $\psi_{x}(z)$, $\psi_{y}(z)$ in terms of their Fourier transforms $\widehat{\psi_{x}}(k)$, $\widehat{\psi_{y}}(k)$, and evaluate the integrals for each value of $k$. The calculations are elementary except for that of the term $\delta_{y}K_{I}^{\rm skew}$, which is presented in Appendix A with the following result:222In our previous stability analysis (Leblond et al., 2011), $F$ was considered to be a function of $|k|/\lambda$ instead of $\xi={\lambda}/{|k|}$, and therefore expressed as $\sqrt{\frac{|k|/\lambda}{|k|/\lambda+1}}$ rather than $1/\sqrt{1+{\lambda}/{|k|}}$.

[TABLE]

With the result (11), the following expressions of the Fourier transforms of the variations $\delta K_{p}$ ($p=I,II,III$) of the SIFs are obtained:

[TABLE]

3.3 Application of the principle of local symmetry

Goldstein and Salganik (1974)’s PLS, written in “Fourier’s form” $\widehat{\delta K_{II}}(x,k)=0$, then yields

[TABLE]

where

[TABLE]

denotes the unperturbed mode mixity ratio and ${\rm sgn}(x)$ the sign of $x$.

3.4 Application of Griffith (1920)’s criterion

Using then equations (12, 13), Irwin’s formula $G=\frac{1-\nu^{2}}{E}(K_{I}^{2}+K_{II}^{2})+\frac{1+\nu}{E}K_{III}^{2}$ ($E$: Young’s modulus) and the definition (6)2 of the mode mixity ratio, one can get the Fourier transforms $\widehat{\delta G}$, $\widehat{\delta\rho}$ of the perturbations of the energy-release-rate and mode mixity ratio, as functions of the Fourier transform $\widehat{\psi_{x}}$ of the sole in-plane perturbation of the crack front. The results, obtained after straightforward but heavy calculations, are as follows:

[TABLE]

where

[TABLE]

denotes the unperturbed energy-release-rate, and $f$ and $g$ the functions defined by

[TABLE]

To enforce Griffith (1920)’s criterion $G(x,z)=G_{c}(x,z)$ at every point of the crack front and every instant, we use the first-order expansions of $G(x,z)$ and $G_{c}(x,z)$ in the pair $(\phi_{x},\phi_{y})$,

[TABLE]

where $\delta G$ and $\delta\rho$ are given, in Fourier’s form, by equations (15). Equating successive terms of the expansions of $G(x,z)$ and $G_{c}(x,z)$, we get:

• •

At order 0:

[TABLE]

This condition specifies the intensity of the loading inducing propagation of the crack, for the given value of the unperturbed mode mixity ratio $\rho^{0}$.

• •

At order 1:

[TABLE]

[TABLE]

There are then two possibilities for each value of the wavenumber $k$: either the Fourier component $\widehat{\psi_{x}}(k)$ is zero and this relation is trivially satisfied; or it is not, and then necessarily

[TABLE]

where equation (19) has been used. This condition specifies the value of the normalized growth rate $\xi={\lambda}/{|k|}$ of the instability mode as a function of the unperturbed mode mixity ratio $\rho^{0}=K_{III}^{0}/K_{I}^{0}$.

3.5 Growth rate of the instability mode and critical mode mixity ratio

Equation (20) may be put in an alternative format by “solving” it with respect to $\xi$, formally considering the value of $F(\xi)$ as known; the function $G_{c}(\rho)$ being assumed to be of the form (7), we get upon use of equations (17):

[TABLE]

The terms $3(2-\nu)-4\sqrt{2}(1-2\nu)F(\xi)$ and $4-\sqrt{2}\frac{1-2\nu}{1-\nu}F(\xi)$ in the expression of the numerator $N(\rho;\xi)$ here vary modestly with $\xi$,333For a standard value of $\nu$ of $0.3$, when $\xi$ goes from [math] to $+\infty$, the first term varies from $2.84$ to $5.1$, and the second from $3.19$ to $4$. In the limit $\nu\rightarrow 1/2$, they both become constant. so that they may be considered as approximately constant. Thus equation (21) “almost” directly provides the value of the normalized growth rate $\xi$, which makes it convenient for both qualitative physical analyses of its predictions, and simple methods of numerical solution such as the fixed-point algorithm. The form (21) is also convenient for a mathematical analysis of the solutions in $\xi$; such an analysis is presented in Appendix B, with the following conclusions:

The equation $N(\rho^{0};0)=0$ admits a unique positive solution $\rho^{\rm cr}$. This solution represents the critical value of the unperturbed mode mixity ratio for which a neutrally stable mode (having a zero normalized growth rate $\xi$) exists. 2. 2.

Assume that $\kappa\leq 3$. Then for values of the unperturbed mode mixity ratio $\rho^{0}$ larger than $\rho^{\rm cr}$, the equation $\xi={N(\rho^{0};\xi)}/{D(\rho^{0})}$ on $\xi$ admits a unique positive solution. This means that an instability mode exists, with normalized growth rate given by this solution. 3. 3.

Assume that $\kappa\leq 3$ and $\nu\geq\frac{2}{31}(9-5\sqrt{2})\simeq 0.124$.444This value may not be optimal, since the proof of Appendix B is based on strict inequalities which may be liable to improvement. Then for values of $\rho^{0}$ smaller than $\rho^{\rm cr}$, the equation $\xi={N(\rho^{0};\xi)}/{D(\rho^{0})}$ on $\xi$ does not admit any positive solution. Hence no instability mode exists.

It is finally worth noting that when the critical energy-release-rate $G_{c}$ is a constant ($\gamma=0$), equation (21) for the normalized growth rate $\xi$ of the instability mode, and the equation $N(\rho^{0};0)=0$ defining the critical mode mixity ratio $\rho^{\rm cr}$, are identical to those obtained in our previous stability analysis, equations (20) and (21) of (Leblond et al., 2011). They are however new in general.

3.6 Geometry of the instability mode

Assuming now that $\rho^{0}\geq\rho^{\rm cr}$ so that equation (21) on $\xi$ admits a solution, we discuss the geometrical shape of the corresponding instability mode. Two cases must be distinguished.

• •

Case 1: $\rho^{0}>\rho^{\rm cr}$.As established in Appendix B, equation (21) admits a single positive solution in $\xi=\frac{\lambda}{|k|}$. For a fixed (positive) value of the growth rate $\lambda$, this value corresponds to a pair of distinct solutions $(k,-k)$ with $k>0$. Then the Fourier transform $\widehat{\psi_{x}}(k^{\prime})$ of the in-plane perturbation of the crack front is nonzero only for $k^{\prime}=k$ or $k^{\prime}=-k$. Since the perturbation $\psi_{x}$ is real, $\widehat{\psi_{x}}(-k^{\prime})=\overline{\widehat{\psi_{x}}(k^{\prime})}$, and it follows that $\widehat{\psi_{x}}(k^{\prime})$ must be of the form $\frac{1}{2}\left[C_{x}\,\delta(k^{\prime}-k)+\overline{C_{x}}\,\delta(k^{\prime}+k)\right]$ where $C_{x}$ is a complex number and $\delta$ Dirac’s function, or equivalently

[TABLE]

where $A_{x}$ and $\theta$ are real numbers. Then by equation (9)1, the in-plane perturbation $\psi_{x}$ of the crack front is given by

[TABLE]

By equation (13), the expression of the Fourier transform $\widehat{\psi_{y}}(k^{\prime})$ of the out-of-plane perturbation of the crack surface then reads

[TABLE]

This implies that the out-of-plane perturbation $\psi_{y}$ of the crack surface is given by

[TABLE]

Combination of equations (8), (23) and (25) yields the final expression of the instability mode:

[TABLE]

Equations (26) show that the instability mode corresponds to a perturbed crack front having the shape of an elliptic helix, of axis coinciding with the present unperturbed crack front and “amplitudes” $A_{x}e^{\xi kx}$, $A_{y}e^{\xi kx}$ in the directions $x$ and $y$ respectively, growing in proportion, and exponentially with the distance of propagation $x$. This conclusion, and even the precise relation (26)3 between the amplitudes $A_{x}$ and $A_{y}$, are exactly the same as in the case of a critical energy-release-rate independent of the mode mixity ratio, see (Leblond et al., 2011).555This is because the equations derived in this Subsection result only from the PLS which bears no relationship to the energy-release-rate. Fig. 4 illustrates the geometry of the instability mode.

• •

Case 2: $\rho^{0}=\rho^{\rm cr}$.Then the solution of equation (21) is $\xi=\frac{\lambda}{|k|}=0$ so that $\lambda=0$ and $k$ is arbitrary. The out-of-plane perturbation $\phi_{y}$ is independent of the distance $x$ of propagation, and since the components of its Fourier transform in the direction $z$ of the crack front are totally unconstrained, it varies arbitrarily in this direction. This means that coplanar propagation is (neutrally) unstable versus arbitrary out-of-plane deviations from coplanarity. On the other hand, once the out-of-plane perturbation $\phi_{y}$ is fixed, the in-plane perturbation $\phi_{x}$ is also fixed and determined by equation (13) (where $\xi=0$), which connects the two perturbations via the PLS.

3.7 Numerical illustrations

In this Section we provide numerical examples of predictions of the values of the critical mode mixity ratio and the normalized growth rate of unstable perturbations. A wide range of possible material parameters is considered, without any reference to a specific material.

Fig. 5 illustrates the dependence of the critical mode mixity ratio $\rho^{\rm cr}$ - obtained by numerically solving the equation $N(\rho^{\rm cr};0)=0$ - upon Poisson’s ratio $\nu$, for several values of $\gamma$ and $\kappa$.

It clearly appears that the dependence of $G_{c}$ upon $\rho$ has a major impact upon the value of the critical mode mixity ratio: the larger the value of $\gamma$, the smaller that of $\rho^{\rm cr}$, whereas the effect of $\kappa$ is more complex and depends upon the value of $\nu$.

More precisely, Fig. 5 seems to show that $\rho^{\rm cr}\ll 1$ for large values of $\gamma$. Quite remarkably, an approximate analytical expression of this critical value may be found in this case. Indeed when $\rho\ll 1$, the expression (21)2 of $N(\rho;\xi)$ becomes for $\xi=0$, using the approximations $-(1-\nu)(2-3\nu)+[...]\rho^{2}\simeq-(1-\nu)(2-3\nu)$, $1-\nu+\rho^{2}\simeq 1-\nu$, $4-5\nu+2(...)\rho^{2}\simeq 4-5\nu:$

[TABLE]

The solution $\rho^{\rm cr}$ of the equation $N(\rho^{0};0)=0$ is therefore given by

[TABLE]

Expression (27) applies provided that $(4-5\nu)\kappa-(2-3\nu)$ is positive, which is true for $\kappa\geq 1/2$ (with $0<\nu<1/2$). It confirms that $\rho^{\rm cr}$ goes to zero when $\gamma$ goes to infinity, as was anticipated from Fig. 5. It also implies that when both $\gamma$ and $\kappa$ are large, $\rho^{\rm cr}$ depends only marginally on $\nu$ since the exponent $1/\kappa$ becomes small.

Fig. 6 illustrates the dependence of the normalized growth rate $\xi$ of the instability mode - obtained by numerically solving equation (21) - upon the unperturbed mode mixity ratio $\rho^{0}$, for a few values of the material parameters. One clearly sees that $\xi$ takes positive values (implying an instability) only for values of $\rho^{0}$ larger than the material-dependent critical value $\rho^{\rm cr}$.

Finally, we make a tentative comparison with the experimental measurements of Lin et al. (2010). Assuming a value of $\kappa=2$, the value of $\gamma$ leading to a best fit (for small $\rho$) of the experimental crack initiation curve is $\gamma\approx 25$ (see Fig. 1). The critical threshold of instability for these material parameters is $\rho^{\rm cr}\sim 0.1$, which is in the range of loading considered in the experiments.

4 Perturbations of the energy-release-rate and the mode mixity ratio for stationary perturbations

The case of stationary perturbations (independent of $x$) satisfying the PLS (equation (13)) is of special interest since they correspond to instability modes at the onset of the instability ($\rho^{0}=\rho^{\rm cr}$, $\lambda=0$). It so happens that simple and explicit expressions, in the physical space rather than that of Fourier, may be obtained in this special case for the perturbations $\delta G$ and $\delta\rho$ of the energy-release-rate and mode mixity ratio.

Perturbations independent of $x$ are of the form (8) with $\lambda\equiv 0$, $\psi_{x}\equiv\phi_{x}$, $\psi_{y}\equiv\phi_{y}$; hence provided they satisfy the PLS, equations (15) are applicable with $\xi\equiv{\lambda}/{|k|}=0$ and yield

[TABLE]

where

[TABLE]

But equation (13) expressing the PLS implies that

[TABLE]

hence $\widehat{\delta G}$ and $\widehat{\delta\rho}$ may be rewritten in the form

[TABLE]

The factors $\frac{2-3\nu}{4\rho^{0}}\,f(\rho^{0};0)$ and $\frac{2-3\nu}{4\rho^{0}}\,g(\rho^{0};0)$ here are independent of the wavenumber $k$. Therefore Fourier inversion of these relations yields directly

[TABLE]

These remarkable equations directly relate the perturbations of the energy-release-rate and the mode mixity ratio to the local slope $d\phi_{y}/dz$ of the crack front in the plane perpendicular to the direction of propagation. They have interesting consequences. For instance, consider (Fig. 7) the case of a crack loaded in mode I+III (with $K_{III}^{0}>0$), with a surface shaped like a “factory roof”, that is consisting of “facets of type A” obtained through rotation of the original crack plane $Ozx$ by an angle $\theta_{A}>0$ about the direction $x$ of propagation, and “facets of type B” obtained through a similar rotation but with an angle $-\theta_{B}<0$.666The denominations “facets of type A”, “facets of type B” are due to Hourlier and Pineau (1979). Assume that these facets have reached a stationary state and satisfy the PLS. Then equations (29) apply, and imply that the perturbations of the energy-release-rate and the mode mixity ratio (and hence the energy-release-rate and the mode mixity ratio themselves) take constant values over facets of each type, given with obvious notations by

[TABLE]

The case of predominant mode I is of special interest. In such a case the mode mixity ratio is much smaller than unity so that the terms proportional to its square in the expressions (28) of $f(\rho;0)$ and $g(\rho;0)$ may be disregarded, and the expressions (29) of $\delta G$ and $\delta\rho$ take the following wonderfully simple form, for arbitrary out-of-plane perturbations $\phi_{y}$:

[TABLE]

A simple and interesting application of formula (312) consists in approximately predicting the inclination angle $\theta_{A}$ of facets of type A on a factory-roof-shaped crack surface (Fig. 7) having reached its stationary state. If we adopt, following Cooke and Pollard (1996), the assumption that facets of this type are shear-free, the resulting condition $\rho_{A}=\rho^{0}+\delta\rho_{A}=0$, combined with equation (31)2, yields

[TABLE]

This expression calls for a few comments.

• •

Since Movchan et al. (1998)’s formulae used in the derivation of equation (32) are accurate only to first order in the perturbations of the crack surface and front, the prediction of this equation is applicable only for small inclinations of facets of type A, resulting from small mode mixity ratios $\rho^{0}$.

• •

The inclination angle $\theta_{A}$ given by equation (32) differs from that predicted by Cooke and Pollard (1996)’s well-known formula, $\theta_{A}^{\rm CP}=\rho^{0}/(1-2\nu)$ (for small $\rho^{0}$), due to different hypotheses and methods of calculation. Cooke and Pollard (1996)’s formula was for incipient facets instead of stationary ones, and was based on an elementary calculation involving only the asymptotic stress field prior to any segmentation, rather than evaluation of the distribution of the local mode III SIF along an already fragmented crack front.

• •

Equation (32) provides a lower estimate of the inclination of facets of type A than Cooke and Pollard (1996)’s formula, agreeing better with the results of both experiments and numerical computations (based on a phase-field model) performed by Chen et al. (2015).

• •

Formula (32) eliminates the physically groundless divergence of the inclination angle predicted by Cooke and Pollard (1996)’s formula in the limit of incompressible materials ($\nu\rightarrow 1/2$).

A final remark pertains to the practical relevance of the factory roof geometry. This geometry has been observed in many experiments. It represents an asymptotically exact stationary solution of the linearized equations of crack propagation for $\rho^{0}=\rho^{\rm cr}$; but, as remarked in Subsection 3.6, the same is in fact true of all smooth small-amplitude out-of-plane perturbations.777This basically stems from the fact that all Fourier modes bifurcate at the same critical value $\rho^{\rm cr}$ for which the instability growth rate vanishes (that is, the bifurcation is infinite dimensional). Hence for $\rho^{0}=\rho^{\rm cr}$ any out-of-plane perturbation, being a sum of Fourier modes of different amplitudes, is a stationary solution. Understanding why the factory roof shape is dynamically selected in experiments and phase-field simulations, rather than any other shape compatible with linear stability analysis, is beyond the scope of this work. It would most likely require to consider both nonlinear effects and the role of the process zone scale, which has been proposed by Pons and Karma (2010) to modify the PLS by making the local mode II SIF dependent on out-of-plane curvature of the crack front, in a way analogous to surface tension in other interfacial pattern forming instabilities.

5 Summary of results and perspectives

In this paper, we presented an extension of our previous linear stability analysis of coplanar crack propagation in mixed-mode I+III (Leblond et al., 2011), including a heuristic hypothesis of dependence of the critical energy-release-rate $G_{c}$ upon the ratio $\rho\equiv K_{III}^{0}/K_{I}^{0}$ of the unperturbed mode III to mode I SIFs. The essential aim was to see whether accounting for a shear dependence of the fracture energy could lead to a significant decrease of the theoretical threshold $\rho^{\rm cr}$ above which coplanar propagation becomes unstable, and generate values low enough to become compatible with those generally observed experimentally.

Section 2 was first devoted to general hypotheses and notations and some fundamental formulae used throughout the paper, providing the SIFs along the front of a slightly perturbed semi-infinite crack in an infinite body. These formulae were taken from the work of Gao and Rice (1986) for the in-plane perturbation of the crack front, and the work of Movchan et al. (1998) (duly completed by that of Leblond et al. (2011)) for the out-of-plane perturbation of the crack surface.

Section 3 then presented the extended stability analysis. Apart from the additional hypothesis of dependence of $G_{c}$ upon $\rho$, the new analysis relied on the same basic elements as our previous one (Leblond et al., 2011); that is, in addition to the technical elements explained in the preceding Section, use of a double propagation criterion combining Griffith (1920)’s energetic condition and Goldstein and Salganik (1974)’s principle of local symmetry, enforced at all instants and all points along the perturbed crack front. The geometrical features of instability modes above the threshold were found to be the same as for a constant $G_{c}$, the perturbed crack front still assuming the shape of an elliptic helix with size growing exponentially with the distance of propagation. More importantly, the dependence of $G_{c}$ upon $\rho$ was found to generate a significant decrease of the theoretical threshold $\rho^{\rm cr}$, as was hoped. Furthermore, an estimate of this threshold obtained by assuming a simple functional form of $G_{c}(\rho)$ with parameters fitted to experimentally measured crack initiation curves of Lin et al. (2010) (cf. Fig. 1) was found to be compatible with the range of $\rho$-values for which fragmentation was observed in these experiments.

Section 4 was devoted to a complementary study of some remarkable properties of the SIFs and energy-release-rate along the perturbed crack front, in the special case of stationary perturbations satisfying the principle of local symmetry. Such perturbations are interesting in that they correspond to instability modes in near-threshold conditions, with negligible rate of growth. The perturbations of the SIFs and the energy-release-rate were shown to be both proportional to the local slope of the perturbed front in the plane perpendicular to the direction of propagation. For instance, for a crack surface having the shape of a “factory roof”, the SIFs of mode I and III and the energy-release-rate are all constant on every flat portion of the “roof”. It is hoped that the formulae derived will be useful in the future to analyze and understand the inclinations of the facets of type A and B observed experimentally, as functions of the mode mixity ratio and the material properties.

This work opens several perspectives:

• •

Theoretical predictions should be compared with the results of more fracture experiments performed under mixed mode I+III conditions. The first aim of experiments would be to study the dependence of the fracture energy upon the mode mixity ratio, and if possible precisely quantify this dependence independently of any observation of the threshold value of this ratio. One could then deduce the theoretical value of the threshold from this observed dependence, and compare it to the experimental value.

• •

The theoretical value of the threshold should be compared to the results of phase-field simulations analogous to those of Pons and Karma (2010) and Chen et al. (2015), but including an extra dependence of the fracture energy upon the mode mixity ratio. The aim of such simulations would be to account for nonlinear geometrical effects and the role of the process zone scale disregarded in the theoretical linear stability analysis, but which might be of primary importance, especially for disconnected facets.

• •

It would be interesting to investigate the question of the physical mechanisms responsible for the dependence of the fracture energy upon the mode mixity ratio, left completely open by the present work where this dependence was simply postulated as a heuristic hypothesis.

• •

Finally the stability analysis should be extended to general mixed-mode I+II+III conditions, since experiments almost always include some small, uncontrolled mode II loading component. The extension will inevitably be limited to small values of the mode II SIF applied remotely, since large values would result in a large general kink of the crack prohibiting the use of Movchan et al. (1998)’s geometrically linearized formulae.

Acknowledgements

The authors wish to express their sincere thanks to Professor K. Ravi-Chandar of the University of Austin, Texas for several very stimulating discussions. A.K. acknowledges support of Grant DEFG02-07ER46400 from the US Department of Energy, Office of Basic Energy Sciences. L.P. and A.V. acknowledge support of the City of Paris through the Emergence Program.

Appendix A Appendix : calculation of the Fourier transform of $\delta_{y}K_{I}^{\rm skew}$

To calculate the Fourier transform $\widehat{\delta_{y}K_{I}^{\rm skew}}(x,k)$ for an out-of-plane perturbation of the crack $\phi_{y}(x,z)$ of the form (8)2, the first task is to evaluate $\delta_{y}K_{I}^{\rm skew}(x,z)$ for the perturbation $e^{\lambda x}\,\widehat{\psi_{y}}(k)e^{ikz}$. For this perturbation equation (5) yields

[TABLE]

where the change of variable $v\equiv x-x^{\prime}$ has been used. The integral over $z^{\prime}$ here has been evaluated by Leblond et al. (2011) using complex analysis, and the result is:

[TABLE]

where $H$ is Heaviside’s function; it follows that

[TABLE]

Now use of the change of variable $w\equiv\sqrt{(|k|+\lambda)v}$ yields

[TABLE]

where $\xi$ and $F(\xi)$ are defined by equations (10) and (11)2 of the text. Therefore

[TABLE]

For the value $-k$, one has $\widehat{\psi_{y}}(-k)=\overline{\widehat{\psi_{y}}(k)}$ since $\psi_{y}(z)$ is real, so that

[TABLE]

For the perturbation $e^{\lambda x}.\,\frac{1}{2}\left[\widehat{\psi_{y}}(k)e^{ikz}+\widehat{\psi_{y}}(-k)e^{-ikz}\right]$, it then follows from use of the relation $H(k)+H(-k)=1$ that

[TABLE]

where the property $\widehat{\psi_{y}}(-k)=\overline{\widehat{\psi_{y}}(k)}$ has been used again.

Since the real perturbation defined by equation (8)2 may be written as

[TABLE]

it follows that for this real perturbation

[TABLE]

Comparison with the definition (9)1 of the Fourier transform then shows that

[TABLE]

which is equation (11) of the text.

Appendix B Appendix : qualitative mathematical analysis of equation (21) on $\xi$

It is assumed in this Appendix that $0<\nu<1/2$ - these inequalities being always satisfied in practice although thermodynamic stability requirements warrant only that $-1<\nu<1/2$.

Equation (21) on $\xi$ is written in the form

[TABLE]

B.1 Study of the function $N(\rho;0)$

The function $N(\rho;0)$ is given by equation (21)2 of the text with $\xi=0$, that is:

[TABLE]

where $X(\rho)$ is given by equation (21)4. The following properties follow:

• •

At the point $\rho=0$, $N(\rho;0)=-(1-\nu)(2-3\nu)$ is negative.

• •

In the limit $\rho\rightarrow+\infty$, $X(\rho)$ becomes constant so that $N(\rho;0)$ is positively proportional to $\rho^{4}$ and therefore goes to $+\infty$.

• •

The quantities

[TABLE]

are all positive and increase with $\rho$, so that $N(\rho;0)$ is an increasing function of $\rho$.

B.2 Consequences for the equation $N(\rho^{0};0)=0$

From the first two properties established in Subsection B.1, combined with the continuity of the function $N(\rho;0)$, one deduces that the equation $N(\rho^{0};0)=0$ admits at least one solution $\rho^{\rm cr}$ in the interval $(0,+\infty)$. From the third, one concludes that this solution is unique, and that $N(\rho;0)<0$ for $\rho<\rho^{\rm cr}$, $N(\rho;0)>0$ for $\rho>\rho^{\rm cr}$.

B.3 Positiveness of the denominator $D(\rho)$

The denominator $D(\rho)$ is given by equation (21)3 of the text. It follows from this expression, that of the quantity $X(\rho)$, equation (21)4, and the obvious inequality $\frac{\gamma\rho^{\kappa}}{1+\gamma\rho^{\kappa}}<1$, that

[TABLE]

For the last quantity to be positive for all values of $\rho$, it suffices that $2-\nu-\kappa\nu$ be positive ($2+\nu-\kappa\nu$ is then also automatically positive). This is equivalent to $\nu<\frac{2}{1+\kappa}$. Since $\nu<\frac{1}{2}$, it suffices that $\frac{1}{2}\leq\frac{2}{1+\kappa}$ or equivalently $\kappa\leq 3$ for this to be true. Such a condition, which seems experimentally reasonable (Lin et al., 2010), will be assumed to be fulfilled in the sequel. Then $D(\rho)$ is positive for all values of $\rho$.

B.4 Study of the function $\Phi(\rho;\xi)$ for fixed $\rho$

The following properties follow from the definition of the function $\Phi(\rho;\xi)$, the expressions (21)2,3 of $N(\rho;\xi)$ and $D(\rho)$, and the properties established in Subsections B.2 and B.3:

• •

At the point $\xi=0$, $\Phi(\rho;0)=-{N(\rho;0)}/{D(\rho)}$ is positive if $\rho<\rho^{\rm cr}$, and negative if $\rho>\rho^{\rm cr}$.

• •

In the limit $\xi\rightarrow+\infty$, $\Phi(\rho;\xi)\sim\xi$ so that this function goes to $+\infty$.

• •

The derivative $\partial\Phi/\partial\xi$ is given by

[TABLE]

since the quantities $D(\rho)$ and $[...]$ are positive, this derivative increases with $\xi$.

B.5 Consequences for the equation $\Phi(\rho^{0};\xi)=0$ for $\rho^{0}>\rho^{\rm cr}$

In the case where $\rho^{0}>\rho^{\rm cr}$, $\Phi(\rho^{0};\xi)$ is negative at the point $\xi=0$ and goes to $+\infty$ in the limit $\xi\rightarrow+\infty$; hence the equation $\Phi(\rho^{0};\xi)=0$ on $\xi$ admits at least one solution in the interval $(0,+\infty)$.

Let $\xi^{1}$ denote the first (smallest) solution. In the interval $(0,\xi^{1})$, $\Phi(\rho^{0};\xi)<0$; hence necessarily $\frac{\partial\Phi}{\partial\xi}(\rho^{0};\xi^{1})\geq 0$. Since $\frac{\partial\Phi}{\partial\xi}(\rho^{0};\xi)$ is an increasing function of $\xi$, this implies that $\frac{\partial\Phi}{\partial\xi}(\rho^{0};\xi)>0$ for $\xi>\xi^{1}$. Hence $\Phi(\rho^{0};\xi)$ increases with $\xi$ over the interval $(\xi^{1},+\infty)$ and cannot vanish there. Thus the solution of the equation $\Phi(\rho^{0};\xi)=0$ on $\xi$ in the interval $(0,+\infty)$ is unique.

B.6 Positivity of the derivative $\frac{\partial\Phi}{\partial\xi}(\rho;0)$ for $\rho<\rho^{\rm cr}$

We now assume that $\rho<\rho^{\rm cr}$ and look for conditions ensuring the positiveness of the derivative $\frac{\partial\Phi}{\partial\xi}(\rho;0)$. In order to avoid any ambiguity, we shall add to all desired, but not necessarily fulfilled inequalities, an interrogation mark “ ? ” after the inequality sign.

The above expression of $\frac{\partial\Phi}{\partial\xi}(\rho;\xi)$ shows that the inequality $\frac{\partial\Phi}{\partial\xi}(\rho;0)>?\ 0$ is equivalent to

[TABLE]

that is, account being taken of the expression (21)3 of $D(\rho)$:

[TABLE]

This inequality cannot be true for all values of $\rho$ since in the limit $\rho\rightarrow+\infty$, the right-hand side goes to $+\infty$ proportionally to $\rho^{4}$ ($X(\rho)$ goes to a constant), more quickly than the left-hand side which is asymptotically proportional to $\rho^{2}$. But we only wish it to be true for values of $\rho$ smaller than $\rho^{\rm cr}$. For such values we have seen that $N(\rho;0)<0$, which implies by the above expression of $N(\rho;0)$ that

[TABLE]

Hence to fulfill the desired inequality for values of $\rho$ smaller than $\rho^{\rm cr}$, it suffices to fulfill the inequality obtained by replacing the quantity $X(\rho)\left(1-\nu+\rho^{2}\right)$ by the upper bound just derived, that is:

[TABLE]

To study this inequality, rewrite it in the form

[TABLE]

where ${\mathcal{H}}(\rho)$ is the homographic function of $\rho^{2}$ defined by

[TABLE]

We wish to show that $0<{\mathcal{H}}(\rho)<1$. Since the denominator $4-5\nu+2\left(4-\sqrt{2}\,\frac{1-2\nu}{1-\nu}\right)\rho^{2}$ does not vanish, ${\mathcal{H}}$ is a monotone function of $\rho^{2}$, so that it suffices to consider its extremal values

[TABLE]

Both of these quantities are homographic, monotone functions of $\nu$; they vary between the extremal values (obtained for $\nu=0$ and $\nu=1/2$ respectively) [math] and $1/3$, $\frac{\sqrt{2}}{2(4-\sqrt{2})}\simeq 0.273$ and [math], respectively. All of these values are in the interval $(0,1)$, so the same is true of the quantities ${\mathcal{H}}(0)$ and ${\mathcal{H}}(+\infty)$, and of ${\mathcal{H}}(\rho)$ itself for all values of $\rho$.

Now rewrite the above inequality in the form

[TABLE]

Since $0<{\mathcal{H}}(\rho)<1$, the term $2-\nu-(2-3\nu){\mathcal{H}}(\rho)$ is larger than $2-\nu-(2-3\nu)=2\nu$ and the term $2+\nu-2\sqrt{2}\,(1-2\nu)+\left[3(2-\nu)-4\sqrt{2}(1-2\nu)\right]{\mathcal{H}}(\rho)$ larger than $2+\nu-2\sqrt{2}\,(1-2\nu)$. Hence for this inequality to be true, it suffices that the inequality

[TABLE]

be satisfied. For this to be true for all values of $\rho$, it suffices that

[TABLE]

This reasonable condition is assumed to be satisfied in the sequel.888It has already been noted by Leblond et al. (2011), in the case of a fracture energy independent of mode mixity, that the case of small positive values of $\nu$ is special. Then $\frac{\partial\Phi}{\partial\xi}(\rho;0)>0$ for $\rho<\rho^{\rm cr}$.

B.7 Consequences for the equation $\Phi(\rho^{0};\xi)=0$ for $\rho^{0}<\rho^{\rm cr}$

Let us assume that $\rho^{0}<\rho^{\rm cr}$. Then:

• •

At the point $\xi=0$, $\Phi(\rho^{0};0)>0$, see Subsection B.4.

• •

At the same point $\frac{\partial\Phi}{\partial\xi}(\rho^{0};0)>0$, see Subsection B.6.

• •

The function $\frac{\partial\Phi}{\partial\xi}(\rho^{0};\xi)$ increases with $\xi$, see Subsection B.4; hence $\frac{\partial\Phi}{\partial\xi}(\rho^{0};\xi)>0$ for $\xi>0$.

One deduces from these elements that when $\xi$ increases from [math] to $+\infty$, $\Phi(\rho^{0};\xi)$ increases from some positive value; hence the equation $\Phi(\rho^{0};\xi)=0$ on $\xi$ does not admit any solution in the interval $(0,+\infty)$.