Analysis and numerical calculation of a coupled creep and strain-softening model for soft rock tunnels

Proper mechanical model selection is critical in tunnel support design and stability analysis, especially to reflect the creep and strain-softening behavior of soft rock. We present a coupled nonlinear Burgers strain-softening (NBSS) model and numerical calculation method to investigate the coupled effects of creep and strain-softening of soft rock tunnels. The nonlinear elastic-viscous model is used to simulate the steady creep behavior of mudstone, and the nonlinear viscoplastic strain-softening model is used to simulate the accelerated creep behavior and post-peak strength attenuation behavior. The experimental results show that the viscoplastic parameters and post-peak softening parameters of mudstone are highly sensitive to confining pressure and exhibit nonlinear characteristics. The accelerated creep curve obtained by the numerical calculation is consistent with the experiments, which verifies the model reliability. We use the NBSS and nonlinear Burgers Mohr-Coulomb (NBMC) models to calculate the plastic zone distribution characteristics and deformation law. The distribution of the plastic zone calculated by the NBSS model is larger with more localized fractures. The NBSS model is useful for studying the evolution of stress and displacement fields of complex surrounding rock mass, which provides important theoretical guidelines for the support design and stability analysis of soft rock tunnel engineering.


Introduction
Rock strength attenuation and creep characteristics control the stability of supporting structures during tunnel excavation and support engineering [1][2][3][4][5]. The material strength characteristics of surrounding rock strength gradually decreases with time, particularly in soft rock settings [6,7]. Soft rock creep prior to reaching the peak value typically undergoes decay and steady creep. Under loading, the micro-cracks in soft rock expand over time and the rock strength gradually decreases with increasing viscoplastic strain. When the strength is reduced to the load stress level, the soft rock will enter the accelerated creep stage. The strength attenuation process of soft rock over time can be characterized by creep damage and strain-softening models [8][9][10][11]. Previous rock creep damage studies focused mainly on how to determine creep damage variables and applied mechanical models to describe accelerated rock creep. Because the rheological process of surrounding rock is very complex in practical engineering and creep parameters tend to change with the stress environment, some previous studies developed an unsteady fractional order creep model to investigate the characteristics of accelerated rock creep [12][13][14][15][16]. The unsteady fractional order creep model well reflects the creep acceleration stage but is impractical for model development and calculations and still requires basic elasticviscoplastic theory as a support [17][18][19]. Gioda and Sterpi [20] considered that rock viscoplastic strain is the root cause of rheological instability in tunnel engineering. Other scholars began from rock elasto-viscoplasticity theory and constructed rock creep damage models. For example, Zhao et al. [21] proposed an elasto-viscoplastic rheological creep damage model and data processing algorithm using experiments and data processing methods of rock visco-elastic and visco-plastic strain separation, and quantified the creep model parameters by linear interpolation [22]. Fu et al. [23] built a three-dimensional discrete element grain-based stress corrosion model incorporating the theories of subcritical crack growth and chemical reaction rate to explore the time-dependent behavior of damage evolution and fracture patterns of brittle rocks on a mesoscopic scale. Manica et al. [24] proposed a rheological constitutive model of mudstone developed within an elastic-plastic framework, which completely describes the rheological characteristics of mudstone. Zhang et al. [25] proposed a novel classifier ensemble to improve prediction accuracy for tunnel squeezing problems by aggregating seven individual classifiers using the weighted voting method.
Within the study of post-peak rock strength attenuation characteristics, previous studies have mainly constructed different types of rock strain-softening models. For example, Alejano et al. [26] constructed a mechanical model of post-peak rock strength attenuation based on elastoplastic mechanics theory using a geological strength index of the surrounding rock rating system [27]. Krajcinovic and Silva [28] established a statistical constitutive model of rock strain-softening damage using a random distribution of internal defects and rock microelement strength following a Weibull distribution [29]. Sun et al. [30] determined the relationship between post-peak softening modulus and confining pressure of granite using nonlinear fitting methods and constructed a strain-softening mechanical model of granite according to the geometric relationship of stress-strain space. Feng et al. [31] developed a strain-softening damage model of the rocks with defect growth based on damage evolution. Wang et al. [32] proposed a new method for the numerical calculation of rock strain-softening processes.
Abundant achievements have been made in the study of damage creep and strength attenuation of rock. However, an approach to study the coupled law of creep and strain-softening in soft rock remains poorly defined, as well as how to numerically reflect the accelerated creep characteristics of soft rock and strength attenuation law after accelerated creep.
To address this issue, we analyzed the coupled mechanism of creep and strainsoftening of soft rock surrounding a tunnel using theoretical methods and performed laboratory experiments to study the strain-softening and creep damage characteristics of mudstone. We construct a nonlinear Burgers creep damage model and Mohr-Coulomb strain-softening model based on the nonlinear exponential fitting method, and propose numerical calculation methods using a nonlinear Burgers Mohr-Coulomb (NBMC) model and nonlinear Burgers strain-softening (NBSS) model. We calculate the coupled effect law of creep and strain-softening of rock surrounding a circular tunnel using numerical simulations, and analyze the distribution characteristics of the plastic zone and deformation evolution law of surrounding rock.

Strain-softening behavior of rock mass
A hypothetical two-dimensional circular soft rock tunnel is shown in Fig 1A. The surrounding rock is in a state of hydrostatic pressure, the original rock stress is σ 0 , and the excavation radius is r 0 . Without considering rock creep, the tangential stress σ θ at point A approaches 2σ 0 and the radial stress σ r approaches 0 owing to the unloading of the surrounding rock at time T 1 . At this time, the surrounding rock at point A is within a state similar to uniaxial compression. If σ θ − σ r is greater than the surrounding rock strength at point A, a plastic strain-softening zone gradually forms in the surrounding rock. At T 2 , if σ θ − σ r is greater than the rock strength at point B, the post-peak rock strength will attenuate with increasing cumulative plastic strain. By analogy, after T 3 and T 4 , the plastic zone tends to stabilize after points C and D reach the yield strength and the surrounding rock between points A and B reaches the residual strength. At this time, the surrounding rock of the tunnel is composed of three parts: an elastic zone; a plastic strain-softening zone; and a plastic residual zone. The stress-strain curves of the surrounding rock at points A, B, C, and D can be approximately expressed as a simplified three-line segment, as shown in Fig 1B. The process by which the strength of rock surrounding a tunnel gradually attenuates with post-peak cumulative plastic strain can be characterized by the plastic softening coefficient η, which has a nonlinear relationship with the plastic strain. It is assumed that the stress component expression of the rock mass yield condition considering damage is as follows: When the surrounding rock reaches the peak strength, its strength parameters gradually change with increasing η. The η value is assumed to decrease to η � and the rock strength reaches the residual strength. This relationship can be expressed as [27]: : where ω is the rock strength parameter, ω p is the peak strength parameter, and ω r is the residual strength parameter. The variation law is shown in Fig 1C. The process by which η decreases from 0 to η � is called surrounding rock strain-softening.

Coupled behavior of rock creep and strain-softening
Under rock creep conditions, the unloading of rock surrounding a tunnel concentrate stress concentration at point a at T 1 , as shown in Fig 1A. The concentrated stress σ θ at point a is considered to be less than the peak strength of the surrounding rock. The creep damage deformation of the surrounding rock at point a increases with time owing to the rock creep damage characteristics. The creep trajectory of the surrounding rock at point a gradually extends to the post-peak stress-strain curve of point A. Accelerated creep failure occurs when the strength of the surrounding rock is less than σ θ owing to creep damage. The creep trajectory is shown in Fig 1B and the creep curve of the corresponding surrounding rock is shown in Fig 1D. After creep failure of the surrounding rock at point a, its strength will continue to attenuate. According to the stress-strain curve characteristics at point A, the post-peak strength of the surrounding rock continues to attenuate to the residual strength with increasing cumulative plastic strain. At T 2 , the surrounding rock at point b is affected by σ θ − σ r and its creep trajectory gradually extends to the post-peak stress-strain curve of point B. After accelerated creep failure, the

PLOS ONE
strength of the surrounding rock at point b continues to decrease to the residual strength with increasing cumulative plastic shear strain. After T 3 and T 4 , the σ θ − σ r stress level of the surrounding rock at point d is relatively small. The relationship between the creep trajectory and post-peak stress-strain curve at point D is shown in Fig 1B. Creep deformation of the surrounding rock tends to be stable at this stress level, and creep damage does not accelerate creep failure of the surrounding rock at point d. The creep curve is shown in Fig 1D. At this time, the plastic zone of the surrounding rock is essentially stable. The creep damage variable D is used to characterize the attenuation process of surrounding rock strength with time and is related to the cumulative plastic strain. With increasing cumulative rock viscoplastic strain, damage cracks appear and the strength gradually attenuates. Accelerated creep failure occurs when the cumulative rock plastic strain at point b reaches Δε 1 . The creep damage variable D is assumed to be η D , as shown in Fig 1C. After accelerated rock creep damage, the surrounding rock strength continues to decrease with increasing cumulative plastic strain until reaching the residual strength. As the rock creep damage changes to accelerated creep, η increases from 0 to η D . After accelerated creep, η increases from η D to η � while approaching the residual strength. The transition from pre-peak creep damage to post-peak residual strength is called rock strain-softening.

Nonlinear behavior of rock failure
When discussing the creep damage and strain-softening process in Fig 1A, the basic assumption is that the σ r of each point in the surrounding rock remains constant. While the confining pressure of the surrounding rock remains constant, the σ r of the surrounding rock in the unloading path of the surrounding rock in the same area changes constantly. The unloading path of point D is shown in Fig 1E. When η = 0, the surrounding rock at point D is in the pre-peak state; when η = η � , the surrounding rock at point D is in the residual strength state. After reaching the peak strength, the stress at point D changes from II to III and ultimately reaches the residual strength IV. The variation range of the confining pressure can be expressed by Δσ 3 , as shown in Fig 1B. In practical engineering, the creep damage and strain-softening of tunnel soft rock show complex nonlinear characteristics. A nonlinear mechanical model should therefore be established in accordance with the creep damage and strength attenuation characteristics of soft rock to better analyze its deformation and failure characteristics in tunnel settings and ultimately propose reasonable control methods.

Rock specimen
Specimens were collected from the surrounding rock face of the DK77+684 section of the Milin tunnel on the Sichuan-Tibet railway. The content of clay minerals is more than 50%, the content of carbonate minerals is less than 25%, and the content of felsic minerals is less than 25%. The yellow mudstone is relatively dense, with an average porosity of 3.29% and water content of 15.7%. The surrounding rock was mainly yellow mudstone with distinct rheological properties. To study the deformation mechanism, we performed conventional triaxial compression tests to determine the rock strain-softening mechanical properties and triaxial creep tests to investigate the rock creep characteristics.

Triaxial compression tests
Experimental protocol and results. A British GDS high precision soft rock rheometer was used in the test, which included a computer and full servo motor control that automatically collected data. The system has good dynamic response function, including a 250-kN motor driven digital load rack, 32-MPa pressure/volume control system, local strain sensor, and multi-function test module. The rheometer has high testing accuracy and meets the requirements of conventional triaxial and creep tests. The experimental loading device and control system are shown in Fig 2A, and the characteristics of the rock specimens after loading failure are shown in Fig 2B. The displacement control method was used in the triaxial compression tests and the confining pressures were set to 0, 1, 3, 7, and 10 MPa. The full stress-strain curves of the specimens under different confining pressures were obtained by applying an axial load at the loading rate of 0.1 MPa/min, as shown in  Table 1.

PLOS ONE
We use M to represent the post-peak stress-strain of the mudstone, i.e. post-peak softening modulus, as shown in Fig 4. The experimental analysis shows a nonlinear correlation between the post-peak softening modulus and confining pressure, which is obtained by nonlinear exponential fitting. The fitted equation is as follows: where Y is the target parameter, which represents softening modulus M in this equation, σ 3 is confining pressure, and a, b, and c are the constants obtained by the fitting ( Table 2).

Creep tests
Experimental protocol. Mudstone creep tests were carried out using a British GDS high precision soft rock rheometer. We adopted hierarchical loading and applied a maximum load of 80%-90% of the peak strength from the conventional triaxial compression tests under the same confining pressure. The tests were divided into 6-8 loading gradients, as shown in Table 3. During the creep tests, the confining pressure was loaded to the target value using a stress loading rate is 0.05 MPa/s and then held fixed. The axial compression was loaded to the target value of the first stage creep and held constant. The following load was applied every 48 h until the specimen failed.
Creep curve characteristics and fitted parameters. Fig 5 shows the axial creep curves of the mudstone under different confining pressures. The creep curve exhibits the typical threestage creep characteristics of rock. When σ 3 = 0, the mudstone shows characteristics of decay creep and steady creep while loading from stages 1 to 5. When the deviatoric stress reached 3 MPa, the mudstone entered the accelerated creep stage and the rock specimen was destroyed. The mudstone creep curves characteristics under different confining pressures show that the steady creep curve under high deviatoric stress gradually tends to flatten with the increasing confining pressure. This shows that confining pressure exerts a strong influence on the viscoplastic rheological properties of mudstone.
The Burgers model was used to fit and analyze the creep experimental curve of the mudstone. The Burgers model is composed of a Maxwell model and Kelvin model in series, which well reflects the decay creep and steady creep stage of rock. The specific methods are as follows: The GDS file obtained from the experiment was imported into an Excel file and the experimental data were processed by the Boltzmann superposition principle. The expression is as follows [33]: where J is the creep compliance, t is the creep time, and i is the number of hierarchical loading stages, which are taken as 1, 2,. . ., i − 1. The processed data were then imported into Origin software. Because the tests were performed in stages, the data were fitted step by step. The obtained nonlinear fitted equation of the Burgers model suitable for the hierarchical loading mode is given as: where t is creep time, y is creep strain,

PLOS ONE
parameters η m , η k , E m , and E k were inversely calculated according to the fitted results. The obtained parameters are listed in Table 4 and the fitted and experimental curves are compared in Fig 6. The effect of confining pressure on the viscosity coefficients based on the Maxwell model and Kelvin model is shown in Fig 6. The viscosity coefficient η k decreases strongly with increasing confining pressure in the low pressure range (0-3 MPa) and tends to stabilize in the higher pressure range (3-7 MPa). The viscosity coefficient η m increases gradually with the increase of confining pressure. The relationship between the viscosity coefficient and confining pressure is obtained by exponential fitting analysis of the viscosity coefficients η m and η k from the Maxwell and Kelvin models, respectively. The data are fitted following Eq (3) and the obtained parameters are listed in Table 5. The R 2 of the fitted equations all exceed 0.95, which indicates high accuracy.

NBSS model and numerical implementation
Irreversible viscoplastic deformation occurs during the steady mudstone creep that leads to gradual strength attenuation. Accelerated creep failure occurs when the strength is less than or equal to the load stress level [34,35]. In practical engineering, mudstone strength continues to decrease to the residual strength after creep failure.

Strain-softening model
In the plastic strain-softening model of rock, the Mohr-Coulomb yield criterion can be expressed as follows [36]: f ðs y ; s r ; ZÞ ¼ s y À N φ ðZÞs r À s c ð6Þ

PLOS ONE
where σ c is the uniaxial compressive strength of the rock and: If c and φ in the Mohr-Coulomb constant are assumed to vary with η, the Mohr-Coulomb strain-softening relation in Eq (2) can be obtained and c and φ can be replaced by ω. Strainsoftening in mudstone is controlled by the softening modulus M. The critical value that mudstones reaches in the residual strength stage, η � , can be expressed in the form of internal variables. The plastic parameter η � can be defined as the plastic shear strain, which is obtained by the difference between the maximum and minimum principal plastic strains as follows [27]: where ε p 1 and ε p 3 are the maximum and minimum plastic strains, respectively. According to the geometric relationship of ε p 1 and ε p 3 in Fig 7, we obtain: where ε p;e 1 is the maximum elastic principal strain prior to reaching the peak value, ε drop 1 is the

PLOS ONE
softening strain after reaching the peak value, and ε e 1 is the maximum elastic principal strain. When σ 3 is constant, the parameters in Eq (9) can be expressed as follows: where s p 1 is the peak principal stress and s r 1 is the residual principal stress. Considering the dilatancy angle ψ and form of the strain increment, we obtain [37]: where ε& p j (j = 1, 3) is the plastic principal strain rate. Eq (11) can be transformed into: Under static conditions, the relationship between ε p 3 and ε p 1 can be constructed by holding ψ fixed according to Eq (12) as follows: The nonlinear expression of η � is then simultaneously obtained by Eqs (8)-(14) as follows: where:

NBSS creep damage model
In The Burgers creep model can reasonably describe the decay creep and steady creep processes, and the stress-strain deviator constitutive relation is [38]: The axial rock creep strain under deviator stress (-) and confining pressure σ 3 can be deduced according to the Burgers model as: where: .. ..
Based on the Mohr-Coulomb criterion, f = 0, the shear yield and tensile yield are obtained in the principal axis stress space as follows [36]: where σ t is the tensile strength.
Let c p and φ p be the initial cohesion and internal friction angle of rock, respectively. During the creep process, c and φ of rock vary with the plastic parameter as: where η varies with η m . By substituting Eqs (22) and (23) into Eqs (20) and (21), the timedependent shear yield criterion of the nonlinear strain-softening plastic elements can be obtained: In FLAC 3D , for the Kelvin model: and for the Maxwell model: The deviator strain rate of nonlinear strain-softening model is: The volume behavior is given as: In FLAC 3D , the plastic parameter η is expressed by the plastic shear strain, and its increment form is given as: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðDε ps where: where Dε ps j (j = 1, 3) is the main plastic shear strain increment. When considering rock dilatancy behavior, the unrelated flow law corresponding to the plastic shear potential function can be expressed as:

Model establishment
Using FLAC 3D software, the nonlinear creep and strain softening coupling model (NBSS) is used to calculate the creep curve of mudstone. The numerical model and boundary conditions are shown in Fig 9. We use a 3D numerical specimen (50 mm in diameter and 100 mm in length) divided into 3100 elements. The bottom of the specimen is constrained and the upper part is loaded. The experimental uniaxial creep parameters are selected for calculation and the

PLOS ONE
viscosity coefficients η m and η k adopt the fitted parameters of the nonlinear equation in Table 5. The nonlinear Burgers model (η m (σ 3 , η), η k (σ 3 )) is used to calculate the steady creep of mudstone. The nonlinear creep damage variable function η D (σ 3 , η m , η) and constant η D are used to calculate the accelerated creep of mudstone, yielding constant values of 0.061, 0.0635, 0.066 and 0.0685. FISH language is used in the calculation to obtain the nonlinear creep damage variable function η D (σ 3 , η m , η).

Model verification
The results of NBSS numerical model were compared with the experimental data of yellow mudstone uniaxial creep to verify the reliability of the model. The steady creep curve of mudstone calculated by the NBSS numerical model is shown in Fig 10A. The numerical results are in good agreement with the experimental data in the decay creep and steady creep stages, which shows that the nonlinear Burgers model (η m (σ 3 , η), η k (σ 3 )) based on the variation of confining pressure is highly reliable. The accelerated creep curve is obtained by the calculation, as shown in Fig 10B. When the creep damage variable η D is constant, the creep curves obtained by the numerical calculation show different characteristics. When η D � 0.0635, the creep curve obtained by the numerical calculation shows notable acceleration and linear characteristics. When η D � 0.066, the acceleration characteristic of the creep curve obtained by the numerical calculation is not apparent and tends to be in a steady stage. This shows that the accelerated creep characteristics of mudstone are highly sensitive to changes of η D . The creep curve calculated by the nonlinear function η D (σ 3 , η m , η) shows clear acceleration and nonlinear characteristics. The creep acceleration start-up time of the numerical calculation is longer than that of the experimental data; however, the trends are very close, which demonstrates the high accuracy of the calculations. The NBSS numerical model can meet the accuracy requirements of the steady state and accelerated creep calculations for mudstone, and its nonlinear creep characteristics are more in line with actual project situations.

Numerical simulation method
Numerical calculation method for coupled effect of creep and strain-softening. We constructed a typical numerical model of tunnel surrounding rock to analyze the coupled action law of creep and strain-softening. The size and boundary conditions of the numerical model are shown in Fig 11. The model is divided into 11,000 numerical elements and the rock surrounding the tunnel is in the state of hydrostatic pressure. The original rock stress is 2.5 MPa and the rock mechanical parameters are listed in Table 1. We apply two models in the calculation. The first is the NBMC, which does not take into account the strain-softening characteristics of mudstone after accelerated creep and the Burgers model parameters (η m , η k ) have a nonlinear relationship with σ 3 . The second model is the NBSS, which takes into account the creep and strain-softening characteristics of the surrounding rock. The relationship between η and σ 3 is nonlinear.
The numerical calculation process of the coupled effect of creep and strain-softening of rock surrounding a tunnel is as follows: First, the NBSS model and initial mechanical parameters are given to the numerical element to calculate the initial stress balance. An excavation simulation is then carried out and the deviator stress and minimum principal stress of each model element are obtained by the traversal method. The Mohr-Coulomb yield criterion is used to judge the elastic-plastic state of the numerical element. If the element is in an elastic state, the nonlinear creep parameters of the element are reassigned. The viscoplastic strain (η m (σ 3 , η)) of the Maxwell component in the Burgers model is converted into plastic shear strain, and the c and φ values vary with increasing plastic shear strain. The Mohr-Coulomb yield criterion is then used to assess the elastic-plastic state of the numerical element. If the numerical element enters the plastic state, the viscoplastic strain is not converted into plastic shear strain. The stress is updated according to the plastic shear strain increment reaching after the peak value, and c and φ continue to vary with increasing plastic shear strain. In the numerical

Rheological failure characteristics
Plastic zone distribution characteristics using the NBMC model. The nonlinear Burgers model is combined with Mohr-Coulomb model without considering the strain-softening characteristics of mudstone after accelerated creep. When the creep damage of the mudstone causes the surrounding rock to reach the Mohr-Coulomb yield strength, the surrounding rock will enter the plastic state. But the strength of the plastic surrounding rock will not attenuate with increasing plastic shear strain. The balance calculation is carried out after the excavation of the surrounding rock. The concentrated stress in the excavation boundary of the cavern does not exceed the mudstone strength. If the creep characteristics are not considered, there is no plastic zone in the surrounding rock. Under the conditions of no support, the evolution law of the plastic zone is calculated as shown in Fig 13. The plastic zone of the surrounding rock gradually expands over time. In the first 200 h, the plastic zone of the surrounding rock is small; between 300 and 500 h, the plastic zone notably increases obviously; and between 600 and 800 h, the expansion range of the plastic zone tends to gradually stabilize with a plastic zone thickness of 2.6 m.
Distribution characteristics of the plastic zone using the NBSS model. The nonlinear Burgers model and strain-softening model are combined and consider the strain-softening characteristics of mudstone. When the creep damage causes the rock to reach the Mohr-Coulomb yield strength, the rock enters the plastic strain-softening state and the post-peak strength decreases with increasing plastic shear strain. The balance calculation is carried out  Fig 14. The plastic zone of the surrounding rock gradually expands with time. In the first 300 h, the range of the plastic zone is small, the maximum plastic zone thickness is 0.4 m, and the distribution is more uniform and complete. After 400 h, the plastic zone substantially increases and hosts the formation of localized fracture zones. When entering the strain-softening stage, the rock mechanical parameters in the local fracture zone strongly decrease, as shown in the red area in the Fig 14A. After 500 h, the plastic zone of the surrounding rock continues to increase. Localized fracturing is more apparent and the fracture zone gradually extends to the interior of the surrounding rock. After 700 h, the range of the plastic zone tends to gradually stabilize and the expanding thickness of the local fracture zone is 4.0 m.
The above analysis shows that the plastic zone distribution characteristics obtained by the NBSS and NBMC models differ substantially. Fig 15 shows the variation of plastic zone thickness with time calculated by the two models. When the NBMC model is used to calculate the plastic zone, the distribution is relatively uniform and the development range is small, which is consistent with theoretical analysis. When the NBSS model is used, the creep damage causes the surrounding rock to enter the plastic state, the strength of the surrounding rock continues to decline with the plastic shear strain, and the plastic zone of the surrounding rock is larger with more local fracture characteristics. Plastic zone distribution characteristics under anchor rod support conditions (NBSS). Anchor rods are used to support the surrounding rock and their stress characteristics and influence on the plastic zone are analyzed under these conditions using the NBSS model [39]. We adopted full-length bonded anchor rods with a diameter of 18 mm, length of 2.2 m, and spacing of 1.0 m as the support structure. The calculation result is shown in Fig 17. Under anchor rod support conditions, the expansion range of the plastic zone is well controlled. Within 300 h, the expansion range of the plastic zone is small; after 400 hours, the surrounding rock in the anchorage zone appears to locally fracture. However, with increasing time, the range of the plastic zone does not further expand. The anchor rods therefore play a superior supporting role and effectively control the development of plastic zone. Fig 18 shows the evolution process of the anchor rod axial force with time. The anchor rod axial force is small in the first 300 h but gradually increases after 400 h even though the extension range of the plastic zone is effectively controlled by the anchor rod. After 500 h, the growth rate of the anchor rod axial force gradually decreases. The maximum anchor rod axial force is 240 kN at 1000 h, which is concentrated at an anchor rod length of approximately 1.5 m.

Discussion
In order to study the deformation law of soft surrounding rock and provide a theoretical basis for soft rock tunnel support engineering, this paper firstly analyzes the creep characteristics of yellow mudstone through laboratory tests, and proposes a nonlinear creep and strain softening coupling model (NBSS). Then FLAC 3D finite difference software was used to numerically solve the NBSS model to verify the reliability of the model. Finally, using FLAC 3D finite difference In rock creep damage research, many scholars' research focuses on how to determine the creep damage variable, and put forward the expression of accelerated creep of rock mechanics model, but there is no in-depth study intensity attenuation laws, and in terms of the new model is put forward and the design is not enough comprehensive has limitations, their calculations often don't meet a complex rheological properties of the soft rock in actual engineering. The nonlinear creep and strain softening coupling model (NBSS) proposed in this paper is different from other research results by avoiding the above research shortcomings and completing the research loopholes. (1) The research of soft rock creep and strain softening coupling rule. (2) The strength attenuation law of soft rock after accelerated creep is reflected in the numerical calculation. (3) The unique result is that the creep damage model of NBSS is more suitable for the analysis of the supporting system of weak surrounding rock in practical engineering because of its larger plastic zone and obvious local fracture characteristics.
Although the proposed NBSS model can reflect the coupled characteristics of rock mass creep and strain-softening, some limitations remain. For instance, the model is more suitable for underground tunnel engineering with relatively intact rock mass but not for jointed or broken rock. The NBSS model can be developed in further research based on discrete element software (e.g., UDEC or 3DEC). The strength damage characteristics of the structural plane can be reflected by constructing the model with a strength parameter that decays with time. The NBSS model can be used to reflect the creep and strain-softening characteristics of intact rock blocks and thus improves the applicability of the numerical method to practical engineering.

Conclusions
We used theoretical analysis methods to reveal the coupled mechanism of creep and strainsoftening of rock surrounding a tunnel. The strain-softening and creep damage characteristics of mudstone were obtained by laboratory experiments. We construct a nonlinear Burgers model and Mohr-Coulomb strain-softening model, and propose a nonlinear Burgers strain- softening model. The numerical calculation of the coupled effects of creep and strain-softening is detailed and the rheological failure characteristics of a circular tunnel surrounding rock are analyzed using FLAC 3D . The following conclusions are obtained from the analysis.
1. The brittle-ductile characteristics of mudstone are apparent after reaching the peak value, and the confining pressure σ 3 strongly influences the post-peak softening modulus M. The nonlinear exponential equation of σ 3 and M obtained by fitting and nonlinear strain-softening mechanical model reflect the post-peak strength attenuation characteristics of mudstone and can also be used as pre-peak creep damage variable parameters. The fitted creep curves at different stress levels show high accuracy. The viscoplastic creep parameters (η m , η k ) of mudstone are highly sensitive to σ 3 and show notable nonlinear characteristics. We establish the nonlinear exponential equation of viscoplastic creep parameters and σ 3 of mudstone. The proposed nonlinear creep damage model is more suitable for numerical development and calculation and more consistent with the evolution characteristics of the stress and displacement fields of complex rock environments in practical engineering.
2. The steady creep curve calculated by the nonlinear Burgers model is consistent with the experimental data, but this model does not reflect the accelerated creep characteristics of mudstone. By taking the softening coefficient η of mudstone as a creep damage variable, η is a constant when variable σ 3 is not considered. η is a nonlinear time-varying function when considering variable σ 3 . Two kinds of creep damage variables are used to calculate the accelerated creep characteristics of mudstone. The accelerated creep curve calculated using the nonlinear creep damage variable shows good agreement with the experimental data. When a constant creep damage variable is used, the accelerated creep characteristics can be obtained but the acceleration curve shows linear characteristics. Therefore, the nonlinear creep and strain softening coupling model (NBSS) of mudstone is more suitable for analyzing the rheological law of surrounding rock in actual engineering.
3. The nonlinear Burgers Mohr-Coulomb model (NBMC) and nonlinear Burgers strain-softening (NBSS) model are used to calculate the plastic zone and deformation law of rock surrounding a tunnel. There are substantial differences between the two models in the distribution characteristics of the plastic zone. When using the NBMC creep damage model, the distribution of plastic zone is uniform and the plastic zone range is small; whereas the plastic zone calculated by the NBSS creep damage model is larger with notable localized fracture characteristics. The extension range of the plastic zone can be effectively controlled by anchor rod supports, but the axial force of the anchor rod increases with time. Anchor rod support does not meet the stability requirements of surrounding rock, thus it is necessary to combine supports to control the rheological deformation of surrounding rock.
The study in this paper is more suitable for analyzing the rheological law of complex and weak surrounding rock in practical engineering and provides a theoretical analysis basis for tunnel support of complex and weak surrounding rock.