Effect of packing characteristics on the discrete element simulation of elasticity and buckling

Abstract Discrete Element Method (DEM) simulations are used to model the elastic properties of a continuous material. The preparation route of the particle packing is shown to have a significant effect on the macroscopic properties. We propose simple relations, generalized from the mean field solution, that are able to fit DEM results. These relations introduce only basic microstructural features such as the coordination number and the packing density. When the tangential to normal stiffness ratio increases above unity, the material becomes potentially auxetic. Buckling is also explored with DEM, and results on cylindrical bars are compared to the classic Euler results for critical stress.


Introduction
Originally designed for granular materials [6], the Discrete Element Method (DEM) has been employed in the recent years to model continua. DEM, when used to model continuous materials, bears some resemblance with lattice methods, which model a solid as a set of nodes (with no volume or mass) connected by beam elements through which both attractive and repulsive forces and torques can be transmitted [27]. This type of method has proven very efficient for simulating fracture, although crack closure may be an issue since nodes carry no volume. A comparable route has been adopted by models to simulate the behavior of cohesive materials with bonded particles under a static or kinematic assumption [9,12,23,15]. Peridynamics is also a similar method proposing a theory of continuum mechanics with a nonlocal model that can simulate fractures and discontinuities [29]. However, the simulation of Poissons ratios other than 1/4 is not straightforward with Peridynamics [18].
In DEM , the solid is represented by bonded particles that rearrange under applied deformation according to Newton's second law. In general, these particles are spherical to ease distance calculation and to benefit from the existing literature on spherical contact forces [35,31,2,13]. In specific applications, it is more interesting to use other shapes, like tubes for example [16]. Nevertheless, spheres have the advantage of ensuring isotropy of the contact network when packed randomly [32,2]. Also, they may be used to obtain a mesh for complex geometries, for example when starting from three-dimensional images of real objects [10].
Most of the above-cited examples of DEM applications for continuum are actually interested in breaking the initial continuum material into discrete entities to model damage and fracture. This is because DEM has shown a great ability to model fracture in continuous material. Indeed, the inher-ent discrete nature of DEM provides an effective formulation with in-built natural discontinuities. DEM has the additional advantage, as compared to the lattice or to the finite element methods, of dealing naturally with post-fracture stages where crack closure may play a role.
In this paper, we aim to provide relations that link the bulk elastic properties of a solid simulated by a packing of bonded particles to the microscopic and macroscopic characteristics of the packing. These relations introduce macroscopic basic parameters that define the structure of the packing, such as the coordination number (the average number of bonds per particle) and the packing density. The microscopic properties of the bonds, i.e. their normal stiffness and their tangential to normal stiffness ratio, are also included.
These relations are generalization of the classical mean field solutions applied to a random packing of spherical particles. They should provide simple and useful guidelines for the generation of particle packings. We explore also the capability for DEM to reproduce auxetic (negative Poisson's ratio) materials. Finally, the possibility to use DEM to model buckling in quasistatic conditions is demonstrated. DEM-type simulations were carried out on half-spherical shell under dynamic conditions, but were not compared quantitatively to analytical results [15]. Here we attempt a direct comparison with Euler's formulation under quasi-static conditions on a cylindrical bar.

Model description
Spherical particles with radius R are connected to each other through solid bonds that transmit normal and tangential forces and resisting moments ( Fig. 1). An equivalent radius at the bond between two spheres of radii R p and R q is defined as: The normal and tangential forces acting on the bond are given by simple linear elastic laws: where δ N and δ T are the normal and tangential relative displacements between the two particle centers and n and t are the unit vectors normal and parallel to the contact plane, respectively. The normal and tangential stiffnesses K N and K T are size dependent. Thus, we prefer using material parameters Σ N = K N 2R * and Σ T = K T 2R * , which have the unit of stress and allow macroscopic elastic properties to be independent of the sphere size. Noting: we restrict our investigation to the range 0 < α ≤ 10. The standard case is given by 0 < α < 1, while we will show that α ≥ 1 potentially leads to auxetic properties (negative Poisson's ratio). Forces are taken positive in tension while the tangential force opposes the accumulated tangential displacement.
Damping forces are included for computational convenience. They oppose the normal relative displacement velocity at the bond: where η n is the damping coefficient, which is chosen as a fixed fraction ξ of the critical damping: A value of 0.05 for ξ ensures underdamped behavior together with negligible effect on macroscopic stresses in quasi-static conditions. Similarly, damping forces oppose the tangential relative displacement velocity with ξ = 0.025.
Bonds transmit resisting moments, M N and M T , in the normal and tangential directions with a formulation originally proposed by Potyondy and Cundall [24]: where θ N and θ T are the accumulated relative rotations in the normal and tangential directions in the contact framework ( Fig. 1). Note that in a packing with multimodal size particles, the stiffness and resisting moments from one bond to another are distributed due to the use of the equivalent radius R * in Eqs. (2) to (8).
The macroscopic stress tensor in the packing is calculated from Love's formulation [36,5]: where the summation is made on all contacts and where V is the sample volume, F i is the i th component of the total contact force (with normal and tangential terms), and l pq,j is the j th component of the l pq vector connecting the centers of two particles p and q (Fig. 1). The simulations described hereafter have been carried out on our in-house code dp3D.

Sample preparation
Three dimensional cubic samples were prepared by packing 5000 spherical particles within a box bounded by periodic conditions. Under periodic boundary conditions, when a particle protrudes outside the periodic cell through a given face, it interacts with the particles on the opposite face. The first preparation step consists of creating a gas of particles. This is attained by locating particles randomly into the volume with the only constraint that they do not contact each other. Using a packing of nearly monomodal particles (R ± 5%) with an average radiusR, a relative density of D 0 = 0.3 is initially obtained. Ten such initial samples were generated using different random seeds to compute standard deviations for each set. As described in the next section, the standard deviation on the elastic properties was small enough to consider that ten samples for each set of conditions are sufficient.
Following this first step, a relative density of 0.65 is sought for, using three 6 different preparation routes. The lengths of the cubic simulation box, after densification to 0.65, are L x = L y = L z = 32R. We found that, under periodic conditions, this size together with 5000 particles is enough to obtain statistically meaningful results.
The first route simply prescribes an affine displacement to all particles to impose a macroscopic isostatic densification until the preset relative density is attained (D = 0.65). This densification technique leads to significant relative interpenetration between particles. We define this packing preparation route as unjammed.
In the second method, starting from the initial relative density of 0.3, the initial gas of particles is isostatically densified by moving inwards the periodic bounds. No normal tensile force nor tangential friction nor resisting moment is transmitted (i.e. only normal repulsive forces are used) at this stage. Weak jamming is obtained by imposing a control pressure P c to the packing (using Eq. (9)) with a proportional controller that dictates the hydrostatic strain-rate to the periodic bounds as described elsewhere [20].
The control pressure is several orders of magnitude smaller than the stiffness (P c /Σ N = 1.0E −07 ). The stiffness of particles is thus large enough to ensure that relative interpenetration between particles is negligible at this stage.
Force equilibrium is imposed in a quasi-static configuration by resolving explicitly the second law of Newton as classically done in the DEM and by ensuring that the total force applied on particles is smaller than where N max is the maximum contact force on the particle. This scheme is imposed until a preset relative density of D 0 = 0.5 is attained. This relative density is much below the Random Close Pack relative density (0.63-0.64).
Such a low value is still interesting because it ensures a fast sample preparation since particles rearrange easily. This packing is further densified, using affine displacement, up to D = 0.65. We define this packing preparation as weakly jammed.
The third method aims at obtaining a fully jammed packing by continuing the densification scheme described above at the control pressure P c until the densification rate attains very low values [20]. At this point, we consider that further densification is prevented by the full jamming of the particle assembly. A few particles (rattlers) with less than 4 contacts are still observed (< 2%). All other particles cannot be displaced without generating an increase in contact forces with neighboring particles, thus qualifying this packing as jammed [33]. The relative density of this jammed packing is D = 0.636 ± 0.0004 and the coordination number is Z = 6.03 ± 0.011. As with the two other methods, this packing is further slightly densified with affine displacement of particles to attain D = 0.65.
The final affine displacement imposed to particles to attain the final packing density 0.65 allows all three packings to have the exact same density. It modifies only very slightly the microstructure of the jammed packing. We have studied the fabric tensor [25] of the packings generated through this methodology and found that they did not exhibit any sign of anisotropy.
We have thus three preparation methods: unjammed, weakly jammed, and fully jammed, which lead, after the last affine densification to three packings that exhibit the same number of particles and the same relative density  interparticle distances as compared to the jammed packing, it still presents the same trends as the jammed packing. In short, the unjammed packing is characterized by a small number of contacts that are very broadly distributed in size, while the jammed packings have a large number of smaller contacts.
From these three initial packings, bonds are installed between two particles p and q with radii R p and R q when the following criterion is met: where |l p,q | is the center to center distance between particles p and q (Fig. 1 (Fig. 4). This feature provides the possibility to adjust the degree of interlocking of the particles forming the numerical medium as proposed by [13] or [28]. The average numbers of bonds per particle, denoted as Z b , are listed in Table 1. Z b is strictly equal to Z for κ = 1.000 and increases with increasing interaction range κ.

Elastic properties
The macroscopic elastic properties of the cubic numerical samples, of which the preparation has been described in the preceding section, have been obtained by testing them uniaxially under fully periodic conditions. Eq.
(.1) in the appendix is used to determine the macroscopic stress response to the imposed strain and the two macroscopic elastic constants, Young's modulus E and Poisson's ratio ν. The normal to tangential stiffness ratio (α = Σ T /Σ N ) has been varied from 0.005 to 10 to test the range of elastic constants attainable with potentially auxetic materials (ν < 0 for α > 1).
The two parameters discussed in the preceding section (D and Z b ), which describe the average microstructure of the packing of particles, can be advan-tageously used to propose a simple normalization of the Young's modulus: The rationale for Eq. (11) is detailed in the appendix. It derives from the mean field assumption, which has been used to model particulate packings [5,17,19] The effect of the ratio between the tangential and normal stiffnesses (α = Σ T /Σ N ) may be captured for the normalized Young's modulus by the following simple equation:Ẽ = κ n 2π where n is a rational number, and a i are fitted parameters. Similarly, for the Poisson's ratio: where b i are fitting parameters. The simple form of Eqs. (12) and (13) originate from a generalization of the equations that describe the elastic behaviour of an aggregate of bonded particles under the approximations of a mean field or of a static solution [17], as detailed in the Appendix.
We first specialize this study to α values in between 0 and 1, which is the For all packings, the normalized Young's modulus increases and the Poisson's ratio decreases as the tangential stiffness increases as predicted qualitatively by the mean field model (see appendix) and by the static hypothesis model from Liao et al. [17]. Still, the normalized Young's moduli of unjammed and weakly jammed prepared packings are significantly smaller than those predicted by the mean field and static models. The fully jammed packing is much closer to the mean field model. In particular, the slope is very similar in a large range of α values. This is because the mean field model has been derived under the assumption of point contacts between particles that are distributed with uniform probability over their surface [34]. The unjammed and weakly jammed packings hardly satisfy these two assumptions at all. Discrete element simulations consistently lead to a softer response as compared to the mean field model. This is expected since quasi-static DEM ensures force equilibrium, in contrast with the mean field method, which does not. It is expected that the mean field method, because of its imposed assumption on uniform strain, leads to a stiffer response than the true behavior. Table 2 lists the values of the fitting parameters a i and b i . These equations allow the macroscopic elastic properties of a packing of particles to be predicted with simple forms using important microstructural parameters Table 2: Parameters a i , n and b i of Eqs. (12) and (13)  (bond coordination number, interaction range of these bonds, and density).
Of particular interest is the value of the a 1 parameter, which dictates the behavior of the packing when the bond tangential stiffness vanishes. As discussed by Liao et al. [17], if a 1 tends to zero (as in the static hypothesis or for unjammed packings), a packing with smooth particles (zero bond tangential stiffness) has no shear modulus but retains bulk modulus if its Poisson's ratio tends to 0.5. In that case, the packing resists isotropic bulk stresses but no deviatoric stress. It thus behaves as a fluid. Reversely, if a 1 is non-zero, packings do retain shear modulus and behave solid-like.
The same methodology is used to obtain the behaviour of packings in the range 1 ≤ α ≤ 10. These packings potentially lead to auxetic materials (negative Poisson's ratio) as shown previously in 2 dimensions from analytical considerations similar to those developed in the appendix [3] or by discrete element simulations [8]. The same tendencies are observed in this range of α values than in the standard 0 < α ≤ 1 range (Fig. 6). The normalized Young's modulus is in relatively good accordance with the mean field model

Application example: buckling of a cylindrical bar
Using the sample preparation method described in section 3, we have built cylindrical samples with length l and radius r. However, instead of periodic conditions, a rigid cylinder and two planes were used as boundary conditions. Once jammed, the cylinder and the planes were removed and bonds were installed. Fixed-fixed conditions are applied as boundary conditions by blocking translation and rotation of the particles in contact with the two planes.
We investigated l/r ratios in between 10 and 120 with the same number of particles per unit length. Thus, for l/r = 10, 10 000 particles were used while 100 000 particles were used for l/r = 100. Three samples were generated for a given l/r value to estimate dispersion. This was simply obtained by using different random seeds for the initial gas of particles. We have limited the sample number due to the CPU intensiveness of the simulations involving buckling. Simulations with 100 000 particles and large strains (to attain buckling), although heavily parallelized, are too long to allow for good statistics with a larger number of samples.
As in the preceding section, the same preparation routes were tested (unjammed, weakly jammed and jammed). A value α = 0.2 was chosen for the tangential to normal stiffness and the interaction range κ was set to unity. These values lead to a Poisson's ratio of the order of 0.2. No initial imperfection in the geometry, in the material properties or in the loading direction were introduced. This in contrast with the method classically used in Finite Element Method [4,26,22], where a small imperfection in the element nodal points or in the load application is imposed to ensure the initiation of buckling. However, it should be clear that the randomness of the packing provides some imperfection at the microscale.
Euler's critical stress is given by [7]: where k = 0.5 for the fixed-fixed conditions used here. Alternative solutions have been proposed. For example, Mazzilli [21] has shown that the analytical compares better to the critical load computed by numerical integration of the exact equation of equilibrium than Euler's formulation. However, the difference is only noticeable for bars with small slenderness. Fig. 7 shows typical stress-strain curves. Axial stresses and strains are normalized by σ Euler and ǫ Euler respectively (ǫ Euler = σ Euler /E). Fig. 7 indicates that the critical strain is larger than the critical strain given by Euler's formulation. Also, the post-buckling behaviour depends on the l/r ratio with a more abrupt decrease of stress as the l/r ratio increases. As 20 shown in Fig. 7, buckling is characterized by a clear maximum in the stressstrain curve. However, a non-linear region exists prior to this maximum.
Still, we compare the critical Euler buckling stress to this maximum stress. Indeed, dynamics lead to buckling stresses that exceed Euler critical stress [14]. We have investigated this issue by calculating the normalized kinetic energy per particle:Ẽ kin = E kin n max (N R) (16) where E kin is the kinetic energy of the system, n the number of particles, N and R the normal force (Eq. (2)) and the particle radius, respectively.Ẽ kin should stay reasonably below 10 −08 to ensure quasi-static conditions [1]. Fig.   9 indicates that this condition is satisfied before buckling for the two lowest

Concluding remarks
The elastic properties of numerical samples made of particles linked by elastic bonds have been investigated. The preparation route used for the samples clearly has an effect on these properties. We have shown that when the ratio of tangential to normal stiffness (α) is smaller than unity, the Young's modulus and the Poisson's ratio are well described by simple equations (Eqs. (12) and (13)), which introduce microstructural properties of the packing of spheres that has served to generate the numerical sample. These microstructural properties are the average coordination number and the packing density.
The proposed equations are generalizations of the mean field and static solutions. They give approximate values of the elastic properties that should be expected for a given numerical sample, knowing its packing history.
Our results also indicate that, with the standard bond model used here, it is difficult to attain Poisson's ratio above 0.3. This is inherent to the fact that the Poisson effect is not introduced at the micro level in the bonds, since discrete element simulations are based only on pair interactions. Taking into account a more accurate strain field would require a more complex coupling between bonds. This has been proposed analytically [12] and developed numerically with an explicit calculation of bond interaction on the elastic field [11].
The three preparation routes investigated here involve very different effort in terms of CPU time to create the packing that meshes the numerical sample.
The unjammed route is a very rapid method since it only requires the random deposition of spheres into the volume. The fully jammed preparation route is much more CPU demanding if a density close to the Random Close Pack (RCP) is sought for. The weakly jammed route is an interesting compromise.
It allows for a wider span of attainable Poisson's ratio (also in the auxetic domain) than the unjammed route but is much less CPU demanding since densification rates in the low density domain are much faster than close to the RCP.
For bonds with tangential to normal stiffness greater than unity, we have shown that hypothetical auxetic materials built from discs [3,8], may also be generated in 3 dimensions. The requirement that bonds be softer in the normal direction than in the tangential direction is difficult to satisfy with standard materials. However, some interesting designs have been proposed in the literature for granular-like materials [8] .
As an application example, we have shown that the simple bond model proposed here, with no built-in buckling mechanism, is capable of simulating buckling. A buckling application example of DEM-type simulations was proposed on half-spherical shell under dynamic conditions [15]. Here, we have attempted to ensure quasi-static conditions to enable quantitative comparison with Euler's formulation. It should be clear that the DEM is not the most appropriate method for simulating buckling. The finite element method, supported by decades of research on the subject, is a far more effective method. However, when post-buckling is accompanied by fracture, the DEM may be an interesting option. This is because, the topological modifications (branching, bifurcation, and new surface generation) that come with fracture are easier to capture with DEM than with the finite element method.
Appendix. Elastic behaviour of an aggregate of bonded particles under the approximation of a mean field solution The spherical particle kinematics may be simply treated under the approximation of a mean field solution (Voigt hypothesis) for a monomodal packing. In that case, the motion of each particle center relative to a reference particle is dictated by the macroscopic deformation gradient and force equilibrium is not enforced. We start from the formulation derived by [17], which considers the case of bonds transmitting both normal and tangential forces (but no resisting moment). The model is similar to those proposed in [5,34,30,19]. Starting from Eq. (9) under the assumption of a random packing of particles of uniform radius R, it leads to the stress-strain relationship under the mean field assumption (Eq. (38) in [17]): The Poisson's ratio is then readily given by while the normalized Young's modulus simplifies tõ An alternative expression has been derived by [17], based on the static hypothesis (or best-fit model). It leads to: Thus, for the static hypothesis the Poisson's ratio is