Peridynamics simulation of the comminution of particles containing micro-craks

. In this study, we rely on a ’bond-based’ peridynamic approach to investigate the strength and failure of 2D particles containing a collection of 1D microcracks. The mechanical tests were performed on disks under diametral compression. In an extensive parametric study, the distribution of microcracks was varied for di ﬀ erent particle sizes. The evolution of yield stress with diameter and the probability of failure in terms of Weibull distributions are investigated in detail. Finally, by means of a ﬂoodﬁll algorithm, we analyze the variation of the mean fragment size as a function of the density of defects.


Introduction
Comminution consists in reducing the average size of particles by different processes such as crushing, grinding or cutting.It is employed in many industrial applications and consumes large amounts of energy in applications such as cement manufacturing [1] and cereal milling [2].Several models, mostly based on empirical parameters, have been developed to estimate the energy consumption of such processes.Population balance models [3,4] involve the probability of failure and distribution of fragments obtained for a single particle of a given size.Brazilian tests (diametral compression) have been used for the determination of particle strength [5,6] and in studies of the influence of the microstructure on the fracture and fragmentation [7,8].Numerical simulations based on the Discrete Element Method (DEM) have also been used for the investigation of particle crushing by modeling each grain as a cohesive agglomerate of smaller (primary) particles [9][10][11][12].Such simulations can now reproduce the evolution of cracks and complex geometries of fragments in 2D [13,14] and 3D [15][16][17].But they fundamentally replace a continuum by a discrete medium with pre-existing potential fragments.
In this paper, we use a peridynamic approach applied to particles which can break into fragments of arbitrary sizes and shapes (up to the discretization limit) in the presence of microcracks.We investigate the crack patterns of 2D rounded particles containing defects of given density and distribution loaded between two horizontal plates in compression.We are interested in their yield stress and distribution of fragments as a function of the density of defects for a range of particle sizes.e-mail: nicolas.blanc@supagro.inra.fr

Bond-based peridynamics
The peridynamics method was introduced by Silling [18] as an alternative to the description of classical continuum mechanics by differential equations.Peridynamics has been successfully used for modeling the fracture of materials such as concrete [19], nickel nanowires [20], wood [21]. . . .We use a simple implementation of peridynamics called bond-based approach which describes the mechanics of deformable bodies as a mass-spring system.A major interest of this model is its weak dependence on the geometry of the mesh contrary to other bond-based methods such as Fuse Models [22] and Lattice Element Method (LEM) [23].The reason is that the computation of mechanical properties relies on a multinodal interactions within a distance h called the horizon.As other bondbased models, this method is allowing inclusion of discontinuities such as pores, cracks or stiffness and damage gradients without requiring special case-sensitive treatment of each problem [18].
The bonds are 1D elements characterized by a bond stiffness k = c (δx) 4 , where c is the elastic micromodulus and δx the space step.The damage is introduced using a yield elongation s 0 above which the bond is broken.A mass m i is associated to each node i.The stress tensor is upscaled from the forces connecting i to neighboring nodes j using the virial expansion where V i is a volume attributed to the node i and F i j and the r i j are the force and length of the bond i j, respectively.The Young modulus E = αc of the material is proportional to the micromodulus with α = πh 3 (1−ν)/6 [24], where ν is Poisson's ratio.For isotropic linear elastic materials with central interactions, the value of Poisson's ratio is 1/3 in 2D.The fracture energy G is easily deduced from the crit- ical elongation s 0 as G = (9Ehs 0 2 )/4π [24].The computation of the evolution of the system is similar to Molecular Dynamics or Discrete Element Method simulations.We use the velocity-Verlet algorithm [25] to integrate the equations of motion with a time step δt.

Sample building
The particles are 2D disks meshed by a regular grid.We consider defects in the form of 1D linear segments.For all bonds crossing these zones we set the elastic modulus to a value below E. The defects are randomly distributed within the particle with number density n.Their lengths and Young moduli follow Gaussian distributions: ¯ ∼ N(μ¯ , σ¯ ), and Ē ∼ N(μ Ē , σ Ē ).
The particle is then subjected to a quasistatic diametral compression test.The bottom plate is kept fixed while the upper plate moves downward.The contacts between the particle and the plates are modeled by adding a normal repulsive force to the particle's nodes touching the plates.A regularized Coulomb friction law with a friction coefficient of 0.5 and a viscosity of 0.8 was used.The compression continues until the primary cracks fully propagate and before secondary fragments appear.As we rely on an explicit resolution of the equations of dynamics, we use a global viscous damping in order to dissipate lowfrequency elastic waves.
The fracture of particle leads to the failure of bonds.From these failures, we determine at each node a damage parameter d i defined as the ratio of the number of broken bonds to the initial number of bonds attached to this node.Fig. 1a shows this damage field after fracture.The color level indicates the magnitude of the local damage, which has its maximum value in the middle of the cracks.As the peridynamics is a non-local approach, the crack-path dependence on the grid orientation is limited, compared with purely local lattice-based frameworks, like LEM.
To identify the fragments, we use a Floodfill algorithm (see Fig. 1b, c and d) [26].This algorithm paints connected zones made of the pixels which have damage values lower than a given threshold T .This threshold affects the width of the cracks and thus the number of potential fragments produced.Fig. 2 shows the evolution of the number of fragments as a function of T .In the following we set T to 30%, which corresponds to a value for which the number of fragments is roughly constant.

Parametric study
The main objective of this parametric study is to investigate the influence of the particle size on the fracture process.It depends on the ratio of particle diameter D to the characteristic size μ¯ of the defects.For a reference diameter D 0 , we introduce n0 = 100 defects uniformly distributed in a square area (see Fig. 3).The particle extracted from this area has thus an average number n0 × 4 π 78 of defects.The average size of the defects was set to μ ¯ 0 D 0 /30 with a standard deviation of σ ¯ 0 D 0 /200.The particle was then meshed with a spatial resolution of  256 nodes per particle diameter.To reduce computational costs, we did not directly vary the particle size but the values of parameters related to the defects: In all cases, the Young moduli of the defects are chosen to be narrowly distributed with an average μ Ē /E= 0.2 and a standard deviation σ Ē = 0.005.Hence, the variability of particle strength mainly depends on the distribution of the lengths and the random positions of the defects.Table 1 shows the seven values of D/D 0 used in our study and the corresponding defect parameters.In order to obtain meaningful statistics on particle fracture, we broke 100 particles for each diameter.

Results and discussion
Figure 4 shows the average size of the fragments as a function of dimensionless particle diameter.For a diameter above 0.2, the mean size of the fragments decreases with particle size.The increase in the number of defects with particle size leads to an increase of potential crack paths and thus to a larger number of fragments.For the smallest size (diameter of 0.005), the density of defects is too small so that the particle breaks homogeneously into numerous small pieces.Assuming that the fracture occurs for the maximum value of tensile stress in the middle of the particle [27], we get: where F is the yield force at failure.In figure 5 we plot as a function of D/D 0 , the average value of σ t /σ th where σ th is the yield stress in tension of the material.We see that the strength of the particle decreases as the diameter increases, and follows a power law with an exponent of -0.101.This behavior is in agreement with experimental results [27].
In order to study the dispersion of the values of the yield stress, we consider the probability for a particle to survive under a given loading.For N particles, we assign the probability P S (i) = 1 − i N to the i th smallest yield stress value.This distribution is then fitted by the Weibull function [28]:  where σ w corresponds to the stress for which 63% of particles are broken and m is the Weibull modulus characterizing the width of the distribution.Figure 6 shows ln(ln(1/P S )) as a function of ln(σ t /σ th ).If the data collapse on a straight line, the slope equals the Weibull modulus of the distribution.Note that a value of 0 on the x axis corresponds to a yield stress equal to σ th and thus to a material without defect.Except for smallest particle diameters (0.05 and 0.2), the distributions are well fitted by the Weibull law.The Weibull modulus globally increases with particle size (m = 4.4 to m = 41).As the number of defects is increasing with diameter, the probability of finding similar critical defects between particles increases.As a consequence, the distributions get narrower.On the contrary, the dispersion is smaller for 0.05 compared to 0.2.As the cracks are initiated in tension, the presence of a defect in the vicinity of the particle center has a dramatic impact.Because of the very low density of defects, many particles have no defect at their center.Thus, the probability distribution becomes narrow as the particles yield stress reaches the theoretical value of the material.

Conclusion
In this paper, we developed a methodology to investigate the fracture of 2D particles containing defects.More than 700 compression tests were performed with a peridynamic bond-based approach to characterize the evolution of yield stress and average fragment size as a function of particle diameter.We found that the Weibull modulus evolves with particle size and the average yield stress decreases as a power law with diameter.
In future work, we would like to analyze in detail the size distributions of the fragments.We would also like to analyze crack patterns.The simulations allow us also to evaluate the distributions of crack lengths and the interactions between the growing crack and defects.As the cracks tend to cross the defects, this leads to a lower number of defects in the fragments.We also plan to investigate the influence of defect lengths (Gaussian or power law distribution) on the probability distributions of the yield stress.

Figure 1 .
Figure 1.a): Damage field after particle fracture.The color level represents the magnitude of the local damage.b), c) and d): Results of the Floodfill algorithm with T =1%, T =33% and T =66% respectively.

Figure 2 .
Figure 2. Average number of fragments of 100 particles as a function of damage threshold.

Figure 3 .
Figure 3. Population of defects in the particle for different diameters.From left to right: D/D 0 =2.85,D/D 0 =1 and D/D 0 =0.2.

Figure 4 .
Figure 4. Mean fragment size normalized by particle area as a function of dimensionless diameter.

Figure 5 .
Figure 5. Average yield stress in tension normalized by the yield stress of the material without defects.

Figure 6 .
Figure 6.Surviaval probability of particles as a function of the applied stress for different values of relative particle size D/D 0 .The straight lines represent Weibull fits to the data.

Table 1 .
Parameters of defect distribution as a function of dimensionless particle diameter.