# Internal states, stress-strain behavior and elasticity in oedometrically compressed model granular materials

M.H. Khalili, Jean-Noël Roux, Jean-Michel Pereira, Matthieu Vandamme, Sébastien Brisard, Michel Bornert

► To cite this version:

M.H. Khalili, Jean-Noël Roux, Jean-Michel Pereira, Matthieu Vandamme, Sébastien Brisard, et al. Internal states, stress-strain behavior and elasticity in oedometrically compressed model granular materials. Particles 2015, Sep 2015, Barcelona, Spain. hal-01214733

HAL Id: hal-01214733

https://enpc.hal.science/hal-01214733

Submitted on 12 Oct 2015

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
# INTERNAL STATES, STRESS-STRAIN BEHAVIOR AND ELASTICITY IN OEDOMETRICALLY COMPRESSED MODEL GRANULAR MATERIALS

M.H. KHALILI, J.-N. ROUX, J.-M. PEREIRA, M. VANDAMME, S. BRISARD and M. BORNERT

Université Paris-Est, Laboratoire Navier (UMR 8205), CNRS, ENPC, IFSTTAR, F-77455 Marne-la-Vallée

**Key words:** Granular Materials, DEM, Oedometric compression, Elasticity

**Abstract.** The behaviour of a model granular material (an assembly of slightly polydisperse spherical beads, with Hertz-Mindlin elastic and frictional contacts) subjected to one dimensional (oedometric) compressions is studied by DEM simulations. We systematically investigate the influence of the (idealized) packing process on the microstructure and stresses in the initial, weakly confined equilibrium state. Such characteristics as density (ranging from maximally dense to moderately loose), coordination number (which might vary independently of solid fraction, especially in dense systems), fabric and stress anisotropies are monitored in oedometric loading cycles in which the major principal stress varies by up to 5 orders of magnitude. The evolution of the solid fraction (or the void ratio) versus the imposed vertical (principal) stress as observed in the loading and unloading paths, like in the case of isotropic compression [2] and unlike laboratory tests on sands, the behaviour shows only very limited plastic strain and is very nearly reversible in dense samples (which tend nevertheless to lose contacts in a loading cycle if the initial coordination number was large). The irreversibility observed in sands should thus be attributed to plasticity or damage within inter granular contacts. The anisotropy of the microstructure is described by the angular distributions of contacts and forces. It is explicitly linked to the stresses in the loading history, by semi-quantitative relations. One of the important characteristics measured during the compression is the ratio of lateral to controlled ('vertical') stress, $K_0$. We discuss conditions in which $K_0$ might be regarded as constant. We calculate, via a static (matrix) method [1], the complete tensor of elastic moduli, expressing response to very small stress increments about the transversely isotropic equilibrium states along the loading path.

## 1 INTRODUCTION

Granular materials are usually studied in soil mechanics at macroscale as a continuum medium. Recent works tend to investigate more the microstructural mechanism (contact
forces, fabric anisotropy, coordination number...) and relate them to macroscopic behavior. Discrete Element Method (DEM) [4], is a powerful tool to investigate the internal state of packing, and provides relevant variables such as fabric and coordination number [5]. In this paper we propose to study slightly polydisperse spherical beads under uni-axial oedometric compression. The evolution of structural parameters of the different assembly is studied to characterize the internal state during oedometric compression and probe the effect of the initial packing process. Strain reversibility and elasticity are also discussed. Note that other works have covered the oedometric compression by studying the effect of bead size distribution on internal states [6] or stress anisotropy compared to sand experiments [7].

## 2 DEM MODEL

In this study the sample is modeled by an assembly of 4000 spherical beads slightly polydisperse (from $d_1 = 1$ to $d_2 = 1.2$), which interact through Hertz-Mindlin elastic and frictional contacts. Glass beads were chosen with Young modulus $E = 70$ GPa, Poisson ratio $\nu = 0.3$ and friction coefficient $\mu = 0.3$. Beads are contained in a periodic cuboidal cell. Both grain positions and cell size are simultaneously determined on solving dynamical equations, implemented as in Ref. [1] for isotropic compression. In our case the applied stress and strain tensors describing the oedometric compression are $\sigma_1 = \sigma_0$ and $\varepsilon_2 = \varepsilon_3 = 0$. The strain rate is controlled during the compression and set to a constant value until the desired stress is reached. Then we search for the equilibrium state by allowing the strain rate to vary without exceeding a maximal value. The value of strain rate is set through the inertia parameter $I = \dot{\varepsilon}\sqrt{m/d\sigma_0}$, with $m$ the average particle mass. We also define a stiffness parameter $\kappa = (\tilde{E}/\sigma_0)^{2/3}$, with $\tilde{E} = E/(1 - \nu^2)$, such that the typical contact deflection is of order $\kappa^{-1}$.

## 3 SAMPLE PREPARATIONS

We investigate the effect of the initial state of the bead assembly on the oedometric compression response. Different sample preparation methods provide different microstructures characterized by their density (or solid fraction $\phi$), coordination number $z$ and fabric anisotropies. Bead diameters are randomly chosen according to a distribution uniform by volume:

$$p(d) = \frac{(d_1 d_2)^2}{2(d_2^2 - d_1^2)} \frac{1}{d^3}, \qquad (1)$$

Random, low density ($\Phi = 0.55$) configurations of non-contacting grains (granular gases) are then subjected to different packing assembly processes, involving either isotropic or oedometric (one-dimensional) compression to a solid state in which contact forces balance the applied stress.
**Moderately loose packing.** A weakly confined equilibrium state is obtained by applying a stress $\sigma_0$ of 10 kPa to the granular gas, until an equilibrium configuration is obtained with contact forces balancing the external load. Samples prepared in these conditions are moderately loose and will be noted Li or Lo according to whether they underwent an isotropic compression (all three principal stresses $\sigma_1 = \sigma_2 = \sigma_3$ set to $\sigma_0$) or an oedometric one ($\sigma_1 = \sigma_0$, no strain in directions 2 and 3).

**Maximally dense packing with high coordination number.** We again confine the granular gas at a confining pressure of 10 kPa. But this time, to maximize density [1], we set the friction coefficient $\mu$ to zero. This process gives an unrealistically high initial coordination number ($z$ about 6). In what follows, these configurations will be denoted DHi and DHo, for isotropic and oedometric compression.

**Maximally dense packing with low coordination number.** To obtain a more realistic coordination number in a dense packing process, we apply a small isotropic dilation to the confined maximally dense sample (the high coordination number one prepared without friction), so that all contacts open, and then confine the dilated configuration under pressure 10 kPa. A mixing step before the final compression phase [1] is necessary in the isotropic case. The final coordination numbers are about 4. We will denote these configurations as DLi (isotropic) and DLo (oedometric).

Tab. 1 characterizes the different initial states, giving values of $K_0$, the ratio of “horizontal” stress $\sigma_2 = \sigma_3$ to “vertical” stress $\sigma_1$ (denotes the applied stress in this paper), fraction $x_0$ of rattlers (grains carrying no force in equilibrated packs), coordination number $z$, and coordination number of non-rattler grains, $z^* = \frac{z}{1-x_0}$.

<table><caption>**Table 1:** Structural parameters of different initial states. Values are averaged over 3 configurations (differences for $\Phi$, $z$, $K_0$ lie below given accuracy level).</caption><thead><tr><th rowspan="2"></th><th colspan="2">Loose</th><th colspan="2">Dense High Coordination</th><th colspan="2">Dense Low Coordination</th></tr><tr><th>oedometric</th><th>isotropic</th><th>oedometric</th><th>isotropic</th><th>oedometric</th><th>isotropic</th></tr></thead><tbody><tr><th>K<sub>0</sub></th><td>0.72</td><td>1.0</td><td>0.94</td><td>1.0</td><td>0.51</td><td>1.0</td></tr><tr><th>φ</th><td>0.584</td><td>0.589</td><td>0.639</td><td>0.638</td><td>0.634</td><td>0.637</td></tr><tr><th>z</th><td>4.22</td><td>4.14</td><td>5.98</td><td>5.99</td><td>4.06</td><td>4.17</td></tr><tr><th>z<sup>*</sup></th><td>4.63</td><td>4.63</td><td>6.07</td><td>6.07</td><td>4.54</td><td>4.65</td></tr><tr><th>x<sub>0</sub>(%)</th><td>8.8</td><td>10.3</td><td>1.5</td><td>1.3</td><td>10.4</td><td>10.37</td></tr></tbody></table>

# 4 OEDOMETRIC COMPRESSION

## 4.1 Loading and unloading path

The different preparations of Sec. 3, are subjected to oedometric compressions, in which equilibrium states are obtained at different axial strain $\sigma_1$: 31.62 kPa, 100 kPa,
316.2 kPa, 1 MPa, 3.162 MPa, 10 MPa, 31.62 MPa (equidistant on a logarithmic scale). The calculation involves two steps. First, the system is compressed at a constant strain rate $\dot{\epsilon}_1$, maintaining inertia parameter $I$ to $10^{-5}$. Then, once the desired stress level $\sigma_1$ is obtained, it is kept constant until equilibrium criteria are satisfied for particles and for internal stress ($I$ is meanwhile requested to stay below $10^{-5}$). Upon unloading, the corresponding value of $I$ is set to $10^{-6}$.

## 4.2 Evolution of solid fraction $\phi$ and coordination number $z$

Fig. 1 shows the evolution of solid fraction during the oedometric compression of the different sample types. The curves present similar behaviours and curves of samples with the same initial solid fraction almost coincide on this scale – although highly coordinated systems are less deformable in the initial compression stages. Coordination number $z$

<figure><img src="image_2.png" /><figcaption>**Figure 1:** Evolution of $\phi$ with the vertical stress during the oedometric compression for samples of different preparation routines.</figcaption></figure>

(Fig. 2) slightly decreases in the first compression steps for samples DHi and DHo, for which its initial value is high (~ 6) due to the absence of friction in the assembling process. Friction stabilizes grain packs at lower coordination levels. $z$ increases again at higher pressures, as in all systems: new contacts are formed by closing the gaps between neighbours in denser structures [2]. Under high pressures, coordination numbers increase with density. The evolution of the coordination number depends on its initial value and the value of the solid fraction. Rattler fractions follow the opposite evolution to that of $z$ except for samples DHi and DHo for which the initial decrease of $z$ does not translate
into an increase of $x_0$.

<figure><img src="image_3.png" /><figcaption>**Figure 2:** Coordination number $z$ and rattler fraction $x_0$ versus vertical stress $\sigma_1$ along oedometric loading path for the different sample preparations.</figcaption></figure>

## 4.3 Variation of coefficient of lateral pressure $K_0$

In experiments on sands [8, 9, 10], $K_0$ is regarded as constant. As shown in Fig. 3 the variation in $K_0$ in our simulated bead packings strongly depends on the initial state. Obviously, samples initially confined by an isotropic compression of the granular gas are first devoid of structural anisotropy. Therefore, $K_0$ varies as the sample becomes anisotropic during the axial compression. This is the case of the samples Li, DHi and
DLi. Within this group, the denser samples reach higher anisotropy, the faster the larger the number of contacts (enabling different force distributions in contact networks).

The samples prepared under oedometric compression of the granular gas (such as DLo and Lo) are anisotropic from the beginning. $K_0$ is then observed to remain almost constant. As to DHo, prepared with no friction, its very small initial anisotropy is typical of frictionless systems [12]. The anisotropy of contacts and forces also depends on the preparation and will be detailed in the next section.

<figure><img src="image_4.png" /><figcaption>Figure 3: $K_0$ versus "vertical" stress $\sigma_1$ in oedometric compression.</figcaption></figure>

## 4.4 Anisotropy characterization

$K_0$ is characteristic of the overall degree of stress anisotropy of the system. At the microscopic level, stresses are determined by contact forces, which may differently contribute to $\sigma_1$ and $\sigma_2$, depending on their directions and their amplitudes. Anisotropies of the force network may thus be characterized by the distribution of contact orientations, and also by the angular dependence of the force distribution. We only consider here the normal component of the contact forces, as the tangential force contribution to normal stress components $\sigma_1$, $\sigma_2 = \sigma_3$ is found negligible. We denote as $\bar{n}$ the normal unit vector at contacts.

**Contact distribution.** Due to rotational invariance about the oedometric compression axis, the density distribution function $P(\bar{n})$ on the unit sphere $\Sigma$ depends on $|n_1|$ only, and
it and can be expanded in a series of Legendre polynomials of even order. The simplest non-isotropic form is truncated at order 2:

$$P(|n_1|) \approx 1 + \frac{5}{4}(3c_2 - 1)(3n_1^2 - 1) \quad (\text{with } c_2 = \langle n_1^2 \rangle) \qquad (2)$$

**Angular distribution of normal force amplitudes.** The distribution of the average force amplitude in the contacts with direction $\underline{n}$ depends on $|n_1|$ as well, and its expansion to second order in Legendre polynomials is analogous to (2). We define

$$f_2 = \frac{1}{4\pi} \int_{\Sigma} \mathcal{F}(|n_1|) n_1^2 d^2 \underline{n}$$

with notation $\mathcal{F}(|n_1|) = \langle F_N \rangle(|n_1|)/\langle F_N \rangle$. Here $\langle F_N \rangle$ is the overall average normal contact force and $\langle F_N \rangle(|n_1|)$ is the average normal force among contacts with normal direction $\underline{n}$.

## 4.5 Evolution of $c_2$ and $f_2$

Fig. 4 is a plot, versus $\sigma_1$, in the oedometric compression, of the second order anisotropy parameters $c_2$ and $f_2$ respectively characterizing the contact distribution and the force distribution. $c_2$ and $f_2$, similarly to $K_0$, depend on initial state. The large variation of $K_0$ observed in isotropically assembled samples appears to be due to the evolution of $f_2$, rather than that of $c_2$. The relation of $K_0$ to both anisotropy parameters $c_2$ and $f_2$ can be made more explicit, as shown in the following section.

## 4.6 Estimation of $K_0$ from anisotropy parameters

The normal force contribution to principal stresses is given by:

$$\sigma_{\alpha} = \frac{3z\phi <d>}{\pi <d^3>} <F_N> \int_S P(\underline{n}) \mathcal{F}(\underline{n}) n_{\alpha} n_{\alpha} d\Omega \qquad (3)$$

Using second order expansion (2), and neglecting the products of anisotropic terms, we get the following approximation

$$\sigma_{\alpha} \approx \frac{3z\phi <d>}{\pi <d^3>} <F_N> \left[ <n_{\alpha}^2> + \langle\langle F_N \rangle(|n_{\alpha}|)n_{\alpha}^2\rangle - \frac{1}{3} \right]. \qquad (4)$$

Knowing that $<n_2^2> = <n_3^2> = \frac{1}{2}(1 - <n_1^2>)$, one may finally estimate $K_0$ as:

$$K_0 = \frac{\sigma_{22}}{\sigma_{11}} \approx \frac{4 - 3(c_2 + f_2)}{6(c_2 + f_2) - 2}. \qquad (5)$$

This estimate of $K_0$ agrees well with measured values (see Fig. 5).
<figure><img src="image_5.png" /><figcaption>Figure 4: Anisotropy parameters $c_2$ and $f_2$ versus applied stress $\sigma_1$.</figcaption></figure>

## 4.7 Elastic moduli

The elastic moduli are computed in equilibrium configurations via the static method of Ref. [3]. Moduli $C_{11}$ and $C_{12}$, respectively relating increments of strain $\varepsilon_1$ to increments in stresses $\sigma_1$ and $\sigma_2$, are plotted in Fig. 6, as functions of axial stress $\sigma_1$. At low pressure $C_{11}$ depends on both density and coordination number. At high pressure, the coordination tends itself to be determined by the density, and the curves separate in two groups according to solid fraction. For $C_{12}$ only the dependence on the coordination number is observed. The results along unloading paths appear to depend on solid fraction only for $C_{11}$, and nearly coincide for $C_{12}$. This is to be correlated to the evolution of coordination
<figure><img src="image_6.png" /><figcaption>Figure 5: $K_0$ versus $\sigma_1$. Dots joined by solid line: measured values. Dots joined by dashed line: Eq. 5.</figcaption></figure>

number upon unloading, as shown below in Sec. 5.

# 5 REVERSIBILITY

## 5.1 Void ratio

The void ratio $e = -1 + 1/\Phi$, as plotted in Fig. 7 versus axial stress in the compression cycle, appears to be very nearly reversible, especially in dense samples. This behaviour is quite different from the typical experimental ones. In sands, the initial and final values of $e$ in such a cycle is of the order of 0.05, and the final value is almost equal to the lowest one reached in the loading phase [11]. The irreversibility observed in sands should thus be attributed to phenomena (plasticity or damage at contacts) which are not present in the present model. Reversiblity, in the simulations, is only apparent, and should not to be confused with elasticity. The right plot of Fig. 7 shows that coordination numbers change irreversibly in the compression cycle: samples with high initial $z$ values end up with much lower final values (as in the isotropic compressions of Ref. [2]). This also justifies the change in the elastic moduli in unloading observed above.

## 5.2 Elastic strain

From the curves in Fig. 6 we can fit a function $C_{11}(\sigma_1)$ which allows to estimate a total elastic deformation as:

$$ \varepsilon^{\text{el}} = \int_{\sigma_{\text{min}}}^{\sigma_{\text{max}}} \frac{d\sigma_1}{C_{11}(\sigma_1)} \qquad (6) $$
<figure><img src="image_7.png" /><figcaption>**Figure 6:** $C_{11}$ (top) and $C_{12}$ (bottom) versus $\sigma_1$ in oedometric compression (doubly logarithmic plots), along loading (left graph) and unloading (right graph) paths. Black dashed line has slope 1/3.</figcaption></figure>

The value of this elastic strain represents the expected strain in the sample if all the occurring evolution were governed by elastic mechanisms. Such elastic strains are smaller than measured strains, as reported in Tab. 2, with larger differences in looser states – thereby providing new evidence of nonelastic strain mechanisms. In the loose case (Lo) the difference between observed and elastic strain is somewhat reduced by considering the unloading stress path, but still of the order as the one obtained in denser systems.
<figure><img src="image_8.png" /><figcaption>**Figure 7:** Void ratio *e* (left) and coordination number *z* (right) versus *σ*₁ in oedometric compression cycle.</figcaption></figure>

<table><caption>**Table 2:** Total strains computed from elastic moduli (ε<sup>el</sup>) and measured in simulations along loading (ε<sup>ld</sup>) and unloading (ε<sup>uld</sup>) paths.</caption><thead><tr><th></th><th>ε<sup>el</sup></th><th>ε<sup>ld</sup></th><th>ε<sup>uld</sup></th></tr></thead><tbody><tr><th>Lo</th><td>1.72 %</td><td>3.54 %</td><td>2.58 %</td></tr><tr><th>DHo</th><td>1.27 %</td><td>2.28 %</td><td>2.38 %</td></tr><tr><th>DLo</th><td>1.38 %</td><td>2.22 %</td><td>2.14 %</td></tr></tbody></table>

# 6 CONCLUSION

This numerical study of granular media under oedometric compression shows that the evolution of $K_0$ depends on the packing process. It was observed that when care is taken to establish an initial anisotropy in the sample (thus mimicking such processes as air pluviation) the coefficient $K_0$ can be regarded as constant during loading. $K_0$ can be simply related to fabric (contact orientation) and normal force anisotropy parameters. We also showed that the irreversibility usually obtained for sands in real experiments is not present in our DEM simulations and hence should be attributed to other phenomena that were not included in the model such as grains damage and plasticity.

**Acknowledgment:** This work has benefited from a French government grant managed by ANR within the frame of the national program Investments for the Future ANR-11-LABX-022-01
## REFERENCES

[1] I. Agnolin and J.-N. Roux (2007) *Internal states of model isotropic granular packings*. *I. Assembling process, geometry, and contact networks*. Physical Review E, 76(6), 061302.

[2] I. Agnolin and J.-N. Roux (2007) *Internal states of model isotropic granular packings*. *II. Compression and pressure cycles*. Physical Review E, 76(6), 061303.

[3] I. Agnolin and J.-N. Roux (2007) *Internal states of model isotropic granular packings*. *III. Elastic properties*. Physical Review E, 76(6), 061304.

[4] P. A. Cundall and O. D. L. Strack (1979) *A discrete numerical model for granular assemblies*. Geotechnique 29,47.

[5] F. Radjaï and F. Dubois (2011). *Discrete-element modeling of granular materials* (pp. 425-p). Wiley-Iste.

[6] J. Wiacek and M. Molenda (2014) *Effect of particle size distribution on micro-and macromechanical response of granular packings under compression*. International Journal of Solids and Structures, 51(25), 4189-4195.

[7] Y. Gao and Y. H. Wang (2014) *Experimental and DEM Examinations of K0 in Sand under Different Loading Conditions*. Journal of Geotechnical and Geoenvironmental Engineering, 140(5).

[8] O. Daramola (1980) *On estimating Ko for overconsolidated granular soils*. Geotechnique, 30(3).

[9] J. A. Yamamuro, P. A. Bopp and P. V. Lade (1996) *One-dimensional compression of sands at high pressures*. Journal of Geotechnical Engineering, 122(2), 147-154.

[10] J. Feda (1984) *Ko-Coefficient of Sand in Triaxial Apparatus*. Journal of Geotechnical Engineering, 110(4), 519-524.

[11] A. Sanzeni and A. Whittle (2012) *Compression and Creep of Venice Lagoon Sands*. Journal of Geotechnical, (October), 1266-1276.

[12] P.-E. Peyneau and J.-N. Roux (2008) *Solidlike behavior and anisotropy in rigid frictionless bead assemblies*. Physical Review E, 78(4), 041307.