# Study on Nonlinear Damage Creep Model for Rocks under Cyclic Loading and Unloading

Authors: Liangliang Zhang; Xiaojian Wang
Journal: Advances in Materials Science and Engineering
 (2021)

---

## Abstract

To determine the nonlinear creep characteristics of rocks under cyclic loading and unloading conditions, a nonlinear Kelvin model and damage viscoplastic model are proposed. The models are connected in series with a linear elastic body to establish a nonlinear damage creep model. The differential damage constitutive equations of the proposed creep model under one-dimensional and three-dimensional stress states are derived based on the creep mechanics and elasticity theory. The damage and unloading creep equations are then obtained based on the superposition principle, and a simple and feasible method for determining the model parameters is determined. Finally, the step cyclic loading and unloading creep test data for lherzolite and limestone are used to verify the rationality and feasibility of the nonlinear damage creep model. The results show that the theoretical creep curves of the nonlinear damage creep model are consistent with the experimental curves which indicates that the proposed model can not only determine the creep properties of lherzolite and limestone under cyclic loading and unloading but also determine the nonlinear characteristics of rocks in the transient and steady-state creep stages and particularly within the accelerating creep stage.


---

## Body

## 1. Introduction

Creep characteristics are intrinsic mechanical properties of rocks, which are significant for the analysis of various rock engineering failure problems, such as dam foundation stability, reservoir subsidence, and tunnel support design [1]. Particularly, deep underground excavation activities frequently affect the stability of surrounding rocks. Moreover, the perturbed surrounding rocks, whose creep characteristics cannot be ignored, are subjected to loading and unloading alternatively because of excavation [2, 3].Currently, there are hundreds of creep models, which can be classified into three types: empirical [4], component [5], and damage models [6, 7]. The empirical creep model is a mathematical expression established by fitting the test results. The empirical creep equation has the advantages of a simple form and few creep parameters. However, predicting the long-term deformation of salt rocks using empirical models tends to produce large errors. This is because the models cannot reflect the creep mechanical mechanism of salt rocks by purely fitting mathematically, and the parameters of the creep model lack physical significance [8, 9]. The component models are constructed by connecting the components with basic functions (including elastic, plastic, and viscosity components) in series and parallel. These models can directly reflect the complex mechanical properties of rocks, their parameters have clear physical significance, and the creep equations are easy to establish. However, the model parameters of component models are constant under different stress levels; thus, they cannot reflect the nonlinear characteristics of rocks under cyclic loading and unloading. In most cases, the theoretical results of component models are not consistent with the experimental data [10, 11]. The damage creep model can determine the deterioration damage effect of rocks during cyclic loading and unloading, and its creep model parameters vary with stress level and time. Therefore, the model can reflect the nonlinear creep characteristics of the rock.In this study, a novel nonlinear Kelvin body and viscoplastic damage body are established and are then connected in series with the elastic body to obtain a novel nonlinear damage creep model. This model effectively determines the loading and unloading creep properties of the lherzolite and limestone under different stress levels and also favorably represents the nonlinear characteristics of accelerating creep, overcoming the shortcomings of the empirical creep models and the component models.


## 2. Establishment of the Nonlinear Damage Creep Model

In the rock step cyclic loading and unloading compression creep test, the total strain can be divided into instantaneous elastic strain, instantaneous plastic strain, viscoelastic strain, and viscoplastic strain [12, 13] (Figure 1). The total strain can then be expressed as follows:(1)ε=εm+εv=εme+εve+εrpεm=εme+εmpεv=εve+εvpεrp=εmp+εvp,where ε is the total strain of rock under uniaxial compression creep state, εmandεv are the instantaneous strain and creep strain, respectively, and εme,εmp,εve,εvp,andεrp are the instantaneous elastic strain, instantaneous plastic strain, viscoelastic strain, viscoplastic strain, and residual plastic strain, respectively.Figure 1
Typical creep curve during the loading and unloading test.The instantaneous plastic strain only accounts for 15%–20% of the instantaneous strain, which is smaller than the instantaneous elastic strain. Moreover, the increment of the instantaneous plastic strain per unit stress decreases with the increase in stress level. Simultaneously, if the effect of the instantaneous plastic strain is considered, the creep mechanical model and creep equation of the rock will be more complex, which is not conducive for application. Therefore, the novel nonlinear damage creep model established in this study does not consider the effect of instantaneous plastic strain [14]. The relationship between the stress and instantaneous elastic strain can be described by the linear elastic body.


### 2.1. Nonlinear Kelvin Model

The transient and unloading creep curves have obvious nonlinear characteristics, which are difficult to accurately determine using the traditional Kelvin model [15]. This is because the parameters of the Kelvin model are constant and do not change with creep time. Therefore, it is necessary to establish a nonlinear Kelvin model to simultaneously determine the viscoelastic creep characteristics of the rocks during loading and unloading. To analyze the nonlinear creep characteristics of rocks, many researchers have considered the viscous coefficient of the Kelvin model as a time-dependent function [16]. In this study, a nonlinear Kelvin model is established based on the assumption that the viscous coefficient satisfies a power function with time during the creep test. The mechanical model is shown in Figure 2.Figure 2
Nonlinear Kelvin model.The differential constitutive equation of the nonlinear Kelvin model can be expressed as(2)σ=Enkεve+ηnkt1−λε˙ve,where σ is the applied stress; Enkandηnk represent the elastic modulus and viscous coefficient of the nonlinear Kelvin model, respectively; λ is the material constant, which is related to the stress level; and t is the creep time.By solving equation (2), the creep equation of the nonlinear Kelvin model is obtained as follows:(3)εve=σEnk1−exp−Enkηnkλtλ.Assuming the applied stress is unloaded att0, equation (2) can be expressed as(4)ε˙ve+Enkηnkt01−λεve=0.Equation (5) can be obtained by solving the differential equation (4) as follows:(5)εve=Aexp−Enkηnkt01−λt,where A is the integral constant.Whent=t0, εve=σ/Enk1−exp−Enk/ηnkλt0λ, and this boundary condition is substituted into equation (5) to obtain the integral constant as follows:(6)A=σEnk1−exp−Enkηnkλt0λexpEnkηnkt01−λt0.By substituting equation (6) into equation (5), the creep strain εul during the unloading stage can be expressed as(7)εul=σEnk1−exp−Enkηnkλt0λexpEnkηnkt01−λt0−t.


### 2.2. Damage Viscoplastic Model

When the applied stress exceeds the long-term strength of the rock, it undergoes unstable creep, in which the strain and strain rate increases rapidly with the creep time. The rock is destroyed within a short period. Because the unstable creep has no unloading process, the strain generated in the unstable creep stage is regarded as viscoplastic strain [17]. The law governing the change between the stress and viscoplastic strain in the accelerated creep stage can be determined by the damage viscoplastic model. The viscosity coefficient, ηvp, in the nondamaged state is replaced by the effective viscosity coefficient, ηvp1−D, to determine the deterioration process of the rock. The mechanical model is shown in Figure 3.Figure 3
Damage viscoplastic model.Based on Lemaitre’s strain equivalence principle [18, 19], the stress-strain relation of the damage viscous body can be expressed as follows:(8)ε˙vp=σ˜ηvp=σηvp1−D,where σ˜ is the effective stress and D is the damage variable.The damage evolution equation proposed by Kachanov [20] is widely used to determine the damage behavior of rocks during the accelerated creep stage.(9)D˙=Cσ1−Dn,where C and n represent the rock material constants.The creep damage critical failure time can be derived using equation (9) as follows:(10)tf=C1+nσn−1,where tf represents the critical failure time.The relationship between the damage variableD and time can be obtained using equations (9) and (10) as follows:(11)D=1−1−ttf1/1=n.We can obtain the constitutive equation of the damaged viscous body by substituting equation (11) into equation (8) as follows:(12)ε˙vp=σηvp1−ttf−1/1+n.By integrating equation (12), we can obtain the following equation:(13)εvp=−σηvptf1+nn1−ttfn/1+n+B,where B is another integral constant.Whent=0, ε=0, and this boundary condition is substituted into equation (13) to obtain the integral constant as follows:(14)B=σηvptf1+nn.By substituting equation (14) into equation (13), we obtain(15)εvp=σηvptf1+nn1−1−ttfn/1+n=σηvptf1+nn1−1−Dn.By substituting equation (10) into equation (15), we obtain(16)εvp=σ1−nnCηvp1−1−Dn.In previous research, the creep equation of the damage viscous body was obtained using the Kachanov creep damage rate as follows:(17)εvp=σηvp1−Dt.By calculating the derivative of equation (17) with respect to time, the relationship between the stress and strain rate of the damage viscous body can be expressed as follows:(18)ε˙vp=σηvp11−Dσ,t+Cσn1−Dσ,tn+2t≠σηvp1−D.Comparing equation (18) and equation (9), the relationship between the stress and strain rate does not satisfy the Lemaitre strain equivalence principle; therefore, it is less rigorous to directly replace the stress in the creep equation with the effective stress in theory.For a damage viscoplastic model under the one-dimensional stress state, the creep equation is obtained by replacing stressσ with stress σ−σs.(19)εvp=σ−σsηvptf1+nn1−1−ttfn/1+n,where σs is the long-term strength, which can be obtained by the creep test.


### 2.3. Nonlinear Damage Creep Model

The instantaneous elastic body determining the instantaneous strain, the nonlinear Kelvin model determining the viscoelastic strain, and the damage viscoplastic model determining the viscoplastic strain are connected in series to form a nonlinear damage creep model, as shown in Figure4.Figure 4
Nonlinear damage creep model.Based on the stress-strain relationship of the model in the series-parallel connection, we can obtain the following equations:(20)σ=σe=σve=σvpε=εe+εve+εvp,where σandε are the total stress and strain, respectively, σe,σve,andσvp and εe,εve,andεvp represent the three-part stress and strain of the nonlinear damage creep model, respectively.The constitutive equations for each part are expressed as(21)σ=σe=Eeεeσ=σve=Enkεve+ηnkt1−λε˙veσ=σvp=σS+ηvp1−Dε˙vp,where EeandEnk are the elastic modulus, ηnkandηvp are the viscosity coefficients, and ε˙veandε˙vp are the strain rates.When it is one-dimensional, the constitutive equation is as follows:(1)
whenσ<σs,(22)Ee+EnkEeηnkt1−λσ+1Eeσ˙=ε˙+Enkηnkt1−λε.(2)
whenσs≤σ,(23)σ¨Ee+Ee+Enk+ηnk1−λt−λEeηnkt1−λσ˙+Enk+ηnk1−λt−ληnkηvpt1−λ1−Dσ−σs+σ˙ηvp1−D+D˙σ−σsηvp1−D2=ε¨+Enk+ηnk1−λt−ληnkt1−λε˙.The axial creep equations of the nonlinear damage rheological model can be obtained based on the superposition principle.(1)
Ifσ<σS, then(24)ε=σEe+σEnk1−exp−Enkηnkλtλ.(2)
IfσS≤σ, then(25)ε=σEe+σEnk1−exp−Enkηnkλtλ+σ−σsηvptf1+nn1−1−ttfn/1+n.


### 2.4. Three-Dimensional Creep Equation

Because actual rock mass engineering often involves rocks in a complex three-dimensional stress condition, it is necessary to establish the axial creep and unloading equations of the nonlinear damage creep model under three-dimensional stress conditions to reflect the creep properties of rocks more accurately. Assuming that the rock specimen is a continuous and uniform material, the total strain of rock under triaxial compression creep state,εij, of the nonlinear damage creep model under three-dimensional stress conditions can be expressed as(26)εij=εije+εijve+εijvp.For an instantaneous elastic body, the strain under three-dimensional stress conditions is obtained according to the general Hooke's law [21] as(27)εije=12GeSij+13Keσmδij,where GeandKe are the shear modulus and bulk modulus, respectively, Sijandσmδij are the stress tensor and spherical stress tensor, respectively, and σm is the mean stress.For the nonlinear Kelvin model, the creep equation under three-dimensional stress conditions can refer to the creep equation under one-dimensional stress conditions as follows:(28)εijve=Sij2Gnk1−exp−Gnkηnkλtλ,where Gnk is the shear modulus of the nonlinear Kelvin model.The three-dimensional constitutive relationship of the damage viscoplastic model is(29)εijvp=tf1+n1−1−Dn2ηvpnΦFF0∂Q∂σij,where • is the switch function, F is the rock yield function, F0 is the initial value of the rock yield function, which is typically 1 for simplicity [22], and ΦF/F0 is generally an exponential or power function. Here, we assume a power function, and the power exponent is set to 1 for the rock materials [23]. ΦF/F0=F/F0. Q is the plastic potential function. Moreover, for convenience of calculation, the associated flow rule is used; therefore, F=Q [24].The Mohr–Coulomb and Von Mises yield criteria are the most applicable rock yield criteria. However, the Mohr–Coulomb criterion does not consider the influence of the intermediate principal stress on the yield failure of the rocks, and the Von Mises yield criterion ignores the influence of the spherical stress on the creep properties of the rocks, especially soft rocks. Therefore, the Drucker–Prager yield criterion [25] is adopted in this study. This yield criterion can effectively compensate for the shortcomings of the first two yielding criteria and is expressed as(30)F=J2−αI1−k,where J2 is the second deviatoric stress invariant, I1 is the first invariant of the stress tensor, and αandk are the material parameters expressed as follows:(31)α=sinφ33+sin2φk=3ccosφ3+sin2φ.Here,φ is the internal friction angle of the rock material and c is the sticky cohesion of the rock material.The following equations can be obtained in the conventional triaxial step cycle loading and unloading compression creep test:(32)σ1>σ2=σ3σm=σ1+σ2+σ33=σ1+2σ33J2=σ1−σ33I1=σ1+σ2+σ3=σ1+2σ3S11=2σ1−σ33∂F∂σ11=33−α.According to equations (26)–(32), the three-dimensional creep equations of the nonlinear damage creep model are expressed as follows:(33)ε=σ1+2σ39Ke+σ1−σ33Ge+σ1−σ33Gnk1−exp−Gnkηnkλtλσ<σsσ1+2σ39Ke+σ1−σ33Ge+σ1−σ33Gnk1−exp−Gnkηnkλtλ+3−3α3tf1+n1−1−Dn2ηvpnσ1−σ3−σsσs≤σ.Assuming the applied stress is unloaded att0, the unloading equations of the nonlinear damage creep model are expressed as(34)ε=σ1−σ33Gnk1−exp−Gnkηnkλt0λexpGnkηnkt01−λt0−t.


## 2.1. Nonlinear Kelvin Model

The transient and unloading creep curves have obvious nonlinear characteristics, which are difficult to accurately determine using the traditional Kelvin model [15]. This is because the parameters of the Kelvin model are constant and do not change with creep time. Therefore, it is necessary to establish a nonlinear Kelvin model to simultaneously determine the viscoelastic creep characteristics of the rocks during loading and unloading. To analyze the nonlinear creep characteristics of rocks, many researchers have considered the viscous coefficient of the Kelvin model as a time-dependent function [16]. In this study, a nonlinear Kelvin model is established based on the assumption that the viscous coefficient satisfies a power function with time during the creep test. The mechanical model is shown in Figure 2.Figure 2
Nonlinear Kelvin model.The differential constitutive equation of the nonlinear Kelvin model can be expressed as(2)σ=Enkεve+ηnkt1−λε˙ve,where σ is the applied stress; Enkandηnk represent the elastic modulus and viscous coefficient of the nonlinear Kelvin model, respectively; λ is the material constant, which is related to the stress level; and t is the creep time.By solving equation (2), the creep equation of the nonlinear Kelvin model is obtained as follows:(3)εve=σEnk1−exp−Enkηnkλtλ.Assuming the applied stress is unloaded att0, equation (2) can be expressed as(4)ε˙ve+Enkηnkt01−λεve=0.Equation (5) can be obtained by solving the differential equation (4) as follows:(5)εve=Aexp−Enkηnkt01−λt,where A is the integral constant.Whent=t0, εve=σ/Enk1−exp−Enk/ηnkλt0λ, and this boundary condition is substituted into equation (5) to obtain the integral constant as follows:(6)A=σEnk1−exp−Enkηnkλt0λexpEnkηnkt01−λt0.By substituting equation (6) into equation (5), the creep strain εul during the unloading stage can be expressed as(7)εul=σEnk1−exp−Enkηnkλt0λexpEnkηnkt01−λt0−t.


## 2.2. Damage Viscoplastic Model

When the applied stress exceeds the long-term strength of the rock, it undergoes unstable creep, in which the strain and strain rate increases rapidly with the creep time. The rock is destroyed within a short period. Because the unstable creep has no unloading process, the strain generated in the unstable creep stage is regarded as viscoplastic strain [17]. The law governing the change between the stress and viscoplastic strain in the accelerated creep stage can be determined by the damage viscoplastic model. The viscosity coefficient, ηvp, in the nondamaged state is replaced by the effective viscosity coefficient, ηvp1−D, to determine the deterioration process of the rock. The mechanical model is shown in Figure 3.Figure 3
Damage viscoplastic model.Based on Lemaitre’s strain equivalence principle [18, 19], the stress-strain relation of the damage viscous body can be expressed as follows:(8)ε˙vp=σ˜ηvp=σηvp1−D,where σ˜ is the effective stress and D is the damage variable.The damage evolution equation proposed by Kachanov [20] is widely used to determine the damage behavior of rocks during the accelerated creep stage.(9)D˙=Cσ1−Dn,where C and n represent the rock material constants.The creep damage critical failure time can be derived using equation (9) as follows:(10)tf=C1+nσn−1,where tf represents the critical failure time.The relationship between the damage variableD and time can be obtained using equations (9) and (10) as follows:(11)D=1−1−ttf1/1=n.We can obtain the constitutive equation of the damaged viscous body by substituting equation (11) into equation (8) as follows:(12)ε˙vp=σηvp1−ttf−1/1+n.By integrating equation (12), we can obtain the following equation:(13)εvp=−σηvptf1+nn1−ttfn/1+n+B,where B is another integral constant.Whent=0, ε=0, and this boundary condition is substituted into equation (13) to obtain the integral constant as follows:(14)B=σηvptf1+nn.By substituting equation (14) into equation (13), we obtain(15)εvp=σηvptf1+nn1−1−ttfn/1+n=σηvptf1+nn1−1−Dn.By substituting equation (10) into equation (15), we obtain(16)εvp=σ1−nnCηvp1−1−Dn.In previous research, the creep equation of the damage viscous body was obtained using the Kachanov creep damage rate as follows:(17)εvp=σηvp1−Dt.By calculating the derivative of equation (17) with respect to time, the relationship between the stress and strain rate of the damage viscous body can be expressed as follows:(18)ε˙vp=σηvp11−Dσ,t+Cσn1−Dσ,tn+2t≠σηvp1−D.Comparing equation (18) and equation (9), the relationship between the stress and strain rate does not satisfy the Lemaitre strain equivalence principle; therefore, it is less rigorous to directly replace the stress in the creep equation with the effective stress in theory.For a damage viscoplastic model under the one-dimensional stress state, the creep equation is obtained by replacing stressσ with stress σ−σs.(19)εvp=σ−σsηvptf1+nn1−1−ttfn/1+n,where σs is the long-term strength, which can be obtained by the creep test.


## 2.3. Nonlinear Damage Creep Model

The instantaneous elastic body determining the instantaneous strain, the nonlinear Kelvin model determining the viscoelastic strain, and the damage viscoplastic model determining the viscoplastic strain are connected in series to form a nonlinear damage creep model, as shown in Figure4.Figure 4
Nonlinear damage creep model.Based on the stress-strain relationship of the model in the series-parallel connection, we can obtain the following equations:(20)σ=σe=σve=σvpε=εe+εve+εvp,where σandε are the total stress and strain, respectively, σe,σve,andσvp and εe,εve,andεvp represent the three-part stress and strain of the nonlinear damage creep model, respectively.The constitutive equations for each part are expressed as(21)σ=σe=Eeεeσ=σve=Enkεve+ηnkt1−λε˙veσ=σvp=σS+ηvp1−Dε˙vp,where EeandEnk are the elastic modulus, ηnkandηvp are the viscosity coefficients, and ε˙veandε˙vp are the strain rates.When it is one-dimensional, the constitutive equation is as follows:(1)
whenσ<σs,(22)Ee+EnkEeηnkt1−λσ+1Eeσ˙=ε˙+Enkηnkt1−λε.(2)
whenσs≤σ,(23)σ¨Ee+Ee+Enk+ηnk1−λt−λEeηnkt1−λσ˙+Enk+ηnk1−λt−ληnkηvpt1−λ1−Dσ−σs+σ˙ηvp1−D+D˙σ−σsηvp1−D2=ε¨+Enk+ηnk1−λt−ληnkt1−λε˙.The axial creep equations of the nonlinear damage rheological model can be obtained based on the superposition principle.(1)
Ifσ<σS, then(24)ε=σEe+σEnk1−exp−Enkηnkλtλ.(2)
IfσS≤σ, then(25)ε=σEe+σEnk1−exp−Enkηnkλtλ+σ−σsηvptf1+nn1−1−ttfn/1+n.


## 2.4. Three-Dimensional Creep Equation

Because actual rock mass engineering often involves rocks in a complex three-dimensional stress condition, it is necessary to establish the axial creep and unloading equations of the nonlinear damage creep model under three-dimensional stress conditions to reflect the creep properties of rocks more accurately. Assuming that the rock specimen is a continuous and uniform material, the total strain of rock under triaxial compression creep state,εij, of the nonlinear damage creep model under three-dimensional stress conditions can be expressed as(26)εij=εije+εijve+εijvp.For an instantaneous elastic body, the strain under three-dimensional stress conditions is obtained according to the general Hooke's law [21] as(27)εije=12GeSij+13Keσmδij,where GeandKe are the shear modulus and bulk modulus, respectively, Sijandσmδij are the stress tensor and spherical stress tensor, respectively, and σm is the mean stress.For the nonlinear Kelvin model, the creep equation under three-dimensional stress conditions can refer to the creep equation under one-dimensional stress conditions as follows:(28)εijve=Sij2Gnk1−exp−Gnkηnkλtλ,where Gnk is the shear modulus of the nonlinear Kelvin model.The three-dimensional constitutive relationship of the damage viscoplastic model is(29)εijvp=tf1+n1−1−Dn2ηvpnΦFF0∂Q∂σij,where • is the switch function, F is the rock yield function, F0 is the initial value of the rock yield function, which is typically 1 for simplicity [22], and ΦF/F0 is generally an exponential or power function. Here, we assume a power function, and the power exponent is set to 1 for the rock materials [23]. ΦF/F0=F/F0. Q is the plastic potential function. Moreover, for convenience of calculation, the associated flow rule is used; therefore, F=Q [24].The Mohr–Coulomb and Von Mises yield criteria are the most applicable rock yield criteria. However, the Mohr–Coulomb criterion does not consider the influence of the intermediate principal stress on the yield failure of the rocks, and the Von Mises yield criterion ignores the influence of the spherical stress on the creep properties of the rocks, especially soft rocks. Therefore, the Drucker–Prager yield criterion [25] is adopted in this study. This yield criterion can effectively compensate for the shortcomings of the first two yielding criteria and is expressed as(30)F=J2−αI1−k,where J2 is the second deviatoric stress invariant, I1 is the first invariant of the stress tensor, and αandk are the material parameters expressed as follows:(31)α=sinφ33+sin2φk=3ccosφ3+sin2φ.Here,φ is the internal friction angle of the rock material and c is the sticky cohesion of the rock material.The following equations can be obtained in the conventional triaxial step cycle loading and unloading compression creep test:(32)σ1>σ2=σ3σm=σ1+σ2+σ33=σ1+2σ33J2=σ1−σ33I1=σ1+σ2+σ3=σ1+2σ3S11=2σ1−σ33∂F∂σ11=33−α.According to equations (26)–(32), the three-dimensional creep equations of the nonlinear damage creep model are expressed as follows:(33)ε=σ1+2σ39Ke+σ1−σ33Ge+σ1−σ33Gnk1−exp−Gnkηnkλtλσ<σsσ1+2σ39Ke+σ1−σ33Ge+σ1−σ33Gnk1−exp−Gnkηnkλtλ+3−3α3tf1+n1−1−Dn2ηvpnσ1−σ3−σsσs≤σ.Assuming the applied stress is unloaded att0, the unloading equations of the nonlinear damage creep model are expressed as(34)ε=σ1−σ33Gnk1−exp−Gnkηnkλt0λexpGnkηnkt01−λt0−t.


## 3. Validation of the Nonlinear Damage Creep Model

### 3.1. Parameters Determination

#### 3.1.1. Determination ofGe and Ke

Ge and Ke can be determined by the elastic modulus Ee and Poisson's ratio μ according to the following equation:(35)Ge=Ee21+μKe=Ee31−2μ.


#### 3.1.2. Determination ofGnk,ηnk, and λ

Whenσ<σs and t⟶∞, equation (33) can be converted into the following:(36)εt⟶∞=σ1+2σ39Ke+σ1−σ33Ge+σ1−σ33Gnk.Subtracting equation (36) from equation (33) and taking the logarithmic operation on both sides of equation (36), we can obtain the following:(37)lnεt⟶∞−ε=−Gnkηnkλtλ+lnσ1−σ33Gnk.Supposeft=lnεt⟶∞−εt,a0=−Gnk/ηnkλ,b0=λ,andc0=lnσ1−σ3/3Gnk, then equation (37) becomes a nonlinear equation about a0,b0, and c0 as follows:(38)ft=a0tb0+c0.The fitting parametersa0,b0,andc0 can be obtained by nonlinear fitting based on the experimental data during the stable creep stage. Further, Gnk,ηnk, and λ can be derived from the fitting parameters a0,b0, and c0 as follows:(39)λ=b0,Gnk=σ1−σ33ec0,ηnk=−σ1−σ33a0b0ec0.


#### 3.1.3. Determination ofηvp and n

Numerous studies have been conducted to determine the creep model parameters; these studies can be divided into two types. The first is a graphical method, which determines the creep parameters based on the relation between the creep test curve and physical significance of the creep parameters. The second is an optimization analysis algorithm method, such as the regression analysis method and the least squares method (LSM). In this study, the widely used Levenberg–Marquardt nonlinear least squares method [26] is adopted to determine ηvp and n because of its simplicity and applicability.


### 3.2. Model Verification

To verify the rationality and feasibility of the nonlinear damage creep model, the proposed model was evaluated by multilevel loading and unloading cycles compression creep test data for lherzolite and limestone [27, 28]. The basic mechanical parameters of the rock samples are shown in Table 1.Table 1
Basic parameters of the rock specimens.
Rock specimensConfining pressure (MPa)Compressive strength (MPa)Elastic modulus (GPa)Poisson's ratioSticky cohesion (MPa)Internal friction angle (°)Lherzolite378.6319.070.2278.6252.64Limestone459.207.790.2518.1235.10As shown in Table2, the loading level was six and the stress at each loading level was 25.2%, 42.0%, 50.4%, 54.6%, 58.84%, and 75.6% of the peak strength, respectively. During the unloading process, the confining pressure was kept constant, the axial stress was unloaded to the same level as the confining pressure, and the deviatoric stress was maintained at zero. Except for the last stage, the creep duration of each stage was 90 h, and the unloading rest time was 20–30 h. The rock specimen exhibited creep failure until the last load level.Table 2
Test scheme for triaxial creep.
Rock specimenConfining pressure (MPa)Deviatoric stress level (MPa)First levelSecond levelThird levelFourth levelFifth levelSixth levelLherzolite313.322.126.628.831.039.9Limestone429.633.838.042.246.450.6The model parameters of lherzolite and limestone can be obtained based on the experimental data and parameter identification method (as shown in Tables3 and 4). By substituting the model parameters into equations (33) and (34), the theoretical curves of the nonlinear damage creep model are obtained, as shown in Figures 5 and 6.Table 3
Model parameters under triaxial creep for lherzolite.
Stress levels (MPa)Ke (GPa)Ge (GPa)Gnk (GPa)ηnk (GPa·hλ)λαtf (h)nηvp (GPa·h)13.311.647.774.6819.690.750.24———22.13.1713.811.09———26.62.9117.871.35———28.82.6822.121.57———31.02.176.120.67———39.92.899.940.736554.8220.45Table 4
Model parameters under triaxial creep for limestone.
Stress levels (MPa)Ke (GPa)Ge (GPa)Gnk (GPa)ηnk (GPa·hλ)λαtf (h)nηvp (GPa·h)29.65.213.1112.5532.270.910.18———33.811.04118.891.60———38.09.6694.041.78———42.26.0868.581.87———46.44.1114.940.49———50.64.9424.380.786763.3525.18Figure 5
Comparisons between the nonlinear damage creep model and experimental curves at various deviatoric stresses for lherzolite specimen: (a) first stress (13.3 MPa) and second stress (22.1 MPa); (b) third stress (26.6 MPa) and fourth stress (28.8 MPa); (c) fifth stress (31.0 MPa) and sixth stress (39.9 MPa).
(a)(b)(c)Figure 6
Comparisons between the nonlinear damage creep model and experimental curves at various deviatoric stresses for limestone specimen.A comparison between the experimental and theoretical curves of the nonlinear damage creep model (Figures5 and 6) indicates that the experimental values agree well with the theoretical values. Particularly, the theoretical value is very close to the experimental value in the accelerating creep stage, which shows that the proposed model can not only reflect the loading and unloading creep characteristics of lherzolite and limestone at different stress levels but also characterize the nonlinear characteristics of the attenuation, steady-state, and accelerating creep stages. Therefore, the model established in this study can be used to determine the nonlinear creep characteristics of rocks during loading and unloading, which overcomes the shortcomings of the classical linear element combination model.


## 3.1. Parameters Determination

### 3.1.1. Determination ofGe and Ke

Ge and Ke can be determined by the elastic modulus Ee and Poisson's ratio μ according to the following equation:(35)Ge=Ee21+μKe=Ee31−2μ.


### 3.1.2. Determination ofGnk,ηnk, and λ

Whenσ<σs and t⟶∞, equation (33) can be converted into the following:(36)εt⟶∞=σ1+2σ39Ke+σ1−σ33Ge+σ1−σ33Gnk.Subtracting equation (36) from equation (33) and taking the logarithmic operation on both sides of equation (36), we can obtain the following:(37)lnεt⟶∞−ε=−Gnkηnkλtλ+lnσ1−σ33Gnk.Supposeft=lnεt⟶∞−εt,a0=−Gnk/ηnkλ,b0=λ,andc0=lnσ1−σ3/3Gnk, then equation (37) becomes a nonlinear equation about a0,b0, and c0 as follows:(38)ft=a0tb0+c0.The fitting parametersa0,b0,andc0 can be obtained by nonlinear fitting based on the experimental data during the stable creep stage. Further, Gnk,ηnk, and λ can be derived from the fitting parameters a0,b0, and c0 as follows:(39)λ=b0,Gnk=σ1−σ33ec0,ηnk=−σ1−σ33a0b0ec0.


### 3.1.3. Determination ofηvp and n

Numerous studies have been conducted to determine the creep model parameters; these studies can be divided into two types. The first is a graphical method, which determines the creep parameters based on the relation between the creep test curve and physical significance of the creep parameters. The second is an optimization analysis algorithm method, such as the regression analysis method and the least squares method (LSM). In this study, the widely used Levenberg–Marquardt nonlinear least squares method [26] is adopted to determine ηvp and n because of its simplicity and applicability.


## 3.1.1. Determination ofGe and Ke

Ge and Ke can be determined by the elastic modulus Ee and Poisson's ratio μ according to the following equation:(35)Ge=Ee21+μKe=Ee31−2μ.


## 3.1.2. Determination ofGnk,ηnk, and λ

Whenσ<σs and t⟶∞, equation (33) can be converted into the following:(36)εt⟶∞=σ1+2σ39Ke+σ1−σ33Ge+σ1−σ33Gnk.Subtracting equation (36) from equation (33) and taking the logarithmic operation on both sides of equation (36), we can obtain the following:(37)lnεt⟶∞−ε=−Gnkηnkλtλ+lnσ1−σ33Gnk.Supposeft=lnεt⟶∞−εt,a0=−Gnk/ηnkλ,b0=λ,andc0=lnσ1−σ3/3Gnk, then equation (37) becomes a nonlinear equation about a0,b0, and c0 as follows:(38)ft=a0tb0+c0.The fitting parametersa0,b0,andc0 can be obtained by nonlinear fitting based on the experimental data during the stable creep stage. Further, Gnk,ηnk, and λ can be derived from the fitting parameters a0,b0, and c0 as follows:(39)λ=b0,Gnk=σ1−σ33ec0,ηnk=−σ1−σ33a0b0ec0.


## 3.1.3. Determination ofηvp and n

Numerous studies have been conducted to determine the creep model parameters; these studies can be divided into two types. The first is a graphical method, which determines the creep parameters based on the relation between the creep test curve and physical significance of the creep parameters. The second is an optimization analysis algorithm method, such as the regression analysis method and the least squares method (LSM). In this study, the widely used Levenberg–Marquardt nonlinear least squares method [26] is adopted to determine ηvp and n because of its simplicity and applicability.


## 3.2. Model Verification

To verify the rationality and feasibility of the nonlinear damage creep model, the proposed model was evaluated by multilevel loading and unloading cycles compression creep test data for lherzolite and limestone [27, 28]. The basic mechanical parameters of the rock samples are shown in Table 1.Table 1
Basic parameters of the rock specimens.
Rock specimensConfining pressure (MPa)Compressive strength (MPa)Elastic modulus (GPa)Poisson's ratioSticky cohesion (MPa)Internal friction angle (°)Lherzolite378.6319.070.2278.6252.64Limestone459.207.790.2518.1235.10As shown in Table2, the loading level was six and the stress at each loading level was 25.2%, 42.0%, 50.4%, 54.6%, 58.84%, and 75.6% of the peak strength, respectively. During the unloading process, the confining pressure was kept constant, the axial stress was unloaded to the same level as the confining pressure, and the deviatoric stress was maintained at zero. Except for the last stage, the creep duration of each stage was 90 h, and the unloading rest time was 20–30 h. The rock specimen exhibited creep failure until the last load level.Table 2
Test scheme for triaxial creep.
Rock specimenConfining pressure (MPa)Deviatoric stress level (MPa)First levelSecond levelThird levelFourth levelFifth levelSixth levelLherzolite313.322.126.628.831.039.9Limestone429.633.838.042.246.450.6The model parameters of lherzolite and limestone can be obtained based on the experimental data and parameter identification method (as shown in Tables3 and 4). By substituting the model parameters into equations (33) and (34), the theoretical curves of the nonlinear damage creep model are obtained, as shown in Figures 5 and 6.Table 3
Model parameters under triaxial creep for lherzolite.
Stress levels (MPa)Ke (GPa)Ge (GPa)Gnk (GPa)ηnk (GPa·hλ)λαtf (h)nηvp (GPa·h)13.311.647.774.6819.690.750.24———22.13.1713.811.09———26.62.9117.871.35———28.82.6822.121.57———31.02.176.120.67———39.92.899.940.736554.8220.45Table 4
Model parameters under triaxial creep for limestone.
Stress levels (MPa)Ke (GPa)Ge (GPa)Gnk (GPa)ηnk (GPa·hλ)λαtf (h)nηvp (GPa·h)29.65.213.1112.5532.270.910.18———33.811.04118.891.60———38.09.6694.041.78———42.26.0868.581.87———46.44.1114.940.49———50.64.9424.380.786763.3525.18Figure 5
Comparisons between the nonlinear damage creep model and experimental curves at various deviatoric stresses for lherzolite specimen: (a) first stress (13.3 MPa) and second stress (22.1 MPa); (b) third stress (26.6 MPa) and fourth stress (28.8 MPa); (c) fifth stress (31.0 MPa) and sixth stress (39.9 MPa).
(a)(b)(c)Figure 6
Comparisons between the nonlinear damage creep model and experimental curves at various deviatoric stresses for limestone specimen.A comparison between the experimental and theoretical curves of the nonlinear damage creep model (Figures5 and 6) indicates that the experimental values agree well with the theoretical values. Particularly, the theoretical value is very close to the experimental value in the accelerating creep stage, which shows that the proposed model can not only reflect the loading and unloading creep characteristics of lherzolite and limestone at different stress levels but also characterize the nonlinear characteristics of the attenuation, steady-state, and accelerating creep stages. Therefore, the model established in this study can be used to determine the nonlinear creep characteristics of rocks during loading and unloading, which overcomes the shortcomings of the classical linear element combination model.


## 4. Conclusions

(1)
In this study, a novel nonlinear damage creep model for rocks incorporating instantaneous elasticity, viscoelasticity, and viscoplasticity was proposed based on the elastic body, nonlinear Kelvin model, and damage viscoplastic model.(2)
The creep equations and unloading equations for rocks under one-dimensional and three-dimensional stress conditions were derived by introducing the creep mechanics theory, respectively. The parameters of the nonlinear damage creep model were determined using the nonlinear least squares method.(3)
The proposed nonlinear damage creep model was used to fit lherzolite and limestone creep test data; the results showed that the model accurately determines the loading and unloading creep properties of lherzolite and limestone under different stress levels and also favorably represents the nonlinear characteristics of accelerating creep, thus validating our model.



---

*Source: 5512972-2021-04-12.xml*