Creep damage model of rock with varying-parameter under the step loading and unloading conditions

The creep characteristics of rock under step loading and unloading conditions were investigated in this study. Based on the generalized Burgers model, the total strain of rock was decomposed into elastic, viscoelastic, varying-parameter viscoelastic, and viscoplastic strains considering the damage. The four strains were connected in series to establish a new varying-parameter creep damage model that can characterize the creep characteristics of rock under step loading and unloading conditions as well as identify and verify the model parameters. The study results showed that the varying-parameter creep damage model could better describe the creep characteristics of rock under step loading and unloading conditions, significantly the non-linear both the strain and time of attenuation creep and accelerated creep. The model fitting curve was highly consistent with the experimental data, and the correlation coefficient R2 was greater than 0.98, which thoroughly verified the accuracy and rationality of the model. These findings can provide theoretical support for analyzing the deformation and long-term stability of rock and soil.

Rock creep has been a prominent research hotspot for experts and scholars in geotechnical engineering. Rock creep is closely related to rock and soil deformation and long-term stability 1 . In the development process of slope excavation and tunnel excavation, the surrounding rock is typically in loading and unloading states. The loading and unloading effects seriously weakens the strength of the rock mass and induces secondary disasters 2,3 . Therefore, understanding the rock creep model under step loading and unloading conditions is theoretically significant and practically valuable for slope engineering design, landslide warning, and tunnel deformation prevention.
Scholars have achieved fruitful advancements in studying rock creep models, which are mainly categorized into empirical and component combination models [4][5][6] . Empirical models are based on experimental phenomena and laws [7][8][9][10] ; component combination models are based on essential components, such as elasticity, viscosity, and plasticity, and are formed by series and parallel methods 11,12 . Wang et al. 13 established a six-element composite creep model by analyzing the loading and unloading creep curves of limestone. Song et al. 14 used saturated red sandstone as the step loading and unloading creep test object and constructed a freeze-thaw-damage creep model. Deng et al. 15 utilized Q3 loess as a research object, introduced a fractional element model based on the Nishihara model, and deduced the constitutive equation of rock under the cyclic loading and unloading conditions. Yang et al. 16 conducted graded loading and unloading uniaxial creep tests on sandstone and established a creep damage model. Based on the triaxial cyclic loading and unloading creep phenomenon of rock, Yang et al. 17 separated rock strain by viscoelastic-plasticity and constructed a creep model considering the damage. Wang et al. 18 conducted triaxial compression tests and unloading confining pressure tests to determine the mechanical properties of sandstone under loading and unloading conditions. Despite these efforts, no study has investigated the rock creep model under step loading and unloading conditions, which considers the damage of the accelerated creep stage and the nonlinear function of the viscosity coefficient and time.
Based on previous research results, this study utilized the sandstone of a tunnel floor as the test object, introduced the time-dependent viscosity coefficient into the viscoelastic body, and introduced the damage variable into the viscoelastic body. The total strain of the rock was decomposed into elastic, viscoelastic, variable parameter viscoelastic, and viscoplastic strains considering the damage. The four strains were connected in series to

Analysis of sandstone creep test data under step loading and unloading conditions
Yang et al. 19 utilized the sandstone of a tunnel floor as the test object: the sample size was Φ50 mm × 100 mm, and the RLW-2000 triaxial rheological test system was used to conduct a step loading and unloading creep test on the sandstone sample. The confining pressures were 4 MPa, 8 MPa, and 12 MPa, respectively. The levels of deviator stress during the test are shown in Table 1, and the creep test curve of graded loading and unloading sandstone is shown in Fig. 1. Figure 1 indicates that instantaneous elastic strain occurs initially at low-stress levels (level 1 to level 4). As the loading time increases, the axial strain rate of sandstone gradually decreases, and the axial strain increases nonlinearly, all tending to a particular value. With an increase in the unloading time, the axial strain rate and axial strain of the sandstone gradually decrease, and both tend to a particular value. The sandstone sample exhibits attenuation and constant velocity creep stages under the step loading and unloading conditions. At high-stress levels (level 5), the sandstone has a nonlinear accelerated creep phase. In summary, the creep stages of sandstone under step loading and unloading conditions include attenuation creep, constant velocity creep, and accelerated creep stages.

Establishment of the varying-parameter creep damage model
The classic Nishihara and Burgers creep models can describe the deformation characteristics of rock attenuation creep and constant velocity creep. However, the models cannot describe the nonlinear law of strain and time in the accelerated creep stage. Therefore, to accurately describe the full-stage creep characteristics of the rock under step loading and unloading conditions, a varying-parameter creep damage model is established based on

Elastomer. Lin et al. 23 established a one-dimensional creep equation of an elastometer,
where σ is the stress of the elastometer, k 1 is the elastic modulus of the elastometer, and ε 1 is the strain of the elastometer.
The creep equation of an elastometer under three-dimensional stress can be written as follows 24 : where v is the Poisson's ratio of rock, S ij is the deviation stress tensor, σ m is the average stress, δ ij is the Kronecker tensor, G 1 is the shear modulus of the elastometer, and K is the bulk modulus of the elastometer. The deviatoric stress tensor in the principal stress space can be expressed as Viscoelastic body. The viscoelastic body is also known as a Kelvin's material, and its constitutive equation where ε 2 is the strain of the Kelvin's material, k 2 is the elastic modulus of the Kelvin's material, and η 1 is the viscosity coefficient of the viscoelastic body (Kelvin's material). A one-dimensional creep equation of the viscoelastic body is established as follows 26 : By extending to the three-dimensional space analysis, the three-dimensional creep equation of the viscoelastic body can be obtained as 27 where G 2 is the shear modulus of the viscoelastic body.
Assume that the unloading occurs at t 0 , Eq. (4) can be determined using where A is the integral constant. When t = t 0, Eq. (5) can be written as follows: (1)  Varying-parameter viscoelastic body. Zhang et al. 28 assumed that the viscosity coefficient is a power function related to time; function η 2 (t) of the viscosity coefficient related to time is given as follows: where η 20 is the initial viscosity coefficient of the varying-parameter viscoelastic body, and λ is a constant.
According to the one-dimensional constitutive equation of a Kelvin's material, a one-dimensional creep equation of the varying-parameter viscoelastic body is established as follows: where ε 3 is the strain of the varying-parameter viscoelastic body, and k 3 is the elastic modulus.
A one-dimensional creep equation of the varying-parameter viscoelastic body is established by integrating Eq. (13) as follows: By extending to the three-dimensional space analysis, the three-dimensional creep equation of the varyingparameter viscoelastic body can be obtained as where G 3 is the shear modulus of the varying-parameter viscoelastic body.
If t 0 is the unloading time, Eq. (13) can be written as follows: Integrate Eq. (16) as follows: where B is the integral constant. When t = t 0 , Eq. (14) can be written as follows: Based on Eq. (17) and Eq. (18), the integral constant B can be written as follows: Based on Eq. (17) and Eq. (19), the one-dimensional creep equation of the varying-parameter viscoelastic body under unloading conditions can be established as where ε 4 is the strain rate of the damaged viscoplastic body, and η 3 (t) is the function of the viscosity coefficient related to time.
Based on the results of numerous rock creep damage tests, the damage variable D took the form of a negative exponential function related to time during rock creep [29][30][31][32] . In this study, the damage variable is reduced to Eq. (23).
where α is the coefficient related to the properties of rock materials. where ε 4 is the strain of the damaged viscoplastic body, and η 3 is the initial viscosity coefficient of the damaged viscoplastic body.
The three-dimensional creep equation of the damaged viscoplastic body can be obtained as: where F 0 is the initial reference value of the rock yield function, F is the rock yield function, Q is the plastic potential function, and φ(·) is the power function form. When F ≥ 0, based on Eq. (25) and associated flow standards, the following relationship can be obtained: where and σ s is the yield strength of the rock. Under the conventional triaxial compression test conditions (σ 2 = σ 3 ), J 2 can be obtained as Assume that the initial reference value of the rock yield function F 0 = 1, then the three-dimensional creep equation of the damaged viscoplastic body can be written as follows: Establishment of the varying-parameter creep damage model. According to the series character of components, the following relationship can be obtained: where σ 1 is axial compression, σ 2 is confinement pressure, σ 3 is confinement pressure, and S 11 is axial deviatoric stress tensor.

Model parameter identification and validation
To verify the accuracy and rationality of the variable parameter creep damage model, the author used MATLAB to identify the parameters of the varying-parameter creep damage model based on the sandstone creep test data under step loading and unloading conditions. To simplify the parameter identification process, letters are used to replace the complicated items in the theoretical model; the following relationship can be obtained:  Tables 2 and 3. Owing to many images, only the model curve and test data when the confining pressure is 4 MPa are listed (as shown in Fig. 4).

Scientific Reports
It can be seen from Fig. 4 that the varying-parameter creep damage model can better simulate the three stages of rock creep Moreover, the scattered points are approximately on both sides of the fitting curve, indicating that the model can achieve results under different loading stress conditions and different lithologies. The better fit results demonstrate the correctness and applicability of the model.
As shown in Fig. 4a-d, regardless of the stage (loading or unloading), the model curve is highly consistent with the changing trend of the test data and can better describe the creep mechanical behavior of rocks. Hence, the model is also suitable for studying the law of stable creep deformation of rocks.
It can be observed from Fig. 4e that the accelerated creep phase occurs quickly, and the creep curve rises rapidly. The model can also capture the creep deformation law adequately, especially as the description of the inflection point is more accurate. The model can describe the nature of the rock strain changing drastically with time in the accelerated creep stage, which fully reflects the superiority of the established model for describing the nonlinear accelerated creep characteristics of rock.
Overall, the varying-parameter creep damage model has a better identification effect and a higher fitting accuracy under different confining pressures and deviator stresses. The model can better characterize the creep properties of rock during loading and unloading, especially the nonlinear changes of strain and time in the attenuation creep and accelerated creep stages. The fitting curve is highly consistent with the experimental data, the correlation coefficient R 2 is greater than 0.98, and the parameters conform to knowledge of physics www.nature.com/scientificreports/ www.nature.com/scientificreports/ and mechanics; therefore, the accuracy and rationality of the model are thoroughly verified. consequently, the varying-parameter creep damage model developed in this study can be used to predict the creep deformation of rock under long-term loading and unloading stress.

Conclusion
The total strain of rock is decomposed into elastic, viscoelastic, varying-parameter viscoelastic, and viscoplastic strains considering the damage. The four strains were connected in series to establish a new varying-parameter creep damage model that can characterize the creep characteristics of rock under step loading and unloading conditions. The varying-parameter creep damage model could better describe the creep characteristics of rock under step loading and unloading conditions, significantly the non-linear both the strain and time of attenuation creep and accelerated creep. The model fitting curve was highly consistent with the experimental data, and the correlation coefficient R 2 was greater than 0.98, which thoroughly verified the accuracy and rationality of the model.

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