# A constitutive model for granular materials considering grain breakage

Wei Hu, Zhenyu Yin, Christophe Dano, Pierre-Yves Hicher

► To cite this version:

Wei Hu, Zhenyu Yin, Christophe Dano, Pierre-Yves Hicher. A constitutive model for granular materials considering grain breakage. SCIENCE CHINA Technological Sciences, 2011, 54 (8), pp.2188-2196. 10.1007/s11431-011-4491-0. hal-01004692

HAL Id: hal-01004692

https://hal.science/hal-01004692

Submitted on 12 Feb 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>
# A constitutive model for granular materials considering grain breakage

HU Wei¹, YIN ZhenYu², DANO Christophe³ & HICHER Pierre-Yves³

¹ State Key Laboratory of Geohazard Prevention and Geoenvironment Protection, Chengdu University of Technology, Chengdu 610059, China;

² Center for Marine Geotechnics Research, Department of Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China;

³ Research Institute in Civil and Mechanical Engineering, UMR CNRS 6183, Ecole Centrale de Nantes, 44321 Nantes, France

A particle breakage has a significant impact on the mechanical behavior of granular materials. In this paper, we present an elasto-plastic model with two yield surfaces to which the influence of particle breakage has been introduced. The main feature of this model is to incorporate the change in the critical state line (CSL) consequent to the grain breakage induced by isotropic and deviatoric stresses during loading. For this purpose we propose a breakage function which connects the evolution of the CSL to the energy consumed. Results from earlier studies on drained and undrained compression and extension triaxial tests were used to calibrate and validate the model. Comparison between earlier results and our simulations indicates that the model can reproduce with good accuracy the mechanical behavior of crushable granular materials and predict the evolution of the grain size distribution during loading.

granular materials, grain breakage, elasto-plasticity, critical state, grain size distribution

## 1 Introduction

A grain breakage commonly occurs when a granular material undergoes compression and shearing, especially under high confining stress (e.g., earth dams, deep well shafts). Its impact on the mechanical behavior of granular materials has been widely studied in the past decades [1–13], and different modeling methods have been developed. Daouadji et al. [14] connected the position of the critical state line (CSL) to the amount of energy needed for grain breakage, showing that the CSL in the e-logp' (void ratio versus mean effective stress in log scale) descends according to the evolution of gradation. Muir Wood et al. [12] confirmed the change of

the position of CSL with grain gradation. Based on this result, they connected the CSL with a grading state index $I_G$ which is a state parameter that evaluates the evolution of the gradation as a result of grain breakage. Russell et al. [15] used a three-segment type CSL within a boundary surface constitutive model to describe the behavior of crushable granular materials. Salim et al. [16] formulated a ratio between the deviatoric and the mean stresses as a function of dilation and derived a new plastic flow rule from this formulation to take into account the effects of particle breakage. Sun et al. [17] and Yao et al. [18] modified the plastic hardening parameter of their models in order to take into account the effect of particle breakage.

The change of the CSL with gradation can also be found when fines are added to sands (Thevanayagam et al. [19]), which in turn can be evidence of the evolution of the CSL
with gradation changes due to particle breakage. Therefore, the models by Daouadji et al. [14], Muir Wood et al. [20], Daouadji et al. [21] can be better justified from a physical point of view. The model proposed in this paper situates itself along this line. However, departing from the models of Daouadji et al. [14, 21], we propose a simple two-yield surface plastic model using the evolution of the CSL with gradation, and departing from the model of Muir Wood et al. [20], we connect the evolution of the gradation to the plastic work during loading. Since the gradation is an important factor in the proposed model, we have also made it possible to predict its evolution at each stage of loading.

In the first part of the paper, we present an analysis of the connection between the CSL, particle breakage and energy consumption based on experimental results. Then, we formulate an elasto-plastic model within the framework of critical state soil mechanics under triaxial condition. Finally, results of numerical simulations of triaxial tests performed on Cambria Sand under different loading conditions at high confining stresses are compared to experimental data.

# 2 Analysis of breakage

## 2.1 Definition of breakage index

Hardin [2] suggested a breakage index $B_r$ in order to quantify the amount of particle breakage. The index is based on the changes in particle size. Einav [11] modified the definition of this factor, taking into account the changes on the overall grain size distribution and assuming a fractal rule for particle breakage [13].

$$B_r = \frac{B_p}{B_t} = \frac{\int_{d_m}^{d_M} (F(d) - F_0(d)) d(\log(d))}{\int_{d_m}^{d_M} (F_u(d) - F_0(d)) d(\log(d))}, \quad (1)$$

where $B_p$ is the area between the original and the present grain size distributions; $B_t$ is the total area between the original and the ultimate fractal grain size distributions; $F_0(d)$ and $F_u(d)$ represent respectively the initial gradation before grain breakage and the ultimate fractal distribution; $F(d)$ is the present gradation; $d$ is the present grain size; $d_M$ and $d_m$ are the maximum and minimum grain sizes of the material. The present gradation can be expressed as

$$F(d) = (d / d_M)^{\alpha}, \quad (2)$$

where $\alpha$ is a material constant. For $F_0(d)$ the value of $\alpha$ can be measured from the initial grain size distribution, for example $\alpha = 5.5$ for Cambria Sand. For $F_u(d)$, the value of $\alpha = 0.4$, proposed by Coop et al. [3], is adopted in order to obtain the ultimate fractal grain size distribution. Thus, for a given grain size distribution, $\alpha$ can be obtained by fitting the grading curve with eq. (2). Then, the breakage index $B_r$

can be obtained by eq. (1). In turn, for a given $B_r$, the present grain size distribution can be determined.

## 2.2 Influence of plastic work on evolution of gradation

In this section, results from drained triaxial compression tests on Cambria sand [5, 6] are analyzed to investigate the relation between the breakage index $B_r$ and the modified plastic work $w_p$ representing energy consumption. Cambria sand is a coarse and uniform sand consisting of sub-angular to well-rounded grains. The grain diameter varies from 0.83 to 2 mm. The maximum and minimum void ratios are 0.792 and 0.503, respectively. The specific gravity is 2.69.

The modified plastic work is expressed as follows

$$w_p = \int p' \langle d\varepsilon_v^p \rangle + q d\varepsilon_d^p, \quad (3)$$

where $p'$ is the mean effective stress: $p' = (\sigma'_1 + 2\sigma'_3) / 3$; $q$ is the deviatoric stress: $q = \sigma'_1 - \sigma'_3$; $d\varepsilon_v^p$ and $d\varepsilon_d^p$ are volumetric and deviatoric plastic strain increments, respectively: $d\varepsilon_v^p = d\varepsilon_1^p + 2d\varepsilon_3^p$ and $d\varepsilon_d^p = 2(d\varepsilon_1^p - d\varepsilon_3^p) / 3$; $\langle F \rangle$ is the MacCauley function: $\langle F \rangle = 0$ for $F < 0$ and $\langle F \rangle = F$ for $F > 0$. The MacCauley function in eq. (3) implies that the shear induced dilation ($d\varepsilon_v^p < 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 based on drained triaxial tests with confining stresses less than 2.1 MPa.

The plastic strain increments were calculated from the total strain increments by subtracting the elastic strain increments, using the following elastic law:

$$d\varepsilon_v^e = \frac{dp}{K}, \quad d\varepsilon_d^e = \frac{dq}{3G}, \quad (4)$$

where $G$ and $K$ are the hypo-elastic shear and bulk modulus, respectively, defined as follows (Richart et al. [22]):

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

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

where $G_0$, $K_0$ and $n$ are elastic parameters; $p_{at}$ is the atmospheric pressure used as reference pressure ($p_{at} = 101$ kPa). For Cambria sand, $K_0 = 26.3$ MPa and $n = 0.4$ were determined from isotropic compression test, and $G_0 = 35$ MPa from the initial slope of the stress-strain curve (e.g., $\varepsilon_1 < 0.1\%$) of drained triaxial compression tests.

The values of the breakage index $B_r$ as well as the modified plastic work were measured for different tests. $B_r$ is plotted versus $w_p$ in Figure 1(a) which shows that the value of $B_r$ increases with the modified plastic work. Based on
these results, we suggest a hyperbolic function between $B_r$ and $w_p$

$$B_r = \frac{w_p}{\chi + w_p}, \qquad (7)$$

where $\chi$ is a material constant controlling the evolution rate of the gradation. For Cambria sand, $\chi = 15000$ was obtained.

## 2.3 Influence of gradation on the position of the CSL

One of the important elements to be considered in soil modeling is the critical state concept. The critical state void ratio $e_c$ is a function of the mean effective stress $p'$. The relationship is traditionally written in the $e-\log p'$ plane as

$$e_c = e_{ref} - \lambda \ln\left(\frac{p'}{p_{ref}}\right), \qquad (8)$$

where $e_{ref}$ and $p_{ref}$ determine a particular point of the critical-state in the $e-\log p'$ plane, $\lambda$ determines the slope of the CSL. Thus, the CSL can be determined by the two parameters $e_{ref}$ and $\lambda$ when $p_{ref}$ is assigned (in this study, $p_{ref} = 100$ kPa).

According to Biarez et al. [4], the position of the CSL moves down in the $e-\log p'$ plan with the increase of the coefficient of uniformity $C_u = d_{60}/d_{10}$. The position of the CSL as a function of $C_u$ has also been demonstrated by Muir Wood et al. [12] through simulations by the discrete element method. Their results agree with the concept developed by Biarez and Hicher. However, up to now the studies on the relation between the position of the CSL and the gradation are not based on experimental results. The concept of the CSL is based on the assumption that at the critical state the material remains at a constant volume while being subjected to continuous distortion. If ever the CSL is able to move, this concept becomes invalid. This paper extends the concept of critical state by defining the position of the CSL at a given loading stage, corresponding to the present gradation.

The drained compression tests performed by Yamamuro

et al. [5] on Cambria sand were used to investigate the evolution of the CSL with the gradation. For each drained compression test, the void ratio at failure was measured and the state $e, \log p'$ was considered as the critical state corresponding to the gradation at the final stage of the test. The value of $e_{ref}$ representing the position of the CSL was obtained by using eq. (8). Based on all drained compression tests, $e_{ref}$ is plotted versus the breakage index $B_r$, as shown in Figure 1(b), from which a hyperbolic relation can be derived as

$$e_{ref} = e_{ref0} + (e_{refu} - e_{ref0}) \frac{B_r}{\rho + B_r}, \qquad (9)$$

where $e_{ref0}$ and $e_{refu}$ are the initial and ultimate reference critical state void ratios, respectively; $\rho$ is a material constant controlling the evolution rate of the CSL with particle breakage. For the Cambria sand, $e_{ref0} = 0.59$ and $\lambda = 0.006$ were obtained from drained triaxial compression tests under low confining stresses (less than 1 MPa) for which Yamamuro et al. [5] indicated that very limited grain breakage occurred. $e_{refu} = 0.13$ and $\rho = 0.16$ were obtained from Figure 1(b).

# 3 Constitutive model

An isotropic hypo-elasticity was assumed for the elastic part of the model (see eqs. (4) to (6)). Thus, three parameters are required for the elastic behavior: $G_0$, $n$, $\nu$. For the plastic behavior, the proposed approach uses two yield surfaces, one for shear sliding and one for normal compression. Thus, the framework of the proposed approach is similar to that of the double-hardening model developed by Vermeer [23].

## 3.1 Shear sliding criterion

As in many models for sand [20, 23], the shape of the yield surface for the shear component is linear in $p' - q$ plot, written as follows:

$$f_s = \eta - H, \qquad (10)$$

<figure><img src="image_2.png" /><figcaption>**Figure 1** Particle crushing effect. (a) Evolution of breakage index versus modified plastic work, and (b) evolution of reference critical state void ratio versus breakage index.</figcaption></figure>
where $\eta = q/p'$; H is the hardening variable defined by a hyperbolic function in the H-$\varepsilon_d^p$ plane, similar to the one proposed by Yin et al. [24]

$$H = \frac{M_p \varepsilon_d^p}{\frac{M_p p'}{G_p K} + \varepsilon_d^p}, \qquad (11)$$

where $G_p$ is used to control the initial slope of the hyperbol-ic curve $\eta-\varepsilon_d^p$. Eq. (11) guarantees that the stress ratio will approach the peak value of stress ratio $M_p$.

In order to take into account dilation or contraction dur-ing shear sliding, a Roscoe-type stress dilatancy equation is used

$$\frac{d\varepsilon_v^p}{d\varepsilon_d^p} = D(M_{pt} - \eta), \qquad (12)$$

where D is a soil parameter. $M_{pt}$ is the slope of the phase transformation line for sand as defined by Ishihara et al. [25] or the characteristic line defined by Luong [26].

## 3.2 Normal Compression Criterion

In order to describe the compressible behavior of breakable granular materials, we added a second yield surface. The second yield function is assumed to be as follows

$$f_N = p' - p_y, \qquad (13)$$

where $p_y$ is the hardening variable controlling the size of the yield surface. The yield surface expands with the plastic volumetric strain. The hardening rule of the modified Cam Clay model has been adopted:

$$d p_y = p_y \frac{d \varepsilon_v^p}{c_p}. \qquad (14)$$

An associated flow rule has been adopted for the normal compression. The initial value of compression yield stress $p_{y0}$ for $\varepsilon_v^p = 0$ is also needed for the model.

In order to interpolate the slope of critical state line in $p' - q$ plan “M” between its values $M_c$ (for compression) and $M_e$ (for extension) by means of the Lode angle $\theta$ (see Sheng et al. [27]), the following relation is adopted:

$$M = M_c \left[ \frac{2c^4}{1+c^4+(1-c^4)\sin 3\theta} \right]^{\frac{1}{4}}, \qquad (15)$$

where $c=(3-\sin\phi_\mu)/(3+\sin\phi_\mu)$, assuming the same friction angle for compression and extension; $-\frac{\pi}{6} \le \theta = \frac{1}{3}\sin^{-1}\left(\frac{-3\sqrt{3}J_3}{2J_2^{3/2}}\right) \le \frac{\pi}{6}$; $J_2 = (I_1^2 + 2I_3)/3$ and $J_3 = 2I_1^3 + 9I_1I_2 + 27I_3/2$

27 are the second and the third invariants of the deviatoric stress tensor with $I_1 = \sigma'_1 + \sigma'_2 + \sigma'_3$, $I_2 = \sigma'_1\sigma'_2 + \sigma'_2\sigma'_3 + \sigma'_3\sigma'_1$, $I_3 = \sigma'_1\sigma'_2\sigma'_3$.

## 3.3 Density state effect

The material’s density state is defined by the ratio $e_c/e$, where $e$ is the present void ratio and $e_c$ is the critical void ratio for the present value of $p'$ obtained by eq. (8).

According to Biarez et al. [4], the peak friction angle $\phi_p$ (related to $M_p = 6\sin\phi_p/(3-\sin\phi_p)$ for triaxial compression) is linked to the intrinsic friction angle $\phi_\mu$ (related to the critical state value $M=6\sin\phi_\mu/(3-\sin\phi_\mu)$ for triaxial compression) and the density state $e_c/e$

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

$M_p$ is then obtained through the critical state $M$ and the den-sity state $e_c/e$. Eq. (16) means that in a loose assembly the peak friction angle $\phi_p$ is smaller than $\phi_\mu$. On the other hand, a dense state provides a higher degree of interlocking. Therefore, the peak friction angle $\phi_p$ is greater than $\phi_\mu$. When the stress state reaches the phase transformation line, the dense assembly dilates and the degree of interlocking decreases. As a consequence, the peak friction angle is reduced, which results in a strain-softening phenomenon.

$M_{pt}$ is the slope of the phase transformation line for sand which we assume to be a function of the void ratio

$$e_c \tan \phi_{pt} = e \tan \phi_\mu. \qquad (17)$$

Eq. (17) indicates that a dense packing has a smaller phase transformation angle than a loose packing, producing the same effect as in the formulation used by Muir Wood et al. [20].

From eqs. (7), (8) and (9), the particle breakage directly influences the position of CSL, which results in the change of the density state $e_c/e$. All the terms related to $e_c/e$ (e.g., $\phi_p$, $\phi_{pt}$, H, etc.) are then influenced by breakage, which allows us to incorporate the influence of particle breakage into the model.

# 4 Test simulations

## 4.1 Determination of model parameters

The drained triaxial tests on Cambria sand performed by Yamamuro et al. [5] and the undrained triaxial tests on the same material by Bopp et al. [9] were used to calibrate and validate the model. All selected samples were isotropically consolidated with an initial void ratio of 0.52 at different stress levels before shearing. The determination of the mod-el parameters is based on one isotropic compression test and drained triaxial tests in compression as follows.
• $K_0 = 26.3$ MPa and $n = 0.4$ were calibrated from an isotropic compression test ($p' < 12$ MPa, see Figure 2(a)), and $G_0 = 35$ MPa was calibrated from the initial slope of the stress strain curve (e.g., $\epsilon_1 < 0.1$ %) of drained triaxial compression tests.

• $p_{y0} = 12$ MPa and $c_p = 0.028$ were obtained by the curve fitting of the isotropic compression test (see Figure 2(a)). In Figure 2(a), simulations were also carried out for different values of $c_p$ which controls the slope of the compression line under high stress levels.

• $G_p = 3.5$ MPa was obtained by fitting the initial slope of the curve $q-\epsilon_a$ (for $\epsilon_a < 1\%$) of the drained test at the confining stress 26 MPa.

• $\phi_\mu = 37.5^\circ$ was determined from drained triaxial tests at lower confining stresses (2.1, 4, 5.8 MPa). $e_{ref0} = 0.59$ and $\lambda = 0.006$ were obtained from drained triaxial compression tests under low confining stresses (less than 1 MPa) for which Yamamuro et al. [5] indicated very slight grain breakage.

• $\chi = 15000$, $\rho = 0.16$ and $e_{refu} = 0.13$ were determined on the basis of the breakage analysis (see Figures 1(a), (b), 2(b) and (c)). Numerical simulations with different values of $\chi$ agree with the experimental results on carbonate sands presented by Coop [3]: particle breakage increases the contraction and decreases the peak friction angle of the material.

All the determined values of the model parameters are summarized in Table 1, and are used for simulating the tests with different stress paths.

## 4.2 Simulations of drained triaxial tests in compression and extension

Figure 3 shows comparison between experimental results and numerical simulations for drained triaxial tests in compression with confining stresses varying from 2.1 to 52 MPa. A good agreement was achieved for all comparisons. With

the model, we are able to reproduce the main features of the mechanical behavior of sand influenced by particle breakage.

(i) Under the lowest confining stress (2.1 MPa), the material exhibits a dilative behavior.

(ii) For higher confining stresses, the material becomes contractive. The disappearance of the dilation is linked to the increase of particle breakage occurring under high stresses (2.1–26 MPa). For tests under confining stresses from 2.1 to 26 MPa, the volumetric strain increases with the increase of the confining stress.

(iii) For tests under confining stresses from 26 to 52 MPa, the volumetric strain decreases with the increase of the confining stress. Yamamuro et al. [5] indicated that this effect is caused by increasing amounts of volumetric contraction and particle breakage during the isotropic consolidation stage, and that lower void ratios are obtained by increasing the confining stress, which in turn generates less volumetric contraction during shearing. This trend was well reproduced by the model incorporating particle breakage. The increase of particle breakage during isotropic loading results in a smaller amount of breakable grains left over during the shearing stage (see Figure 1(a) for the evolution of $B_r$ which becomes stable for high plastic work). As a result, the material becomes less contractive.

Using the method presented in the section concerning the breakage analysis, we can predict the evolution of the gradation. The test results and the model's predictions of the gradation are presented in Figure 4, which demonstrates that the model's capacity to predict gradation changes during loading.

Using the same set of parameters (Table 1), the model has also been applied to simulate drained extension tests on Cambria sand with confining stresses varying from 6 to 42 MPa (Bopp et al. [9]). A good agreement was also achieved between experimental results and numerical simulations, as shown in Figure 5. With the model, we are capable of re-

<figure><img src="image_3.png" /><figcaption>**Figure 2** Parametric study of particle breakage related parameters. (a) For isotropic compression test and (b)–(c) deviatoric stress and volumetric strain versus major principle strain, respectively, of drained compression test with constant confining stress 26 MPa.</figcaption></figure>

<table><caption>**Table 1** Values of model parameters for Cambria sand</caption><thead><tr><th>G<sub>0</sub>(MPa)</th><th>K<sub>0</sub>(MPa)</th><th>n</th><th>G<sub>p</sub></th><th>φ<sub>μ</sub></th><th>p<sub>y0</sub>(MPa)</th><th>c<sub>p</sub></th><th>e<sub>ref0</sub></th><th>λ</th><th>χ</th><th>ρ</th><th>e<sub>refu</sub></th></tr></thead><tbody><tr><td>35</td><td>26.3</td><td>0.4</td><td>3.5</td><td>37.5°</td><td>12</td><td>0.028</td><td>0.59</td><td>0.006</td><td>15000</td><td>0.16</td><td>0.13</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 3** Comparison between experimental data and numerical simulations for drained triaxial compression tests. (a-c) Deviatoric stress versus major principle strain; (b-d) volumetric strain versus major principle strain.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4** Grain size distributions for drained compression tests. (a) Experiments, and (b) simulations.</figcaption></figure>

producing the stress-strain relation and the volumetric strain response for different confining stress levels.

The evolution of the gradation for all the selected drained extension tests was also predicted. Again, good agreement was achieved between experimental data and simulations as shown in Figure 6.

## 4.3 Simulations of undrained triaxial tests in compression and extension

The model with the same set of parameters was also used to

simulate the undrained compression and extension behaviors of Cambria sand under high consolidated stresses varying from 16.7 to 68.9 MPa for compression tests and from 12 to 52 MPa for extension tests. A good agreement was achieved between experimental and numerical results as presented in Figures 7 and 9.

(i) The initial stress strain curves increases with increasing consolidation pressure.

(ii) Different from drained tests, the peak deviatoric stress is reached at a very low strain level, followed by a distinct reduction in the deviatoric stress.
<figure><img src="image_6.png" /><figcaption>**Figure 5** Comparison between experimental results and simulations for drained triaxial extension tests. (a) Deviatoric stress versus major principle strain, and (b) volumetric strain versus major principle strain.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6** Grain size distributions for drained extension tests. (a) Experiments, and (b) simulations.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7** Comparison between experimental results and simulations for undrained triaxial compression tests. (a) Deviatoric stress versus major principle strain, and (b) effective stress paths.</figcaption></figure>

(iii) The axial strain corresponding to the peak deviatoric stress increases slightly with an increase of the initial consolidation pressure.

(iv) The effective confining stress decreases as a result of the rapidly increasing pore pressure, and the loading resistance of the material becomes reduced, which corresponds to an unstable state also observed in similar tests

performed on loose specimens at low confining stresses.

The experimental and predicted evolutions of the gradation are presented in Figures 8 and 10, which show that the model with parameters determined from drained compression tests can predict the evolution of the gradation during undrained shearing in compression and extension.
<figure><img src="image_9.png" /><figcaption>**Figure 8** Grain size distributions for undrained compression tests. (a) Experiments, and (b) simulations.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9** Comparison between experimental results and simulations for undrained triaxial extension tests. (a) Deviatoric stress versus major principle strain, and (b) effective stress paths.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10** Grain size distributions for undrained extension tests. (a) Experiments, and (b) simulations.</figcaption></figure>

# 5 Conclusions

We analyzed the evolution of gradation as a function of the amount of plastic work and the evolution of the position of CSL with the gradation. Based on these analyses, we suggested two constitutive equations for the relation between the breakage index, the modified plastic work, and the reference critical state void ratio. A double-yield surface mod-

el accounting for the influence of particle breakage was developed, which includes the constitutive equations controlling grain breakage.

Triaxial tests on Cambria sand were used to calibrate and validate the model. The parameters can be easily determined from one isotropic compression test and several drained compression tests. Using the set of determined parameters, several other tests including drained tests in com-
pression and extension and undrained tests in compression and extension were simulated. The grain size distributions at the end of each test were also predicted. All comparisons between the experimental data and the numerical simulations demonstrate that the model reproduces with good accuracy the mechanical behavior of granular materials with particle breakage along various loading paths, as well as the evolution of the grain size distribution during loading.

1. Vesic A S, Clough G W. Behaviors of granular materials under high stresses. ASCE J Soil Mech Found Div, 1968, 94: 661-668

2. Hardin B O. Crushing of soil particles. ASCE J Geotech Eng, 1985, 111: 1177-1192

3. Coop M R. The mechanics of uncemented carbonate sands. Geotech, 1990, 40: 607-626

4. Biarez J, Hicher P Y. Elementary Mechanics of Soil Behaviors. Balkema, ROTTERDAM. 1994

5. Yamamuro J A, Lade P V. Drained sand behavior in axisymmetric tests at high pressures. ASCE J Geotech Eng, 1996, 122: 109-119

6. Lade P V, Yamamuro J A, Bopp P A. Significance of particle crushing in granular materials. ASCE J Geotech Eng, 1996, 122: 309-316

7. Biarez J, Hicher P Y. Influence de la granulométrie et de son évolution par ruptures de grains sur le comportement mécanique de matériaux granulaires. Revue Française de Génie Civil, 1997, 1: 607-631

8. Coop M R, Sorensen K K, Freitas T B, et al. Particle breakage during shearing of a carbonate sand. Geotech, 2004, 54: 157-163

9. Bopp P A, Lade P V. Relative density effects on undrained sand behavior at high pressures. Soils Found, 2005, 45: 15-26

10. Huang W, Ren Q W, Sun D A. A study of mechanical behaviour of rock-fill materials with reference to particle crushing. Sci China Ser E-Tech Sci, 2007, 50(Supp I): 125-135

11. Einav I. Breakage mechanics- Part I: Theory. J Mech Phys Solids, 2007, 55: 1274-1297

12. Muir Wood D, Maeda K. Changing grading of soil: Effect on critical state. Acta Geotechnica, 2008, 3: 3-14

13. McDowell G, Bolton M, Robertson D. The fractal crushing of granular materials. J Mech Phys Solids, 1996, 44: 2079-2101

14. Daouadji A, Hicher P Y, Rahma A. An elastoplastic model for granular materials taking into account grain breakage. Eur J Mech A-Solids, 2001, 20: 113-137

15. Russell A R, Khalili N. A bounding surface plasticity model for sands exhibiting particle crushing. Can Geotech J, 2004, 41: 1179-1192

16. Salim W, Indraratna B. A new elasto-plastic constitutive model for coarse granular aggregates incorporating particle breakage. Can Geotech J, 2004, 41: 657-671

17. Sun D, Huang W, Sheng D, et al. An elastoplastic model for granular materials exhibiting particle crushing. Key Eng Mater, 2006, 340-341: 1255-1248

18. Yao Y P, Yamamoto H, Wang N D. Constitutive model considering sand crushing. Soils Found, 2008, 48: 601-608

19. Thevanayagam S, Shenthan T, Mohan S, et al. Undrained fragility of clean sands, silty sands and sandy silts. ASCE J Geotech Geoenviron Eng, 2002, 128: 849-859

20. Muir Wood D, Kikumoto M, Russell A R. Particle crushing and deformation behavior. In: Proceedings of Prediction and Simulation Methods for Geohazard Mitigation. Oka F, Murakami A, Kimoto S, eds. CRC Press/Balkema, 2009

21. Daouadji A, Hicher P Y. An enhanced constitutive model for crushable granular materials. Int J Numer Anal Methods Geomech, 2010, 34: 555-580

22. Richard F E, Hall J R, Woods R D. Vibration of Soils and Foundations. International Series in Theoretical and Applied Mechanics. Prentice-Hall: Englewood Cliffs, NJ, 1970

23. Vermeer P A. A double hardening model for sand. Geotech, 1978, 28: 413-433

24. Yin Z Y, Chang C S, Hicher P Y. Micromechanical modelling for effect of inherent anisotropy on cyclic behaviour of sand. Int J Solids Struct, 2010, 47: 1933-1951

25. Ishihara K, Towhata I. Cyclic behavior of sand during rotation of principal axes. Mech Granular Mater, Elsevier, 1983. 55-73

26. Luong M P. Stress-strain aspects of cohesionless soils under cyclic and transient loading. In: Proc of Int Symp on Soils under Cyclic and Transient Loading, Swansea, 1980. 353-376

27. Sheng D, Sloan S W, Yu H S. Aspects of finite element implementation of critical state models. Comput Mech, 2000, 26: 185-196