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

. 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 identiﬁed depending on the toughness ratio.


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. e-mail: xavier.frank@inra.fr

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 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 Under a deformation B 0 → B t the material points x undergo a displacement u x, t . Thus, a bond ξ = x − x becomes after deformation ξ+ η where η = u x , t − u x, t is the relative displacement ( Figure 1).
The force exerted on a material point x results from the contribution of all bonds in H( x). Here, we use the so-called bond-based approach in which we only consider radial pairwise forces f that depend on both ξ and η [19]. The value of the Poisson ratio of a homogeneous sample is fixed to ν = 1/3 [8]. The forces can be expressed as where f ξ , η is the bond force magnitude.
In this work, we use harmonic pairwise forces defined as: where s is the bond elongation and c is the elastic micromodulus. The value of the effective Young modulus of the material is given by [9]: which linearly depends on the micromodulus c. The damage of the material is introduced using a history-dependent scalar function μ as a prefactor modifying c to μ ξ, t c, and defined as 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 ξ , η = 0. The fracture energy G is given by [9] G = (9Ehs 0 2 )/4π For a heterogeneous material, each phase α is characterized by its elastic modulus E α and its toughness K α = √ E α G α , where G α is the fracture energy of the phase α.
If a bond connects two nodes i and j in the same phase α the mechanical properties of the bond are directly deduced from Equations 4 and 6. If the bond connects two different phases α and β, we choose E αβ = min (E α , E β ) and we set G αβ so that K αβ = 1 2 (K α + K β ). More generally, the interphase properties can be independent of bulk properties.

Implementation
The discretization of the material is performed in 2D according to a regular grid of nodes x i with a spatial resolution δx. A mass m i = (δx) 2 ρ x i is attributed to each node. The mechanical system is similar to a mass-spring lattice in which each bond in H x i represents a spring. The evolution of the system is computed by the discretized equations of motion: 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 quasistatic conditions. For the dynamic part we use the classical velocity-Verlet algorithm [20] with a time step δ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 κ. For a sufficiently small value, the process is quasistatic but the computation time increases dramatically. We fixed the value of κ as the energy required to propagate a crack of distance equal to 1% of the sample width.

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 ϕ were used (0.3, 0.5 and 0.7) with six values of grainmatrix toughness ratio ψ = 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 parallezied 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( x) associated with the nodes. When a bond connecting two nodes x i and x j breaks, the values of n( x i ) and n( 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 developped and we compute the grain damage parameter d = N g /N tot (8) where N g = x i ∈G n( x i ), G is the subset of B 0 corresponding to the grain phase and N tot = x i ∈B 0 n( x i ).  Figure 3 shows two typical examples of crack patterns for ϕ = 0.5 and a toughness ratio ψ of 0.25 and 4.00, corresponding respectively to soft and hard inclusions. In these figures, we clearly observe different failure behaviors. For ψ = 0.25 the failure crosses the particles whereas for ψ = 4 the crack remains inside the matrix and bypasses the particles. This typical behavior was already observed for an single inclusion in a cement matrix [22].

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 ϕ,ψ . Figure 4 shows for three values of the grain volume fraction ϕ, the average grain damage parameter d (Equation 8) as a function of toughness ratio ψ. The amount of broken bonds in the grain phase increases with ϕ. It is interesting to note that  the value of ψ below which no failure occurs in grains is the same (ψ 3) whatever the value of ϕ.
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 ψ = 0.5, and ψ = 2.0, respectively. As the grains do not have the same size, each crack trajectory was normalized by the corresponding 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. We see that, for a low value of the grain-matrix toughness ratio ψ, the grains tend to focus the cracks like a convergent optical lens. On the contrary, for ψ > 1, the grains deflect the cracks. This convergent-divergent transition occurs for a value of ψ = 1 which corresponds to a regime of straight failure.

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