# Experiment and Numerical Simulation on Creep Mechanical Behaviors of Mudstone under Unloading Condition

Authors: Xudong Feng
Journal: Geofluids
 (2022)

---

## Abstract

The deep soft rock with complex occurrence conditions is easy to produce significant creep deformation under the influence of high in situ stress, high temperature, and disturbance, which has a great influence on the safety and long-term stability of structures with long service life, such as underground chambers and nuclear waste storage rooms. To study the creep effect of surrounding rock of deep chamber under the condition of excavation and unloading, creep tests of deep mudstone over 1000 meters in unloading confining pressure were carried out in this work. The experimental results show that high confining pressure has a significant inhibition effect on the creep behavior of rock, while high axial stress aggravates the creep behavior, which is manifested by the increase of creep deformation and the shortening of total creep time. When the axial stress is higher than the yield strength of the rock mass, the rock mass enters the tertiary creep stage (which is characterized by large deformation and instability failure in a short period of time). An improved nonlinear viscoelastic-plastic rock creep model, which can simulate three creep stages, is established (well verified). In addition, the creep deformation of surrounding rock after excavation is simulated based on this model, and the simulation results are consistent with the engineering practice: after 165 days of excavation, there is still a large deformation rate.


---

## Body

## 1. Introduction

Creep property of rock is one of the mechanical properties that cannot be ignored in engineering. In the design of underground structures, not only the instantaneous deformation under load but also the creep deformation over time should be considered, especially in the construction of deep soft rock with complex occurrence conditions [1]. Under the influence of high in situ stress, high temperature, and human disturbance, deep soft rock is prone to creep deformation, and joints and fissures are usually well developed under long-term action [2, 3]. A large number of studies have proved that surrounding rock deformation has a significant time effect [4–8]. How to predict and solve the deformation problem caused by rock creep of underground structures has become an urgent problem to be solved in underground support engineering. Therefore, the study of creep mechanical properties of deep soft rock is crucial to the safety and long-term stability of deep underground structures with long service life, such as underground chambers and nuclear waste storage rooms [9, 10].Scholars have carried out a large number of studies on the creep behavior of rock, which can be mainly divided into two categories: (i) indoor creep tests under single and triaxial graded loading and unloading [11, 12]. (ii) Theoretical research on the combination model of the corresponding components (elastic components, viscous components, and plastic components) was established according to the test results [13–16]. In terms of experimental research, Yang et al. [17] conducted a triaxial rheological study on diabase under different stress paths, believing that there would be three stages of deceleration, constant velocity, and accelerated creep under fracture stress level. Heap et al. [18] carried out triaxial creep test on Darley Dale sandstone under saturated water condition, and the test results showed that the new damage generated in the primary and secondary creep process was mainly dominated by axial crack growth parallel to the maximum principal stress. In addition, the axial strain rate is largely dependent on the applied differential stress, and a reduction of only 10% of the differential stress can result in a decrease in the strain rate of more than two orders of magnitude. Brantut et al. [19, 20] conducted triaxial compression creep test of granite and basalt and discussed the influence of compressive deviational stress and pore pressure on creep duration and deformation rate. In addition, by comparing the relationship between the constant strain rate and the applied stress and inelastic strain in creep test, it is found that the mechanical work required by creep deformation is small. Through creep test, Fabre and Pellet [21] found that there is a certain stress threshold between second-order creep and third-order creep. When the axial loading stress is less than the threshold stress, only the initial rock creep occurs. When the axial loading stress is greater than the threshold stress, secondary and tertiary creep stages occur. In terms of theoretical model research, scholars have constructed different nonlinear creep models to describe the three stages of rock creep completely: Boukharov G. et al. [22] analyzed three typical process in rheological curves of brittle rock, and the stress and strain trigger hook body trigger nonlinear sticky pot of components in series is used to describe unsteady creep process. Sterpi and Gioda [23] connected a viscoplastic element in series on the basis of the Kelvin-Voigt model to form a nonlinear elastoviscoplastic creep model. Fahimifar A. et al. [24] aimed at the Sterpi model’s deficiency in predicting steady-state creep under low stress condition; they connected plastic elements that also met Mohr-Coulomb criterion in series on its basis.Stress redistribution caused by underground excavation leads to the emergence of unloading zone in surrounding rock [25, 26]. The mechanical behavior of rock under unloading state is different from that under loading state [27–30]. Unloading often destroys the integrity of the original rock mass, leading to the fragmentation of some rock masses and providing access to groundwater and weathering stress. For example, Kimberley et al. [31] observed that cracks appeared and gradually expanded during unloading of the sample. After further study, it is concluded that the residual stress of inelastic strain gradually accumulates during loading process. If the load is removed slowly, the residual stress is released by creep, but if the load is removed quickly, the residual stress is released by crack initiation and propagation. Huang D. et al. [32] conducted the creep test of marble unloaded under high pressure confining pressure and constant axial pressure decomposition and obtained the relationship between graded unloading degree and deviant stress at failure, elastic modulus, fracture angle, and steady creep time. Wang Y. et al. [33] conducted creep test research on muddy siltstone under different stress levels under confining pressure unloading, and the results showed that the rheological process of rock was anisotropic, with obvious axial flow deformation under high confining pressure, and lateral deformation was larger than axial deformation under low confining pressure.In conclusion, a large number of meaningful studies on rock creep have been carried out in the past, but most of these studies are limited to rocks subjected to constant axial and confining pressures. In fact, in the process of roadway excavation, the stress field is redistributed due to excavation unloading, and the stress of surrounding rock generally decreases. At present, most creep tests are based on loading mechanics, and unloading tests mainly focus on loading and unloading cycles, especially on unloading creep in deep soft rock with high stress. In this paper, triaxial creep tests of mudstone at depth of more than 1000 m under the condition of stage unloading lateral stress and constant axial stress are carried out. By studying the axial instantaneous strain, creep strain, and axial strain rate during unloading, the creep characteristics of rock under different axial stresses and confining pressures are analyzed. Based on the experimental analysis, a new creep damage model of mudstone is established, which can accurately describe the whole process of mudstone creep with a unified function. In addition, FLAC3D software (based on the new creep damage model) was used to simulate the surrounding rock deformation within 165 days after excavation and analyze the surrounding rock deformation law.


## 2. Unloading Creep Test of Mudstone

### 2.1. Test Method

Mudstone samples were taken from the horizontal development roadway of -1350 m underground of Huafeng coal mine. Rock samples are processed into cylinders 50 mm in diameter and 100 mm in height, and triaxial creep test was carried out with MTS 815.02 rock mechanics test system (Figure1) under the condition of step unloading lateral stress and constant axial stress.Figure 1
MTS 815.02 rock mechanics test system.The test is divided into three stages: (1) the confining pressure is increased toσ3=25 MPa and the loading rate is 0.2 MPa/s; (2) when the confining pressure reaches 25 MPa and the deformation tends to be stable, maintain σ3=25 MPa and apply the axial stress at a rate of 0.5 MPa/s. After σ1 is applied to the predetermined test value, observe the axial strain until the deformation of sample tends to be stable; (3) keep σ1 unchanged, and unload σ3 by stages according to the unloading path in Figure 2. The unloading capacity Δσ3 of each stage is 5 MPa, and the unloading rate is 0.05 MPa/s. For each unloading stage, the stress level of this stage shall be maintained unchanged, and the next unloading stage shall be conducted after the given creep time of confining pressure, until the unloading of confining pressure is completed or the test is stopped as the specimen is damaged. Stress-strain curves of rock specimens in triaxial compression tests are shown in Figure 3.Figure 2
Variety of stress during step by step unloading creep.Figure 3
σ1−ε curves of triaxial compression tests.


### 2.2. Analysis of Test Results

According to the creep test, the axial strain-time curve of rock unloading creep under each constant axial stress condition is obtained, as shown in Figure4. The specimens for the creep test are listed in Table 1.Figure 4
The relationship between axial strain and creep time.Table 1
The grouping of specimens for creep test.
50% of peak strength (σ1=65 MPa)60% of peak strength (σ1=73 MPa)70% of peak strength (σ1=81 MPa)1-12-13-11-22-23-21-32-33-32-4As illustrated in Figure4, strain-time curves are generated by performing creep experiments under the condition of graded unloading, in which the strain increases gradually as the confining pressure is reduced. The rock samples utilized in the test are all from coal mines, where rock stress and mining have a significant impact. There are numerous structural planes and joints; therefore, the results of creep tests vary. To gather more creep information, select test blocks with large creep orders for examination when examining the results. The analysis of samples 1-2, 2-4, and 3-2 will be the emphasis of this article.According to the relationship between axial strain and time, the creep process of mudstone can be divided into three stages: (1) in the first stage,σ1 and σ3 have just been loaded to the predetermined stress. At this time, the strain of the sample changes rapidly with time. The ε1 increases continuously, while the axial strain rate gradually decreases to a certain value. This stage is called as the “primary stage.” (2) As the strain rate tends to be constant, the specimen enters secondary stage creep, and the strain-time curve in this stage is close to a straight line with a long duration. (3) At the end of secondary stage creep, the slope of the strain-time curve increases obviously, and the strain rate increases continuously which does not tend to a certain value. The fracture of the sample develops rapidly, even making the material lose stability in a way similar to brittle failure. This failure phenomenon is the result of the strength weakening or unstable deformation of rock mass during creep deformation. It only exists near the failure point of the test, which is called tertiary stage creep. It should be noted that tertiary stage creep does not always occur. In fact, there is almost no tertiary stage creep occurring in hard rock [34–36]. It is still likely to lose stability and destroy directly during unloading even in soft rock. The complete creep curve consists of three stages, as shown in Figure 5.Figure 5
Creep deformation diagram.In addition, the axial strain rate for each axial load condition is determined based on the test data as shown in Figure6 (the initial test section is the axial loading section, and the axial strain rate is quite big, which has no bearing on the study in this paper). The initial section’s axial strain rate is excluded in this article for convenience of analysis.Figure 6
The relationship between axial creep rate and creep time.The changes in creep orders under different confining pressures may be noted in Figures4–6. Creep curves only comprise the primary and secondary stages, in addition to the last level load. In the last order, there are usually three stages of creep. Under certain extreme circumstances, the rock mass may lose stability during the unloading process, causing the fracture to penetrate quickly and lose bearing capacity. The curve indicates that it has progressed to the tertiary stage of creep.The axial strain under each stage of confining pressure is also different. The summary of creep test results is shown in Table2. As the axial load is constant, the instantaneous strain ε0, primary stage creep strain ε1, and secondary stage creep strain ε2 increase step by step along with the unloading behavior of confining pressure, which shows that high confining pressure has a significant restraining effect on rock creep. When the confining pressure decreases, the bearing capacity of rock mass decreases, the fracture development accelerates, and the deformation increases. When the confining pressure drops to a certain value, axial load σ1 is higher than rock mass ultimate bearing capacity under current confining pressure, and then, rock mass enters tertiary stage creep, manifested as large displacement and instability damage in a short time. According to analyzing the proportion of creep deformation in the total deformation, the proportion of creep deformation increases with the decrease of confining pressure. At the last level of load, creep deformation accounts for more than 80% of the total deformation, and creep has become the main factor of rock deformation.Table 2
Creep deformation statistics.
No.σ1 (MPa)σ3 (MPa)ε0/10−4ε1/10−4ε3/10−4ε4/10−4ε1+ε2+ε3ε0ε1+ε2+ε3ε1-265251.780.370.15/0.290.23202.210.780.35/0.510.34152.520.980.57/0.620.38103.031.220.68/0.630.3954.04//17.254.270.812-473252.250.380.32/0.310.24202.690.490.68/0.430.30152.871.561.28/0.990.50103.414.472.1918.477.370.883-281252.460.650.68/0.540.35202.870.911.26/0.760.43153.331.692.8/1.350.57103.672.043.5524.288.140.89In this paper, the effect of axial stress on creep is also studied. By observing Table2 and Figures 3–6, under the same confining pressure, a sample with a larger axial stress tends to have a larger creep deformation, a larger deformation rate, and a larger creep deformation proportion in the total deformation. In addition, it is not difficult to find that the average creep time of σ1=81 MPa sample is 5.84 h, which is far less than the 6.2 h of σ1=73 MPa sample and the 7.89 h of σ1=65 MPa sample. From the point of view of the pressure relief order, σ1=81 MPa sample is generally in the third order; that is, when σ3=15 MPa, it has lost its stability, while σ1=65 MPa sample often reaches the end of σ3=10 MPa stage. Under the condition of unloading confining pressure, the high axial stress can accelerate rock creep deformation significantly. This is because the sample with high axial stress can always reach the axial ultimate bearing strength of its current confining pressure earlier. Under the effect of creep damage and microcrack development, the final instability failure of the sample occurred naturally.To summarize, high axial stress and low confining pressure can expedite creep behavior, resulting in increased creep deformation and a shorter overall creep time. The two have different outward shapes, but the mechanism of action is the same for samples. Samples will be brought closer to the ultimate bearing state of their current stress level as axial stress and confining pressure increase. That is, the creep behavior of sample is determined by its proximity to the ultimate bearing condition under the current stress level.


## 2.1. Test Method

Mudstone samples were taken from the horizontal development roadway of -1350 m underground of Huafeng coal mine. Rock samples are processed into cylinders 50 mm in diameter and 100 mm in height, and triaxial creep test was carried out with MTS 815.02 rock mechanics test system (Figure1) under the condition of step unloading lateral stress and constant axial stress.Figure 1
MTS 815.02 rock mechanics test system.The test is divided into three stages: (1) the confining pressure is increased toσ3=25 MPa and the loading rate is 0.2 MPa/s; (2) when the confining pressure reaches 25 MPa and the deformation tends to be stable, maintain σ3=25 MPa and apply the axial stress at a rate of 0.5 MPa/s. After σ1 is applied to the predetermined test value, observe the axial strain until the deformation of sample tends to be stable; (3) keep σ1 unchanged, and unload σ3 by stages according to the unloading path in Figure 2. The unloading capacity Δσ3 of each stage is 5 MPa, and the unloading rate is 0.05 MPa/s. For each unloading stage, the stress level of this stage shall be maintained unchanged, and the next unloading stage shall be conducted after the given creep time of confining pressure, until the unloading of confining pressure is completed or the test is stopped as the specimen is damaged. Stress-strain curves of rock specimens in triaxial compression tests are shown in Figure 3.Figure 2
Variety of stress during step by step unloading creep.Figure 3
σ1−ε curves of triaxial compression tests.


## 2.2. Analysis of Test Results

According to the creep test, the axial strain-time curve of rock unloading creep under each constant axial stress condition is obtained, as shown in Figure4. The specimens for the creep test are listed in Table 1.Figure 4
The relationship between axial strain and creep time.Table 1
The grouping of specimens for creep test.
50% of peak strength (σ1=65 MPa)60% of peak strength (σ1=73 MPa)70% of peak strength (σ1=81 MPa)1-12-13-11-22-23-21-32-33-32-4As illustrated in Figure4, strain-time curves are generated by performing creep experiments under the condition of graded unloading, in which the strain increases gradually as the confining pressure is reduced. The rock samples utilized in the test are all from coal mines, where rock stress and mining have a significant impact. There are numerous structural planes and joints; therefore, the results of creep tests vary. To gather more creep information, select test blocks with large creep orders for examination when examining the results. The analysis of samples 1-2, 2-4, and 3-2 will be the emphasis of this article.According to the relationship between axial strain and time, the creep process of mudstone can be divided into three stages: (1) in the first stage,σ1 and σ3 have just been loaded to the predetermined stress. At this time, the strain of the sample changes rapidly with time. The ε1 increases continuously, while the axial strain rate gradually decreases to a certain value. This stage is called as the “primary stage.” (2) As the strain rate tends to be constant, the specimen enters secondary stage creep, and the strain-time curve in this stage is close to a straight line with a long duration. (3) At the end of secondary stage creep, the slope of the strain-time curve increases obviously, and the strain rate increases continuously which does not tend to a certain value. The fracture of the sample develops rapidly, even making the material lose stability in a way similar to brittle failure. This failure phenomenon is the result of the strength weakening or unstable deformation of rock mass during creep deformation. It only exists near the failure point of the test, which is called tertiary stage creep. It should be noted that tertiary stage creep does not always occur. In fact, there is almost no tertiary stage creep occurring in hard rock [34–36]. It is still likely to lose stability and destroy directly during unloading even in soft rock. The complete creep curve consists of three stages, as shown in Figure 5.Figure 5
Creep deformation diagram.In addition, the axial strain rate for each axial load condition is determined based on the test data as shown in Figure6 (the initial test section is the axial loading section, and the axial strain rate is quite big, which has no bearing on the study in this paper). The initial section’s axial strain rate is excluded in this article for convenience of analysis.Figure 6
The relationship between axial creep rate and creep time.The changes in creep orders under different confining pressures may be noted in Figures4–6. Creep curves only comprise the primary and secondary stages, in addition to the last level load. In the last order, there are usually three stages of creep. Under certain extreme circumstances, the rock mass may lose stability during the unloading process, causing the fracture to penetrate quickly and lose bearing capacity. The curve indicates that it has progressed to the tertiary stage of creep.The axial strain under each stage of confining pressure is also different. The summary of creep test results is shown in Table2. As the axial load is constant, the instantaneous strain ε0, primary stage creep strain ε1, and secondary stage creep strain ε2 increase step by step along with the unloading behavior of confining pressure, which shows that high confining pressure has a significant restraining effect on rock creep. When the confining pressure decreases, the bearing capacity of rock mass decreases, the fracture development accelerates, and the deformation increases. When the confining pressure drops to a certain value, axial load σ1 is higher than rock mass ultimate bearing capacity under current confining pressure, and then, rock mass enters tertiary stage creep, manifested as large displacement and instability damage in a short time. According to analyzing the proportion of creep deformation in the total deformation, the proportion of creep deformation increases with the decrease of confining pressure. At the last level of load, creep deformation accounts for more than 80% of the total deformation, and creep has become the main factor of rock deformation.Table 2
Creep deformation statistics.
No.σ1 (MPa)σ3 (MPa)ε0/10−4ε1/10−4ε3/10−4ε4/10−4ε1+ε2+ε3ε0ε1+ε2+ε3ε1-265251.780.370.15/0.290.23202.210.780.35/0.510.34152.520.980.57/0.620.38103.031.220.68/0.630.3954.04//17.254.270.812-473252.250.380.32/0.310.24202.690.490.68/0.430.30152.871.561.28/0.990.50103.414.472.1918.477.370.883-281252.460.650.68/0.540.35202.870.911.26/0.760.43153.331.692.8/1.350.57103.672.043.5524.288.140.89In this paper, the effect of axial stress on creep is also studied. By observing Table2 and Figures 3–6, under the same confining pressure, a sample with a larger axial stress tends to have a larger creep deformation, a larger deformation rate, and a larger creep deformation proportion in the total deformation. In addition, it is not difficult to find that the average creep time of σ1=81 MPa sample is 5.84 h, which is far less than the 6.2 h of σ1=73 MPa sample and the 7.89 h of σ1=65 MPa sample. From the point of view of the pressure relief order, σ1=81 MPa sample is generally in the third order; that is, when σ3=15 MPa, it has lost its stability, while σ1=65 MPa sample often reaches the end of σ3=10 MPa stage. Under the condition of unloading confining pressure, the high axial stress can accelerate rock creep deformation significantly. This is because the sample with high axial stress can always reach the axial ultimate bearing strength of its current confining pressure earlier. Under the effect of creep damage and microcrack development, the final instability failure of the sample occurred naturally.To summarize, high axial stress and low confining pressure can expedite creep behavior, resulting in increased creep deformation and a shorter overall creep time. The two have different outward shapes, but the mechanism of action is the same for samples. Samples will be brought closer to the ultimate bearing state of their current stress level as axial stress and confining pressure increase. That is, the creep behavior of sample is determined by its proximity to the ultimate bearing condition under the current stress level.


## 3. Constitutive Model

### 3.1. Burgers Creep Constitutive Model

In general, creep curves usually include primary stage and secondary stage, with almost no tertiary stage creep. The Burgers model is generally used to describe the creep behavior of rock [16, 37]. The one-dimensional creep constitutive equation of the Burgers model is as follows:
(1)ε=1E1+1η1t+1E21−e−E2/η2tσ,where ε is the axial strain; E1 is the instantaneous elastic modulus; t is the creep time; E2 is the viscoelastic modulus; η1 and η2 are all viscosity coefficients; and σ is the axial stress.The three-dimensional creep constitutive equation of the Burgers model is as follows:(2)eij=Sij2G1+Sij2η1t+Sij2G21−e−G2/η2t,εm=σm3K,where eij represents deviatoric strain, Sij represents deviatoric stress, K represents bulk modulus, G1 and G2 represent shear modulus, εm represents spherical strain, and σm represents spherical stress. K, G1, G2, η1, and η2 are all model parameters, which are determined by fitting test results, where K, G1, and G2 are defined as
(3)K=E131−2μ,G1=E121+μ,G2=E221+μ.In which,μ is Poisson’s ratio, which is the physical characteristic parameter of a rock sample. Therefore, five model parameters, such as K, G1, G2, η1, and η2, can be simplified into four parameters, i.e., E1, E2, η1, and η2. As a result, for the Burgers model, whether in one-dimensional or three-dimensional conditions, the solution of equation parameters can always be simplified to E1, E2, η1, and η2.


### 3.2. An Improved Nonlinear Viscoelastoplastic Rheological Model Based on the Burgers Model

According to the previous analysis, the creep curve of mudstone can be divided into three stages (primary, secondary, and tertiary), and the instantaneous strain will appear during the unloading of the confining pressure, so the axial strain valueε of the whole creep process can be expressed as
(4)ε=ε0+ε1t+ε2t+ε3t,where ε denotes the total axial strain of the mudstone, ε0 denotes the instantaneous strain, ε1t denotes the primary stage strain, ε2t denotes the secondary stage strain, and ε3t denotes the tertiary stage strain.The Burgers model has a good fitting effect on the primary and secondary stage creep. However, it cannot simulate the whole creep process completely because of the lock of tertiary stage simulation components. In order to fully simulate the three stages of creep, considering the influence of time on creep deformation, an improved nonlinear viscoelastoplastic rheological model is constructed by connecting nonlinear viscoplastic bodies (NVPB) in series on the basis of burger bodies [12, 38, 39].As shown in Figure7, NVPB is composed of a switch and a nonlinear dashpot in parallel. When σ0≤σs, the switch is connected. At this time, the nonlinear dashpot is short circuited and does not participate in the deformation of the model, and the model is equivalent to the Burgers model, which can realize the simulation of instantaneous strain and the first two stages of creep. When σ0>σs, the constitutive equation of NVPB is as follows:
(5)εt=Hσ0−σSηtnt0n−1=Hσ0−σSηtn,where n is creep index, which reflects the speed of accelerated creep of rock, as shown in Figure 8. t0 is reference time, which is taken as 1; σs is yield strength or long-term strength [40]. σ0 is deviatoric stress, and H is shown in the following formula:
(6)Hσ0−σS=0,σ0≤σS,σ0−σS,σ0>σS.Figure 7
Nonlinear viscoplastic body (NVPB).Figure 8
Creep curve of NVPB model.From the creep characteristic curve, it can be seen that the strain rate of rock increases gradually in the tertiary stage, so takingn>1. At this time, NVPB connects with the Burgers model and begins to participate in the deformation of the model, as shown in Figure 9. The constitutive equation is
(7)ε=1E1+1η1t+1E21−e−E2/η2tσ+Hσ−σSη3tn,where ε is the axial strain, E1 is the instantaneous elastic modulus, t is the creep time, E2 is the viscoelastic modulus, η1, η2, and η3 are all viscosity coefficients, σ is the axial stress, and σs is the yield strength or long-term strength.Figure 9
An improved nonlinear viscoelastoplastic rheological model based on burgers model.Based on the experimental data and the results of constitutive derivation, the parameters of the creep constitutive model are solved by the least square method. The creep parameters of the Burger model under different stress conditions are shown in Tables3–5.Table 3
Parameters of improved nonlinear viscoelastoplastic rheological model based on burgers model under 50% axial ultimate load (the creep duration under 5 MPa confining pressure is short, and there is too little data to be used for formula fitting, so it is omitted).
σ3 (MPa)E1 (GPa)η1 (GPa·h)E2 (GPa)η2 (GPa·h)R2256967.831669.64222.9018.670.98205737.031363.63203.5216.250.97154359.091267.78185.1014.470.98103825.261112.59178.6312.090.97Table 4
Parameters of improved nonlinear viscoelastoplastic rheological model based on burgers model under 60% axial ultimate load.
σ3 (MPa)E1 (GPa)η1 (GPa·h)E2 (GPa)η2 (GPa·h)R2255694.571343.39177.4216.880.95205292.271003.04154.0815.040.96153898.54762.28136.5413.070.94103014.65701.13132.4611.460.98Table 5
Parameters of improved nonlinear viscoelastoplastic rheological model based on burgers model under 70% axial ultimate load.
σ3 (MPa)E1 (GPa)η1 (GPa·h)E2 (GPa)η2 (GPa·h)R2254945.23906.10108.9210.240.93204076.14824.6895.0211.120.94153688.78547.6378.008.480.97It can be seen from Tables3–5 that the elastic modulus increases significantly with the decrease of confining pressure, and the change trend is positively related to the change of axial stress. The elastic modulus E represents the elastic deformation capacity of the model, and the larger the elastic modulus is, the smaller the elastic deformation is. In this model, the amplitude of variation of E1 under each axial stress decreases with the decrease of confining pressure, and the variation amplitude of E1 under the low axial stress is larger than that under the high axial stress. It shows that the instantaneous strain increases with the increase of confining pressure, which is consistent with the previous observation. Under the condition of unloading confining pressure, the viscosity coefficient has a similar trend to the elastic modulus. The parameters of the viscosity coefficient represent the rate of model creep deformation, and the larger the value is, the smaller the rate is. It can be seen from Tables 3–5 that the viscosity coefficient has a similar trend to elastic modulus under the condition of unloading confining pressure, which indicates that the reduction of confining pressure under the condition of constant axial force will not only increase the instantaneous strain during unloading but also have the same effect on the creep value of the test block. With the unloading process, the deformation of the test block accumulates and the cracks develop gradually. The samples reach their ultimate bearing state and lose stability finally.After about 1 hour of secondary stage creep under 10 MPa confining pressure, sample 3-2’s axial strain rate and creep displacement increase consistently, and sample 3-2 loses its stability quickly in a short period, exhibiting distinct tertiary stage creep features. The least square approach is utilized to calculate the tertiary stage creep parameters of mudstone obtained by entering the time-strain data received from the test into the constitutive equation of the enhanced nonlinear viscoelastoplastic rheological model based on the burgers model (in Table6).Table 6
Parameters of tertiary stage creep of the improved nonlinear viscoelastoplastic rheological model based on the burgers model.
σ3 (MPa)E1 (GPa)η1 (GPa·h)E2 (GPa)η2 (GPa·h)η3 (GPa·h)nR2102308.575216.35229.135.853.72e+1043.740.99It can be seen from Figure10 that before the sample enters the tertiary stage, the crack development has been accelerated and the deformation increases rapidly. This process makes the damage to the sample accumulate continuously, until finally it loses the bearing capacity and leads to instability and failure. The improved nonlinear viscoelastoplastic rheological model based on the burgers model can accurately simulate the time-displacement relationship of mudstone in tertiary stage creep.Figure 10
Comparison between calculated displacement and actual displacement.


## 3.1. Burgers Creep Constitutive Model

In general, creep curves usually include primary stage and secondary stage, with almost no tertiary stage creep. The Burgers model is generally used to describe the creep behavior of rock [16, 37]. The one-dimensional creep constitutive equation of the Burgers model is as follows:
(1)ε=1E1+1η1t+1E21−e−E2/η2tσ,where ε is the axial strain; E1 is the instantaneous elastic modulus; t is the creep time; E2 is the viscoelastic modulus; η1 and η2 are all viscosity coefficients; and σ is the axial stress.The three-dimensional creep constitutive equation of the Burgers model is as follows:(2)eij=Sij2G1+Sij2η1t+Sij2G21−e−G2/η2t,εm=σm3K,where eij represents deviatoric strain, Sij represents deviatoric stress, K represents bulk modulus, G1 and G2 represent shear modulus, εm represents spherical strain, and σm represents spherical stress. K, G1, G2, η1, and η2 are all model parameters, which are determined by fitting test results, where K, G1, and G2 are defined as
(3)K=E131−2μ,G1=E121+μ,G2=E221+μ.In which,μ is Poisson’s ratio, which is the physical characteristic parameter of a rock sample. Therefore, five model parameters, such as K, G1, G2, η1, and η2, can be simplified into four parameters, i.e., E1, E2, η1, and η2. As a result, for the Burgers model, whether in one-dimensional or three-dimensional conditions, the solution of equation parameters can always be simplified to E1, E2, η1, and η2.


## 3.2. An Improved Nonlinear Viscoelastoplastic Rheological Model Based on the Burgers Model

According to the previous analysis, the creep curve of mudstone can be divided into three stages (primary, secondary, and tertiary), and the instantaneous strain will appear during the unloading of the confining pressure, so the axial strain valueε of the whole creep process can be expressed as
(4)ε=ε0+ε1t+ε2t+ε3t,where ε denotes the total axial strain of the mudstone, ε0 denotes the instantaneous strain, ε1t denotes the primary stage strain, ε2t denotes the secondary stage strain, and ε3t denotes the tertiary stage strain.The Burgers model has a good fitting effect on the primary and secondary stage creep. However, it cannot simulate the whole creep process completely because of the lock of tertiary stage simulation components. In order to fully simulate the three stages of creep, considering the influence of time on creep deformation, an improved nonlinear viscoelastoplastic rheological model is constructed by connecting nonlinear viscoplastic bodies (NVPB) in series on the basis of burger bodies [12, 38, 39].As shown in Figure7, NVPB is composed of a switch and a nonlinear dashpot in parallel. When σ0≤σs, the switch is connected. At this time, the nonlinear dashpot is short circuited and does not participate in the deformation of the model, and the model is equivalent to the Burgers model, which can realize the simulation of instantaneous strain and the first two stages of creep. When σ0>σs, the constitutive equation of NVPB is as follows:
(5)εt=Hσ0−σSηtnt0n−1=Hσ0−σSηtn,where n is creep index, which reflects the speed of accelerated creep of rock, as shown in Figure 8. t0 is reference time, which is taken as 1; σs is yield strength or long-term strength [40]. σ0 is deviatoric stress, and H is shown in the following formula:
(6)Hσ0−σS=0,σ0≤σS,σ0−σS,σ0>σS.Figure 7
Nonlinear viscoplastic body (NVPB).Figure 8
Creep curve of NVPB model.From the creep characteristic curve, it can be seen that the strain rate of rock increases gradually in the tertiary stage, so takingn>1. At this time, NVPB connects with the Burgers model and begins to participate in the deformation of the model, as shown in Figure 9. The constitutive equation is
(7)ε=1E1+1η1t+1E21−e−E2/η2tσ+Hσ−σSη3tn,where ε is the axial strain, E1 is the instantaneous elastic modulus, t is the creep time, E2 is the viscoelastic modulus, η1, η2, and η3 are all viscosity coefficients, σ is the axial stress, and σs is the yield strength or long-term strength.Figure 9
An improved nonlinear viscoelastoplastic rheological model based on burgers model.Based on the experimental data and the results of constitutive derivation, the parameters of the creep constitutive model are solved by the least square method. The creep parameters of the Burger model under different stress conditions are shown in Tables3–5.Table 3
Parameters of improved nonlinear viscoelastoplastic rheological model based on burgers model under 50% axial ultimate load (the creep duration under 5 MPa confining pressure is short, and there is too little data to be used for formula fitting, so it is omitted).
σ3 (MPa)E1 (GPa)η1 (GPa·h)E2 (GPa)η2 (GPa·h)R2256967.831669.64222.9018.670.98205737.031363.63203.5216.250.97154359.091267.78185.1014.470.98103825.261112.59178.6312.090.97Table 4
Parameters of improved nonlinear viscoelastoplastic rheological model based on burgers model under 60% axial ultimate load.
σ3 (MPa)E1 (GPa)η1 (GPa·h)E2 (GPa)η2 (GPa·h)R2255694.571343.39177.4216.880.95205292.271003.04154.0815.040.96153898.54762.28136.5413.070.94103014.65701.13132.4611.460.98Table 5
Parameters of improved nonlinear viscoelastoplastic rheological model based on burgers model under 70% axial ultimate load.
σ3 (MPa)E1 (GPa)η1 (GPa·h)E2 (GPa)η2 (GPa·h)R2254945.23906.10108.9210.240.93204076.14824.6895.0211.120.94153688.78547.6378.008.480.97It can be seen from Tables3–5 that the elastic modulus increases significantly with the decrease of confining pressure, and the change trend is positively related to the change of axial stress. The elastic modulus E represents the elastic deformation capacity of the model, and the larger the elastic modulus is, the smaller the elastic deformation is. In this model, the amplitude of variation of E1 under each axial stress decreases with the decrease of confining pressure, and the variation amplitude of E1 under the low axial stress is larger than that under the high axial stress. It shows that the instantaneous strain increases with the increase of confining pressure, which is consistent with the previous observation. Under the condition of unloading confining pressure, the viscosity coefficient has a similar trend to the elastic modulus. The parameters of the viscosity coefficient represent the rate of model creep deformation, and the larger the value is, the smaller the rate is. It can be seen from Tables 3–5 that the viscosity coefficient has a similar trend to elastic modulus under the condition of unloading confining pressure, which indicates that the reduction of confining pressure under the condition of constant axial force will not only increase the instantaneous strain during unloading but also have the same effect on the creep value of the test block. With the unloading process, the deformation of the test block accumulates and the cracks develop gradually. The samples reach their ultimate bearing state and lose stability finally.After about 1 hour of secondary stage creep under 10 MPa confining pressure, sample 3-2’s axial strain rate and creep displacement increase consistently, and sample 3-2 loses its stability quickly in a short period, exhibiting distinct tertiary stage creep features. The least square approach is utilized to calculate the tertiary stage creep parameters of mudstone obtained by entering the time-strain data received from the test into the constitutive equation of the enhanced nonlinear viscoelastoplastic rheological model based on the burgers model (in Table6).Table 6
Parameters of tertiary stage creep of the improved nonlinear viscoelastoplastic rheological model based on the burgers model.
σ3 (MPa)E1 (GPa)η1 (GPa·h)E2 (GPa)η2 (GPa·h)η3 (GPa·h)nR2102308.575216.35229.135.853.72e+1043.740.99It can be seen from Figure10 that before the sample enters the tertiary stage, the crack development has been accelerated and the deformation increases rapidly. This process makes the damage to the sample accumulate continuously, until finally it loses the bearing capacity and leads to instability and failure. The improved nonlinear viscoelastoplastic rheological model based on the burgers model can accurately simulate the time-displacement relationship of mudstone in tertiary stage creep.Figure 10
Comparison between calculated displacement and actual displacement.


## 4. Numerical Simulation

At present, no accurate theoretical formula has been established for the influence of rock rheology on deep underground engineering. As an analysis method, computer numerical simulation is widely used in practical engineering construction and also provides an important reference value for solving the problem of surrounding rock stability under the condition of in situ rock stress [41–45]. To study the influence of creep on geological engineering and explore the law of the roadway deformation, FLAC3D software is used in this section to simulate the deformation behavior of the roadway surrounding rock over a long time, obtain the displacement evolution of roadway surrounding rock, and summarize the law of roadway deformation.


### 4.1. Establish Numerical Simulation Model

The displacement and stress constraint boundary conditions were chosen to fix the displacement of the model’s edge and bottom, and a uniform load was given to the top of the model to replicate in situ stress. The road will be 1350 meters long, 3.6 meters wide in cross-section, 1.5 meters high on the side walls, and 1.8 meters high in the arch. The model’s dimensions are50m×50m×60 m. The parameters utilized in the entire model are based on real-world engineering characteristics that approximate nearby rock displacement 165 days after the motorway was excavated. Figure 11 depicts the model diagram.Figure 11
Numerical simulation model.


### 4.2. Numerical Simulation Results

The FLAC3D software is used to simulate the whole period of 165 days after the excavation of permanent main roadway of Huafeng mine. The vertical and horizontal displacement contours of surrounding rock 1 d, 30 d, 60 d, 90 d, 120 d, and 165 d after the excavation is shown in Figures 12 and 13.Figure 12
Horizontal displacement contours.Figure 13
Vertical displacement contours.


### 4.3. Roadway Deformation Analysis

The deformation of the roadway can be immediately determined by evaluating the displacement of adjacent rock around the route. Figure14 depicts the data of the roadway’s deformation following the excavation. The highest displacement is the change of roof position, which is 1129.50 mm; the second largest displacement is the approach of the two sides, which is 884.75 mm; and the lowest displacement is the change of floor position, which is 743.23 mm. The highest horizontal displacement occurs from the middle of the side wall to the arch foot, and the horizontal displacement is primarily represented in the internal displacement of the side wall. In addition, the maximum vertical displacement is located in the center of the roof, the position with the largest floor heave is in the middle of the floor, and the vertical displacement of the side wall is very small.Figure 14
Statistics of roadway deformation.The roadway would distort fast without support following excavation due to high in situ stress, especially in the vertical direction. The roof lowers quickly, as does the abscission layer, and the floor heave is severe. The inward shrinkage of the side wall is the principal manifestation of the horizontal displacement, and the approach amount of the higher section of the side wall is much larger than the lower half. In general, the deformation of the roadway is mostly manifested in the deformation of the roadway’s upper structure. The surrounding rock is influenced by the simultaneous forces of its own gravity and in situ stress when the roadway is excavated, generating a caving arch, and the upper structure would quickly lose its stability.In addition, according to the numerical simulation of roadway deformation, the deformation of surrounding rock around the roadway presents significant rheological characteristics. At the beginning of excavation, the roadway was severely deformed. The deformation rates of the roof, side walls, and floor in 90 d were 7.83 mm/d, 5.96 mm/d, and 4.89 mm/d, respectively. At this time, the displacement is mainly caused by the reduction of in situ stress. After the initial rapid deformation, the redistribution of in situ stress has stabilized at this time, and the deformation rate of the roadway has gradually slowed down. The displacement rates of the roof, side wall, and floor in 90~165 days are 5.66 mm/d, 4.65 mm/d, and 4.04 mm/d, respectively. It is not difficult to know that in the late stage of roadway deformation, although the deformation rate tends to decrease, the deformation of roadway is far from convergence, and it still has a large deformation rate after 165 days of excavation. In fact, roadway often takes months or even years to enter into stable deformation after undergoing severe initial deformation. As for the convergence of the deformation of the roadway, it depends on the support and the distribution of in situ stress. In fact, many roadways will continue to deform after many years, which will pose great safety risks to production.


## 4.1. Establish Numerical Simulation Model

The displacement and stress constraint boundary conditions were chosen to fix the displacement of the model’s edge and bottom, and a uniform load was given to the top of the model to replicate in situ stress. The road will be 1350 meters long, 3.6 meters wide in cross-section, 1.5 meters high on the side walls, and 1.8 meters high in the arch. The model’s dimensions are50m×50m×60 m. The parameters utilized in the entire model are based on real-world engineering characteristics that approximate nearby rock displacement 165 days after the motorway was excavated. Figure 11 depicts the model diagram.Figure 11
Numerical simulation model.


## 4.2. Numerical Simulation Results

The FLAC3D software is used to simulate the whole period of 165 days after the excavation of permanent main roadway of Huafeng mine. The vertical and horizontal displacement contours of surrounding rock 1 d, 30 d, 60 d, 90 d, 120 d, and 165 d after the excavation is shown in Figures 12 and 13.Figure 12
Horizontal displacement contours.Figure 13
Vertical displacement contours.


## 4.3. Roadway Deformation Analysis

The deformation of the roadway can be immediately determined by evaluating the displacement of adjacent rock around the route. Figure14 depicts the data of the roadway’s deformation following the excavation. The highest displacement is the change of roof position, which is 1129.50 mm; the second largest displacement is the approach of the two sides, which is 884.75 mm; and the lowest displacement is the change of floor position, which is 743.23 mm. The highest horizontal displacement occurs from the middle of the side wall to the arch foot, and the horizontal displacement is primarily represented in the internal displacement of the side wall. In addition, the maximum vertical displacement is located in the center of the roof, the position with the largest floor heave is in the middle of the floor, and the vertical displacement of the side wall is very small.Figure 14
Statistics of roadway deformation.The roadway would distort fast without support following excavation due to high in situ stress, especially in the vertical direction. The roof lowers quickly, as does the abscission layer, and the floor heave is severe. The inward shrinkage of the side wall is the principal manifestation of the horizontal displacement, and the approach amount of the higher section of the side wall is much larger than the lower half. In general, the deformation of the roadway is mostly manifested in the deformation of the roadway’s upper structure. The surrounding rock is influenced by the simultaneous forces of its own gravity and in situ stress when the roadway is excavated, generating a caving arch, and the upper structure would quickly lose its stability.In addition, according to the numerical simulation of roadway deformation, the deformation of surrounding rock around the roadway presents significant rheological characteristics. At the beginning of excavation, the roadway was severely deformed. The deformation rates of the roof, side walls, and floor in 90 d were 7.83 mm/d, 5.96 mm/d, and 4.89 mm/d, respectively. At this time, the displacement is mainly caused by the reduction of in situ stress. After the initial rapid deformation, the redistribution of in situ stress has stabilized at this time, and the deformation rate of the roadway has gradually slowed down. The displacement rates of the roof, side wall, and floor in 90~165 days are 5.66 mm/d, 4.65 mm/d, and 4.04 mm/d, respectively. It is not difficult to know that in the late stage of roadway deformation, although the deformation rate tends to decrease, the deformation of roadway is far from convergence, and it still has a large deformation rate after 165 days of excavation. In fact, roadway often takes months or even years to enter into stable deformation after undergoing severe initial deformation. As for the convergence of the deformation of the roadway, it depends on the support and the distribution of in situ stress. In fact, many roadways will continue to deform after many years, which will pose great safety risks to production.


## 5. Conclusion

To study the influence of creep of seep soft rock on the deformation behaviors of surrounding rock under unloading conditions, experimental tests, constitutive models, and numerical simulation were combined to illuminate the deformation laws. The following conclusions can be drawn:(1)
Under the same axial load, the deformation of mudstone increases with the decrease of confining pressure, which is not only reflected in the increase of instantaneous strain caused by the reduction of confining pressure but also in the creep deformation and creep rate of mudstone. Under the same confining pressure, the greater the axial load is, the more severe the deformation of the specimen is (specifically, the greater the instantaneous strain and creep rate are), and the less the unloading stage is, the earlier the instability and failure of the specimen are. In addition, the severity of the creep behavior of samples is closely related to the current stress level. The closer the external load is to the ultimate bearing capacity under the current stress level, the stronger the creep behavior is(2)
The Burgers model and creep curve are analyzed. NVPB is introduced to simulate the tertiary stage creep. The influence of time on creep parameters is considered. The improved model is simple in structure and has fewer calculation parameters. It can simulate the whole creep process well, and the calculated value has a high correlation with the test value. The creep behavior of rock under other stress conditions can be simulated only by modifying creep parameters. Based on the above constitutive model, FLAC3D software was used to simulate the deformation of the unsupported roadway after excavation. During the simulation period of 165 d, the deformation of the roadway shows distinct timeliness. The deformation law of the roadway was summarized, which has a certain guiding effect on the actual production



---

*Source: 7175774-2022-05-04.xml*