# Grading-Dependent Behavior of Granular Materials: From Discrete to Continuous Modeling

Gang Li, Ying-Jing Liu, Christophe Dano, Pierre-Yves Hicher

► To cite this version:

Gang Li, Ying-Jing Liu, Christophe Dano, Pierre-Yves Hicher. Grading-Dependent Behavior of Granular Materials: From Discrete to Continuous Modeling. Journal of Engineering Mechanics - ASCE, 2014, 141 (6), 10.1061/(ASCE)EM.1943-7889.0000866. hal-01006919

HAL Id: hal-01006919

https://hal.science/hal-01006919

Submitted on 11 May 2017

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
# Grading-Dependent Behavior of Granular Materials: From Discrete to Continuous Modeling

Gang Li<fnref n="1" />; Ying-Jing Liu<fnref n="2" />; Christophe Dano<fnref n="3" />; and Pierre-Yves Hicher<fnref n="4" />

**Abstract:** The aim of this paper is to investigate the effect of soil grading on the stress-strain and the critical-state behavior of granular materials from idealized spheres to natural soils as well as from discrete to continuous modeling. The three-dimensional discrete-element method has been applied to study the mechanical behavior of idealized granular materials. The results confirm that the critical-state line (CSL) shifts downward as grading broadens with an increase of the coefficient of uniformity $C_u$. An exponential relationship between the critical-state parameters and the coefficient of uniformity $C_u$ can then be established. Furthermore, experimental investigations on an artificial material (glass beads) and a natural material (Hostun sand) were carried out. The experimental results confirm the grading dependency of the critical-state behavior, and a similar relationship between CSL and $C_u$, for both glass beads and natural sand, extending previous literature results for high $C_u$ values up to 20. Moreover, the experimental results show that the CSL tends to stabilize when the coefficient of uniformity $C_u$ becomes higher than 10. Based on all test results, a simple elastoplastic model accounting for the influence of the grain size distribution (GSD) on the mechanical behavior of granular materials has been developed within the framework of critical-state soil mechanics. Comparisons between experimental results and simulations show that the model can accurately reproduce the mechanical behavior of granular materials with different GSDs.

**keywords:** Granular materials; Grain size distribution (GSD); Critical state; Discrete-element modeling (DEM); Triaxial testing; Elastoplastic model.

## Introduction

The mechanical behavior of granular materials is dependent on the properties of the grains and on the distribution of sizes [grain size distribution (GSD)] as well as on the geometrical arrangement of the grains. The change in the GSD upon loading attributable to grain breakage occurs in many geotechnical applications, such as within high rockfill dams (Indraratna et al. 1993; Ghanbari et al. 2008; Carrera et al. 2011), at the tip of foundation piles (Datta et al. 1979; Yasufuku and Hyde 1995; Yang et al. 2010), or in ballasted railway structures (Indraratna et al. 2005; Lobo-Guerrero and Vallejo 2006; Lu and McDowell 2010).

Engineering practice and research findings reveal that granular crushing, which induces GSD changes, has a significant effect on the mechanical behavior of granular materials: grain crushing generates a more contractive volumetric behavior concomitant with a decrease of the peak friction angle compared with the behavior of an intact material (Coop 1990; Verdugo and de la Hoz 2007).

Biarez and Hicher (1997), Daouadji et al. (2001), and Hu et al. (2011) suggested that, as the amount of grain breakage increases, the critical-state line (CSL) shifts downward in the e-log $p'$ plane according to the evolution of GSD possibly described by the increase of the coefficient of uniformity. Based on the numerical analysis of two-dimensional (2D) discrete-element modeling (DEM) simulations, Muir Wood and Maeda (2008) confirmed that the CSL of a well-graded material lies below that of a poorly graded material. Muir Wood et al. (2010) proposed a continuum model incorporating a relationship between a grading index $I_g$ and the CSL in the e-$p'$ plane. Yan and Dong (2011) performed three-dimensional (3D) DEM simulations on sphere assemblies with four different GSDs to study the GSD effect on the stress-strain and volumetric behavior of granular materials, and calibrated the model parameters for the continuum-based model proposed by Li and Dafalias (2000). They suggested a linear relationship between the critical-state parameters and the coefficient of uniformity. However, the study focused only on small values of the uniformity coefficient $C_u$ (smaller than 2.2), and the evidence presented is not so systematic. Furthermore, the calibrated parameters were found to be different for each grading.

The purpose of this study is to extend the previous research studies (Muir Wood and Maeda 2008; Yan and Dong 2011) and to render the results more general and more reliable by adopting different types of materials (DEM spheres, glass beads, and Hostun sand) with a wider range of GSDs ($C_u$ = approximately 1–20). In addition, the materials selected in this paper are made of perfectly rounded (glass balls) or highly self-similar (Hostun sand) grains, which eliminates or at least minimizes the influence of the grain shape on the critical state (Cho et al. 2006).

<footnote n="1">¹Engineer, LUNAM Univ., Ecole Centrale de Nantes, Research Institute in Civil and Mechanical Engineering (GeM), Unité Mixte de Recherche (UMR) 6183 of the Centre National de la Recherche Scientifique (CNRS), 6183 Nantes, France. E-mail: engineer.ligang@gmail.com</footnote>

<footnote n="2">²Engineer, Dept. of Civil Engineering, Shanghai Jiao Tong Univ., Shanghai 200240, P.R. China; formerly, Ph.D. Candidate, LUNAM Univ., Ecole Centrale de Nantes, Research Institute in Civil and Mechanical Engineering (GeM), Unité Mixte de Recherche (UMR) 6183 of the Centre National de la Recherche Scientifique (CNRS), 6183 Nantes, France. E-mail: liuyj301@gmail.com</footnote>

<footnote n="3">³Associate Professor, LUNAM Univ., Ecole Centrale de Nantes, Research Institute in Civil and Mechanical Engineering (GeM), Unité Mixte de Recherche (UMR) 6183 of the Centre National de la Recherche Scientifique (CNRS), 6183 Nantes, France (corresponding author). E-mail: christophe.dano@ec-nantes.fr</footnote>

<footnote n="4">⁴Professor, LUNAM Univ., Ecole Centrale de Nantes, Research Institute in Civil and Mechanical Engineering (GeM), Unité Mixte de Recherche (UMR) 6183 of the Centre National de la Recherche Scientifique (CNRS), 6183 Nantes, France. E-mail: pierre-yves.hicher@ec-nantes.fr</footnote>
From these results, the authors seek to establish a more general relationship between the critical-state parameters and the coefficient of uniformity $C_u$ for different types of materials. Then, this relationship will be introduced into an elastoplastic model developed within the framework of critical-state soil mechanics for predicting the behavior of the tested materials under triaxial loading. The aim is to propose, for any given material, a unique set of parameters for specimens with different GSDs.

It should be noted that the authors consider in this paper only granular materials of the sandy type. Investigating the influence of nonplastic fines content on the mechanical behavior of granular material is beyond the scope of this study.

## Grading Dependency of Idealized Granular Materials

### Principle of Numerical Simulations

The DEM numerical analysis was performed using Particle Flow Code in 3 Dimensions (PFC3D) 3.10 software (Cundall and Strack 1979). For the sake of simplicity, a linear contact model between particles was adopted. One of the long-term objectives of this study is to accurately describe the behavior of crushable soils and rockfills. For such applications, a Talbot-type GSD is usually assumed because it leads to compact earth fills in the field. This Talbot-type distribution has also been used to fit the fractal distribution when grain crushing occurs, as was suggested by Einav (2007). The GSDs are therefore characterized by the following function (Einav 2007):

$$F(d) = \left(\frac{d}{d_M}\right)^{\alpha} \qquad (1)$$

where $\alpha$ = dimensionless variable; and $d_M$ = maximum grain size. Changing the value of $\alpha$ produces different GSDs. The GSDs of the specimens used for the analysis are shown in Fig. 1. The relationship between $C_u$ and $\alpha$ is the following:

$$C_u = \frac{d_{60}}{d_{10}} = 6^{1/\alpha} \qquad (2)$$

### Specimen Creation

A specified number of spheres for each predetermined GSD (Fig. 1) was placed randomly in a parallelepiped container enclosed with six rigid walls (not shown in Fig. 2). An initial high porosity (high void ratio $e$ as listed in Table 1) of the specimen was chosen to obtain a very loose assembly with almost zero initial forces at the particle contacts. The spheres were then expanded to their final sizes. The

<figure><img src="image_2.png" /><figcaption>Fig. 1. Grain size distributions for DEM analyses</figcaption></figure>

system was iterated to equilibrium, whereas the wall location was maintained constant. The system was considered to be in quasi-static equilibrium when the ratio of the average unbalanced force to the average contact force became smaller than a target tolerance value (e.g., $10^{-3}$). A Cundall-type local damping ratio of 0.628 (Ng 2006) was used in this study to accelerate the equilibrium process. The tests were simulated under conditions of zero gravity. Detailed information for each specimen is listed in Table 1. The values of the parameters used in the DEM model are summarized in Table 2.

<figure><img src="image_3.png" /><figcaption>Fig. 2. (Color) Initial specimen of spheres</figcaption></figure>

<table><caption>Table 1. Information for DEM Specimens</caption><thead><tr><th>Grading</th><th>Void ratio <i>e</i></th><th>Specimen size (cm<sup>3</sup>)</th><th>Particle number</th></tr></thead><tbody><tr><td>C<sub>u</sub> = 1.0</td><td>0.835</td><td>96 × 192 × 96</td><td>3,588</td></tr><tr><td>C<sub>u</sub> = 1.2</td><td>0.835</td><td>96 × 192 × 96</td><td>4,633</td></tr><tr><td>C<sub>u</sub> = 1.4</td><td>0.835</td><td>88 × 176 × 88</td><td>5,433</td></tr><tr><td>C<sub>u</sub> = 1.6</td><td>0.835</td><td>80 × 160 × 80</td><td>5,274</td></tr><tr><td>C<sub>u</sub> = 1.8</td><td>0.835</td><td>80 × 160 × 80</td><td>8,170</td></tr><tr><td>C<sub>u</sub> = 3.6</td><td>0.667</td><td>48 × 96 × 48</td><td>28,676</td></tr><tr><td>C<sub>u</sub> = 6.0</td><td>0.667</td><td>36 × 72 × 36</td><td>28,196</td></tr></tbody></table>

<table><caption>Table 2. Properties of DEM Specimens</caption><thead><tr><th>Property</th><th>Sphere</th><th>Wall</th></tr></thead><tbody><tr><td>Density (kg/m<sup>3</sup>)</td><td>2,630</td><td>0.0</td></tr><tr><td>Normal contact stiffness (N/m)</td><td>1.5 × 10<sup>8</sup></td><td>1.5 × 10<sup>7</sup></td></tr><tr><td>Tangential contact stiffness (N/m)</td><td>1.0 × 10<sup>8</sup></td><td>1.5 × 10<sup>7</sup></td></tr><tr><td>Friction coefficient</td><td>0.5</td><td>0.0</td></tr><tr><td>Damping constant α</td><td>0.628</td><td>—</td></tr></tbody></table>
# Numerical Simulation Process

## Isotropic Compression Stage

The specimen was loaded through moving the six walls inwards or outwards and brought to an isotropic stress. A stress-controlled mode was used through a servocontrol mechanism by adjusting the wall movements. The stress was estimated from the summation of all contact forces acting on a parallel pair of walls, divided by the current function area of the walls. All stress components were determined.

The specimen in the container was first subjected to an isotropic confining stress of 15 kPa with an initial interparticle friction coefficient $\mu_{ini}$. The system was iterated to equilibrium. The initial interparticle friction coefficient $\mu_{ini}$ was carefully adjusted by trial and error to obtain the specified void ratio $e_0 = 0.687$ (in DEM simulations, $e_0$ means the initial void ratio after consolidation at a confinement pressure of 500 kPa) with a maximum error of 0.002 for all the specimens with different $C_u$ (the values of $\mu_{ini}$ are, respectively 0.0265, 0.0375, 0.085, 0.14, and 0.291 for the five first samples with $C_u$ between 1.0 and 1.8). Then, the interparticle friction coefficient was changed instantaneously from $\mu_{ini}$ to $\mu = 0.5$, and some cycles were run until the same applied isotropic stress state of equilibrium was reached. Isotropic consolidation to different confining levels can then be carried out. For any other simulation, $\mu_{ini} = 0.5$.

## Drained Triaxial Tests

After isotropic compression, the specimen was compressed in the vertical direction in a strain-controlled mode by specifying the

velocities of the top and bottom walls, giving a strain rate of 2%/s. The horizontal stress was maintained as nearly constant by continuously adjusting the width and length of the assembly as the vertical compression proceeded. The tests were conducted over a range of confining pressures from approximately 100–2,500 kPa, with an average relative overlap between particles $\Delta u/(r_1 + r_2)$ in the range of approximately 0.0020–0.0039 ($\Delta u$ is the overlap at the contact between two particles, and $r_1, r_2$ are the radii of the contacting spheres). Before presenting the simulation results, some issues relative to DEM have to be addressed.

# Some DEM Simulation Issues

## Number of Particles

Three different sample sizes with gradings $C_u = 1.0$ and $C_u = 6.0$ were performed to study the influence of the particle number on the stress-strain response, in particular at the critical state (Ng 2008). The results are shown in Fig. 3. It could be observed that the evolution of the deviatoric stress is very similar for each specimen of the same grading. On the other hand, differences can be seen in the void ratio evolution. The specimen size should be big enough for the result to converge toward similar behavior. Thus, the minimum size for the specimens can be determined through Fig. 3(b): 3,588 particles for the grading $C_u = 1.0$, and 28,196 particles for the grading $C_u = 6.0$. The oscillations in stresses are smaller for an increasing number of particles as shown in Figs. 3(c and d).

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** (Color) Influence of particle number on stress-strain response: (a and b) response for specimen with $C_u = 1.0$; (c and d) response for specimen with $C_u = 6.0$</figcaption></figure>
## Variability of Specimen Generation

Because the DEM specimens were randomly generated, different generations can produce different initial fabrics for the same initial porosity and confinement. These differences can affect the

mechanical behavior (Yimsiri and Soga 2011). To study the influence of the variability of specimen generation on the stress-strain response, two extreme GSDs ($C_u$ = 1.0 and $C_u$ = 6.0) were generated with three different random seeds (g1 = 10,000; g2 = 50,000;

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** (Color) Shear behavior comparison between different random seed generations: (a) deviatoric stress $q$ versus axial strain $\varepsilon_1$ for grading of $C_u$ = 1.0 and $C_u$ = 6.0; (b) volumetric strain $\varepsilon_v$ versus axial strain $\varepsilon_1$ for grading of $C_u$ = 1.0 and $C_u$ = 6.0</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Measurement of porosity distribution in specimen at different shearing stages (grading index $C_u$ = 1.0; initial void ratio = 0.687; cell pressure = 500 kPa): (a) front view of measurement spheres; (b) plan view of measurement spheres; (c) $\varepsilon_1$ = 2%; (d) $\varepsilon_1$ = 4%; (e) $\varepsilon_1$ = 8%; (f) $\varepsilon_1$ = 16%; (g) $\varepsilon_1$ = 24%; (h) $\varepsilon_1$ = 28%; (i) $\varepsilon_1$ = 32%; (j) $\varepsilon_1$ = 34%</figcaption></figure>
and g3 = 70,000, as defined in PFC3D notice) but with the same specimen size and particle numbers. They were then sheared under a confining pressure of 500 kPa.

The stress-strain responses are plotted in Fig. 4. For both specimens with C_u = 1.0 or C_u = 6.0, the stress-strain responses are almost identical for different generations, although the results present some fluctuations for high C_u specimens. The authors can conclude that, for different random generations, the specimen sizes used in this paper can be regarded as accurate enough for the following analysis. Considering the preceding results, the sensitivity to the random seed can be disregarded. Therefore, the generation random seed g1 was chosen in this study for all the specimen generations, which is also the default value in PFC3D.

## Potential Presence of Strain Localization

Researchers who used DEM to study the strain localization phenomenon adopted most of the time the condition of flexible boundaries. Under this condition, strain localization can occur, as shown, for example, by Iwashita and Oda (2000), O'Sullivan et al. (2003), or Jiang et al. (2011). However, a reflection of localized strain bands at rigid boundary was sometimes observed. To check the potential presence of strain localization, measurement spheres were used to characterize the void ratio distribution within the specimens (Jiang et al. 2011). Fig. 5 presents an example of the measured void ratios within a specimen. The locations of the measurement spheres are shown in Figs. 5(a and b), where no evident strain localization is observed. Each sphere represents a data point, and the measured results are projected to the middle section of the specimen. It is shown that, at large strains, the particle rearrangement

occurs almost over the entire granular assembly, which gives a rather homogeneous distribution of the void ratio within the specimen. Therefore, it seems reasonable to use the average stresses and strains to quantify the critical state of the assembly.

## Grading Dependent Stress-Strain Behavior

After these preliminary considerations, isotropically consolidated drained triaxial compression (CIDC) simulations were performed at the same consolidation stress of 500 kPa. Five specimens with different GSDs were created with an almost identical initial void ratio (0.687) after consolidation. Two additional specimens with different void ratios and gradings were also simulated to verify the absence of influence of the initial void ratio on the critical state. It is shown in Fig. 6 that at relatively large strains (higher than 25%), all specimens reach a stable state (constant stress and constant void ratio), which corresponds to the so-called critical state.

As in Yan and Dong (2011), Muir Wood and Maeda (2008), or Biarez and Hicher (1997), it is observed that the peak stress is particularly sensitive to the changes of GSD, decreasing as the C_u value increases [Fig. 6(a)]. For the specimens with a lower degree of uniformity (C_u = 1.0 and 1.2), a higher strength is reached at small strain levels and thereafter drops remarkably to a unique state. These specimens also need much larger strains to stabilize the volumetric strain [Fig. 6(b)]. Inversely, the constant volumetric strain state is reached sooner for specimens with higher degrees of uniformity (C_u = 1.8, 3.6, and 6.0). These results have to be analyzed in the framework of the state parameter concept $\psi = e - e_c$ (Been and Jefferies 1985), where $\psi$ is the state parameter, $e$ is the current void ratio, and $e_c$ is the critical-state void ratio. Considering the

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** (Color) Stress-strain responses of DEM analyses of specimens with different GSDs</figcaption></figure>
distribution characterized by $C_u = 1.57$ as an example [Fig. 6(c)], ($p'_0$, $e_0$) is the typical initial state point, and ($p'_c$, $e_c$) is the critical-state point for this specific grading. The sign of the state parameter determines the volumetric behavior of the specimens; hence, for the specimen with $C_u = 1.57$, the initial state parameter $\psi_0 = e_0 - e_c$ is negative, and during loading, the void ratio remains smaller than the $e_c$ until the critical state is reached (keeping $\psi < 0$). Therefore, the specimen exhibited dilatancy upon shearing. For the same reason, the specimen with $C_u = 1.0$ is also dilative. On the contrary, the specimen with $C_u = 3.6$ shows contraction because $\psi_0 = e_0 - e_c$ is positive. Hence, it can be concluded that the occurrence of a peak can be attributed to the material dilatancy, which depends on the relative position of the initial void ratio and the critical void ratio associated with a given grading, rather than to the GSD of the material. Besides, the results also show that specimens with different GSDs reached the same stress level at critical state regardless of the initial void ratio but at different critical void ratios.

## Grading Dependency of CSL

As previously mentioned, the deviatoric stress and void ratio became stable after a large strain level was reached, so the final state of shearing could be defined as the critical state. By measuring the critical states for the specimens of different GSDs and plotting them on the q-p' plane and e-p' plane, it was possible to define the CSLs as presented in Fig. 7. All simulations led to the conclusion that the critical stress ratio $\eta = q/p' = M$ is constant whatever the GSD may be, provided that the critical-state friction angles are similar [Fig. 7(a)]. An average value of $M = 0.737$ (equivalent to a friction angle of $19.2^{\circ}$) was determined. Because the stress-strain curves were not smooth, the stress plateau was averaged over the final 5% axial strains.

Unlike the final stress state in the q-p' plane, the location of the CSL on the e-p' plane depends on the GSD of the material. As demonstrated in Fig. 7(b), the CSL shifts downward on the compression e-p' plane, and its gradient decreases slightly when the $C_u$ value of the specimen increases.

Li and Wang (1998) have suggested the following form for the location of the CSL on the e-p' plane:

$$e_c = e_{\text{ref}} - \lambda \left( \frac{p'}{p_{\text{at}}} \right)^{\xi} \qquad (3)$$

where $e_c$ = critical-state void ratio at the mean effective stress $p'$; reference void ratio $e_{\text{ref}}$ allocates the position of the CSL; $\lambda$ controls the slope of the CSL; and $\xi$ = material constant. Thus, the CSL for a given material can be determined by the two parameters, $e_{\text{ref}}$ and $\lambda$, because $p_{\text{at}}$ and $\xi$ are assigned; in this study, $p_{\text{at}} = 101.3$ kPa, and $\xi = 0.9$. This value of $\xi$, calibrated from the experimental results, is similar to the value proposed by Yan and Dong (2011).

Based on the CSL family, the parameters $e_{\text{ref}}$ and $\lambda$ are plotted versus the GSD index $C_u$, respectively (Fig. 8). A similar exponential relation with $C_u$ is derived for both $e_{\text{ref}}$ and $\lambda$ as follows:

$$e_{\text{ref}} = a_{e,\text{ref}} + b_{e,\text{ref}}e^{(-c_{e,\text{ref}}C_u)} \qquad (4a)$$

$$\lambda = a_{\lambda} + b_{\lambda}e^{(-c_{\lambda}C_u)} \qquad (4b)$$

where $a_{e,ref}$, $b_{e,ref}$, $c_{e,ref}$, $a_{\lambda}$, $b_{\lambda}$, and $c_{\lambda}$ = fitting parameters.

The numerical results confirm the general scheme empirically proposed by Biarez and Hicher (1997) in that an almost exponential decrease of the critical void ratio with coefficient of uniformity up to

<figure><img src="image_8.png" /><figcaption>Fig. 7. CSLs of DEM analysis: (a) in q-p' plane; (b) in e-p' plane</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Relationship between CSL and $C_u$: (a) $e_{\text{ref}}$ and $C_u$; (b) $\lambda$ and $C_u$</figcaption></figure>
$C_u = 10$ was established. For $C_u$ values greater than 10, the critical void ratio remains almost constant. The values of $e_{ref}$ for $C_u = 1.0$ are higher when the grains are more angular (Biarez and Hicher 1994).

Experimental Evidence of Grading Dependency of Glass Beads and Natural Sand

It is widely accepted that DEM simulations can reproduce qualitatively the mechanical behavior of granular materials. Hence, more

<figure><img src="image_10.png" /><figcaption>Fig. 9. Grain shape of Hostun sand</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Grain size distributions for experimental testing</figcaption></figure>

experimental results supporting the computational predictions are needed. So far, however, only very few experimental data have been able to confirm the previous numerical findings. In this section, two different granular materials were chosen to verify the assumptions deduced from the preceding simulations: artificial spherical glass beads and natural Hostun sand of similar particle shapes for grains with different sizes, as shown by scanning electron microscope (SEM) observations in Fig. 9. The GSDs adopted for both glass beads and Hostun sand are plotted in Fig. 10, with almost the same mean grain size $d_{50}$, and the maximum grain size $d_M$ is close to 2 mm in all cases.

All specimens were prepared with the moist tamping method (Ishihara 1993) to obtain very loose specimens, so as to avoid early strain localization as much as possible. The CIDC tests were performed on specimens of 100 mm in diameter and 200 mm in height, under confining pressures of 100, 200, and 400 kPa. The strength of two material grains are relatively high so that grain breakage could be negligible under the loading applied (Biarez and Hicher 1997), and hence the same GSD could be maintained along the stress paths.

Test results illustrated a common mechanical behavior for loose specimens, namely a monotonic increase of the shear stress up to a plateau and a mainly contractive volumetric behavior up to the critical state. The CSLs obtained for different $C_u$ in the $p'-q$ plane are presented in Fig. 11(a) for glass beads and in Fig. 12(a) for Hostun sand. As previously stated, the tested specimens exhibited a unique critical stress ratio ($M = q/p'$) upon compression regardless of GSD ($M = 0.81$ for glass beads and $M = 1.13$ for Hostun sand). The uniqueness of the critical stress ratio confirms that the shear resistance is not influenced by the change in GSD.

The CSL on the $e-p'$ plane for specimens with different $C_u$ under different confining pressures is presented in Fig. 11(b) for glass beads and Fig. 12(b) for Hostun sand. It has been found that the location of the CSL varies as $C_u$ changes, the higher the $C_u$ value and the lower the location of the CSL. In other words, a well-graded material has a smaller void ratio at critical state than a more uniform material under the same stress level.

The findings from experiments on glass beads and Hostun sand have confirmed that the CSL is grading dependent, moving downward when the $C_u$ value of the material increases, which strongly confirms the DEM results (Muir Wood et al. 2010; Yan and Dong 2011). Similar relationships as those derived from the DEM results for both $e_{ref}$ and $\lambda$ with $C_u$ are obtained and plotted in Fig. 13.

Constitutive Modeling

From the preceding discussion, both DEM simulations and experimental test results have confirmed that the CSL shifts downward in

<figure><img src="image_12.png" /><figcaption>Fig. 11. CSLs of glass beads: (a) in $q-p'$ plane; (b) in $e-p'$ plane</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12. CSLs of Hostun sand: (a) in q-p' plane; (b) in e-p' plane</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Relationship between (a) e ref and uniformity coefficient C u ; (b) λ and uniformity coefficient C u</figcaption></figure>

the e-p' plane as the grading index C_u increases. Although DEM analysis can reflect the grading effect on the mechanical behavior of granular materials, its application to full-scale physical boundary-value problems is limited because of the heavy computational effort required. Thus, soil constitutive models under the continuum solid mechanics framework are still one of the most preferred tools for solving geotechnical problems. However, few constitutive models have considered the initial GSD effect (Muir Wood et al. 2010; Yan and Dong 2011). In this study, a specific constitutive model for granular materials was used to illustrate that the physical origin of the model constants could be better understood through a comparison with DEM simulation results and experimental test results. The constitutive model was developed in the same framework as the Muir Wood et al. (2010) model. However, some differences were inserted in the formulation of the CSL and in the choice of the parameters that governs the evolution of the GSD. The grading index used in the Muir Wood model assumes an ultimate fractal GSD, which, although some elements have been presented in the literature, remains very difficult to identify experimentally. Therefore, the calibration of this ultimate GSD requires the assumption of some parameters. Fortunately, granular materials exhibit common characteristics, the fractal parameter (2.6, corresponding to a C_u of approximately 88) among them. In the model presented in this paper, there is no need for such an assumption.

Elastoplastic Model

The elastic behavior is assumed to be isotropic and nonlinear

$$
d\epsilon_v^e = \frac{dp}{K}; \quad d\epsilon_d^e = \frac{dq}{3G} \qquad (5)
$$

where dε<sup>e</sup><sub>v</sub> = elastic volume strain increment; dε<sup>e</sup><sub>d</sub> = elastic deviatoric strain increment; and G and K = hypoelastic shear modulus and bulk modulus, respectively, defined as follows (Richart et al. 1970):

$$
G = G_0 \frac{(2.97 - e)^2}{(1 + e)} \left(\frac{p'}{p_{at}}\right)'' \quad (6)
$$

$$
K = K_0 \frac{(2.97 - e)^2}{(1 + e)} \left(\frac{p'}{p_{at}}\right)'' \quad (7)
$$

where G₀, K₀, and n = elastic parameters; and pₐt = atmospheric pressure used as reference pressure (pₐt = 101.3 kPa).

For the plastic behavior, the proposed model has two yield sur-
faces, one for shear sliding and the other for normal compression.
Thus, the framework of the proposed approach is similar to that of
the double-hardening model developed by Vermeer (1978).

As in many constitutive models for sand, the shape of the yield
surface for the shear component is linear in the p'-q plane, written as
follows:

$$
f_s = \eta - H \tag{8}
$$

where η = q/p'; and H = hardening variable defined by a hyperbolic function in the H-ε<sup>p</sup><sub>d</sub> plane, similar to the following one proposed by Yin et al. (2010):

$$
H = \frac{M_p \epsilon_d^p}{(M_p p' / G_p K) + \epsilon_d^p} \quad (9)
$$
where $\epsilon_d^p$ = deviatoric plastic strain; and $G_p$ is used to control the initial slope of the hyperbolic curve $\eta-\epsilon_d^p$. The Eq. (9) guarantees that the stress ratio will approach peak value at the stress ratio $M_p$.

To take into account dilation or contraction during shear, the following Roscoe-type stress dilatancy equation is used:

<table><caption>**Table 3.** Values of Model Parameters</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>DEM</th><th>Glass beads</th><th>Hostun sand</th></tr></thead><tbody><tr><td>Friction angle</td><td>φ<sub>μ</sub> (degrees)</td><td>19.2</td><td>21.0</td><td>28.4</td></tr><tr><td rowspan="3">Elasticity</td><td>G<sub>0</sub> (kPa)</td><td>263</td><td>60</td><td>34</td></tr><tr><td>K<sub>0</sub> (kPa)</td><td>350</td><td>80</td><td>45</td></tr><tr><td>N</td><td>0.5</td><td>0.5</td><td>0.5</td></tr><tr><td>Plasticity</td><td>G<sub>p</sub></td><td>0.004</td><td>0.004</td><td>0.004</td></tr><tr><td>Dilatancy</td><td>D</td><td>1.0</td><td>1.0</td><td>0.8</td></tr><tr><td rowspan="3">GSD effect on e<sub>ref</sub></td><td>a<sub>e,ref</sub></td><td>0.354</td><td>0.338</td><td>0.590</td></tr><tr><td>b<sub>e,ref</sub></td><td>0.624</td><td>0.42</td><td>0.181</td></tr><tr><td>c<sub>e,ref</sub></td><td>0.356</td><td>0.201</td><td>0.123</td></tr><tr><td rowspan="3">GSD effect on λ</td><td>a<sub>λ</sub></td><td>5.61 × 10<sup>-4</sup></td><td>4.82 × 10<sup>-4</sup></td><td>4.62 × 10<sup>-3</sup></td></tr><tr><td>b<sub>λ</sub></td><td>4.39 × 10<sup>-3</sup></td><td>2.04 × 10<sup>-3</sup></td><td>5.79 × 10<sup>-3</sup></td></tr><tr><td>c<sub>λ</sub></td><td>1.85</td><td>0.265</td><td>0.139</td></tr></tbody></table>

$$ \frac{d\epsilon_v^p}{d\epsilon_d^p} = D(M_{pt} - \eta) \qquad (10) $$

where $D$ = material parameter; and $M_{pt}$ = slope of the phase transformation line for sand as defined by Ishihara et al. (1975) or the characteristic line defined by Luong (1982).

The material's density state is defined by the ratio $e_c/e$, where $e$ is the current void ratio and $e_c$ is the critical void ratio for the current value of $p'$ obtained by Eq. (3). According to Biarez and Hicher (1997), the peak friction angle $\phi_p$ [related to $M_p = 6 \sin \phi_p / (3 - \sin \phi_p)$ for triaxial compression test] is linked to the intrinsic friction angle $\phi_μ$ [related to the critical-state value $M = 6 \sin \phi_μ / (3 - \sin \phi_μ)$ for triaxial compression test] and the density state $e_c/e$

$$ e \tan \phi_p = e_c \tan \phi_\mu \qquad (11) $$

The $M_p$ is then obtained through the critical-state $M$ and the density state $e_c/e$. The Eq. (11) means that in a dense state, the peak friction angle $\phi_p$ is greater than $\phi_μ$ owing to a higher degree of interlocking. When the stress state reaches the phase transformation line, the dense assembly dilates and the degree of interlocking decreases. As

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** (Color) Model parameter calibration of drained compression tests on DEM sphere, glass beads, and Hostun sand specimens: (a–c) isotropic compression tests; (d–f) deviatoric stress versus axial strain; (g–i) void ratio versus mean effective stress</figcaption></figure>
a consequence, the peak friction angle is reduced, which results in a strain-softening phenomenon.

The $M_{pt}$ [related to $6 \sin \phi_{pt}/(3 - \sin \phi_{pt})$] is the slope of the phase transformation line, which is assumed to be a function of the void ratio

$$e_c \tan \phi_{pt} = e \tan \phi_{\mu} \quad (12)$$

The Eq. (12) indicates that a dense packing has a smaller phase transformation angle than does a loose packing. The Eqs. (11) and (12) were successfully applied by Yin and Chang (2011).

## Calibration of Model Constants

The proposed model requires the calibration of 12 parameters. The material parameters can be obtained mainly from standard laboratory tests. The first part of the calibration requires triaxial compression tests at different values of the initial void ratio and different confining pressures to calibrate different features of hardening/softening and dilatancy/contractancy for the model that are associated with the plastic strain rate contribution attributable to stress-ratio changes. The other part of the calibration incorporates variables for describing the critical-state response under different GSDs. The determination of the model parameters for each material requires only one drained or undrained triaxial test up to failure with an initial isotropic consolidation stage. The model parameters are determined as follows:

* $K_0$ and $n$ were calibrated from an isotropic compression test, and $G_0$ can be calibrated from the initial slope of the stress-strain curve (e.g., $\epsilon_a < 0.1\%$) of drained triaxial compression tests (in this study, the authors assumed a Poisson's ratio $\nu = 0.2$ to calculate $G_0$ from $K_0$);

* $G_p$ was obtained by fitting the initial slope of the curve $q-\epsilon_a$ (for $\epsilon_a < 1\%$) in drained triaxial tests;

* $\phi_{\mu}$ was determined from drained triaxial tests; and

* $a_{e,ref}$, $b_{e,ref}$, $c_{e,ref}$, $a_{\lambda}$, $b_{\lambda}$, and $c_{\lambda}$ were obtained by curve fitting using Eqs. (4a) and (4b), which account for the effect of gradings in terms of $C_u$ value.

All the parameters in Table 3 were calibrated from the results of three drained triaxial tests (Fig. 14). The DEM specimen was sheared at a confining stress of 500 kPa after an isotropic consolidation stage; the glass beads and the Hostun sand specimens were sheared at the confining stress of 400 kPa after an isotropic consolidation stage.

## Test Simulations of Stress-Strain Responses

Using the sets of calibrated parameters in Table 3, the predicted results for DEM spheres, glass beads, and Hostun sand are shown in the following sections.

### Idealized Material

Fig. 15 compares DEM results and model simulations for drained triaxial tests in compression with a confining stress equal to 500 kPa. As shown in Fig. 15(a), the simulations reproduced the GSD effect

<figure><img src="image_16.png" /><figcaption>Fig. 15. (Color) Comparisons between DEM results and model simulations for drained triaxial compression tests: (a) deviatoric stress versus axial strain; (b) void ratio versus effective stress</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. (Color) Comparisons between glass beads and simulations for drained triaxial compression tests: (a) deviatoric stress versus axial strain; (b) void ratio versus effective stress</figcaption></figure>
on the stress-strain responses of granular materials: the peak deviatoric stress decreases as $C_u$ increases. The model also reproduces the dilatancy behavior obtained in DEM results [Fig. 15(b)]. For specimens with the same initial void ratio, the dilatancy decreases as $C_u$ increases: for large $C_u$ and specimens with different initial void ratios, the contractive behavior is well reproduced. The model, as shown in the comparison, describes with good accuracy the stress-strain responses of idealized granular materials.

## Glass Beads

Fig. 16 compares the experimental data and simulation results for drained triaxial compression tests on isotropically consolidated specimens of glass beads under a confining stress of 100 kPa. The model is able to reproduce the main features of the mechanical behavior of glass beads influenced by GSD. The small $C_u$ specimens

present a slight dilatancy, whereas the high $C_u$ specimens are contractive.

## Hostun Sand

Fig. 17 compares the data and simulation results for drained triaxial compression tests on isotropically consolidated specimens of Hostun sand with confining stresses varying from 100 to 400 kPa. Because the specimens are initially loose, the stress-strain responses show a contractive behavior. The specimens reach a stable stress state at approximately 20% in axial strain. In the e-p' plane, the void ratios decrease during compression and tend to become constant after 20% in axial strain. The numerical simulations show a good agreement with the experimental results, especially for the samples with a relatively small $C_u$, whereas the model is not able to reproduce in a good manner the peak observed in the stress-strain behavior for

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** (Color) Comparisons between experimental results of Hostun sand and simulations for stress-strain behavior of drained triaxial compression tests for (a and b) $C_u = 1.1$; (c and d) $C_u = 1.4$; (e and f) $C_u = 2.5$; (g and h) $C_u = 5$; (i and j) $C_u = 10$; (k and l) $C_u = 20$</figcaption></figure>
$C_u = 20$. Nevertheless, the proposed model can accurately simulate the granular material behavior, from idealized spheres to natural sand, assuming a Talbot-type GSD.

## Undrained Behavior

Simulations were carried out to evaluate the model's performance for predicting the undrained behavior of granular materials. Fig. 18 compares experimental data and numerical simulations of undrained triaxial compression tests on isotropically consolidated specimens of glass beads and Hostun sand at the confining stress of 400 kPa. The results show an overall good agreement between the experimental and numerical results. Therefore, the proposed model can accurately simulate the granular material behavior, from idealized spheres to natural sand of a Talbot-type GSD, taking into account the characteristics of the GSD.

## Discussion

To demonstrate the predictive ability of the proposed model, comparisons between simulations by the Muir Wood model (Muir Wood et al. 2010), denoted as Model-$I_g$, and the proposed model, denoted as Model-$C_u$, were performed. The main differences between the two models are the equation of the CSL and the effect of grading on the position of CSL

$$v_{cs} = (v_{\text{ref}} - \Delta v) + \Delta v \exp \left[ - \left( \frac{p'}{p_{cs}} \right)^{\beta} \right] \quad (13)$$

$$v_{\text{ref}} = v_0 - v_g I_g + \Delta v \quad (14)$$

where $v_{cs}$ = specific critical volume; $v_{\text{ref}}$ = reference position of the CSLs in the $v-p'$ plane when $p' = 0$ as defined in Eq. (14); $I_g$ = grading index; $v_0$ = value of the minimum void ratio for

a monodisperse material $I_g = 0$; $v_g \approx v_0 - 1$; $\Delta v$ = material parameter; $p_{cs}$ = reference stress; and $\beta$ controls the shape of the CSL. Because $\Delta v$ is constant, $v_{\text{ref}}$ has a linear relationship with the grading index $I_g$.

Taking tests on glass beads as an example, Eq. (13) is used to fit the critical-state data of glass beads with the following parameters: $p_{cs} = 506.5$ kPa, $\beta = 0.145$, and $\Delta v = 0.3$, referring to Muir Wood et al. (2010) and Kikumoto et al. (2010) [see Fig. 19(a)]. Then, the relationship between $e_{\text{ref,MW}}$ and $I_g$ was fitted by $v_{\text{ref}} = 1.57 - 0.57I_g + \Delta v$ with $v_0 = 1.57$, $v_g = 0.57$ [see Fig. 19(b)].

Beside the critical-state parameters, the other parameters of Model-$I_g$ were determined from simulating a drained triaxial test with $C_u = 1.1$ or $I_g = 0.0043$ sheared under a confining stress of 400 kPa: $a = 0.0016$, $M = 0.81$, $A = 0.9$, $k_R = 2$, and $k_D = 0.1$, also referring to Gajo and Muir Wood (1999) and Kikumoto et al. (2010), as illustrated in Fig. 20. Then, using these parameters, tests on three different gradings of glass beads and two confining pressures of 100 and 400 kPa were simulated and compared with simulations by the proposed model as shown in Fig. 21. All results show that both models have captured the trend of the stress-strain response of glass beads with different gradings. The difference mainly comes from the concise but also crude linear $v_{\text{ref}}-I_g$ relationship implemented in the Model-$I_g$ and the exponential, so nonlinear, $e_{\text{ref}}-C_u$ relationship in the Model-$C_u$.

## Conclusions

Three series of drained triaxial compression tests were performed to investigate the GSD effect on the stress-strain and the critical-state behavior of granular materials, from idealized spheres to natural sand as well as from discrete to continuous modeling. A wide range of gradings, all of Talbot-type, for different materials was studied: for DEM simulations, $C_u$ values were up to 6, whereas for glass beads and Hostun sand, the $C_u$ was up to 20.

<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Comparison between experimental results and model predictions for undrained tests on glass beads and Hostun sand: (a and b) deviatoric stress versus axial strain; (c and d) effective stress path</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 19. Fitted locations of CSLs using Model-$I_g$ for glass beads: (a) CSL of glass beads; (b) relationship between $v_{ref}$ and $I_g$</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 20. (Color) Determination of model parameters for Model-$I_g$: (a) deviatoric stress versus axial strain; (b) void ratio versus mean effective stress</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 21. (Color) Comparison between Model-$I_g$ and Model-$C_u$ on simulating glass beads results: (a) deviatoric stress versus axial strain for $C_u = 5.0$; (b) void ratio versus mean effective stress for $C_u = 5.0$; (c) deviatoric stress versus axial strain for $C_u = 20$; (d) void ratio versus mean effective stress for $C_u = 20$</figcaption></figure>
The results confirm that for the same initial void ratio, specimens with a wider GSD show more contractive behavior and strain hardening upon shearing. The GSD does influence the critical state of the granular materials. On the one hand, a larger $C_u$ value leads to a lower location of the CSL in the e-p' plane; on the other hand, the test specimens with different GSDs under the same confining stress reach almost the same stress ratio at critical state. In other words, the GSD affects the void ratio at critical state, but not the ultimate shear strength or friction angle. Furthermore, the results demonstrate that when the coefficient of uniformity $C_u$ is higher than 10, the position of the CSL tends to become stable in the e-p' plane. A new exponential relationship between the critical-state parameters and the coefficient of uniformity $C_u$ was established for the different materials, which extends the scope of previous studies to incorporate higher $C_u$ values.

A simple elastoplastic model accounting for the influence of GSD changes was developed within the framework of critical-state soil mechanics. Drained triaxial compression tests on DEM materials, glass beads, and Hostun sand were used to calibrate and validate the model. The model unifies all the gradings with only one group of parameters for one material and can predict drained and undrained triaxial compression tests. All comparisons between the experimental results (including DEM results) and the simulations demonstrate the capacity of the model to reproduce with good accuracy the mechanical behavior of granular materials with different GSDs under drained and undrained triaxial loading.

## Acknowledgments

The authors gratefully acknowledge the financial support provided by the European Program TANDEM, the Natural Science Foundation of China/French National Research Agency (NSFC/ANR) Joint Research Grant (No. 51161130523), the China Scholarship Council, and the Natural Science Foundation of China (Grant Nos. 41372285, 41330633, 41172251, 41002095, and 41272317). The authors thank Professor Zhen-Yu Yin of Shanghai Jiao Tong University and Dr. Jidong Zhao of Hong Kong University of Science and Technology for their advice and support.

## References

Been, K., and Jefferies, M. G. (1985). "A state parameter for sands." *Géotechnique*, 35(2), 99-112.

Biarez, J., and Hicher, P. Y. (1994). *Elementary mechanics of soil behaviour: Saturated remoulded soils*, Balkema, Rotterdam, Netherlands.

Biarez, J., and Hicher, P. Y. (1997). "Influence de la granulométrie et de son évolution par ruptures de grains sur le comportement mécanique de matériaux granulaires." *Rev. Fr. Genie Civ.*, 1(4), 607-631 (in French).

Carrera, A., Coop, M., and Lancellotta, R. (2011). "Influence of grading on the mechanical behaviour of Stava tailings." *Géotechnique*, 61(11), 935-946.

Cho, G.-C., Dodds, J., and Santamarina, J. C. (2006). "Particle shape effects on packing density, stiffness, and strength: Natural and crushed sands." *J. Geotech. Geoenviron. Eng.*, 10.1061/(ASCE)1090-0241(2006)132:5(591), 591-602.

Coop, M. (1990). "The mechanics of uncemented carbonate sands." *Géotechnique*, 40(4), 607-626.

Cundall, P. A., and Strack, O. D. L. (1979). "A discrete numerical model for granular assemblies." *Géotechnique*, 29(1), 47-65.

Daouadji, A., Hicher, P.-Y., and Rahma, A. (2001). "An elastoplastic model for granular materials taking into account grain breakage." *Eur. J. Mech. A, Solids*, 20(1), 113-137.

Datta, M., Gulhati, S., and Rao, G. (1979). "Crushing of calcareous sands during shear." *Proc.*, 11th Annual Offshore Technology Conf., Offshore Technology Conference, Houston.

Einav, I. (2007). "Breakage mechanics—Part I: Theory." *J. Mech. Phys. Solids*, 55(6), 1274-1297.

Gajo, A., and Muir Wood, D. (1999). "Severn-Trent sand: A kinematic-hardening constitutive model: The q-p formulation." *Géotechnique*, 49(5), 595-614.

Ghanbari, A., Sadeghpour, A. H., Mohamadzadeh, H., and Mohamadzadeh, M. (2008). "An experimental study on the behavior of rockfill materials using large scale tests." *Electron. J. Geotech. Eng.*, 13(G), 1-16.

Hu, W., Yin, Z. Y., Dano, C., and Hicher, P.-Y. (2011). "A constitutive model for granular materials considering grain breakage." *Sci. China, Ser. E: Technol. Sci.*, 54(8), 2188-2196.

Indraratna, B., Lackenby, J., and Christie, D. (2005). "Effect of confining pressure on the degradation of ballast under cyclic loading." *Géotechnique*, 55(4), 325-328.

Indraratna, B., Wijewardena, L. S. S., and Balasubramaniam, A. S. (1993). "Large-scale triaxial testing of greywacke rockfill." *Géotechnique*, 43(1), 37-51.

Ishihara, K. (1993). "Liquefaction and flow failure during earthquakes." *Géotechnique*, 43(3), 351-415.

Ishihara, K., Tatsuoka, F., and Yasuda, S. (1975). "Undrained deformation and liquefaction of sand under cyclic stresses." *Soils Found.*, 15(1), 29-44.

Iwashita, K., and Oda, M. (2000). "Micro-deformation mechanism of shear banding process based on modified distinct element method." *Powder Technol.*, 109(1-3), 192-205.

Jiang, M. J., Yan, H. B., Zhu, H. H., and Utili, S. (2011). "Modeling shear behavior and strain localization in cemented sands by two-dimensional distinct element method analyses." *Comput. Geotech.*, 38(1), 14-29.

Kikumoto, M., Muir Wood, D., and Russell, A. (2010). "Particle crushing and deformation behaviour." *Soils Found.*, 50(4), 547-563.

Li, X. S., and Dafalias, Y. F. (2000). "Dilatancy for cohesionless soils." *Géotechnique*, 50(4), 449-460.

Li, X. S., and Wang, Y. (1998). "Linear representation of steady-state line for sand." *J. Geotech. Geoenviron. Eng.*, 10.1061/(ASCE)1090-0241(1998)124:12(1215), 1215-1217.

Lobo-Guerrero, S., and Vallejo, L. E. (2006). "Discrete element method analysis of railtrack ballast degradation during cyclic loading." *Granul. Matter*, 8(3-4), 195-204.

Lu, M., and McDowell, G. R. (2010). "Discrete element modelling of railway ballast under monotonic and cyclic triaxial loading." *Géotechnique*, 60(6), 459-467.

Luong, M. P. (1982). "Stress-strain aspects of cohesiveless soils under cyclic and transient loading [C]." *Proc., Int. Symp. on Soils under Cyclic and Transient Loading*, Wiley, New York, 315-324.

Muir Wood, D., and Maeda, K. (2008). "Changing grading of soil: Effect on critical states." *Acta Geotech.*, 3(1), 3-14.

Muir Wood, D., Maeda, K., and Nukudani, E. (2010). "Modelling mechanical consequences of erosion." *Géotechnique*, 60(6), 447-457.

Ng, T.-T. (2006). "Input parameters of discrete element methods." *J. Eng. Mech.*, 10.1061/(ASCE)0733-9399(2006)132:7(723), 723-729.

Ng, T.-T. (2008). "Effect of sample size and shape on critical state." *Proc., GeoCongress 2008: Characterization, Monitoring, and Modeling of GeoSystems*, Geotechnical special publication 179, A. N. Alshawabkeh, K. R. Reddy, and M. V. Khire, eds., ASCE, Reston, VA, 950-957.

O'Sullivan, C., Bray, J. D., and Li, S. (2003). "A new approach for calculating strain for particulate media." *Int. J. Numer. Anal. Methods Geomech.*, 27(10), 859-877.

*Particle Flow Code in 3 Dimensions (PFC3D) 3.10* [Computer software]. Minneapolis, Itasca Consulting Group.

Richart, F. E., Hall, J. R., and Woods, R. D. (1970). *Vibrations of soils and foundations: International series in theoretical and applied mechanics*, Prentice Hall, Upper Saddle River, NJ.

Verdugo, R., and de la Hoz, K. (2007). "Strength and stiffness of coarse granular soils." *Proc., Geotechnical Symp.: Soil Stress-Strain Behavior: Measurement, Modeling and Analysis*, Springer, Rotterdam, Netherlands, 243-252.

Vermeer, P. (1978). "A double hardening model for sand." *Géotechnique*, 28(4), 413-433.
Yan, W. M., and Dong, J. (2011). "Effect of particle grading on the response of an idealized granular assemblage." *Int. J. Geomech.*, 10.1061/(ASCE) GM.1943-5622.0000085, 276-285.

Yang, Z. X., Jardine, R. J., Zhu, B. T., Foray, P., and Tsuha, C. H. C. (2010). "Sand grain crushing and interface shearing during displacement pile installation in sand." *Géotechnique*, 60(6), 469-482.

Yasufuku, N., and Hyde, A. F. L. (1995). "Pile end-bearing capacity in crushable sands." *Géotechnique*, 45(4), 663-676.

Yimsiri, S., and Soga, K. (2011). "Effects of soil fabric on behaviors of granular soils: Microscopic modeling." *Comput. Geotech.*, 38(7), 861-874.

Yin, Z.-Y., and Chang, C. S. (2011). "Stress-dilatancy behavior for sand under loading and unloading conditions." *Int. J. Numer. Anal. Methods Geomech.*, 37(8), 855-870.

Yin, Z.-Y., Chang, C. S., and Hicher, P.-Y. (2010). "Micromechanical modelling for effect of inherent anisotropy on cyclic behaviour of sand." *Int. J. Solids Struct.*, 47(14-15), 1933-1951.