# Numerical study of the failure of materials embedding soft to hard particles

Xavier Frank, Jean-Yves Delenne, Farhang Radjai

► To cite this version:

Xavier Frank, Jean-Yves Delenne, Farhang Radjai. Numerical study of the failure of materials embedding soft to hard particles. International Conference on Micromechanics of Granular Media (Powders & Grains), Jul 2017, Montpellier, France. 10.1051/epjconf/201714002029. hal-01595485

HAL Id: hal-01595485

https://hal.science/hal-01595485

Submitted on 26 Sep 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical study of the failure of materials embedding soft to hard particles

Xavier Frank<fnref n="1" />,*, Jean-Yves Delenne<fnref n="1" />, and Farhang Radjai<fnref n="2" />,<fnref n="3" />

**Abstract.** In this study, we use a bond-based peridynamic approach to investigate the mechanical strength and cracking of composite materials with spherical inclusions. The total volume fraction of particles and the particle-matrix toughness ratio were varied to cover a range of soft to hard inclusions. The mean particle damage was characterized together with crack patterns at a sub-particle scale. Three types of crack patterns are identified depending on the toughness ratio.

## 1 Introduction

Damage and failure of materials including a granular phase is crucial in numerous natural processes and industrial applications. The case of tough inclusions was studied in detail by several authors for various applications such as the failure of concrete [1], the milling of wheat [2] (whose internal micro-structure can be seen as a cemented granular material composed of starch granules glued by a protein matrix), the strength of particle-enriched composites [3]... However, in many examples the particles can have a lower toughness than the matrix. This is the case of some biocomposites used in packaging [4] and in food products such as granola bars or rice chocolate in which the taste largely depends on the textural properties [5].

In this paper, we use a bond-based peridynamic approach to investigate this transition from fragile to tough inclusions in particle-laden composites. A major advantage of this approach is its versatile nature in application to discontinuities such as pores, cracks and jumps in mechanical properties... allowing one to simulate crack patterns [6, 7]. The bond-based peridynamic approach has been successfully applied for modeling elastic brittle materials [7] including the dynamic effects and fracture (for example, crack patterns in concrete [8], dynamic cracks branching [9], fracture of polycrystals [10], nanoscale mechanics [11, 12]...). In this paper, we investigate the crack patterns in 2D composites including a suspension of disk particles. We performed an extensive parametric study by varying the density of particles in the matrix and the yield stress of particles in order to simulate a range of very soft to very hard inclusions.

## 2 Bond-based peridynamics

In recent years, various numerical approaches have been designed for the computation of crack patterns such as cohesive zone models [13], the cohesive discrete elements method [14, 15], the contact dynamics-based bonded-cell models [16] and the lattice element method [17, 18]. The discrete peridynamic approach was introduced by Silling [6]. In this approach the mechanical behavior is based on integro-differential equations rather than partial differential equations. The main advantage is to simplify the treatment of the discontinuities in the material. We use this method for its versatile and simple implementation to investigate the cracking of composites with soft and hard inclusions. In the following, we describe the basics of the method in two dimensions.

Lets consider a deformable body $B_t$ defined at the time $t$ as a subset of $\mathbb{R}^2$. For each point in this domain, the strain, stress, damage... can be directly upscaled from the behaviour of bonds defined at the microscale. Figure 1 shows, for the reference configuration $B_0$ (at time $t = 0$) a bond $\xi = \vec{x}' - \vec{x}$ where $\vec{x}$ and $\vec{x}' \in B_0$. At each point $\vec{x}$ a set of bonds $H(\vec{x}) = \{\xi \in \mathbb{R}^2 | \vec{x} + \xi \in B_0 \wedge \|\xi\| < h\}$ is defined, where $h$ is the cutoff distance usually called horizon in peridynamics.

Under a deformation $B_0 \to B_t$ the material points $\vec{x}$ undergo a displacement $\vec{u}(\vec{x}, t)$. Thus, a bond $\xi = \vec{x}' - \vec{x}$ becomes after deformation $\xi+\vec{\eta}$ where $\vec{\eta} = \vec{u}(\vec{x}', t) - \vec{u}(\vec{x}, t)$ is the relative displacement (Figure 1).

The force exerted on a material point $\vec{x}$ results from the contribution of all bonds in $H(\vec{x})$. Here, we use the so-called bond-based approach in which we only consider radial pairwise forces $\vec{f}$ that depend on both $\xi$ and $\vec{\eta}$ [19]. The value of the Poisson ratio of a homogeneous sample

<footnote n="1">¹IATE, CIRAD, INRA, Montpellier SupAgro, Université de Montpellier, F-34060, Montpellier, France</footnote>

<footnote n="2">²Laboratoire de Mécanique et Génie Civil (LMGC), Université de Montpellier, CNRS, Montpellier, France</footnote>

<footnote n="3">³*MSE*², UMI 3466 CNRS-MIT, CEE, MIT, 77 Massachusetts Avenue, Cambridge CA 02139, USA</footnote>

<footnote>*e-mail: xavier.frank@inra.fr</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Initial configuration $B_0$: schematic representation of a bond $\vec{\xi}$ and the bond family $H(\vec{x})$ of a point $\vec{x}$ in a peridynamic body, definition of the horizon $h$. Schematic representation of the deformation of a peridynamic body from initial configuration $B_0$ to $B_t$ and definition of both displacement $\vec{u}$ and relative displacement $\vec{\eta}$.</figcaption></figure>

is fixed to $v = 1/3$ [8]. The forces can be expressed as

$$ \vec{f}(\vec{\xi}, \vec{\eta}) = f(\vec{\xi}, \vec{\eta}) \frac{\vec{\xi} + \vec{\eta}}{\|\vec{\xi} + \vec{\eta}\|} \quad (1) $$

where $f(\vec{\xi}, \vec{\eta})$ is the bond force magnitude.

In this work, we use harmonic pairwise forces defined as:

$$ s(\vec{\xi}, \vec{\eta}) = \frac{\|\vec{\xi} + \vec{\eta}\| - \|\vec{\xi}\|}{\|\vec{\xi}\|} \quad (2) $$

$$ f(\vec{\xi}, \vec{\eta}) = c s(\vec{\xi}, \vec{\eta}) \quad (3) $$

where $s$ is the bond elongation and $c$ is the elastic micro-modulus. The value of the effective Young modulus of the material is given by [9]:

$$ E = c(\pi h^3 (1 - v)/6) \quad (4) $$

which linearly depends on the micromodulus $c$. The damage of the material is introduced using a history-dependent scalar function $\mu$ as a prefactor modifying $c$ to $\mu(\vec{\xi}, t)c$, and defined as

$$ \mu(\vec{\xi}, t) = \begin{cases} 1 & s(\vec{\xi}, \vec{\eta}(t')) < s_0 \quad \forall t' < t \\ 0 & \exists t' < t \mid s(\vec{\xi}, \vec{\eta}(t')) \ge s_0 \end{cases} \quad (5) $$

where $s_0$ is the critical elongation of the bond. In other words, when at a time $t'$ the elongation $s$ exceeds $s_0$, the bond is definitely broken and $f(\vec{\xi}, \vec{\eta}) = 0$. The fracture energy $G$ is given by [9]

$$ G = (9Eh s_0^2)/4\pi \quad (6) $$

For a heterogeneous material, each phase $\alpha$ is characterized by its elastic modulus $E_\alpha$ and its toughness $K_\alpha = \sqrt{E_\alpha G_\alpha}$, where $G_\alpha$ is the fracture energy of the phase $\alpha$.

If a bond connects two nodes $i$ and $j$ in the same phase $\alpha$ the mechanical properties of the bond are directly deduced from Equations 4 and 6. If the bond connects two different phases $\alpha$ and $\beta$, we choose $E_{\alpha\beta} = \min(E_\alpha, E_\beta)$ and we set $G_{\alpha\beta}$ so that $K_{\alpha\beta} = \frac{1}{2}(K_\alpha + K_\beta)$. More generally, the inter-phase properties can be independent of bulk properties.

### 3 Implementation

The discretization of the material is performed in 2D according to a regular grid of nodes $\vec{x}_i$ with a spatial resolution $\delta x$. A mass $m_i = (\delta x)^2 \rho(\vec{x}_i)$ is attributed to each node. The mechanical system is similar to a mass-spring lattice in which each bond in $H(\vec{x}_i)$ represents a spring. The evolution of the system is computed by the discretized equations of motion:

$$ m_i \ddot{u}_i(t) = \sum_{\vec{x}_j \in H(\vec{x}_i)} k s(\vec{\eta}_{ij}, \vec{\xi}_{ij}) \frac{\vec{\xi}_{ij} + \vec{\eta}_{ij}}{\|\vec{\xi}_{ij} + \vec{\eta}_{ij}\|} + \vec{b}_i(t) \quad (7) $$

where $\vec{u}_i(t) = \vec{u}(\vec{x}_i, t)$, $\vec{\xi}_{ij} = \vec{x}_j - \vec{x}_i$, $\vec{\eta}_{ij} = \vec{u}_j - \vec{u}_i$, $\vec{b}_i(t) = (\delta x)^2 \vec{b}(\vec{x}_i, t)$ and $k = c(\delta x)^4$.

External forces or displacements are applied at the boundaries of the samples. The stress-strain evolution of the system can be computed either in dynamics or quasi-static conditions. For the dynamic part we use the classical velocity-Verlet algorithm [20] with a time step $\delta t$. In order to efficiently damp low frequency modes, we use a global viscous term at all nodes. To ensure the quasi-static evolution of the system, we use a snap-back algorithm which consists in relaxing the loading strain periodically when the cumulative energy due to crack growth reaches a given value $\kappa$. For a sufficiently small value, the process is quasi-static but the computation time increases dramatically. We fixed the value of $\kappa$ as the energy required to propagate a crack of distance equal to 1% of the sample width.

### 4 Sample building & mechanical tests

We prepared a series of cemented granular samples. These samples are composed of a collection of disks introduced into a matrix. For the granular phase, six samples were first prepared using a DEM code [21] and packed under a weak confining pressure. The particle radii are slightly dispersed in order to avoid local crystallization. These samples are then meshed using a rectilinear grid with a spatial resolution of approximately 50 nodes in the grain's mean diameter and the left boundary of each sample was notched to initiate crack propagation. To study the influence of particle density, we apply a geometrical shrinkage factor to each grain. Three particle solid fractions $\varphi$ were used (0.3, 0.5 and 0.7) with six values of grain-matrix toughness ratio $\psi = K_g/K_m$ (0.25, 0.375, 0.5, 2.0, 3.0 and 4.0). Finally, 108 samples were meshed. As the computation effort is relatively high (with approximately 500 000 nodes and 7 000 000 bonds for each sample), we paralleled the code using Message Passing Interface (MPI).
All the samples were subjected to quasi-static tensile tests until full failure occurs. Figure 2 shows an example of the snap-back procedure and the stress-strain envelope. For each sample we consider the damage after full failure. We consider a scalar field $n(\vec{x})$ associated with the nodes. When a bond connecting two nodes $\vec{x}_i$ and $\vec{x}_j$ breaks, the values of $n(\vec{x}_i)$ and $n(\vec{x}_j)$ are incremented. The damage at a node can then evolve between 0 and the maximum number of bonds connected to a node which depends on the horizon. In the following, we consider only the total damage at the end of the simulation when the failure is fully developed and we compute the grain damage parameter

$$d = N_g/N_{tot} \qquad (8)$$

where $N_g = \sum_{x_i \in G} n(\vec{x}_i)$, G is the subset of $B_0$ corresponding to the grain phase and $N_{tot} = \sum_{x_i \in B_0} n(\vec{x}_i)$.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** (color online) Evolution of the stress normalized by the yield stress of the matrix $\sigma_m$ as a function of the strain $\epsilon$ for a sample with grain volume fraction $\varphi = 0.7$ and a toughness ratio $\psi = 3.00$.</figcaption></figure>

Figure 3 shows two typical examples of crack patterns for $\varphi = 0.5$ and a toughness ratio $\psi$ of 0.25 and 4.00, corresponding respectively to soft and hard inclusions. In these figures, we clearly observe different failure behaviors. For $\psi = 0.25$ the failure crosses the particles whereas for $\psi = 4$ the crack remains inside the matrix and bypasses the particles. This typical behavior was already observed for a single inclusion in a cement matrix [22].

# 5 Parametric study

In this section we discuss the results of the parametric study. For a quantitative description, we computed the total number of broken bonds in each phase cumulated on all samples for each parameter set $\varphi, \psi$. Figure 4 shows for three values of the grain volume fraction $\varphi$, the average grain damage parameter $d$ (Equation 8) as a function of toughness ratio $\psi$. The amount of broken bonds in the grain phase increases with $\varphi$. It is interesting to note that

<figure><img src="image_4.png" /><figcaption>**Figure 3.** (color online) Crack path for grain volume fraction $\varphi = 0.5$. A: $\psi = 0.25$, B: $\psi = 4.0$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** (color online) Grain damage parameter $d$ corresponding to the average fraction of broken bonds inside grains (cumulated on all samples) as a function of toughness ratio $\psi$ for three values of the grain volume fraction $\varphi$.</figcaption></figure>

the value of $\psi$ below which no failure occurs in grains is the same ($\psi \approx 3$) whatever the value of $\varphi$.

For the same samples, Figure 5a and Figure 5b show the superposition of all crack paths in the vicinity and inside the damaged grains for a toughness ratio $\psi = 0.5$, and $\psi = 2.0$, respectively. As the grains do not have the same
size, each crack trajectory was normalized by the corre-
sponding grain radii. Note also that in the simulations, as
the failure is initiated by a notch on the left hand side of
the sample, the propagation occurs from left to right.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** (color online) Superposition of all crack patterns (crack propagation from the left to the right) in the neighborhood of broken grains for grain volume fraction 0.4. Toughness ratio is A: 0.5, and B:2.0.</figcaption></figure>

We see that, for a low value of the grain-matrix tough-
ness ratio ψ, the grains tend to focus the cracks like a con-
vergent optical lens. On the contrary, for ψ > 1, the grains
deflect the cracks. This convergent-divergent transition oc-
curs for a value of ψ = 1 which corresponds to a regime of
straight failure.

**6 Conclusion**

In this paper we introduced a bond-based quasi-static peri-
dynamic approach for the simulation of granular compos-
ites. This method is quite well suited to the computation
of damage and crack patterns in highly heterogeneous me-
dia. An extensive parametric study was performed to in-
vestigate the effects of toughness ratio and grain volume
fraction on the damage and fracture behavior across the
grain-matrix interfaces.

Acknowledgment

We are grateful to the genotoul bioinformatics platform
Toulouse Midi-Pyrenees (Bioinfo Genotoul) for providing
computing resources.

**References**

[1] M. Issa, M. Issa, M. Islam, A. Chudnovsky, Eng. Fract. Mech. **70**, 125 (2003)

[2] V. Topin, F. Radjaï, J.Y. Delenne, A. Sadoudi, F. Ma-
bille, J. Cereal Sci. **47**, 347 (2008)

[3] R. Eckner, M. Krampf, C. Segel, L. Krüger, Mechan-
ics of Composite Materials **51**, 707 (2016)

[4] M.A. Berthet, H. Angellier-Coussy, V. Guillard, N. Gontard, Journal of Applied Polymer Science **133**, 42528 (2016)

[5] T. van Vliet, Rheology and fracture mechanics of foods (CRC Press, 2013)

[6] S. Silling, J. Mech. Phys. Solids **48**, 175 (2000)

[7] P. Perré, G. Almeida, M. Ayouz, X. Frank, Annals of Forest Science **73**, 147 (2016)

[8] W. Gerstle, N. Sau, S. Silling, Nuclear Eng. & De-
sign **237**, 1250 (2007)

[9] Y. Ha, F. Bobaru, Eng. Fract. Mech. **78**, 1156 (2011)

[10] E. Askari, F. Bobaru, R. Lehoucq, M. Parks, S. Silling, O. Weckner, J. of Phys.: Conf. Series **125**, 012078 (2008)

[11] F. Bobaru, Modelling Simul. Mater. Sci. Eng. **15**, 397 (2007)

[12] E. Celik, I. Guven, E. Madenci, Theor. Appl. Fract. Mech. **55**, 185 (2011)

[13] X. Ren, J. Li, Phys. Rev. E **85**, 055102 (2012)

[14] J.Y. Delenne, M.S. El Youssoufi, F. Cherblanc, J.C. Bénet, International Journal for Numerical and Analytical Methods in Geomechanics **28**, 1577 (2004)

[15] H.A. Carmona, A. Guimarães, J. Andrade, I.N. Jr., F. Wittel, H. Herrmann, Phys. Rev. E **91**, 012402 (2015)

[16] D.H. Nguyen, E. Azéma, P. Sornay, F. Radjaï, Phys. Rev. E **91**, 022203 (2015)

[17] V. Topin, J.Y. Delenne, F. Radjaï, L. Brendel, F. Ma-
bille, Eur. Phys. J. E **23**, 413 (2007)

[18] R. Affes, J.Y. Delenne, Y. Monerie, F. Radjaï, V. Topin, Eur. Phys. J. E **35**, 117 (2012)

[19] S. Silling, M. Epton, O. Weckner, J. Xu, E. Askari, J. Elasticity **88**, 151 (2007)

[20] M. Allen, D. Tildesley, *Computer Simulation of Liquids* (Oxford University Press, Oxford, 1986)

[21] F. Radjaï, F. Dubois, *Discrete Numerical Modeling of Granular Materials* (Wiley ISTE, 2011)

[22] O. Buyukozturk, B. Hearing, Int. J. Solids Struct. **35**, 4055 (1998)