#### 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.

#### 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:where is the total strain of rock under uniaxial compression creep state, are the instantaneous strain and creep strain, respectively, and are the instantaneous elastic strain, instantaneous plastic strain, viscoelastic strain, viscoplastic strain, and residual plastic strain, respectively.

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.

The differential constitutive equation of the nonlinear Kelvin model can be expressed aswhere is the applied stress; 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 , 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 , , and this boundary condition is substituted into equation (5) to obtain the integral constant as follows:

By substituting equation (6) into equation (5), the creep strain during the unloading stage can be expressed as

##### 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, , in the nondamaged state is replaced by the effective viscosity coefficient, , to determine the deterioration process of the rock. The 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.

The 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.

The creep damage critical failure time can be derived using equation (9) as follows:where represents the critical failure time.

The 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 , , 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 .where 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 Figure 4.

Based on the stress-strain relationship of the model in the series-parallel connection, we can obtain the following equations:where are the total stress and strain, respectively, and represent the three-part stress and strain of the nonlinear damage creep model, respectively.

The constitutive equations for each part are expressed aswhere are the elastic modulus, are the viscosity coefficients, and are the strain rates.

When it is one-dimensional, the constitutive equation is as follows:(1)when ,(2)when ,

The axial creep equations of the nonlinear damage rheological model can be obtained based on the superposition principle.(1)If , then(2)If , then

##### 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, , of the nonlinear damage creep model under three-dimensional stress conditions can be expressed as

For an instantaneous elastic body, the strain under three-dimensional stress conditions is obtained according to the general Hooke's law [21] aswhere are the shear modulus and bulk modulus, respectively, are the stress tensor and spherical stress tensor, respectively, and 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:where is the shear modulus of the nonlinear Kelvin model.

The three-dimensional constitutive relationship of the damage viscoplastic model iswhere is the switch function, is the rock yield function, is the initial value of the rock yield function, which is typically 1 for simplicity [22], and 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]. . is the plastic potential function. Moreover, for convenience of calculation, the associated flow rule is used; therefore, [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 aswhere is the second deviatoric stress invariant, is the first invariant of the stress tensor, and are the material parameters expressed as follows:

Here, is the internal friction angle of the rock material and 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:

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 , the unloading equations of the nonlinear damage creep model are expressed as

#### 3. Validation of the Nonlinear Damage Creep Model

##### 3.1. Parameters Determination

###### 3.1.1. Determination of and

and can be determined by the elastic modulus and Poisson's ratio according to the following equation:

###### 3.1.2. Determination of and

When and , 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:

Suppose , then equation (37) becomes a nonlinear equation about , and as follows:

The fitting parameters can be obtained by nonlinear fitting based on the experimental data during the stable creep stage. Further, , and can be derived from the fitting parameters , and as follows:

###### 3.1.3. Determination of and

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 and 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.

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. The rock specimen exhibited creep failure until the last load level.

The 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)**

**(b)**

**(c)**

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. 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.#### Data Availability

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

#### Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was supported by the National Natural Science Foundation of China (51374010 and 51474004).