Nonlinear Elasto-Visco-Plastic Creep Behavior and New Creep Damage Model of Dolomitic Limestone Subjected to Cyclic Incremental Loading and Unloading

: For many rock engineering projects, the stress of surrounding rocks is constantly increasing and decreasing during excavating progress and the long-term operation stage. Herein, the triaxial creep behavior of dolomitic limestone subjected to cyclic incremental loading and unloading was probed using an advanced rock mechanics testing system (i.e., MTS815.04). Then, the instantaneous elastic strain, instantaneous plastic strain, visco-elastic strain, and visco-plastic strain components were separated from the total strain curve, and evolutions of these different types of strain with deviatoric stress increment were analyzed. Furthermore, a damage variable considering the proportion of irrecoverable plastic strain to the total strain was introduced, and a new nonlinear multi-element creep model was established by connecting the newly proposed damage viscous body in series with the Hookean substance, St. Venant body, and Kelvin element. The parameters of this new model were analyzed. The ﬁndings are listed as follows: (1) When the deviatoric stress is not more than 75% of the compressive strength, only instantaneous deformation, transient creep, and steady-state creep deformation occur, rock deformation is mainly characterized by the instantaneous strain, whereas the irrecoverable instantaneous plastic strain accounts for 38.02–60.27% of the total instantaneous strain; (2) Greater deviatoric stress corresponds to more obvious creep deformation. The visco-elastic strain increases linearly with the increase of deviatoric stress, especially the irrecoverable visco-plastic strain increases exponentially with deviatoric stress increment, and ﬁnally leads to accelerated creep and delayed failure of the sample; (3) Based on the experimental data, the proposed nonlinear creep model is veriﬁed to describe the full creep stage perfectly, particularly the tertiary creep stage. These results could deepen our understanding of the elasto-visco-plastic deformation behavior of dolomitic limestone and have theoretical and practical signiﬁcance for the safe excavation and long-term stability of underground rock engineering.


Introduction
The time-dependent deformation characteristics of rock, especially the creep behaviors, are related to the long-term stability of underground rock engineering, such as tunnels and roadways [1][2][3][4][5][6][7][8][9][10][11]. Many scholars have carried out extensive research on creep mechanical behaviors of different types of rocks. Maranini and Brignoli [12], Yang et al. [11], Fujii et al. [13], and Wang et al. [14] have obtained a lot of research achievements in uniaxial and triaxial compression creep characteristics of soft rocks such as limestone, salt, clasticrock, weathered sandstone, and coal rock. The creep mechanical properties of hard rocks, such as greenschist [15], diabase [16], marble [17], and red sandstone [18] under stepwise Sustainability 2021, 13, 12376 2 of 15 loading, have been systematically studied. Furthermore, the constitutive model theory has been developed and applied to describe rock creep behavior. There are two main expression types for rheological models: One is the empirical model, the curtain function, exponential function, and logarithmic function are commonly used to fit the creep test data. Generally, the fitting effect is perfect, but the physical meaning of the model is not clear [14]. The other is the component combined model. The classical multi-component models include the Maxwell model, Kelvin model, Bingham model, and Burgers model [9]. However, since these models are composed of linear elements, they cannot describe the nonlinear accelerated rheological behavior [9,14]. Thus, more and more nonlinear rheological models have been established to describe the whole creep process reasonably. There are mainly two modeling methods: one is to propose a new nonlinear rheological element, and then connect it with the original linear model to build a nonlinear creep model that can describe the three creep stages [16,17]; The other is to establish a nonlinear creep damage model based on fracture or damage mechanics, which could reasonably describe the damage evolution during three creep stages of rocks [18,19].
However, the triaxial compression creep tests are still limited due to the limitation of test apparatus and strict requirements for creep test conditions; Previous studies have paid much attention to the creep behavior of rock under a constant load or multi-stage load, with little being paid to the time-dependent deformation behavior of rock during unloading. Previous creep experiments seldom have been carried out to separate the instantaneous elastic, instantaneous plastic strain, visco-elastic strain, and visco-plastic strain components from creep strain curves [20]. Furthermore, the previous constitutive models have not sufficiently considered the difference between visco-elastic and visco-plastic strains. For a lot of rock engineering, the stress of surrounding rocks is constantly changing and adjusting between increasing and decreasing during the excavating progress and long-term operation stage [20][21][22]. Therefore, the time-dependent deformation behavior of surrounding rocks subjected to cyclic loading and unloading is a scientific problem that needs to be focused on, which has important theoretical and practical significance for the safe construction and long-term operation of rock engineering.
In this work, the chosen dolomitic limestone specimens were drilled from Yangzong Tunnel in Yunnan Province, southwest China. The triaxial compression creep experiments were performed to capture the visco-elastic-plastic deformation behaviors of the dolomitic limestone under cyclic incremental loading and unloading. Creep data was analyzed, and the total creep strain was separated into the instantaneous elastic strain, instantaneous plastic strain, visco-elastic strain, and visco-plastic strain components. Then, a new nonlinear visco-elasto-plastic creep model was proposed to describe the full creep stages, especially the tertiary stage. Finally, the creep parameters were identified by employing the Levenberg-Marquardt optimization algorithm and the accuracy of the creep model was verified based on the experimental data.

Experimental Method
The creep experiment was carried out using a servo-controlled high rigidity rock mechanics testing system (MTS815.04) at Shaoxing University in China. The specific parameters of this apparatus: maximum axial pressure is 4600 kN; maximum confining pressure is 140 MPa; adjustable temperature is in the range of 0-200 • C The duration, axial pressure, confining pressure, axial and lateral strain can be displayed on the computer monitor during the testing progress. The size of cylindrical samples used in both the conventional compression experiment and the creep experiment are Φ50 mm × 100 mm. The test apparatus and the deformation monitoring sensor system are shown in Figure 1. A conventional triaxial compression test was carried out first to reveal the short-term mechanical properties of the studied rock and provide a basis for determining the stress level in the creep test. According to in-situ stress measurement, the horizontal stress is about 9 MPa, which is designed as the confining stress in conventional test and creep tests. An axial displacement rate of 0.001 mm/s and a confining stress loading rate of 0.1 MPa/s were used for conventional compression experiments. Based on testing results, the short-term mechanical parameters could be calculated, as shown in Table 1. The confining pressure was applied first at a rate of 0.1 MPa/s to achieve 9 MPa, which then was maintained unchanged throughout creep testing until the specimen failed. Six levels of axial stress were set in the creep test, as listed in Table 1, and the first axial loading stress was 35% of the axial compressive strength, i.e., the first stress level was 27 MPa. The loading and unloading process is shown in Figure 2, an axial loading rate of 0.5 kN/s was adopted to reach the target stress and then remained constant for 48 h, after which, the deviatoric stress (i.e., the difference between axial stress and confining stress) was unloaded to 0 at a rate of 0.5 kN/s and then maintained for 20-30 h until the visco-elastic strain was recovered. After this unloading period, the loading and unloading of the next axial stress level were performed, until the specimen was damaged, the experiment was completed. The axial and radial strain data of the specimen were acquired in real-time by the testing machine and stored in the computer during the whole experiment.

Data Processing
Three samples were tested in the creep experiment, and the results show little scattering due to the natural heterogeneities and discreteness of the rock samples. In this paper, only the most representative test results are selected and analyzed. Figure 3 shows the total axial and radial creep curves of dolomitic limestone subjected to cyclic incremental loading and unloading. Taking the axial strain-time curve (as shown in Figure 4) under loading and unloading of a certain deviatoric stress level (i.e., cycle i) for example, the elasto-visco-plastic strain separation is analyzed as follows.  The total strain could be separated into the instantaneous elastic, instantaneous plastic strain, visco-elastic strain, and visco-plastic strain components.
where ε ve during the unloading stage, according to the principles of elastic and plastic mechanics [23]. However, the plastic strain can only be obtained by considering the history of loading [24], i.e., ε (i) mp and ε (i) vp can be expressed as: where ∆ε It should be pointed out that the evolution law of radial creep strain is similar to that of axial creep strain, the elasto-visco-plastic strain characteristics of radial creep behavior would not be analyzed here.

Elasto-Visco-Plastic Strain Analysis
The strain curves (in Figure 3), during different deviatoric stress levels, indicate that instantaneous deformation occurred firstly, after which creep deformation occurred. When the loading deviatoric stress was not more than 51 MPa (i.e., the fifth deviatoric stress level), only transient creep stage and steady-state creep stage occurred. However, the rock underwent accelerated creep deformation as the deviatoric stress reached 56 MPa, which accounts for 82.35% of the conventional compressive strength.
According to the strain separation method shown in Figure 4, the axial elasto-viscoplastic strain data of the specimen subjected to cyclic incremental loading and unloading could be obtained, as listed in Table 2. Additionally, the variation curves of instantaneous elastic strain, instantaneous plastic strain, visco-elastic strain, and visco-plastic strain with deviatoric stress are illustrated in Figure 5.
where mp ε Δ and vp ε Δ are the instantaneous plastic strain increment and visco-plastic strain increment caused by deviatoric stress increment It should be pointed out that the evolution law of radial creep strain is similar to that of axial creep strain, the elasto-visco-plastic strain characteristics of radial creep behavior would not be analyzed here.

Elasto-Visco-Plastic Strain Analysis
The strain curves (in Figure 3), during different deviatoric stress levels, indicate that instantaneous deformation occurred firstly, after which creep deformation occurred. When the loading deviatoric stress was not more than 51 MPa (i.e., the fifth deviatoric stress level), only transient creep stage and steady-state creep stage occurred. However, the rock underwent accelerated creep deformation as the deviatoric stress reached 56 MPa, which accounts for 82.35% of the conventional compressive strength.
According to the strain separation method shown in Figure 4, the axial elasto-viscoplastic strain data of the specimen subjected to cyclic incremental loading and unloading could be obtained, as listed in Table 2. Additionally, the variation curves of instantaneous elastic strain, instantaneous plastic strain, visco-elastic strain, and visco-plastic strain with deviatoric stress are illustrated in Figure 5.  where mp ε Δ and vp ε Δ are the instantaneous plastic strain increment and visco-plasti strain increment caused by deviatoric stress increment

Confining Pressure
It should be pointed out that the evolution law of radial creep strain is similar to tha of axial creep strain, the elasto-visco-plastic strain characteristics of radial creep behavio would not be analyzed here.

Elasto-Visco-Plastic Strain Analysis
The strain curves (in Figure 3), during different deviatoric stress levels, indicate tha instantaneous deformation occurred firstly, after which creep deformation occurred When the loading deviatoric stress was not more than 51 MPa (i.e., the fifth deviatori stress level), only transient creep stage and steady-state creep stage occurred. Howeve the rock underwent accelerated creep deformation as the deviatoric stress reached 5 MPa, which accounts for 82.35% of the conventional compressive strength.
According to the strain separation method shown in Figure 4, the axial elasto-visco plastic strain data of the specimen subjected to cyclic incremental loading and unloadin could be obtained, as listed in Table 2. Additionally, the variation curves of instantaneou elastic strain, instantaneous plastic strain, visco-elastic strain, and visco-plastic strain wit deviatoric stress are illustrated in Figure 5.  It is observed from Figure 5a,b that both instantaneous elastic and instantaneous plastic strain increase linearly with deviatoric stress. The irrecoverable instantaneous plastic strain accounts for 37.97-60.27% of instantaneous strain (using the axial elasto-viscoplastic strain data listed in Table 2), yet it would often belong to recoverable instantaneous elastic deformation in the creep tests without unloading process [14,16], which may cause errors in engineering design. Besides, the creep strain, which consists of visco-elastic strain and visco-plastic strain, increases with deviatoric stress increment. Specifically, as the deviatoric stress increases from 18 MPa to 56 MPa, the proportion of creep strain in total strain increases from 5.06% to 33.56%, which means creep deformation is more and more obvious as the deviatoric stress increases. In addition, it should be noted that the viscoplastic strain increases nonlinearly and exponentially with the increase of deviatoric stress, as shown in Figure 5d. Moreover, the proportion of irrecoverable visco-plastic strain in total creep strain increases from 35.00% to 63.89% with the deviatoric stress increasing from 18 MPa to 51 MPa. It suggests that the recoverable visco-elastic strain accounts for the main part of the total creep strain when the deviatoric stress is not more than 51 MPa. However, the creep deformation is mainly characterized by visco-plastic strain as the deviatoric stress is greater than 51 MPa, especially during the accelerated creep stage, the total creep deformation could be considered as the development of irrecoverable visco-plastic strain, which leads to delayed failure of the specimen. As illustrated in Figure 6, the specimen shows shear failure, an inclined failure surface runs through the whole specimen, which is similar to the creep failure mode of hard rock found by Zhao et al. [20].

Development of Nonlinear Creep Damage Model
In Figure 7, a new multi-element rheological model is developed to describe the nonlinear visco-elasto-plastic creep behavior of dolomitic limestone under cyclic incremental loading and unloading. This multi-element model consists of the Hookean substance, St. Venant substance, Kelvin body, and a newly proposed damage viscous body, which are used to exhibit instantaneous elastic strain, instantaneous plastic strain, visco-elastic strain, and nonlinear visco-plastic strain, respectively. Their constitutive relations and creep equations are derived as follows.

Instantaneous Strain Model
According to the constitutive relation of Hooke body [25], axial instantaneous elastic strain ε me is given by where E 1 is the elastic modulus of the single Hooke body 1. Similarly, the axial instantaneous plastic strain could be expressed by [26] where P is the plastic deformation modulus of the St. Venant substance, and I(σ 1 − σ 3 − σ s ) is a unit jump function, namely,

Visco-Elastic Strain Model
Based on the constitutive relationship of Kevin body [25,27], the axial visco-elastic strain of the specimen during the loading stage could be presented as where E 2 and η 1 are elastic modulus and viscous coefficient of Kelvin substance, respectively, and t is creep time.
At the time point t = t 1 , the deviator stress σ 1 − σ 3 is unloaded to zero, the unloading creep equation of visco-elastic strain ε ve is given by Equation (8) [28].

Visco-Plastic Strain Model
In terms of damage mechanics, rock failure is the joint action result of external load and evolution of internal defects of rock materials, which is a gradual creep accumulation process [9,18,29]. Under the long-term influence of large deviatoric stress, the damage of rock materials will continue to accumulate, resulting in fissures increase in the amount and crack expansion and penetration. In this process, the damage degree of rock will develop rapidly, the viscous coefficient will continue to decrease. Macroscopically, it shows the accumulation and nonlinear increase of unrecoverable visco-plastic strain (as shown in Figures 3 and 5d). Finally, the rock deformation increases sharply, and the specimen is destroyed. Therefore, similar to the elastic-plastic theory, a new damage variable could be defined by considering the development of the plastic strain of the rock material, i.e., where ε e and ε p are the elastic and plastic strain of rock material, respectively. The former is a summation of instantaneous elastic strain and visco-elastic strain, the latter consists of instantaneous plastic strain and visco-plastic strain. At the initial loading time point t = 0, instantaneous plastic deformation occurs, and the damage variable, D = 0, indicates that the damage variable can also reflect the initial damage, which is an essential feature of rock due to the existence of joints and micro-fractures. Then, based on the constitutive relationship of dashpot element [27], the constitutive equation of this newly proposed damaged viscous body could be given as where η 2 is the initial viscosity coefficient of visco-plastic element, . ε vp is the creep rate of visco-plastic strain.
Then, the following Equation (10) could be obtained. .
where q = (σ 1 − σ 3 ). Besides, the visco-elastic strain could be omitted due to the fact that the visco-elastic strain of rock is more than one order of magnitude smaller than the instantaneous elastic strain during the triaxial compression creep process (as proven in Table 2), then, Equation (11) can be simplified to Equation (12).
Equation (12) is a first-order linear non-homogeneous differential equation concerning time. The creep equation of visco-plastic strain can be obtained by solving this differential equation, as shown in Equation (13). During creep testing, when the loading stress is constant, the instantaneous elastic strain and instantaneous plastic strain of the sample are constant, so Equation (13) can be expressed as Equation (14).
where A is a variable related to instantaneous strain, B is a variable related to instantaneous elastic strain, C is the primary integration constant. Based on the above contents, the axial creep equation of the new nonlinear creep model can be obtained as follows.

Model Parameters Identification
In general, there are two methods for model parameter identification [14], i.e., the curve analysis method and the numerical iteration method. This work uses the direct curve analysis method to identify the instantaneous parameters, and then uses the numerical iteration method (also named numerical fitting) to determine the creep parameters. The parameters of E 1 , P, and σ s could be identified using the instantaneous strain data in Table 2. Then, E 2 and η 1 could be obtained by fitting the unloading creep curves using Equation (8). Finally, η 2 , A, B, and C can be solved by fitting the visco-plastic creep curve with Equation (14) based on the Levenberg-Marquardt optimization algorithm. All the model parameters of the specimen under six deviatoric stress levels are listed in Table 3. Figure 8 presents the experimental curves and theoretical curves of the proposed model at every stress level. Additionally, the evolution curves of the damage variable with time under two typical deviatoric stress levels are plotted in Figure 9. Table 3. Model parameters of the testing specimen under six deviatoric stress levels.  It could be found that the theoretical curves fit well with the experimental curves, correlation coefficients are between 0.72-0.93. The damage variables are in good agreement with the increasing trend of strain. These suggest that the proposed nonlinear creep damage model could legitimately describe the elasto-visco-plastic strain characteristics of dolomitic limestone under cyclic incremental loading and unloading.

Model Parameters Analysis
In order to expand the application scope of the new model and deepen the understanding of this model, a parameter analysis is required. As the Hookean substance, St. Venant substance, and Kelvin body are well-known models, herein, only the analysis of the parameter A, B, and C within the newly proposed viscous body is carried out, the variation curves of visco-plastic strain with time under different model parameters are drawn in Figure 10. Here, η 2 and σ 1 − σ 3 are taken as 300 GPa·h and 50 MPa, respectively.
It could be found that the value of parameter A only affects the initial value of viscoplastic strain (see Figure 10a,b), while parameters B and C have an important impact on the evolution law of visco-plastic creep curve. Specifically: (1) When the value of parameter C is among 1−4 × 10 −4 , the model curve is characterized by an accelerated creep process, and the creep rate decreases with the increase of parameter B, as presented in Figure 10c. While parameter B is constant, the accelerated creep rate increases with the increase of parameter C (see Figure 10e). (2) When the value of C is among −4−(−1) × 10 −4 , the model curves are featured by attenuation creep or steady-state creep. As the value of parameter C is relatively large, the model curve mainly shows an attenuation creep process, and the creep rate increases with the increase of parameter B, i.e., the duration of the attenuation creep stage decreases, as revealed in Figure 10d. While the value of parameter C is small and the value of parameter B is large, the visco-plastic model curves mainly show a steady-state creep process, and the rate decreases with the increase of parameter C, as presented in Figure 10f.
Moreover, the theoretical model in this paper is suitable for the creep law description of many rocks under the condition of loading and unloading in engineering practice. The creep parameters of different rocks could be obtained by using the creep testing results combined with one of the model parameter identification methods. Then, a more accurate creep theoretical model with definite creep parameters of the studied rocks could be applied to the long-term stability analysis and engineering design.

Advantages and Disadvantages of the Developed Model
The classical multi-component models, such as the Maxwell model, Kelvin model, Bingham model, Burgers model, as well as Xiyuan model, cannot describe the nonlinear accelerated rheological behavior [9,14]. Moreover, they cannot reveal the progressive damage properties of rock during the creep process, especially the accelerated creep stage. The developed model includes a newly proposed damage viscous body, which could describe the evolution of nonlinear visco-plastic strain in the accelerated creep stage.
Besides, the damage variable is usually defined as an equation concerning the elastic modulus in previous studies [21,29]. However, the elastic modulus changes with time, and the initial elastic modulus is difficult to define during the creep process [21]. Here, during compressive creep test of rock subjected to cyclic loading and unloading, elastic strain and plastic strain could be measured and calculated, which is consistent with the important principle of the definition of damage variable: the parameters should be easy to measure and easy to establish contact with macro mechanical properties [29].
In addition, it can be seen from Figure 10 that the new visco-plastic model would reflect different evolutionary trends of creep strain under different testing conditions and has wide applicability. Thus, the multi-element creep model can be used universally for describing the creep behaviors of rock subjected to cyclic incremental loading and unloading.
However, the disadvantage of this model is that, as shown in Figure 10a, there was a short-lived fluctuation (i.e., the value of the damage variable decreased) in the initial stage of the damage variable curve. This is mainly because that the visco-elastic strain curve is always attenuated creep process, while the visco-plastic strain curve shows steady-state creep (i.e., red curve Figure 10f) under the condition of deviatoric stress being 51 MPa, the initial creep rate of the former is relatively large, whereas the initial creep rate of the latter is relatively small, causing the growth rate of visco-elastic strain will be greater than that of visco-plastic stain in a short time during the beginning stage of creep process. In fact, as shown in Figure 10d, this phenomenon would not occur when the visco-plastic strain curve also belongs to the type of attenuation creep. Subsequent studies should pay more attention to and explore the definition of more satisfactory damage variables.

Conclusions
(1) The triaxial compression creep behaviors of dolomitic limestone subjected to cyclic incremental loading and unloading manifest that the instantaneous deformation occurs firstly during every deviatoric stress level and then creep deformation occurs. When the loading deviatoric stress is not more than 75% of the conventional triaxial compression strength, only transient creep stage and steady-state creep stage occur successively. However, the rock undergoes accelerated creep deformation while the loading stress reaches 82.35% of the instantaneous compressive strength. (2) The separation of elasto-visco-plastic strain component shows that rock deformation is mainly characterized by instantaneous strain while the deviatoric stress is not more than 75% of the conventional compressive strength, especially the irrecoverable instantaneous plastic strain accounts for 37.97-60.27% of total instantaneous strain and should not be ignored during rock engineering design. Greater deviatoric stress contributes to more obvious creep deformation, the visco-elastic strain increases linearly with deviatoric stress, More importantly, the irrecoverable visco-plastic strain increases exponentially with the deviatoric stress increasing from 26.47% to 82.35% of the compressive strength, and finally leads to the occurrence of accelerated creep and rock failure. The obtained findings can help our understanding of the elasto-visco-plastic deformation behavior of dolomitic limestone and provide a more accurate creep theoretical model for long-term stability analysis and engineering design.
Author Contributions: Conceptualization, X.W. prepared the manuscript and conducted the investigation, L.S. and C.X. were responsible for the methodology, G.H. carried out the data processing, L.S. acquired the funding, C.X. and Z.Z. reviewed and edited the manuscript. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.