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

AbstractThe 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 Cu. An exponential relationship between the critical-state parameters and the coefficient of uniformity Cu 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 Cu, for both glass beads and natural sand, extending previous literature results for high Cu values up to 20. Mo...


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 p9 plane according to the evolution of GSD possibly described by the increase of the coefficient of uniformity. Based on the numerical analysis of twodimensional (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  proposed a continuum model incorporating a relationship between a grading index I g and the CSL in the e-p9 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 5 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 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.

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): where a 5 dimensionless variable; and d M 5 maximum grain size.
Changing the value of a produces different GSDs. The GSDs of the specimens used for the analysis are shown in Fig. 1. The relationship between C u and a is the following: 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 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 23 ). 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.

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 m ini . The system was iterated to equilibrium. The initial interparticle friction coefficient m ini was carefully adjusted by trial and error to obtain the specified void ratio e 0 5 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 m 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 m ini to m 5 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, m ini 5 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 Du=ðr 1 1 r 2 Þ in the range of approximately 0.0020-0.0039 (Du 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.

Number of Particles
Three different sample sizes with gradings C u 5 1:0 and C u 5 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 5 1:0, and 28,196 particles for the grading C u 5 6:0. The oscillations in stresses are smaller for an increasing number of particles as showed in Figs. 3(c and d).

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 5 1:0 and C u 5 6:0) were generated with three different random seeds (g1 5 10,000; g2 5 50,000; (g) ɛ 1 5 24%; (h) ɛ 1 5 28%; (i) ɛ 1 5 32%; (j) ɛ 1 5 34% and g3 5 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 5 1:0 or C u 5 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 5 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 5 1:8, 3:6, and 6:0). These results have to be analyzed in the framework of the state parameter concept c 5 e 2 e c (Been and Jefferies 1985), where c is the state parameter, e is the current void ratio, and e c is the critical-state void ratio. Considering the distribution characterized by C u 5 1:57 as an example [Fig. 6(c)], (p 0 9, e 0 ) is the typical initial state point, and (p c 9, 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 5 1:57, the initial state parameter c 0 5 e 0 2 e c is negative, and during loading, the void ratio remains smaller than the e c until the critical state is reached (keeping c , 0). Therefore, the specimen exhibited dilatancy upon shearing. For the same reason, the specimen with C u 5 1:0 is also dilative. On the contrary, the specimen with C u 5 3:6 shows contraction because c 0 5 e 0 2 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-p9 plane and e-p9 plane, it was possible to define the CSLs as presented in Fig. 7. All simulations led to the conclusion that the critical stress ratio h 5 q=p9 5 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 5 0:737 (equivalent to a friction angle of 19.2°) 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-p9 plane, the location of the CSL on the e-p9 plane depends on the GSD of the material. As demonstrated in Fig. 7(b), the CSL shifts downward on the compression e-p9 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-p9 plane: where e c 5 critical-state void ratio at the mean effective stress p9; reference void ratio e ref allocates the position of the CSL; l controls the slope of the CSL; and j 5 material constant. Thus, the CSL for a given material can be determined by the two parameters, e ref and l, because p at and j are assigned; in this study, p at 5 101:3 kPa, and j 5 0:9. This value of j, calibrated from the experimental results, is similar to the value proposed by Yan and Dong (2011). Based on the CSL family, the parameters e ref and l are plotted versus the GSD index C u , respectively (Fig. 8). A similar exponential relation with C u is derived for both e ref and l as follows: where a e,ref , b e,ref , c e,ref , a l , b l , and c l 5 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   (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 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 p9-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 5 q=p9) upon compression regardless of GSD (M 5 0:81 for glass beads and M 5 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-p9 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 Yan and Dong 2011). Similar relationships as those derived from the DEM results for both e ref and l 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   the e-p9 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 boundaryvalue 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 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 where dɛ e v 5 elastic volume strain increment; dɛ e d 5 elastic deviatoric strain increment; and G and K 5 hypoelastic shear modulus and bulk modulus, respectively, defined as follows (Richart et al. 1970): where G 0 , K 0 , and n 5 elastic parameters; and p at 5 atmospheric pressure used as reference pressure (p at 5 101:3 kPa). For the plastic behavior, the proposed model has two yield surfaces, 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 p9-q plane, written as follows: where h 5 q=p9; and H 5 hardening variable defined by a hyperbolic function in the H-ɛ p d plane, similar to the following one proposed by Yin et al. (2010): where ɛ p d 5 deviatoric plastic strain; and G p is used to control the initial slope of the hyperbolic curve h-ɛ p d . 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: where D 5 material parameter; and M pt 5 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 p9 obtained by Eq. (3). According to Biarez and Hicher (1997), the peak friction angle f p [related to M p 5 6 sin f p =ð3 2 sin f p Þ for triaxial compression test] is linked to the intrinsic friction angle f m [related to the critical-state value M 5 6 sin f m =ð3 2 sin f m Þ for triaxial compression test] and the density state e c =e e tan f p ¼ e c tan f m 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 f p is greater than f m 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  (d-f) deviatoric stress versus axial strain; (g-i) void ratio versus mean effective stress a consequence, the peak friction angle is reduced, which results in a strain-softening phenomenon. The M pt [related to 6 sin f pt =ð3 2 sin f pt Þ] is the slope of the phase transformation line, which is assumed to be a function of the void ratio 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., ɛ a , 0:1%) of drained triaxial compression tests (in this study, the authors assumed a Poisson's ratio n 5 0:2 to calculate G 0 from K 0 ); • G p was obtained by fitting the initial slope of the curve q-ɛ a (for ɛ a , 1%) in drained triaxial tests; • f m was determined from drained triaxial tests; and • a e,ref , b e,ref , c e,ref , a l , b l , and c l 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. 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  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. 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. 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-p9 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 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 5 1:1; (c and d) C u 5 1:4; (e and f) C u 5 2:5; (g and h) C u 5 5; (i and j) C u 5 10; (k and l) C u 5 20 C u 5 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 , 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 where v cs 5 specific critical volume; v ref 5 reference position of the CSLs in the v-p9 plane when p9 5 0 as defined in Eq. (14); I g 5 grading index; v 0 5 value of the minimum void ratio for a monodisperse material I g 5 0; v g v 0 2 1; Dv 5 material parameter; p cs 5 reference stress; and b controls the shape of the CSL. Because Dv is constant, v ref has a linear relationship with the grading index I g .
Beside the critical-state parameters, the other parameters of Model-I g were determined from simulating a drained triaxial test with C u 5 1:1 or I g 5 0:0043 sheared under a confining stress of 400 kPa: a 5 0:0016, M 5 0:81, A 5 0:9, k R 5 2, and k D 5 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 ref -I g relationship implemented in the Model-I g and the exponential, so nonlinear, e 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.   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-p9 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-p9 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 criticalstate 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.