# An Improved Nishihara Model for Frozen Loess considering the Influence of Temperature

Authors: Xin Li; Enlong Liu; Bingtang Song; Xingyan Liu
Journal: Advances in Materials Science and Engineering
 (2018)

---

## Abstract

A series of triaxial creep tests under the constant confining pressure are performed on frozen loess specimens, and the creep behavior of the frozen loess with respect to variations in both temperature and deviator stress is examined. Experimental results illustrate that the frozen loess specimens present the attenuation creep at the lower deviatoric stress, whereas the nonattenuation creep under the higher deviatoric stress level, and with a drop in the temperature, the deviator stress value which the exhibition of nonattenuation creep needs will increase under the constant confining pressure condition. According to the microscopic analysis on deformation characteristics in the creep process of frozen soil, both temperature and external stress will cause the hardening and weakening effects, affecting the creep properties of frozen loess. By introducing the hardening variable and damage variable to consider the hardening and weakening effects of the frozen loess, an improved Nishihara model is proposed. The correlations between model parameters and the temperature as well as deviator stress are determined. The comparisons between model predictions and experimental results show that the improved creep constitutive model proposed here can not only describe the whole creep process well, but also reveal the influences of the temperature and deviator stress on the creep behavior of frozen loess, which demonstrate its accuracy and usefulness.


---

## Body

## 1. Introduction

Over the years, with more and more development efforts made in the Qinghai-Tibet Plateau, such as railways, highways, and dams, many engineering problems concerned with frozen soils have arisen [1–4]. Creep is one of the fundamental mechanical behaviors of frozen soil, which has significant impacts on the stability and safety of engineering construction in cold regions. Hence, it is essential to conduct researches on the creep property of frozen soils and establish the proper creep constitutive model so as to provide a reference for engineering design and guarantee the normal service life of the constructions in cold regions.As a kind of four-phase mixture composed of mineral particles, liquid water, ice inclusions, and gaseous inclusions, the creep property of frozen soils is heavily affected by the composition of grains, the temperature, the contents of unfrozen water, and the stress conditions [5]. In view of the external load and the temperature, investigators have carried out some experiments to study the creep behavior of frozen soils [6–10]. As for the microscopic analysis in creep process of frozen soil, Vyalov [11] proposed that the external load led to two opposite phenomena, namely, the hardening and the weakening effects of the frozen soils, whose relative intensity determined the type of creep. Wu et al. [12] dynamically monitored the submicrostructure of frozen soils in uniaxial and triaxial creep tests by computerized tomography and confirmed the theory of Vyalov.In order to properly predict frozen soil deformation over time, investigators have put forward numerous creep constitutive models. Fish [13] viewed the creep deformation of frozen soils as the result of entropy change and then presented an exponential creep constitutive model. Yang et al. [14] conducted triaxial compression tests and triaxial creep tests of frozen silt, found that the increase in confining pressure can result in the deformation suppression, and established the disturbed-state elastoplastic constitutive model and a new creep constitutive model. Hou et al. [15] introduced the effects of hardening and weakening induced from stress on the creep process of frozen soils with different contents of coarse grains and improved the Nishihara model. Based on thermodynamics and damage mechanics, He and Cheng [16] proposed a viscoelastic-plastic constitutive model considering damage. In summary, creep models of geotechnical materials were divided into three types, that is, empirical models, rheological models, and general stress-strain-time models [17]. Among them, the empirical models are widely used for practical engineering due to their simple forms, but their applications are strictly limited to specific boundary and loading conditions; the rheological models have the finite parameters, clear physical meaning, and available three-dimensional expression, and besides, they could be applied to finite element analysis, all of which make them the top priority for investigators while they establish creep constitutive models. However, the number of rheological models with good description for the whole creep process is limited, and the Nishihara model is the one that are relatively general, but the strain rate predicted by the Nishihara model during the tertiary stage is inconsistent with experimental results of frozen soils, for its model prediction value remaining constant. Moreover, there are fewer creep constitutive models available considering the effects of temperature and stress, which reduce the validity and practicability of models.Therefore, it is necessary to acquire information on creep behavior of frozen soils under various temperatures and stress conditions and develop an efficient creep constitutive model which could depict the impacts of those two variations. This paper is organized as follows. First, the triaxial creep tests of frozen loess specimens at the same confining pressure, various temperatures, and deviator stresses are shown and analyzed. Subsequently, the hardening variable and damage variable are introduced to consider the effects of the hardening and weakening in creep process, and then an improved Nishihara model is presented which can produce the reasonable prediction about three creep stages of frozen loess. Furthermore, the model parameters are determined as the functions of temperature and deviator stress. Finally, the improved Nishihara model is verified against the aforementioned triaxial creep tests results of frozen loess, and the dependencies of hardening variable and weakening variable on temperature and deviator stress are discussed as well.


## 2. The Triaxial Creep Test of the Frozen Loess

### 2.1. Test Conditions

The material used in the tests was the loess taken from the Lanzhou province in China, whose physical parameters are listed in Table1. The nature moisture content of the soil was 1.35%, and the liquid limit and plastic index were 27.58% and 19.37%, respectively. The desired amount of loess is put into a cylindrical mold with 61.8 mm in diameter and 125 mm in height. Following this, the specimens were placed into another mold with three same parts and saturated for 12 h after keeping under the vacuum conditions for 3 hours; the measured moisture content of the specimens was 22.35% at that time. After the specimens were laid in a refrigerator to freeze for over 24 h at the temperature of −30°C, they were taken out from the refrigerator and then were put into an incubator at the preset temperature for 24 h. The specimens were then placed into a pressure cell of MTS-810 which can apply the confining pressure and the axial load; the given temperature of the triaxial pressure cell was identical with the incubator. After the confining pressure and axial pressure were applied successively, the triaxial creep test began. The testing temperatures were −15°C, −10°C, and −1.5°C with a precision of ±0.1°C, the confining pressure σ3 was a constant with the value of 0.3 MPa, and the way of applying axial pressure was stress controlled. Before the triaxial creep test, the triaxial compression test should be conducted to determine the axial stress value of the triaxial creep test (ranging from 40% to 100% of compression strength).Table 1
The physical parameters of loess.
Composition of grains (%)
Dry unit weight (g/cm3)
>0.25 mm
0.1–0.25 mm
0.05–0.10 mm
0.0005–0.05 mm
<0.0005 mm
2
3.9
18.1
73.7
2.3
1.7


### 2.2. Test Results and Analyses

A set of creep curves from the tests at various temperatures and deviator stresses are exhibited in Figures1(a)–1(c), from which some conclusions can be reached:(1)
The specimens show both attenuation creep and nonattenuation creep at three temperatures as the deviator stress varies. When the deviator stress is low, the specimens undergo transient stage and steady-state stage, and finally the attenuation creep occurs, whereas the specimens with higher stress exhibit the nonattenuation creep after suffering the transient stage, the steady-state stage, and the tertiary stage in succession. The phenomenon is attributed to the hardening and weakening effects induced by the external load. The external load would result in viscoplastic flow and thawing of ice, followed by movements of mineral particles. On the one hand, those movements would bring about the cracks and structural flaws of frozen loess which mean the weakening effect; on the other hand, they give rise to the recovery of destroyed combination between mineral particles which means the hardening effect. While the hardening effect is more ascendant than the weakening effect all the time, the specimens are manifested as the attenuation creep, whereas they show the nonattenuation creep under the higher stress condition where the weakening effect plays a predominant role in the tertiary stage. The threshold stress is called the long-term strength. In the aforementioned triaxial creep tests of frozen loess, the specimens will present the nonattenuation creep once their applied deviator stress exceeds the corresponding long-term strength at different temperatures.(2)
With an increase in temperature, the dramatic strain appears and it is more likely to suffer from the destructive nonattenuation creep for specimens of frozen loess under the lower deviator stress condition. It is possibly because the temperature is another crucial factor that influences the hardening and weakening effects except for the external load. The pore ice works as the cement between the mineral particles of frozen loess, and as pore ice thaws at warmer temperature, the portion of combination between mineral particles caused by pore ice will fail and the strength of frozen loess will be reduced, which could be regarded as the weakening effect induced by the temperature. Hence, the weakening effect is greater than hardening effect at the lower deviator stress as the temperature rises, which makes the nonattenuation creep for specimens of frozen loess.Figure 1
Axial strain versus time at different temperatures. (a)T=−15°C. (b) T=−10°C. (c) T=−1.5°C.
(a)
(b)
(c)


## 2.1. Test Conditions

The material used in the tests was the loess taken from the Lanzhou province in China, whose physical parameters are listed in Table1. The nature moisture content of the soil was 1.35%, and the liquid limit and plastic index were 27.58% and 19.37%, respectively. The desired amount of loess is put into a cylindrical mold with 61.8 mm in diameter and 125 mm in height. Following this, the specimens were placed into another mold with three same parts and saturated for 12 h after keeping under the vacuum conditions for 3 hours; the measured moisture content of the specimens was 22.35% at that time. After the specimens were laid in a refrigerator to freeze for over 24 h at the temperature of −30°C, they were taken out from the refrigerator and then were put into an incubator at the preset temperature for 24 h. The specimens were then placed into a pressure cell of MTS-810 which can apply the confining pressure and the axial load; the given temperature of the triaxial pressure cell was identical with the incubator. After the confining pressure and axial pressure were applied successively, the triaxial creep test began. The testing temperatures were −15°C, −10°C, and −1.5°C with a precision of ±0.1°C, the confining pressure σ3 was a constant with the value of 0.3 MPa, and the way of applying axial pressure was stress controlled. Before the triaxial creep test, the triaxial compression test should be conducted to determine the axial stress value of the triaxial creep test (ranging from 40% to 100% of compression strength).Table 1
The physical parameters of loess.
Composition of grains (%)
Dry unit weight (g/cm3)
>0.25 mm
0.1–0.25 mm
0.05–0.10 mm
0.0005–0.05 mm
<0.0005 mm
2
3.9
18.1
73.7
2.3
1.7


## 2.2. Test Results and Analyses

A set of creep curves from the tests at various temperatures and deviator stresses are exhibited in Figures1(a)–1(c), from which some conclusions can be reached:(1)
The specimens show both attenuation creep and nonattenuation creep at three temperatures as the deviator stress varies. When the deviator stress is low, the specimens undergo transient stage and steady-state stage, and finally the attenuation creep occurs, whereas the specimens with higher stress exhibit the nonattenuation creep after suffering the transient stage, the steady-state stage, and the tertiary stage in succession. The phenomenon is attributed to the hardening and weakening effects induced by the external load. The external load would result in viscoplastic flow and thawing of ice, followed by movements of mineral particles. On the one hand, those movements would bring about the cracks and structural flaws of frozen loess which mean the weakening effect; on the other hand, they give rise to the recovery of destroyed combination between mineral particles which means the hardening effect. While the hardening effect is more ascendant than the weakening effect all the time, the specimens are manifested as the attenuation creep, whereas they show the nonattenuation creep under the higher stress condition where the weakening effect plays a predominant role in the tertiary stage. The threshold stress is called the long-term strength. In the aforementioned triaxial creep tests of frozen loess, the specimens will present the nonattenuation creep once their applied deviator stress exceeds the corresponding long-term strength at different temperatures.(2)
With an increase in temperature, the dramatic strain appears and it is more likely to suffer from the destructive nonattenuation creep for specimens of frozen loess under the lower deviator stress condition. It is possibly because the temperature is another crucial factor that influences the hardening and weakening effects except for the external load. The pore ice works as the cement between the mineral particles of frozen loess, and as pore ice thaws at warmer temperature, the portion of combination between mineral particles caused by pore ice will fail and the strength of frozen loess will be reduced, which could be regarded as the weakening effect induced by the temperature. Hence, the weakening effect is greater than hardening effect at the lower deviator stress as the temperature rises, which makes the nonattenuation creep for specimens of frozen loess.Figure 1
Axial strain versus time at different temperatures. (a)T=−15°C. (b) T=−10°C. (c) T=−1.5°C.
(a)
(b)
(c)


## 3. An Improved Nishihara Model for Frozen Loess considering the Hardening and Weakening Effect

### 3.1. Hardening Variable and Damage Variable

The variations of the hardening and weakening effects of geological material over time are shown in Figure2 [18], from which it can be clearly noticed that the hardening effect increases gradually with a decreased rate during the transient stage and the steady-state stage and finally stays constant; however, the weakening effect increases linearly at first and evolves into an exponential growth trend during the tertiary creep stage. In this paper, the hardening variable and the damage variable are introduced to take account of the hardening and weakening effects of frozen loess, respectively. The hardening variable stands for the value of hardening effect in the creep process, while the damage variable represents the reduction ratio of material parameters induced by the weakening effect.Figure 2
The illustration of the hardening and weakening curves.The hardening effect can be recognized as the plastic property, so the hardening variableH works in the viscoplastic body. On the basis of Figure 2, there are two following requirements for H should be met: (a) t⟶0, H⟶0 and (b) t⟶∞, H⟶a certain constant. The hardening variable H then can be assumed as(1)Hσ,T=1−2et+1A,where A denotes the level of the hardening effect, and it is closely linked to the temperature and stress level and σ, T, and t represents the stress, the temperature, and the creep duration, respectively.It is generally accepted that the weakening effect has no threshold, and the weakening effect would occur once the external load is applied. It should be noted that the weakening effect becomes remarkable at the end of the steady-state stage and further controls the tertiary stage. Therefore, the damage variable can be viewed as zero for simplification during the transient stage and the steady-state stage. The variableD should meet two requirement as follows: (a) t⟶0, D⟶0 and (b) t⟶∞, D⟶1. The damage variable D then can be expressed as the form of piecewise function, which can be written as(2)Dσ,T=0,q≤σu,1−e−Bt,q>σu,where B reflects the varying speed of weakening effect, which is strongly sensitive to temperature and stress and σu and q are the long-term strength and the deviator stress, respectively.


### 3.2. The Improved Nishihara Model

The Nishihara model [19] consists of the Hook body, the viscoelastic body, and the viscoplastic body as shown in Figure 3(a), where E0 stands for the elastic modulus of the Hook body, E1 and η1 denote, respectively, the elastic modulus and the viscosity coefficient of the viscoelastic body, and σ0 and η2 present, respectively, the yield stress and the viscosity coefficient of the viscoplastic body. It can model the transient and steady-state stages of geomaterial well but not for the tertiary stage. Based on the microscopic analysis on deformation characteristics in the creep process of frozen soil, the hardening variable H and the damage variable D are introduced into the viscoplastic body of the Nishihara model to modify it. The viscosity coefficient of the dashpot in the viscoplastic body of the initial Nishihara model is represented as η2; after the hardening variable Hσ,T and the damage variable Dσ,T are taken into account, it actually reflects effective value which is expressed as η2′. In this way, for the improved Nishihara model as shown in Figure 3(b), the hardening variable is intended for describing the transient and steady-state stages, whereas the damage variable is added with the aim of producing a good description for the rapid deformation growth during the tertiary stage.Figure 3
(a) The Nishihara model. (b) The improved Nishihara model.Combining Equation (1) with Equation (2), the effective viscosity coefficient in the viscoplastic body can be obtained as(3)η2′H,D=η2s1+H1−D=η2s1+1−2et+1A,q≤σu,η2s1+1−2et+1Ae−Bt,q>σu,where η2s is the initial viscosity coefficient of the viscoplastic body, which is assumed to be purely thermodependent.


#### 3.2.1. Formulation of 1-D Creep Constitutive Equations

Whenσ<σu, only the Hook and viscoelastic body are engaged in creep, and the stress-strain relations can be expressed as the following equation:(4a)ε=εe+εve,σ=σe=σveε˙=ε˙e+ε˙ve,σe=E1εve+η1ε˙veσe=E0εe,and furthermore, it could be derived as(4b)σ+η1E0+E1σ˙=E0E1E0+E1ε+E0η1E0+E1ε˙.Whenσ≥σu, all of the three parts of improved model, including the Hook body, the viscoelastic body and viscoplastic body as well, participate in the creep process, and the stress-strain relations can be written as(5a)ε=εe+εve+εvp,σ=σe=σve=σvpε˙=ε˙e+ε˙ve+ε˙vp,σve=E1εve+η1ε˙veσe=E0εe,σvp=σ0+η2′H,Dε˙vp,and it could be further formulated as(5b)E0E1η1η2′H,Dσ−σ0+E0η2′H,D+E0η1+E1η2′H,Dη1η2′H,Dσ˙+σ¨=E0E1η1+E0ε¨.There is such an initial condition for Equation (4b) that ε=0 and ε˙=0 when t=0, and the initial condition is that ε=0, ε˙=0, and ε¨=0 when t=0 for Equation (5b). Solving the differential Equations (4b) and (5b) combined with the initial conditions, the 1-D expression of the improved Nishihara model can be obtained as(6)ε=εe+εve+εvp=σE0+σE11−exp−E1η1t,σ<σ0,σE0+σE11−exp−E1η1t+σ−σ0η2′H,Dt,σ≥σ0.


#### 3.2.2. Formulation of 3-D Creep Constitutive Equations

Letsij=σij−1/3σkkδij, eij=εij−1/3εkkδij, and sij=σij−1/3σkkδij, in which δij is the Kronecker symbol.According to the generalized Hook’s law, 3-D constitutive equation for the Hook body is given by(7)εije=eije=sij2G0,where εije, eije, and G0 represent the strain tensor, deviatoric strain tensor, and the shear modulus in the Hook body, respectively.The equation for the viscoelastic body is derived as follows:(8)εijve=sij2G11−exp−2G1η1t,where εijve and G1 are the strain tensor and the shear modulus in the viscoelastic body, respectively.Zienkiewicz and Cormeau [20] proposed that the viscoplastic strain could be expressed as a function of the yield function and the plastic potential. Based on that theory, the strain rate tensor of the viscoplastic body can be formulated as the following equation:(9)ε˙ijvp=0,F<0,1η2′H,DϕFF0∂Q∂σij,F≥0,where εijvp is the strain tensor of the viscoplastic body; F is the yield function; F0 denotes the reference value of F, which equals to 1.0 MPa for frozen soils; Q is the plastic potential; and ϕ· can be simply written as the function of power whose power equals to 1.By gaining the antiderivative of Equation (9) and then combining it with Equations (7) and (8), the general relation of the improved Nishihara model when the associated flow law is adapted can be obtained finally as(10)εij=εije+εijve+εijvp=sij2G0+sij2G11−exp−2G1η1t,F<0,sij2G0+sij2G11−exp−2G1η1t+1η2′H,DF∂F∂σijt,F≥0.When the hardening and weakening effects are not considered, Equation (10) degenerates into the initial Nishihara model, which is expressed as(11)εij=εije+εijve+εijvp=sij2G0+sij2G11−exp−2G1η1t,F<0,sij2G0+sij2G11−exp−2G1η1t+1η2F∂F∂σijt,F≥0.There are several strength criteria widely used for unfrozen soils, such as Mises criterion, Mohr–coulomb criterion, and Drucker–Prager criterion. However, the failure property of frozen soils is a little bit different from other geomaterials, and it is the four-phase mixture composition of frozen soils that makes their strength surface strongly affected by pressure and phase transition. Fish [21] put forward a parabolic strength criterion which is applicable to frozen soils under a wider range of confining pressure conditions. It is assumed for frozen loess that the initial yield surface is discrepant in size but similar in shape. Therefore, the yield function of frozen loess adopts the similar parabolic form to judge whether the initial yield happens in the improved Nishihara model, which can be expressed as(12)F=3J2−c−σmb+b2aσm2,where σm and J2 represent the first principal stress invariant and the second deviatoric stress invariant and a, b, and c are the material parameters.


#### 3.2.3. Verification Formulation in the Triaxial Creep Tests

Letσ1 denotes the axial stress, σ3 the confining pressure, q the shear stress, and ε1 the axial stress, in triaxial creep tests, we have(13)ε2=ε3,σ2=σ3,q=3J2=σ1−σ3,=σ1−σ3,σm=13σ1+2σ3,s11=23σ1−σ3.The yield function can be rewritten as the following equation:(14)F=σ1−σ3−c−13σ1+2σ3b+b18cσ1+2σ32.Thus, the axial strain is obtained finally as(15)ε1=σ1−σ33G0+σ1−σ33G11−exp−2G1η1t,F<0,σ1−σ33G0+σ1−σ33G11−exp−2G1η1t+1η2′H,DFMt,F≥0,where M=1−1/3b+b/9aσ1+2σ3.


### 3.3. Determination of Model Parameters

The model parameters of the yield functionF can be gotten by fitting the turning point between the line and the curve of stress-strain diagram in the triaxial compression tests, which is recognized as the initial yield point, and the parameters associated with the yield function at various temperatures are listed in Table 2. The long-term strength is fitted via the method of the isochronal curve, which is also listed in Table 2.Table 2
The improved Nishihara model parameters.
Temperature (°C)
c(MPa)
b
a(MPa)
σ
u (MPa)
G
0
/
10
2 (MPa)
−1.5
0.3
0.01
2.5
1.82
0.09
−10
3.86
0.06
6.08
6.04
0.45
−15
5.25
0.09
7.74
8.3
0.75The shear modulus of the Hook body can be gained by triaxial creep tests. According to Equation (10), when t=0, ε11=s11/2G0. Let ε0 denotes the instantaneous axial strain in triaxial creep tests, and the instantaneous deformation is regarded as elasticity for simplification, we get(16)ε0=s112G0=q3G0,and then the shear modulus of the Hook body is equal to the slope of one-third of q and ε0, which is also listed in Table 2.The shear modulus of the viscoelastic body can be assumed as follows:(17)G1=κexp−αq−βT,where κ and α are the material constants; βT is the function with respect to temperature; for frozen loess, it is expressed as βT=0.2009T−T0/1°C−0.0186, and T0 is the reference temperature.The viscosity coefficient of the viscoelastic body can be assumed to be as follows:(18)η1=γexp−δq−θT,where γ and δ are the material constants; θT is the function of temperature; for frozen loess, it can be determined as θT=0.0182T−T0/1°C2+0.3208T−T0/1°C.The model parameterA depicting the level of the strengthening effect can be assumed as the following equation:(19)A=uexp−χq−λT,where u and χ are the material constants; for frozen loess,(20)λT=0.0225T−T01°C2+0.3824T−T01°C.The model parameterB that stands for the varying speed of weakening effect can be assumed as(21)B=ξq−σuσum,where m=1.7+0.41T/T0−0.8T/T02 is based on the aforementioned triaxial creep tests for frozen loess.ξ in Equation (21) can be determined as the following equation:(22)ξ=υexpw¯T,where υ and w¯ are the material constants.For the initial viscosity coefficientη2s in Equation (3), it varies with the temperature and can be assumed as(23)η2s=η20ωT,where η20 is a constant and ωT=0.8245T/T02−0.6767T/T0+0.6983 for frozen loess.


## 3.1. Hardening Variable and Damage Variable

The variations of the hardening and weakening effects of geological material over time are shown in Figure2 [18], from which it can be clearly noticed that the hardening effect increases gradually with a decreased rate during the transient stage and the steady-state stage and finally stays constant; however, the weakening effect increases linearly at first and evolves into an exponential growth trend during the tertiary creep stage. In this paper, the hardening variable and the damage variable are introduced to take account of the hardening and weakening effects of frozen loess, respectively. The hardening variable stands for the value of hardening effect in the creep process, while the damage variable represents the reduction ratio of material parameters induced by the weakening effect.Figure 2
The illustration of the hardening and weakening curves.The hardening effect can be recognized as the plastic property, so the hardening variableH works in the viscoplastic body. On the basis of Figure 2, there are two following requirements for H should be met: (a) t⟶0, H⟶0 and (b) t⟶∞, H⟶a certain constant. The hardening variable H then can be assumed as(1)Hσ,T=1−2et+1A,where A denotes the level of the hardening effect, and it is closely linked to the temperature and stress level and σ, T, and t represents the stress, the temperature, and the creep duration, respectively.It is generally accepted that the weakening effect has no threshold, and the weakening effect would occur once the external load is applied. It should be noted that the weakening effect becomes remarkable at the end of the steady-state stage and further controls the tertiary stage. Therefore, the damage variable can be viewed as zero for simplification during the transient stage and the steady-state stage. The variableD should meet two requirement as follows: (a) t⟶0, D⟶0 and (b) t⟶∞, D⟶1. The damage variable D then can be expressed as the form of piecewise function, which can be written as(2)Dσ,T=0,q≤σu,1−e−Bt,q>σu,where B reflects the varying speed of weakening effect, which is strongly sensitive to temperature and stress and σu and q are the long-term strength and the deviator stress, respectively.


## 3.2. The Improved Nishihara Model

The Nishihara model [19] consists of the Hook body, the viscoelastic body, and the viscoplastic body as shown in Figure 3(a), where E0 stands for the elastic modulus of the Hook body, E1 and η1 denote, respectively, the elastic modulus and the viscosity coefficient of the viscoelastic body, and σ0 and η2 present, respectively, the yield stress and the viscosity coefficient of the viscoplastic body. It can model the transient and steady-state stages of geomaterial well but not for the tertiary stage. Based on the microscopic analysis on deformation characteristics in the creep process of frozen soil, the hardening variable H and the damage variable D are introduced into the viscoplastic body of the Nishihara model to modify it. The viscosity coefficient of the dashpot in the viscoplastic body of the initial Nishihara model is represented as η2; after the hardening variable Hσ,T and the damage variable Dσ,T are taken into account, it actually reflects effective value which is expressed as η2′. In this way, for the improved Nishihara model as shown in Figure 3(b), the hardening variable is intended for describing the transient and steady-state stages, whereas the damage variable is added with the aim of producing a good description for the rapid deformation growth during the tertiary stage.Figure 3
(a) The Nishihara model. (b) The improved Nishihara model.Combining Equation (1) with Equation (2), the effective viscosity coefficient in the viscoplastic body can be obtained as(3)η2′H,D=η2s1+H1−D=η2s1+1−2et+1A,q≤σu,η2s1+1−2et+1Ae−Bt,q>σu,where η2s is the initial viscosity coefficient of the viscoplastic body, which is assumed to be purely thermodependent.


### 3.2.1. Formulation of 1-D Creep Constitutive Equations

Whenσ<σu, only the Hook and viscoelastic body are engaged in creep, and the stress-strain relations can be expressed as the following equation:(4a)ε=εe+εve,σ=σe=σveε˙=ε˙e+ε˙ve,σe=E1εve+η1ε˙veσe=E0εe,and furthermore, it could be derived as(4b)σ+η1E0+E1σ˙=E0E1E0+E1ε+E0η1E0+E1ε˙.Whenσ≥σu, all of the three parts of improved model, including the Hook body, the viscoelastic body and viscoplastic body as well, participate in the creep process, and the stress-strain relations can be written as(5a)ε=εe+εve+εvp,σ=σe=σve=σvpε˙=ε˙e+ε˙ve+ε˙vp,σve=E1εve+η1ε˙veσe=E0εe,σvp=σ0+η2′H,Dε˙vp,and it could be further formulated as(5b)E0E1η1η2′H,Dσ−σ0+E0η2′H,D+E0η1+E1η2′H,Dη1η2′H,Dσ˙+σ¨=E0E1η1+E0ε¨.There is such an initial condition for Equation (4b) that ε=0 and ε˙=0 when t=0, and the initial condition is that ε=0, ε˙=0, and ε¨=0 when t=0 for Equation (5b). Solving the differential Equations (4b) and (5b) combined with the initial conditions, the 1-D expression of the improved Nishihara model can be obtained as(6)ε=εe+εve+εvp=σE0+σE11−exp−E1η1t,σ<σ0,σE0+σE11−exp−E1η1t+σ−σ0η2′H,Dt,σ≥σ0.


### 3.2.2. Formulation of 3-D Creep Constitutive Equations

Letsij=σij−1/3σkkδij, eij=εij−1/3εkkδij, and sij=σij−1/3σkkδij, in which δij is the Kronecker symbol.According to the generalized Hook’s law, 3-D constitutive equation for the Hook body is given by(7)εije=eije=sij2G0,where εije, eije, and G0 represent the strain tensor, deviatoric strain tensor, and the shear modulus in the Hook body, respectively.The equation for the viscoelastic body is derived as follows:(8)εijve=sij2G11−exp−2G1η1t,where εijve and G1 are the strain tensor and the shear modulus in the viscoelastic body, respectively.Zienkiewicz and Cormeau [20] proposed that the viscoplastic strain could be expressed as a function of the yield function and the plastic potential. Based on that theory, the strain rate tensor of the viscoplastic body can be formulated as the following equation:(9)ε˙ijvp=0,F<0,1η2′H,DϕFF0∂Q∂σij,F≥0,where εijvp is the strain tensor of the viscoplastic body; F is the yield function; F0 denotes the reference value of F, which equals to 1.0 MPa for frozen soils; Q is the plastic potential; and ϕ· can be simply written as the function of power whose power equals to 1.By gaining the antiderivative of Equation (9) and then combining it with Equations (7) and (8), the general relation of the improved Nishihara model when the associated flow law is adapted can be obtained finally as(10)εij=εije+εijve+εijvp=sij2G0+sij2G11−exp−2G1η1t,F<0,sij2G0+sij2G11−exp−2G1η1t+1η2′H,DF∂F∂σijt,F≥0.When the hardening and weakening effects are not considered, Equation (10) degenerates into the initial Nishihara model, which is expressed as(11)εij=εije+εijve+εijvp=sij2G0+sij2G11−exp−2G1η1t,F<0,sij2G0+sij2G11−exp−2G1η1t+1η2F∂F∂σijt,F≥0.There are several strength criteria widely used for unfrozen soils, such as Mises criterion, Mohr–coulomb criterion, and Drucker–Prager criterion. However, the failure property of frozen soils is a little bit different from other geomaterials, and it is the four-phase mixture composition of frozen soils that makes their strength surface strongly affected by pressure and phase transition. Fish [21] put forward a parabolic strength criterion which is applicable to frozen soils under a wider range of confining pressure conditions. It is assumed for frozen loess that the initial yield surface is discrepant in size but similar in shape. Therefore, the yield function of frozen loess adopts the similar parabolic form to judge whether the initial yield happens in the improved Nishihara model, which can be expressed as(12)F=3J2−c−σmb+b2aσm2,where σm and J2 represent the first principal stress invariant and the second deviatoric stress invariant and a, b, and c are the material parameters.


### 3.2.3. Verification Formulation in the Triaxial Creep Tests

Letσ1 denotes the axial stress, σ3 the confining pressure, q the shear stress, and ε1 the axial stress, in triaxial creep tests, we have(13)ε2=ε3,σ2=σ3,q=3J2=σ1−σ3,=σ1−σ3,σm=13σ1+2σ3,s11=23σ1−σ3.The yield function can be rewritten as the following equation:(14)F=σ1−σ3−c−13σ1+2σ3b+b18cσ1+2σ32.Thus, the axial strain is obtained finally as(15)ε1=σ1−σ33G0+σ1−σ33G11−exp−2G1η1t,F<0,σ1−σ33G0+σ1−σ33G11−exp−2G1η1t+1η2′H,DFMt,F≥0,where M=1−1/3b+b/9aσ1+2σ3.


## 3.2.1. Formulation of 1-D Creep Constitutive Equations

Whenσ<σu, only the Hook and viscoelastic body are engaged in creep, and the stress-strain relations can be expressed as the following equation:(4a)ε=εe+εve,σ=σe=σveε˙=ε˙e+ε˙ve,σe=E1εve+η1ε˙veσe=E0εe,and furthermore, it could be derived as(4b)σ+η1E0+E1σ˙=E0E1E0+E1ε+E0η1E0+E1ε˙.Whenσ≥σu, all of the three parts of improved model, including the Hook body, the viscoelastic body and viscoplastic body as well, participate in the creep process, and the stress-strain relations can be written as(5a)ε=εe+εve+εvp,σ=σe=σve=σvpε˙=ε˙e+ε˙ve+ε˙vp,σve=E1εve+η1ε˙veσe=E0εe,σvp=σ0+η2′H,Dε˙vp,and it could be further formulated as(5b)E0E1η1η2′H,Dσ−σ0+E0η2′H,D+E0η1+E1η2′H,Dη1η2′H,Dσ˙+σ¨=E0E1η1+E0ε¨.There is such an initial condition for Equation (4b) that ε=0 and ε˙=0 when t=0, and the initial condition is that ε=0, ε˙=0, and ε¨=0 when t=0 for Equation (5b). Solving the differential Equations (4b) and (5b) combined with the initial conditions, the 1-D expression of the improved Nishihara model can be obtained as(6)ε=εe+εve+εvp=σE0+σE11−exp−E1η1t,σ<σ0,σE0+σE11−exp−E1η1t+σ−σ0η2′H,Dt,σ≥σ0.


## 3.2.2. Formulation of 3-D Creep Constitutive Equations

Letsij=σij−1/3σkkδij, eij=εij−1/3εkkδij, and sij=σij−1/3σkkδij, in which δij is the Kronecker symbol.According to the generalized Hook’s law, 3-D constitutive equation for the Hook body is given by(7)εije=eije=sij2G0,where εije, eije, and G0 represent the strain tensor, deviatoric strain tensor, and the shear modulus in the Hook body, respectively.The equation for the viscoelastic body is derived as follows:(8)εijve=sij2G11−exp−2G1η1t,where εijve and G1 are the strain tensor and the shear modulus in the viscoelastic body, respectively.Zienkiewicz and Cormeau [20] proposed that the viscoplastic strain could be expressed as a function of the yield function and the plastic potential. Based on that theory, the strain rate tensor of the viscoplastic body can be formulated as the following equation:(9)ε˙ijvp=0,F<0,1η2′H,DϕFF0∂Q∂σij,F≥0,where εijvp is the strain tensor of the viscoplastic body; F is the yield function; F0 denotes the reference value of F, which equals to 1.0 MPa for frozen soils; Q is the plastic potential; and ϕ· can be simply written as the function of power whose power equals to 1.By gaining the antiderivative of Equation (9) and then combining it with Equations (7) and (8), the general relation of the improved Nishihara model when the associated flow law is adapted can be obtained finally as(10)εij=εije+εijve+εijvp=sij2G0+sij2G11−exp−2G1η1t,F<0,sij2G0+sij2G11−exp−2G1η1t+1η2′H,DF∂F∂σijt,F≥0.When the hardening and weakening effects are not considered, Equation (10) degenerates into the initial Nishihara model, which is expressed as(11)εij=εije+εijve+εijvp=sij2G0+sij2G11−exp−2G1η1t,F<0,sij2G0+sij2G11−exp−2G1η1t+1η2F∂F∂σijt,F≥0.There are several strength criteria widely used for unfrozen soils, such as Mises criterion, Mohr–coulomb criterion, and Drucker–Prager criterion. However, the failure property of frozen soils is a little bit different from other geomaterials, and it is the four-phase mixture composition of frozen soils that makes their strength surface strongly affected by pressure and phase transition. Fish [21] put forward a parabolic strength criterion which is applicable to frozen soils under a wider range of confining pressure conditions. It is assumed for frozen loess that the initial yield surface is discrepant in size but similar in shape. Therefore, the yield function of frozen loess adopts the similar parabolic form to judge whether the initial yield happens in the improved Nishihara model, which can be expressed as(12)F=3J2−c−σmb+b2aσm2,where σm and J2 represent the first principal stress invariant and the second deviatoric stress invariant and a, b, and c are the material parameters.


## 3.2.3. Verification Formulation in the Triaxial Creep Tests

Letσ1 denotes the axial stress, σ3 the confining pressure, q the shear stress, and ε1 the axial stress, in triaxial creep tests, we have(13)ε2=ε3,σ2=σ3,q=3J2=σ1−σ3,=σ1−σ3,σm=13σ1+2σ3,s11=23σ1−σ3.The yield function can be rewritten as the following equation:(14)F=σ1−σ3−c−13σ1+2σ3b+b18cσ1+2σ32.Thus, the axial strain is obtained finally as(15)ε1=σ1−σ33G0+σ1−σ33G11−exp−2G1η1t,F<0,σ1−σ33G0+σ1−σ33G11−exp−2G1η1t+1η2′H,DFMt,F≥0,where M=1−1/3b+b/9aσ1+2σ3.


## 3.3. Determination of Model Parameters

The model parameters of the yield functionF can be gotten by fitting the turning point between the line and the curve of stress-strain diagram in the triaxial compression tests, which is recognized as the initial yield point, and the parameters associated with the yield function at various temperatures are listed in Table 2. The long-term strength is fitted via the method of the isochronal curve, which is also listed in Table 2.Table 2
The improved Nishihara model parameters.
Temperature (°C)
c(MPa)
b
a(MPa)
σ
u (MPa)
G
0
/
10
2 (MPa)
−1.5
0.3
0.01
2.5
1.82
0.09
−10
3.86
0.06
6.08
6.04
0.45
−15
5.25
0.09
7.74
8.3
0.75The shear modulus of the Hook body can be gained by triaxial creep tests. According to Equation (10), when t=0, ε11=s11/2G0. Let ε0 denotes the instantaneous axial strain in triaxial creep tests, and the instantaneous deformation is regarded as elasticity for simplification, we get(16)ε0=s112G0=q3G0,and then the shear modulus of the Hook body is equal to the slope of one-third of q and ε0, which is also listed in Table 2.The shear modulus of the viscoelastic body can be assumed as follows:(17)G1=κexp−αq−βT,where κ and α are the material constants; βT is the function with respect to temperature; for frozen loess, it is expressed as βT=0.2009T−T0/1°C−0.0186, and T0 is the reference temperature.The viscosity coefficient of the viscoelastic body can be assumed to be as follows:(18)η1=γexp−δq−θT,where γ and δ are the material constants; θT is the function of temperature; for frozen loess, it can be determined as θT=0.0182T−T0/1°C2+0.3208T−T0/1°C.The model parameterA depicting the level of the strengthening effect can be assumed as the following equation:(19)A=uexp−χq−λT,where u and χ are the material constants; for frozen loess,(20)λT=0.0225T−T01°C2+0.3824T−T01°C.The model parameterB that stands for the varying speed of weakening effect can be assumed as(21)B=ξq−σuσum,where m=1.7+0.41T/T0−0.8T/T02 is based on the aforementioned triaxial creep tests for frozen loess.ξ in Equation (21) can be determined as the following equation:(22)ξ=υexpw¯T,where υ and w¯ are the material constants.For the initial viscosity coefficientη2s in Equation (3), it varies with the temperature and can be assumed as(23)η2s=η20ωT,where η20 is a constant and ωT=0.8245T/T02−0.6767T/T0+0.6983 for frozen loess.


## 4. Validation of the Improved Model

The improved Nishihara model is verified through the aforementioned experimental results obtained on frozen loess, whose parameters areT0=−10°C, κ=190MPa, α=0.31MPa−1, γ=1473.3MPa·h, δ=0.529MPa−1, u=439.25, χ=0.64MPa−1, v=772.9, w¯=0.362°C−1, and η20=1300MPa·h. Figures 4(a)–4(c) show the tested and predicted results of triaxial creep cases conducted at −15°C, −10°C, and −1.5°C, respectively.Figure 4
Comparison between tested results and predicted results. (a)T=−15°C. (b) T=−10°C. (c) T=−1.5°C.
(a)
(b)
(c)It can be drawn from Figures4(a)–4(c) that the improved model is able to simulate both attenuation creep and nonattenuation creep and there is a good agreement with the experimental results, which demonstrate that the model can reproduce effects of the temperature and deviator stress on creep behavior of frozen loess.Through performing the derivative of Equation (11), the strain rate of the initial Nishihara model is derived as(24)ε˙ij=sijη1exp−2G1η1t,F<0,sijη1exp−2G1η1t+1η2F∂F∂σij,F≥0.Based on the Equation (24), the strain rate during the tertiary stage will approach to a constant that equals to 1/η2F∂F/∂σij when t⟶∞.Taking the case of frozen loess at the constant temperature of −1.5°C and the constant deviator stress of 1.84 MPa as an example, the compared result between the improved Nishihara model and the initial one (as shown in Figure 5) implies that the improved model is more suitable for describing the tertiary stage of frozen loess during which the strain grows exponentially.Figure 5
Compared result between the improved model and the initial one.


## 5. Discussion

Figures6 and 7 exhibit the varying trends of the hardening variable and the damage variable of frozen loess incurring under various temperatures and deviator stresses. Observations can be made from Figures 6 and 7: for the attenuation creep and the nonattenuation creep, the hardening variable H standing for the value of the hardening effect grows to a constant over time, and the constant is negatively correlated with temperature and deviator stress; the damage variable D presenting the reduction ratio of material parameters induced by the weakening effect is considered only when the nonattenuation creep occurs, and it will increase to 1 over time; the higher the temperature and deviator stress, the more rapidly the damage variable will reach to 1. These indicate that the hardening effect will drop as the deviator stress becomes higher and the temperature turns warmer, whereas the weakening effect will increase with the rises in temperature and deviator stress; that is, for frozen loess, when the temperature and deviator stress increase, the hardening effect will decrease and the weakening effect will increase; thus the deformation is bigger, and it will be more likely to result in the nonattenuation creep, which is identical to strain-time curve of frozen loess in the triaxial creep tests under the conditions where the confining pressure is the same and the temperature and the deviator stress are various.Figure 6
The changes of the hardening variable with temperature and stress.
(a)
(b)Figure 7
The changes of the damage variable with temperature and stress.
(a)
(b)


## 6. Conclusions

On the basis of the triaxial creep test results of frozen loess, the improved Nishihara model considering the combining influences of temperature and deviator stress is proposed and validated, and the following conclusions can be drawn:(1)
The frozen loess specimens show both attenuation creep and nonattenuation creep at different experimental temperatures, which implies that the deviator stress heavily affect the hardening and weakening effects; therefore, the two phenomena that the hardening effect is always greater than the weakening effect and the weakening effect is less than the hardening effect during the tertiary stage will happen ignoring the testing temperature, which leads to two kinds of creep appearance. When the temperature is higher, the bigger deformation is more likely to occur for specimens under the lower deviator stress condition, and the nonattenuation creep arises, it is because the rise in temperature causes the ice to melt and the weakening effect is enhanced.(2)
According to the concept of hardening and weakening effects, the hardening variable and the damage variable are introduced to improve the Nishihara model, and the model parameters are determined as functions of deviator stress and temperature. The verification result confirms that the improved model can simulate the whole creep process well, including the transient stage, the steady-state stage, and the tertiary stage. There is the application for the improved model to predict the creep deformation of frozen loess subjected to the constant confining pressure and different temperatures and various deviator stresses.



---

*Source: 9073435-2018-11-01.xml*