#### Abstract

In order to study the creep behavior of the surrounding rock of Hengda coal mine in Fuxin under different temperatures, the triaxial creep test of sandstone is carried out by the MTS815.02 test system. The relationship between damage variables and temperature is constructed based on the Weibull distribution of the meso-probability voxel intensity. Aiming at the nonlinear characteristics of rock creep, a nonlinear viscous pot element and a nonlinear spring element are proposed. The two linear viscous pot elements and one linear spring element in the Nishihara model can be replaced separately. Thus, an unsteady parameter creep model is established. The comparison between the Nishihara model curve and the model and the experimental curves in this article has been added to the article. Furthermore, the superiority of this model can be proved. The results show that the established variable-time aging creep model not only can describe the rock attenuation creep and stable creep deformation characteristics but also can make up for the shortcomings of the traditional creep model that cannot describe the accelerated creep characteristics. Moreover, it predicts the development law of creep deformation well. The model is in good agreement with the test curve, which shows the correctness and rationality of the model. It has guiding significance for actual engineering support and prediction of long-term deformation of surrounding rock.

#### 1. Introduction

After the deep coal mining, the surrounding rock of the roadway undergoes time-dependent deformation under the action of the supporting body, and this time-related deformation is generally called rheology [1]. For the study of rheological properties of the rock under external loads, a rheological model describing rheological properties is established, which has practical guiding significance for predicting surrounding rock deformation and design support scheme [2, 3]. As the mining depth increases, the mechanical properties of rock are no longer the linear variation law when the shallow part is present, and the rock mechanical properties show nonlinear variation characteristics [4]. Not only soft rock produces rheology but also hard rock (such as sandstone, marble, and granite) has obvious rheological properties [5]. In order to further study the creep deformation characteristics and failure mechanism of rock under high confining pressure, high temperature, high ground stress, and strong disturbance condition, it is necessary to carry out an indoor creep test on the rock. Furthermore, a suitable creep constitutive model can be constructed by experimental conditions and data, which is of great significance for studying the creep mechanism of materials [6–8].

In recent years, scholars had done a lot of research on the rheological properties of rocks, the mechanism of deformation and failure, and how to construct a suitable creep model [9–11]. The establishment of the model includes the following methods: (1) a nonlinear creep model was constructed by the series-parallel connection of basic components and the nonlinearization of components; (2) the rock damage creep constitutive model was constructed by combining thermodynamic and damage mechanics theories; (3) through a large number of experimental research results, the empirical creep model formula was summarized, and the model parameter determination method was proposed; (4) a fractional-order creep constitutive model was established by transforming the integer-order component model into a fractional-order component model. Based on the multiload creep test of sandstone with different initial damage levels, a new nonlinear creep damage model was proposed. The model not only described three typical creep stages but also showed the effect of initial damage on creep failure stress [12]. Wang et al. [13] constructed a micromechanical damage-friction coupling model and applied the model to the triaxial compression test and triaxial creep test of simulated granite. The numerical data could be used to predict the experimental data well. The fractional differential element was used to replace the viscous element in the traditional Nishihara model, and the nonlinear creep parameter creep model of rock was established [14]. The shear creep behavior of the key unit rock of the potential sliding surface of the landslide was studied by the shear creep test. A new plastic nonlinear model (PFY model) was proposed to characterize the progressive failure creep characteristics reflected in the limiting shear creep process of rock [15]. Hadiseh and Rassoul [16] conducted axial compression and triaxial creep tests on the salt rock in salt wells in southern Iran. The mechanical properties of rocks were studied, and a nonlinear creep model of rock was constructed. Using the joint inclination angle *α* as a variable, Liu and Zhang [17] established numerical models of different joint inclination angles using discrete element numerical methods and studied the shear mechanical properties of the sawtooth joint under creep conditions. Aiming at the state of the direct shear test, Lin et al. [18] used Kachanov’s creep damage law to describe the time characteristics of the shear strength of the rock in the accelerated creep stage. A nonlinear viscoplastic element based on time-based shear strength is established. This model can reflect the mechanical properties of the rock in the three creep stages.

Although the model in the above study can describe the acceleration of creep deformation of the rock into the cerebral blood vessels, the creep equation was finally formed to have a low degree of fitting for accelerated creep deformation. Therefore, it was very important to choose a simple transformation of creep parameters into time-dependent functions and then derived nonlinear creep equations. In this paper, the creep deformation characteristics of the surrounding rock of the roadway in Hengda coal mine of Liaoning Province under different temperatures were analyzed. The Nishihara model creep parameters were all converted into time-dependent functions. The rheological theory was used to recalculate the relationship between temperature and stress and time. A time-dependent creep constitutive model of deep surrounding rock under temperature-stress coupling was constructed. Finally, the correctness of the creep model was verified by comparing the experimental data with the model curve.

#### 2. Establishment of Time-Dependent Model

##### 2.1. Establishment of Variable Parameter Model in One-Dimensional State

The Nishihara model consists of an elastomer, a viscoelastic body, and a viscoplastic body [19]. The one-dimensional model expression for the Nishihara model is as follows [20–22]:

When *σ* < *σ*_{s}

When *σ* ≥ *σ*_{s}where *σ* is the stress, *E*_{0} is the elastic modulus of the elastomer, *E*_{1} is the elastic modulus of the viscoelastic body, *η*_{1} is the viscosity coefficient of the viscoelastic body, *η*_{2} is the viscosity coefficient of the viscoplastic body, and *σ*_{s} is the yield stress.

The mechanical properties of roadway surrounding rock in a complex deep geological environment cannot be explained by conventional mechanics theory. At this time, the mechanical and the creep properties of surrounding rock have obvious nonlinear characteristics. In order to describe this nonlinear creep deformation characteristic, the creep parameter can no longer be used as a fixed value. It is a function of time variation. The creep parameter and time relationship function can be expressed as follows:where *Z* is a variable representing all creep parameters.

It is assumed that the variation law of damage variable and time product satisfies the following equation [23]where *α* is the parameter that characterizes the internal damage of the rock, and *m* is the nonuniformity of damage distribution of meso-probability voxels.

The distribution of microdefects in rock materials has significant self-similar characteristics, which can be characterized by the fractal dimension. The Weibull distribution function parameter *m* is the main parameter that affects the shape. Chen et al. [24] proposed that the distribution function parameter *m* has a linear relationship with the fractal dimension. The rocks in the deep strata are subjected to high temperatures, and the temperature on the microscopic scale mainly affects the interconnection of the fracture network along with the crystal. Wu [25] found that the fractal dimension describing the connection of the salt rock fracture network has a quadratic function relationship with temperature. It is assumed that the nonuniform distribution parameter *m* inside the rock is connected with the fracture network. Generally, the room temperature is 25°C. The relationship between the nonuniformity coefficient *m* and the temperature is as follows [26]:where *a*, *b*, and *c* are coefficients related to temperature.

The relationship between damage variable and temperature and time is obtained by combining equations (3) and (4) as follows:

Liu and Zhang [17] believe that the creep parameters of rocks in the deep complex address environment are related to time. It is believed that the relationship between the damage variable, time, and creep parameters of the rock satisfies the following equation. It can be defined that the rock creep parameter *Z* has the following relationship with the damage variable *D*.

Equation (5) is substituted into Equation (6). The relationship between rock creep parameters and time is obtained as follows:

The total strain *ε* satisfies the following formula under one-dimensional stress [21]:where *ε*_{e} is the elastic strain, *ε*_{ve} is the viscoelastic strain, and *ε*_{vp} is the viscoplastic strain.

The instantaneous elastic strain is only related to stress, not time. The change law of elastic modulus is only related to stress. The damage variable of elastic modulus satisfies equation (9) [27]. The relationship between the damage variable and time in equation (5) can be transformed into the relationship between the damage variable and stress. Therefore, the relationship between the elastic modulus *E*_{0} and the damage variable can be expressed as follows:where *α*_{0} is the damage influence factor related to stress; *a*_{0}, *b*_{0}, and *c*_{0} are coefficients related to temperature.

The elastic strain of the rock is not affected by time, but it is affected by stress. In the one-dimensional state, the elastic strain *ε*_{e} of the rock satisfies the following equation:

In order to obtain the creep deformation of the viscoelastic body and make the calculation process conform to the laws of mathematics, this article will only improve the viscous pot element in the viscoelastic body and will not improve the spring element. In the one-dimensional state, the viscoelastic strain *ε*_{ve} of the rock satisfies the following equation:where *α*_{1} is the time coefficient of influence of the viscoelastic body; *a*_{1}, *b*_{1}, and *c*_{1} are coefficients related to temperature.

Separate the variables to calculate the definite integral of Equation (11), and the creep deformation of the viscoelastic body can be obtained as follows:

The exponential function in Equation (12) can be expanded by the Taylor series as follows:

By integrating Equations (13) and (14), the following can be obtained:

The viscoelastic deformation is as follows:

In the one-dimensional state, the viscoelastic strain *ε*_{vp} of the rock satisfies the following equation:where *α*_{2} is the time influence coefficient of the viscoplastic body, and *a*_{2}, *b*_{2}, and *c*_{2} are coefficients related to temperature.

The exponential function is expanded by using the Taylor series. The integral solution of Equation (16) is solved to obtain a viscoplastic strain expression as follows:

##### 2.2. Establishment of Variable Parameter Model in Three-Dimensional State

In the actual project, the surrounding rock of the roadway is in a three-direction stress state, and the above model cannot describe the multidirectional force creep characteristics. This requires transforming the one-dimensional model into a three-dimensional model [28, 29].

However, the total strain *ε*_{11} of the Nishihara model in the three-direction stress state satisfies the following equation:where is the elastic strain in the three-dimensional state, is the viscoelastic strain in the three-dimensional state, and is the viscoplastic strain in the three-dimensional state.

The elastic strain in the three-dimensional state is not affected by time but by the stress state. The elastic strain can be expressed as a function of the stress state by the elastic model in the one-dimensional state as follows:where *K* is the bulk modulus, and *G*_{0} is the elastomer shear modulus.

The expression of viscoelastic strain affected by time in the three-dimensional state is as follows:

*A*'(*t*) is the corresponding term of *A*(*t*) in the three-dimensional stress state. The function expression is as follows:where is the time coefficient of influence of the viscoelastic body under the three-dimensional stress state.

The viscoplastic strain in the three-dimensional state cannot be directly converted by analogy. It is also affected by the plastic potential and the yield functions. Therefore, the viscoplastic strain in the three-dimensional state can be expressed as follows:

When *σ* ＜ *σ*_{s}

When *σ* ≥ *σ*_{s}where *Q* is the plastic potential function, and *n* is a constant. Generally, *n* can be taken as 1 [30]. *F*_{0} is the initial reference value of the yield function of the rock.

It can be assumed that the initial yield function value of the rock is 1. Moreover, according to the flow law in plasticity theory, Equation (21b) can be expressed as follows:

In general, the yield function selects the generalized Drucker–Prager yield function. The function expression is as follows:where *I*_{1} is the first invariant of stress, *J*_{2} is the second invariant of stress bias, and *ξ* and *k* is the test parameter.

The test parameters *ξ* and *k* are functions of internal friction angle and cohesion. The function expressions are as follows:where *φ* is the angle of internal friction, and *C* is the cohesion.

The internal friction angle and cohesion will also deteriorate under the action of time and temperature. The function expressions are as follows:

Therefore, the viscoplastic strain of the rock can be obtained as follows:

The time-dependent constitutive equation for obtaining rock in a three-dimensional state is as follows:

When *σ* < *σ*_{s}

When *σ* ≥ *σ*_{s}

#### 3. Creep Test

##### 3.1. Creep Test Plan

In this paper, the surrounding rock (sandstone) of the roadway in Hengda coal mine of Fuxin is shown in Figure 1. The selected sandstone sample has a buried depth of 500–750 m. The appearance of the rock sample is dark gray; the structure is relatively uniform; and the texture is relatively hard. There are no microcracks and bedding visible to the naked eye. Due to the high in situ stress everywhere in the rock, the maximum horizontal in situ stress reached 25 MPa, and the vertical in situ stress was 10.56 MPa after actual measurement, so this area belongs to a high in situ stress environment. In order to facilitate the experiment, the confining pressure value is selected as 10 MPa. In order to simulate the effect of high temperature on the creep characteristics of rock, the temperature is set to 100°C and 200°C.

According to the requirements of the International Rock Mechanics Society for standard test pieces, the surrounding rock was made into a cylindrical sample with a height of 100 mm and a diameter of 50 mm. It must be ensured that the nonparallelism and unevenness of both ends of the specimen are less than 0.05 mm. Prior to the test, Vaseline was evenly applied to both end faces of the rock to eliminate the end effect during the test. The test equipment used was the MTS815.02 rock test system (shown in Figure 2). Its maximum confining pressure was 100 MPa, and the precise range of force measurement was −1%∼+1%, which met the requirements of this test.

In this paper, the indoor triaxial creep test was carried out by the single specimen gradual loading method. First, the confining pressure was applied to a predetermined value, and the confining pressure was selected to be 10 MPa. The temperature of the creep test is set to 100 and 200°C. The stress levels were 50, 60, 70, and 80 MPa. After the confining pressure was stabilized, the axial pressure was applied. The applied load rate was set to 500 N/s. The temperature is loaded to the predetermined temperature at a rate of 0.5°C/s. When applying axial stress, it must be ensured that the confining pressure had been changed within a controllable range of the predetermined value. After the stress level creep deformation entered the stable creep, the next level of load application began. This cycle was repeated until the rock sample was destroyed. Finally, the test data was saved at intervals of 3 s. After unloading the confining pressure and the axial pressure, the sample is taken out and stored.

##### 3.2. Analysis of Creep Characteristics

The single test piece is gradually loaded to complete the creep test, and the creep data under each load is affected by the historical load. Therefore, the data need to be processed by Chen’s superposition method. The axial creep deformation-time curve of the surrounding rock of the roadway under different temperatures is shown in Figure 3.

**(a)**

**(b)**

It can be seen from Figure 3 that the creep deformation curves of rock under different temperatures are basically the same. Under the action of the temperature of 200°C, the deformation failure occurs after the fourth stage of load. Under the action of the temperature of 100°C, the deformation failure occurs after the fifth stage of load. The instantaneous and creep strains increase with the increase of stress level, and the ratio of instantaneous strain to total strain first decreases and then increases. This is due to the compaction of the internal voids of the rock under the initial stress level. Taking the temperature of 200°C as an example, under the action of low stress (50 MPa), the rock deformation has only transient strain and decay creep deformation. Finally, the rock creep deformation rate decays to zero. Under medium stress (60 and 70 MPa), rock creep has two kinds of attenuation creep and stable creep. At this time, the creep rate does not decay to zero but decays to a stable value and enters the stage of creep to be stabilized. The creep deformation of the rock keeps this creep rate continuously deformed. Under high stress (80 MPa), the rock appears to accelerate creep deformation. After the first two creep stages, the rock enters the accelerated creep stage. The creep rate accumulation increases, and finally, the internal crack of the rock penetrates to form a fracture surface. The last creep time of the rock decreases with the increase in temperature. When the temperature is 200°C, the last stage creep time is 6.32 h, and the final stage creep time is 7.89 h when the confining pressure is 100°C. This shows that the decrease of temperature effectively increases the rock bearing capacity and increases the creep damage time.

#### 4. Verification of Creep Aging Constitutive Model

Before verifying the rock creep aging model, the long-term strength of the rock under different confining pressures needs to be determined [31]. The long-term strength of the rock can be determined by isochronous stress-strain curves. The isochronous stress-strain curve under the action of the temperature of 200°C is shown in Figure 4. The isochronous stress-strain curve is the curve drawn by the strain at different periods under the same stress load. Moreover, every isochronous stress-strain curve has an obvious turning point. Before the turning point, the curve basically showed a linear increasing trend. After the turning point, the curve basically shows a nonlinear increasing trend. Generally, the stress value corresponding to the turning point can be used as the long-term strength of the rock.

It can be seen from Figure 4 that the isochronous stress-strain curve of the rock is a cluster of divergent broken line segments. The rock stress and strain show a linear change before the divergence point. After the divergence point, the rock stress and strain basically show nonlinear linear variation. Therefore, the divergence point corresponding stress value can be used as the long-term strength value of the rock. Under the action of the temperature of 200°C, the long-term strength of the rock is 60 MPa.

The improved creep model was fitted by the least-squares method, and the creep parameters are obtained as shown in Table 1 (taking 200°C as an example).

By substituting parameters of different stress levels into the model, the model curves of the rock under different stress levels are obtained. The comparison between the model curve and the experimental data is shown in Figure 5(a). In the same way, the model curve and test curve of rock under the temperature of 100°C can be obtained as shown in Figure 5(b).

**(a)**

**(b)**

It can be seen from Figure 5 that the rock model curve and the experimental data have good fitness, and the correlation coefficients are all above 0.90. It is shown that the variable parameter time-dependent creep model established in this paper not only can describe the rock attenuation creep and stable creep deformation characteristics but also can make up for the shortcomings of the traditional creep model that cannot describe the accelerated creep characteristics. Moreover, it is a good predictor of the development of creep deformation. At the same time, the model is in good agreement with the experimental curve, which indicates the correctness and rationality of the model. It has guiding significance for actual engineering support and prediction of long-term deformation of surrounding rock.

In order to verify the correctness of the model established in this paper better, the creep test curve under different temperature conditions and the curve of the traditional Nishihara model and the model curve of this paper are compared and analyzed as shown in Figure 6.

**(a)**

**(b)**

It can be seen from Figure 6 that the parameters of the traditional creep model are unsteady based on the degradation characteristics of the strength parameters and creep parameters. Compared with the traditional Nishihara model, the improved creep model can better describe the accelerated creep deformation characteristics of rocks.

#### 5. Conclusion

In this paper, the creep behavior of the surrounding rock of Hengda coal mine in Fuxin under different temperatures is analyzed. The creep parameters of the Nishihara model are transformed into time-dependent functions. Furthermore, a time-dependent creep constitutive model of deep surrounding rock under temperature-stress coupling is constructed.

The instantaneous strain and creep strain increases with the increase of stress level, and the ratio of instantaneous strain to total strain first decreases and then increases. This is due to the compaction of the internal voids of the rock under the initial stress level.

The rock model curve and the experimental data had good fitness, and the correlation coefficients are all above 0.90. It is shown that the variable parameter time-dependent creep model established in this paper not only can describe the rock attenuation creep and stable creep deformation characteristics but also can make up for the shortcomings of the traditional creep model that cannot describe the accelerated creep characteristics.

It is a good predictor of the development of creep deformation. The model is in good agreement with the experimental curve, which indicates the correctness and rationality of the model. It has guiding significance for actual engineering support and prediction of long-term deformation of surrounding rock.

#### Data Availability

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

#### Disclosure

The authors declare that this paper has been presented as a preprint in Research Square.

#### Conflicts of Interest

The authors declare that they have no financial and personal relationships.

#### Acknowledgments

The authors gratefully acknowledge the National Natural Science Foundation of China (51774173) for funding this work. Special thanks are due to Professor Xiangzhi Yin for providing technical support.