Compaction of granular materials composed of deformable particles

. In soft particle materials such as metallic powders the particles can undergo large deformations without rupture. The large elastic or plastic deformations of the particles are expected to strongly a ﬀ ect the mechanical properties of these materials compared to hard particle materials more often considered in research on granular materials. Herein, two numerical approaches are proposed for the simulation of soft granular systems: (i) an implicit formulation of the Material Point Method (MPM) combined with the Contact Dynamics (CD) method to deal with contact interactions, and (i) Bonded Particle Model (BPM), in which each deformable particle is modeled as an aggregate of rigid primary particles using the CD method. These two approaches allow us to simulate the compaction of an assembly of elastic or plastic particles. By analyzing the uniaxial compaction of 2D soft particle packings, we investigate the e ﬀ ects of particle shape change on the stress-strain relationship and volume change behavior as well as the evolution of the microstructure.


Introduction
Granular materials made of soft deformable particles represent a broad class of materials. Besides granular materials such as metallic powders, in which the particles can undergo large plastic deformations, a large number of cosmetic and food products, emulsions and biological cells may be included in this class. The deformation of the individual soft particles has been a subject of experimental and theoretical investigations at low level of deformation and more recently for large deformations [1,2]. Although the behavior of individual particles has been studied experimentally with tools designed for this purpose or inside a shear flow, the behavior of a dense packing of soft particles under stress remains largely unexplored. Unlike weakly deformable particle materials, the mechanical behavior of an assembly of soft particles is governed by both rearrangements and shape change of the particles, and also partially by their compressibility. Particle shape change modifies drastically the collective behavior of the particles, allowing the packing to reach high packing fraction beyond Random Close Packing (RCP) state.
A major obstacle to the numerical study of soft particle materials has been the lack of a proper discrete modeling for the soft particles. Discrete Elements Methods (DEM) (like Contact Dynamics (CD) [3,4]) often used for hard particles, are based on the approximation of small strains of contacts between particles. Thus, the only degrees of e-mail: nthai@dut.udn.vn e-mail: saeid.nezamabadi@umontpellier.fr e-mail: jean-yves.delenne@supagro.inra.fr e-mail: franck.radjai@umontpellier.fr freedom considered in these techniques are those of particles (displacements and rotations). For soft particles, in contrast, it is necessary to take into account the internal degrees of freedom of each individual particle as well as the contact interactions.
In this paper, we consider the Material Point Method (MPM) and CD method, both recently implemented and applied to small assemblies of particles [5]. We also consider a Bonded Particle Model (BPM) in which each particle is represented by an aggregate of primary rigid particles with appropriate long-range interactions such that the aggregate can undergo plastic deformations without rupture of the aggregate and dispersion of primary particles [6]. The major advantage of this latter approach is that it uses the framework of the DEM. However, as we shall see below, the resulting behavior is plastic. In the following, we first briefly introduce the applied numerical approaches. Then, we present the behavior of a single deformable particle. Finally, we study the compaction of an assembly of soft circular particles. We conclude with a brief summary and scope of this work.

Numerical procedures
There are two general possible representations of ultrasoft particles: 1) as an aggregate of rigid primary particles with attractive interactions (BPM) and 2) as a set of material points (MPM). In the MPM, the information carried by material points is projected onto a background mesh, where equations of motion are solved. The mesh solution is then used to update the material points. The important point for the treatment of soft particles with MPM is the necessity of coupling it with the CD approach for a rigorous treatment of frictional contacts. More details of this procedure can be found in [5].
The BPM approach has been mainly used for the simulation of breakable particles [7]. Its application to soft particle materials involves the introduction of a long-range attraction force between primary particles. Hence, the interaction force between primary particles is composed of two forces: i) a long-range center-to-center attraction force F a , as in colloidal systems, and ii) a short-range repulsion force F r [8]. In this manner, each BPM particle is an aggregate of primary rigid particles, which can move and rearrange according to the external forces acting on its boundaries by other aggregates.
To define the interaction force between primary particles, let us consider two rigid circular particles of the radii a 1 et a 2 separated by a distance δ. We use a Lennard-Jones attractive force: where a 0 = a 1 + a 2 , and F 0 is the maximum value of the attraction force at contact between particles (δ = 0) [8,9]. The choice of the exponent γ is not crucial for this study. It has been shown that this law with different values of this exponent may represent the potential of the mean force in complex materials such as cement and clay [9]. The repulsive part of the interaction between primary particles can be defined in different ways. But, as we use the CD method for the simulation of primary particles, we consider that the particles are perfectly rigid so that the repulsion force F r is activated only when the two particles touch (δ = 0). Therefore, the interaction force F = F a + F r as a function of δ consists of only attraction at arbitrary distance δ > 0 and takes a value ∈ [−F 0 , +∞] at contact (δ = 0). In the CD method, this force is calculated by using nonsmooth an implicite integration scheme [3,4,10]. To avoid strain shear localization and volume change within aggregates, the coefficient of friction between primary particles belonging to a particle is set to zero. Detailed descriptions of the BPM method can be found in [6].

Diametral compression of one particle
In this section, we study the accuracy and efficiency of the proposed models by considering the behavior of a single particle by using the MPM and BPM approaches. To do so, we compressed a soft particle of radius R = 5 mm between two rigid walls by fixing the bottom wall and moving downwards the top wall at a constant velocity of 0.2 m/s. The BPM particle is composed of 1750 rigid frictionless primary particles with a weak size polydispersity (diameters ∈ [0.16, 0.26] mm). Their interactions are governed by equation (1) with F 0 = 100 N and γ = 7. For an efficient implementation of this law, a cut-off distance δ = a 0 is introduced. Indeed, for γ = 7, the attraction force is negligibly small beyond δ a 0 , in which case only the first and second neighbors of each primary particle are involved in the interactions with other primary particle.
Two-dimensional MPM simulations in plane strain conditions were performed in the same conditions. The computation domain was meshed with four-node quadrangular elements, and an initial distribution of four material points per element was used. Young's modulus, Poisson's ratio and density of the particles were set to E = 10 MPa, ν = 0.45 and ρ = 990 kg/m 3 , respectively. We chose a high value of Poisson's ratio in order to get nearly constant volumes of the particles as in BPM particles. The spatial relative resolution is Δr R = 0.029, where Δr is the mean distance between material points.   Figure 1 shows the vertical stress σ = F L , where L is the horizontal diameter of the particle at any time step, as a function of the cumulative axial strain ε = ln(1 + d/R), with d being the displacement of the center of the particle. The MPM particle shows a linear elastic behavior at low strain as predicted by the Hertz analysis for a disk [5]: Deviation from linear behavior is observed for ε > 0.05. This range of the Hertz scaling (with an exponent 3/2 in 3D) is generally used in molecular dynamics simulations of granular materials composed of spherical particles by assuming small contact strains [11,12]. A transition to a linear mode with a lower slope is observed beyond ε = 0.15. The BPM particle is characterized by a linear response at small deformations of the order of 2% and a plastic behavior beyond. The plastic strain threshold remains constant and is almost equal to a characteristic stress defined by σ I = F 0 <a 0 > = 0.5 MPa, where < a 0 > 0.2 mm, is the mean distance between primary particles. Moreover, to demonstrate the plastic behavior of BPM particle, we applied an unloading at two strain levels (10% and 30%); see Fig. 1. The stress decreases linearly during unloading with a slope that appears to increase slightly with strain.

Uni-axial compression of a soft particle system
In this section, we use the BPM and MPM procedures to analyze the compaction of a packing of soft particles. It is composed of 300 soft circular particles confined inside a rectangular box of width L in which only the top wall is mobile and moves downwards at constant velocity of 5 m/s with a time step of δt = 0.1 μs. In the initial configuration, a small size polydispersity is introduced in order to avoid long-range ordering (diameters ∈ [1.4 , 2.4] mm). The gravitational acceleration and friction between the particles and the walls are set to be zero in order to avoid stress gradients. We consider two cases below: 1) without friction and 2) with a coefficient of friction μ = 0.5 between the particles. In the BPM simulations, each particle is composed of a minimum of 200 primary particles. The same material parameters are considered for the two methods as in previous section. Figure 2 shows the vertical stress σ, calculated from the forces acting on the bottom wall and normalized by the consolidation stress σ 0 , defined for the packing fraction Φ 0 , as a function of the cumulative vertical strain ε = ln(1 + Δh/h 0 ). The couple (σ 0 , Φ 0 ) defines a reference state. We chose Φ 0 = 0.8 which corresponds to the beginning of particle jamming. In the MPM simulations, the relationship between stress and the cumulative strain is linear whereas in the BPM simulations, a nonlinear increase of stress with ε is observed. Three deformation mechanisms control this development: (1) the elastic and/or plastic volume change of particles, (2) the rearrangement of particles and (3) particles shape change. The latter dominates the behavior of the system and allows the packing fraction to exceed considerably the RCP state. It is also worth noting that in the BPM simulations, the vertical stress grows at an increasing rate with ε as in the compaction of metallic powders [13]. Moreover, in both cases, the friction between particles has little effect on the deformation.  The microstructure of an assembly of soft particles can be characterized in terms of the coordination number Z. Above 20% of deformation, the particles are so much deformed that the contacts can no more be modeled as point contacts. In this regime, the particle shapes are well approximated by polygons, and Z then represents the number of sides of these polygons. Figure 3 shows the mean coordination number Z and the packing fraction Φ as a function of cumulative vertical strain ε defined from the beginning of the simulations. Φ is an almost linear function of ε. Indeed, as the volume of the particles does not vary, we have Φ = Φ 0 e ε Φ 0 (1 + ε) (Φ 0 = 0.7). In the case of MPM particles, the behavior is linear elastic, and therefore the remaining small pores for high values of packing fraction can only be filled at much higher stresses (requiring the radii of curvature to be increasingly smaller). However, in the BPM simulations, the pores can be filled more efficiently by plastic deformation of the particles for a much lower stress level. The porefilling process is similar in both approaches as the pore size between the BPM particles is greater than the size of the primary particles. The limit where the internal pore sizes of the particles are of the same order as the pores at the interface between two particles, corresponds to 40% of deformation and a packing fraction of 0.93 in the BPM simulations.
The mean coordination number Z increases with ε. In the BPM simulations, we observe no difference between frictional and frictionless particles, which is due to the absence of friction between the primary particles. Indeed, the friction force at the contact between two primary particles a t the boundary of two particles in contact, can not be mobilized because these primary particles can rotate freely without dissipation. In MPM simulations, we see slightly higher values of Z for frictionless particles while the packing fraction is the same in all cases with or without friction at all levels of deformation. The maximum value of Z is 5.5 in all simulations. This value reflects the low size polydispersity of particles.
The jamming transition occurs for ε 0.08 where Z 3 and Φ 0.75. From this state, the vertical stress begins to increase. However, the assembly remains fragile (rearrangements of particles) until a critical packing fraction of Φ 0 0.8 for a stress σ 0 , which was used earlier as the the consolidation stress. Figure 4 shows the evolution of Z as a function of packing fraction for both approaches and both cases of with and without friction. We see that by normalizing Z − Z 0 per Z 1 − Z 0 , where Z 0 is the coordination number corresponding to Φ 0.8 and Z 1 for Φ 1, all data points collapse on a function well-fitted by a power law. Beyond Φ 0 , Z increases as (Φ − Φ 0 ) 0.5 as in the case of emulsions and foams [1,14]; see Fig. 4. This dependence of the coordination number with the packing fraction is quite robust and almost independent of the interaction potential, polydispersity and friction [14].

Conclusion
In this paper, we analyzed the uniaxial compression of a soft-particle assembly beyond the RCP state using two numerical techniques: Material Point Method (MPM) and Bonded Particle Model (BPM). We characterized the macroscopic behavior, and simulations by two methods allowed us to highlight several significant differences between the assemblies of elastic and plastic particles. One interesting finding concerns the relationship between the packing fraction and the coordination number by a simple power law as also observed in emulsions.