An elastic-viscoplastic double-yield-surface model for coarse-grained soils considering particle breakage

Computers and Geotechnics 85 (2017) 59–70

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

An elastic-viscoplastic double-yield-surface model for coarse-grained soils considering particle breakage Yufei Kong, Ming Xu, Erxiang Song ⇑ Department of Civil Engineering, Tsinghua University, Beijing 100084, China

a r t i c l e

i n f o

Article history: Received 18 May 2016 Received in revised form 10 December 2016 Accepted 12 December 2016

Keywords: Constitutive model Coarse-grained soil Shear creep Compression creep

a b s t r a c t The proper modelling of time-dependent behavior, such as creep and strain rate effects, of coarse-grained soils (e.g., rockfills) is important to predict the long-term deformation of high-fill projects. In this paper, the primary features of coarse-grained soil are discussed, and an elastic-viscoplastic model is proposed for simulating its time-dependent behaviors in both shear and compression. The stress-dependence of the strength, dilatancy/contractibility and creep is carefully considered based on recent experimental findings. Verifications with various experimental results demonstrate that the model is capable of predicting the creep behavior in various stress states satisfactorily. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction In recent years, high-fill engineering projects have become increasingly common due to the construction of dams or embankments for infrastructures (e.g., airports, highways, etc.) in mountainous regions. The construction of these projects always involves considerable amounts of geo-materials being excavated from high places and filled into valleys, forming a high-fill ground using broken rocks or gravels. For example, the maximum filled height of the Jiuzhai Huanglong Airport is more than 104 m [15]. In contrast to ordinary projects, where stress levels are low and the time-dependent behaviors of filled coarse-grained soils are negligible, the stress level is clearly high in high-fill projects, which may induce a significant instantaneous as well as time-dependent increase in deformation [1,32]. In fact, considerable postconstruction settlements have been observed for some projects [37], which could influence the operation and even cause stability problems [10,21]. Thus, the long-term stability and deformation of high-fill projects must be considered carefully. To describe the time-dependent behaviors of the filling materials and predict the long-term behaviors of high-fills, it is necessary to establish a proper constitutive model for this type of material. To model the deformation behavior of a material, one needs to know the mechanism of the deformation, which is quite complex, especially for the creep of rocks. In modelling the time-dependent ⇑ Corresponding author. E-mail addresses: [email protected] (Y. Kong), mingxu@tsinghua. edu.cn (M. Xu), [email protected] (E. Song). http://dx.doi.org/10.1016/j.compgeo.2016.12.014 0266-352X/Ó 2016 Elsevier Ltd. All rights reserved.

deformation of sedimentary rocks (e.g., shales), Pietruszczak et al. [26] attributed the creep of these rocks to the propagation of micro-cracks as well as the progressive change in their microstructures and proposed a mathematical framework to describe the creep. To model the creep of water saturated chalks, Pietruszczak et al. [27] and Lydzba et al. [16] proposed that the creep of chalks is due to the time-dependent degradation of the intergranular interface, which is caused by the chemo-mechanical interaction, i.e., the dissolution/diffusion within the intergranular interface in the chalks. On that basis they introduced a multiscale framework to describe chalk creep. In these studies, the macroscopic models are established using the mechanism at a microlevel. The philosophy is quite reasonable, but extensive work is still required to determine the quantitative relationship between the macro-phenomenon and the micro-mechanism. Coarse-grained soil consists of broken rock particles, whose short-term and longterm deformations are closely related with the individual grains, and their interactions affected by the environment. The mechanism behind its creep is more complicated. The primary cause of coarse-grained soil creep is commonly recognized as the timedependent breakage of the rockfill grains, especially those with a large size and angular shape [1,17]. There must be also timedependent degradation of inter-grain contacts due to some chemo-mechanical process. Currently, it is still difficult to establish constitutive models intended for the analysis of practical problems directly on the micro-mechanism. Tapias and Alonso [32] and Zhou and Song [51] used DEM (Discrete Element Method) to simulate the time-dependent crack of particles and predict the creep deformation of rockfills. The method proposed may be suitable for

60

Y. Kong et al. / Computers and Geotechnics 85 (2017) 59–70

research but is not yet applicable for practical analysis. Thus, the model here is proposed on the basis of macro-scale creep experiments, and the viscoplastic theory is used as a phenomenological description of creep while considering particle breakage. Based on laboratory tests, it is observed that the mechanical characteristics of coarse-grained soils are prominently different from that of fine grained soils (e.g., clays and silts) due to particle crushing [11,14]. The following 4 features have been observed. (1) Non-linear strength. Its friction angle decreases with the confining pressure. (2) Dilatancy and contractibility. It demonstrates a more pronounced shear contractibility under high pressures and a shear dilatancy under low confining pressures [41]. (3) Isotach behavior. The behavior of this material in drained shearing fits the Isotach hypotheses. A faster loading strain rate corresponds to a higher shear stress and failure strength [2,34]. (4) Stress level dependence of the creep coefficient. It has been discovered through previous studies on crushable sands [6] and recent experiments on soilrock mixtures and rockfills [42] that the creep coefficient (coefficient of secondary compression) for coarse-grained soils varies with the applied load. The first two features have been considered in a few elastoplastic models. For example, Xu and Song [41] and Cao [5] performed modifications to the Mohr-Coulomb failure envelope and Rowe’s dilatancy using these considerations. Wang & Yao [40,48] considered the effects of particle breakage on the timeindependent mechanical behaviors in the UH model with a single yield surface. However, there are few models thus far that have given sufficient consideration to reflect features 3 and 4 in describing the creep behavior of coarse-grained soils. The currently available models simulating the time-dependent behavior of coarse-grained soils are mostly phenomenological, which can be primarily categorized into the following three types. The first type is the single-yield-surface model. These models were originally proposed for soft soils based on oedometer tests or isotropic compression tests, and they are widely accepted nowadays, e.g., the Soft Soil Creep model [36] and the Hunter Clay model [43], for the deformation analysis of soft soils. Some researchers assumed that the creep behavior of coarse-grained soil is similar to that of soft soil. For instance, Wahls (1965) [37] and Parkin [22] once concluded from certain oedometer tests that the compression creep of coarse-grained soils fit a time-logarithm pattern similar to the secondary compression of soft soils. Based on this view, single-yield-surface creep models for soft soils were also applied in the coarse-grained soil calculation [23]. Although the viscoplastic framework is applicable to coarse-grained soils, the application above has at least two defects: (1) The loading rate dependence of the material strength, i.e., feature 3 of the coarsegrained soil mentioned above, cannot be reflected by the existing single-surface models. Indeed, they can simulate the creep and rate-dependent behavior in compression. However, for shear deformation, the drained shear strength simulated is irrelevant to the loading rates because the calculated effective stress of the soil at failure will lie on the critical state line, whose location does not change with time. (2) The stress dependence of the compression creep coefficient, i.e., feature 4 of the coarse-grained soil mentioned above, cannot be reflected. Most of the currently available models use the compression creep law of soft soils, where the creep coefficient is independent of the stress level. This is different from what has been observed from creep tests on coarse-grained soils, as previously mentioned. In the literature, there are a few creep models for soft soils formulated with double yield surfaces [20,38] that can overcome certain disadvantages of the current models and better describe the drained shear creep behavior using a yield surface different from that for compression creep. The second type is the time-hardening model [13,44]. In timehardening models, the creep strain calculated by the model for a

given time merely depends on the current stress state and the time duration from the start to the current moment. Thus, the creep strain is irrelevant to any previous stress states in the time history, i.e., the effects of the loading path and loading history are neglected. The third type is the component model. A component model uses a combination of springs, dampers and plastic elements to compose a mechanical model. Tatsuoka et al. (2006) [33] proposed an Isotach component model for granular materials, in which the strain-stress-time relationship in shearing was simulated using the combination of several components. It is a shear model which can represent the Isotach shear behavior of granular materials well, and is validated by several constant rate loading tests, varying rate tests and creep tests. However, consideration is not sufficiently given for such models to simulate the compression behaviors of soil, such as compression hardening and compression creep. Much work is needed to ensure such models suitably simulate the soil behaviors for both shearing and compression. It can be seen from the above discussion that the currently available creep models are insufficient to describe the primary characteristics of coarse-grained soils, and improvements are required. Some of them considered the strength features but failed to suitably describe creep, and some other models considered creep behaviors but neglected the time effects on drained strength. In constructing a constitutive model for coarse-grained soils, the material’s distinctive characteristics should be included so that better simulations of the deformation for high-fills can be achieved. Moreover, because soils in the high-filled ground may be subjected to various stress states, both shear creep and compression creep can be significant during the lifecycle of the filled body, which can be better described using double yield surfaces. The aim of this paper is to present a practical model for the stress-strain-time behaviors of coarse-grained filling materials, especially rockfills, which are consistent with previous experiences and recent experimental observations. The proposed model is a double-yield-surface hardening model based on the elasticviscoplastic (EVP) framework. The primary features that are different from the previous EVP models consider both shear creep and compression creep, suitability for creep behavior in various stress states, stress dependence and the effects of particle crushing on the material parameters relative to strength and deformation. The model is established using available experimental findings from literature as well as the authors’ group. Even though the underlying physico-chemical mechanisms are not directly considered in the model, the phenomena induced by the time-dependent breakage of grains are properly considered. Then, the proposed model is implemented into a finite element program and validated with various experimental results.

2. Model description The proposed model is based on the elastic-viscoplastic (EVP) framework. In the EVP framework [24,25], the total strain rate tensor fe_ g can be decomposed into an elastic strain rate tensor fe_ e g and a time-dependent viscoplastic strain rate tensor fe_ cr g as follows:

fe_ g ¼ fe_ e g þ fe_ cr g

ð1Þ

where fe_ e g is calculated based on the generalized Hook’s law, and the elastic modulus takes the value of the unloading-reloading modulus [50], which is dependent on the minor principal effective m

0 ref stress Eur ¼ Eref ur ðr3 =p Þ . In the following sections, the direction and magnitude of tensor fe_ cr g will be derived from the strainstress-time relationships observed in relevant experiments.

61

Y. Kong et al. / Computers and Geotechnics 85 (2017) 59–70

2.1. Shear hardening and flow rule The relationship between the vertical strain e1 and the deviatoric stress q from the triaxial compression with a constant confining pressure can be simplified as a hyperbolic curve [8]:

e1 ¼

1 q Ei 1  q=qa

for q < qf

ð2Þ

where q ¼ r01  r03 , and qa is the asymptotic value of q; Ei is the initial stiffness and related to E50 , which is dependent on the confining stress r03 and can be given by the following equations:

Ei ¼

2E50 2  Rf

g 13 g 23

ð11Þ

2.2. Shear creep

c cos u þ r03 sin u c cos u þ pref sin u

m ð4Þ

where pref is a reference pressure; Eref 50 is the corresponding stiffness; and m is a power index reflecting the stress dependence of the stiffness. The asymptotic value of the deviator stress is related with the peak deviatoric stress by the following equation:

qa ¼

ðr1  r2 Þ ðr1 þ r2 Þ  sin wm 2 2 ðr1  r3 Þ ðr1 þ r3 Þ  sin wm ¼ 2 2 ðr2  r3 Þ ðr2 þ r3 Þ  sin wm ¼ 2 2

g 12 ¼

ð3Þ


E50 ¼ Eref 50

where um is the mobilized friction angle, sin um ¼ ðr1  r3 Þ=ðr1 þ r3 Þ; and uv 0 is the mobilized friction angle at the minimum point of the ev  e1 curve and can be obtained from experimental results at different confining pressures. The plastic potential function can be given as follows:

qf Rf

ð5Þ

whereas the peak deviatoric stress qf can be expressed as:

qf ¼ ðc cot u þ r03 Þ

2 sin u 1  sin u

ð6Þ

where Rf is the failure ratio with a range of 0.75–1.00, which is essentially independent of the confining pressure [8]. For coarse-grained soils, c ¼ 0. The hardening law for shear was derived in [30] as follows:

According to Tatsuoka et al. [34] and Anh Dan et al. [2], coarsegrained soils used in high-fill projects can typically be categorized as the Isotach type, which signifies that the location of the q  ea curve measured in a constant-rate loading triaxial compression test is positively correlated to the applied strain rate, and when there is a sudden increase in the loading rate, the transiting q  ea curve will monotonically increase before joining the curve corresponding to the higher loading rate. The schematic relationship for this type of material as well as the three other types (including the TESRA type, P & N type and Combined type) are depicted in Fig. 1. To describe the time-dependent behavior of coarse-grained soil in the EVP framework, the shear behavior of the Isotach type will be expressed using a differential overstress function.

where cp is the plastic shear strain and can be used as the hardening parameter and can be related with the axial plastic strain in triaxial test by: cp ¼ ep1  ep2  ep3  2ep1 . The confining pressure has a significant influence on the peak friction angle u for crushable soils. As more particles become crushed under higher pressures, the friction angle tends to decrease. A non-linear strength envelope can be determined by setting the friction angle as [9]:

2.2.1. Stress-strain curves for constant rate loading For triaxial compression tests with a constant loading rate, the peak strength increases with the axial loading rate e_ a . The deviatoric stress is a unique function of the axial strain and the loading rate, which indicates that different strain rates correspond to different stress-strain curves. Any q  ea curve with a particular constant axial strain rate e_ a is located above the qR  ea curve, which is the q  ea curve when e_ a ! 0. This qR  ea curve is called the static curve in the overstress theory [24] or the reference stress-strain relationship [33], where R is the abbreviation for reference. In the triaxial compression (TC) tests, ea ¼ e1 . The reference stress-strain relations correspond to the relationships described in Section 2.1. The qR  ea relationship (e_ a ! 0) can be expressed by a hyperbolic curve as follows:

u ¼ u0  Du lgðr03 =pref Þ

qR ¼ min

2 q 2q   cp Ei 1  q=qa Eur

ð7Þ

ð8Þ

where u0 is the friction angle when r03 ¼ pref . A non-associated flow rule is assumed for shear yielding. The plastic volumetric strain epv can be related with the plastic shear strain cp in rate form by

e_ pv ¼ sin wm c_ p

ð9Þ

The mobilized dilatancy angle wm is calculated based on Rowe’s stress dilatancy theory [28,29]. However, the particles in coarsegrained soils may have highly irregular shapes and are crushable, which are different from the idealized rods or spheres originally considered in Rowe’s stress dilatancy theory. During shearing, particle breakage, rotation and re-arrangement contribute to the plastic volume change. Coarse-grained soils exhibit shear dilation under low confining pressures and shear contraction under high confining pressures. Xu and Song modified Rowe’s stress dilatancy equation and recommended the following expression of wm for rockfills [41]:

sin um  sin uv 0 sin wm ¼ 1  sin um sin uv 0

ð10Þ




ea

a0 þ b0 ea



ð12Þ

; qRf

Deviatoric stress q

fs ¼

Isotach Combined TESRA P&N

Continuous constant rate loading

Axial Strain εa Fig. 1. Schematic stress-strain curves for four types of viscosity in triaxial compression.

62

Y. Kong et al. / Computers and Geotechnics 85 (2017) 59–70

where a0 ¼ 1=ERi ; b0 ¼ 1=qRa ¼ Rf =qRf ; qRf ¼ 2 sin uR r03 =ð1  sin uR Þ; ERi is the initial stiffness; qRa is the asymptotic value of qR ; and uR is the peak friction angle for the reference curve, which is independent of the loading rate. Furthermore, the q  ea curves for a constant rate loading (e_ a > 0) can be expressed in a similar form:


q ¼ min

ea

a þ bea

 ; qf

ð13Þ

where qf is the peak stress, which increases with the loading rate e_ a ; and b ¼ Rf =qf , where Rf is assumed to be independent of the strain rate. However, the change in the a, b and qf parameters with the strain rate e_ a must still be determined. A few relationships can be concluded from experiments as follows:

2.2.2. Differential law for shear creep At a given time, if the deviatoric stress is no higher than the reference relationship (q 6 qR ) depicted in Fig. 2, i.e., if the current stress lies inside the shear yield surface, then the shear yield surface will not change with time, and no creep will occur. If the deviatoric stress is greater than the reference relationship (q > qR ), then shear creep develops, and e_ cr a – 0. In a strain-hardening model, if the stress and strain states of the soil (i.e., q and cp ) at a given time are known, then the axial viscoplastic strain rate in triaxial shearing, e_ cr a , is a unique function. The function that uses the current stress and strain states as input variables will be explained in this section and can be derived from the constant rate triaxial compression tests. p p Based on the relationships in Section 2.1, ecr a ¼ ea ¼ c =2. Thus, the total axial strain can be expressed by q and cp as

ea ¼ ecra þ eea ¼

1. Several equations are recommended by Perzyna [24] for describing the rate dependence of strength, including

"

cp 2

þ

q Eur

ð21Þ

Using Eq. (20), we can obtain




  q  qR  1 fqR


f # e_ a R qf ¼ qf 1 þ  e_

ð14Þ


 e_ a qf ¼ qRf 1 þ  e_

ð15Þ


 e_ a qf ¼ qRf 1 þ f ln 1 þ  e_

where qR can be calculated from Eq. (12), and q and qR correspond to the same ea . For the constant rate loading process in Fig. 2, the strain rate e_ a includes an elastic strain rate e_ ea ¼ ð1=Eur Þðdq=dtÞ and a viscoplastic strain rate e_ cr a , which can be further expressed as

ð16Þ

e_ cra ¼ e_ a 

where e_  and f are the material constants. These equations are used to fit the results of experiments performed by [2], as indicated in Fig. 3. Eqs. (14) and (16) demonstrate a reasonably good agreement with the experimental results. Eq. (16) is used in the proposed model. 2. By substituting the qf  e_ a relationship into b ¼ Rf =qf , the b  e_ a relationship can be derived. 3. When considering the rate dependence of a, the assumption of a / b is made, which indicates a fair accuracy in fitting the experimental data from [2]. Based on the discussions above, the equations can be used to express the parameters a, b and qf for loading with a constant strain rate e_ a as follows:




 1 dq 1 dq dea 1 dq ¼ e_ a  ¼ e_ a 1  Eur dt Eur dea dt Eur dea

ð23Þ


 
 dq e_ a dqR e_ a ¼ 1 þ f ln 1 þ  ¼ 1 þ f ln 1 þ  ERt dea e_ dea e_

ð24Þ

where ERt can be calculated using

(

ERt

¼

a0 ða0 þb0 ea Þ2

; for qR < qRf

ð25Þ

for qR ¼ qRf

0;

From Eqs. (20) and (24), we know that dq=dea ¼ ðq=qR ÞERt . By substituting this equation and Eq. (22) into Eq. (23), we can obtain



e_ cra ¼ e_  exp

!
  q  qR ERt q  1 1  fqR Eur qR

ð26Þ

ð17Þ

b¼

Rf Rf b0 ¼    ¼   qf qRf 1 þ f ln 1 þ ee__ a 1 þ f ln 1 þ ee__ a

ð18Þ

a¼

a0 a0 b¼   b0 1 þ f ln 1 þ ee__ a

ð19Þ

By substituting (12), (17)–(19) into (13), the deviatoric stress q can be expressed by a unique function of the instantaneous strain ea and strain rate e_ a by


 e_ a q ¼ qR 1 þ f ln 1 þ  e_

ð22Þ

According to Eq. (20), dq=dea is proportional to dqR =dea or ERt , and the tangential modulus on the reference curve can be given by the following equation:

Deviatoric stress q


 e_ a qf ¼ qRf 1 þ f ln 1 þ  e_

e_ a ¼ e_  exp

q f | a3

a1

a2

a3 a3

Strain rate increase Strain rate decrease

a2 a1

q Rf

a

0

Reference Stress-strain relation

ð20Þ

where qR is the stress on the reference stress-strain curve, as expressed by Eq. (12); and f and e_  are the material constants reflecting the extent of the rate effect.

Axial Strain Fig. 2. Schematic stress-strain curves for the Isotach type.

a

63

Y. Kong et al. / Computers and Geotechnics 85 (2017) 59–70

(a) 2500 Experimental data of Chiba grave (old batch)

2450 2400

10

2000

Eq. (14) Eq. (15)

2350

10

Eq. (16)

q (kPa)

Failure deviatoric stress q f (kPa)

2500

2300 2250

0.1

0

0

0

0

1500 0

1000 0.1

2200

Experiment of Chiba gravel (old batch) Simulating CRL 10 0

0

Simulating CRL

500

2150

0.006%/min

0

Simulating CRL 0.1

0

0

Simulating Changing Rate Loading 0

2100

0

2

4

2050 10 -8

10 -7

10 -6

10 -5

Rate of axial loading

10-4 a(



ð27Þ

D=1.1 D=1.05 D=1 D=0.9 D=0.8 D=0.4

(%)

8

6

4

2

0

0

200

400

Therefore, when q 6 q , the rate of shear creep will be zero.

800

1000

(c)

0.108

0.106

(%)

v

'

1MPa

0.104

v

2.2.3. Extension to 3-D conditions By referencing the generalized form of the functions recommended by Zienkiewicz and Cormeau [52] and the method for generalizing constitutive models by Yao et al. [45–47], the differential law for shear creep in the drained TC tests can be extended to 3-D conditions. The isotropic assumption can be used in the extension as follows:

600

t (min)

R

0 3;

12

a

To consider the complex conditions, q  q is replaced by hq  qR i, which signifies that hxi ¼ x for x > 0 and hxi ¼ 0 for x 6 0. Then, Eq. (26) is now written as R

!


hqqR i  ERt q 1 @g s cr  R fq _ _ fes g ¼ e e 1 1 0 R Eur q @g s =@ r1 @ r0

10

(b) 10

/min )

!
  hq  qR i ERt q  _ _ecr 1 1 a ¼ e exp fqR Eur qR

8

(%)

10 -3

Fig. 3. Comparison of the three equations in fitting the experimental data for the qf  e_ a relationship.

0 1

6 a

0.102 Simulating creep behavior in oedometer apparatus with the shear creep model

ð28Þ 0.1

f _ cr s g

where q ¼ r  r e is a vector denoting the shear creep strain rate; fr0 g is a vector representing the current effective stress; e_  and R

ERt

f are the material parameters; and q and are calculated using Eqs. (12) and (25) respectively. Because the plastic deformation in this model is a time-dependent viscoplastic deformation, the plastic potential surface is the same as the viscoplastic potential surface. The expression for the plastic (viscoplastic) potential function is defined in Eq. (11). 2.2.4. Discussion on the shear creep model Eq. (28) describes the law of shear creep that is validated with a drained triaxial compression test on Chiba gravel. The test was conducted by Anh Dan et al. [2] under a constant confining stress r03 ¼ 490 kPa. The stress-strain relations of the constant strain rate loading and the varying rate loading are simulated and presented in u0 = 43.53°, Du = 5°, Fig. 4a. The parameters include Eref m = 0.5, tur = 0.2, uv 0 = 40°, Eref ur = 1370 MPa, 50 = 176 MPa, _e = 9.58  106/min, f = 0.2364, Rf = 0.87, and pref = 0.49 MPa. The deviatoric stress level D is defined as the ratio of the deviatoric stress to the reference peak deviatoric stress qRf and can be larger than 1.0. For the Chiba gravel test, the reference failure strength is qRf = 2168 kPa, and the applied maximum deviatoric stress level is Dmax ¼ qmax =qRf = 1.10. The model can suitably simulate the effects

0.098 -1 10

0

10

1

10

2

10

3

10

t (min) Fig. 4. (a) Numerical simulation of the triaxial compression test on Chiba gravel (CRL is the abbr. for constant rate loading, and the reference rate e_ 0 ¼ 0:006%=min); (b) simulation of creep behavior under different stress levels; and (c) model with a singular shear yield surface used to simulate an oedometer creep test lasting 1000 min, in which the creep rate reduces to zero in 30 min, thus differing from experimental observations.

of the rate on the behavior of this material under a high stress level. In addition, the proposed model is capable of simulating transient creep and stationary creep for the condition in which the shear stress is higher than the long-term shear strength. By assuming that the shear stress q increases to a certain level instantly and then remains constant, the creep process can be calculated using the model, as indicated in Fig. 4b. The shear stress levels considered include q ¼ 867:6 kPa (D ¼ q=qRf ¼ 0:4), q ¼ 1735:2 kPa ðD ¼ 0:8Þ, q ¼ 1952:1 kPa ðD ¼ 0:9Þ, q ¼ 2169 kPa ðD ¼ 1Þ, q ¼ 2277:5 kPa ðD ¼ 1:05Þ, and q ¼ 2385:9 kPa ðD ¼ 1:1Þ. As D becomes smaller, the shear creep converges faster.

64

Y. Kong et al. / Computers and Geotechnics 85 (2017) 59–70

Then, this set of parameters is used to predict the compressive creep in an oedometer test, where the vertical load is r0v = 1 MPa and applied at t ¼ 0. The vertical strain ev versus time t is plotted in Fig. 4c. In the simulation of a 1000 min period, the predicted creep rate reaches zero in 30 min. This result is notably inconsistent from the phenomenon commonly observed in oedometer compression tests using coarse-grained materials, in which the vertical strain increases with log time in the long term. It is revealed that the model, when only considering shear creep, will fail in predicting the creep behavior of soil under low deviatoric stress levels, especially the volumetric creep under isotropic compression or oedometer compression. This limitation indicates the necessity of introducing a cap yield surface in the model. 2.3. Cap hardening and flow rule The quarter elliptical cap yield surface suggested by Schanz et al. [30] is used here to describe the compression hardening, with the associated flow rule. Thus, the plastic potential function and the cap yield surface can be given as

gc ¼ f c ¼

~2 q

a2

þ p02  ðpp Þ

2

ð29Þ ~ ¼ r01 þ ðd  1Þr02  dr03 ; q qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi nc a ¼ 3 ð1  K nc 0 Þ=ð1 þ 2K 0 Þ

p0 ¼ ðr01 þ r02 þ r03 Þ=3;

where

d ¼ ð3 þ sin uÞ=ð3  sin uÞ;

and

(default: K nc 0 ¼ 1  sin u). The magnitude of the yield cap is determined by the isotropic pre-consolidation pressure pp . As is shown in Fig. 5, pp indicates the intersection of the cap yield surface and the p-axis. For a given vertical effective stress r0v , a contour surface of the cap yield function crosses the stress point for the current state nc 0 0 ~ (p0 ¼ ð1 þ 2K nc 0 Þrv =3, q ¼ ð1  K 0 Þrv ), which is known as the equivalent surface. The magnitude of the equivalent surface can be denoted by peq as follows:

peq ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~2 q þ p02 ¼ vr0v 2

where

a

ð30Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi nc 2 ð1 þ 2K nc 0 Þ=9 þ ð1 þ 2K 0 Þ =9. Similarly, the isotropic

v¼

pre-consolidation pressure pp can be related to the vertical preconsolidation pressure in oedometer tests r0p v by

pp ¼ vr0p v

ð31Þ

The compressive stiffness for the oedometer tests Eoed can be expressed by the vertical effective stress r0v ð¼ r01 Þ or confining pressure r03 and the oedometer stiffness at the reference pressure as follows: Eref oed




Eoed ¼ Eref oed

r0v

pref

m


¼ Eref oed

0 3 ref K nc p 0

r

m ð32Þ

For silicon sands, it is reported that 0:5 < m < 1 [3]. However, for soils with crushable grains, 0 < m < 1 (e.g., after fitting results from isotropic compression tests on rockfills, m ¼ 0:4 [31]). The relationship between the hardening parameter pp and the volumetric cap strain epc v can be derived from the above equations as follows:

epc v ¼


p 1m b p þC 1  m pref ref




where b ¼ vp1m

p ' (1 2K 0nc) v '/ 3 q (1 K 0nc ) v '

pp v

'

v

'p

ne K0-li Cap Yield Surface

Equivalent Surface

p eq

v

p ' p

Fig. 5. Current stress, yield surface and equivalent surface.



e_ pc v ¼b

ð1þmur Þð12mur Þ

pp pref

m

p_ p pref

ð33Þ

 .

m ref Eur ð1mur ÞðK nc 0 Þ

2.4. Compression creep Compression creep tests typically refer to two types: (1) oedometer compression tests with a vertical stress rv ¼ const and (2) isotropic compression tests with a constant confining pressure. Because the oedometer (i.e., one-dimensional) stress state is more common in large-area high-fill projects, the results of the oedometer compression tests are investigated here. In this section, the experimental observations are summarized from literature as well as some tests conducted by the authors. 2.4.1. Creep behavior with a constant vertical stress in oedometer tests It has been widely accepted that a logarithmic relationship is practical to describe the long-term volumetric change of clays in creep tests [35]. There is evidence that the compression creep of coarse grained soils under constant load also follows such a rule [18]:

ev ¼ C e lnðt þ tc Þ þ const

ð34Þ

where C e is the coefficient of secondary compression [36], and in oedometer tests, it is denoted as C e;oed . The coefficient C e is essentially stress-independent [7] for clays. However, experimental studies reveal that C e is stress-dependent for coarse-grained soils [19,42]. It should be noted that because the permeability of coarsegrained soils is high, for a saturated sample in the oedometer test, the excess pore pressure dissipates quickly, and the consolidation time is negligible. Thus, the applied load immediately induces an increase in the corresponding effective stress. A conceptual graph describing the deformation increasing with time is illustrated in Fig. 6, where s is an artificially assigned time (generally a unit time, e.g., 1 h or 1 d [49]), at which the state of the soil is regarded as normally consolidated. With the normallyconsolidated state strains (ev ;NC ) under different vertical loads and m can be calculated by fit(r0v ), the material parameters Eref oed ting the stress-strain data with the equation describing their rela ref m p 1 . tionship. According to Eq. (32), we know E 1 ¼ ref r0v oed E oed

ref m de dev ;NC p 1 Meanwhile, E 1 ¼ dvr;NC ¼ ref . By integrating this 0 . Hence, dr0 r0 oed

q

1 ref E oed




or

v

v

E

oed

v

equation from zero to the current value of ev ;NC  r0v relationship as follows:

r0v , we can obtain the

m

ev ;NC ¼

ðpref Þ ðr0v Þ1m ð1  mÞEref oed

ð35Þ

2.4.2. Stress dependence of compression creep Several large oedometer tests were performed on the rock-soil mixture from Chongqing Airport, Chongqing, China, and the results are presented in Fig. 7. The samples have a diameter D ¼ 200 mm and a height h ¼ 190 mm. The vertical load was increased gradually. After each load increment, the load was kept constant for 12 h. The loading process was 0 MPa ? 0.34 MPa

65

Y. Kong et al. / Computers and Geotechnics 85 (2017) 59–70

ε

1

ε

1

Cε

Cε

1

τ

tc

tc τ

ln(t )

t

Fig. 6. Conceptual strain-time relationship under a constant effective stress (after Vermeer and Neher [36]), in which the loading time is t ¼ 0.




0.08

C e ¼ C e;oed ¼ C ref e;oed

0.07

Δε v (%)

0.04 0.03

0.01

2

10

5

t (h), log scale Fig. 7. Vertical creep strain versus time for each loading step, where creep is assumed to start at t ¼ 1 h, i.e., Dev ¼ ev ðtÞ  ev ð1 hÞ.

(12 h) ? 0.54 MPa (12 h) ? 0.74 MPa (12 h) ? 0.94 MPa (12 h). Because the slopes of the curves are observed to become nearly constant after 1 h of loading, s is assigned to be 1 h. The coefficient C e corresponding to each r0v is calculated by fitting the data for t > s. It is observed from this experiment that C e is positively correlated with r0v . The value of C e is plotted with r0v in Fig. 8a in addition to the results from other tests performed by the authors. Furthermore, the stress dependence of C e is supported by previous 1-D compression tests on coarse-grained soils [42]. Based on the test results, an empirical relationship for C e and r0v in oedometer tests can be proposed as

(b)

4.5 4

Chongqing Sample 1

3.5

Chongqing Sample 2 Kunming Sample

0

×1

⎞

ref K nc 0 p

2.5 2

0.7569

⎛ σ v ⎞⎟ 0.360 × ⎜⎝ pref ⎠

1.5 1

8 22 0.6

70

v ⎛ σ ref ⎟ ⎠ ×⎜ p 41 -4 3 ⎝ 3 0.75 0.6 × 10 ⎞ v σ ⎛ ref ⎟ 8× ⎜ p ⎠ ⎝ 0.56

3

Cε

60

-4

× 10

p

9 0.547

Cε

50

p

40 30 Calcareous Sand(dense) Calcareous Sand(loose)

20

0.5 0

mc

80

-4

7 27 0.7

Cε (×10-4 )

-4

Cε (×10 )

(a)

r03

2.4.3. Differential law for 1-D creep Eqs. (34) and (36) provide the creep laws for granular soils under 1-D conditions with a constant vertical effective stress, from which the differential equation for 1-D creep is formulated in this section to implement this law in a numerical model. In oedometer compression, the state of the soil can be represented by (i) its vertical effective pressure r0v and vertical preeq consolidation pressure r0p and the v or (ii) the equivalent stress p isotropic pre-consolidation stress pp . Their relationships are expressed in Eqs. (30) and (31). p With an increase in the plastic strain, r0p v and p increase with time. The time-hardening process is illustrated in Fig. 9. The pre-loading stress state is presented in Fig. 9a, which is located inside the yield surface. A stress increment is then applied, as indicated in Fig. 9b. In an EVP framework [4], the material behaves elastically during a sudden application of the stress increment, and the instant elastic stress-strain relationship is governed by Hooke’s law. The stress state is brought beyond the cap yield surface by the loading, which is not allowed in traditional elastoplastic mechanics but can be achieved in the EVP framework. As indicated in Fig. 6, the strain continues developing when t > 0, although the strain rate decreases with time. The expanding cap yield surface is located inside the equivalent surface in Fig. 9c,

0.02

1

pref


or C e ¼ C e;oed ¼ C ref e;oed

where mc is a material parameter. This power function relationship can be also approximately applied to other loading conditions in addition to oedometer compression. For example, in isotropic compression tests, C e / pmc (see Fig. 8b) [6].

0.05

0

mc

ð36Þ

0.94MPa 0.74MPa 0.54MPa 0.34MPa

0.06

r0v

10 0

0.2

0.4

0.6

0.8

σ v (MPa)

1

1.2

1.4

0

0

5

10

15

Confining pressure p (MPa)

Fig. 8. (a) C e -rv relationship in the oedometer tests (pref ¼ 0:1 MPa); (b) C e -p relationship in the isotropic compression tests on Calcareous Sand, which is derived from data in [6].

66

Y. Kong et al. / Computers and Geotechnics 85 (2017) 59–70

(a)

q

q

Pre-loading state

q (c) 0 < t <τ OCR < 1

(b) State at the time of loading t=0 σ v 'p

σv '

σ v0 '

Cap yield surface keeps expanding

cap yield surface

p

(d)

q

q

t =τ

OCR = 1

pp

(e)

p

p eq = χσ v '

t >τ OCR > 1

pp

p

p eq

Cap yield surface

σv '

p

Equivalent surface

p

p p = p eq

Current state

p

p eq p p

~ plane (q ~ ¼ r01 þ ðd  1Þr02  dr03 ) at different times: (a) the stress state before loading; (b) the stress is increased after a Fig. 9. Position of cap yield surfaces on the p  q stepwise loading; (c) the cap yield surface expands with time, and the plastic deformation increases; (d) when t ¼ s, the sample reaches a normally-consolidated state (OCR ¼ 1); and (e) under a constant effective stress, the creep and hardening process continue.

if 0 < t < s. When the loading time t ¼ s (Fig. 9d), the stress point 0 p eq is located on the yield surface, and r0p v ¼ rv , p ¼ p , and OCR ¼ 1. At this time, the material is considered to be normally consolidated. For t > s, the cap yield surface moves beyond the current stress state with a decreasing rate of hardening (Fig. 9e). It is assumed that the initial vertical effective stress is r0v0 , and the initial vertical pre-consolidation pressure is r0p v0 . The vertical effective stress is increased to r0v during loading. The vertical strain e and the vertical pre-consolidation pressure r0pv at any given time t > 0 can be derived in the following. According to Eq. (35), ev  ðr0v Þ1m has a linear relationship in oedometer compression, and the slope for this line is B ¼

ðpref Þ

m

ð1mÞE

ref oed

,

as indicated in Fig. 10. For unloading conditions, the unloading modulus is Eur , which is also stress-dependent, i.e., ref m

0 Eur ¼ Eref ur ðr3 =p Þ . It should be noted that Eur is Young’s modulus 0 instead of the oedometer modulus. By assuming r03  K nc 0 rv in the unloading process, the oedometer unloading modulus can be esti 0 m ref r3 mur 1mur mated as Eur;oed ¼ Eur ð1þmur1Þð12 . Thus, mur Þ ¼ Eur pref ð1þmur Þð12mur Þ

ev  ðr0v Þ1m has a linear relationship in the oedometer unloading process with the slope A ¼ ðp

ref Þm ð1þm Þð12m Þ ur ur ref ð1mÞEur ð1mur Þ

.

Based on the expressions of the loading/unloading modulus, the vertical strain at t can be expressed using the stresses and state parameters (see Fig. 10) as follows:

e ¼ ee þ ecr

ð37Þ

ee ¼ A½ðr0v Þ1m  ðr0v0 Þ1m 

ð38Þ

ecr ¼ ðB  AÞ½ðr0pv Þ1m  ðr0pv0 Þ1m 

ð39Þ

e

where e is time-independent during the creep process; B is the slope for the normal consolidation line (NCL); and A is the slope for the unloading-reloading relationship in the oedometer tests. 0 As a special case of Eq. (37), the strain for t ¼ s (r0p v ¼ rv ) can be given by

es ¼ ejt¼s ¼ ee þ ecr jt¼s

ð40Þ

ecr jt¼s ¼ ðB  AÞ½ðr0v Þ1m  ðr0pv0 Þ1m 

ð41Þ

On the other hand, the strain for time t is

e ¼ C e lnðt=sÞ þ es

Therefore, based on Eqs. (39), (40) and (42), the strain increment between s and t (t > t c ) can be expressed in two forms as follows:

e  es ¼ ðB  AÞ½ðr0pv Þ1m  ðr0v Þ1m  ¼ C e lnðt=sÞ

σ v0 ' σ v0 '

σv '

σv '

p

σ

Eq. (43) implies the relationship between the hardening parameter

From Eq. (43), we can obtain the relationships as 1m

1− m v

t

s

¼ exp


 Ce t 1m þ ðr0v Þ ln BA s

ð44Þ

 B  A 0p 1m 1m ½ðrv Þ  ðr0v Þ  Ce

ð45Þ

¼

Under a constant effective vertical pressure 0p 1m vÞ

time. The time derivative of ðr

ε εv

ε

1

cr

εe

1

A

0p 1m vÞ

B

Cε ln(t / τ )

Fig. 10. Idealized stress-strain curve for the oedometer test.

ð43Þ

r0pv and time t (t > tc ). r0v0 and r0pv0 do not appear in this equation. Based on Eq. (36), the current stress state r0v has an effect on C e . ðr0p vÞ

p

ð42Þ

dðr dt

¼

r0v , r0pv varies with

can be given by

Ce 1 Ce ¼ B  A t sðB  AÞ

 B  A 0p 1m 1m  exp ½ðrv Þ  ðr0v Þ  Ce

ð46Þ

Then, from Eqs. (41) and (46), we can obtain the formula for the plastic strain rate as

67

Y. Kong et al. / Computers and Geotechnics 85 (2017) 59–70 1m

decr dðr0p vÞ ¼ ðB  AÞ dt dt

 Ce B  A 0p 1m 1m exp ½ðrv Þ  ðr0v Þ  ¼ Ce s

ð47Þ

2.4.4. 3-D differential equation for oedometer compression For the analysis of practical problems, it is necessary to formulate a differential law for general stress states. With reference to the methods in [36], the notations in Eq. (47) are changed from r0v and r0pv to peq and pp , respectively. Then, a differential equation for oedometer compression can be formulated in vector form as follows:

f _ cr c g

e

"

1m #) 1m pp b peq ¼ exp  C e ð1  mÞ pref s pref

 1 @g c  P3 0 @ r0 i¼1 @g c =@ ri Ce

Eref oed = 50 MPa,

Eref ur = 320 MPa,

m = 0.5,

0.15

0.1 D=0.5 Experiment D=0.75 Experiment D=0.5 Simulation D=0.75 Simulation Half of line

0.05

0

20

40

60

80

100

120

Creep time (h)

ð48Þ

2.4.5. Discussion on the compression creep model As the compression creep model expressed by Eq. (48) is derived from oedometer compression, it is capable of simulating creep processes in oedometer tests and triaxial tests that are at approximately the same deviatoric stress level. In oedometer compression, r03  ð1  sin uÞr01 , q  r01  r03 ¼ 0 ðr3 sin uÞ=ð1  sin uÞ ¼ 0:5qf , and D  0:5. The single cap yield surface model is validated using triaxial compression creep tests conducted by Wang and Yin [39]. The parameters include C ref e;oed

0.2

0

(

where fe_ cr c g is a vector denoting the plastic strain rate due to cap hardening; and C e is the coefficient of secondary compression determined using the oedometer tests, as expressed in Eq. (36). The viscoplastic potential function is the same as the plastic potential function expressed by Eq. (29).

u = 40°,

Volumetric creep strain ε vcr (%)

0.25

e_ cr ¼

tur = 0.2,

4

¼ 1:66  10 , mc = 0.5, and pref = 0.1 MPa. In the experiment, the volumetric creep strain rate for D ¼ 0:75 is approximately half of that for D ¼ 0:5. However, if the equation for the oedometer test is used directly for D ¼ 0:75, it would predict nearly the same ev  t relationship as D ¼ 0:5 (see Fig. 11). Thus, the cap model described by Eq. (48) needs a correction when considering the deviatoric stress level. Moreover, as two important components of the double surface model, the quarter elliptic yield surface and the shear yield surface cannot be separated from each other to form a suitable creep model. If only the quarter elliptic yield surface is used to simulate shear behavior, the shear strain calculated would be too small to reflect the true deformation of a soil. 2.5. Formulation of the hardening soil creep model It can be concluded from the above analysis that to model the creep behaviors under different stress states properly, it is necessary to involve both shear and compression creep mechanisms in the constitutive model. However, for the condition where D ! 0, the rate of shear creep is negligible in the total creep strain, and the compression creep primarily contributes to the creep. Conversely, when D is increased, the effects of the time-dependent hardening of the cap surface will decrease as C e decreases. In the double-yield-surface model proposed here, the total strain rate is the sum of the elastic strain rate and the creep strain rate fe_ g ¼ fe_ e g þ fe_ cr g. By considering the deviatoric stress level, the creep strain rate fe_ cr g can be expressed as a combination of

Fig. 11. Simulation of the drained triaxial compression creep tests at different stress levels (r03 ¼ 200 kPa) using the compression creep model expressed by Eq. (48).

the shear yield surface creep and the cap yield surface creep as follows:


hqqR i   @g  ER 1 s fe_ cr g ¼ De_  e fqR  1 1  Eurt qqR @gs =@ r01 @ r0 þ

 

1m  eq 1m p C P3  prefp ð1  DÞ es;iso exp ð1DÞC b ð1mÞ ppref e;iso

i¼1

1 @g c =@ r0i

@gc  @ r0

;

for 0 6 D 6 1
hqqR i   @g  ER 1 s fe_ cr g ¼ e_  e fqR  1 1  Eurt qqR @gs =@ r0 @ r0 ; for D > 1 1

ð49Þ Details regarding the notations and their relationships with the model parameters have been explained in previous sections. For D ¼ 0:5, the volumetric strain rate over the long-term period described in the equation above is equal to Eq. (48). Thus, ð1  0:5ÞC e;iso ¼ C e;oed and C e;iso ¼ 2C e;oed . For D P 1, this equation is equivalent to Eq. (28). 2.6. Parameter determination There are 13 parameters in the proposed model, including: u0 , friction angle at the reference confining pressure; Du, the parameter reflecting the effects of the confining pressure on the friction ref ref angle; Rf , the failure ratio; Eref 50 , Eoed and Eur , the stiffness parameters; m, the power index for the stress dependence of stiffness; tur , Poisson’s ratio for unloading/reloading; uv 0 , the mobilized friction angle for dilatancy/contraction; e_  and f, the parameters

related to shear creep; C ref e;oed , the coefficient of secondary compression for the oedometer creep test at the reference vertical pressure; and mc , the power index for the stress dependence of C ref e;oed . _ Among the 13 parameters, u0 , Du, Rf , Eref 50 , m, uv 0 , e and f can be determined using constant rate loading triaxial compression tests with different confining pressures and loading rates or triaxial creep tests. Eref ur and Eref oed ,

tur can be determined using unloading tests.

C ref e;oed

and mc can be determined using oedometer creep tests. To determine the shear creep parameters e_  and f, several constant rate shearing triaxial tests are needed to obtain a group of curves, as indicated in Fig. 2. By fitting the peak strength qf from the curves and the corresponding shearing rate e_ a to the relation   ship qf ¼ qRf 1 þ f ln 1 þ ee__ a , e_  and f as well as the reference strength qRf can be evaluated for the given confining pressure. With a group of curves under a given confining pressure r3 , it is not dif-

68

Y. Kong et al. / Computers and Geotechnics 85 (2017) 59–70

ficult to calculate Rf , u and E50 corresponding to the given r3 . After the completion of the loading process, an unloading test can be conducted to measure Eur and tur . For coarse-grained soils, u, E50 and Eur are stress-dependent. By conducting the above experiments under several different confining pressures, their values can be obtained for different confining pressures. Then, by fitting these values to the already known relationship formulas, the

3. Validation of the model The proposed model has been implemented in the finite element program PLAXIS, and validations are performed by comparing the calculated results with the experimental data. 3.1. Triaxial creep test of a rockfill material

ref parameters u0 , Du, Eref 50 , Eur and m, contained in these formulas, can be obtained. The parameter uv 0 is the mobilized friction angle at the minimum point of the ev  ea curve. It is the mobilized friction angle corresponding to the dividing point between shear dilation and contraction during triaxial shear. To obtain this value, we should first determine ea corresponding to the minimum point of the recorded ev  ea curve, and the corresponding qR value can be calculated based on the qR  ea curve. Then, uv 0 can be obtained from qR ¼ 2 sin uv 0 r03 =ð1  sin uv 0 Þ.

The drained triaxial compression creep tests conducted by Wang and Yin [39] are simulated here. In these tests, the confining pressure r3 = 200 kPa was first applied, and r1 was then increased over a short amount of time until the deviatoric stress level reached 0.5 or 0.75, which was followed by the creep process under the applied load. The strain increment starting from 0.5 h after the load was applied was defined as the creep strain in the experiments. The numerical simulation follows the same loading process as that in

0.3

the tests. The parameters include u = 42°, Rf = 0.9, Eref 50 = 40 MPa, Eref m = 0.5, tur = 0.2, uv 0 = 35°, Eref ur = 320 MPa, oed = 50 MPa, 4 _e = 1.5  106/min, f = 0.03, C ref e;oed = 1.66  10 , mc = 0.5, and

(b) Volumetric creep strain ε vcr (%)

0.35

cr

(a) Axial creep strain ε a (%)

ref The parameters Eref oed , C e;oed and mc should be obtained from the staged loading oedometer tests. Their determination has been provided in Section 2.4.

0.25 0.2 0.15 0.1

D=0.5 Experiment D=0.75 Experiment D=0.5 Simulation D=0.75 Simulation

0.05 0

0

20

40

60

80

100

0.25

0.2

0.15

0.1

0

120

D=0.5 Experiment D=0.75 Experiment D=0.5 Simulation D=0.75 Simulation

0.05

0

20

40

Creep time (h)

60

80

100

120

Creep time (h)

Fig. 12. Simulation of the drained triaxial compression creep tests of the rockfill material at different stress levels (r03 ¼ 200 kPa) (a) axial creep strain ecr a versus creep time; (b) volumetric creep strain ecr v versus creep time.

(a)

(b)

2 0.94MPa Experiment 0.94MPa Simulation

0.94MPa Experiment 0.94MPa Simulation 0.06 0.74MPa Experiment 0.74MPa Simulation

0.74MPa Experiment 0.74MPa Simulation

0.5

0.04

v

1

Δ ε (%)

v

Δ ε (%)

1.5

0.02 0.54MPa Experiment 0.54MPa Simulation

0

0.54MPa Experiment 0.54MPa Simulation

0.34MPa Experiment 0.34MPa Simulation

-0.02 0

0

2

4

6

t (h)

8

10

12

0.5

1

2

5

10

t (h),log scale

Fig. 13. Vertical strain versus time in the experiment and simulation: (a) incremental strain after the first loading step, Dev ¼ ev ðtÞ  ev ;0:34 MPa ; (b) creep strain of each loading step, Dev ¼ ev ðtÞ  ev ð1 hÞ.

69

Y. Kong et al. / Computers and Geotechnics 85 (2017) 59–70

400

(b)

Sample after isotropic drained creep for 180 min 300

Sample after isotropic drained creep for 3 min

200

100

Drained creep at q=200kPa for 24h

Experiment

cr cr Creep strain ε 1 − ε 2 (%)

Deviatoric stress q (kPa)

(a)

0.5

Sample after isotropic drained creep for 180 min

0.4

0.3

0.2

Sample after isotropic drained creep for 3 min

0.1

Experiment

Simulation 0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Simulation 0

0

5

p p Plastic strain ε1 − ε 2 (%)

8

10

15

20

25

Creep time (h)

Fig. 14. Simulation of the loading-history-dependence experiment (a) q versus plastic strain

pref = 0.1 MPa. The relationships between the axial creep strain and the volumetric strain versus time are presented in Fig. 12.

ep1  ep2 ; (b) creep strain ecr1  ecr2 versus creep time.

and the same stress states. However, the shear strains are different, which demonstrates that the effects of the loading history can be properly captured by the proposed model.

3.2. Oedometer creep test of rock-soil mixture The staged loading oedometer creep tests on the rock-soil mixture from the Chongqing Airport, performed by the authors, are simulated using the proposed model. By taking the strain at the end of the first loading step (0.34 MPa) ev ;0:34 MPa as the initial value, the incremental strains versus time for the other three higher load levels, i.e., 0.54 MPa, 0.74 MPa and 0.94 MPa, are plotted in Fig. 13a. The creep strain-time relationships for the four load levels are presented in Fig. 13b, where Dev ¼ ev ðtÞ  ev ð1 hÞ represents a relative value to the strains at t ¼ 1 h. It can be seen that the simulation fits the experiment results well and reflects the stress dependence of C e . The parameters are u = 40°, Du = 8°, Rf = 0.9, Eref 50 = 16.6 MPa,

Eref = 20.8 MPa, oed _

Eref ur = 175 MPa,

m = 0.247,

6

5 tur = 0.2, uv 0 = 35°, e = 2  10 /min, f = 0.02, C ref e;oed = 5.86  10 ,

mc = 0.754, and pref = 0.1 MPa. 3.3. Loading-history test of loose sand A set of experiments on loose sand was conducted by Kiyota et al. [12], which demonstrated the double yielding characteristics as well as the loading-history-dependence of coarse-grained soils. The model parameters were first calibrated using the triaxial compression test results under the confining pressures of 100 kPa, 200 kPa, and 400 kPa and then used to simulate two sets of triaxial creep tests. In the tests, two specimens were subjected to isotropic drained creep for either 3 min or 180 min at p0 = 400 kPa. Then, they were subjected to drained constant rate loading (CRL, e_ a ¼ 0:125%= min) triaxial compression until q = 200 kPa followed by drained creep for 24 h. Fig. 14a illustrates the stress-strain relationships from the experiment and the numerical simulation, and Fig. 14b depicts the increase in the creep strain with time. It can be seen from both the tests and simulations that a longer isotropic creep time results in a smaller shear strain. The parameters in the

4. Conclusions In this study, the time-dependent behaviors of coarse-grained soils have been investigated. A double-yield-surface EVP model is proposed for simulating such behaviors. The primary features of the proposed model are the following: 1. A strain-hardening model using plastic strains as the hardening parameters for both yield surfaces has been developed. 2. The stress dependence of the material modulus, strength and dilatancy is considered. 3. The model is capable of modelling the loading rate effects on shearing and the shear creep phenomenon, including transient creep and stationary creep. 4. The simulation of compression creep has considered the particle crush-induced stress level dependence of the creep coefficient, which is a newly found mechanical characteristic of coarse-grained soils. 5. Compared with the original hardening soil model, 6 more parameters are added for the consideration of stress dependence, shear creep and compression creep. The soil behavior at different stress states can be described with one set of parameters. The differential equations that describe creep behaviors are derived from basic phenomena discovered through experiments. By comparing with certain testing results, it has been demonstrated that this model is capable of predicting the timedependent behaviors of coarse-grained soils with fair accuracy. In the future, further research is desired to explore the micro mechanism of the time-dependent behavior of coarse-grained soils in complex situations, which would be helpful for improving the proposed model.

simulation include u = 34.95°, Du = 0.72°, Rf = 0.9, Eref 50 = 12.4 MPa, Eref oed = 15.5 MPa,

Eref ur = 37.2 MPa,

m = 0.5, tur = 0.2, uv 0 = 30°, 4 e_  = 1  10 /min, f = 0.01729, C ref e;oed = 8.5  10 , mc = 0.5, and

Acknowledgements

pref = 0.1 MPa. Points A and B are two points on the simulation curve that correspond to the same time from the start of isotropic compression

This study was financially supported by the National Key Fundamental Research and Development Program of China (Project No. 2014CB047003).

6

70

Y. Kong et al. / Computers and Geotechnics 85 (2017) 59–70

References [1] Alonso EE, Oldecop L, Pinyol NM. In: Kolymbas D, Viggiani G, editors. Long term behaviour and size effects of coarse granular media. Berlin: SpringerVerlag Berlin; 2009. p. 255–81. [2] Anh Dan LQ, Tatsuoka F, Koseki J. Viscous effects on the stress-strain behavior of gravelly soil in drained triaxial compression. Geotech Test J 2006;29 (4):330–40. [3] Augustesen A, Liingaard M, Lade PV. Evaluation of time-dependent behavior of soils. Int J Geomech 2004;4(3):137–56. [4] Bjerrum L. Engineering geology of Norwegian normally-consolidated marine clays as related to settlements of buildings. Geotechnique 1967;17(2):83–118. [5] Cao GX. Study on post-construction settlement of high fill foundation in mountainous airport. Beijing: Tsinghua University; 2012. [6] Colliat-Dangus JL, Desrues J, Foray P. Triaxial testing of granular soil under elevated cell pressure. Advanced triaxial testing of soil and rock 1988;SATM STP 977:290–310. [7] Den Haan EJ. Vertical compression of soils. Delft University of Technology; 1994. [8] Duncan JM, Chang C. Nonlinear analysis of stress and strain in soils. J Soil Mech Found Divis 1970;96(5):1629–53. [9] Duncan JM, Wong KS, Mabry P. Strength, stress-strain and bulk modulus parameters for finite element analyses of stresses and movements in soil masses. University of California; 1980. [10] Eliahu U, Harrell E. Long-term performance of engineered fills Presented. In: Geo-congress 2013@ stability and performance of slopes and embankments III. p. 628–41. [11] Hardin BO. Crushing of soil particles. J Geotech Eng 1985;10(111):1177–92. [12] Kiyota T, Tatsuoka F, Yamamuro J. Drained and undrained creep characteristics of loose saturated sand. In: Proceedings of the sessions of the geo-frontiers 2005 congress, site characterization and modeling (GSP 138). Reston: ASCE; 2005. [13] Lade PV, Liu CT. Modeling creep behavior of granular materials. In: Proceedings of the 10th international conference on computer methods and advances in geomechanics. Tucson, Arizona, USA: CRC Press; 2001. p. 277–83. [14] Lade PV, Yamamuro JA, Bopp PA. Significance of particle crushing in granular materials. J Geotech Eng 1996;122(4):309–16. [15] Liu H. A systematic research on the deformation and stability of high embankment of Jiuzhai-Huanglong Airport, Sichuan China. Earth Environ 2006;33(B10):133–5. [16] Lydzba D, Pietruszczak S, Shao J. Intergranular pressure solution in chalk: a multiscale approach. Comput Geotech 2007;34(4):291–305. [17] Mcdowell GR, Bolton MD. On the micromechanics of crushable aggregates. Geotechnique 1998;48(5):667–79. [18] Mcdowell GR, Khan JJ. Creep of granular materials. Granul Matter 2003;5:115–20. [19] Mesri G, Barames V. Compression of granular materials. Can Geotech J 2009;46 (4):369–92. [20] Morsy MM. Effective stress modeling of creep behaviour of clay. Edmonton, Alberta: University of Alberta; 1994. [21] Oldecop LA, Alonso EE. Theoretical investigation of the time-dependent behaviour of rockfill. Geotechnique 2007;57(3):289–301. [22] Parkin AK. Settlement rate behaviour of some fill dams in Australia. In: Proceedings of 11th ICSMFE. San Francisco. p. 2007–10. [23] Parkin AK. Creep of rockfill. Springer; 1991. p. 221–37. [24] Perzyna P. Fundamental problems in viscoplasticity. Adv Appl Mech 1966;9:243–377. [25] Perzyna P. Thermodynamic theory of viscoplasticity. Adv Appl Mech 1971;11:313–54. [26] Pietruszczak S, Lydzba D, Shao J. Description of creep in inherently anisotropic frictional materials. J Eng Mech 2004;130(6):681–90.

[27] Pietruszczak S, Lydzba D, Shao J. Modelling of deformation response and chemo-mechanical coupling in chalk. Int J Numer Anal Met 2006;30 (10):997–1018. [28] Rowe PW. The stress-dilatancy relation for static equilibrium of an assembly of particles in contact. Roy Soc 1962:500–27. [29] Schanz T, Vermeer PA. Angles of friction and dilatancy of sand. Geotechnique 1996;46(1):145–52. [30] Schanz T, Vermeer PA, Bonnier PG. The hardening soil model: formulation and verification. In: Presented in beyond 2000 in computational geotechnics–10 years of PLAXIS international. Rotterdam: AA Balkema, Publishers; 1999. p. 281–96. [31] Sun DA, Huang WX, Sheng DC, et al. An elastoplastic model for granular materials exhibiting particle crushing. Key Eng Mater 2007;340:1273–8. [32] Tapias M, Gili J, Alonso EE. Scale effects in rockfill behaviour. Géotechn Lett 2012;2(July-September):155–60. [33] Tatsuoka F. Inelastic deformation characteristics of geomaterial. In: Presented in soil stress-strain behavior: measurement, modeling and analysis. Roma: Springer; 2006. p. 1–108. [34] Tatsuoka F, Benedetto HD, Enomoto T, et al. Various viscosity types of geomaterials in shear and their mathematical expression. Soils Found 2008;1 (48):41–60. [35] Taylor DW. Fundamentals of soil mechanics. New York: John Wiley & Sons; 1970. [36] Vermeer PA, Neher HP. A soft soil model that accounts for creep. In: Presented in beyond 2000 in computational geotechnics–10 years of PLAXIS international. Rotterdam: AA Balkema, Publishers; 1999. p. 1–13. [37] Wahls HE. Analysis of primary and secondary consolidation. J Soil Mech Found 1965;91 [Proc. Paper 4336]. [38] Wang HJ, Yin ZZ. Creep tests of rockfill and double-yield surface creep model. Chin J Geotech Eng 2008;30(7):959–63. [39] Wang HJ, Yin ZZ. Experimental study on creep deformation of rockfill considering load action. Hydro-Sci Eng 2008;02:48–53. [40] Wang N, Yao Y. A generalized constitutive model considering sand crushing. Jpn Geotech Soc Spec Publ 2015;1(2):12–5. [41] Xu M, Song EX. Numerical simulation of the shear behavior of rockfills. Comput Geotech 2009;36(8):1259–64. [42] Xu M, Song EX, Cao GX. Compressibility of broken rock-fine grain soil mixture. Geomech Eng 2009;1(2):169–78. [43] Yang C, Carter JP, Sheng DC, et al. An isotach elastoplastic constitutive model for natural soft clays. Comput Geotech 2016;77:134–55. [44] Yao YP, Huang G. Visco-elasto-plastic constitutive model for rockfill considering particle crushing. Ind Constr 2010;03:71–6. [45] Yao YP, Hou W, Zhou AN. UH model: three-dimensional unified hardening model for overconsolidated clays. Geotechnique 2009;59(5):451. [46] Yao YP, Sun DA, Matsuoka H. A unified constitutive model for both clay and sand with hardening parameter independent on stress path. Comput Geotech 2008;35(2):210–22. [47] Yao YP, Wang ND. Transformed stress method for generalizing soil constitutive models. J Eng Mech 2013;140(3):614–29. [48] Yao YP, Yamamoto H, Wang ND. Constitutive model considering sand crushing. Soils Found 2008;48(4):603–8. [49] Yin JH. Constitutive modelling of time-dependent stress-strain behaviour of soils. Winnipeg: The University of Manitoba; 1990. [50] Yin JH, Graham J. Elastic viscoplastic modelling of the time-dependent stressstrain behaviour of soils. Can Geotech J 1999;36(4):736–45. [51] Zhou MJ, Song EX. A random virtual crack DEM model for creep behavior of rockfill based on the subcritical crack propagation theory. Acta Geotech 2016;11(4):827–47. [52] Zienkiewicz OC, Cormeau IC. Visco-plasticity-plasticity and creep in elastic solids—a unified numerical solution approach. Int J Numer Meth Eng 1974;8 (4):821–45.