# Out-of-plane deviation of a mode I+III crack encountering a tougher obstacle

Jean-Baptiste Leblond, Laurent Ponson

► To cite this version:

Jean-Baptiste Leblond, Laurent Ponson. Out-of-plane deviation of a mode I+III crack encountering a tougher obstacle. Comptes Rendus Mécanique, 2016, 344 (7), pp.521-531. 10.1016/j.crme.2016.04.003. hal-01340573

HAL Id: hal-01340573

https://hal.sorbonne-universite.fr/hal-01340573

Submitted on 1 Jul 2016

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Out-of-plane deviation of a mode I+III crack encountering a tougher obstacle

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

CrossMark

## Déviation hors plan d'une fissure chargée en mode I+III et rencontrant un obstacle plus tenace

Jean-Baptiste Leblond *, Laurent Ponson

Sorbonne Universités, UPMC Université Paris-6 et CNRS, UMR 7190, Institut Jean-Le-Rond-d'Alembert, 75005 Paris, France

## ARTICLE INFO

**Article history:**

Received 23 February 2016

Accepted 4 April 2016

Available online 5 May 2016

**Keywords:**

Crack

Mode I+III

Obstacle

Out-of-plane deviation

**Mots-clés :**

Fissure

Mode I+III

Obstacle

Déviation hors-plan

## ABSTRACT

One possible explanation of out-of-plane deviations of cracks loaded in mode I+III was suggested by Gao and Rice in 1986. These authors noted that small in-plane undulations of the crack front, arising from fluctuations of the fracture toughness, should generate a small local mode-II component, causing the crack to depart from planarity. Their analysis is completed here by explicitly calculating the evolution in time of the out-of-plane deviation of a mode-I+III crack encountering a tougher obstacle. The calculation is based on (i) first-order formulae for the stress intensity factors of a crack slightly perturbed within and out of its plane; and (ii) a "double" propagation criterion combining a Griffith condition on the local energy-release rate and a Goldstein-Salganik condition on the local stress intensity factor of mode II. It is predicted that the crack must evolve toward a stationary state, wherein the orthogonal distance from the average fracture plane to the perturbed crack front is constant outside the obstacle and varies linearly across it. We hope that this theoretical prediction will encourage comparison with experiments, and propose a fracture test involving propagation of a mode-I+III crack through a 3D-printed specimen containing some designed obstacle.

© 2016 Académie des sciences. Published by Elsevier Masson SAS. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

## RÉSUMÉ

Une explication possible de la déviation hors plan de fissures chargées en mode I+III a été suggérée par Gao et Rice en 1986. Ces auteurs ont remarqué que de petites ondulations coplanaires du front de fissure, dues à des fluctuations de la tenacité, doivent engendrer une petite composante locale de mode II, forçant ainsi la fissure à sortir de son plan. Leur analyse est complétée ici en calculant explicitement l'évolution au cours du temps de la déviation hors plan d'une fissure chargée en mode I+III et rencontrant un obstacle plus résistant. Le calcul repose sur (a) des formules donnant, au premier ordre, les facteurs d'intensité de contraintes d'une fissure légèrement perturbée dans son plan et en dehors et (b) un critère de propagation «double» combinant une condition de Griffith sur le taux de restitution d'énergie local et une condition de Goldstein-Salganik sur le facteur d'intensité local du mode II. On prédit que la fissure doit évoluer vers un état stationnaire, où la

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: jbl@lmm.jussieu.fr (J.-B. Leblond), ponson@dalembert.upmc.fr (L. Ponson).</footnote>
distance orthogonale du plan moyen de fracture au front de fissure perturbé est constante en dehors de l'obstacle et varie linéairement à l'intérieur de ce dernier. Nous espérons que cette prédiction théorique suscitera une comparaison avec des expériences, et proposons un essai de rupture impliquant la propagation d'une fissure chargée en mode I+III dans une éprouvette fabriquée grâce à une imprimante 3D, contenant un obstacle contrôlé.

© 2016 Académie des sciences. Published by Elsevier Masson SAS. This is an open access

article under the CC BY-NC-ND license

(http://creativecommons.org/licenses/by-nc-nd/4.0/).

# 1. Introduction

Out-of-plane deviations of cracks loaded in mixed mode I+III in elastic brittle materials are well documented experimentally. Such deviations are currently observed in a large variety of materials, as soon as the local stress field near the front of the initially planar crack involves a mode-III component.

Motivated by Pons and Karma's previous numerical study [1] of the instability of coplanar crack propagation in mode I+III, based on a "phase-field" model, Leblond et al. [2] recently performed a linear stability analysis of this type of propagation, based on the standard tools of linear elastic fracture mechanics. This analysis was based on two main elements: (i) first-order formulae due to Gao and Rice [3] and Movchan et al. [4] for the stress intensity factors (SIF) along the front of a crack slightly perturbed both within its plane and out of it; and (ii) a "double" propagation criterion, suggested by Pons and Karma's numerical study [1], enforced at all points of the crack front and at all instants, combining a Griffith condition of constant energy-release rate [5] and a Goldstein-Salganik condition of zero SIF of mode II [6]. The conclusion was that a bifurcation from coplanar to non-coplanar propagation, leading to an instability, exists for values of the ratio $K_{III}^{0}/K_{I}^{0}$ of the unperturbed SIF of modes III and I larger than some "threshold" depending on Poisson's ratio.

However, since this threshold is of the order of 0.5 for standard values of Poisson's ratio, this analysis cannot explain out-of-plane deviations commonly observed for smaller, or even much smaller values of the ratio $K_{III}^{0}/K_{I}^{0}$, see e.g. the recent work of Pham and Ravi-Chandar [7].

A plausible explanation of this discrepancy was proposed very recently by Chen et al. [8], using new phase-field simulations supplementing those of Pons and Karma [1]. These simulations confirmed the value of the theoretical threshold but evidenced a strongly subcritical character of the bifurcation, inaccessible to Leblond et al.'s linear stability analysis [2]. This subcritical character makes the theoretical threshold virtually unobservable in actual experiments.

Another possible explanation lies in the existence of "discontinuous" perturbations of the crack, in the form of disconnected tilted facets. Such facets have been observed experimentally many times, and appear to be a generic nonlinear feature of mixed-mode-I+III fracture (Pham and Ravi-Chandar [7]), again inaccessible to linear stability analysis.

Still another possible explanation (by no means incompatible with the previous two but rather complementary) could reside in a strong influence of imperfections upon the propagation path. An interesting suggestion in this direction was made as early as 1986 by Gao and Rice [3]. These authors pointed out that small, inevitable fluctuations of the fracture toughness within the crack plane are bound to create tiny in-plane undulations of the crack front; because of the coupling of modes II and III for planar crack problems, these undulations must generate a small mode-II component along this front, resulting in future out-of-plane deviations of the crack. The analysis was very recently supplemented by Leblond and Lazarus [9] by studying the stability of these deviations, using Cotterell and Rice's heuristic "directional stability" criterion [10]. It however remained incomplete, insofar as although out-of-plane deviations of the crack were predicted to be inevitable, their development in time was not studied.<fnref n="1" />

The aim of this work is precisely to complete Gao and Rice's analysis [3] by studying the evolution in time of out-of-plane deviations of a mode-I+III crack generated by small in-plane fluctuations of the fracture toughness. It is based on exactly the same elements as Leblond et al.'s bifurcation analysis [2], that is (i) first-order perturbation formulae for the SIF resulting from in-plane and out-of-plane perturbations of the crack, and (ii) a double propagation criterion combining Griffith and Goldstein-Salganik conditions. The only difference is that the fracture toughness is no longer considered to be constant within the initial crack plane. It is assumed instead that the initially planar crack suddenly encounters an obstacle of infinite elongation in the direction of propagation, with the same elastic properties but a slightly larger toughness than the rest of the material. This model obstacle geometry is considered for the sake of simplicity of the resulting analytical calculations, and to hopefully permit comparisons with future experimental observations.

Two calculations are performed. In a first step, the ratio $K_{III}^{0}/K_{I}^{0}$ of the unperturbed SIF of modes III and I (for the planar configuration of the crack) is assumed to be small. In this case the whole evolution in time of the out-of-plane deviations of the crack may be calculated. In a second step, arbitrary values of the same ratio are considered. It is then no longer possible to analytically determine the full evolution of out-of-plane deviations, but the stationary shape of the crack front may still

<footnote n="1">¹ Such a study was impossible at the time of publication of Gao and Rice's paper [3] since Movchan et al.'s formulae [4] for the SIF of a crack slightly perturbed out of its plane were not available yet.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Initial configuration of the crack.</figcaption></figure>

be calculated explicitly. In both cases it is found that the out-of-plane shape of the crack front evolves toward a very simple stationary state, wherein the orthogonal distance from the average fracture plane to the perturbed crack front is constant outside the obstacle, and varies linearly across it.

## 2. First-order expressions of the stress intensity factors for a mode-I+III crack encountering a slightly tougher obstacle

In the initial configuration, the geometry consists of an infinite body made of some isotropic elastic material, containing a planar, semi-infinite crack with a straight front (Fig. 1). Following the usual conventions, the Ox axis is chosen oriented along the direction of propagation of the crack, the Oy axis perpendicular to the crack plane and the Oz axis parallel to the crack front. The crack is loaded in mode I+III through some system of forces independent of the coordinate z, so that the SIF $K_I^0$, $K_{II}^0 = 0$, $K_{III}^0$ are independent of the position z along the crack front. (They may however depend upon its location within the crack plane.) The front is just about to penetrate into some obstacle, painted in grey in the figure, of width $2d$ (in the direction z) and infinite length (in the direction x), of elastic properties identical to those of the rest of the material but slightly larger fracture toughness. The Ox axis is placed outside the obstacle in the figure for legibility, but in the actual calculations will be relocated at the centre of the bounding surface of this obstacle, which will thus occupy the domain $x \ge 0, -\infty < y < +\infty, -d \le z \le d$.

In the present configuration, the crack front has penetrated into the obstacle by a mean distance x, resulting in both a small in-plane perturbation of the crack front, $\phi_x(x, z)$, and a small out-of-plane perturbation of the fracture surface, $\phi_y(x, z)$ (Fig. 2). The resulting perturbations of the p-th SIF ($p=1, II, III$) at the point z of the crack front are denoted $\delta_x K_p(x, z)$ and $\delta_y K_p(x, z)$, respectively. To first order in $(\phi_x, \phi_y)$, these perturbations are additive so that the total perturbation of the p-th SIF is simply $\delta K_p(x, z) = \delta_x K_p(x, z) + \delta_y K_p(x, z)$.

The perturbations of the SIF have been calculated by Gao and Rice [3] for the in-plane perturbation of the crack front, and Movchan et al. [4] for the out-of-plane perturbation of the crack surface. To apply their results, we introduce the reasonable assumption that the half-width d of the obstacle, which obviously governs the typical distance of variation of the perturbations $\phi_x(x, z)$ and $\phi_y(x, z)$, is much smaller than the characteristic length L – of the order of $K_I^0/(dK_I^0/dx)$ or $K_{III}^0/(dK_{III}^0/dx)$ – defined by the loading, in the absence of any length defined by the infinite geometry itself. This permits to retain only, in Gao and Rice's and Movchan, Gao and Willis's formulae [3,4], these terms involving the unperturbed SIF and disregard those involving the unperturbed non-singular stresses and similar higher-order constants.<fnref n="2" /> With this hypothesis Gao and Rice's formulae [3] for the $\delta_x K_p(x, z)$ read

$$ \left\{ \begin{aligned} \delta_x K_I(x, z) &= \frac{K_I^0}{2\pi} PV \int_{-\infty}^{+\infty} \frac{\phi_x(x, z') - \phi_x(x, z)}{(z' - z)^2} dz' \\ \delta_x K_{II}(x, z) &= -\frac{2}{2-\nu} K_{III}^0 \frac{\partial \phi_x}{\partial z}(x, z) \\ \delta_x K_{III}(x, z) &= \frac{2+\nu}{2-\nu} \frac{K_{III}^0}{2\pi} PV \int_{-\infty}^{+\infty} \frac{\phi_x(x, z') - \phi_x(x, z)}{(z' - z)^2} dz' \end{aligned} \right. \quad (1) $$

where $\nu$ denotes Poisson's ratio and the symbol PV a Cauchy principal value. Also, Movchan et al.'s formulae [4] for the $\delta_y K_p(x, z)$ read

<footnote n="2">² Because, K, T and φ denoting typical values of the SIF, the non-singular stresses and the crack perturbations, the former terms are of order $K\phi/d$ whereas the latter ones are, for those involving the non-singular stresses for instance, of order $T\phi/\sqrt{d} \sim K\phi/\sqrt{Ld} \ll K\phi/d$.</footnote>
<figure><img src="image_4.png" /><figcaption>(a) In-plane perturbation of the crack front.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Out-of-plane perturbation of the crack surface.</figcaption></figure>

<caption>**Fig. 2.** Present configuration of the crack, involving in-plane and out-of-plane perturbations of its front and surface.</caption>

$$
\left\{
\begin{aligned}
\delta_y K_1(x, z) &= -2K_{\text{III}}^0 \frac{\partial \phi_y}{\partial z}(x, z) + \delta_y K_1^{\text{skew}}(x, z) \\
\delta_y K_{\text{II}}(x, z) &= \frac{K_1^0}{2} \frac{\partial \phi_y}{\partial x}(x, z) - \frac{2-3\nu}{2-\nu} \frac{K_1^0}{2\pi} PV \int_{-\infty}^{+\infty} \frac{\phi_y(x, z') - \phi_y(x, z)}{(z' - z)^2} dz' \\
\delta_y K_{\text{III}}(x, z) &= \frac{2(1-\nu)^2}{2-\nu} K_1^0 \frac{\partial \phi_y}{\partial z}(x, z)
\end{aligned}
\right.
\qquad (2)
$$

In the first expression the quantity $\delta_y K_1^{\text{skew}}(x, z)$ – which is connected to Bueckner’s skew-symmetric crack-face weight functions [11], whence the notation – is given by<fnref n="3" />

$$
\delta_y K_1^{\text{skew}}(x, z) = \frac{\sqrt{2}}{4\pi} \frac{1 - 2\nu}{1 - \nu} K_{\text{III}}^0 \text{Re} \left[ \int_{0}^{x} dx' \int_{-\infty}^{+\infty} \frac{(\partial \phi_y / \partial z)(x', z')}{(x-x')^{1/2} [x-x' + i(z-z')]^{3/2}} dz' \right] \quad (3)
$$

3. In-plane and out-of-plane perturbations of the crack front and surface

3.1. *Hypotheses and notations*

We shall now determine the in-plane perturbation of the crack front and the out-of-plane perturbation of the crack surface resulting from penetration of the crack into the obstacle. These perturbations will be calculated as functions of (i) the assumedly small “normalized toughness contrast” $\epsilon$ defined by

$$
\epsilon = \frac{G_c^{\text{Obstacle}} - G_c^{\text{Matrix}}}{G_c^{\text{Matrix}}} \qquad (4)
$$

<footnote n="3">³ Movchan et al.’s expression [4] of $\delta_y K_1^{\text{skew}}(x, z)$, applicable to the sole case of a perturbation $\phi_y(x, z)$ independent of $x$, was extended by Leblond et al. [2] to the general case.</footnote>
where $G_{\text{matrix}}$ and $G_{\text{obstacle}}$ denote the critical energy-release rates of the matrix and obstacle, respectively, and (ii) the ratio $\rho$ of the initial SIF of modes III and I,

$$ \rho = \frac{K_{\text{III}}^{0}}{K_{I}^{0}} \qquad (5) $$

In a first step, this ratio will be assumed to be small; in a second one, it will be considered to be arbitrary.

Like in Leblond et al.'s work [2], we shall make use of a double propagation criterion consisting of

- Griffith's condition [5] $G(x, z) = G_c(z)$ where $G(x, z)$ denotes the local energy-release rate and $G_c(z)$ the local critical value of this rate – which depends only on $z$, see Fig. 1;

- Goldstein and Salganik's principle of local symmetry [6] stipulating that $K_{\text{II}}(x, z) = 0$ where $K_{\text{II}}(x, z)$ denotes the local SIF of mode II.

This double criterion will be enforced at all points of the crack front and all instants.<fnref n="4" />

### 3.2. Case of a small ratio of the initial mode III to mode I stress intensity factors

We shall show that when the ratio $\rho$ is small, the perturbations $\phi_x(x, z)$, $\phi_y(x, z)$ may be calculated in a “decoupled” way, one after another: $\phi_x(x, z)$ is determined independently of $\phi_y(x, z)$ from Griffith's condition $G(x, z) = G_c(z)$, and $\phi_y(x, z)$ is then deduced from $\phi_x(x, z)$ using Goldstein and Salganik's condition $K_{\text{II}}(x, z) = 0$.

We first note that the in-plane perturbation $\phi_x(x, z)$ of the crack front is generated by the contrast of toughness between the matrix and the obstacle but exists even in the absence of mode III, whereas the out-of-plane perturbation $\phi_y(x, z)$ of the crack surface arises from both the contrast of toughness and the presence of mode III. This implies that $\phi_x(x, z)$ and $\phi_y(x, z)$ must be of order $O(\epsilon)$ and $O(\epsilon\rho)$ respectively, so that $|\phi_y(x, z)| \ll |\phi_x(x, z)|$.

Using equations (1), (2) and (3), one then sees that:

- $\delta_x K_1(x, z)$ is of order $O(\epsilon)$, but both terms in the expression of $\delta_y K_1(x, z)$ are of order $O(\epsilon\rho^2)$, so that $|\delta_y K_1(x, z)| \ll |\delta_x K_1(x, z)|$ and $\delta K_1(x, z) \simeq \delta_x K_1(x, z)$ is globally of order $O(\epsilon)$;

- both $\delta_x K_{\text{III}}(x, z)$ and $\delta_y K_{\text{III}}(x, z)$ are of order $O(\epsilon\rho)$.

It follows that in the general expression of the perturbation $\delta G(x, z)$ of the energy-release rate,

$$ \delta G(x, z) = 2 \frac{1 - \nu^2}{E} K_I^0 \delta K_1(x, z) + 2 \frac{1 + \nu}{E} K_{\text{III}}^0 \delta K_{\text{III}}(x, z) \qquad (6) $$

(where $E$ denotes Young's modulus), the terms proportional to $K_I^0 \delta K_1(x, z)$ and $K_{\text{III}}^0 \delta K_{\text{III}}(x, z)$ are of order $O(\epsilon)$ and $O(\epsilon\rho^2)$ respectively; therefore

$$ \delta G(x, z) \approx 2 \frac{1 - \nu^2}{E} K_I^0 \delta K_1(x, z) \approx 2 \frac{1 - \nu^2}{E} K_I^0 \delta_x K_1(x, z) \qquad (7) $$

Thus the perturbation of the energy-release rate is, at the order considered, exactly the same as in the absence of mode III. This implies that the calculation of the in-plane perturbation $\phi_x(x, z)$, based on Griffith's condition $G(x, z) = G_c(z)$ enforced at all points of the crack front and all instants, is also exactly the same. This calculation is classical (see, e.g., Chopin et al. [12]) and will not be repeated here; the result for the difference $\phi_x(x, z) - \phi_x(x, 0)$, characterizing the in-plane shape of the deformed crack front, reads

$$ \phi_x(x, z) - \phi_x(x, 0) = \frac{\epsilon}{\pi} \left[ (z+d) \ln\left(|\frac{z}{d}+1|\right) - (z-d) \ln\left(|\frac{z}{d}-1|\right) \right] \qquad (8) $$

where the origin $O$ has been relocated at the centre of the bounding surface of the obstacle.

Next, one may calculate the out-of-plane perturbation $\phi_y(x, z)$ by enforcing Goldstein and Salganik's condition $K_{\text{II}}(x, z) = 0$ at all points of the crack front and all instants. By equations (1) and (2), this condition yields upon multiplication by $2/K_I^0$:

$$ \frac{\partial \phi_y}{\partial x}(x, z) - \frac{2 - 3\nu}{2 - \nu} \cdot \frac{1}{\pi} P V \int_{-\infty}^{+\infty} \frac{\phi_y(x, z') - \phi_y(x, z)}{(z' - z)^2} dz' = \frac{4}{2 - \nu} \rho \frac{\partial \phi_x}{\partial z}(x, z) \qquad (9) $$

<footnote n="4">⁴ This means disregarding the short initial transient period during which $G(x, z) < G_c(z)$ at some points, because the front has not fully penetrated into the obstacle yet. This is justified by the observation that the propagation distance needed to reach stationarity of the in-plane perturbation, being obviously of the order of $\epsilon d$, is much shorter than that needed to reach stationarity of the out-of-plane perturbation, of the order of $d$ as will be apparent below.</footnote>
This integro-differential equation on the function $\phi_y(x, z)$ may be transformed into an ordinary differential equation (ODE) by taking its Fourier transform in the direction $z$ of the crack front. The Fourier transform $\hat{\psi}(k)$ of an arbitrary function $\psi(z)$ is defined here by the equivalent formulae:

$$ \psi(z) = \int_{-\infty}^{+\infty} \hat{\psi}(k) e^{ikz} dk \quad \Leftrightarrow \quad \hat{\psi}(k) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} \psi(z) e^{-ikz} dz \quad (10) $$

With this definition, the Fourier transform of equation (9) reads

$$ \frac{\partial \hat{\phi}_y}{\partial x}(x, k) + \frac{2-3\nu}{2-\nu} |k| \hat{\phi}_y(x, k) = -i \frac{4\epsilon\rho}{2-\nu} \frac{\sin(kd)}{\pi|k|} \quad (11) $$

where account has been taken of the expression (8) of $\phi_x(x, z)$. For fixed $k$, equation (11) is an ODE on $\hat{\phi}_y(x, k)$, considered as a function of $x$. Integration of this ODE yields

$$ \hat{\phi}_y(x, k) = -i \frac{4\epsilon\rho}{2-3\nu} \frac{\sin(kd)}{\pi k^2} + C \exp\left(-\frac{2-3\nu}{2-\nu} |k|x\right) $$

where $C$ is a constant. The initial condition $\hat{\phi}_y(0, k) = 0$, resulting from the fact that the crack is unperturbed prior to encountering the obstacle, yields the value of this constant, and the final result for $\hat{\phi}_y(x, k)$ reads

$$ \hat{\phi}_y(x, k) = -i \frac{4\epsilon\rho}{2-3\nu} \frac{\sin(kd)}{\pi k^2} \left[ 1 - \exp\left(-\frac{2-3\nu}{2-\nu} |k|x\right) \right] \quad (12) $$

It remains to take the inverse Fourier transform of this function to get $\phi_y(x, z)$. This inverse Fourier transform is given by

$$ \phi_y(x, z) = \int_{-\infty}^{+\infty} \hat{\phi}_y(x, k) e^{ikz} dk \\ = \frac{8\epsilon\rho}{\pi(2-3\nu)} \int_{0}^{+\infty} \sin(kd) \sin(kz) \left[ 1 - \exp\left(-\frac{2-3\nu}{2-\nu}kx\right) \right] \frac{dk}{k^2} $$

where account has been taken of the fact that $\hat{\phi}_y(x, k)$ is an odd function of $k$. The integral here is of the type $\int_0^{+\infty} \sin(a\nu) \sin(b\nu)(1 - e^{-p\nu}) du/u^2$ (with $p \ge 0$); its calculation is explained in Appendix A, and the final result for $\phi_y(x, z)$ reads:

$$ \phi_y(x, z) = \frac{4}{\pi} \epsilon \rho \left\{ \frac{x}{2(2-\nu)} \ln \left[ \frac{(2-3\nu)^2 x^2 + (2-\nu)^2 (z+d)^2}{(2-3\nu)^2 x^2 + (2-\nu)^2 (z-d)^2} \right] + \frac{z+d}{2-3\nu} \arctan \left[ \frac{(2-3\nu)x}{(2-\nu)(z+d)} \right] - \frac{z-d}{2-3\nu} \arctan \left[ \frac{(2-3\nu)x}{(2-\nu)(z-d)} \right] \right\} \quad (13) $$

The stationary state of the crack surface, obtained after propagation over a very large distance, is of special interest; this state corresponds to the limit $x \to +\infty$ in the preceding expression, $z$ being fixed. In this limit the term $x\ln[\dots]$ goes to zero and the terms $\arctan[\dots]$ go to $+\infty$ or $-\infty$ depending on the position of $z$ with respect to $d$ and $-d$; the result is

$$ \phi_y(+\infty, z) = \begin{cases} -\frac{4\epsilon\rho}{2-3\nu} d & \text{if } z \le -d \\ \frac{4\epsilon\rho}{2-3\nu} z & \text{if } -d \le z \le d \\ \frac{4\epsilon\rho}{2-3\nu} d & \text{if } z \ge d. \end{cases} \quad (14) $$

Therefore, in the stationary state, the orthogonal distance $\phi_y(x, z)$ from the original crack plane to the perturbed crack front is constant and equal to $\pm \frac{4\epsilon\rho}{2-3\nu} d$ outside the obstacle, and varies linearly with $z$ inside it. The angle of rotation of the front about the $Ox$ axis within the obstacle is thus $-\frac{4\epsilon\rho}{2-3\nu}$; its sign is opposite to that of the initial SIF $K_{III}^0$ of mode III (for a positive contrast $\epsilon$).
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Normalized in-plane perturbation of the crack front $[\phi_x(x, z) - \phi_x(x, 0)]/(\epsilon d)$ – Case of a small ratio $\rho = K_{\text{III}}^0/K_1^0$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Normalized out-of-plane perturbation of the crack surface $\phi_y(x, z)/(\epsilon\rho d)$ – $\nu = 0.3$ – Case of a small ratio $\rho = K_{\text{III}}^0/K_1^0$.</figcaption></figure>

Figs. 3, 4 and 5 illustrate the results obtained:

- Fig. 3 shows the in-plane perturbation of the crack front, measured in reference to the centre of the obstacle and normalized by the toughness contrast $\epsilon$ and the half-width $d$ of this obstacle, $[\phi_x(x, z) - \phi_x(x, 0)]/(\epsilon d)$, as a function of $z/d$. Note that this perturbation is independent of the distance of propagation $x$, and that the centre of the obstacle is (quite naturally) the most retarded point of the crack front.

- Fig. 4 shows, for the value $\nu = 0.3$ of Poisson's ratio, the out-of-plane perturbation of the crack surface, normalized by the toughness contrast $\epsilon$, the ratio $\rho$ of the unperturbed mode III to mode I SIF and the half-width $d$ of the obstacle, $\phi_y(x, z)/(\epsilon\rho d)$; this perturbation is plotted as a function of $z/d$ for various values of $x/d$. Note that at each point $z$ of the crack front, when the distance of propagation $x$ increases from zero to infinity, the crack gradually evolves from an initial flat configuration to a stationary non-planar one obeying equation (14). Obvious dimensional considerations dictate that the typical distance of propagation required to reach this stationary configuration must be of the form $d\Phi(z/d)$ where $\Phi$ is a certain function, so for every $z/d$ it scales with $d$ (the only characteristic distance in the problem). But clearly the larger the value of $z/d$, the slower the evolution toward the stationary state; for a fixed distance $x$ of propagation, at a sufficiently large distance $z$ from the centre of the obstacle, the crack goes back to its original plane.

- Finally Fig. 5 shows a full 3D view of the developing crack surface. The crack propagates toward the upper left corner of the figure.

### 3.3. Case of an arbitrary ratio of the initial mode III to mode I stress intensity factors

When the ratio $\rho$ of the initial SIF of modes III and I takes arbitrary values, it becomes impossible to calculate analytically the full evolution in time of the perturbations $\phi_x(x, z)$, $\phi_y(x, z)$ of the crack front and surface. The reason is that the term $\delta_y K_1^{\text{skew}}(x, z)$ in the expression (3) of $\delta_y K_1(x, z)$ is no longer negligible. This term creates a very complex dependence of the perturbed SIF of mode I upon the past history of the crack, which leads to intractable calculations: $\phi_y(x, z)$ may be expressed simply as a function of $\phi_x(x, z)$ using Goldstein and Salganik's condition, but Griffith's condition then yields an evolution equation for $\phi_x(x, z)$ which is of integrodifferential type and remains so even after Fourier transform in the direction of the crack front.

However, we shall show that it remains possible to calculate analytically the stationary state of the crack, corresponding to the limit of a very large distance of propagation, $x \to +\infty$. The reason is that in this limit the perturbation $\phi_y(x', z')$ may
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Geometry of the fracture surface, defined by the sequence of successive crack fronts – $v = 0.3$ – Case of a small ratio $\rho = K_{\text{III}}^0 / K_I^0$. The x-axis is represented in a logarithmic scale since a long distance of propagation is necessary to reach stationarity of the out-of-plane perturbation.</figcaption></figure>

be considered as independent of $x'$ in the expression (3) of $\delta_y K_I^{\text{skew}}(x, z)$, because the contribution to the integral of the region where this function varies with $x'$ becomes negligibly small. The quantity $\delta_y K_I^{\text{skew}}(x, z)$ then takes the following very simple form (Movchan et al. [4]):

$$ \delta_y K_I^{\text{skew}}(z) = \frac{1-2v}{\sqrt{2}(1-v)} K_{\text{III}}^0 \frac{\partial \phi_y}{\partial z}(z) \quad (15) $$

where unnecessary indications of dependence upon $x$ have been discarded.

Taking the Fourier transforms of equations (1), (2) and (15) and expressing $K_{\text{III}}^0$ as $\rho K_I^0$, one gets

$$ \left\{ \begin{aligned} \widehat{\delta K}_I(k) &= K_I^0 \left[ -\frac{|k|}{2} \hat{\phi}_x(k) + i \left( -2 + \frac{1-2v}{\sqrt{2}(1-v)} \right) \rho k \hat{\phi}_y(k) \right] \\ \widehat{\delta K}_{II}(k) &= K_I^0 \left[ -i \frac{2}{2-v} \rho k \hat{\phi}_x(k) + \frac{2-3v}{2(2-v)} |k| \hat{\phi}_y(k) \right] \\ \widehat{\delta K}_{III}(k) &= K_I^0 \left[ -\frac{2+v}{2(2-v)} \rho |k| \hat{\phi}_x(k) + 2i \frac{(1-v)^2}{2-v} k \hat{\phi}_y(k) \right] \end{aligned} \right. \quad (16) $$

From there and equation (6) follows the expression of the Fourier transform $\widehat{\delta G}(k)$ of the perturbation of the energy-release rate $\delta G(z)$:

$$ \widehat{\delta G}(k) = \frac{G^0}{1 + \frac{\rho^2}{1-\nu}} \left\{ - \left[ 1 + \frac{2+\nu}{(1-\nu)(2-\nu)} \rho^2 \right] |k| \hat{\phi}_x(k) + 2i \left[ -\frac{2}{2-\nu} + \frac{1-2\nu}{\sqrt{2}(1-\nu)} \right] \rho k \hat{\phi}_y(k) \right\} \quad (17) $$

where

$$ G^0 = \frac{1 - \nu^2}{E} (K_I^0)^2 + \frac{1 + \nu}{E} (K_{\text{III}}^0)^2 = \frac{1 - \nu^2}{E} (K_I^0)^2 \left( 1 + \frac{\rho^2}{1 - \nu} \right) \quad (18) $$

is the unperturbed energy-release rate.

Then Goldstein and Salganik's condition written in "Fourier's form", $\widehat{\delta K}_{II}(k) = 0$, yields

$$ \hat{\phi}_y(k) = i \frac{4}{2 - 3\nu} \rho \operatorname{sgn}(k) \hat{\phi}_x(k) \quad (19) $$

where $\operatorname{sgn}(k)$ denotes the sign of $k$. Also, Griffith's condition reads $G(z) = G^0 + \delta G(z) = G_c(z) = G_c^{\text{Matrix}} + \delta G_c(z)$ where $\delta G_c(z)$ is equal to zero in the matrix and to $G_c^{\text{Obstacle}} - G_c^{\text{Matrix}} = \epsilon G_c^{\text{Matrix}}$ in the obstacle; identifying terms of order 0 and 1 in $\epsilon$, one gets:

- at order 0: $G^0 = G_c^{\text{Matrix}}$;

- at order 1: $\delta G(z) = \delta G_c(z)$, or equivalently in Fourier's form, $\widehat{\delta G}(k) = \widehat{\delta G}_c(k) = \epsilon G_c^{\text{Matrix}} \frac{\sin(kd)}{\pi k}$.

The last condition yields upon use of equations (17) and (19), after some algebraic manipulations:

$$ \hat{\phi}_x(k) = -\epsilon A \frac{\sin(kd)}{\pi k |k|} \quad \text{where} \quad A = \frac{1 + \frac{\rho^2}{1-\nu}}{1 - \frac{3(2-\nu) - 4\sqrt{2}(1-2\nu)}{(1-\nu)(2-3\nu)} \rho^2} \quad (20) $$
The case envisaged previously of a small $\rho$ corresponds to a value of A of unity. Equation (20) shows that the effect of a finite $\rho$ is simply to multiply the stationary value of $\hat{\phi}_x(k)$ by the “amplification factor” A, which is independent of k. It follows that the stationary value of the in-plane perturbation $\phi_x(z)$ of the crack front is also obtained from that corresponding to a small $\rho$, equation (8), through multiplication by the same factor<fnref n="5" />:

$$ \phi_x(z) - \phi_x(0) = \frac{\epsilon A}{\pi} \left[ (z+d) \ln\left(\left|\frac{z}{d} + 1\right|\right) - (z-d) \ln\left(\left|\frac{z}{d} - 1\right|\right) \right] \quad (21) $$

Furthermore, the relation (19) between the stationary values of $\hat{\phi}_x(k)$ and $\hat{\phi}_y(k)$ is the same no matter whether $\rho$ is small or arbitrary. Therefore the stationary value of the out-of-plane perturbation $\phi_y(z)$ of the crack surface is also obtained from that corresponding to a small $\rho$, equation (14), through multiplication by the amplification factor A:

$$ \phi_y(z) = \begin{cases} \frac{4\epsilon\rho A}{2 - 3\nu} d & \text{if } z \le -d \\ \frac{4\epsilon\rho A}{2 - 3\nu} z & \text{if } -d \le z \le d \\ \frac{4\epsilon\rho A}{2 - 3\nu} d & \text{if } z \ge d \end{cases} \quad (22) $$

Equation (22) means that in the stationary state, the angle of rotation of the front about the Ox axis within the obstacle is $-\frac{4\epsilon\rho A}{2-3\nu}$. By the expression (20) of the amplification factor A, this stationary angle diverges to infinity when $\rho$ approaches the value $\rho_c$ given by

$$ \rho_c = \sqrt{\frac{(1 - \nu)(2 - 3\nu)}{3(2 - \nu) - 4\sqrt{2}(1 - 2\nu)}} \quad (23) $$

This “critical” value of the ratio $\rho = K_{III}^0/K_I^0$ is interestingly exactly the same as that derived by Leblond et al. [2] for the onset of a bifurcation leading to an instability in a homogeneous material (no obstacle).<fnref n="6" />

The conclusion that the stationary angle of rotation of the front about the Ox axis within the obstacle should diverge when $\rho$ approaches $\rho_c$, is to be taken with caution for two reasons: (i) because of the issues raised by the unfavorable comparison of the theoretical and experimental values of $\rho_c$ (see the Introduction); (ii) because the perturbation approach employed breaks down anyway when the in-plane and out-of-plane perturbations of the crack front and surface become large. It is however probably reasonable to retain the qualitative conclusion that this stationary angle of rotation should increase more quickly than the ratio of the initial SIF of modes III and I.

## 4. Summary and perspectives

The aim of this work was to study the evolution in time of the out-of-plane deviation of a crack loaded in mode I+III, and suddenly encountering a slightly more resistant obstacle of infinite length in the direction of propagation. Two calculations were performed. First, the ratio of the unperturbed SIF of modes III and I (for the planar configuration of the crack) was assumed to be small; the whole evolution in time of the out-of-plane deviation of the crack was then determined. Second, arbitrary values of this ratio were considered; the sole stationary state of the propagating crack front could then be calculated. In both cases it was found that in this stationary state, the orthogonal distance from the average fracture plane to the perturbed crack front is constant outside the obstacle and varies linearly across it.

It would be instructive to compare these theoretical predictions to actual experimental observations. In order to do so, it is proposed to use a multimaterial 3D printer to design specimens containing obstacles of controlled geometry and toughness. Mode-I+III loading conditions may be achieved through 3-point or 4-point bending of beams with an initial crack inclined at some angle over the cross section. Various difficulties may be anticipated, some technological, some of more basic nature:

- it will be necessary to use, for the matrix and the obstacle, materials obeying strong, maybe conflicting constraints: (i) be compatible with present-day technology of multimaterial 3D printers; (ii) possess similar or nearly similar elastic properties; (iii) possess different toughness values;

- the requirement of an *inclined* initial crack will make it difficult to create it by fatigue. It will probably be necessary to generate it through 3D printing, but what will then be obtained will rather be a notch of small but finite root radius than a real crack;

<footnote n="5">⁵ Of course, the same conclusion can also be reached by directly taking the inverse Fourier transform of equation (20).</footnote>

<footnote n="6">⁶ This coincidence arises from the fact that the divergence of A is due to the vanishing of $\delta\tilde{G}(k)$ in the limit $\rho \to \rho_c$; this vanishing is independent of whether $\delta\tilde{G}(k)$ is equated to 0, as in Leblond et al.'s previous work [2], or to $\delta\tilde{G}_c(k)$, as in the present one.</footnote>
- most importantly, the presence in the loading of a mode-III component will inevitably result in the creation of fracture facets inclined over the mean surface of the crack. This will complicate the observation of the rotation of the crack front about the direction of propagation induced by the obstacle. An interesting observation in this respect is that the sign of the rotation of the front due to the obstacle is predicted to be opposite to that of the rotation of the fracture facets.<fnref n="7" /> This feature should facilitate the distinction between the two rotations.

It is hoped that in spite of these foreseeable difficulties, experimental assessment of the deformed crack front shapes predicted theoretically will be possible.

The basic aim of this study was to confirm and complete Gao and Rice [3]'s analysis of the influence of fluctuations of the fracture toughness upon the instability of coplanar crack propagation in mode I+III, which essentially expresses itself through the formation of fracture facets. One might object that this aim has not really been fulfilled for two reasons:

- the analysis does not shed light on the formation of fracture facets, since the obstacle considered has been shown to induce a rotation of the crack front in the “wrong” direction, opposite to that of the facets;

- the geometry considered is too idealized, since it involves a single obstacle instead of a more or less random distribution of zones with different toughness values.

The answer to the first objection is that the reason for considering a region of slighter higher toughness was only that experimental verification of the theoretical predictions should be easier in such a case, as explained above. But considering a region of slightly lower toughness only requires to change the sign of the normalized toughness contrast $\epsilon$ in the above formulae. For such a region the crack front is predicted to rotate in the same direction as the commonly observed facets.

The answer to the second objection will be fully provided in a future paper considering a periodic (thus still idealized but more realistic), two-step distribution of fracture toughness.

## Appendix A

To calculate the integral

$$I(a, b; p) = \int_{0}^{+\infty} \sin(au) \sin(bu)(1 - e^{-pu}) \frac{du}{u^2} \quad (p \ge 0)$$

the simplest solution is probably to differentiate with respect to $p$ and get

$$\frac{\partial I}{\partial p}(a, b; p) = \int_{0}^{+\infty} \sin(au) \sin(bu) e^{-pu} \frac{du}{u} = \frac{1}{4} \ln \left[ \frac{p^2 + (a+b)^2}{p^2 + (a-b)^2} \right]$$

where Gradsteyn and Ryzhik's formula (3.947.1) [13] has been used; and then to integrate with respect to $p$, accounting for the fact that $I(a, b; 0) = 0$. One thus gets

$$I(a, b; p) = \frac{p}{4} \ln \left[ \frac{p^2 + (a+b)^2}{p^2 + (a-b)^2} \right] + \frac{a+b}{2} \arctan \left( \frac{p}{a+b} \right) - \frac{a-b}{2} \arctan \left( \frac{p}{a-b} \right)$$

**Remark.** Gradsteyn and Ryzhik's formula (3.947.2) [13] directly provides the value of the integral looked for, but this formula contains an error, as can be seen by considering the special case $p = 0$. In contrast formula (3.947.1) has been checked to be correct.

## References

[1] A.J. Pons, A. Karma, Helical crack-front instability in mixed-mode fracture, Nature 464 (2010) 85–89.

[2] J.B. Leblond, A. Karma, V. Lazarus, Theoretical analysis of crack front instability in mode I+III, J. Mech. Phys. Solids 59 (2011) 1872–1887.

[3] H. Gao, J.R. Rice, Shear stress intensity factors for planar crack with slightly curved front, ASME J. Appl. Mech. 53 (1986) 774–778.

[4] A.B. Movchan, H. Gao, J.R. Willis, On perturbations of plane cracks, Int. J. Solids Struct. 35 (1998) 3419–3453.

[5] A.A. Griffith, The phenomena of rupture and flow in solids, Philos. Trans. R. Soc. Lond. Ser. A 221 (1920) 163–198.

[6] R.V. Goldstein, R.L. Salganik, Brittle fracture of solids with arbitrary cracks, Int. J. Fract. 10 (1974) 507–523.

[7] K.H. Pham, K. Ravi-Chandar, Further examination of the criterion for crack initiation under mixed-mode I+III loading, Int. J. Fract. 189 (2014) 121–138.

[8] C.H. Chen, T. Cambonie, V. Lazarus, M. Nicoli, A. Pons, A. Karma, Crack front segmentation and facet coarsening in mixed-mode fracture, Phys. Rev. Lett. 115 (2015) 265503, 1–5.

<footnote n="7">⁷ The sign of the rotation of the front due to the obstacle has been noted above to be opposite to that of the initial SIF of mode III, whereas the sign of the rotation of the facets is experimentally known to be identical.</footnote>
[9] J.B. Leblond, V. Lazarus, On the strong influence of imperfections upon the quick deviation of a mode I+III crack from coplanarity. Special Issue in Memoriam: Huy Dong Bui, J. Mech. Mater. Struct. 10 (2015) 299–316.

[10] B. Cotterell, J.R. Rice, Slightly curved or kinked cracks, Int. J. Fract. 16 (1980) 155–169.

[11] H.F. Bueckner, Weight functions and fundamental fields for the penny-shaped and the half-plane crack in three-space, Int. J. Solids Struct. 23 (1987) 57–93.

[12] J. Chopin, A. Prevost, A. Boudaoud, M. Adda-Bedia, Crack front dynamics across a single heterogeneity, Phys. Rev. Lett. 107 (2011) 144301.

[13] I.S. Gradshteyn, I.M. Ryzhik, Table of Integrals, Series, and Products, Academic Press, New York, 1980.