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

To determine the nonlinear creep characteristics of rocks under cyclic loading and unloading conditions, a nonlinear Kelvin model and damage viscoplastic model are proposed. ,e models are connected in series with a linear elastic body to establish a nonlinear damage creep model. ,e 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. ,e 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. ,e 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.


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]. e empirical creep model is a mathematical expression established by fitting the test results. e 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. is 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]. e component models are constructed by connecting the components with basic functions (including elastic, plastic, and viscosity components) in series and parallel. ese 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]. e 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. erefore, 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. is 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.

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). e total strain can then be expressed as follows: where ε is the total strain of rock under uniaxial compression creep state, ε m and ε 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. e 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. erefore, the novel nonlinear damage creep model established in this study does not consider the effect of instantaneous plastic strain [14]. e relationship between the stress and instantaneous elastic strain can be described by the linear elastic body.

Nonlinear Kelvin Model.
e transient and unloading creep curves have obvious nonlinear characteristics, which are difficult to accurately determine using the traditional Kelvin model [15].
is is because the parameters of the Kelvin model are constant and do not change with creep time.
erefore, 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. e mechanical model is shown in Figure 2. e differential constitutive equation of the nonlinear Kelvin model can be expressed as where σ is the applied stress; E nk and η 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: Assuming the applied stress is unloaded at t 0 , equation (2) can be expressed as Equation (5) can be obtained by solving the differential equation (4) as follows: where A is the integral constant. When t � t 0 , ε ve � (σ/E nk )[1 − exp(− (E nk /η nk λ)t λ 0 )], and this boundary condition is substituted into equation (5) to obtain the integral constant as follows:

Advances in Materials Science and Engineering
By substituting equation (6) into equation (5), the creep strain ε ul during the unloading stage can be expressed as

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. e 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]. e law governing the change between the stress and viscoplastic strain in the accelerated creep stage can be determined by the damage viscoplastic model. e viscosity coefficient, η vp , in the nondamaged state is replaced by the effective viscosity coefficient, η vp (1 − D), to determine the deterioration process of the rock. e mechanical model is shown in Figure 3.
Based on Lemaitre's strain equivalence principle [18,19], the stress-strain relation of the damage viscous body can be expressed as follows: where σ is the effective stress and D is the damage variable. e damage evolution equation proposed by Kachanov [20] is widely used to determine the damage behavior of rocks during the accelerated creep stage.
where C and n represent the rock material constants. e creep damage critical failure time can be derived using equation (9) as follows: where t f represents the critical failure time. e relationship between the damage variable D and time can be obtained using equations (9) and (10) as follows: We can obtain the constitutive equation of the damaged viscous body by substituting equation (11) into equation (8) as follows: By integrating equation (12), we can obtain the following equation: where B is another integral constant. When t � 0, ε � 0, and this boundary condition is substituted into equation (13) to obtain the integral constant as follows: By substituting equation (14) into equation (13), we obtain By substituting equation (10) into equation (15), we obtain In previous research, the creep equation of the damage viscous body was obtained using the Kachanov creep damage rate as follows: 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: 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 .

Advances in Materials Science and Engineering
where σ s is the long-term strength, which can be obtained by the creep test.

Nonlinear Damage Creep
Model. e 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 Figure 4.
Based on the stress-strain relationship of the model in the series-parallel connection, we can obtain the following equations: 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. e constitutive equations for each part are expressed as where E e and E nk are the elastic modulus, η nk and η vp are the viscosity coefficients, and _ ε ve and _ ε vp are the strain rates. When it is one-dimensional, the constitutive equation is as follows: (1) when σ < σ s , (2) when σ s ≤ σ, e axial creep equations of the nonlinear damage rheological model can be obtained based on the superposition principle.

ree-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 For an instantaneous elastic body, the strain under threedimensional stress conditions is obtained according to the general Hooke's law [21] as where G e and K e are the shear modulus and bulk modulus, respectively, S ij and σ 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: 4 Advances in Materials Science and Engineering where G nk is the shear modulus of the nonlinear Kelvin model. e three-dimensional constitutive relationship of the damage viscoplastic model is where 〈•〉 is the switch function, F is the rock yield function, F 0 is the initial value of the rock yield function, which is typically 1 for simplicity [22], and Φ(F/F 0 ) 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/F 0 ) � F/F 0 . Q is the plastic potential function. Moreover, for convenience of calculation, the associated flow rule is used; therefore, F � Q [24]. e 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. erefore, the Drucker-Prager yield criterion [25] is adopted in this study. is yield criterion can effectively compensate for the shortcomings of the first two yielding criteria and is expressed as where J 2 is the second deviatoric stress invariant, I 1 is the first invariant of the stress tensor, and α and k are the material parameters expressed as follows: Here, φ is the internal friction angle of the rock material and c is the sticky cohesion of the rock material. e following equations can be obtained in the conventional triaxial step cycle loading and unloading compression creep test: According to equations (26)-(32), the three-dimensional creep equations of the nonlinear damage creep model are expressed as follows: Assuming the applied stress is unloaded at t 0 , the unloading equations of the nonlinear damage creep model are expressed as  Figure 4: Nonlinear damage creep model.

Advances in Materials Science and Engineering
3.1.2. Determination of G nk , η nk , and λ. When σ < σ s and t ⟶ ∞, equation (33) can be converted into the following: Subtracting equation (36) from equation (33) and taking the logarithmic operation on both sides of equation (36), we can obtain the following: , a 0 � − (G nk /η nk λ) ,b 0 � λ, and c 0 � ln(σ 1 − σ 3 /3G nk ), then equation (37) becomes a nonlinear equation about a 0 , b 0 , and c 0 as follows: e fitting parameters a 0 , b 0 , and c 0 can be obtained by nonlinear fitting based on the experimental data during the stable creep stage. Further, G nk , η nk , and λ can be derived from the fitting parameters a 0 , b 0 , and c 0 as follows:

Determination of η vp and n.
Numerous studies have been conducted to determine the creep model parameters; these studies can be divided into two types. e 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. e 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.

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]. e basic mechanical parameters of the rock samples are shown in Table 1.
As shown in Table 2, 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. e rock specimen exhibited creep failure until the last load level. e model parameters of lherzolite and limestone can be obtained based on the experimental data and parameter identification method (as shown in Tables 3 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.
A comparison between the experimental and theoretical curves of the nonlinear damage creep model (Figures 5 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. erefore, 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.       (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) e creep equations and unloading equations for rocks under one-dimensional and three-dimensional stress conditions were derived by introducing the creep mechanics theory, respectively. e parameters of the nonlinear damage creep model were determined using the nonlinear least squares method. (3) e 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.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding the publication of this paper.