Elasto-plastic behaviour of frozen soil subjected to long-term low-level repeated loading, Part II: Constitutive modelling

Cold Regions Science and Technology 122 (2016) 58–70

Contents lists available at ScienceDirect

Cold Regions Science and Technology journal homepage: www.elsevier.com/locate/coldregions

Elasto-plastic behaviour of frozen soil subjected to long-term low-level repeated loading, Part II: Constitutive modelling Qionglin Li a,b,c, Xianzhang Ling a,c,⁎, Daichao Sheng b,d a

School of Civil Engineering, Harbin Institute of Technology, Heilongjiang, Harbin 150090, China Centre of Excellence for Geotechnical Science and Engineering, The University of Newcastle, Callaghan, NSW 2308, Australia State Key Laboratory of Frozen Soil Engineering, Chinese Academy of Science, Gansu, Lanzhou 730000, China d National Engineering Laboratory for High-Speed Railway Construction, Central South University, Changsha, Hunan 410075, China b c

a r t i c l e

i n f o

Article history: Received 16 July 2015 Received in revised form 4 November 2015 Accepted 18 November 2015 Available online 27 November 2015 Keywords: Frozen soil Empirical equations Long-term low-level repeated stress Constitutive model

a b s t r a c t Based on the conclusions from the experimental investigation in Part I, empirical equations for the accumulated shear strain, accumulated direction ratio and elastic modulus of frozen soil are derived. The basic issues observed in the test results can also be observed in these empirical equations. Subsequently, an elaborate constitutive model for frozen soil subjected to long-term low-level repeated loading is produced by combining the empirical equations and classical elasto-plastic theory. This model accounts for the dependency of the accumulated behaviour on the initial stress state, repeated stress amplitude and frozen soil strength. In addition, the evolution behaviour of the elastic modulus with accumulated strain is also considered. The model is verified for frozen soil with the aid of the triaxial test results represented in Part I. This study is the first attempt to model the elasto-plastic behaviour under long-term low-level repeated loading for frozen soil. © 2015 Elsevier B.V. All rights reserved.

1. Introduction In Part I, the elasto-plastic behaviour, including the accumulated amount, direction and elastic modulus, of frozen soil subjected to long-term low-level repeated loading was presented based on experimental investigations. For the further development of deformation analyses and design methods for frozen soil foundations, an efficient numerical approach that can capture the long-term accumulated and elastic behaviour of frozen soil under a large number of loading cycles is indispensable. Thus, the purpose of the present study is to develop a constitutive model by combining the framework of time-dependent elasto-plastic theory and long-term mechanical behaviours of frozen soil. Although many elasto-plastic models or creep models on frozen soils under static loading have been proposed (Lai et al., 2009, 2014; Li et al., 2011; Wang et al., 2014; Yang et al., 2010), comparatively little research is done in the field of cyclic elasto-plastic models of frozen soils. Accordingly, in this subject not much experiences is available on the aspect regarding constitutive modelling of frozen soil under long-term repeated loading. Moreover, the studies on modelling the long-term elasto-plastic behaviour of the unfrozen sands or clays are also not investigated sufficiently. Previously, a number of elasto-plastic models, including multisurface types and bounding surface types (Dafalias, 1986; Mroz et al., ⁎ Corresponding author at: School of Civil Engineering, Harbin Institute of Technology, Heilongjiang, Harbin, 150090, China. Tel.: +86 189 461 52163. E-mail address: [email protected] (X. Ling).

http://dx.doi.org/10.1016/j.coldregions.2015.11.009 0165-232X/© 2015 Elsevier B.V. All rights reserved.

1979; Prevost, 1985; Yang et al., 2003), have been developed for sand or clay with kinematic hardening. These models can predict the details of stress–strain hysteresis loops, which require hundreds of load increments per cycle. Therefore, they have generally been used in the prediction of soil behaviour under earthquake vibrations. Compared with earthquake vibrations, these long-term soil dynamic problems represent a large number of loading cycles. These models are of limited use in problems induced by long-term vibration because of the high computing cost and uncontrollable cumulative calculation errors. Another efficient and economic countermeasure against long computing times and error accumulation is the utilization of empirical approaches. The power equation first proposed by Monismith (1975) is frequently used to predict the permanent deformation of soil under repeated or cyclic loading. Li and Selig (1996) introduced the ratio of the cyclic deviator stress over the static deviator stress at failure into the constants of the power equation. In this method, the level of cyclic loading and physical soil state are considered. Chai and Miura (2002) further considered the effect of the initial static deviator stress and shear strength to develop a new power equation for calculating the cumulative axial strain of soft cohesive soil under repeated loading. Various types of logarithmic formulations are used to describe the evolution of strain accumulation in relation to the number of loading cycles. The parameters in most of these empirical models depend on the stress state (static and dynamic stress state), and soil properties were determined by additional experimental results (Behzadi and Yandell, 1996; Karg et al., 2010; Li et al., 2013; Niemunis et al., 2005; Sawicki and Swidzinski, 1989; Sweere, 1990). Certain researchers

Q. Li et al. / Cold Regions Science and Technology 122 (2016) 58–70

recommend the use of accumulated strain at a certain cycle (εacc(N = 1), εacc(N = 100) or εacc(N = 1000)) as a reference to predict the long-term accumulated strain, with the effects of stress state or soil properties considered directly or indirectly (Guo et al., 2013). Approaches formulated in other structures have also been proposed (Drabkin et al., 1996; François et al., 2010; Lekarp and Dawson, 1998; Paute et al., 1996). These empirical equations, however, are not frequently used, and their details are not represented in this study. Most of the above-mentioned empirical models were proposed for accumulated axial strain, and cyclic or repeated stress was assumed to be constant during the long-term loading process. However, for the boundary-value problem, the cyclic or repeated stress distribution must be determined by a general constitutive model (elastic or elastoplastic model) before the settlement or residue deformation prediction. This method is referred to as the “mechanistic–empirical method” and has been frequently used in previous studies (Puppala et al., 1999). However, stress redistribution induced by the accumulated strain is inevitable during long-term loading processes, and the assumption that the cyclic or repeated stress is constant is irresponsible. Therefore, an advanced model suitable for implementation in a finite element framework must be developed, and the deviatoric and volumetric portions of the accumulated strain must be considered. Marr and Chrisitian (1981), Bouckovalas et al. (1984) and Kaggwa et al. (1991) first proposed separate empirical equations for the accumulated deviatoric and volumetric strain in relation to the number of loading cycles, which are then coupled by the accumulated direction observed in recent studies. Wichtmann et al. (2006, 2014) found that the accumulated direction of sand was consistent with the flow rule of the modified Cam Clay model, which was integrated by the average stress ratio. The mechanism of accumulation for granular soil was decomposed into frictional sliding and volumetric compaction by François et al. (2010) and Karg et al. (2010), and these phenomenological laws were advanced for the volumetric portion and deviatoric portion, respectively. These authors suggested that a portion of the accumulated volumetric strain resulted from dilation induced by the deviatoric deformation, and the accumulated direction of frictional sliding is derived from the Drucker–Prager or Mohr–Coulomb criterion. Furthermore, by combining the coupled formulations for accumulated strain with a classical plasticity framework, Niemunis et al. (2005), François et al. (2010), and Karg et al. (2010) proposed accumulation models. In previous studies, a large number of models for the resilient modulus derived from the axial stress–strain relationship were proposed to evaluate the repetitive behaviour of subgrade soils. Certain significant influencing factors, such as the stress state (static and repeated) and soil physical properties, were integrated into these models (Drumm et al., 1997; Fall et al., 2008; Guo et al., 2013; Hicks and Monismith, 1971). Poisson's ratio is another elastic property index for material, and it was assumed to be constant in the above models. However, studies (Brown and Hyde, 1975; Hicks and Monismith, 1971; Kolisoja, 1997; Sweere, 1990) have shown that Poisson's ratio is not a constant and varies with applied stress. Thus, separate mathematical formulations for defining the Poisson's ratio were proposed by certain researchers (Boyce, 1980; Hicks and Monismith, 1971). The resilient modulus and Poisson's ratio can also be replaced by the bulk and shear moduli. Boyce (1976) and Sweere (1990) developed generalized formulas for the bulk and shear modulus to define their dependence on the stress state. Nevertheless, the evolution of the elastic modulus represented in certain experimental results during cyclic or repeated loading is not considered in the above models. In the accumulation models proposed by Niemunis et al. (2005), François et al. (2010), and Karg et al. (2010), elastic behaviour is defined by the pressure-dependent stiffness and a constant Poisson's ratio, but the evolution behaviour of the elastic modulus during long-term loading processes, which significantly contributes to the stress redistribution, was not considered. Moreover, most of the current research focuses on unfrozen soil, and a systematic attempt to model the long-term dynamic behaviour of

59

frozen soil has not been performed. Distinct differences have been found in the long-term elasto-plastic behaviour between frozen and unfrozen soil, which are presented in Part I. Therefore, it is necessary to acquire additional information on long-term elasto-plastic behaviour and develop an efficient and economic constitutive model for frozen soil. This paper is organized as follows. First, empirical equations are proposed for accumulated shear strain, accumulated direction ratio and elastic modulus based on the experimental data in Part I. Based on these empirical equations, an efficient strategy is presented that incorporates a constitutive model for simulating the elasto-plastic behaviour of frozen soil under long-term low-level repeated loading. Subsequently, the numerical integration of this constitutive model is conducted by employing a fully explicit Euler-forward algorithm, and a reasonable integration step is determined. In the final section, the proposed constitutive model and model parameters are verified against the test results presented in Part I. In addition, the simulation capacity and validity and shortcomings of this constitutive model are discussed. 2. Empirical equations 2.1. Empirical equation for accumulated shear strain In Part I, the accumulated shear strain increased with the number of repeated loading cycles and was significantly influenced by the initial stress ratio, initial mean stress, repeated stress amplitude and strength of the frozen soil were detected As a result, the above-mentioned influencing factors should be directly or indirectly integrated into the empirical equation for accumulated shear strain. The test results indicate that at a given number of loading cycles, additional plastic shear strain will accumulate for samples with a higher initial stress ratio. Thereafter, the function g(η0) (shown in Eq. (1)) is proposed in a simple method used to consider the influence of initial stress ratio on the magnitude of accumulated shear strain.

 g η0 ¼ a1  exp b1  η0

ð1Þ

where a1 and b1 are material parameters, with a1 = 0.74 and b1 = 1.22; and η0 = 0.25 is used as a reference initial stress ratio η0,ref, for which g(η0,ref) = 1 holds. Eq. (1) is fitted by the normalized data from the test results of Series 1 presented in Part I and shown in Fig. 1. The magnitude of the accumulated shear strain is significantly dependent on the initial mean stress, and additional plastic shear strain is accumulated for specimens with a higher initial mean stress, which was demonstrated in Part I. Furthermore, Eq. (2) is formulated to

Fig. 1. Normalized accumulated shear strain versus initial stress ratio.

60

Q. Li et al. / Cold Regions Science and Technology 122 (2016) 58–70

consider the dependence on the initial mean stress as follows:   p g ðp0 Þ ¼ a2  exp b2 0 : pre f

ð2Þ

The material parameters a2 = 0.66 and b2 = 0.83 are generally independent of the number of loading cycles. A default value of 1 MPa for the reference pressure pref is used in this paper to make the pressure ratio dimensionless. This value for the reference pressure will also be used in certain empirical formulations. The reference initial mean stress p0,ref is set to 0.6 MPa, for which g(p0,ref) = 1 is true. The function g(p0) is derived from the fitting curve of the normalized data of the test in Series 2 as shown in Fig. 2. The influence of repeated stress amplitudes on the magnitude of the accumulated shear strain is significant, and additional plastic shear strain will accumulate under a higher repeated stress. The following function g(σd) is formulated to describe this dependence:     σ g ðσ d Þ ¼ a3  exp b3  d −1 pre f

ð3Þ

where a3 and b3 are two material parameters. The normalized accumulated shear strains from test Series 3 are plotted versus the repeated stress amplitude as shown in Fig. 3. The fitting curve from Eq. (3) with a3 = 0.19 and b3 = 1.22 in this figure demonstrates that the loss of accuracy caused by the assumption a3,b3 = const is acceptable. When the reference repeated stress amplitude σd,ref is set to 1.5 MPa, then the function g(σd,ref) = 1 is true. Any other influence factors of frozen soil, including the temperature and initial moisture content, have a significant effect on the magnitude of the accumulated shear strain. To account for the influence of them, the strength s0.3 obtained from the triaxialtriaxial compression test at a confining pressure of 0.3 MPa is introduced into this study to represent variations in the physical behaviour. From the test results in Part I, specimens with a higher strength will accumulate less plastic shear strain up to a given number of loading cycles. The dependence on strength s0.3 is approached well by the function g(s0.3) as follows: g ðq0:3 Þ ¼

a4 −b4 : s0:3 =pre f

ð4Þ

Fig. 3. Normalized accumulated shear strain versus repeated stress amplitude.

In brief, the reference states are η0,ref = 0.25, p0,ref = 0.6 MPa, σd,ref = 1.5 MPa and q0.3,ref = 4.0 MPa. The accumulated shear strains from test Series 1 to 5 are normalized by the functions g(η0), g(p0), g(σd) and g(s0.3) to the values at the reference state as shown in Fig. 5. All of the normalized test data fall into a narrow band, which can be approximated by the reference accumulation curve as follows: εacc q;re f ðN Þ ¼ C N1 ln ð1 þ C N2 N Þ

ð5Þ

where CN1 and C N2 are material parameters. For the fitting curve, CN1 = 0.15 and CN2 = 0.26. For the bottom and top curves, 0.26 is used for CN2, whereas a range from 0.08 to 0.20 is used for C N1. In this reference accumulation curve, the behaviour in which the accumulated shear strain increases but the accumulated rate decreases with an increasing number of loading cycles can be well described. The influence functions g(η0), g(p0), g(σd), g(s0.3) and reference accumulation curve εacc q,ref(N) have been presented, and they can be integrated into an explicit empirical equation for the accumulated shear strain as follows (Eq. (6)):

The material parameters a4 and b4 in Eq. (4) are independent of the number of loading cycles. The fitting curve with a4 = 7.67 and b4 = 0.92 from this equation, which is shown in Fig. 4, indicates that this function can be confirmed by the normalized accumulated shear strain of test Series 4 and 5. The reference strength corresponding to g(s0.3,ref) = 1 is chosen as s0.3,ref = 4.0 MPa.

Eq. (6) can be used to predict the accumulated shear strain of frozen soil with different initial stress states, repeated stress amplitudes and

Fig. 2. Normalized accumulated shear strain versus initial mean stress.

Fig. 4. Normalized accumulated shear strain versus common strength.


 acc εacc q ðN Þ ¼ ε q;re f ðN Þ  g ðσ d Þ  g η0  g ðp0 Þ  g ðs0:3 Þ:

ð6Þ

Q. Li et al. / Cold Regions Science and Technology 122 (2016) 58–70

61

Fig. 5. Normalized accumulated shear strain versus the number of repeated loading cycles. Fig. 6. Normalized fitting parameter αd form Eq. (7) versus initial stress ratio.

strengths. However, the test program used in this investigation does not provide a sufficient range of influencing factors.

According to the description above, α d can be approximated by Eq. (11):

2.2. Empirical equation for accumulated direction ratio


 α d ¼ α re f d η0 dðp0 Þdðσ d Þ:

From the test results in Part I, the accumulated direction for frozen soil under long-term low-level repeated loading is significantly affected by repeated loading processes, initial stress ratios, initial mean stress, repeated stress amplitudes and frozen soil strength. A reasonable empirical equation for the accumulated direction ratio should integrate such effects. The accumulated direction ratio proposed in Part I is rewritten in Eq. (7) to incorporate the evolution during repeated loading process. dg ¼

dεacc v;N dεacc q;N

 ¼ dg0 exp α d εacc q;N

ð7Þ

According to the analysis in Part 1, the initial accumulated direction ratio dg0 was found to be relatively insensitive to changes in the initial stress ratio, initial mean stress and repeated stress amplitude. The effect on the accumulated direction from the previously mentioned influencing factors can be characterized by the effect on αd. With increases of the initial stress ratio, initial mean stress and repeated stress amplitude, αd clearly decreases. To account for these effects, three functions are proposed as follows:

 d η0 ¼ c1 exp d1 η0 :

ð8Þ

  p dðp0 Þ ¼ c2 exp d2 0 : pre f

ð9Þ

  σ dðσ d Þ ¼ c3 exp d3 d : pre f

ð10Þ

In Eqs. (8)–(10), the material parameters c1 = 1.48, d1 = − 1.58, c2 = 2.92, d2 = − 1.82, c3 = 13.07 and d3 = − 1.71 are independent of the number of loading cycles, and they are fitted well by the normalized test results in Figs. 6–8. In accordance with the explicit empirical equation for the accumulated shear strain, η0,ref = 0.25, p0,ref = 0.6 MPa and σd,ref = 1.5 MPa are also used as the reference states, wherein d(η 0,ref) = 1, d(p0,ref) = 1 and d(σd,ref) = 1. The normalized αd from test Series 1, 2, 3, 6 and 7 presented in Part I by functions d(η 0 ), d(p0 ), and d(σd ) fall into a relatively small band. The values of α ref used in this paper are in the range of 3.5 to 5.3.

ð11Þ

In addition, the initial accumulated direction ratio dg0 is detected to be significantly affected by the strength s0.3 in this section. With increases of strength s0.3, the initial accumulated direction ratio increases nonlinearly as shown in Fig. 9. The initial accumulated direction ratio normalized by the pre-set reference strength s0.3,ref = 4.0MPa can be described by Eq. (12):   s0:3 dðs0:3 Þ ¼ c4 exp d4 pre f

ð12Þ

where the material parameters are c4 = 0.07 and d4 = 0.66. When the strength s0.3 is set as the reference strength s0.3,ref = 4.0 MPa, then this function at d(s0.3,ref) = 1.0 is true. According to Eq. (13), the initial

Fig. 7. Normalized fitting parameter αd form Eq. (7) versus initial mean stress.

62

Q. Li et al. / Cold Regions Science and Technology 122 (2016) 58–70

repeated loading. In addition, the test results from Part I also indicate that the elastic modulus is independent of the initial stress ratio and repeated stress amplitudes but significantly dependent on the initial mean stress. The evolution behaviour of the elastic shear modulus and bulk modulus can be approximated by Eqs. (15) and (16), which were proposed in Part I.

Fig. 8. Normalized fitting parameter αd form Eq. (7) versus repeated stress amplitude.

accumulated direction ratio dg0 can be approximated by Eq. (13): dg0

   s0:3 ¼ dg0;re f c4 exp d4 pre f

ð13Þ

in which dg0,ref is the initial accumulated direction ratio at the reference strength s0.3 = 4 MPa. dg0,ref can be assumed to be equal to the mean value of the normalized test results at 0.16. Thus far, empirical equation configured for the accumulated direction ratio can be described in Eq. (14): h i
 dg ¼ dg0;re f dðs0:3 Þ exp α re f d η0 dðp0 Þdðσ d Þεacc q;N

 GN ¼ G0 1 þ α 1 εacc v;N

ð15Þ

 K N ¼ K 0 1 þ α 2 εacc v;N

ð16Þ

According to the related analysis in Part I, the dependence of the initial mean stress on the elastic modulus can be characterized by the effect on the initial elastic modulus G0 or K0. For further study, the initial elastic moduli G0 and K0 are plotted against the initial mean stress in Fig. 10(a) and (b), respectively. Subsequently, Eqs. (17) and (18) are formulated to approach the pressure-dependent behaviour:  n1 p G0 ¼ k1 p0;re f 1 þ 0 p0:re f K 0 ¼ k2 p0;re f 1 þ

p0

!n2

p0;re f

ð18Þ

where n 1 and n 2 are non-dimensional material constants, with n1 = 0.4032 and n2 = 0.2866; and k1 and k2 are non-dimensional parameters dependent on the physical behaviour of the frozen soil.

ð14Þ

Eq. (14) is applicable for predicting the evolving behaviour with accumulated strain and effects of the initial stress state, repeated stress amplitude and frozen soil strength. However, the limited test program used in this investigation does not provide sufficient ranges of the influencing factors discussed above. 2.3. Empirical equation for elastic modulus The test results from Part I indicate that the elastic modulus, which is defined by the resilient modulus, increases linearly with accumulating compression plastic volumetric strain under long-term low-level

Fig. 9. Normalized accumulated direction ratio versus strength of frozen soil.

ð17Þ

Fig. 10. The initial modulus (shear and bulk) versus initial mean stress.

Q. Li et al. / Cold Regions Science and Technology 122 (2016) 58–70

63

stress ratio, initial mean stress and repeated stress amplitude. To obtain reasonable values for these parameters, the normalized curves acc GN/G0 ~ εacc v,N and KN/K0 ~ εv,N from test Series 1, 2, 3, 6 and 7 are plotted in Fig. 11(a) and (b), respectively. Fig. 11(a) and (b) demonstrate that the normalized curves acc GN/G0 ~ εacc v,N and KN/K0 ~ εv,N fall into their bands respectively. However, to keep the number of material constants manageable, using values from the fitting curves, which are α1 = 0.05 and α2 = 0.07, is proposed. In addition, the loss of accuracy caused by the assumption that α1 and α2 are constant values is acceptable. 3. Constitutive model for long-term deformation 3.1. Stress and strain variables Tensor quantities are denoted by bold-faced letters in direct notation, and a colon denotes an inner product. Stress, strain and their increments are positive in compression. The mean stress and volumetric strain are calculated as follows: p¼

1 trσ 3

εv ¼ trε

ð19Þ ð20Þ

where σ is the stress tensor and ε is the strain tensor. The deviatoric stress and strain tensors are defined as follows: s ¼ σ−pI e ¼ ε−

εv I 3

ð21Þ ð22Þ

where I = δij denotes the identity tensor of rank two, or the Kronecker delta δij(1 for i = j, 0 for i ≠ j). Moreover, the following invariants will frequently be used in this model definition: q¼ Fig. 11. The normalized elastic modulus (shear and bulk) versus accumulated volumetric strain.

For the frozen soil specimens at freezing temperature T = − 5 °C and initial moisture content W = 14.5% used in Part I, k1 = 681.66 and k2 = 474.86. The default value of 1.0 MPa is used as the reference initial pressure p0,ref in these formulations. The test results presented in Part I also indicate that the material constants α1 and α2 in Eqs. (15) and (16) are independent of the initial

rffiffiffiffiffiffiffiffiffiffiffiffi 3 s:s 2

εq ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffi 2 e:e 3

ð23Þ

ð24Þ

in which q denotes the deviatoric stress invariant (shear stress) and εq is the deviatoric strain invariant (shear strain). 3.2. Basic simplification for simulation Based on previous statements and experimental findings, simulating the long-term response is challenging because of the large number of

Fig. 12. The equivalent process of long-term low-level repeated loading.

64

Q. Li et al. / Cold Regions Science and Technology 122 (2016) 58–70

loading cycles and complex stress states for different frozen soil elements at different loading cycles. To address the above difficulties, three reasonable simplifications are proposed. (a): As previously discussed, accumulated and elastic deformation, rather than details of the hysteresis loop are concerned for simulations of the long-term response. Therefore, long-term low-level repeated loadings can be reduced to a constant static loading equivalent to the amplitude of repeated loading. Moreover, the loading time t of the repeated loadings is replaced by a continuous number of repeated loading cycles N. This equivalent loading process is called “pseudo-repeated loading versus N.” (b): In the simulation of boundary values problem for a frozen soil foundation system, the static stress induced by self-weight or static loading must be completed before applying the longterm low-level repeated loading. For simplication, in the process of long-term low-level repeated loading, the static stress redistribution induced by accumulated strain is ignored. Then, the initial stress state of each element is assumed to be constant during the long-term low-level repeated loading process. (c): Compared with the stress state of samples under a triaxialtriaxial test, the stress state of any frozen soil element in a foundation system is more complex. In the present study, the repeated stress amplitude for each element can be calculated by Eq. (25): rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 σ d ¼ qd ¼ sN : sN 2

ð25Þ

in which sN is the deviatoric stress tensor induced by pseudo-repeated loading at the Nth cycle of repeated loading. Based on the above simplifications and discussions, a numerical scheme is developed herein to analyse a frozen soil geotechnical system experiencing long-term low-level repeated loading. In this scheme, the two loading steps in the traditional simulation of bounding values for a frozen soil foundation system can be transformed to three steps, which is shown in Fig. 12: Step 1, initial static loading process versus loading time t, in which the static stress distribution of the computational field can be obtained; Step 2, pseudo-repeated loading process versus loading time t, in which the initial repeated stress amplitude of each element can be obtained; Step 3, pseudo-repeated loading process versus continuous number of repeated loading cycles N, in which the accumulated strain, elastic strain and dynamic stress redistribution can be obtained. Thus, the first two steps provide the initial values of the constitutive model used in step 3. 3.3. General formulation

1 acc 3st acc 1 acc acc acc ε_ N ¼ e_ N þ ε_ v;N I ¼ þ ε_ I ε_ 3 2qt q;N 3 v;N

ð28Þ

in which st and qt are the total deviatoric stress tensor and total shear stress (total deviatoric stress invariant), respectively, which include the initial static and repeated elements induced by pseudo-repeated loading, which can be calculated by Eq. (29): rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 st : st : ð29Þ 2 Upon substitution of the definition for the accumulated direction ratio dg as shown in Eq. (7), an alternative formulation for the accumulated strain rate is formulated:

st ¼ s0 þ sN ; q t ¼



 3st 1 acc þ dg I ε_ q;N 2qt 3

ð30Þ

acc

where ε_ q;N is the accumulated shear strain rate versus loading cycles N, which can be obtained from the empirical equation for the accumulated shear strain as shown Eq. (31). acc

acc ε_ q;N ¼

dεacc dεq;N
 q;re f ¼ g ðσ d Þ  g η0  g ðp0 Þ  g ðs0:3 Þ  dN dN

ð31Þ

ð26Þ

e acc where ε_ N ; ε_ N ; ε_ N denote the strain tensor rate versus loading cycles N and its elastic (resilient) and accumulated parts, respectively. The elastic strain can be determined by a generalized Hooke's law as shown in Eq. (27):

1 1 1 e e s_ N þ ε_ N ¼ e_ N þ ε_ ev I ¼ p_ I 3 2GN 3K N N

in which GN and KN represent the shear and bulk modulus, respectively, at loading cycle N. As discussed in Section 2, the dependence of these parameters on the initial mean stress and accumulated volumetric strain after N loading cycles can be described by Eqs. (15) to (18). The accumulated strain rate versus loading cycles N can analogously be decomposed into a deviatoric volumetric elements. In the present study, the deviatoric accumulated strain tensor rate is considered to be dependent on the total deviatoric stress tensor, and it can be written as follows:

acc ε_ N ¼

A robust constitutive model for frozen soil subjected to long-term low-level repeated loading can be framed within a finite-element formulation that can reflect the salient behaviour trends that were previously described. As described in Part I, the strain induced by long-term low-level repeated loading can be decomposed into the unrecovered accumulated strain and resilient strain. Therefore, the strain rate versus loading cycles N can be compiled with the additive decomposition: e acc ε_ N ¼ ε_ N þ ε_ N

Fig. 13. The accumulated shear strain with different integration step size versus the number of repeated loading cycles.

ð27Þ

Table 1 Summary of model parameters. Accumulation parameters a1 b1 a2 b2 a3

0.74 1.22 0.66 0.83 0.19

b3 a4 b4 CN1 CN2

1.22 7.67 0.92 0.08–0.20 0.26

Direction parameters c1 d1 c2 d2 c3

1.48 −1.58 2.92 −1.82 13.07

d3 αref c4 d4 dg0,ref

Elastic parameters −1.71 3.5–5.3 0.07 0.66 0.16

n1 n2 α1 α2

0.40 0.29 0.05 0.07

Q. Li et al. / Cold Regions Science and Technology 122 (2016) 58–70

Fig. 14. Comparison for accumulated shear strain between the test and model.

65

66

Q. Li et al. / Cold Regions Science and Technology 122 (2016) 58–70

Fig. 14 (continued).

From Eq. (5), the reference accumulated shear strain rate versus loading cycles N can be described as follows: dεacc q;re f dN

¼ C N1 C N2 exp −

εacc q;re f C N1

Upon the substitution of Eqs. (31), (32) and (33) into Eq. (30), the accumulated strain rate tensor can be obtained as follows:

! :

ð32Þ

acc ε_ N ¼

dεacc N ¼ dN

!  εacc 3st 1 q;re f þ dg I κ 1 exp − g ðσ d Þ 2qt 3 C N1



ð35Þ

in which

According to assumption (b) in Section 3.2, functions g(η0), g(p0)0 and g(s0.3) have been determined in the first two steps and remain constant during the pseudo-repeated loading process versus N. Similarly, the functions d(η0), d(p0) and d(s0.3) can be determined in the first two steps. For simplification, alternative variations κ1 and κ2 are defined as follows:

The dependence of the initial accumulated direction ratio dg0 on the strength can be described by Eq. (13).


 κ 1 ¼ C N1  C N2  g η0  g ðP 0 Þ  g ðs0:3 Þ:

ð33Þ

3.4. Numerical algorithm


 κ 2 ¼ α re f d η0 dðp0 Þ:

ð34Þ

A numerical scheme framed in the finite element method is developed herein to analyse frozen soil systems that experience long-term

h i dg ¼ dg0 exp κ 2 dðσ d Þεacc q;N :

ð36Þ

Q. Li et al. / Cold Regions Science and Technology 122 (2016) 58–70

Fig. 15. Comparison for accumulated volumetric strain between the test and model.

67

68

Q. Li et al. / Cold Regions Science and Technology 122 (2016) 58–70

Fig. 15 (continued).

low-level repeated loading, and formulations for the constitutive model used in this scheme are also described in detail. In this section, the algorithm of this scheme is presented. According to the three loading steps described previously, this algorithm of the numerical scheme involves three modules.

3.4.1. Module#1: initial stress conditions As previously described, when using this model to calculate large-scale boundary value problems, the initial stress state at the material point level must first be determined. The stress field induced by self-weight and initial static loads are calculated using a standard from the FEM with a selected mechanical constitutive model. Moreover, the static model parameters should be used in the selected constitutive model. At every integration point, the calculated initial stress is σ0, and the corresponding strain can be reset to zero. Thereafter, the explicit functions related to the initial stress at every integration point can be obtained.

3.4.2. Module#2: initial dynamic conditions This module corresponds to loading step 2, which has been defined. As presented in Fig. 12, the pseudo-repeated loading process versus t is monotonous, and the total time size of this loading process is a relative value set as 1.0 in this study. The elastic constitutive model with elastic parameters G0 and K0 is employed to calculate the initial dynamic stress conditions. At every integration point, the stress state at t = 1.0 and repeated stress amplitude calculated from Eq. (25) by integrating the deviatoric stress at t = 1.0 will provide the initial values for the stress integration of the pseudo-repeated loading process versus N.

3.4.3. Module#3: stress integration in the long-term process In this module, the details of the stress integration in a pseudorepeated loading process versus the number of loading cycles N will be presented. The accumulated strain for each element must satisfy compatibility requirements, and the system must be at equilibrium throughout the domain, regardless of whether this frozen soil element

Q. Li et al. / Cold Regions Science and Technology 122 (2016) 58–70

has experienced strain accumulation. As previously mentioned, the functions related to the initial stress state and strength will not be updated in this process. The variables updated from cycle N to cycle N + ΔN are described in Eq. (37):   acc σN ; εN ; εacc N ; ψN → σNþΔN ; εNþΔN ; εNþΔN ; ψNþΔN

ð37Þ

where subscript N designates the material state after the application of N loading cycles (i.e., the last converged state) and subscript N + ΔN designates the new state after the application of N + ΔN load cycles yet to be determined. Vector Ψ contains two components: elastic modulus (i.e., K and G) and accumulated direction ratio (i.e., dg). Based upon the displacement method of the FEM program, Eq. (38) is governed by the update of the total strain: εNþΔN ¼ εN þ Δε

ð38Þ

where the total strain increment Δε is calculated from displacements by geometric equations. The update, Eq. (37), is performed using a fully explicit Euler forward algorithm in which the pressure-accumulationdependent elastic stiffness is included via its resilient values. Eq. (35) shows that the incremental accumulated strain tensor has this form: !

Δεacc ¼

ð39Þ

From Eqs. (20) and (22), the incremental deviatoric accumulated strain Δeacc and incremental volumetric accumulated strain Δεacc v can be obtained. Thereafter, the stress increment induced by the accumulated strain is defined as follows:
 Δσ ¼ Δs þ ΔpI ¼ 2GN ðΔe−Δeacc Þ þ K N Δεv −Δεacc I v

ð40Þ

where KN and GN are the elastic modulus at N loading cycles. The updated accumulated strain tensor and stress tensor are calculated by Eqs. (41) and (42), respectively. acc acc εacc NþΔN ¼ εNþΔN þ Δε

ð41Þ

σNþΔN ¼ σN þ Δσ

ð42Þ

From assumption (c) in Section 3.2, the updated repeated stress amplitude is defined as follows: rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 sNþΔN sNþΔN 2

ð43Þ

where the updated deviatoric stress vector sN + ΔN can be obtained from Eq. (21) with an updated stress vector σN + ΔN. The updated shear modulus and bulk modulus are defined as follows:  GNþΔN ¼ G0 1 þ α 1 εacc v;NþΔN

ð44Þ

 K NþΔN ¼ K 0 1 þ α 2 εacc v;NþΔN

ð45Þ

4. Model verification For this exploratory study, an advanced constitutive model was implemented in the commercial finite element code ABAQUS using UMAT capabilities. The triaxial tests presented in Part I have been considered as numerical examples used to verify and evaluate the performance of this model and its numerical implementation. In these triaxial tests, the frozen soil specimen was subjected to an axisymmetric compressive stress state, including a constant initial stress state and an additional vertical repeated stress. To minimize inaccuracies caused by the step-size dependency, the maximum integration step size was set to ΔN = 1. The parameters that have been determined in Section 2 are summarized in Table 1 and shown to be independent of the number of loading cycles. Comparisons of the accumulated shear strain and accumulated volumetric strain in frozen soil between the experimental results and model results are plotted in Figs. 14 and 15, respectively. The calibrated values of CN1 and αref were within the range listed in Table 1, and each value is presented in the figures. The comparison results from Figs. 14 and 15 indicate that sufficient accuracy is reached in this prediction despite the limited amount of test data. In addition, the accumulated behaviour for frozen soil under long-term low-level repeated loading can be well described by this model using the recommended parameters. However, despite the number of parameters incorporated into this model, most are constant parameters for this frozen soil, and only three other parameters (the parameters k1 in Eq. (17), k2 in Eq. (18) and frozen soil strength at confining pressure of 0.3 MPa s0.3) must be obtained through laboratory testing. 5. Conclusions and discussion

From Eq. (36), the updated accumulated direction ratio is calculated by the updated accumulated shear strain and updated repeated stress amplitude: h i
 dg;NþΔN ¼ dg0 exp κ 2 d σ d;NþΔN εacc q;NþΔN

induced by the accumulated strain. The accuracy of this integration algorithm is dependent on the size of the integration step, with a smaller step size resulting in a more accurate result. However, the amount of required calculations will increase significantly with decreasing integration step sizes. Hence, a balance must be attained between accuracy and computational cost. Therefore, a special procedure focused on solving constitutive equations has been developed. The accumulated shear strain calculated with different integration step sizes can be achieved by this procedure, and the accumulated shear strain can be accurately calculated by the explicit model as shown in Eq. (6), with the accuracy evaluated by comparing the results. Moreover, the computational cost can be evaluated by the number of steps. A test case with an initial stress ratio of 0.5, initial mean stress of 1.5 MPa, repeated stress amplitude of 1.5 MPa and reference strength of 4.25 MPa has been selected. The results of the accumulated shear strain up to 10,000 are presented in Fig. 13. This figure indicates that results with a smaller integration step size are closer to the accurate curve, and ΔN = 1 can provide sufficient accuracy. The computational efficiency with integration step size ΔN = 1 is acceptable.

!

εacc
 3st;N 1 q;re f;N þ dg;N I κ exp − g σ d;N ΔN: 2qt;N 3 C N1

σ d;NþΔN ¼ qd;NþΔN ¼

69

ð46Þ

3.5. Integration step size For simplicity and ease of implementation, a fully explicit integration algorithm (fully Euler-forward scheme) is applied to calculate the stress

A constitutive model is presented here for the elastic and accumulated behaviour of frozen soil under long-term low-level repeated loading. This model is formulated using empirical equations derived from experimental results and classical elasto-plastic theory. In this model, the coupled empirical formulations for accumulated strain and evolution formulations for the elastic modulus as a function of accumulated strain are developed based on experimental evidence in which the significant influence factors presented in Part I are considered. Subsequently, a numerical strategy integrating the presented constitutive model is proposed to predict the permanent deformation or dynamic response under long-term vibration, wherein the repeated loading process is simplified to a constant pseudo-static loading process with respect to a continuous number of repeated loading cycles and the magnitude of

70

Q. Li et al. / Cold Regions Science and Technology 122 (2016) 58–70

constant quasi-static loading is regarded as the amplitude of repeated loading. This model has been implemented in a three-dimensional finite element framework by employing a fully Euler-forward scheme with the commercial finite element code ABAQUS and considering the triaxial tests presented in Part I as numerical examples to evaluate the capability of this advanced constitutive model. The essential features of this model are validated by comparing the simulation results with experimental data. Despite the use of several simplifications, consistent results are obtained between the simulation results and test data for all considered cases. This model is the first attempt to predict the long-term elasto-plastic behaviour under repeated loading for frozen soil. The essential advantages of using this approach include accurate representations of the accumulated and elastic behaviours of frozen soil observed from the experimental data and simulations of stress redistribution behaviour in a long-term repeated loading process. Furthermore, the simple formulations allow for easy comprehension, extension and general applicability. The easy implementation in the finite element code supports the use of this model in boundary value problems. Thus, the presented constitutive model is semi-empirical. If additional influence factors were integrated into this model, then additional parameters would be required, although determining many of these parameters would require a substantial amount of laboratory work for different types of frozen soil. Furthermore, the link between laboratory-determined parameters and in situ conditions is not always indisputable. Consequently, the applicability of such variables in models would be limited in practical applications. Acknowledgments The work presented in this paper is performed within the framework of the National Basic Research Program of China (Grant No. 2012CB026104), National Nature Science Foundation of China (Grant No. 41430634 and 51174261) and the Open Research Fund Program of State Key Laboratory of Frozen Soil Engineering of China (Grant No. SKLFSE201216). The authors are grateful for this financial support. References Behzadi, G., Yandell, W., 1996. Determination of elastic and plastic subgrade soil parameters for asphalt cracking and rutting prediction. Transp. Res. Rec. J. Transp. Res. Board (1540), 97–104. Bouckovalas, G., Whitman, R.V., Marr, W.A., 1984. Permanent displacement of sand with cyclic loading. J. Geotech. Eng. 110 (11), 1606–1623. Boyce, J.R., 1976. The Behaviour of a Granular Material Under Repeated Loading. University of Nottingham . Boyce, H., 1980. A non-linear model for the elastic behaviour of granular materials under repeated loading. Proceedings International Symposium on Soils under Cyclic and Transient Loading, Swansea, UK, pp. 285–294. Brown, S., Hyde, A., 1975. Significance of cyclic confining stress in repeated-load triaxial testing of granular material. Transp. Res. Rec. (537). Chai, J.-C., Miura, N., 2002. Traffic-load-induced permanent deformation of road on soft subsoil. J. Geotech. Geoenviron. 128 (11), 907–916. Dafalias, Y.F., 1986. Bounding surface plasticity. I: mathematical foundation and hypoplasticity. J. Eng. Mech. 112 (9), 966–987.

Drabkin, S., Lacy, H., Kim, D., 1996. Estimating settlement of sand caused by construction vibration. J. Geotech. Eng. 122 (11), 920–928. Drumm, E.C., Reeves, J.S., Madgett, M.R., Trolinger, W.D., 1997. Subgrade resilient modulus correction for saturation effects. J. Geotech. Geoenviron. 123 (7), 663–670. Fall, M., Sawangsuriya, A., Benson, C.H., Edil, T.B., Bosscher, P.J., 2008. On the investigations of resilient modulus of residual tropical gravel lateritic soils from Senegal (West Africa). Geotech. Geol. Eng. 26 (1), 13–35. François, S., Karg, C., Haegeman, W., Degrande, G., 2010. A numerical model for foundation settlements due to deformation accumulation in granular soils under repeated small amplitude dynamic loading. Int. J. Numer. Anal. Methods Geomech. 34 (3), 273–296. Guo, L., Wang, J., Cai, Y., Liu, H., Gao, Y., Sun, H., 2013. Undrained deformation behavior of saturated soft clay under long-term cyclic loading. Soil Dyn. Earthq. Eng. 50, 28–37. Hicks, R.G., Monismith, C.L., 1971. Factors influencing the resilient response of granular materials. Highw. Res. Rec. (345), 15–31. Kaggwa, W.S., Booker, J.R., Carter, J., 1991. Residual strains in calcareous sand due to irregular cyclic loading. J. Geotech. Eng. 117 (2), 201–218. Karg, C., François, S., Haegeman, W., Degrande, G., 2010. Elasto-plastic long-term behavior of granular soils: modelling and experimental validation. Soil Dyn. Earthq. Eng. 30 (8), 635–646. Kolisoja, P., 1997. Resilient Deformation Characteristics of Granular Materials. Tampere University of Technology. Lai, Y., Jin, L., Chang, X., 2009. Yield criterion and elasto-plastic damage constitutive model for frozen sandy soil. Int. J. Plast. 25 (6), 1177–1205. Lai, Y., Xu, X., Yu, W., Qi, J., 2014. An experimental investigation of the mechanical behavior and a hyperplastic constitutive model of frozen loess. Int. J. Eng. Sci. 84, 29–53. Lekarp, F., Dawson, A., 1998. Modelling permanent deformation behaviour of unbound granular materials. Constr. Build. Mater. 12 (1), 9–18. Li, D., Selig, E.T., 1996. Cumulative plastic deformation for fine-grained subgrade soils. J. Geotech. Eng. 122 (12), 1006–1013. Li, D.-W., Fan, J.-H., Wang, R.-H., 2011. Research on visco-elastic-plastic creep model of artificially frozen soil under high confining pressures. Cold Reg. Sci. Technol. 65 (2), 219–225. Li, Q.-l., Ling, X.-z., Wang, L.-n., Zhang, F., Wang, J.-h., Xu, P.-j., 2013. Accumulative strain of clays in cold region under long-term low-level repeated cyclic loading: experimental evidence and accumulation model. Cold Reg. Sci. Technol. 94, 45–52. Marr, W.A., Chrisitian, J.T., 1981. Permanent displacements due to cyclic wave loading. J. Geotech. Eng. Div. 107 (8), 1129–1149. Monismith, C.L., Ogawa, N., Frume, C.R., 1975. Permanent deformation characteristics of subsoil due to repeated loading. Transp. Res. Rec. 537, 1–17. Mroz, Z., Norris, V., Zienkiewicz, O., 1979. Application of an anisotropic hardening model in the analysis of elasto-plastic deformation of soils. Geotechnique 29 (1), 1–34. Niemunis, A., Wichtmann, T., Triantafyllidis, T., 2005. A high-cycle accumulation model for sand. Comput. Geotech. 32 (4), 245–263. Paute, J., Hornych, P., Benaben, J., 1996. Repeated load triaxial testing of granular materials in the French Network of Laboratories des Ponts et Chaussées, FLEXIBLE PAVEMENTS. Proceedings of the European Symposium Euroflex 1993 Held in Lisbon, Portugal, 20–22 September 1993. Prevost, J.H., 1985. A simple plasticity theory for frictional cohesionless soils. Int. J. Soil Dyn. Earthquake Eng. 4 (1), 9–17. Puppala, A.J., Mohammad, L.N., Allen, A., 1999. Permanent deformation characterization of subgrade soils from RLT test. J. Mater. Civ. Eng. 11 (4), 274–282. Sawicki, A., Swidzinski, W., 1989. Mechanics of a sandy subsoil subjected to cyclic loadings. Int. J. Numer. Anal. Methods Geomech. 13 (5), 511–529. Sweere, G.T.H., 1990. Unbound Granular Bases for Roads. TU Delft, Delft University of Technology. Wang, S., Qi, J., Yin, Z., Zhang, J., Ma, W., 2014. A simple rheological element based creep model for frozen soils. Cold Reg. Sci. Technol. 106, 47–54. Wichtmann, T., Niemunis, A., Triantafyllidis, T., 2006. Experimental evidence of a unique flow rule of non-cohesive soils under high-cyclic loading. Acta Geotech. 1 (1), 59–73. Wichtmann, T., Niemunis, A., Triantafyllidis, T., 2014. Flow rule in a high-cycle accumulation model backed by cyclic test data of 22 sands. Acta Geotech. 9 (4), 695–709. Yang, Z., Elgamal, A., Parra, E., 2003. Computational model for cyclic mobility and associated shear deformation. J. Geotech. Geoenviron. 129 (12), 1119–1127. Yang, Y., Lai, Y., Dong, Y., Li, S., 2010. The strength criterion and elastoplastic constitutive model of frozen soil under high confining pressures. Cold Reg. Sci. Technol. 60 (2), 154–160.