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

Rishi Kumar, Sarshad Rommel, David Jauffrès, Pierre Lhuissier, Christophe

L. Martin

► To cite this version:

Rishi Kumar, Sarshad Rommel, David Jauffrès, Pierre Lhuissier, Christophe L. Martin. Effect of packing characteristics on the discrete element simulation of elasticity and buckling. International Journal of Mechanical Sciences, 2016, 110, pp.14-21. 10.1016/j.ijmecsci.2016.02.009 . hal-01367851

HAL Id: hal-01367851

https://hal.science/hal-01367851

Submitted on 17 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Effect of packing characteristics on the discrete element simulation of elasticity and buckling

Rishi Kumar, Sarshad Rommel, David Jauffrès, Pierre Lhuissier, Christophe L. Martin*

Univ. Grenoble Alpes, CNRS, SIMAP, F-38000 Grenoble, France

## 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.

**Keywords:** Discrete Element Method; buckling; elasticity

## 1. Introduction

Originally designed for granular materials [6], the Discrete Element Method (DEM) has been employed in the recent years to model continua. DEM,

<footnote>*Corresponding author. Email: christophe.martin@grenoble-inp.fr</footnote>
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 non-local 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 proper-
ties 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 nor-
mal stiffness and their tangential to normal stiffness ratio, are also included.
These relations are generalization of the classical mean field solutions ap-
plied 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 quasi-
static conditions is demonstrated. DEM-type simulations were carried out on
half-spherical shell under dynamic conditions, but were not compared quan-
titatively to analytical results [15]. Here we attempt a direct comparison
with Euler's formulation under quasi-static conditions on a cylindrical bar.

## 2. 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$
<figure><img src="image_2.png" /><figcaption>Figure 1: Sketch of the bond between two spherical particles p and q transmitting normal and tangential forces and resisting moments.</figcaption></figure>

and $R_q$ is defined as:

$$R^* = \frac{R_p R_q}{R_p + R_q} \qquad (1)$$

The normal and tangential forces acting on the bond are given by simple linear elastic laws:

$$\mathbf{N} = -K_N \delta_N \mathbf{n} = -2 \Sigma_N R^* \delta_N \mathbf{n} \qquad (2)$$

$$\mathbf{T} = -K_T \delta_T \mathbf{t} = -2 \Sigma_T R^* \delta_T \mathbf{t} \qquad (3)$$

where $\delta_N$ and $\delta_T$ are the normal and tangential relative displacements between the two particle centers and $\mathbf{n}$ and $\mathbf{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 $\Sigma_N = \frac{K_N}{2R^*}$ and $\Sigma_T = \frac{K_T}{2R^*}$, which have the unit of stress and allow
macroscopic elastic properties to be independent of the sphere size. Noting:

$$ \alpha = \frac{\Sigma_T}{\Sigma_N} \qquad (4) $$

we restrict our investigation to the range $0 < \alpha \le 10$. The standard case is given by $0 < \alpha < 1$, while we will show that $\alpha \ge 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:

$$ \mathbf{N}_v = -\eta_n \frac{d\delta_N}{dt} \mathbf{n} \qquad (5) $$

where $\eta_n$ is the damping coefficient, which is chosen as a fixed fraction $\xi$ of the critical damping:

$$ \eta_n = 2\xi\sqrt{2m\Sigma_N R^*} \qquad (6) $$

A value of 0.05 for $\xi$ ensures underdamped behavior together with negligible effect on macroscopic stresses in quasi-static conditions. Similarly, damping forces oppose the tangential relative displacement velocity with $\xi = 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]:

$$ M_N = -4\Sigma_T (R^*)^3 \theta_N \qquad (7) $$

$$ M_T = -2\Sigma_N (R^*)^3 \theta_T \qquad (8) $$

where $\theta_N$ and $\theta_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]:

$$ \Sigma_{ij} = \frac{1}{V} \sum_{\text{contacts}} F_i l_{pq,j}, \qquad (9) $$

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.

### 3. 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 ($\bar{R} \pm 5\%$) with an average radius $\bar{R}$, 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
different preparation routes. The lengths of the cubic simulation box, after densification to 0.65, are $L_x = L_y = L_z = 32\bar{R}$. 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/\Sigma_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 $10^{-3}N_{max}$, 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 \pm 0.0004$ and the coordination number is $Z = 6.03 \pm 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 ($D = 0.65$), but their microstructures differ markedly as demonstrated by the increasing coordination number with jamming extent: $Z = 4.52 \pm 0.02, 5.39 \pm$
<figure><img src="image_3.png" /><figcaption>Figure 2: Typical packings for the three preparation routes investigated, starting from the same gas of particles. Packing density is 0.65 after affine densification. Colors indicate the number of contacts per particle (coordination number).</figcaption></figure>

0.02, 6.56 ± 0.02, respectively. Fig. 2 shows the typical packings associated with each preparation route. It is difficult to discern differences between the three packings. However, the coordination number of each particle is a clear indication of the marked characteristics of the jammed microstructure as compared to the unjammed and weakly jammed.

The differences between the preparation routes is also clear when examining Fig. 3, which plots the distribution of interparticle distances between contacting particles. Fig. 3 shows that the unjammed packing exhibit a very large distribution, entirely resulting from the affine densification process imposed to the initial gas of particles. On the opposite, the jammed packing presents a very narrow range of interparticle distance and these distances are close to the diameter of particles (small interpenetration of particles). The weakly jammed packing shares important qualitative features with the jammed packing. Although it exhibits a broader distribution and smaller
<figure><img src="image_4.png" /><figcaption>Figure 3: Distributions of interparticle distance $l_{pq}$ normalized by the average diameter of particles, for the three preparation routes investigated. Packing density is 0.65 after affine densification. Average coordination are indicated for each preparation route.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Effect of the interaction range on the bonding stage between particles. (a) $\kappa = 1.000$ (b) $\kappa > 1$.</figcaption></figure>

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:

$$ \kappa (R_p + R_q) \geq |l_{p,q}| \qquad (10) $$

where $|l_{p,q}|$ is the center to center distance between particles $p$ and $q$ (Fig. 1), and $\kappa$ is the interaction range ($\kappa \geq 1.0$). Four values for $\kappa$ were tested: $\kappa = 1.000, 1.025, 1.050, 1.075$. $\kappa = 1.000$ leads to the classic granular description of the packing where only particles in contact are bonded. For $\kappa > 1$, an interaction range bonds particles that are not geometrically in contact but are
<table><caption>Table 1: Average number of bonds per particle for the 3 preparation routes at $D = 0.65$ for the 4 values of the interaction range $\kappa$</caption><thead><tr><th></th><th colspan="4">Z<sub>b</sub></th></tr><tr><th>preparation</th><th>κ = 1.000</th><th>κ = 1.025</th><th>κ = 1.050</th><th>κ = 1.075</th></tr></thead><tbody><tr><th>unjammed</th><td>4.52±0.02</td><td>5.00±0.01</td><td>5.49±0.02</td><td>5.99±0.01</td></tr><tr><th>weakly jammed to 0.5</th><td>5.39±0.02</td><td>5.81±0.01</td><td>6.23±0.01</td><td>6.56±0.01</td></tr><tr><th>jammed to 0.637</th><td>6.56±0.02</td><td>7.46±0.01</td><td>8.06±0.01</td><td>8.57±0.02</td></tr></tbody></table>

neighbors (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 $\kappa = 1.000$ and increases with increasing interaction range $\kappa$.

## 4. 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 $\nu$. The normal to tangential stiffness ratio ($\alpha = \Sigma_T/\Sigma_N$) has been varied from 0.005 to 10 to test the range of elastic constants attainable with potentially auxetic materials ($\nu < 0$ for $\alpha > 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:

$$ \tilde{E} = \frac{E}{Z_b D \Sigma_N} \qquad (11) $$

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 ($\alpha = \Sigma_T/\Sigma_N$) may be captured for the normalized Young's modulus by the following simple equation:

$$ \tilde{E} = \frac{\kappa^n a_1 + a_2 \alpha}{2\pi \left( 4 + a_3 \alpha \right)} \qquad (12) $$

where $n$ is a rational number, and $a_i$ are fitted parameters. Similarly, for the Poisson's ratio:

$$ \nu = \frac{1 - b_1 \alpha}{b_2 + b_3 \alpha} \qquad (13) $$

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 $\alpha$ values in between 0 and 1, which is the standard case. For this range, Fig. 5 shows the normalized Young's modulus multiplied by $2\pi\kappa^{-n}$ and the Poisson's ratio $\nu$ for the four values of interaction ranges $\kappa$ and with the three preparation methods. The value of $n$ is fitted to obtain an approximate master curve for all simulation points. Error bars are calculated from the standard deviations on ten simulations. The standard deviations are less than 0.05 and 0.01 for the normalized Young's modulus
and Poisson's ratio, respectively. They are thus barely discernible in Fig. 5, indicating small dispersion from one random sample to another. The fact that all curves fall on a master curve indicates that our simple normalization of the Youngs modulus by the coordination number and by the relative density (Eq. (11)) is effective.

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 $\alpha$ 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><caption>Table 2: Parameters $a_i$, $n$ and $b_i$ of Eqs. (12) and (13) fitted for the three preparation routes ($0 < \alpha \le 1$).</caption><thead><tr><th>preparation</th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>3</sub></th><th>n</th><th>b<sub>1</sub></th><th>b<sub>2</sub></th><th>b<sub>3</sub></th></tr></thead><tbody><tr><td>unjammed</td><td>0.18</td><td>5.45</td><td>9.64</td><td>5</td><td>0.24</td><td>3.38</td><td>6.91</td></tr><tr><td>weakly jammed to 0.5</td><td>0.42</td><td>7.45</td><td>8.68</td><td>2</td><td>0.55</td><td>3.02</td><td>5.75</td></tr><tr><td>jammed to 0.637</td><td>1.09</td><td>6.72</td><td>5.12</td><td>1</td><td>0.80</td><td>3.07</td><td>4.07</td></tr><tr><td>mean field</td><td>2</td><td>3</td><td>1</td><td>0</td><td>1</td><td>4</td><td>1</td></tr><tr><td>static hypothesis</td><td>0</td><td>10</td><td>6</td><td>0</td><td>1</td><td>2</td><td>3</td></tr></tbody></table>

(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 \le \alpha \le 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 $\alpha$ values than in the standard $0 < \alpha \le 1$ range (Fig. 6). The normalized Young's modulus is in relatively good accordance with the mean field model
<figure><img src="image_6.png" /><figcaption>Figure 5: Evolution of (a) the normalized Young's modulus $\tilde{E}$ and (b) the Poisson's ratio for the three preparation routes and for various values of the interaction range $\kappa$, in the $0 \le \alpha \le 1$ range.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Evolution of (a) the normalized Young's modulus $\tilde{E}$ and (b) the Poisson's ratio for the three preparation routes and for various values of the interaction range $\kappa$, in the $1 \le \alpha \le 10$ range.</figcaption></figure>
<table><caption>Table 3: Parameters $a_i$, $n$ and $b_i$ of Eqs. (12) and (13) fitted for the three preparation routes ($1 \le \alpha \le 10$, potentially auxetic materials).</caption><thead><tr><th>preparation</th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>3</sub></th><th>n</th><th>b<sub>1</sub></th><th>b<sub>2</sub></th><th>b<sub>3</sub></th></tr></thead><tbody><tr><td>unjammed</td><td>1.10</td><td>0.86</td><td>0.72</td><td>6</td><td>0.54</td><td>6.41</td><td>1.10</td></tr><tr><td>weakly jammed to 0.5</td><td>1.62</td><td>1.33</td><td>0.72</td><td>3</td><td>0.75</td><td>6.00</td><td>1.07</td></tr><tr><td>jammed to 0.637</td><td>2.06</td><td>2.07</td><td>0.78</td><td>1.5</td><td>0.91</td><td>5.09</td><td>0.99</td></tr></tbody></table>

for the fully jammed packing, while it deviates markedly from this model and from the static hypothesis model for the two other preparations. Also, the Poisson's ratio is further decreasing with increasing stiffness ratio $\alpha$ values down to negative values. Poisson's ratio of the order of -0.5 are obtained for $\alpha = 10$. Introducing the interaction range parameter $\kappa$ into the fitting Eq. (12) allows a master curve $\tilde{E}\kappa^{-n}$ to be obtained. However, this was not possible with the Poisson's ratio, and we prefer to present directly the Poisson's ratio values as a function of $\alpha$ without proposing a complex correlation between Poisson's ratio and $\alpha$.

## 5. 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 $\alpha = 0.2$ was chosen for the tangential to normal stiffness and the interaction range $\kappa$ 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]:

$$ \sigma_{Euler} = \frac{\pi^2 E}{4 (k \frac{l}{r})^2} \qquad (14) $$

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
<figure><img src="image_8.png" /><figcaption>Figure 7: Stress-strain curve examples for three values of the ratio $l/r$.</figcaption></figure>

solution:

$$ \sigma_{\text{Mazzilli}} = \left( 1 - \sqrt{1 - \frac{\pi^2}{\left(k \frac{l}{r}\right)^2}} \right) \frac{E}{2} \qquad (15) $$

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 $\sigma_{Euler}$ and $\epsilon_{Euler}$ respectively ($\epsilon_{Euler} = \sigma_{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
shown in Fig. 7, buckling is characterized by a clear maximum in the stress-strain curve. However, a non-linear region exists prior to this maximum. Still, we compare the critical Euler buckling stress to this maximum stress.

Fig. 8 summarizes the simulation results for the three preparation routes investigated here. The critical stress obtained from the simulation is normalized by $\sigma_{Euler}$. The randomness of the granular structure inevitably results in some dispersion illustrated by error bars in Fig. 8. The three preparation routes lead to larger maximum stresses than the Euler's formulation for large $l/r$. However, our definition of the critical stress (see typical stress-strain curves Fig. 7) should lead in any case to an upper bound. Also, the overestimation of the critical stress would suggest that dynamics may be at play in our simulations which attempt to approach the quasi-static limit. 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:

$$ \tilde{E}_{kin} = \frac{E_{kin}}{n \max(NR)} \qquad (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. $\tilde{E}_{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 values of $l/r$. For the largest slenderness, dynamics may explain partially the stress overestimation exhibited by DEM simulations.
<figure><img src="image_9.png" /><figcaption>Figure 8: Critical buckling stresses normalized by Euler's stress versus the $l/r$ ratio for the three preparation routes investigated. Mazzilli's solution (Eq. (15) (— —) is plotted for completeness [21].</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Normalized kinetic energy evolution versus normalized strain for the three values of the ratio *l/r*, as in Fig. 7.</figcaption></figure>
## 6. 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 ($\alpha$) 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]):

$$ \sigma_{ij} = \frac{Z_b D}{10\pi} (2 + 3\alpha) \left[ \varepsilon_{ij} + \frac{1 - \alpha}{2 + 3\alpha} \varepsilon_{kk} \delta_{ij} \right] \Sigma_N \quad (1) $$

The Poisson's ratio is then readily given by

$$ \nu = \frac{1 - \alpha}{4 + \alpha} \quad (2) $$

while the normalized Young's modulus simplifies to

$$ \tilde{E} = \frac{E}{Z_b D \Sigma_N} = \frac{1}{2\pi} \frac{2 + 3\alpha}{4 + \alpha} \quad (3) $$

An alternative expression has been derived by [17], based on the static hypothesis (or best-fit model). It leads to:

$$ \sigma_{ij} = \frac{Z_b D}{2\pi} \frac{5\alpha}{3 + 2\alpha} \left[ \varepsilon_{ij} + \frac{1 - \alpha}{5\alpha} \varepsilon_{kk} \delta_{ij} \right] \Sigma_N \quad (4) $$

Thus, for the static hypothesis the Poisson's ratio is

$$ \nu = \frac{1 - \alpha}{2 + 3\alpha} \quad (5) $$
while the normalized Young's modulus is given by

$$ \tilde{E} = \frac{5}{2\pi} \frac{\alpha}{2 + 3\alpha} \qquad (6) $$

We use generalized forms of Eqs. ((.2), (.3), (.5) and (.6)) in section 4 to fit the numerical results for the Poisson's ratio and the Young's modulus.
References

[1] Agnolin, I., Roux, J.N.. Internal states of model isotropic granular packings. I. Assembling process, geometry, and contact networks. Phys Rev E 2007;76:61302.

[2] André, D., Iordanoff, I., Charles, J.l., Néauport, J.. Discrete element method to simulate continuous material by using the cohesive beam model. Computer Methods in Applied Mechanics and Engineering 2012;213-216:113-125.

[3] Bathurst, R.J., Rothenburg, L.. Note on a random isotropic granular material with negative poisson's ratio. International Journal of Engineering Science 1988;26:373-383.

[4] Chajes, A.. Stability Behavior Illustrated by Simple Models. Journal of the Structural Division 1969;95:1153-1172.

[5] Christoffersen, J., Mehrabadi, M.M., Nemat-Nasser, S.. A micromechanical description of granular material behavior. Journal of Applied Mechanics 1981;48:339-344.

[6] Cundall, P.A., Strack, O.D.L.. A discrete numerical model for granular assemblies. Géotechnique 1979;29:47-65.

[7] Euler, L.. Methodus inveniendi lineas curvas maximi minimive proprietate gaudentes, sive solutio problematis isoperimetrici latissimo sensu accepti. Additamentum I. De curvis elasticis. MarcumMichaelem Bousquet 1744.
[8] Gaspar, N.. A granular material with a negative Poissons ratio. Mechanics of Materials 2010;42:673-677.

[9] Griffiths, D.V., Mustoe, G.G.W.. Modelling of elastic continua using a grillage of structural elements based on discrete element concepts. International Journal for Numerical Methods in Engineering 2001;50:1759-1775.

[10] Hedjazi, L., Martin, C.L., Guessasma, S., Della Valle, G., Dendievel, R.. Experimental investigation and discrete simulation of fragmentation in expanded breakfast cereals. Food Research International 2014;55:28-36.

[11] Jauffrès, D., Martin, C.L., Lichtner, A., Bordia, R.K.. Simulation of the elastic properties of porous ceramics with realistic microstructure. Modell Simul Mater Sci Eng 2012;20:45009.

[12] Jefferson, G., Haritos, G.K., McMeeking, R.M.. The elastic response of a cohesive aggregate - a discrete element model with coupled particle interaction. J Mech Phys Solids 2002;50:2539-2575.

[13] Jerier, J., Molinari, J.. Normal contact between rough surfaces by the Discrete Element Method. Tribology International 2012;47:1-8.

[14] Karagiozova, D., Alves, M.. Dynamic Elastic-Plastic Buckling of Structural Elements: A Review. Applied Mechanics Reviews 2008;61:040803.

[15] Kuzkin, V.A., Asonov, I.E.. Vector-based model of elastic bonds for simulation of granular solids. Phys Rev 2012;86:51301.
[16] Leblicq, T., Vanmaercke, S., Ramon, H., Saeys, W.. Discrete element modelling of bendable tubes. International Journal of Mechanical Sciences 2015;94-95:75-83.

[17] Liao, C., Chang, T., Young, D., Chang, C.. Stress-strain relationship for granular materials based on the hypothesis of best fit. International Journal of Solids and Structures 1997;34:4087-4100.

[18] Liu, W., Hong, J.W.. Discretized peridynamics for linear elastic solids. Computational Mechanics 2012;50:579-590.

[19] Martin, C.L.. Elasticity, fracture and yielding of cold compacted metal powders. J Mech Phys Solids 2004;52:1691-1717.

[20] Martin, C.L., Bordia, R.K.. Influence of adhesion and friction on the geometry of packings of spherical particles. Phys Rev E 2008;77:31307.

[21] Mazzilli, C.E.N.. Buckling and post-buckling of extensible rods revisited: A multiple-scale solution. International Journal of Non-Linear Mechanics 2009;44:200-208.

[22] Meimand, V.Z., Graham-Brady, L., Schafer, B.W.. Imperfection Sensitivity and Reliability Using Simple Bar-Spring Models for Stability. International Journal of Structural Stability and Dynamics 2013;13:1250075.

[23] Misra, A., Yang, Y.. Micromechanical model for cohesive materials based upon pseudo-granular structure. International Journal of Solids and Structures 2010;47:2970-2981.
[24] Potyondy, D.O., Cundall, P.A.. A bonded-particle model for rock. Int J Rock Mech Min Sci 2004;41:1329-1364.

[25] Satake, M.. Fabric tensor in granular materials. In: Vermeer, P., Luger, H., editors. Deformation and Failure of Granular Materials. Balkema; 1982. p. 63-68.

[26] Schenk, C., Schuëller, G.. Buckling analysis of cylindrical shells with random geometric imperfections. International Journal of Non-Linear Mechanics 2003;38:1119-1132.

[27] Schlangen, E., Garboczi, E.J.. Fracture Simulations of concrete Using Lattice Models: Computational Aspects. Engineering Fracture Mechanics 1997;57:319-332.

[28] Scholtès, L., Donzé, F.V.. A DEM model for soft and hard rocks: Role of grain interlocking on strength. Journal of the Mechanics and Physics of Solids 2013;61:352-369.

[29] Silling, S.A., Askari, E.. A meshfree method based on the peridynamic model of solid mechanics. Computers and Structures 2005;83:1526-1535.

[30] Storakers, B., Fleck, N.A., McMeeking, R.M.. The viscoplastic compaction of composite powders. J Mech Phys Solids 1999;47:785-815.

[31] Tan, Y., Yang, D., Sheng, Y.. Discrete element method (DEM) modeling of fracture and damage in the machining process of polycrystalline SiC. Journal of the European Ceramic Society 2009;29:1029-1037.
[32] Tavarez, F.A., Plesha, M.E.. Discrete element method for modelling solid and particulate materials. Int J Numer Meth Engng 2007;70:379-404.

[33] Torquato, S., Stillinger, F.H.. Multiplicity of Generation, Selection, and Classification Procedures for Jammed Hard-Particle Packings. J Phys Chem B 2001;105:11849-11853.

[34] Walton, K.. The effective elastic moduli of a random packing of spheres. Journal of the Mechanics and Physics of Solids 1987;35:213-226.

[35] Wang, Y., Alonso-Marroquin, F.. A finite deformation method for discrete modeling: Particle rotation and parameter calibration. Granular Matter 2009;11:331-343.

[36] Weber, J.. Recherches concernant les contraintes intergranulaires dans les milieux pulvérulents. Bulletin de liaison des Ponts et Chaussées 1966;1-20.