Modeling Mechanical Behavior of Very Coarse Granular Materials

A novel approach has been developed to predict the mechanical behavior of very coarse granular materials with a constitutive model, which considers both grain breakage and size effect. The behavior of granular assemblies is significantly influenced by particle breakage. A critical-state double-yield-surface model incorporating the change in the critical-state line and in elastic stiffness caused by grain breakage during loading has been adopted. The amount of grain breakage was estimated by extending the size effect theory on individual grains to granular assemblies. The results from earlier studies on granular materials with parallel gradations have been usefully exploited to calibrate and to validate the model. Comparisons between experiments and simulations suggest that this approach can predict the mechanical behavior of very coarse granular materials from test results performed on a finer fraction with a homothetic gradation. keywords: Coarse granular material; Critical state; Grain breakage; Plastic work; Weibull’s theory; Elastoplastic model.


Introduction
The construction of large civil engineering works involving coarse grained materials, such as rockfill dams, is increasing worldwide. However, despite the steady efforts of geotechnical engineers and scientists to study the mechanical behavior of coarse granular materials, the design methods of such geotechnical structures need still to be significantly improved to avoid severe accidents. Accidents occur because it is difficult to construct apparatus large enough for testing such materials, even for small granular materials, and the cost can be prohibitive. For example, for a 0-250-mm rockfill, a representative cylindrical probe for a triaxial test should measure approximately 1.5 m in diameter by 3 m high, and weigh more than 10 t.
Pioneering experimental work in large-scale testing of a wide range of different materials was developed in the 1960s, together with the development of large rockfill dam construction (Marsal 1967;Marachi et al. 1969). Smaller devices have also been used to provide useful results concerning the behavior of coarse granular materials (Chavez and Alonso 2003;Frossard et al. 2012). In these tests, grain breakage was generally observed and shown to depend on the grain sizes. As a consequence, the size effect of granular assemblies was observed, but its quantification remains difficult to assess.
In this paper, the authors suggest a novel modeling approach to evaluate the mechanical properties of assemblies of large particles from the mechanical behavior of assemblies of smaller particles on the basis of size effect in granular materials. For this purpose, a critical-state-based elastoplastic model is adopted. The influence of the particle size on the mechanical behavior of granular assemblies is related to the amount of particle breakage during loading and the size influence on particle crushing strength. The effect of grain breakage is incorporated into the model by introducing the evolution of the grain-size distribution in the constitutive equations. The model has been finally validated by simulating the previously done triaxial tests on Pyramid dam rockfill.

Basic Elastoplastic Model for Granular Material
For building an efficient constitutive model that is able to reproduce the main features of coarse granular materials, one needs at first to consider the main aspects of behavior of a given granular material. For this purpose, the developed model, similar to that by Yin et al. (2013a), contains two yield surfaces: one for shear sliding (f 1 ) and one for compression (f 2 ). The plastic strain increment can be expressed as follows: where the scripts 1 and 2 = shear sliding and compression components, respectively. For one component with f < 0, the corresponding dλ is equal to zero.

Nonlinear Elasticity
The nonlinear hypoelastic model by Richart et al. (1970) was selected. The bulk modulus K and the shear modulus G are power functions of the mean effective stress p 0 and are dependent on the material void ratio e: ð2.97 − eÞ 2 ð1 þ eÞ p 0 p at n ; where G 0 , K 0 , and n = material constants; and p at = atmospheric pressure used as reference pressure (p at ¼ 101.3 kPa). A constant Poison's ratio ν ¼ ð3K 0 − 2G 0 Þ=2ð3K 0 þ G 0 Þ was implied by Eq.
(2). Thus, the input parameters can be two of G 0 , K 0 , and ν with n.

Plastic-Sliding Behavior
The shear sliding yield function, f 1 , has the following expression: where the deviatoric plastic strain ε p d = hardening variable; q = deviatoric stress; G pR = relative plastic stiffness controlling the hardening rate, which is similar to micromechanical models of Yin et al. (2009Yin et al. ( , 2010Yin et al. ( , 2011Yin et al. ( , 2013bYin et al. ( , 2014; and M p = peak stress ratio, which can be expressed as where ϕ p = "peak" friction angle with more details described subsequently.
The plastic potential function is given by where g 1 is different from f 1 ; therefore, the model is nonassociated. M pt can be expressed as where ϕ pt = phase transformation angle or characteristic angle, which corresponds to the transition from the contracting to the dilating domain, with more details described subsequently.

Critical-State Line-Related Density State Effect
Critical state is one of the most important concepts in soil mechanics. At critical state, the material will remain constant in both volume and stress while subjected to a continuous distortion. The void ratio corresponding to this state is e c , named as the critical void ratio, which is a function of the mean effective stress p 0 as follows: The friction angle at critical state ϕ u is considered constant for the material. However, the peak friction angle ϕ p is dependent on the stress state and the void ratio where β = material constant. A common value is β ¼ 1, according to Biarez and Hicher (1994). This formulation indicates that, in a loose structure, in which particles have more freedom to move, the interparticle shear force can be fully mobilized during sliding. Under these conditions, the peak friction angle ϕ p is equal to ϕ u . Conversely, in dense packing, in which there is a higher degree of interlocking, more effort is needed to mobilize the particles in contact. In this case, the peak friction angle ϕ p is greater than ϕ u . When the stress state passes the phase transformation state, the granular assembly begins to dilate. The void ratio increases. As a consequence, the peak friction angle is reduced, which results in a strain-softening phenomenon. The phase transformation angle ϕ pt depends on the packing density and can be simply expressed as follows by adopting the same form as the peak strength: which indicates that a dense packing has a smaller phase transformation angle than does a loose packing. This equation as shown in Yin et al. (2010Yin et al. ( , 2014 gives the same trend as experiments and exponential forms used in Jefferies and Been (2006) and Kikumoto et al. (2010).

Compression Behavior
To describe the compression behavior of granular materials, a second yield surface was added. The function of the breakage limit by Kikumoto et al. (2010) was adopted as the compression yield surface in this study where p m = hardening variable controlling the size of the yield surface. This yield surface expands with the plastic volumetric strain, as in the modified Cam-clay model where κ ¼ ð1 þ e 0 Þp 0 =K is used with K by Eq. (2); and λ takes the same value as that of the CSL.

Considering Grain Breakage
Estimating the Amount of Grain Breakage Grain breakage commonly occurs when a granular material undergoes compression and shearing. The amount of grain breakage caused by mechanical loading is usually evaluated on the basis of the change of the grain-size distribution (GDS). Different measures have been proposed in previous studies (Lee and Farhoomand 1967;Marsal 1967;Hardin 1985;Lade et al. 1996;Nakata et al. 1999;Einav 2007;Kikumoto et al. 2010). Because the GSD tends to become fractal as shown by McDowell et al. (1996) and Coop et al. (2004), the breakage index B Ã r suggested by Einav (2007) varies from 0 to 1 with the change in the GSD. Therefore, B Ã r should be a pertinent measure of the amount of grain breakage during mechanical loading.
For a given material, the amount of grain breakage increases when stresses and strains increase (Biarez and Hicher 1997). Therefore, the amount of grain breakage can be considered as a function of the plastic work during loading (Daouadji et al. 2001 where a = material constant controlling the evolution of the breakage amount. In this study, the expression of the plastic work is slightly modified where dε p v and dε p d = volumetric and deviatoric plastic strain increments, respectively. The MacCauley function in Eq. (13) implies that the shear induced dilation (dε p v < 0) is not accounted for in the modified plastic work. As a result, the evolution of the gradation is not influenced by shear-induced dilation (Biarez and Hicher 1997).
Once the amount of grain breakage is determined, the actual GSD can be obtained by the formulation proposed by Einav (2007) where F 0 ðdÞ and F u ðdÞ = initial and fractal GSDs, respectively, with F u ðdÞ ¼ ðd=d M Þ 0.3 according to Coop et al. (2004).

Dynamic Critical-State Line
Recent experimental and numerical works by Liu et al. (2013) and Li et al. (2014) showed that the position of the critical-state line in the elog p 0 plane depends on the grain-size distribution represented by the coefficient of uniformity C u ¼ d 60 =d 10 . If grain breakage occurs during loading, the grain-size distribution evolves, and the value of d 60 =d 10 increases.  Biarez and Hicher (1997). A simple way to consider the evolution of e ref or p ref is, therefore, to make it dependent of the amount of grain breakage.
In this study, p ref corresponding to a fixed value of e ref ¼ 0.5 was adopted to describe the evolution of the critical-state line caused by grain breakage; the following expression is suggested: where parameter b controls the rate of the CSL shifting caused by grain breakage. Because there exists an ultimate CSL corresponding to the fractal GSD, the value of b ¼ lnðp ref0 =p refu Þ can be obtained by the ultimate p refu when B r ¼ 1, if tests for the material with fractal GSD are available. From Eq. (15) combined with Eqs. (7)-(9), ruptures of particles will reduce the dilatancy of the granular assembly and increase its contractancy. As a consequence, the material will develop more deformation when subjected to oedometric compression and the maximum strength along deviatoric loading paths will reduce.

Influence of Grain Breakage on Compression Behavior
The evolution of p m caused by grain breakage was added in Eq. (11) as follows: Because the same evolution is used for p m and p ref , the same amount of shifting is imposed for the isotropic compression curve and the CSL.
For simplicity, an associated flow rule has been adopted for the normal compression. Therefore, only the initial value of p m0 for ε p v ¼ 0 and B Ã r ¼ 0 is needed for input.

Influence of GSD on Elastic Stiffness
According to Iwasaki and Tatsuoka (1977), the shear modulus of natural sands at small strain decreases with the broadening of the grain-size distribution. Because the change of GSD is estimated by the amount of grain breakage, the elastic stiffness can be modified as follows: where the parameter ρ controls the degradation rate and also the final degradation percentage of elastic stiffness when B r ¼ 1.
For instance, on the basis of the results of Iruma sands shown by Iwasaki and Tatsuoka (1977), ρ ¼ 0.94 was estimated for a final degradation of up to 60%. Because the plastic stiffness for η-ε p d in Eq. (3) is controlled by G pR G, the aforementioned modification will also lead to the degradation to the plastic stiffness, in agreement with the experimental observations.

Size Effect in Individual Grains
The statistical theory of the strength of brittle materials proposed by Weibull (1939) gives the following distribution for the survival probability of a material of size d subjected to a tensile stress σ where σ o = characteristic strength (P s ¼ 37% for a sample of size d o ); and m controls the amplitude of data scatter on crushing strength. According to Bažant and Planas (1997), n d is the geometric similarity of the mechanical problem: n d ¼ 1, 2 or 3, for linear, surface or volume similarity, respectively. For a given survival probability and known empirical parameters, a size-effect relation is obtained between the induced tensile stress at failure (σ f ) and the characteristic size of the sample (d) Marsal (1973) suggested the following empirical expression between the crushing force (f f ) and the rock particle size for a crushing compression test between two stiff parallel platens: where d = characteristic size of the particle (distance between platens before crushing); and η and ζ = experimentally fitted parameters. f f induces a tensile stress σ f , which causes the brittle failure of the particle and can be expressed as (Jaeger 1967) σ From Eqs. (20) and (21), a power law expressing the size effect on particle crushing strength can be obtained (where a ¼ ζ − 2 is an empirical parameter) Empirical validation of Weibull's theory is possible by confirming the validity of the relation α ¼ −n d =m [from Eqs. (18) and (22)]. For soil particles and rock aggregates under compression between two stiff parallel platens, it has been generally assumed that failure occurs under bulk-induced tensile stress (i.e., in Mode I), which gives n d ¼ 3 (Lobo-Guerrero and Vallejo 2006). Values of m for different materials have been reported in the literature (Ovalle et al. 2014), i.e., approximately 3-4, corresponding to sound heterogeneous materials as quartz sands, biotite and ballast; as for feldspars, limestones, and carbonates, m varies between 1 and 3.
McDowell and Bolton (1998) examined the application of Weibull's theory on weak limestone sand by checking the ability of Eq. (18) for n d ¼ 3, obtaining a good agreement. Similar conclusions for silica sand and for rock fragments of biotite gneiss and grey quartzite can also be found. Conversely, geometric similarities of n d ¼ 2 was also obtained to represent size effects on granite and railway ballast particles, respectively, in which failure was generally induced on weak surface flaws. In general, if no information is available on material flaws distribution or on internal stress distribution on particles, the physical sense of n d might be missed. As aforementioned, Weibull's theory assumed that the structure is equivalent to a uniaxial stressed bar. Afterward, if one considers an alternative structure, all information about the mode of failure is lost and becomes irrelevant. This could be the case of soil grains and rock aggregates with significant data scatter in parameters such as particle size, particle shape, grain anisotropy, and flaws distribution. After several crushing tests, failure could occur in different geometric similarities for different aggregates. Therefore, for practical applications, n d could just be fitted to use Weibull's distribution as a statistical tool (Ovalle et al. 2013).

Size Effect in Granular Assemblies
The method proposed by Frossard et al. (2012) considers two homothetic granular materials (say, G1 as the finer and G2 as the coarser) under the same loading condition and with the same mineralogy and grain shape. Thus, if an identical survival probability is assumed for two grains under compression of characteristic sizes d 1 and d 2 in G1 and G2, respectively, the following relation for the particle crushing strength (σ Gi ) can be obtained from Eq. (18): Then, according to Eq. (21), the crushing forces of each particle (f Gi ) are related by Now, if both granular materials present the same amount of grain breakage after shearing, the stress magnitude on G2 should be lower because coarse grains have lower strength. By focusing on contact micromechanics, a condition for the intensity of the contact forces between grains n and p (f ðn=pÞ ) is obtained In parallel, the geometric scaling gives a relation between branch vectors (l ðn=pÞ ) and grains volumes The stress tensors σ ij for G2 and G1 are given by the following expressions (Rothenburg and Selvadurai 1981): Hence, combining Eqs. (26) and (27), the stress tensor of G2 as a function of G1 is obtained Therefore, the modified plastic works w pG1 and w pG2 necessary to create the same amount of grain breakage in materials G1 and G2 along the same loading path are linked by the relation And, as a consequence, the material constants a G1 and a G2 are linked by a similar relation Therefore, if one can assume that all the other parameters are identical for materials G1 and G2, it is possible to calibrate the parameters of coarse material G2 from experimental testing performed on fine material G1. This assumption is reasonable if materials G1 and G2 are made of the same grains with similar shapes for different sizes and homothetic grain-size distributions. Then, to model the behavior of a given coarse granular material, it is possible to select a finer homothetic grain-size distribution to prepare a representative finer material whose grain sizes are compatible with the available experimental means of the laboratory.

Summary of the Model Parameters
The parameters of the model can be classified into four groups: (1) elastic parameters K 0 , G 0 and n, shown in Eq.
On the basis of all aforementioned constitutive equations and according to conventional elasto-plasticity and the framework of double-yield-surface model (Yin et al. 2013a), the stress-strain relationship can be solved for test simulations.
To predict the mechanical behavior of materials made of large particles, it is possible to determine all the aforementioned parameters from tests on samples made of small particles. Parameters of Groups 1 and 2 can be determined in a conventional way [see details in Yin et al. (2013b)] or by an optimization procedure    Fig. 4. Comparison between experimental and predicted results of drained triaxial tests on rockfill material with a maximum grain size of 4.9 cm (Yin and Hicher 2008). The size effect-related parameters can be determined separately from crushing tests on individual particles of different sizes. Thus, only the value of d M is needed for modeling large assemblies.

Experimental Validation
Selected Material and Experimental Testing Marachi et al. (1969) performed drained triaxial tests on samples of Pyramid dam granular material with parallel gradations and respective maximum sample sizes of 7.1, 30.5, and 91.4 cm (with maximum grain size d M of 1.18, 4.9, and 14.5 cm). All specimens had the same initial compactness. Relatively high values between 200 and 4,500 kPa of the confining pressure were applied. The material was selected because of the homogeneity of the grains over the different granular fractions in terms of mineralogy and grain angularity, resulting from the crushing of a sedimentary rock. Fig. 1 presents the grain-size distributions of the tested materials and the grain-size distribution of the material used to build the dam.

Determination of Model Parameters
The determination of the model parameters was based on four triaxial tests, together with the isotropic compression stage and measurements of the GSD before and after testing performed on the finest material (d max ¼ 1.18 cm), shown in Figs. 2 and 3.
The elastic parameters were determined from the isotropic compression curve [IC stage of tests in Fig. 2(a)]. One obtained K 0 ¼ 70 and n ¼ 0.4. The shear modulus G 0 was calculated by assuming a commonly adopted constant Poisson's ratio ν ¼ 0.25.
The plastic parameters were determined from the triaxial test results. The relative plastic stiffness G pR ¼ 16 was determined from four ε a -q curves of an axial strain level of up to 2% [ Fig. 3(a)].
(d) The critical friction angle ϕ μ ¼ 40 was calculated from the slope of the CSL in the p 0 -q plane measured directly from the final stress states [ Fig. 3(b)]. Because only the triaxial test performed at the smallest confining stress did not produce any significant grain breakage, the slope of the CSL in the eln p 0 plane had to be assumed: λ ¼ 0.04 was reasonably selected corresponding to different granular materials (Biarez and Hicher 1994;Jefferies and Bean 2006). For the grain breakage related parameters, the breakage rate parameter a ¼ 3 MPa was determined on the basis of the evolution of GSD [Figs. 3(e and f)]. The yield stress of the breakage limit P m0 ¼ 1.5 MPa and the degradation of the yield stress b ¼ 4 were determined from the isotropic compression curve up to a high stress level [ Fig. 2(b)]. Because no measurement at small strain is available for the selected material, the parameter of the elastic stiffness degradation ρ ¼ 0.94 was assumed on the basis of the test results on Iruma sands shown by Iwasaki and Tatsuoka (1977).
Because no statistical data on crushing tests on grains were available, no direct determination of the Weibull parameters could be made. The values of the size effect-related parameters, n d ¼ 3 and m ¼ 4, were assumed according to Lobo-Guerrero and Vallejo (2006) and Frossard et al. (2012).
All determined parameters are summarized in Table 1. The triaxial tests on coarser material were simulated by using the same set of parameters for the validation of the model.

Test Simulations
In Fig. 3(a-d), the simulations without grain breakage effect by taking n d ¼ 0 were plotted together with those with grain breakage (n d ¼ 3). The model appears able to reproduce the main features of the mechanical behavior of granular material influenced by particle breakage: (1) under the lowest confining stress (200 kPa), the material shows a dilative behavior; (2) for higher confining stresses, the material becomes contractive, and the disappearance of the dilation is caused partly by the increase of particle breakage occurring under high stresses (from 207 to 4,480 kPa); (3) the grain-size distribution broadens more significantly for higher stress levels [ Fig. 3(e versus f)].
Figs. 4 and 5 show the comparisons between experimental results and model predictions for drained triaxial tests in compression with confining stresses varying from 207 to 4,482 kPa on samples with a maximum grain size of 4.9 and 14.5 cm., taking into account the size effect by setting n d ¼ 3. A relatively good agreement was achieved for all comparisons. The comparison shows that the grain breakage effect on the mechanical behavior is well predicted by the model for samples with larger grain sizes. In addition to this, the simulations without considering the size effect by giving n d ¼ 0 were also plotted in Figs. 4 and 5. The comparisons demonstrate that the behavior of coarse granular materials can be reproduced by the model if the size effect is correctly estimated.

Discussions
The influence of the grain size is illustrated in Fig. 6, in which the numerical simulations of the triaxial tests at the confining pressures of 207 to 4,480 kPa are presented for the four different materials including the field gradation with the maximum grain size varying from 1.18 to 37.2 cm according to Fig. 1. The simulations for the maximum grain size of 37.2 cm corresponding to the GSD of dam rockfills were also carried out and plotted in dashed line. The results demonstrate that at higher stress levels the amount of grain breakage is greater, and, as a consequence, the size effect is more pronounced on both the shear stiffness and the maximum strength. These results agree well with the experimental observations on coarse granular materials. Fig. 7 illustrates the evolution of the GSDs for the different material grain sizes. One can see that this evolution is more pronounced at a given confining stress for coarser materials. The computed GDSs agree well with the experimental ones measured at the end of each triaxial test. The dashed lines show again the computed evolution for the field gradation.

Conclusions
A novel approach for modeling the mechanical behavior of very coarse granular materials has been proposed. First of all, the authors constructed a constitutive model able to reproduce the main features of granular materials. This model incorporates the critical state and the elastic stiffness which, in our approach, are both considered dependent on the amount of grain breakage developed within the granular assembly during loading. The amount of grain breakage is made a function of the plastic work developed along a given loading path. The constitutive relation includes two model parameters, which depend on the individual strength of the grains. A size effect is then introduced, which takes into account the decrease of the grain strength with the grain size. To validate the approach, drained triaxial tests under different levels of confining stresses on rockfill material of Pyramid dam with parallel grain-size distributions were simulated. The same set of parameters was used for all simulations. Comparisons between experiments and simulations demonstrate that the model can predict the mechanical behavior of very coarse granular materials from test results performed on a finer fraction with a homothetic gradation.
From the results of this study, one can suggest a procedure for modeling very coarse granular materials for which experimental behavior is not accessible because of the lack of adapted testing. Standard tests can be performed on a specimen with a smaller grain size and a grain-size distribution homothetic to the one corresponding to the real material. The grain-size distribution after the tests should be measured to analyze the amount of grain breakage. Moreover, crushing tests on individual grains of different sizes should be performed to obtain a quantification of the size effect. By combining these results, the calibration of the model parameters can be made, and numerical simulations can be performed to obtain the stress-strain behavior of the coarse material.