Bulk modulus of soft particle assemblies under compression

Using a numerical approach based on the coupling of the discrete and finite element methods, we explore the variation of the bulk modulus K of soft particle assemblies undergoing isotropic compression. As the assemblies densify under pressure-controlled boundary conditions, we show that the non-linearities of K rapidly deviate from predictions standing on a small-strain framework or the Equivalent Medium Theory (EMT). Using the granular stress tensor and extracting the bulk properties of single representative grains under compression, we propose a model to predict the evolution of K as a function of the sample's solid fraction and a reference state as the applied pressure P tends to zero. The model closely reproduces the trends observed in our numerical experiments confirming the behavior scalability of soft particle assemblies from the individual particle scale. Finally, we present the effect of the interparticle friction on K's evolution and how our model easily adapts to such a mechanical constraint.


Introduction
Soft granular particulate materials such as powders, gels, bubbles, rubber chunks, and even cells are challenging materials to characterize and model due to the large deformation they can undergo. Besides, their discrete nature calls for adequate interaction laws for multi-contact systems. Amongst the different numerical approaches to model highly deformable particle assemblies, the discreteelement approach is one of the most frequently used. By means of smoothed interaction laws between bodies, simulations may reproduce some elasticity due to a virtual contact deflection, although the bodies themselves do not undergo strains. Using this approach, the elastic properties of particle assemblies have been studied, however restricting considerations of the small-strain domain of deformations [1][2][3][4]. More recently, the development of more advanced methods coupling discrete and finite elements [5][6][7][8][9][10] or meshless methods [11,12] have permitted to explore the compression behavior of soft granular media beyond jamming. Nonetheless, the study of the evolution of elastic properties of particulate assemblies undergoing large deformation is still challenging to characterize.
In this paper, we simulate assemblies of 2D circular particles under isotropic compression using the contact dynamics method and the finite-element method to account for the large deformation of meshed bodies. Systematically increasing the pressure on assemblies of disks, we have access to the macroscopic stress-strain relation letting us deduce the evolution of the bulk modulus during the compaction. Then, using the granular stress tensor and * e-mail: david.cantor@polymtl.ca the behavior of individual representative particles, we propose an analytical equation for the bulk modulus evolution fitting very well our numerical experiments.
In Sec. 2, we present the details of the numerical approach and the isotropic test procedures. Then, in Sec. 3, we show the stress-strain relation measured at the sample scale and the bulk modulus evolution. In order to look for the origins of the macroscopic bulk modulus, we explore in Sec. 4 the behavior of individual particles under compression. In Sec. 5, we introduce an analytical approach based on the stress tensor decomposition from microstructural parameters letting us deduce an equation for the macroscopic bulk behavior. Finally, in Sec. 6, we present the effect of the interparticle coefficient of friction on the evolution of the bulk modulus and how our analytical model easily adapts to such mechanical constraint. Section 7 concludes this work with a summary and perspectives.

Numerical procedures
To simulate assemblies of soft particles, we used the coupling of the discrete method known as contact dynamics (CD) [13,14] and classic finite elements in the framework coined as non-smooth contact dynamics (NSCD) by M. Jean [15]. With this method, we were capable of building assemblies of circular meshed bodies that interact using unilateral contacts and dry friction.
We built samples composed of N p = 1500 circular bodies slightly disperse in size (the ratio between the maximal over the minimal particle diameter is 1.5), meshed using 92 triangular finite elements, and deposited within square boxes by means of an algorithm based on simple geometrical properties. For the finite elements, we used the non-Hookean hyper-elastic material model [16], setting incompressible bulk behavior, plane strain conditions, and an elastic modulus E. We previously tested the mesh resolution concluding that the number of finite elements per grain we ended up using does not compromise the results' quality. Then, we set pressure controlled conditions on the four rigid walls of the boxes in a series of steps varying the relative pressure P/E in the range [1 × 10 −4 , 5 × 10 −1 ]. The interparticle coefficient of friction µ was set to zero and gravity was neglected to avoid pressure gradients. Also, note that the loading was performed using slow gradual steps to avoid dynamic effects and promote a rapid dissipation of elastic waves. Figure 1 shows the sample's configuration at the beginning of the loading and its deformed state at the end of the tests.

Macroscopic bulk modulus
In our simulations, we can easily track the deformations of the servo-controlled boundary walls. For convenience, we relate the volumetric strain of the sample ε v to the solid fraction φ, as ε v = − ln(φ 0 /φ), with φ 0 being the solid fraction as the ratio P/E → 0. This reference state can be understood either as the solid fraction the assembly presents with perfectly rigid bodies or with vanishing external pressure. The macroscopic bulk modulus K can be then computed as a function of the solid fraction as K = (dP/dφ)(dφ/dε v ). Figure 2 presents the evolution of K as a function of φ. Different strategies have been used to predict the nonlinear evolution of the bulk modulus, being most of then introduced in the small-strain framework or the Equivalent Medium Theory (EMT) [17][18][19][20]. That approach adopts an analogue model considering a set of springs joining the center of mass of bodies in contact and whose deformation represents the relative approaching of their centers as they deform. Using this type of approach, it is then possible to integrate the set of springs' deformation and deduce a stress-strain relation for the whole system; thus, a bulk equation can also be deduced. Nonetheless, the definition of the spring behavior is, in that vein, at the origin of the macroscopic strains and bulk evolution. Similarly, smooth discrete-element methods can consider the forceoverlapping relations to be the spring characteristics for an equivalent medium analogy. If we use a similar strategy, we can first consider to be the average spring length and ε = ln( /d) to be the average strain, with d being the average particle diameter (our numerical experiments also show that ε v = 4ε ). Second, the homogeneous field of springs allows us to suppose that the contact level pressure is P = Eε . The granular stress tensor can be written as σ i j = n c f i j , with the contact number density being n c = N c /V, N c is the number of contacts, V is the sample volume, and . . . the average of contact forces f and branch vectors (i.e., the inter-center vector of particles in contact). We can then deduce a micromechanical definition of pressure as with the coordination number Z = 2N c /N p and σ = f /d 2 . Note that to deduce the expression above, we consider P = (σ 1 + σ 2 )/2, where σ 1 and σ 2 are the principal stresses of tensor σ i j . If we consider that the coordination number evolves as a power-law of the solid fraction (as consistently shown in previous studies as [5,12,21], and also verified in our simulations), in the form (Z −Z 0 ) = k(φ − φ 0 ) α , with α = 0.5, Z 0 the coordination number in the reference state, and k a proportionality parameter easily deduced knowing that, when φ tends to unity, the particle structure tends to a hexagonal-like arrangement and Z → 6. So, k 5. Finally, considering that σ = P , we can deduce a microscopic definition of the bulk modulus upon the derivative of Eq. (1) as This equation is displayed in Fig. 2 with a solid blue line. We can observe that the predictions with such an expression can be considered relatively good for the first part of the compression, and up to φ 0.85. Beyond that value, the bulk modulus starts to increase more rapidly and is expected to diverge as the solid fraction tends to unity (i.e., the sample ends up behaving as an incompressible solid). The equivalent medium composed of springs cannot capture the diverging behavior of the bulk. It is necessary then to track the evolution of K for the multi-particle system using an alternative approach.

Single-particle scale
Let us consider the system composed of a single soft particle inside a square box and following the preset boundary conditions as undertaken with the assembly. Figure 3 presents the evolution of the corresponding solid fraction (φ p ) as a function of the applied pressure P p . Note that the volumetric strain for this case is akin to the multi-particle system but denoted ε v,p . The compression behavior of the single-particle test can be considered analogous to the collapse of a cavity inside a circular body employing a homogeneous external pressure [22]. For that case, elastic solutions allow us to write the relation P p − φ p as with φ p,max the solid fraction as P p → ∞, and φ p,0 = π/4 the solid fraction at the reference state (i.e., as P p → 0). Finally, parameter b is found out to be 0.12 after fitting Eq. (3) to our data. We can observe that this expression fits very well the behavior of the single-particle compression as it is plotted with a red dashed line in the same Fig. 3. We can then easily deduce an analytical bulk modulus equation by using the derivative of the last expression, yielding to K p /E = (bφ p )/(φ p,max − φ p ). In sum, at singleparticle scale, we can comprehensibly identify analytical equations for both the compression and bulk evolution.

Conciliating scales
Let us consider again Eq. (2) on the macroscopic bulk modulus. In the limit of φ → φ 0 , and in agreement with many previous works [17][18][19][20], that equation shows that K ∝ ZφE. This means that, for small deformations, the bulk modulus scales via the Young modulus of a single particle and the structural parameters Z and φ of an assembly of particles. Naively, we may wonder whether the single-particle bulk modulus K p is more representative than the Young modulus E for such scaling. In other words, is there a more general scaling between the bulk modulus of the multi-particle and the single-particle system? In fact, by comparing these two systems at equivalent deformation (i.e., ε v = ε v,p ), our numerical simulations reveal that Equation (4) points out that the single-particle configuration can indeed be considered the smallest representative scale in our multi-particle system. Now, by replacing K p by its analytical form identified before and by mapping φ with φ p (since ε v = ε v,p ), we get a general micro-mechanically based constitutive equation for the bulk modulus evolution beyond the jamming as Equation (5) is shown in Fig. 2 showing an excellent agreement with the measured K in our tests. In this case, the equation follows the nonlinear increase of K for large deformations and its divergence.

Effect of the coefficient of friction
The inter-particle coefficient of friction is a mechanical constraint rapidly limiting the particle reorganization and, thus, their capacity of filling voids. We reproduce a series of tests with the same samples as for the frictionless case, but setting this time the interparticle coefficient of friction to µ p = {0.2, 0.4, 0.6, 0.8}. We tracked the evolution of pressure and deformation, and finally, we could compute the evolution of K for these cases. Figure 4 gathers the results for the evolution of K as a function of the assembly's density for the different values of µ p . We can observe that for small increments of solid fraction after jamming, the evolution of K is indistinguishable for the different coefficients of friction. Nonethe-less, for higher values of solid fraction (i.e., larger deformations), the data gathered from the simulations slightly shifts towards the left as µ p increases. This means that the additional constraint added by friction limits the deformability of the assembly. Many authors have shown that friction modifies the jammed state of rigid grains assemblies and the maximal solid fraction a soft assembly of grains can reach [23,24]. Consequently, the reference state, as P → 0, and the limit state, as P → ∞, are modified with respect to the frictionless case. These facts indeed affect the proposed model in Eq. (5). While for the frictionless case Z 0 3.9, φ 0 0.805, and our fitting procedure resulted in φ max 0.998, the case with µ = 0.8 shows that Z 0 3.6, φ 0 0.795, and fitting of the data ends up estimating φ max 0.978. Using these parameters, we plot in the same figure the predictions of bulk evolution using our model with dashed lines. We can observe that the fitting is still excellent, and we can follow the macroscopic evolution of the bulk modulus up to densities very close to unity.

Conclusions
We studied the bulk modulus behavior of two-dimensional particle grain assemblies undergoing isotropic compression using a coupled discrete-finite element simulation platform. By finely describing the assemblies' compaction behavior, we presented the bulk modulus evolution as a function of the solid fraction.
We first introduced a model approach for the bulk modulus evolution standing on the ideas of the equivalent medium approach and considering the soft particle assembly as a series of interconnected springs governed by the particles elastic modulus E. We observed that such an approach leads to a formulation of the bulk modulus evolution that is acceptable only for relatively small strains and does not reproduce the divergence of K as the density approaches unity. To overcome this challenge, we first studied the compression of a single particle under isotropic compression in a box. We observed that the compression behavior and bulk evolution for such a case could be analytically obtained based on elastic solutions of a circular geometry under compression. Then, standing on the definition of the granular stress tensor and supposing that the most elemental scale of the assemblies is the singleparticle, we conciliated the single and multi-particle scales and deduced a model equation for K. Our model turned out to reproduce the evolution of the bulk modulus very well, highlighting its strong non-linearities and asymptotes. We finally showed that it is straightforward to introduce friction between particles into our model and then consider the impact this mechanical constraint has on limit configurations the system can reach.
The reader can also deduce that our modeling approach lets us describe the compaction evolution the writing of a relationship P vs. φ. This constitutive description can be found in Ref. [25] by the same authors of this work, or in Ref. [26] taking into account mixtures of rigid and de-formable particles. Nonetheless, to enrich the constitutive modeling of granular materials and structures upon discrete systems, it is still necessary to explore the shear behavior of soft particle media. It is also important to extend this work to varied shape and size particles, as well as to the 3D case.