A State-Dependent Constitutive Model for Unsaturated Rockfill Materials

Abstract

This chapter presents a state-dependent elastoplastic constitutive model for both saturated and unsaturated rockfill materials. The model, which is developed within an extended critical-state framework, uses two independent stress state variables: total stress and total suction. The loading-collapse (LC) curve proposed by Oldecop and Alonso for unsaturated rockfills is used herein. A unified hardening parameter, which could consider the effects of stress level, internal state (density) and relative humidity, is introduced to describe the state-dependent dilatancy of saturated and unsaturated rockfill materials. The details of the model formulation and parameters determination are described and reported. Numerical simulations on the triaxial tests, such as the drained shear tests on the saturated specimens with different initial dry densities, shear tests on the specimens with different relative humidity and wetting deformation tests under constant vertical strain rate, have been carried out using the proposed model. The numerical results show that the stress-strain relationships at both loose and dense, saturated and unsaturated states can be properly modelled with a single set of parameters. Additionally, the proposed model can also capture some other key features such as the strain-softening behaviour at the dense state and low confining stress, the sudden stress relaxing subjected to the flooding under a constant vertical strain.

1. Introduction

It is well recognised from the experimental studies and engineering practice that the influence of water on the mechanical behaviour of rockfill materials is significant. A lot of laboratory experiments carried out in the past showed that the rockfill materials may undergo significant amounts of strain upon flooding [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. This phenomenon is usually called ‘collapse deformation’, referring to a strain increment not related to the load increment. On the other hand, a remarkable influence of the compacted density and confining pressure on the strength and deformation behaviours of rockfill materials was observed in some large-scale triaxial tests [11, 12, 13, 14, 15]. The experimental results on rockfill materials at a dense state showed the remarkable strain-softening and dilatancy behaviours, while the rockfill materials at a loose state exhibited the strain hardening and volumetric contraction behaviours. This state-dependent behaviour was also observed in the triaxial tests on rockfill materials used in Guanyinyan dam that will be introduced in the next section. Furthermore, compared with the fully saturated rockfill materials, a more obvious state-dependent behaviour was observed for unsaturated rockfills in the relative humidity-controlled triaxial tests [7, 16]. The rockfills subjected to the larger suction exhibited more remarkable strain-softening and dilatancy at a low confining pressure or with a small initial void ratio. Thus, presenting a constitutive model which could consider the collapse deformation and state-dependent behaviour for unsaturated rockfill materials is necessary.

Rockfill materials have been widely used in the constructions of rockfill dams, and they should qualitatively be defined very well before beginning to the design stage [17, 18, 19, 20]. The collapse settlement of upstream shell for clay core rockfill dam occurs during the impounding, which has attracted extensive attentions from both engineers and scientists since 1970s [3], for it may cause cracks on the dam crest and result in damage to impermeable systems. Thereby, proposing a satisfactory technique for the prediction of collapse settlement is clearly important. Nobari and Duncan [3] first introduced collapse effects into finite element (FE) analysis to predict collapse settlement in rockfill dams. This calculation starts by performing the FE analysis using a set of material parameters corresponding to the dry conditions; then the collapse is numerically simulated in two stages. First, the stress change caused by the saturation at the constant strain is determined using the experimental data. In the second stage, the nodal forces that restore the equilibrium are applied, and the computed strains reproduce the collapse strains. Naylor et al. [21] generalised this method for arbitrary constitutive models, which required the knowledge of two sets of constitutive parameters for dry and saturated conditions. However, numerous experiments have invalidated the implicit assumption that the sequence of loading and wetting does not affect the final state of materials and led to several attempts to establish other approaches [22, 23, 24, 25]. Shen and Wang [22] expressed the magnitudes of the wetting-induced volumetric strain and deviatoric strain using empirical functions, including both confining pressure and deviatoric stress level during the wetting process. By virtue of the assumption of coaxiality of the strain rate tensor and the stress tensor, the strain components are estimated and translated into nodal forces, which are further applied to the finite elements to obtain the wetting deformation of the dam. Recently, this so-called single curve method is widely used in dam engineering in China due to its practicality. Nevertheless, the point to underline here, in the above two procedures, is essentially a computational device that is not necessarily associated with a physical reality [26, 27]. The mechanical properties of rockfill materials are closely related to the breakage properties of rock particles [4, 28] and these breakage properties are sometimes remarkably affected by the degree of saturation [1, 10, 29]. Hence, the physical mechanism underlying the collapse deformation seems to be the difference of breakage rates between saturated and unsaturated rock particles. Therefore, introducing the effect of relative humidity on the rockfill performance is by no means an academic one to model the collapse deformation. Oldecop and Alonso [1] first proposed an elastoplastic model that was valid to express compression properties of unsaturated rockfills considering the influence of relative humidity. Thereafter, within the framework of the Barcelona Basic Model [30] for unsaturated soils, they extended this compression model to the three-dimensional stress state and applied it in the analyses of Beliche and Lechago rockfill dams [31, 32]. Furthermore, Chávez and Alonso [7] developed a work-hardening model to describe the effect of particle breakage at different relative humidities. Besides, Kohgo et al. [33] proposed a three-dimensional elastoplastic model for unsaturated rockfills according to a modification of their generalised elastoplastic model for unsaturated soils. Bauer [34] extended his hypoplastic model by introducing a moisture-dependent solid hardness and taking the stress relaxation caused by the change of solid hardness into consideration. The state-dependent behaviour was not considered in these models.

On the other hand, due to the obvious state-dependent behaviour of sands, many state-dependent models have been proposed [35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45]. Among them, the state-dependent model proposed by Li and Dafalias [41], with only a set of model parameters, can precisely represent the complex stress-strain relationship of Toyoura sand over a wide range of densities and pressures. Within the framework of this model, Xiao et al. [46] developed a generalised elastoplastic model to describe the state-dependent dilatancy for rockfill materials. Sun et al. [47] introduced the fractional order method to describe the state-dependent behaviour of granular materials without using the plastic potential function. However, to the best of our knowledge, these models cannot describe the combination effect of the initial dense state and prevailing relative humidity on the stress-strain relationship of rockfill materials.

The purpose of this chapter is to present an elastoplastic model for unsaturated rockfill materials by introducing a unified hardening parameter which could well capture the essential physical phenomena behind the observed collapse deformation and state-dependent dilatancy for saturated and unsaturated rockfill materials. The conceptual bases for the model and the basic formulation are given in Section 2. Thereafter, the proposed model is incorporated into the coupled flow-deformation analysis FE code [48, 49, 50, 51] using the explicit stress integration algorithm for the elastoplastic model of unsaturated soils [52, 53]. Finally, the comparisons are conducted using the experimental data from a series of triaxial tests on saturated and unsaturated rockfill materials, and the capability of the proposed model to capture the state-dependent behaviour and collapse deformation is validated.

2. Constitutive framework and formulations

2.1 Basic constitutive variables

Due to the large permeability for rockfills, water never fills the large rockfill voids unless structures are submerged. Environmental changes and rainfall can only modify the relative humidity. The extreme case of rockfill flooding is a concern in the rockfill structures which become inundated, such as the upstream rockfill shells of the zoned earth and rockfill dams. Hence, a suitable stress space to describe the isotropic compression states of unsaturated rockfill materials is (p, RH), where p is the mean total stress and RH is the relative humidity. In fact, the relative humidity is the ratio of the vapour pressure present in the air to the vapour pressure when the air is saturated with water vapour. The relative humidity in the gas phase and the matric suction, s, in the rockfill pore water are related by the psychrometric relationship [52]:

$$\psi =s+\pi =-\frac{\mathit{RT}{\rho }_w}{M_w}ln\left(\mathit{RH}\right)$$E1

where R is the gas constant, T is the absolute temperature of the reference system, M_w is the molecular mass of water and ρ_w is the density of water at the reference temperature. π is called the osmotic suction, which is due to the presence of solutes in the rockfill pore water. The sum ψ = s + π is called the total suction. Total and matric suction would be equal with each other in the case when the rockfill pores contained only pure water with no solutes, for instance, rockfill materials applied in the dams. The total stress and total suction are chosen as the basic constitutive variables herein.

2.2 Normal compression model

Considering an isotropic compression test on unsaturated rockfills, in which a rockfill sample at a given relative humidity (subsequently maintained) is subjected to p-load increments along virgin states, which is in full correspondence with the compression behaviour of saturated rockfills, the specific volume will be given by

$$v=1+e=N-\lambda \left(\psi \right)p$$E2

where ν is the specific volume, e is the void ratio and N is the intercept of the normal compression lines with the ν-axis when p = 0. The compression index λ(ψ) is assumed to be a function of the total suction and can be interpolated from the compression indices at the fully saturated state (ψ = 0) and very dry state. For rockfill materials, the very dry state is usually defined as the water content close to zero, generally corresponding to the water content less than 0.45% in practice. A schematic representation of Eq. (1) can be discretized into a series of normal compression lines, as shown in Figure 1.

Figure 1.

Normal compression lines for saturated and unsaturated rockfill materials.

As mentioned by McDowell and Bolton [53], the plastic deformation of the granular materials can be attributed to the two deformation mechanisms. Under low stress, the plastic deformation is due to the particle rearrangement. The second mechanism is called clastic yielding, and it is attained when the applied stress causes the onset of particle breakage. According to this concept and corresponding experimental results, an elastoplastic compressibility model for rockfill was developed by Oldecop and Alonso [1]. Under isotropic compression, a mean threshold total stress, p_y, which marks the beginning of particle breakage, was introduced in the model. Below this threshold stress, the influences of suction on the compression index were not measured in the compression tests. On the other hand, the effect of water action on the compression behaviour was strong when beyond this threshold stress, which is the source of the collapse deformation. Note that, the difference of compression index below and beyond this threshold stress is vanished when the material is in the very dry state. Therefore, the compression index for unsaturated rockfills can be expressed as

$$\begin{array}{l}p\le p_y\lambda \left(\psi \right)={\lambda }^i\\ p>p_y\begin{array}{cc}& \lambda \left(\psi \right)={\lambda }^i+{\lambda }^d\left(\psi \right)\end{array}\end{array}$$E3

where λⁱ is the slope of the normal compression line when particle rearrangement is active only, and (λⁱ + λ^d) is the slope of the normal compression line when both plastic deformation mechanisms (particle rearrangement and particle breakage) are active. According to the previous literature [1], the compressibility index accounting for the clastic yielding is linearly related to the logarithm of total suction through a material parameter:

$${\lambda }^d\left(\psi \right)={\lambda }_0^d-{\alpha }_{\psi }ln\left(\frac{\psi +p_{\mathit{atm}}}{p_{\mathit{atm}}}\right)$$E4

where ${\lambda }_0^d$ is the maximum clastic compressibility index in the saturated condition, p_atm is the atmospheric pressure and α_ψ is a model parameter. The elastic volume change for rockfill materials within unloading-reloading (URL) paths are given by

$$\mathit{dv}=-\mathit{\kappa dp}$$E5

where the slope κ is assumed to be independent of the water action. In addition, there was a moderate swelling behaviour observed in the experiments on the unsaturated rockfills, which was caused by the increase in water content. This swelling is assumed to be a reversible elastic deformation, and a linear relationship is assumed between the swelling volume change and the logarithm of total suction:

$${\mathit{dv}}^{\psi }=-{\kappa }_{\psi }\frac{\mathit{d\psi }}{\left(p_{\mathit{atm}}+\psi \right)}$$E6

where κ_ψ is the total suction-based expansion/compression index. For simplicity, κ_ψ is assumed to be independent of the stress level.

2.3 Loading-collapse yield curve

For the isotropic compression condition, the yield stress, p₀, was defined by the following expression such as [1]:

$$\left\{\begin{array}{l}{p}_0^{\ast }\le p_y\Rightarrow p_0\left(\psi \right)=p_0^{\ast }\\ p_0^{\ast }>p_y\Rightarrow p_0\left(\psi \right)=p_y+\frac{\left({\lambda }^i-\kappa \right)\left(p_0^{\ast }-p_y\right)}{{\lambda }^i+{\lambda }^d\left(\psi \right)-\kappa }\end{array}\right.$$E7

Eq. (7), which is referred to as loading-collapse (LC) yield surface, describes the relationship between the equivalent yield stress ($p_0^{\ast }$ defines the position of the yield curve, and it was identified as the yield stress for a very dry rockfill) and the yield stress (p₀) in unsaturated states, which is controlled by the total suction. By introducing the compressibility index for the clastic yielding expressed in Eq. (4) into Eq. (7), the LC yield surface is then obtained. A plot of the LC yield curves for different values of $p_0^{\ast }$ for Pancrudo slate that were tested by Oldecop and Alonso [1] is given in Figure 2. There are two regions in which the yield stress, p₀, does not depend on the total suction, which correspond to the stress levels below the clastic yield stress py ($p_0^{\ast }$< 0.29 MPa) and to the very dry states (ψ > 67 MPa), respectively.

Figure 2.

Load-collapse curves for different values of $p_0^{\ast }$.

Thereafter, a simple volumetric hardening is used here to follow the evolution of $p_0^{\ast }$:

$${\mathit{dp}}_0^{\ast }=\frac{\mathit{dH}}{{\lambda }^i-\kappa }$$E8

where H is the unified hardening parameter related to the plastic volumetric strain and state-dependent behaviour, which will be introduced in Section 2.5.

2.4 Yield surface

Yao et al. [54] developed a simple but robust constitutive model (the original UH model) to reproduce the mechanical behaviours in the isotropic and triaxial states for both normally consolidated (NC) and overconsolidated (OC) soils, such as non-elastic deformation in reloading, peak strength, shear-dilation and strain-softening. In his model, a unified hardening parameter (H) was proposed, and an associated flow rule was adopted for the robust and convenience of numerical implementation. Hence, in accordance with the aim of simplicity, the original UH model is adopted as the base for developing a new state-dependent constitutive model of unsaturated rockfill materials. Thus, the yield and plastic potential surfaces in the stress space (p, q and ψ) are defined as follows:

$$f=g=\left\{\begin{array}{l}p+\frac{q^2}{M{}^{2}p}-\frac{H}{{\lambda }^i-\kappa }\begin{array}{cc}& \mathit{\text{if}}\begin{array}{c}{p}_0^{\ast }\le p_y\end{array}\end{array}\\ p+\frac{q^2}{M{}^{2}p}-\frac{{\lambda }^d\left(\psi \right)p_y-H}{{\lambda }^i+{\lambda }^d\left(\psi \right)-\kappa }\begin{array}{cc}& \mathit{\text{if}}\begin{array}{c}{p}_0^{\ast }>p_y\end{array}\end{array}\end{array}\right.$$E9

where q is the deviator stress and M is the critical-state slope.

With respect to the rockfill materials, the critical-state slope, M, is no longer a constant. Experimental results on rockfill materials show that the critical state friction angle is dependent on the mean stress [19, 55, 56, 57] that decreases with an increase in p. Moreover, according to the RH-controlled triaxial tests on rockfill materials in Chávez and Alonso [7], M was also found to be dependent on the RH, and the relationship between M and logσ₃ for samples with different RH seemed to be a series of parallel lines. From the experimental data [7, 16], the value of M increases with the increase of total suction (decrease RH). Besides, considering the critical stress ratio is impossible to be negative with the increase of mean stress, the following function is adopted for the sake of simplicity:

$$M\left(p,\psi \right)=\left[M_{\mathit{res}}+\left(M_0-M_{\mathit{res}}\right)e^{-{\alpha }_{M}p}\right]\left[1+{\beta }_{\psi }ln\left(\frac{\psi +p_{\mathit{atm}}}{p_{\mathit{atm}}}\right)\right]$$E10

where M_o and M_res are the initial and residual stress ratios for rockfills in saturated condition under low (0.1 MPa) and high confining pressures, respectively; α_M and β_ψ are material constants. Figure 3 shows that Eq. (10) can well describe the change of the critical stress ratio with mean pressure and relative humidity for unsaturated rockfill materials.

Figure 3.

Variations of critical stress ratio with mean stress and relative humidity: (a) crushed cambric slate and (b) Pancrudo slate.

Using above equations, a three-dimensional view of the yield surfaces in the (p, q and ψ) space is given in Figure 4, where ψ_i is the total suction corresponding to the very dry state.

Figure 4.

Three-dimensional yield surface in stress space (p, q and ψ).

2.5 Unified hardening parameter

The unified hardening parameter was first proposed to provide a unified description of the mechanical behaviours for both clays and sands by Yao et al. [58]. It was then modified and employed to model the hardening process for both normally consolidated and overconsolidated clays [54, 59]. The unified hardening parameter, which is adopted to describe the hardening of the yield surface of rockfill materials can be written as

$$H=\int \mathit{dH}=\int \frac{M_f^4-{\eta }^4}{M^4-{\eta }^4}d{\epsilon }_v^p$$E11

where η is the stress ratio, η = q/p; ${\epsilon }_v^p$ is the plastic volumetric strain; and M_f is the potential failure stress ratio, which represents the potential peak strength that changes with the compaction density and relative humidity. Since dH is always larger than or equal to zero, the following conclusions can be drawn from Eq. (11):

1. 0 ≤ η < M (negative dilatancy condition): $d{\epsilon }_v^p>0$.

2. η = M (characteristic state condition): $d{\epsilon }_v^p=0$.

3. M ≤ η < M_f (positive dilatancy condition): $d{\epsilon }_v^p<0$.

As indicated above, the dilatancy of rockfill materials can be reasonably described by this hardening parameter.

The next is introducing the determination of potential failure stress ratio. For rockfill materials, the peak friction angle, ϕ_p, on a contact plane is dependent on the degree of interlocking by neighbouring particles, which can be related to the state of the packing void ratio [60]:

$$tan{\varphi }_p={\left(\frac{e_c}{e}\right)}^{m}tan{\varphi }_{\mu }$$E12

where ϕ_μ is the internal friction angle corresponding to the critical state, which has the relationship with the critical stress ratio M as shown in Eq. (17); and m is a material constant when granular materials are in the fully saturated condition, but it varies with the confining pressure and suction in unsaturated conditions according to the experimental observations, which can be corrected using the following function:

$$m\left({\sigma }_3,\psi \right)=m_2+\left(m_1-m_2\right){\left(\frac{m_2}{m_1}\right)}^{\psi {\left(\frac{{\sigma }_3}{p_{\mathit{atm}}^2}\right)}^n}$$E13

where σ₃ is the confining pressure; m₁, m₂ and n are material constants. From Eq. (13), we can find that parameter m is kept as m₁ in the saturated condition, but decreases with the increasing confining pressure when in the unsaturated condition. Thus, the nonlinear variation of the peak strength with the changes of suction and confining pressure can be well described using this function.

Furthermore, e_c in Eq. (12) is the void ratio corresponding to the critical state, which is a function of the mean stress. However, for the rockfills in the unsaturated condition, the influence of the total suction on the critical void ratio should be considered as well. From the experimental data in [7, 16], the relation between critical stress ratio and mean stress for saturated and unsaturated rockfills seems to be two parallel lines. Thus, the critical void ratio can be expressed based as [61].

$$e_c\left(p,\psi \right)=e_{\mathit{ref}}\left(\psi \right)-\zeta log\left(\frac{p}{p_{\mathit{ref}}}\right)$$E14

where ζ is the material constant; (e_ref, p_ref) is a reference point on the critical state line e_ref is a function of total suction, which nonlinearly increases with the increase of total suction. For the sake of simplicity, the following function is used to express the variation of reference critical void ratio with total suction:

$$e_{\mathit{ref}}\left(\psi \right)=e_{\mathit{ref}}^0+{\alpha }_{e}ln\left(\frac{\psi +p_{\mathit{atm}}}{p_{\mathit{atm}}}\right)$$E15

where $e_{\mathit{ref}}^0$ means the reference critical void ratio in the saturated condition; and α_e is the material constant, which describes the increase rate of reference critical void ratio with total suction.

In the triaxial compression condition, the stress ratios, M_f and M, can be expressed by using peak friction angle (ϕ_p) and internal friction angle (ϕ_μ) as

$$M_f=\frac{6sin{\varphi }_p}{3-sin{\varphi }_p}$$E16

$$M=\frac{6sin{\varphi }_{\mu }}{3-sin{\varphi }_{\mu }}$$E17

Combining Eqs. (12), (16) and (17), the relation between the potential stress ratio (M_f) and critical state stress ratio (M) can then be obtained.

Substituting Eq. (16) into Eq. (11), the hardening parameter is then obtained, which can reflect the influence of the compacted density and relative humidity on the mechanical behaviour of unsaturated rockfill materials. Considering the peak frictional angle, ϕ_p, is generally greater than ϕ_μ, especially in the dense and very dry states, the dilatancy behaviour for the rockfill materials can thus be well captured using the modified hardening parameter. In addition, once an obvious dilatancy occurs for the rockfill in the low confining pressure, the degree of interlocking and the peak frictional angle are reduced, and the strain-softening phenomenon can also be well described by using this hardening parameter.

2.6 Elastic moduli

Considering the nonlinear changes of the elastic moduli with the changes of degree of compaction and stress level, the elastic shear modulus, G, is calculated by using the following empirical Equation [46]:

$$G=G_0\frac{{\left(2.97-e\right)}^2}{1+e}{\left({\mathit{pp}}_{\mathit{atm}}\right)}^{0.5}$$E18

where G₀ is a material constant. In this equation, the current void ratio, e, is used instead of the initial void ratio. The elastic bulk modulus, K, is equal to

$$K=\frac{2\left(1+\nu \right)}{3\left(1-2\nu \right)}G$$E19

where ν is the Possion’s ratio.

2.7 Determination of model parameters

The application of the model requires the information on the following stress states and parameters.

1. Initial state: initial stresses (p_i, q_i and s_i), initial void ratio, e₀, and initial mean yield stress, $p_0^{\ast }$.

2. Parameters directly associated with the isotropic compression behaviour: p_y, threshold yield mean stress for the onset of clastic phenomena; λⁱ, normal compressibility index when instantaneous deformation mechanism is active only; ${\lambda }_0^d$, maximum clastic compressibility index for fully saturated conditions; κ, compressibility coefficient along elastic (unloading-reloading) stress paths and α_s, compressibility parameter, which controls the rate of stiffness increase with total suction.

3. Parameters directly associated with the critical state: M₀, maximum critical stress ratio for the mean stress approaching to 0.1 MPa; M_res, residual critical stress ratio corresponding to the large mean stress; α_M, controls the rate of critical stress ratio decrease with mean stress; β_ψ, controls the rate of critical stress ratio increase with total suction; ($e_{\mathit{ref}}^0$, p_ref) is a reference point on the critical state line in the saturated conditions; χ, controls the rate of critical void ratio decrease with mean stress and α_e, controls the rate of critical void ratio increase with total suction.

4. Parameters directly associated with elastic moduli: G₀, material constant associated with shear modulus within the elastic domain; ν, Possion’s ratio, which is assumed to be independent of suction and taken as 0.2 ∼ 0.3 for rockfill materials and κ_ψ, compressibility coefficient for changes in suction within the elastic region.

5. Other default parameters: m₁, controls the rate of peak strength change with void ratio in the saturated conditions; m₂, minimum value of parameter, m, corresponding to the high confining stress at very dry state and n, controls the vary of m with suction and confining pressure.

In general, the determination of the model parameters will require relative humidity-controlled testing methods, and suggested stress paths for different sets of parameters are the following:

1. Tests that involve isotropic compression (loading and unloading) for very dry specimen and saturated ones; they provide data to find λⁱ, $p_0^{\ast }$, ${\lambda }_0^d$ and κ.

2. Tests that involve isotropic compression for very dry specimen subsequent with flooding under low applied stress and constant stress p₀; they provide data to find κ_ψ and p_y. As mentioned in the previous study [1], p_y can be obtained using the following function:

$$p_y=p_0-\frac{{\epsilon }_{\mathit{\text{collapse}}}+{\epsilon }_{exp\mathit{\text{ansion}}}}{{\lambda }_0^d}$$E20

where ε_expansion is the measured total expansion strain due to flooding under low applied stress and ε_collapse is the measured total collapse strain due to flooding under applied stress p₀.

1. Tests that involve isotropic compression under different relative humidity from what a linear relationship between collapse strains and the logarithm of rockfill water content can be observed; Oldecop and Alonso [1] provided a method to find α_ψ using the following function

$${\alpha }_{\psi }=\frac{{\kappa }_{\psi }+{\chi }_{\psi }}{p_0-p_y}$$E21

where χ_ψ is the experimental coefficient relating collapse strain with total suction at constant stress p₀ (p₀ > p_y).

1. Triaxial tests (loading, unloading and reloading) at different relative humidities and saturated condition, and they provide data to find M₀, M_res, α_M, β_ψ, $e_{\mathit{ref}}^0$, p_ref, χ, α_e and G₀.

Note that the default material parameters, m₁, m₂ and n, should be determined according to the back analysis using the experimental data by fixing the previous determined parameters. Besides, the above-mentioned test is the minimum experimental programme, so the more tests are required to determine the more reliably model parameters.

3. Model validations

In order to verify the capability and feasibility of the proposed constitutive model, it was incorporated into the hydro-mechanical coupling FEM code [48] that was developed in the Institute of Hydraulic Structures of Hohai University [49, 50]. The explicit stress-integration algorithm for unsaturated soils proposed by Sheng et al. [51, 62] was used for the numerical implementation of the proposed constitutive model. The triaxial tests on the saturated specimens with different initial void ratio, triaxial tests with different relative humidity and triaxial wetting tests under constant vertical strain rate were simulated by means of a hexahedral element. The confining pressure was applied against the upper and lateral boundaries, a constant vertical displacement rate was then applied to the upper boundary for shearing and the axial stresses and volumetric strains were computed finally. For the tests at dry state, a constant high suction, which represents the laboratory conditions, was maintained in the specimen. Flooding was simulated by converting negative pore-water pressures (total suctions) of element nodes into zero.

3.1 Behaviour of rockfill compacted to various initial void ratios

At first, the model was used to predict drained triaxial tests on saturated rockfills with different initial void ratios. The mixture of weak and strong weathered limestone, which was selected as the rockfill material for the Guanyinyan dam (Panzhihua, China), was used in the triaxial tests. The gradation of the tested material was scaled down by making the specimen’s cumulative grain size distributions parallel to the gradation curve of the in situ rockfill material combining with the equivalent substitution method. Figure 5 shows the grain size distribution. The initial void ratio was controlled as 0.24, 0.28 and 0.37 in the sample preparation, respectively. Before shearing, the samples were saturated and isotropically compressed to the confining pressures ranging from 0.2 to 1.2 MPa. The model parameters, which are listed in Table 1, were determined using the procedure introduced in the above section based on the experimental stress-strain relationships. Note that, rockfill specimens for laboratory tests were prepared by compacting the material in thin layers to obtain uniform samples [63]. Thus, a significant pre-consolidation effective mean stress ($p_0^{\ast }$ = 0.3 MPa), associated with compaction, was used in the model to simulate the triaxial tests.

Figure 5.

Grain size distribution curves.

Elastic parameters	Isotropic compression parameters	Critical state parameters	Default parameters
G0 = 125	λ^i-κ = 0.001 MPa^−1	M0 = 1.95	m1 = 0.36
ν = 0.28	${\lambda }_0^d$ = 0.02 MPa^−1	Mres = 1.55	m2 = 0.1
κψ = 0	αψ = 0.001 MPa^−1	αM = 1.2 MPa^−1	n = 3.5
	py = 0.01	βψ = 0.01
	$p_0^{\ast }$ = 0.3 MPa	$e_{\mathit{ref}}^0$ = 0.51
	e0 = 0.24, 0.28 and 0.37	pref = 0.02 MPa
		ζ = 0.105
		αe = 0.03

Table 1.

Model parameters for tested limestone.

The comparisons between the numerical results and experimental measurements at different confining pressures are shown in Figure 6. It can be found that proposed model is satisfied to capture the key features of the loose and dense saturated specimens with one unique set of soil parameters. For instance, the denser of the specimen, the larger the deformation modulus and shear strength. Besides, it clearly shows that the rockfill with a smaller initial void ratio (denser state) dilates more remarkably due to the difference between the peak stress ratio and critical stress ratio is much greater, especially at a low confining pressure. With the increase of initial void ratio, the difference between the peak stress ratio and critical stress ratio decreases, which leads to the decrease of volumetric dilatancy degree. It is worth to note that there is an obvious strain-softening phenomenon observed for the dense specimen under the low confining pressures. To our satisfaction, the proposed model can also predict this strain-softening to some extent, although there is a little difference which may be caused by the difficult to determine the parameters corresponding to the critical state.

Figure 6.

Measured and predicted stress-strain relationship (q vs. εa and εv vs. εa) for rockfill with different initial void ratios under confining pressures: (a) σ₃ = 0.2 MPa, (b) σ₃ = 0.4 MPa, (c) σ₃ = 0.8 MPa and (d) σ₃ = 1.2 Mpa.

Figure 7 shows the predicted relationship between the potential failure stress ratio, M_f, and stress ratio, η, for three different initial void ratios when σ₃ = 0.2 MPa, which can be used to explain the strain-softening phenomenon. From the initial state to the final state, the state η = M_f (i.e., the peak strength state) appears, shown as the intersections of the 45⁰ line and the M_f (η) lines for different initial void ratios. Before the peak strength state is reached, η < M_f, the rockfill material undergoes a hardening process and remains in a stable state. After the peak strength state, η becomes slightly larger than M_f, and then the rockfill material undergoes a softening process that is unstable. It can be found that the part of η > M_f becomes more remarkable with the decrease of the initial void ratio, which explains the appearance of the most significant strain-softening process for the densest specimen (e₀ = 0.24). Additionally, Figure 8 presents the influence of the confining pressure on the changes of potential failure stress ratio with the changes of stress ratio when e₀ = 0.24. It can be found from this figure that M_f (η) line for the larger confining pressure is always in the left side of the 45⁰ line (η < M_f) from the initial state to the final state. Thus, there is no softening phenomenon appearing when there is a large confining pressure applied on the rockfills, and the deviatoric stress only increases up to the strength corresponding to the critical state.

Figure 7.

Relationships between potential failure stress ratio, mf, and stress ratio, η, for three different initial void ratios when σ₃ = 0.2 MPa.

Figure 8.

Relationships between potential failure stress ratio, mf, and stress ratio, η, for different confining pressures when e₀ = 0.24.

Noting that, the difference of stress-strain curve between the measured and predicted is relatively larger with the increase of confining pressure, which is shown in Figure 6(c) and (d). Compared with that experimental stress-strain response, the predicted deviatoric stress increases more slowly with the increase of axial strain. From our analysis, this may be caused by the significant particle breakage under the high confining pressures for the rockfills used in the tests, which led to the considerable increase of compressibility modulus before shearing. However, the influence of the particle breakage on the compressibility modulus of rockfill materials has not been considered in the proposed model. In the construction of model, a linear relationship in the e-p plane is assumed, which means the compression index is kept unchanged throughout the isotropic and triaxial compression process. Even so, the quantitative results are in general in good agreement with the experimental measurements. Finally, from above comparisons, it is indicated that the behaviours of the state-dependent dilatancy and strain-softening for rockfill materials can be well captured by the proposed model.

3.2 Behaviour of rockfill with various RH

In order to illustrate the capability of proposed model to predict the stress-strain relationships of unsaturated rockfills, the results of RH-controlled triaxial tests on the crushed cambric slate were adopted as comparison. These tests were conducted by Chavez and Alonso [7], which involved two series of tests. The first series was performed in the fully saturated condition. After compaction, the samples were confined under pressures of 0.1, 0.3, 0.5 and 0.8 MPa and then sheared under the constant confining pressure. In a second series, samples were compacted in a similar manner as saturated ones and maintained at a constant low relative humidity (RH = 36%, equivalent to the total suction, s = 142 MPa). After that, similar stress paths were applied to these dry samples. The model parameters are given in Table 2. Most of parameters were taken from the Chavez and Alonso [7], while the rest were determined according to the experimental stress-strain curves.

Elastic parameters	Isotropic compression parameters	Critical state parameters	Default parameters
G0 = 80	λ^i-κ = 0.001 MPa^−1	M0 = 2.1	m1 = 0.36
ν = 0.29	${\lambda }_0^d$ = 7.48 × 10^−2 MPa^−1	Mres = 1.45	m2 = 0.1
κψ = 0.4 × 10^−3	αs = 7.86 × 10^−3 MPa^−1	αM = 1.8 MPa^−1	n = 3.5
	py = 0.01 MPa	βψ = 0.01
	$p_0^{\ast }$ = 0.3 MPa	$e_{\mathit{ref}}^0$ = 1.0
	e0 = 0.6	pref = 0.01 MPa
		ζ = 0.1
		αe = 0.03

Table 2.

Model parameters for crushed cambric slate.

Figure 9 shows the predicted and measured stress-strain curves in the conditions of fully saturated and dry (RH = 36%). For the saturated sample under a low confining pressure (0.1 MPa), the deviatoric stress increases with the axial strain up to a constant value for axial strains reaching around 10% and then remains constant, which shows the hardening process. However, for the dry sample with the same confining pressure, a peak strength is observed for axial strain reaching the value between 4 and 6%, and a little strain-softening appears after that. On the other hand, when the sample is sheared under the high confining pressures, the increase of deviatoric stress with the axial strain continues at the end of triaxial tests for both saturated and dry samples. Thus, it is difficult to determine the critical state condition in this case. It can be found a good agreement of deviatoric curves (q vs. ε_a) between modelling and experimental data at both saturated and dry states. Besides, the unloading-reloading process can be well predicted as well with the proposed model. However, some difficulties are observed in reproducing the volumetric response by this model. Computed dilation rates, especially at the low confining pressures, tend to underestimate measured values. Experimental measured dilation rates are significantly larger than that predicted, especially for the unsaturated rockfills under the low confining pressures. There are two reasons accounting for this difference. The first is possibly due to the assumption of critical state conditions in the determination of model parameters. In fact, the attainment of the critical state conditions for a gravel material is probably outside the possibilities of a triaxial test. Very larger strains and a fully broken gravel structure are probably required, but strain localization is also the failure mode in most of tests, which will complicate the interpretation. Thus, the proper observation of critical state conditions would require the larger deformations which were not achieved in the triaxial tests. The second relevant reason is the shape of particles in a gravel. The tested rockfill material has elongated particles with sharp edges, so interlocking and dilation seem very likely to happen even after some significant breakage and deformation, and this also increases the difficulty to modelling the volumetric strains.

Figure 9.

Comparison between model predications and experimental measurements: Triaxial tests on rockfill in the fully saturated and dry states (RH = 36%). (a) σ₃ = 0.1 MPa. (b) σ₃ = 0.3 MPa. (c) σ₃ = 0.5 MPa. (d) σ₃ = 0.8 MPa.

3.3 Behaviour of rockfill during shearing and wetting

Finally, two sets of triaxial wetting tests on rockfills considering the suction reduction [32] were adopted as comparison to verify the capability to predict the collapse deformation for the proposed model. The tested materials were schists and greywacke, respectively, which were used as the rockfill materials for the inner and outer shell of Beliche Dam. Tests were performed on a dry material first (presumably equilibrated with the relative humidity prevailing in the laboratory, for total suction equal to 20 MPa).

At some stage during the triaxial tests, specimens were then flooded. The applied vertical strain rate was maintained during the sample inundation and subsequent straining. The stress-suction path actually experienced by the specimens in the stress space (p, q and ψ) is illustrated in Figure 10. The compacted specimen is initially loaded isotropically under very dry condition (point A). The specimen is then loaded vertically at constant suction until the limiting condition (point B) is reached. Thereafter, it was flooded with the constant of vertical strain rate maintained. The volumetric collapse induces a sudden loss of vertical stress for the reduction of suction to zero, and Point C is reached. As the total suction reaching zero, the yield envelope experiences a reduction in size. As the vertical strain rate is maintained, the specimen starts resisting again along the loading path and eventually reaches a new yield condition at point D, and then a final limiting state at point E. The model parameters for schist and greywacke are listed in Table 3, for most of the parameters are referred from Alonso et al. [32], and the rest are determined based on the corresponding experimental results combined with back-analysis method.

Figure 10.

Sketch showing stress-suction path applied in triaxial tests initially dry and later flooded when limiting conditions are reached.

Elastic parameters	Isotropic compression parameters	Critical state parameters	Default parameters
Inner schists/outer greywacke
G0 = 80/100	λ^i-κ = 0.025/0.01 MPa^−1	M0 = 1.75/1.8	m1 = 0.36/0.4
ν = 0.3/0.3	${\lambda }_0^d$ = 0.028/0.01 MPa^−1	Mres = 1.2/1.6	m2 = 0.16/0.2
κψ = 0/0	αs = 7.0/2.0 (×10^−3 MPa^−1)	αM = 0.1/0.5 MPa^−1	n = 3.5/3.5
	py = 0.01/0.01 MPa	βψ = 0.03/0.01
	$p_0^{\ast }$ = 0.3/0.3 MPa	$e_{\mathit{ref}}^0$ = 0.72/0.87
	e0 = 0.538/0.538	pref = 0.01/0.01 MPa
		ζ = 0.22/0.15
		αe = 0.042/0.01

Table 3.

Model parameters for inner schists and outer greywacke.

The predicted and measured responses of triaxial wetting tests on inner schists and outer greywacke are represented in Figures 11 and 12. It can be found that an overall good agreement is achieved between the experimental observations and the model predictions. The measured stress-strain curves in the dry condition, the transient loss of deviatoric stress associated with flooding and the final recovery of strength in the saturated condition can be well captured with the proposed model. However, the dilatancy is underestimated, especially for samples under the low confining pressures. It is worth to note that the shear dilatancy response in the shearing process after flooding for rockfills of outer greywacke can be reproduced by the proposed model, although the predicted compression volumetric strain is a little larger than that measured. In general, the proposed model can reasonably predict the collapse deformation of the unsaturated rockfill materials subjected to the flooding, although it also has a little shortage to predict the volumetric strain.

Figure 11.

Triaxial tests on compacted schists from inner shell: (a) q vs. ε_a and (b) ε_v vs. ε_a.

Figure 12.

Triaxial tests on compacted greywacke from outer shell: (a) q vs. ε_a and (b) ε_v vs. ε_a.

4. Conclusions

For modelling the strength and deformation behaviour of rockfill materials over various ranges of densities and relative humidity, a state-dependent elastoplastic constitutive model was developed in this study. This model is formulated within an extended critical-state framework by using two independent stress state variables: total stress and total suction. The compressibility model for rockfill materials proposed by Oldecop and Alonso [1] was adopted in this model. The exponential function was used to describe the decrease of the critical stress ratio with the mean stress, while the logarithmic increases of the critical stress ratio and critical void ratio with the total suction were assumed. By substituting these functions into the relation function of peak friction angle with critical friction angle for granular materials proposed by Biarez and Hicher [63], a unified hardening parameter, which could reflect the state-dependent dilatancy, strain-softening and unsaturated effects, was presented for modelling stress-strain relationships of rockfill materials. In addition, the relevant testing procedures have been demonstrated to calibrate and obtain the required model parameters.

The capabilities of the proposed model were illustrated by modelling three series of triaxial tests performed on different rockfill materials (drained triaxial test on saturated samples compacted to various l void ratio, triaxial test with relative humidity controlled and triaxial wetting test). An agreement was obtained between the experimental and numerical results. For the drained triaxial tests on the saturated samples with different initial void ratios, both stress-strain relationship and volumetric strain can be well predicted, which indicates the capability of the proposed model to consider state-dependent behaviour. Through the comparison with the results of relative humidity-controlled triaxial test, the proposed model can be found to well consider the influence of the relative humidity on the stress-strain relationship and volumetric strain of rockfill materials. In addition, the proposed model can well simulate the feature of stress-relaxing according to the comparison with the results of triaxial wetting test. However, the dilatancy is underestimated with the proposed model, especially for the rockfills at the dry state and under the low confining pressures. In general, the model presents an advantage to predict the state-dependent behaviour and collapse deformation of unsaturated rockfills, such as the strain-softening and dilatancy at the dense state. Despite the progress achieved, there is still a room for improvements to be made, for instance, to consider the influence of particle breakage on the changes of compressibility and critical state at different relative humidities.