# Pinning of crack fronts by hard and soft inclusions: A phase field study

Hervé Henry

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

## To cite this version:

Hervé Henry. Pinning of crack fronts by hard and soft inclusions: A phase field study. Physical Review E, 2024, 109 (2), pp.025002. 10.1103/PhysRevE.109.025002. hal-04609803

HAL Id: hal-04609803

https://hal.science/hal-04609803

Submitted on 12 Jun 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Pinning of crack fronts by hard and soft inclusions: A phase field study

Hervé Henry

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

## To cite this version:

Hervé Henry. Pinning of crack fronts by hard and soft inclusions: A phase field study. Physical Review E, 2024, 109 (2), pp.025002. 10.1103/PhysRevE.109.025002. hal-04609803

HAL Id: hal-04609803

https://hal.science/hal-04609803v1

Submitted on 12 Jun 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Pinning of crack fronts by hard and soft inclusions: a phase field study

Hervé Henry

Laboratoire PMC, Ecole Polytechnique, CNRS, IP Paris, 91120 Palaiseau Cedex

Through tridimensional numerical simulations of crack propagating in material with an elastic moduli heterogeneity it is shown that the presence of a simple inclusion can affect dramatically the propagation of the crack. Both the presence of soft and hard inclusions can lead to the arrest of a crack front. Here the mechanism leading to the arrest of the crack are described and shown to depend on the nature of the inclusion. This is also the case in regimes where the presence of the inclusion leads to a slow down of the crack.

## I. INTRODUCTION

While most theoretical work on fracture propagation has been devoted to the study of fracture propagation in homogeneous materials, brittle or ductile, most materials that are actually used are heterogeneous. This is the case at the macroscopic scale for construction materials such as stones or concrete. This is also true for materials such as nacre or polymer blends that are homogeneous at a macroscopic scale but present heterogeneities at the microscopic scale. In all cases, the heterogeneity of the microstructure leads to heterogeneities of both fracture energy and elastic moduli. When considering the elastic properties of intact materials, the later are the only one at play and homogeneization techniques have allowed to propose effective material properties, even in cases where the phase constitutive of the material have dramatically different properties as it is the case for instance with porous materials[46-48].

However, in the context of fracture the process zone is most of the time at a scale that is smaller than a typical representative volume element. As a result contrarily to what happens with elastic properties, the use of effective material properties when describing crack propagation cannot be satisfactory. This is exemplified by Nacre[49-51]. Hence a proper description of the crack front in an heterogeneous material is first needed. In this case the possible effects of the microstructure are many. For instance, heterogeneities, by modifying the stress field locally can affect crack front trajectory. Fracture toughness heterogeneities can lead to crack deflection or arrest and to the propagation of cracks along weak interfaces.

In this context many studies have been devoted to the effects of regions with different fracture energy and the presence of interfaces[52-58]. The effect of a misfitting inclusion (i.e. with eigenstrain) on a crack has also been studied[59-61]. However little attention has been paid on changes elastic moduli which can vary by order magnitude and can strongly affect crack propagation as is evidenced by the deflection of cracks toward voids[62, 63] and deflection by hard inclusions[64, 65]. Recently, some progress on such materials has been made through numerical simulations using phase field models[66-68]. The focus of these work was on the nature of the inclusion itself and on its elastic properties. For instance in [67] the inclusion is lying exactly in the plane of propagation

of the crack front and the effects of both fracture energy and material stiffness are considered. Other studies were focused on 2D geometries while the bridging mechanism discussed in[61] and the results of [57] show clearly that 3D effects must be taken into account.

Here the possible effects of simple heterogeneities on crack propagation are shown and the emphasis is put on three dimensional effects since the microstructure of most materials are tridimensional. More specifically this work focuses on the effects of the position of the inclusion with respect to the crack propagation plane. In addition, since in composite materials, the toughening effects of inclusions are often sought in the context of catastrophic failure where the dynamic aspects of fracture propagation cannot be overlooked[69], this work brings some light on these aspects while being limited to non branching cracks[70].

To this purpose, a phase field model is used and the effect of an inclusion with a different elastic moduli on a crack front is studied. Simulations results indicate that both hard and soft inclusions affect crack propagation and can eventually stop a crack front that is called pinning here. However, the way the crack front is stopped by the inclusion differ. In addition in all cases, if the stress is increased the crack front can resume normal propagation without any additional crack nucleation as would be needed in two dimensional systems. The mechanisms that lead to this depinning are found to depend strongly on the nature of the inclusion:hard with higher elastic moduli or soft with lower elastic moduli.

## II. MODEL, SYSTEM DESCRIPTION AND NUMERICS

The system considered here consists of a simple spherical inclusion embedded in a matrix where a crack is propagating along the *x* axis with an applied strain along the *y* axis as schematized in fig.1.This would correspond to a void or a softer phase in an homogeneous material such as glass or, for instance, a tougher/harder phase in a polymeric material such as carbon black. For the sake of simplicity and to avoid to have to deal with crack front curvature[71, 72], the boundary conditions along the crack front direction, denoted *z*, is periodic (in all simulations presented here, the period is 30 space units
# Pinning of crack fronts by hard and soft inclusions: a phase field study

Hervé Henry

Laboratoire PMC, Ecole Polytechnique, CNRS, IP Paris, 91120 Palaiseau Cedex

Through tridimensional numerical simulations of crack propagating in material with an elastic moduli heterogeneity it is shown that the presence of a simple inclusion can affect dramatically the propagation of the crack. Both the presence of soft and hard inclusions can lead to the arrest of a crack front. Here the mechanism leading to the arrest of the crack are described and shown to depend on the nature of the inclusion. This is also the case in regimes where the presence of the inclusion leads to a slow down of the crack.

## I. INTRODUCTION

While most theoretical work on fracture propagation has been devoted to the study of fracture propagation in homogeneous materials, brittle or ductile, most materials that are actually used are heterogeneous. This is the case at the macroscopic scale for construction materials such as stones or concrete. This is also true for materials such as nacre or polymer blends that are homogeneous at a macroscopic scale but present heterogeneities at the microscopic scale. In all cases, the heterogeneity of the microstructure leads to heterogeneities of both fracture energy and elastic moduli. When considering the elastic properties of intact materials, the later are the only one at play and homogeneization techniques have allowed to propose effective material properties, even in cases where the phase constitutive of the material have dramatically different properties as it is the case for instance with porous materials[46-48].

However, in the context of fracture the process zone is most of the time at a scale that is smaller than a typical representative volume element. As a result contrarily to what happens with elastic properties, the use of effective material properties when describing crack propagation cannot be satisfactory. This is exemplified by Nacre[49-51]. Hence a proper description of the crack front in an heterogeneous material is first needed. In this case the possible effects of the microstructure are many. For instance, heterogeneities, by modifying the stress field locally can affect crack front trajectory. Fracture toughness heterogeneities can lead to crack deflection or arrest and to the propagation of cracks along weak interfaces.

In this context many studies have been devoted to the effects of regions with different fracture energy and the presence of interfaces[52-58]. The effect of a misfitting inclusion (i.e. with eigenstrain) on a crack has also been studied[59-61]. However little attention has been paid on changes elastic moduli which can vary by order magnitude and can strongly affect crack propagation as is evidenced by the deflection of cracks toward voids[62, 63] and deflection by hard inclusions[64, 65]. Recently, some progress on such materials has been made through numerical simulations using phase field models[66-68]. The focus of these work was on the nature of the inclusion itself and on its elastic properties. For instance in [67] the inclusion is lying exactly in the plane of propagation

of the crack front and the effects of both fracture energy and material stiffness are considered. Other studies were focused on 2D geometries while the bridging mechanism discussed in[61] and the results of [57] show clearly that 3D effects must be taken into account.

Here the possible effects of simple heterogeneities on crack propagation are shown and the emphasis is put on three dimensional effects since the microstructure of most materials are tridimensional. More specifically this work focuses on the effects of the position of the inclusion with respect to the crack propagation plane. In addition, since in composite materials, the toughening effects of inclusions are often sought in the context of catastrophic failure where the dynamic aspects of fracture propagation cannot be overlooked[69], this work brings some light on these aspects while being limited to non branching cracks[70].

To this purpose, a phase field model is used and the effect of an inclusion with a different elastic moduli on a crack front is studied. Simulations results indicate that both hard and soft inclusions affect crack propagation and can eventually stop a crack front that is called pinning here. However, the way the crack front is stopped by the inclusion differ. In addition in all cases, if the stress is increased the crack front can resume normal propagation without any additional crack nucleation as would be needed in two dimensional systems. The mechanisms that lead to this depinning are found to depend strongly on the nature of the inclusion:hard with higher elastic moduli or soft with lower elastic moduli.

## II. MODEL, SYSTEM DESCRIPTION AND NUMERICS

The system considered here consists of a simple spherical inclusion embedded in a matrix where a crack is propagating along the *x* axis with an applied strain along the *y* axis as schematized in fig.1.This would correspond to a void or a softer phase in an homogeneous material such as glass or, for instance, a tougher/harder phase in a polymeric material such as carbon black. For the sake of simplicity and to avoid to have to deal with crack front curvature[71, 72], the boundary conditions along the crack front direction, denoted *z*, is periodic (in all simulations presented here, the period is 30 space units
<figure><img src="image_3.png" /><figcaption>FIG. 1. Perspective view of the simulation setup. The crack is propagating from left to right and the inclusion is the sphere seen on the r.h.s. of the simulation domain.</figcaption></figure>

(SU), while the size of the sample in the x and y directions is 480 SU and 240 SU). The loading is imposed through the use of fixed displacement $\Delta Y$ at top and bottom boundaries along the y axis in an infinite strip geometry. In the context of dynamic crack propagation

(see pages 243-247 of [69]), this leads to a crack propagating with a steady velocity that is function of the available elastic energy denoted $G = (0.5\lambda + \mu)(\Delta Y/L)^2$.

Hence the system geometry is given by its length along the y axis denoted $L$, its width (or period) along the z axis $W$. The radius of the spherical inclusion is denoted $R$ and the distance between its center and the midplane of the system is denoted $b$. When $R$ is larger than $b$ the inclusion intersects the midplane that corresponds to the plane of steady crack propagation. The crack is initiated by the introduction of a pre crack that corresponds to a spatial distribution of the phase field with a thickness of the same order of magnitude as the regularization length. It must be noted that the crack propagates freely over a significant distance before its path is affected by the inclusion. Simulations performed using twice smaller values of $L$, keeping $W = 30$ lead to quantitatively similar results.

To model crack propagation a phase field model is used. It relies on the introduction of an additional variable the phase field $\varphi$ that can be seen as $1-d$ where $d$ is a damage variable. The evolution of the phase field derives from a free energy functional that writes:

$$ \mathcal{F} = \int \left( (hV(\varphi) - \varepsilon_c g(\varphi) + \frac{D}{2} |\nabla\varphi|^2) + g(\varphi) \mathcal{E}_{el} \right) dV \quad (1) $$

$$ \mathcal{E}_{el} = \frac{\lambda(\mathbf{x})}{2}(\mathrm{tr}\epsilon)^2 + \mu(\mathbf{x})\mathrm{tr}(\epsilon^2) \quad (2) $$

where $h = 0.25$, $D = 4$. and $\varepsilon_c = 0.25$ are model parameters and $\mathcal{E}_{el}$ is the elastic energy density for an intact material as a function of the strain tensor $\epsilon$. Its expression involves the space dependant coefficients $\lambda(\mathbf{x})$ and $\mu(\mathbf{x})$ that encode the material heterogeneity as will be described below. $V(\varphi) - \varepsilon_c g(\varphi)$ is a tilted double well potential with $V = \varphi^2(1 - \varphi)^2$ and $g(\varphi) = 4\varphi^3 - 3\varphi^4$ is a coupling function chosen so that when $\varphi = 1$ the contribution of the elastic strain corresponds to the intact material while when $\varphi = 0$, $g = 0$ and the contribution of strain to the energy of the system is 0. It has been shown in [73] that with this model the equilibrium solutions in 1D (local minima of $\mathcal{F}$) can be either $\varphi = 1$ and uniform strain or localized solutions where $\varphi$ is equal to 1 except in a small region of thickness $w_\varphi$, the so called regularization length, where it goes toward small values ($\ll 1$). If, as it was the case here, the shape of the potential (i.e. the ratio $h/\varepsilon_c$) is kept constant the regularization length scales like $D/h$. With the values considered here, $w_\varphi$ is of the order of 6. SU if $w_\varphi$ is the thickness of the region where $\varphi < 0.5$. In this case the strain is concentrated in this region and the stress in the material is vanishing[73]. The corresponding fracture energy is $2\Gamma = 1.046$ energy units and will be called in the following $G_c$. This model does not present any strain softening[74]

and has the peculiarity that contrarily to what is observed in other approaches such as the Francfort-Marigo variational approach[75-77], there is no crack nucleation. Hence this model is limited to the modelisation of crack propagation. This limitation comes with the advantage that no unphysical crack nucleation can occur when running macroscopic simulations and using a model interface larger than the microscopic process zone. Indeed it has been shown that if the regularisation length used in the model is chosen larger than the actual size of the process zone, then the material stress at failure is reduced proportionally, leading to the possibility of unphysical nucleation events at low stress values. [78, 79].

In order to describe the evolution of a crack the relaxation evolution equation<fnref n="1" /> that derive from this free energy is:

$$ \partial_t \varphi = -\frac{1}{\beta} \frac{\delta \mathcal{F}}{\delta \varphi} \quad (3) $$

<footnote n="1">¹ More complex kinetic equations[80] may be used, especially for fast crack propagation, however in the context of crack propagation the complexity of the events taking place in the process zone would not allow to estimate them.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Perspective view of the simulation setup. The crack is propagating from left to right and the inclusion is the sphere seen on the r.h.s. of the simulation domain.</figcaption></figure>

(SU), while the size of the sample in the x and y directions is 480 SU and 240 SU). The loading is imposed through the use of fixed displacement $\Delta Y$ at top and bottom boundaries along the y axis in an infinite strip geometry. In the context of dynamic crack propagation

(see pages 243-247 of [69]), this leads to a crack propagating with a steady velocity that is function of the available elastic energy denoted $G = (0.5\lambda + \mu)(\Delta Y/L)^2$.

Hence the system geometry is given by its length along the y axis denoted $L$, its width (or period) along the z axis $W$. The radius of the spherical inclusion is denoted $R$ and the distance between its center and the midplane of the system is denoted $b$. When $R$ is larger than $b$ the inclusion intersects the midplane that corresponds to the plane of steady crack propagation. The crack is initiated by the introduction of a pre crack that corresponds to a spatial distribution of the phase field with a thickness of the same order of magnitude as the regularization length. It must be noted that the crack propagates freely over a significant distance before its path is affected by the inclusion. Simulations performed using twice smaller values of $L$, keeping $W = 30$ lead to quantitatively similar results.

To model crack propagation a phase field model is used. It relies on the introduction of an additional variable the phase field $\varphi$ that can be seen as $1-d$ where $d$ is a damage variable. The evolution of the phase field derives from a free energy functional that writes:

$$ \mathcal{F} = \int \left( (hV(\varphi) - \varepsilon_c g(\varphi) + \frac{D}{2} |\nabla\varphi|^2) + g(\varphi) \mathcal{E}_{el} \right) dV \quad (1) $$

$$ \mathcal{E}_{el} = \frac{\lambda(\mathbf{x})}{2}(\mathrm{tr}\epsilon)^2 + \mu(\mathbf{x})\mathrm{tr}(\epsilon^2) \quad (2) $$

where $h = 0.25$, $D = 4$. and $\varepsilon_c = 0.25$ are model parameters and $\mathcal{E}_{el}$ is the elastic energy density for an intact material as a function of the strain tensor $\epsilon$. Its expression involves the space dependant coefficients $\lambda(\mathbf{x})$ and $\mu(\mathbf{x})$ that encode the material heterogeneity as will be described below. $V(\varphi) - \varepsilon_c g(\varphi)$ is a tilted double well potential with $V = \varphi^2(1 - \varphi)^2$ and $g(\varphi) = 4\varphi^3 - 3\varphi^4$ is a coupling function chosen so that when $\varphi = 1$ the contribution of the elastic strain corresponds to the intact material while when $\varphi = 0$, $g = 0$ and the contribution of strain to the energy of the system is 0. It has been shown in [73] that with this model the equilibrium solutions in 1D (local minimas of $\mathcal{F}$) can be either $\varphi = 1$ and uniform strain or localized solutions where $\varphi$ is equal to 1 except in a small region of thickness $w_\varphi$, the so called regularization length, where it goes toward small values ($\ll 1$). If, as it was the case here, the shape of the potential (i.e. the ratio $h/\varepsilon_c$) is kept constant the regularization length scales like $D/h$. With the values considered here, $w_\varphi$ is of the order of 6. SU if $w_\varphi$ is the thickness of the region where $\varphi < 0.5$. In this case the strain is concentrated in this region and the stress in the material is vanishing[73]. The corresponding fracture energy is $2\Gamma = 1.046$ energy units and will be called in the following $G_c$. This model does not present any strain softening[74]

and has the peculiarity that contrarily to what is observed in other approaches such as the Francfort-Marigo variational approach[75-77], there is no crack nucleation. Hence this model is limited to the modelisation of crack propagation. This limitation comes with the advantage that no unphysical crack nucleation can occur when running macroscopic simulations and using a model interface larger than the microscopic process zone. Indeed it has been shown that if the regularisation length used in the model is chosen larger than the actual size of the process zone, then the material stress at failure is reduced proportionally, leading to the possibility of unphysical nucleation events at low stress values. [78, 79].

In order to describe the evolution of a crack the relaxation evolution equation<fnref n="1" /> that derive from this free energy is:

$$ \partial_t \varphi = -\frac{1}{\beta} \frac{\delta \mathcal{F}}{\delta \varphi} \quad (3) $$

<footnote n="1">¹ More complex kinetic equations[80] may be used, especially for fast crack propagation, however in the context of crack propagation the complexity of the events taking place in the process zone would not allow to estimate them.</footnote>
where $1/\beta$ is kinetic parameter that must be positive or 0. It has for value $\beta_0$ here except in two cases:

* if the r.h.s of eq. 3 is positive then $1/\beta$ is 0. This ensures irreversibility.

* $1/\beta$ is $\max(0, (1/\beta_0)(A + g'(\phi)K_{Lamé}(tr\epsilon)^2)/A)$ with $A = -\delta\mathcal{F}/\delta\phi$ if $tr\epsilon < 0$ so that the compression energy does not contribute to crack growth. However this does not prevent crack faces interpenetration.

These two cases do not occur during normal crack propagation in the setup described here (mode I loading). Then, as has been shown in [70], the energy dissipation at the crack tip is proportional to $\beta_0$, taken equal to 4 here, and the total energy of the system is guaranteed to decrease due to fracture propagation:

$$ \partial_t \mathcal{F} = \int \frac{\delta \mathcal{F}}{\delta \varphi} \partial_t \varphi = -\frac{1}{\beta} \int \left( \frac{\delta \mathcal{F}}{\delta \varphi} \right)^2 = -\beta \int (\partial_t \varphi)^2 \quad (4) $$

In addition, due to the double well potential with minima in $\varphi = 0$ and 1 and the zero slope of $g(\varphi)$ at $\varphi = 0$ and 1, the phase field cannot take values out of the 0-1 range. This model has been proven to reproduce well the Griffith criterion [70, 73].

The evolution equation for the elastic field is the wave equation that derives from the Newton’s law with an additional damping term:

$$ \partial_{tt} u_i = - \frac{\delta F}{\delta u_i} - \eta \partial_t u_i \quad (5) $$

where $u_i$ is the displacement field in the $i$ direction with respect to a reference configuration. The damping term $\eta = 0.02$ is chosen large enough to prevent wave reflexions at the boundaries of the sample. It is also chosen small enough to ensure that the wave equation is not overdamped. Here, the chosen value of $\eta$ corresponds to a damping length of the order of 100.SU for an elastic wave, to be compared with the size of the process zone of 3.SU and typical obstacle sizes of 20SU. This term must be seen as a numerical trick to keep crack propagation slow enough.

In this work the effect of spatial heterogeneity is limited to elastic properties of the material. To this purpose space dependant Lamé coefficients $\lambda(\mathbf{x})$ and $\mu(\mathbf{x})$ are considered. They are considered to obey the following law when considering an inclusion of radius $R$:

$$ \lambda(\mathbf{x}) = \lambda_1(h(r)) + \lambda_0(1 - h(r)) \quad (6) $$

$$ \mu(\mathbf{x}) = \mu_1(h(r)) + \mu_0(1 - h(r)) \quad (7) $$

$$ h(r) = \left( \tanh\left(\frac{r-R}{w}\right) + 1 \right) / 2 \quad (8) $$

where $r$ is the distance from the center of the inclusion and $\mu_{0,1}$, $\lambda_{0,1}$ are the Lamé coefficients inside and outside the inclusion. Here for the sake of simplicity, $\lambda_1 = \mu_1 = 1$ and only two cases are presented: $\lambda_0 = \mu_0 = 5$ for a hard inclusion and $\lambda_0 = \mu_0 = 0.1$ for a soft inclusion. This

corresponds to variation of the Young’s modulus, keeping the poisson ratio constant. Obviously other kind of heterogeneities can be considered [68] and may lead to different behaviour. Here, the choice has been to focus more on geometric aspects than on the nature of inclusions. The parameter $w$ sets the abruptness of the elastic moduli change and the smooth variation of elastic moduli used here allows to have inclusions defined implicitly on a regular cartesian grid. It is chosen to be 0.3 SU. Simulations with $w$ ranging between 0.3 and 1.2 SU do not lead to noticeable differences, indicating that as long as $w$ is much smaller than the characteristic size of the inclusion, no significant effects of its value should be observed.

The model equations are solved using a finite difference scheme for spatial derivatives that ensures that the discretized equations derive from a discretized free energy. An explicit Verlet [81] scheme is used for time stepping using a discretized elastic energy in order to ensure mechanical energy conservation (without the small dissipation). The grid is following the crack tip so that its center is kept close to the center of the grid along the x axis (a tread-mill geometry)$^2$. The use of a simple numerical scheme and geometry allows to use efficiently massively parallel architecture such as GPUs. The in house code (using openACC) is used on a single GPU architecture using either an NVIDIA Tesla P100 with double precision) or a NVIDIA RTX3080 or A40 with single precision. In all instances the simulations lead to the same results indicating that the use of single precision is not detrimental here. Most of the simulations presented here have been performed using single precision.

### III. RESULTS

Here the effects of both hard and soft inclusions are described. Since they differ dramatically both cases are discussed in two separate sections. First when the inclusion is soft, second when the inclusion is hard.

#### A. Soft inclusions

Before discussing the three dimensional case, it is useful to give a comprehensive view of the effects of a two dimensional soft inclusion on crack propagation. As is expected from the behavior of a crack in the presence of a void, the soft inclusion attracts the crack since its softness favors the built up of the singularity at the crack tip. This was for instance exemplified in [65](p 1674-1676) where it is shown that the apparent SIFs are higher when a softer region lies ahead of the crack. As a results

<footnote>$^2$ periodically the crack front (for each $z$ the point with the highest $x$ for which $\phi = 0.5$) is computed and the grid is shifted by an integer number of grid points so that the whole front is behind the middle of the grid along the x axis</footnote>
where $1/\beta$ is kinetic parameter that must be positive or 0. It has for value $\beta_0$ here except in two cases:

* if the r.h.s of eq. 3 is positive then $1/\beta$ is 0. This ensures irreversibility.

* $1/\beta$ is $\max(0, (1/\beta_0)(A + g'(\phi)K_{Lamé}(tr\epsilon)^2)/A)$ with $A = -\delta\mathcal{F}/\delta\phi$ if $tr\epsilon < 0$ so that the compression energy does not contribute to crack growth. However this does not prevent crack faces interpenetration.

These two cases do not occur during normal crack propagation in the setup described here (mode I loading). Then, as has been shown in [70], the energy dissipation at the crack tip is proportional to $\beta_0$, taken equal to 4 here, and the total energy of the system is guaranteed to decrease due to fracture propagation:

$$ \partial_t \mathcal{F} = \int \frac{\delta \mathcal{F}}{\delta \varphi} \partial_t \varphi = -\frac{1}{\beta} \int \left( \frac{\delta \mathcal{F}}{\delta \varphi} \right)^2 = -\beta \int (\partial_t \varphi)^2 \quad (4) $$

In addition, due to the double well potential with minima in $\varphi = 0$ and 1 and the zero slope of $g(\varphi)$ at $\varphi = 0$ and 1, the phase field cannot take values out of the 0-1 range. This model has been proven to reproduce well the Griffith criterion [70, 73].

The evolution equation for the elastic field is the wave equation that derives from the Newton’s law with an additional damping term:

$$ \partial_{tt} u_i = - \frac{\delta F}{\delta u_i} - \eta \partial_t u_i \quad (5) $$

where $u_i$ is the displacement field in the $i$ direction with respect to a reference configuration. The damping term $\eta = 0.02$ is chosen large enough to prevent wave reflexions at the boundaries of the sample. It is also chosen small enough to ensure that the wave equation is not overdamped. Here, the chosen value of $\eta$ corresponds to a damping length of the order of 100.SU for an elastic wave, to be compared with the size of the process zone of 3.SU and typical obstacle sizes of 20SU. This term must be seen as a numerical trick to keep crack propagation slow enough.

In this work the effect of spatial heterogeneity is limited to elastic properties of the material. To this purpose space dependant Lamé coefficients $\lambda(\mathbf{x})$ and $\mu(\mathbf{x})$ are considered. They are considered to obey the following law when considering an inclusion of radius $R$:

$$ \lambda(\mathbf{x}) = \lambda_1(h(r)) + \lambda_0(1 - h(r)) \quad (6) $$

$$ \mu(\mathbf{x}) = \mu_1(h(r)) + \mu_0(1 - h(r)) \quad (7) $$

$$ h(r) = \left( \tanh\left(\frac{r-R}{w}\right) + 1 \right) / 2 \quad (8) $$

where $r$ is the distance from the center of the inclusion and $\mu_{0,1}$, $\lambda_{0,1}$ are the Lamé coefficients inside and outside the inclusion. Here for the sake of simplicity, $\lambda_1 = \mu_1 = 1$ and only two cases are presented: $\lambda_0 = \mu_0 = 5$ for a hard inclusion and $\lambda_0 = \mu_0 = 0.1$ for a soft inclusion. This

corresponds to variation of the Young’s modulus, keeping the poisson ratio constant. Obviously other kind of heterogeneities can be considered [68] and may lead to different behaviour. Here, the choice has been to focus more on geometric aspects than on the nature of inclusions. The parameter $w$ sets the abruptness of the elastic moduli change and the smooth variation of elastic moduli used here allows to have inclusions defined implicitly on a regular cartesian grid. It is chosen to be 0.3 SU. Simulations with $w$ ranging between 0.3 and 1.2 SU do not lead to noticeable differences, indicating that as long as $w$ is much smaller than the characteristic size of the inclusion, no significant effects of its value should be observed.

The model equations are solved using a finite difference scheme for spatial derivatives that ensures that the discretized equations derive from a discretized free energy. An explicit Verlet [81] scheme is used for time stepping using a discretized elastic energy in order to ensure mechanical energy conservation (without the small dissipation). The grid is following the crack tip so that its center is kept close to the center of the grid along the x axis (a tread-mill geometry)$^2$. The use of a simple numerical scheme and geometry allows to use efficiently massively parallel architecture such as GPUs. The in house code (using openACC) is used on a single GPU architecture using either an NVIDIA Tesla P100 with double precision) or a NVIDIA RTX3080 or A40 with single precision. In all instances the simulations lead to the same results indicating that the use of single precision is not detrimental here. Most of the simulations presented here have been performed using single precision.

### III. RESULTS

Here the effects of both hard and soft inclusions are described. Since they differ dramatically both cases are discussed in two separate sections. First when the inclusion is soft, second when the inclusion is hard.

#### A. Soft inclusions

Before discussing the three dimensional case, it is useful to give a comprehensive view of the effects of a two dimensional soft inclusion on crack propagation. As is expected from the behavior of a crack in the presence of a void, the soft inclusion attracts the crack since its softness favors the built up of the singularity at the crack tip. This was for instance exemplified in [65](p 1674-1676) where it is shown that the apparent SIFs are higher when a softer region lies ahead of the crack. As a results

<footnote>$^2$ periodically the crack front (for each $z$ the point with the highest $x$ for which $\phi = 0.5$) is computed and the grid is shifted by an integer number of grid points so that the whole front is behind the middle of the grid along the x axis</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 2. Plots of the crack tip trajectories together with the inclusion boundary for 4 different values of the loading in 2D. In both cases the radius of the obstacle is 13.5, and the value of $Y_{off}$ is 19.5 in (a) and 22.5 in b. The crack is propagating from left to right along the x axis. In the inset, the full simulation domain is represented along the y axis and the rightmost half along the x axis. For $Y_{off} = 19$, further increasing the load does not lead the crack to resume propagation.</figcaption></figure>

the crack deflects toward the inclusion. Depending on the load and on the geometry of the system, different behaviors can be observed.

First, if the crack path is sufficiently far from the inclusion, then two cases arise: if the load is sufficient, then the crack will be deflected toward the inclusion (its tip goes away from the middle of the sample) and slow down. Thereafter, once it has passed the vicinity of the inclusion, it retrieve normal propagation and its tip position relaxes toward the center of the sample. If the load is not sufficient, the crack will stop at some point on the path. It should be noted that in this case the crack path does not intersect the inclusion, is mostly determined by the geometry and is not changed significantly by the load. This behavior are illustrated in fig.2 for $G = 1.5G_c$, $G = 1.75G_c$ and $G = 2,0G - c$ where the position of the crack tip are plotted at regular time intervals together with the inclusion. For such values of $G$, the crack velocity ranges from 0.13 to 0.36 $c_s$ where $c_s$ is the shear wave speed (the Raileigh wave speed is $\approx 0.91c_s$ here)$^3$. Hence

<figure><img src="image_5.png" /><figcaption>FIG. 3. (a1):perspective view of the crack front and inclusion at different times during the pinning of the crack. The crack is propagating from left to right and the color of the line is changed according to the time.(b1) top view of the crack front. Only half of the front is shown. The side of the sample is at the bottom of the graph. (c1) side view of the position of the crack tip in different xy planes (center of the sample, side of the sample and an intermediate plane). Graphs use the actual aspect ratio. a2, b2, c2 : same views in a case where the crack front is pinned but does not intersect the inclusion as it is evidenced in c2</figcaption></figure>

here, the effect of changing $G$ is to change the way the crack goes along a *prescribed path*: for low values of $G$ it slows down and eventually stops and for higher values it slows downs and eventually retrieves normal propagation. One should note that once the crack has stopped, increasing $G$ will make advance and eventually retrieve normal propagation. .

Second, when the inclusion is closer to the crack path, independently of the load, the behavior is dramatically different. The crack deflects toward the inclusion and propagates until it crosses its boundary. There it stops as can be expected from [65] where it is shown that the stress intensity factor at a crack tip in a soft inclusion is much smaller than the stress intensity factor in the bulk. In this case, when the load is increased, once the crack tip has entered the inclusion, it will not resume propagation unless a crack nucleation occurs, for instance at the boundary of the inclusion.

The boundary between the different regimes, for a given value of $G = 1.5G_c$ are represented on fig. 5. It should be noted that in this phase diagram, there is no distinction between the two kinds of crack arrest. This choice is due to the fact that this phase diagram serves the purpose of comparing 2D and 3D results. For the later, as will be seen, such a distinction is not relevant.

<footnote>$^3$ Using eq. 5.4.9 of [69], one can estimate the dissipation at the crack tip induced by the $\beta$ and additional dissipation term. For</footnote>

<footnote>$G = 1.25G_c$ it is $0.238G_c$, for $G = 1.5G_c$ it is $0.446G_c$ and for $G = 1.75G_c$ it is $0.669G_c$.</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 2. Plots of the crack tip trajectories together with the inclusion boundary for 4 different values of the loading in 2D. In both cases the radius of the obstacle is 13.5, and the value of $Y_{off}$ is 19.5 in (a) and 22.5 in b. The crack is propagating from left to right along the x axis. In the inset, the full simulation domain is represented along the y axis and the rightmost half along the x axis. For $Y_{off} = 19$, further increasing the load does not lead the crack to resume propagation.</figcaption></figure>

the crack deflects toward the inclusion. Depending on the load and on the geometry of the system, different behaviors can be observed.

First, if the crack path is sufficiently far from the inclusion, then two cases arise: if the load is sufficient, then the crack will be deflected toward the inclusion (its tip goes away from the middle of the sample) and slow down. Thereafter, once it has passed the vicinity of the inclusion, it retrieve normal propagation and its tip position relaxes toward the center of the sample. If the load is not sufficient, the crack will stop at some point on the path. It should be noted that in this case the crack path does not intersect the inclusion, is mostly determined by the geometry and is not changed significantly by the load. This behavior are illustrated in fig.2 for $G = 1.5G_c$, $G = 1.75G_c$ and $G = 2,0G - c$ where the position of the crack tip are plotted at regular time intervals together with the inclusion. For such values of $G$, the crack velocity ranges from 0.13 to 0.36 $c_s$ where $c_s$ is the shear wave speed (the Raileigh wave speed is $\approx 0.91c_s$ here)$^3$. Hence

<figure><img src="image_5.png" /><figcaption>FIG. 3. (a1):perspective view of the crack front and inclusion at different times during the pinning of the crack. The crack is propagating from left to right and the color of the line is changed according to the time.(b1) top view of the crack front. Only half of the front is shown. The side of the sample is at the bottom of the graph. (c1) side view of the position of the crack tip in different xy planes (center of the sample, side of the sample and an intermediate plane). Graphs use the actual aspect ratio. a2, b2, c2 : same views in a case where the crack front is pinned but does not intersect the inclusion as it is evidenced in c2</figcaption></figure>

here, the effect of changing $G$ is to change the way the crack goes along a *prescribed path*: for low values of $G$ it slows down and eventually stops and for higher values it slows downs and eventually retrieves normal propagation. One should note that once the crack has stopped, increasing $G$ will make advance and eventually retrieve normal propagation. .

Second, when the inclusion is closer to the crack path, independently of the load, the behavior is dramatically different. The crack deflects toward the inclusion and propagates until it crosses its boundary. There it stops as can be expected from [65] where it is shown that the stress intensity factor at a crack tip in a soft inclusion is much smaller than the stress intensity factor in the bulk. In this case, when the load is increased, once the crack tip has entered the inclusion, it will not resume propagation unless a crack nucleation occurs, for instance at the boundary of the inclusion.

The boundary between the different regimes, for a given value of $G = 1.5G_c$ are represented on fig. 5. It should be noted that in this phase diagram, there is no distinction between the two kinds of crack arrest. This choice is due to the fact that this phase diagram serves the purpose of comparing 2D and 3D results. For the later, as will be seen, such a distinction is not relevant.

<footnote n="3">³ Using eq. 5.4.9 of [69], one can estimate the dissipation at the crack tip induced by the $\beta$ and additional dissipation term. For</footnote>

<footnote>$G = 1.25G_c$ it is $0.238G_c$, for $G = 1.5G_c$ it is $0.446G_c$ and for $G = 1.75G_c$ it is $0.669G_c$.</footnote>
In the 3D case, the fact that the inclusion does not
occupy the full width of the sample induces different be-
havior. First, as in 2D it can induce the arrest of the
crack front. However, while in 2D, the crack is simply
stopped, here, the crack front is pinned on the inclusion
and outside of the soft region the crack front is deformed
as has been observed experimentally when the crack front
encounters a tougher region [82]. This is exemplified in
fig. 3 (b1) where one can see a time series of crack front
seen from above, a perspective view(a1) and some cuts
(xy planes) of the crack path(c1). One can see that the
crack front is curved and that the crack surface is bent
toward the inclusion. In this situation the crack front
can be divided into three parts: a central part where the
crack front has reached the inclusion and the two side
parts where the crack front is in the homogeneous do-
main and cannot advance due to the pinned part of the
crack front. This behaviour is similar to the pinning
mechanism described in [83] and was also discussed in
[61] where it was was proven to contribute to toughening
for cracks confined in a plane. As in 2D there exists a
situation where the crack front is stopped while it is not
intersecting the inclusion at all as it is illustrated in fig.
3 (a2, b2 and c2). In these cases, the front cannot be
described as three parts.

In both these pinning situations (3 (c1) and (c2)), the
normal propagation of the crack can be retrieved by sim-
ply increasing the load (once the crack has stopped) with-
out crack nucleation. While this is in contrast with the
2D case where, if the crack front intersect the boundary
of the inclusion, a nucleation event is needed to retrieve
normal crack propagation, it can be expected. Indeed,
since the inclusion does not span the whole sample, there
are always parts of the crack front that do not intersect
the inclusion and for which the *SIF reduction* induced
by the inclusion is absent. For such parts they will ob-
viously advance as the load is increased and the front
motion observed in partial pinning situations (**that is de-
scribed below and illustrated in fig. 4**) is expected.

When considering either higher loads, smaller inclu-
sions or inclusions further apart from the midplane, an
almost unperturbed crack propagation is observed. This
regime is not shown, but at intermediate values of load
there exists a third regime,that is called partial pinning
here, that is worth mentioning and describing. It con-
sists of a pinning of the crack at the inclusion followed by
a slow propagation around it and finally a normal crack
propagation. This is illustrated in figure4(a) and (b)
where the position of the crack front at different equally
spaced times is plotted in perspective and from above
(only half the simulation is shown there). One can see
that the spacing between the front line is dramatically
reduced indicating a slow down of the crack front. It is
then divided into 3 regions. One in the center is pinned
at the inclusion. The other two at the boundary of the
samples are stopped by curvature and by the fact that
since the inclusion is not broken the singularity cannot
fully develop. However they advance slowly around the

<figure><img src="image_6.png" /><figcaption>FIG. 4. **Top** (a):perspective view of the crack front and inclusion at different times during the slow down and unpinning of the crack. The crack is propagating from left to right and the color of the line is changed according to the time.(b) top view of the crack front. Only half of the front is shown. The side of the sample is at the bottom of the graph. (c) side view of the position of the crack tip in different xy planes (center of the sample, side of the sample and an intermediate plane).**Bottom** Sequence of images of the crack (iso surfaceϕ = 0.5 of the phase field) when it goes around the inclusion. The images are labeled by a dimensionless time variable. One can see that between the images 2 and 8, the crack moves very little. This contrasts with the advance of the crack between images 1 and 2 or 8 and 9. It is worth mentioning that the crack isosurface may intersect the sphere while the crack tip is less than a regularization length away from the obstacle</figcaption></figure>

inclusion and once their *free ends* meet the crack resumes
normal fast propagation. Thereafter, the inclusion re-
mains intact **and forms a ligament that lies behind the**
**crack path** for a small time and eventually breaks when
the crack front has advanced enough. This behavior is
also illustrated in figure 4 where the crack surface (the
isosurface *φ* = 0.5) is plotted at different times together
with the spherical inclusion. Once again, it is interesting
to note that this behaviour is similar to the one discussed
in [61] since it implies that for at least some time, intact
inclusions will remain behind the crack acting as *bridges*
and participating potentially to the toughening of the
material. **It must be noted that such ligaments are remi-
niscent of patterns that can be seen during the rupture
of some elastomers [84] at small scales.**

Before turning to the results in the case of an hard
inclusion we represent in figure 5 the different observed
behavior for a given load and different sizes of the inclu-
In the 3D case, the fact that the inclusion does not
occupy the full width of the sample induces different be-
havior. First, as in 2D it can induce the arrest of the
crack front. However, while in 2D, the crack is simply
stopped, here, the crack front is pinned on the inclusion
and outside of the soft region the crack front is deformed
as has been observed experimentally when the crack front
encounters a tougher region [82]. This is exemplified in
fig. 3 (b1) where one can see a time series of crack front
seen from above, a perspective view(a1) and some cuts
(xy planes) of the crack path(c1). One can see that the
crack front is curved and that the crack surface is bent
toward the inclusion. In this situation the crack front
can be divided into three parts: a central part where the
crack front has reached the inclusion and the two side
parts where the crack front is in the homogeneous do-
main and cannot advance due to the pinned part of the
crack front. This behaviour is similar to the pinning
mechanism described in [83] and was also discussed in
[61] where it was was proven to contribute to toughening
for cracks confined in a plane. As in 2D there exists a
situation where the crack front is stopped while it is not
intersecting the inclusion at all as it is illustrated in fig.
3 (a2, b2 and c2). In these cases, the front cannot be
described as three parts.

In both these pinning situations (3 (c1) and (c2)), the
normal propagation of the crack can be retrieved by sim-
ply increasing the load (once the crack has stopped) with-
out crack nucleation. While this is in contrast with the
2D case where, if the crack front intersect the boundary
of the inclusion, a nucleation event is needed to retrieve
normal crack propagation, it can be expected. Indeed,
since the inclusion does not span the whole sample, there
are always parts of the crack front that do not intersect
the inclusion and for which the *SIF reduction* induced
by the inclusion is absent. For such parts they will ob-
viously advance as the load is increased and the front
motion observed in partial pinning situations (**that is de-
scribed below and illustrated in fig. 4**) is expected.

When considering either higher loads, smaller inclu-
sions or inclusions further apart from the midplane, an
almost *unperturbed* crack propagation is observed. This
regime is not shown, but at intermediate values of load
there exists a third regime, *that is called partial pinning
here*, that is worth mentioning and describing. It con-
sists of a pinning of the crack at the inclusion followed by
a slow propagation around it and finally a normal crack
propagation. This is illustrated in figure4(a) and (b)
where the position of the crack front at different equally
spaced times is plotted in perspective and from above
(only half the simulation is shown there). One can see
that the spacing between the front line is dramatically
reduced indicating a slow down of the crack front. It is
then divided into 3 regions. One in the center is pinned
at the inclusion. The other two at the boundary of the
samples are stopped by curvature and by the fact that
since the inclusion is not broken the singularity cannot
fully develop. However they advance slowly around the

<figure><img src="image_6.png" /><figcaption>FIG. 4. **Top** (a):perspective view of the crack front and inclusion at different times during the slow down and unpinning of the crack. The crack is propagating from left to right and the color of the line is changed according to the time.(b) top view of the crack front. Only half of the front is shown. The side of the sample is at the bottom of the graph. (c) side view of the position of the crack tip in different xy planes (center of the sample, side of the sample and an intermediate plane). **Bottom** Sequence of images of the crack (iso surfaceφ = 0.5 of the phase field) when it goes around the inclusion. The images are labeled by a dimensionless time variable. One can see that between the images 2 and 8, the crack moves very little. This contrasts with the advance of the crack between images 1 and 2 or 8 and 9. It is worth mentioning that the crack isosurface may intersect the sphere while the crack tip is less than a regularization length away from the obstacle</figcaption></figure>

inclusion and once their *free ends* meet the crack resumes
normal fast propagation. Thereafter, the inclusion re-
mains intact and forms a ligament that lies behind the
crack path for a small time and eventually breaks when
the crack front has advanced enough. This behavior is
also illustrated in figure 4 where the crack surface (the
isosurface $\phi = 0.5$) is plotted at different times together
with the spherical inclusion. Once again, it is interesting
to note that this behaviour is similar to the one discussed
in [61] since it implies that for at least some time, intact
inclusions will remain behind the crack acting as *bridges*
and participating potentially to the toughening of the
material. *It must be noted that such ligaments are remi-
niscent of patterns that can be seen during the rupture
of some elastomers* [84] *at small scales.*

Before turning to the results in the case of an hard
inclusion we represent in figure 5 the different observed
behavior for a given load and different sizes of the inclu-
<figure><img src="image_7.png" /><figcaption>FIG. 5. Phase diagram in the $R, Y$ plane for simulations for soft inclusions (2D in a and 3D in b) and hard inclusions (2D in c and 3D in d). G is set to $1.5G_c$ for soft inclusions and to $1.25G_c$ for hard inclusions. The criteria chosen to differentiate between freely moving is arbitrary. Here it is chosen to consider that the later was observed when the minimal crack front velocity is smaller than 1/2 its free propagation speed. Obviously, changing the criteria would lead to different shape of the partial pinning region.</figcaption></figure>

sion and offset of its center from the plane of propagation of the crack. One can see that there is a minimal radius below which no pinning of the crack can be observed and solely a slow down (that can be important) is possible. This differs significantly from the 2D case where the transition line between pinning and propagation can be extrapolated down to a 0 radius of the inclusion. This, together with the results described above illustrates clearly the difference between a 2D and a 3D inclusion.

## B. Hard inclusions

We now turn to the case of hard inclusions and proceed similarly. First in 2D systems, the inclusion tends to repel the crack. Indeed, the presence of an hard inclusion tends to prevent large strains in the surrounding material. As a result crack propagation away from the inclusion, in region with nominal elastic moduli, is favored. As a result the crack is repelled. This can, again, be understood using [65] where it was shown that when the material ahead of the crack tip is tougher, the SIF is smaller.

However, the hard inclusion still can pin the crack. Indeed, when the crack tip gets too close to the inclusion, the higher elastic moduli prevents the singularity at the

<figure><img src="image_8.png" /><figcaption>FIG. 6. Snapshots of the crack going around an hard inclusion at different time of propagation The crack surface (isosurface $\phi = 0.5$) is represented in blue while the spherical inclusion is in gray. At $t=4$ and $t=5$ one can see that the crack breaks partly the inclusion. The lower rightmost image corresponds to a crack stopped by an hard inclusion. In all cases, only a portion of the simulation domain is represented.</figcaption></figure>

crack tip to be strong it may lead to the arrest of crack propagation. In this situation, unlike with soft inclusions, the crack never crosses the boundary of the inclusion. As a result, simply increasing the load will lead to crack propagation without the nucleation of a new crack.

In 3D, a spherical hard inclusion can still lead to the arrest of a crack, in a similar manner: most of the crack front is stopped at the inclusion as is shown fig.6 on the bottom right snapshot ($R=12$). In cases where the inclusion is smaller or more off the crack propagation plane, the crack front will go around the inclusion with a significant slow down. This is exemplified in fig.6. Once the crack front is close to the inclusion, it is stopped ($T=2$) and starts to propagate slowly around it ($t=3, t=4$). Once the inclusion is no longer ahead of the front, it retrieves normal propagation ($t=5$). Similar behavior can be seen with smaller inclusions that do not span the whole system. A phase diagram for a fixed value of $G = 1.25G_c$ is presented in fig. 5. It must be noted that the inclusion will affect significantly the crack propagation only if it lies in its path. Moreover, if the center of the spherical inclusion lies in the close vicinity of the crack path it is dramatically more efficient at pinning the crack. The contrast between the *local* effect of hard inclusion and the *non local* effect of soft inclusions lies in the fact that the hard (resp. soft) inclusion tends to repel (resp. attract) the crack, as a result the crack is further away from (resp. closer to ) the inclusion and is less (resp. more) affected by the inclusion.

This shows that the difference in elastic moduli can be sufficient to prevent crack propagation through an hard inclusion and that it can lead to the arrest of the crack. However, it must be noted that it appears to be less efficient. This may be related to the fact that in the case of hard inclusions, they deflect cracks. Indeed it has been shown that crack deflection[57] is less efficient that pinning or bridging[61, 83]. It must also be noted that when there is no crack arrest, the presence of the inclusion can still affect significantly crack propagation. Indeed, the crack front propagates around it and its propagation ve-
<figure><img src="image_7.png" /><figcaption>FIG. 5. Phase diagram in the $R, Y$ plane for simulations for soft inclusions (2D in a and 3D in b) and hard inclusions (2D in c and 3D in d). G is set to $1.5G_c$ for soft inclusions and to $1.25G_c$ for hard inclusions. The criteria chosen to differentiate between freely moving is arbitrary. Here it is chosen to consider that the later was observed when the minimal crack front velocity is smaller than 1/2 its free propagation speed. Obviously, changing the criteria would lead to different shape of the partial pinning region.</figcaption></figure>

sion and offset of its center from the plane of propagation of the crack. One can see that there is a minimal radius below which no pinning of the crack can be observed and solely a slow down (that can be important) is possible. This differs significantly from the 2D case where the transition line between pinning and propagation can be extrapolated down to a 0 radius of the inclusion. This, together with the results described above illustrates clearly the difference between a 2D and a 3D inclusion.

## B. Hard inclusions

We now turn to the case of hard inclusions and proceed similarly. First in 2D systems, the inclusion tends to repel the crack. Indeed, the presence of an hard inclusion tends to prevent large strains in the surrounding material. As a result crack propagation away from the inclusion, in region with nominal elastic moduli, is favored. As a result the crack is repelled. This can, again, be understood using [65] where it was shown that when the material ahead of the crack tip is tougher, the SIF is smaller.

However, the hard inclusion still can pin the crack. Indeed, when the crack tip gets too close to the inclusion, the higher elastic moduli prevents the singularity at the

<figure><img src="image_8.png" /><figcaption>FIG. 6. Snapshots of the crack going around an hard inclusion at different time of propagation The crack surface (isosurface $\phi = 0.5$) is represented in blue while the spherical inclusion is in gray. At $t=4$ and $t=5$ one can see that the crack breaks partly the inclusion. The lower rightmost image corresponds to a crack stopped by an hard inclusion. In all cases, only a portion of the simulation domain is represented.</figcaption></figure>

crack tip to be strong it may lead to the arrest of crack propagation. In this situation, unlike with soft inclusions, the crack never crosses the boundary of the inclusion. As a result, simply increasing the load will lead to crack propagation without the nucleation of a new crack.

In 3D, a spherical hard inclusion can still lead to the arrest of a crack, in a similar manner: most of the crack front is stopped at the inclusion as is shown fig.6 on the bottom right snapshot ($R=12$). In cases where the inclusion is smaller or more off the crack propagation plane, the crack front will go around the inclusion with a significant slow down. This is exemplified in fig.6. Once the crack front is close to the inclusion, it is stopped ($T=2$) and starts to propagate slowly around it ($t=3, t=4$). Once the inclusion is no longer ahead of the front, it retrieves normal propagation ($t=5$). Similar behavior can be seen with smaller inclusions that do not span the whole system. A phase diagram for a fixed value of $G = 1.25G_c$ is presented in fig. 5. It must be noted that the inclusion will affect significantly the crack propagation only if it lies in its path. Moreover, if the center of the spherical inclusion lies in the close vicinity of the crack path it is dramatically more efficient at pinning the crack. The contrast between the *local* effect of hard inclusion and the *non local* effect of soft inclusions lies in the fact that the hard (resp. soft) inclusion tends to repel (resp. attract) the crack, as a result the crack is further away from (resp. closer to ) the inclusion and is less (resp. more) affected by the inclusion.

This shows that the difference in elastic moduli can be sufficient to prevent crack propagation through an hard inclusion and that it can lead to the arrest of the crack. However, it must be noted that it appears to be less efficient. This may be related to the fact that in the case of hard inclusions, they deflect cracks. Indeed it has been shown that crack deflection[57] is less efficient that pinning or bridging[61, 83]. It must also be noted that when there is no crack arrest, the presence of the inclusion can still affect significantly crack propagation. Indeed, the crack front propagates around it and its propagation ve-
locity can dramatically decrease.

IV. CONCLUSION AND PERSPECTIVE

Here we have shown that, when considering crack prop-
agation in heterogeneous materials, elastic effects alone
can have a dramatic effect on crack propagation even
when at a macroscopic scale the heterogeneity is barely
noticeable when considering the outer stress field. The
results of the numerical simulations show that in both
the cases of a soft and of an hard inclusion, it is possible
to pin the crack at the inclusion. However the way the
crack is affected differ qualitatively in two aspects. First
in the case of a soft inclusion, the effects can be signifi-
cant when the inclusion is outside of the crack path while
in the case of hard inclusions, it is not the case. Second,
the unpinning mechanism are dramatically different. In
the case of the soft inclusion, the crack front advances,
moves around it in two parts that eventually reconnect
resume normal propagation and finally once the crack
front is sufficiently far away from the inclusion it even-
tually breaks, as it has already been discussed in [61] in

the case of toughness heterogeneities for a confined crack
front. In the case of the hard inclusion, the crack fronts
keeps it connectivity and the crack surface goes around
the inclusion. In both cases, it must be noted that the
unpinning mechanism does not involve the nucleation of
a new crack contrarily to what is expected in 2D systems.

This works gives some insight on the effects of het-
erogeneities on crack propagation. It paves the way to
further work that would aim at considering more complex
systems such as the effects of either disordered or crys-
talline arrangements of spherical inclusions, more com-
plex inclusions (elliptical, non convex) and also of differ-
ent kind of inclusions (with varying poison ratio, eigen-
strain). Such studies should help creating architectured
materials with higher fracture energy[85]. In addition,
the dramatic slow down of a crack that is observed here,
may be of importance for statistical models of crack front
dynamic in a quenched disorder that are used to model
earthquakes[86].

Another extension of this work that is also needed is to
consider the effects of fracture toughness heterogeneities.
However, in this case, a simple description like the one
presented here is unlikely to be successful[87] and more
complex approaches should be considered[88-90].

[1] A. H.-D. Cheng. *Poroelasticity*. Springer, 2016.

[2] L. Dormieux, A. Molinari, and D. Kondo. Micromechanical approach to the behavior of poroelastic materials. *Journal of the Mechanics and Physics of Solids*, 50(10):2203-2231, 2002.

[3] G. Lu, G.Q. (M.) Lu, and Z.M. Xiao. Mechanical properties of porous materials. *Journal of Porous Materials*, 6:359–368, 1999.

[4] N. Abid, M. Mirkhalaf, and F. Barthelat. Discrete-element modeling of nacre-like materials: Effects of random microstructures on strain localization and mechanical performance. *Journal of the Mechanics and Physics of Solids*, 112:385-402, 2018.

[5] A. G. Evans, Z. Suo, R. Z. Wang, I. A. Aksay, M. Y. He, and J. W. Hutchinson. Model for the robust mechanical behavior of nacre. *Journal of Materials Research*, 16(09):2475-2484, 2001.

[6] K. Okumura and P.-G. de Gennes. Why is nacre strong? Elastic theory and fracture mechanics for biocomposites with stratified structures. *The European Physical Journal E*, 4(1):121–127, 2001.

[7] A. G. Evans, B.J. Dalgleish, M. Y. He, and J. W. Hutchinson. On crack path selection an the interface fracture energy in bimaterial systems. *Acta Metallurgica*, 37(12):3249–3254, 1989.

[8] M. Y. He, A. Bartlett, A. G. Evans, and J. W. Hutchinson. kinking of a crack out of an interface: role of in plane stress. *J. Am Ceram Soc.*, 74(4):767-771, 1991.

[9] M. Y. He, A. G. Evans, and J. W. Hutchinson. Crack Deflection at an interface between dissimilar elastic materials: role of residual stresses. *Int. J. Solids Structures*, 31(24):3443-3455, 1994.

[10] M. Y. He, A. G. Evans, and J. W. Hutchinson. Interface

cracking phenomena in constrained metal layers. *Acta Materialia*, 44(7):2963-2971, 996.

[11] M. Z. Hossain, C.-J. Hsueh, B. Bourdin, and K. Bhattacharya. Effective toughness of heterogeneous media. *Journal of the Mechanics and Physics of Solids*, 71:15-32, 2014.

[12] M. Lebihain, J.-B. Leblond, and L. Ponson. Effective toughness of periodic heterogeneous materials: the effect of out-of-plane excursions of cracks. *Journal of the Mechanics and Physics of Solids*, 137:103876, 2020.

[13] D. Schneider, E. Schoof, Y. Huang, M. Selzer, and B. Nestler. Phase-field modeling of crack propagation in multiphase systems. *Computer Methods in Applied Mechanics and Engineering*, 312:186-195, 2016.

[14] K.P. Herrmann and R. Wang. Interaction of an antiplane shear crack with a transforming inclusion. *Engineering Fracture Mechanics*, 50(3):361-368, 1995.

[15] R. R. Bhargava. A misfitting elastic inclusion in an infinite plane containing a crack. *Journal of Elasticity*, 7(2):201–211, April 1977.

[16] A.F. Bower and M. Ortiz. A three-dimensional analysis of crack trapping and bridging by tough particles. *Journal of the Mechanics and Physics of Solids*, 39(6):815-858, 1991.

[17] F. Erdogan, G. D. Gupta, and M. Ratwani. Interaction Between a Circular Inclusion and an Arbitrarily Oriented Crack. *Journal of Applied Mechanics*, 41(4):1007-1013, 12 1974.

[18] D. Brescakovic, M. Kegl, and O. Kolednik. Interaction of crack and hole: effects on crack trajectory, crack driving force and fracture toughness. *International Journal of Fracture*, 236(1):33-57, July 2022.

[19] E. M. Patton and M. H. Santare. The effect of a rigid el-
locity can dramatically decrease.

IV. CONCLUSION AND PERSPECTIVE

Here we have shown that, when considering crack prop-
agation in heterogeneous materials, elastic effects alone
can have a dramatic effect on crack propagation even
when at a macroscopic scale the heterogeneity is barely
noticeable when considering the outer stress field. The
results of the numerical simulations show that in both
the cases of a soft and of an hard inclusion, it is possible
to pin the crack at the inclusion. However the way the
crack is affected differ qualitatively in two aspects. First
in the case of a soft inclusion, the effects can be signifi-
cant when the inclusion is outside of the crack path while
in the case of hard inclusions, it is not the case. Second,
the unpinning mechanism are dramatically different. In
the case of the soft inclusion, the crack front advances,
moves around it in two parts that eventually reconnect
resume normal propagation and finally once the crack
front is sufficiently far away from the inclusion it even-
tually breaks, as it has already been discussed in [61] in

the case of toughness heterogeneities for a confined crack
front. In the case of the hard inclusion, the crack fronts
keeps it connectivity and the crack surface goes around
the inclusion. In both cases, it must be noted that the
unpinning mechanism does not involve the nucleation of
a new crack contrarily to what is expected in 2D systems.

This works gives some insight on the effects of het-
erogeneities on crack propagation. It paves the way to
further work that would aim at considering more complex
systems such as the effects of either disordered or crys-
talline arrangements of spherical inclusions, more com-
plex inclusions (elliptical, non convex) and also of differ-
ent kind of inclusions (with varying poison ratio, eigen-
strain). Such studies should help creating architectured
materials with higher fracture energy[85]. In addition,
the dramatic slow down of a crack that is observed here,
may be of importance for statistical models of crack front
dynamic in a quenched disorder that are used to model
earthquakes[86].

Another extension of this work that is also needed is to
consider the effects of fracture toughness heterogeneities.
However, in this case, a simple description like the one
presented here is unlikely to be successful[87] and more
complex approaches should be considered[88-90].

[1] A. H.-D. Cheng. *Poroelasticity*. Springer, 2016.

[2] L. Dormieux, A. Molinari, and D. Kondo. Micromechanical approach to the behavior of poroelastic materials. *Journal of the Mechanics and Physics of Solids*, 50(10):2203-2231, 2002.

[3] G. Lu, G.Q. (M.) Lu, and Z.M. Xiao. Mechanical properties of porous materials. *Journal of Porous Materials*, 6:359–368, 1999.

[4] N. Abid, M. Mirkhalaf, and F. Barthelat. Discrete-element modeling of nacre-like materials: Effects of random microstructures on strain localization and mechanical performance. *Journal of the Mechanics and Physics of Solids*, 112:385-402, 2018.

[5] A. G. Evans, Z. Suo, R. Z. Wang, I. A. Aksay, M. Y. He, and J. W. Hutchinson. Model for the robust mechanical behavior of nacre. *Journal of Materials Research*, 16(09):2475-2484, 2001.

[6] K. Okumura and P.-G. de Gennes. Why is nacre strong? Elastic theory and fracture mechanics for biocomposites with stratified structures. *The European Physical Journal E*, 4(1):121–127, 2001.

[7] A. G. Evans, B.J. Dalgleish, M. Y. He, and J. W. Hutchinson. On crack path selection an the interface fracture energy in bimaterial systems. *Acta Metallurgica*, 37(12):3249–3254, 1989.

[8] M. Y. He, A. Bartlett, A. G. Evans, and J. W. Hutchinson. kinking of a crack out of an interface: role of in plane stress. *J. Am Ceram Soc.*, 74(4):767–771, 1991.

[9] M. Y. He, A. G. Evans, and J. W. Hutchinson. Crack Deflection at an interface between dissimilar elastic materials: role of residual stresses. *Int. J. Solids Structures*, 31(24):3443–3455, 1994.

[10] M. Y. He, A. G. Evans, and J. W. Hutchinson. Interface

cracking phenomena in constrained metal layers. *Acta Materialia*, 44(7):2963-2971, 996.

[11] M. Z. Hossain, C.-J. Hsueh, B. Bourdin, and K. Bhattacharya. Effective toughness of heterogeneous media. *Journal of the Mechanics and Physics of Solids*, 71:15-32, 2014.

[12] M. Lebihain, J.-B. Leblond, and L. Ponson. Effective toughness of periodic heterogeneous materials: the effect of out-of-plane excursions of cracks. *Journal of the Mechanics and Physics of Solids*, 137:103876, 2020.

[13] D. Schneider, E. Schoof, Y. Huang, M. Selzer, and B. Nestler. Phase-field modeling of crack propagation in multiphase systems. *Computer Methods in Applied Mechanics and Engineering*, 312:186-195, 2016.

[14] K.P. Herrmann and R. Wang. Interaction of an antiplane shear crack with a transforming inclusion. *Engineering Fracture Mechanics*, 50(3):361-368, 1995.

[15] R. R. Bhargava. A misfitting elastic inclusion in an infinite plane containing a crack. *Journal of Elasticity*, 7(2):201–211, April 1977.

[16] A.F. Bower and M. Ortiz. A three-dimensional analysis of
crack trapping and bridging by tough particles. *Journal*
*of the Mechanics and Physics of Solids*, 39(6):815–858,
1991.

[17] F. Erdogan, G. D. Gupta, and M. Ratwani. Interaction Between a Circular Inclusion and an Arbitrarily Oriented Crack. *Journal of Applied Mechanics*, 41(4):1007–1013, 12 1974.

[18] D. Brescakovic, M. Kegl, and O. Kolednik. Interaction of crack and hole: effects on crack trajectory, crack driving force and fracture toughness. *International Journal of Fracture*, 236(1):33–57, July 2022.

[19] E. M. Patton and M. H. Santare. The effect of a rigid el-
liptical inclusion on a straight crack. *International Journal of Fracture*, 46(1):71–79, November 1990.

[20] H. Gao. Fracture analysis of nonhomogeneous materials via a moduli-perturbation approach. *International Journal of Solids and Structures*, 27(13):1663–1682, 1991.

[21] S. Brach, M.Z. Hossain, B. Bourdin, and K. Bhat-tacharya. Anisotropy of the effective toughness of layered media. *Journal of the Mechanics and Physics of Solids*, 131:96–111, 2019.

[22] J.D. Clayton and J. Knap. A geometrically nonlinear phase field theory of brittle fracture. *International Journal of Fracture*, 189(2):139–148, 2014.

[23] T. Lacondemine, C. Roux-Langlois, and T. Rouxel. Role of poisson's ratio mismatch on the crack path in glass matrix particulate composites. *International Journal of Fracture*, 207(1):73–85, September 2017.

[24] L B Freund. *Dynamic Fracture Mechanics*. Cambridge University Press (UK), 1990.

[25] H. Henry. Study of the branching instability using a phase field model of inplane crack propagation. *Euro-physics Letters*, 83:16004, 2008.

[26] Z.P. Bazant and L.F. Estenssoro. Surface singularity and crack propagation. *International journal of solids and structure*, 15:405–426, 1979.

[27] H. Henry. Study of three-dimensional crack fronts under plane stress using a phase field model. *EPL (Europhysics Letters)*, 92(4):46002, nov 2010.

[28] A. Karma, D.A Kessler, and H. Levine. Phase-Field Model of Mode III Dynamic Fracture. *Phys. Rev. Lett.*, 87:45501, 2001.

[29] C. Kuhn, A. Schlüter, and R. Müller. On degradation functions in phase field fracture models. *Computational Materials Science*, 108:374–384, oct 2015.

[30] G.A. Francfort and J.-J. Marigo. Revisiting brittle fracture as an energy minimization problem. *Journal of the Mechanics and Physics of Solids*, 46(8):1319–1342, aug 1998.

[31] B. Bourdin, G. A. Francfort, and J.-J. Marigo. *The variational approach to fracture*. 2008.

[32] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations. *International Journal for Numerical Methods in Engineering*, 83(10):1273–1311, 2010.

[33] C. Maurini, B. Bourdin, G. Gauthier, and V. Lazarus. Crack patterns obtained by unidirectional drying of a colloidal suspension in a capillary tube: Experiments and numerical simulations using a two-dimensional variational approach. *Int. J. Frac.*, 184:75–91, 2013.

[34] E. Tanné, T. Li, B. Bourdin, J. J. Marigo, and C. Maurini. Crack nucleation in variational phase-field models of brittle fracture. *Journal of the Mechanics and Physics of Solids*, 110:80–99, jan 2018.

[35] Peter Galenko and David Jou. Diffuse-interface model for rapid phase transformations in nonequilibrium systems. *Phys. Rev. E*, 71:046125, Apr 2005.

[36] L. Verlet. Computer "experiments" on classical fluids. i. thermodynamical properties of lennard-jones molecules. *Phys. Rev.*, 159(1):98–103, July 1967.

[37] D. Dalmas, E. Barthel, and D. Vandembroucq. Crack front pinning by design in planar heterogeneous interfaces. *Journal of the Mechanics and Physics of Solids*, 57(3):446–457, 2009.

[38] H. Gao and J.R. Rice. A First-Order Perturbation Analysis of Crack Trapping by Arrays of Obstacles. *Journal of Applied Mechanics*, 56(4):828–836, 12 1989.

[39] S.V. Hainsworth. An environmental scanning electron microscopy investigation of fatigue crack initiation and propagation in elastomers. *Polymer Testing*, 26(1):60–70, 2007.

[40] S. Choukir and C.V. Singh. Role of topology in dictating the fracture toughness of mechanical metamaterials. *International Journal of Mechanical Sciences*, 241:107945, 2023.

[41] E. A. Jagla, François P. Landes, and Alberto Rosso. Viscoelastic effects in avalanche dynamics: A key to earthquake statistics. *Phys. Rev. Lett.*, 112:174301, Apr 2014.

[42] H. Henry. Limitations of the modelling of crack propagating through heterogeneous material using a phase field approach. *Theoretical and Applied Fracture Mechanics*, 104:102384, 2019.

[43] T.T. Nguyen, J. Yvonnet, Q.-Z. Zhu, M. Bornert, and C. Chateau. A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography. *Computer Methods in Applied Mechanics and Engineering*, 2015.

[44] T T Nguyen, J Yvonnet, M Bornert, and C Chateau. Initiation and propagation of complex 3D networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microCT experiments and phase field simulations. *Journal of the Mechanics and Physics of Solids*, 95:320–350, 2016.

[45] L. Schöller, D. Schneider, A. Prahs, and B. Nestler. Phase-field modeling of crack propagation based on multi-crack order parameters considering mechanical jump conditions. *PAMM*, 22(1):e202200039, 2023.

[46] A. H.-D. Cheng. *Poroelasticity*. Springer, 2016.

[47] L. Dormieux, A. Molinari, and D. Kondo. Micromechanical approach to the behavior of poroelastic materials. *Journal of the Mechanics and Physics of Solids*, 50(10):2203–2231, 2002.

[48] G. Lu, G.Q. (M.) Lu, and Z.M. Xiao. Mechanical properties of porous materials. *Journal of Porous Materials*, 6:359–368, 1999.

[49] N. Abid, M. Mirkhalaf, and F. Barthelat. Discrete-element modeling of nacre-like materials: Effects of random microstructures on strain localization and mechanical performance. *Journal of the Mechanics and Physics of Solids*, 112:385–402, 2018.

[50] A. G. Evans, Z. Suo, R. Z. Wang, I. A. Aksay, M. Y. He, and J. W. Hutchinson. Model for the robust mechanical behavior of nacre. *Journal of Materials Research*, 16(09):2475–2484, 2001.

[51] K. Okumura and P.-G. de Gennes. Why is nacre strong? Elastic theory and fracture mechanics for biocomposites with stratified structures. *The European Physical Journal E*, 4(1):121–127, 2001.

[52] A. G. Evans, B.J. Dalgleish, M. Y. He, and J. W. Hutchinson. On crack path selection an the interface fracture energy in bimaterial systems. *Acta Metallurgica*, 37(12):3249–3254, 1989.

[53] M. Y. He, A. Bartlett, A. G. Evans, and J. W. Hutchinson. kinking of a crack out of an interface: role of in plane stress. *J. Am Ceram Soc.*, 74(4):767–771, 1991.

[54] M. Y. He, A. G. Evans, and J. W. Hutchinson. Crack Deflection at an interface between dissimilar elastic materials: role of residual stresses. *Int. J. Solids Structures*,
liptical inclusion on a straight crack. *International Journal of Fracture*, 46(1):71–79, November 1990.

[20] H. Gao. Fracture analysis of nonhomogeneous materials via a moduli-perturbation approach. *International Journal of Solids and Structures*, 27(13):1663–1682, 1991.

[21] S. Brach, M.Z. Hossain, B. Bourdin, and K. Bhat-tacharya. Anisotropy of the effective toughness of layered media. *Journal of the Mechanics and Physics of Solids*, 131:96–111, 2019.

[22] J.D. Clayton and J. Knap. A geometrically nonlinear phase field theory of brittle fracture. *International Journal of Fracture*, 189(2):139–148, 2014.

[23] T. Lacondemine, C. Roux-Langlois, and T. Rouxel. Role of poisson's ratio mismatch on the crack path in glass matrix particulate composites. *International Journal of Fracture*, 207(1):73–85, September 2017.

[24] L B Freund. *Dynamic Fracture Mechanics*. Cambridge University Press (UK), 1990.

[25] H. Henry. Study of the branching instability using a phase field model of inplane crack propagation. *Euro-physics Letters*, 83:16004, 2008.

[26] Z.P. Bazant and L.F. Estenssoro. Surface singularity and crack propagation. *International journal of solids and structure*, 15:405–426, 1979.

[27] H. Henry. Study of three-dimensional crack fronts under plane stress using a phase field model. *EPL (Europhysics Letters)*, 92(4):46002, nov 2010.

[28] A. Karma, D.A Kessler, and H. Levine. Phase-Field Model of Mode III Dynamic Fracture. *Phys. Rev. Lett.*, 87:45501, 2001.

[29] C. Kuhn, A. Schlüter, and R. Müller. On degradation functions in phase field fracture models. *Computational Materials Science*, 108:374–384, oct 2015.

[30] G.A. Francfort and J.-J. Marigo. Revisiting brittle fracture as an energy minimization problem. *Journal of the Mechanics and Physics of Solids*, 46(8):1319–1342, aug 1998.

[31] B. Bourdin, G. A. Francfort, and J.-J. Marigo. *The variational approach to fracture*. 2008.

[32] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations. *International Journal for Numerical Methods in Engineering*, 83(10):1273–1311, 2010.

[33] C. Maurini, B. Bourdin, G. Gauthier, and V. Lazarus. Crack patterns obtained by unidirectional drying of a colloidal suspension in a capillary tube: Experiments and numerical simulations using a two-dimensional variational approach. *Int. J. Frac.*, 184:75–91, 2013.

[34] E. Tanné, T. Li, B. Bourdin, J. J. Marigo, and C. Maurini. Crack nucleation in variational phase-field models of brittle fracture. *Journal of the Mechanics and Physics of Solids*, 110:80–99, jan 2018.

[35] Peter Galenko and David Jou. Diffuse-interface model for rapid phase transformations in nonequilibrium systems. *Phys. Rev. E*, 71:046125, Apr 2005.

[36] L. Verlet. Computer "experiments" on classical fluids. i. thermodynamical properties of lennard-jones molecules. *Phys. Rev.*, 159(1):98–103, July 1967.

[37] D. Dalmas, E. Barthel, and D. Vandembroucq. Crack front pinning by design in planar heterogeneous interfaces. *Journal of the Mechanics and Physics of Solids*, 57(3):446–457, 2009.

[38] H. Gao and J.R. Rice. A First-Order Perturbation Analysis of Crack Trapping by Arrays of Obstacles. *Journal of Applied Mechanics*, 56(4):828–836, 12 1989.

[39] S.V. Hainsworth. An environmental scanning electron microscopy investigation of fatigue crack initiation and propagation in elastomers. *Polymer Testing*, 26(1):60–70, 2007.

[40] S. Choukir and C.V. Singh. Role of topology in dictating the fracture toughness of mechanical metamaterials. *International Journal of Mechanical Sciences*, 241:107945, 2023.

[41] E. A. Jagla, François P. Landes, and Alberto Rosso. Viscoelastic effects in avalanche dynamics: A key to earthquake statistics. *Phys. Rev. Lett.*, 112:174301, Apr 2014.

[42] H. Henry. Limitations of the modelling of crack propagating through heterogeneous material using a phase field approach. *Theoretical and Applied Fracture Mechanics*, 104:102384, 2019.

[43] T.T. Nguyen, J. Yvonnet, Q.-Z. Zhu, M. Bornert, and C. Chateau. A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography. *Computer Methods in Applied Mechanics and Engineering*, 2015.

[44] T T Nguyen, J Yvonnet, M Bornert, and C Chateau. Initiation and propagation of complex 3D networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microCT experiments and phase field simulations. *Journal of the Mechanics and Physics of Solids*, 95:320–350, 2016.

[45] L. Schöller, D. Schneider, A. Prahs, and B. Nestler. Phase-field modeling of crack propagation based on multi-crack order parameters considering mechanical jump conditions. *PAMM*, 22(1):e202200039, 2023.

[46] A. H.-D. Cheng. *Poroelasticity*. Springer, 2016.

[47] L. Dormieux, A. Molinari, and D. Kondo. Micromechanical approach to the behavior of poroelastic materials. *Journal of the Mechanics and Physics of Solids*, 50(10):2203–2231, 2002.

[48] G. Lu, G.Q. (M.) Lu, and Z.M. Xiao. Mechanical properties of porous materials. *Journal of Porous Materials*, 6:359–368, 1999.

[49] N. Abid, M. Mirkhalaf, and F. Barthelat. Discrete-element modeling of nacre-like materials: Effects of random microstructures on strain localization and mechanical performance. *Journal of the Mechanics and Physics of Solids*, 112:385–402, 2018.

[50] A. G. Evans, Z. Suo, R. Z. Wang, I. A. Aksay, M. Y. He, and J. W. Hutchinson. Model for the robust mechanical behavior of nacre. *Journal of Materials Research*, 16(09):2475–2484, 2001.

[51] K. Okumura and P.-G. de Gennes. Why is nacre strong? Elastic theory and fracture mechanics for biocomposites with stratified structures. *The European Physical Journal E*, 4(1):121–127, 2001.

[52] A. G. Evans, B.J. Dalgleish, M. Y. He, and J. W. Hutchinson. On crack path selection an the interface fracture energy in bimaterial systems. *Acta Metallurgica*, 37(12):3249–3254, 1989.

[53] M. Y. He, A. Bartlett, A. G. Evans, and J. W. Hutchinson. kinking of a crack out of an interface: role of in plane stress. *J. Am Ceram Soc.*, 74(4):767–771, 1991.

[54] M. Y. He, A. G. Evans, and J. W. Hutchinson. Crack Deflection at an interface between dissimilar elastic materials: role of residual stresses. *Int. J. Solids Structures*,
31(24):3443-3455, 1994.

[55] M. Y. He, A. G. Evans, and J. W. Hutchinson. Interface cracking phenomena in constrained metal layers. *Acta Materialia*, 44(7):2963-2971, 996.

[56] M. Z. Hossain, C.-J. Hsueh, B. Bourdin, and K. Bhat-tacharya. Effective toughness of heterogeneous media. *Journal of the Mechanics and Physics of Solids*, 71:15-32, 2014.

[57] M. Lebihain, J.-B. Leblond, and L. Ponson. Effective toughness of periodic heterogeneous materials: the effect of out-of-plane excursions of cracks. *Journal of the Mechanics and Physics of Solids*, 137:103876, 2020.

[58] D. Schneider, E. Schoof, Y. Huang, M. Selzer, and B. Nestler. Phase-field modeling of crack propagation in multiphase systems. *Computer Methods in Applied Mechanics and Engineering*, 312:186-195, 2016.

[59] K.P. Herrmann and R. Wang. Interaction of an antiplane shear crack with a transforming inclusion. *Engineering Fracture Mechanics*, 50(3):361-368, 1995.

[60] R. R. Bhargava. A misfitting elastic inclusion in an infinite plane containing a crack. *Journal of Elasticity*, 7(2):201-211, April 1977.

[61] A.F. Bower and M. Ortiz. A three-dimensional analysis of crack trapping and bridging by tough particles. *Journal of the Mechanics and Physics of Solids*, 39(6):815-858, 1991.

[62] F. Erdogan, G. D. Gupta, and M. Ratwani. Interaction Between a Circular Inclusion and an Arbitrarily Oriented Crack. *Journal of Applied Mechanics*, 41(4):1007-1013, 12 1974.

[63] D. Brescakovic, M. Kegl, and O. Kolednik. Interaction of crack and hole: effects on crack trajectory, crack driving force and fracture toughness. *International Journal of Fracture*, 236(1):33-57, July 2022.

[64] E. M. Patton and M. H. Santare. The effect of a rigid elliptical inclusion on a straight crack. *International Journal of Fracture*, 46(1):71-79, November 1990.

[65] H. Gao. Fracture analysis of nonhomogeneous materials via a moduli-perturbation approach. *International Journal of Solids and Structures*, 27(13):1663-1682, 1991.

[66] S. Brach, M.Z. Hossain, B. Bourdin, and K. Bhat-tacharya. Anisotropy of the effective toughness of layered media. *Journal of the Mechanics and Physics of Solids*, 131:96-111, 2019.

[67] J.D. Clayton and J. Knap. A geometrically nonlinear phase field theory of brittle fracture. *International Journal of Fracture*, 189(2):139-148, 2014.

[68] T. Lacondemine, C. Roux-Langlois, and T. Rouxel. Role of poisson's ratio mismatch on the crack path in glass matrix particulate composites. *International Journal of Fracture*, 207(1):73-85, September 2017.

[69] L B Freund. *Dynamic Fracture Mechanics*. Cambridge University Press (UK), 1990.

[70] H. Henry. Study of the branching instability using a phase field model of inplane crack propagation. *Europhysics Letters*, 83:16004, 2008.

[71] Z.P. Bazant and L.F. Estenssoro. Surface singularity and crack propagation. *International journal of solids and structure*, 15:405-426, 1979.

[72] H. Henry. Study of three-dimensional crack fronts under plane stress using a phase field model. *EPL (Europhysics Letters)*, 92(4):46002, nov 2010.

[73] A. Karma, D.A Kessler, and H. Levine. Phase-Field Model of Mode III Dynamic Fracture. *Phys. Rev. Lett.*, 87:45501, 2001.

[74] C. Kuhn, A. Schlüter, and R. Müller. On degradation functions in phase field fracture models. *Computational Materials Science*, 108:374-384, oct 2015.

[75] G.A. Francfort and J.-J. Marigo. Revisiting brittle fracture as an energy minimization problem. *Journal of the Mechanics and Physics of Solids*, 46(8):1319-1342, aug 1998.

[76] B. Bourdin, G. A. Francfort, and J.-J. Marigo. *The variational approach to fracture*. 2008.

[77] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations. *International Journal for Numerical Methods in Engineering*, 83(10):1273-1311, 2010.

[78] C. Maurini, B. Bourdin, G. Gauthier, and V. Lazarus. Crack patterns obtained by unidirectional drying of a colloidal suspension in a capillary tube: Experiments and numerical simulations using a two-dimensional variational approach. *Int. J. Frac.*, 184:75-91, 2013.

[79] E. Tanné, T. Li, B. Bourdin, J. J. Marigo, and C. Maurini. Crack nucleation in variational phase-field models of brittle fracture. *Journal of the Mechanics and Physics of Solids*, 110:80-99, jan 2018.

[80] Peter Galenko and David Jou. Diffuse-interface model for rapid phase transformations in nonequilibrium systems. *Phys. Rev. E*, 71:046125, Apr 2005.

[81] L. Verlet. Computer "experiments" on classical fluids. i. thermodynamical properties of lennard-jones molecules. *Phys. Rev.*, 159(1):98-103, July 1967.

[82] D. Dalmas, E. Barthel, and D. Vandembroucq. Crack front pinning by design in planar heterogeneous interfaces. *Journal of the Mechanics and Physics of Solids*, 57(3):446-457, 2009.

[83] H. Gao and J.R. Rice. A First-Order Perturbation Analysis of Crack Trapping by Arrays of Obstacles. *Journal of Applied Mechanics*, 56(4):828-836, 12 1989.

[84] S.V. Hainsworth. An environmental scanning electron microscopy investigation of fatigue crack initiation and propagation in elastomers. *Polymer Testing*, 26(1):60-70, 2007.

[85] S. Choukir and C.V. Singh. Role of topology in dictating the fracture toughness of mechanical metamaterials. *International Journal of Mechanical Sciences*, 241:107945, 2023.

[86] E. A. Jagla, François P. Landes, and Alberto Rosso. Viscoclastic effects in avalanche dynamics: A key to earthquake statistics. *Phys. Rev. Lett.*, 112:174301, Apr 2014.

[87] H. Henry. Limitations of the modelling of crack propagating through heterogeneous material using a phase field approach. *Theoretical and Applied Fracture Mechanics*, 104:102384, 2019.

[88] T.T. Nguyen, J. Yvonnet, Q.-Z. Zhu, M. Bornert, and C. Chateau. A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography. *Computer Methods in Applied Mechanics and Engineering*, 2015.

[89] T T Nguyen, J Yvonnet, M Bornert, and C Chateau. Initiation and propagation of complex 3D networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microCT experiments and phase field simulations. *Journal of the Mechanics and Physics of Solids*, 95:320-350, 2016.
31(24):3443-3455, 1994.

[55] M. Y. He, A. G. Evans, and J. W. Hutchinson. Interface cracking phenomena in constrained metal layers. *Acta Materialia*, 44(7):2963-2971, 996.

[56] M. Z. Hossain, C.-J. Hsueh, B. Bourdin, and K. Bhat-tacharya. Effective toughness of heterogeneous media. *Journal of the Mechanics and Physics of Solids*, 71:15-32, 2014.

[57] M. Lebihain, J.-B. Leblond, and L. Ponson. Effective toughness of periodic heterogeneous materials: the effect of out-of-plane excursions of cracks. *Journal of the Mechanics and Physics of Solids*, 137:103876, 2020.

[58] D. Schneider, E. Schoof, Y. Huang, M. Selzer, and B. Nestler. Phase-field modeling of crack propagation in multiphase systems. *Computer Methods in Applied Mechanics and Engineering*, 312:186-195, 2016.

[59] K.P. Herrmann and R. Wang. Interaction of an antiplane shear crack with a transforming inclusion. *Engineering Fracture Mechanics*, 50(3):361-368, 1995.

[60] R. R. Bhargava. A misfitting elastic inclusion in an infinite plane containing a crack. *Journal of Elasticity*, 7(2):201-211, April 1977.

[61] A.F. Bower and M. Ortiz. A three-dimensional analysis of crack trapping and bridging by tough particles. *Journal of the Mechanics and Physics of Solids*, 39(6):815-858, 1991.

[62] F. Erdogan, G. D. Gupta, and M. Ratwani. Interaction Between a Circular Inclusion and an Arbitrarily Oriented Crack. *Journal of Applied Mechanics*, 41(4):1007-1013, 12 1974.

[63] D. Brescakovic, M. Kegl, and O. Kolednik. Interaction of crack and hole: effects on crack trajectory, crack driving force and fracture toughness. *International Journal of Fracture*, 236(1):33-57, July 2022.

[64] E. M. Patton and M. H. Santare. The effect of a rigid elliptical inclusion on a straight crack. *International Journal of Fracture*, 46(1):71-79, November 1990.

[65] H. Gao. Fracture analysis of nonhomogeneous materials via a moduli-perturbation approach. *International Journal of Solids and Structures*, 27(13):1663-1682, 1991.

[66] S. Brach, M.Z. Hossain, B. Bourdin, and K. Bhat-tacharya. Anisotropy of the effective toughness of layered media. *Journal of the Mechanics and Physics of Solids*, 131:96-111, 2019.

[67] J.D. Clayton and J. Knap. A geometrically nonlinear phase field theory of brittle fracture. *International Journal of Fracture*, 189(2):139-148, 2014.

[68] T. Lacondemine, C. Roux-Langlois, and T. Rouxel. Role of poisson's ratio mismatch on the crack path in glass matrix particulate composites. *International Journal of Fracture*, 207(1):73-85, September 2017.

[69] L B Freund. *Dynamic Fracture Mechanics*. Cambridge University Press (UK), 1990.

[70] H. Henry. Study of the branching instability using a phase field model of inplane crack propagation. *Europhysics Letters*, 83:16004, 2008.

[71] Z.P. Bazant and L.F. Estenssoro. Surface singularity and crack propagation. *International journal of solids and structure*, 15:405-426, 1979.

[72] H. Henry. Study of three-dimensional crack fronts under plane stress using a phase field model. *EPL (Europhysics Letters)*, 92(4):46002, nov 2010.

[73] A. Karma, D.A Kessler, and H. Levine. Phase-Field Model of Mode III Dynamic Fracture. *Phys. Rev. Lett.*, 87:45501, 2001.

[74] C. Kuhn, A. Schlüter, and R. Müller. On degradation functions in phase field fracture models. *Computational Materials Science*, 108:374-384, oct 2015.

[75] G.A. Francfort and J.-J. Marigo. Revisiting brittle fracture as an energy minimization problem. *Journal of the Mechanics and Physics of Solids*, 46(8):1319-1342, aug 1998.

[76] B. Bourdin, G. A. Francfort, and J.-J. Marigo. *The variational approach to fracture*. 2008.

[77] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations. *International Journal for Numerical Methods in Engineering*, 83(10):1273-1311, 2010.

[78] C. Maurini, B. Bourdin, G. Gauthier, and V. Lazarus. Crack patterns obtained by unidirectional drying of a colloidal suspension in a capillary tube: Experiments and numerical simulations using a two-dimensional variational approach. *Int. J. Frac.*, 184:75-91, 2013.

[79] E. Tanné, T. Li, B. Bourdin, J. J. Marigo, and C. Maurini. Crack nucleation in variational phase-field models of brittle fracture. *Journal of the Mechanics and Physics of Solids*, 110:80-99, jan 2018.

[80] Peter Galenko and David Jou. Diffuse-interface model for rapid phase transformations in nonequilibrium systems. *Phys. Rev. E*, 71:046125, Apr 2005.

[81] L. Verlet. Computer "experiments" on classical fluids. i. thermodynamical properties of lennard-jones molecules. *Phys. Rev.*, 159(1):98-103, July 1967.

[82] D. Dalmas, E. Barthel, and D. Vandembroucq. Crack front pinning by design in planar heterogeneous interfaces. *Journal of the Mechanics and Physics of Solids*, 57(3):446-457, 2009.

[83] H. Gao and J.R. Rice. A First-Order Perturbation Analysis of Crack Trapping by Arrays of Obstacles. *Journal of Applied Mechanics*, 56(4):828-836, 12 1989.

[84] S.V. Hainsworth. An environmental scanning electron microscopy investigation of fatigue crack initiation and propagation in elastomers. *Polymer Testing*, 26(1):60-70, 2007.

[85] S. Choukir and C.V. Singh. Role of topology in dictating the fracture toughness of mechanical metamaterials. *International Journal of Mechanical Sciences*, 241:107945, 2023.

[86] E. A. Jagla, François P. Landes, and Alberto Rosso. Viscoclastic effects in avalanche dynamics: A key to earthquake statistics. *Phys. Rev. Lett.*, 112:174301, Apr 2014.

[87] H. Henry. Limitations of the modelling of crack propagating through heterogeneous material using a phase field approach. *Theoretical and Applied Fracture Mechanics*, 104:102384, 2019.

[88] T.T. Nguyen, J. Yvonnet, Q.-Z. Zhu, M. Bornert, and C. Chateau. A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography. *Computer Methods in Applied Mechanics and Engineering*, 2015.

[89] T T Nguyen, J Yvonnet, M Bornert, and C Chateau. Initiation and propagation of complex 3D networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microCT experiments and phase field simulations. *Journal of the Mechanics and Physics of Solids*, 95:320-350, 2016.
[90] L. Schöller, D. Schneider, A. Prahs, and B. Nestler. Phase-field modeling of crack propagation based on

multi-crack order parameters considering mechanical jump conditions. *PAMM*, 22(1):e202200039, 2023.
[90] L. Schöller, D. Schneider, A. Prahs, and B. Nestler. Phase-field modeling of crack propagation based on

multi-crack order parameters considering mechanical jump conditions. *PAMM*, 22(1):e202200039, 2023.