Introduction

Large scale engineering projects often have a service life of several decades, and the creep characteristics of rock and the creep characteristics of rock significantly impact the long-term operation of engineering structures, especially in areas such as slope engineering, subway tunnels, mine goafs and deep engineering1,2,3. With the continuous advancement of engineering technology and rapid development of infrastructure construction in China, numerous large-scale projects have emerged, rendering stability research increasingly urgent.

A substantial corpus of engineering practice and research has demonstrated that the surrounding rock of tunnels generally does not immediately lose stability and fail after excavation and completion. Under long-term loading, the vast majority of metro tunnels invariably experience creep, which can lead to instability or failure due to the increasing deformation of tunnel surrounding rock over time. This creep behavior can have a significant impact on metro settlement and long-term stability4,5,6,7. The timeliness of rocks can lead to cracking occurs at stresses below instantaneous strength8,9,10. Therefore, evaluating the creep behavior of surrounding rock in tunnel engineering is a necessary indicator for assessing the stability of tunnels11,12,13.

Most of the metro tunnel construction in central China passes through the cretaceous strata, which are mainly composed of muddy siltstone with both sandstone and mudstone characteristics. It is a soft rock with an uneven texture that is prone to instability over time. Deng et al.14 studied the creep characteristics of argillaceous sandstone based on the Dapingshan tunnel and found the creep caused large squeezing deformation. However, the tunnel was a deep-buried soft rock tunnel in areas with high geo-stress. Therefore, there is relatively little research on muddy siltstone in shallow buried tunnels, so it is necessary to pay close attention to the influence of soft rock creep behavior on stability.

To study the deformation of tunnel surrounding rock caused by rheological effects, some scholars use similarity theory and numerical simulations to study tunnel stability15,16,17. Yang et al.18 conducted large-scale model experiments and found that the deformation of the surrounding rock after tunnel excavation in deep-buried composite strata includes instantaneous deformation and continuous creep deformation. Based on the Burgers model, Kou et al.19 evaluated the impact of different step lengths and support times on tunnel stability through model experiments, and compared and analysis on practical engineering monitoring and numerical simulation results. Song et al.20 focus on the surrounding rock in the deep tunnels carried out triaxial creep tests, and found that the simulated deformation of the rock was consistent with the on-site monitoring results. Wang et al.21 developed the Nishihara model to characterize the creep characteristics of surrounding rock and verified the correctness of the numerical simulation by comparing it with the tunnel displacement monitoring results. Zhang et al.22 simulated and analyzed the creep deformation of surrounding rock and supporting structures during tunnel excavation, and established a tunnel model to analyze the creep deformation characteristics during tunnel operation in cold regions.

However, existing numerical simulations of the surrounding rock of deep buried tunnels cannot be used to describe the creep deformation of shallow buried soft rock tunnels. More specifically, due to the unity of the strata in deep buried tunnels, the numerical simulation models established are mostly based on the analysis of a single stratum, and uniformly distributed loads are applied on top to replace geo-stress23,24. Nevertheless, the unevenness of the overlying strata in shallow buried tunnels can affect the numerical simulation results, making it difficult to simply use uniformly distributed loads as boundary conditions, and is also affected by the coupling effect of seepage25. Therefore, the existing numerical models for tunnel creep deformation are not applicable to shallow-buried soft rock tunnels.

In general, the creep behavior of rock can be described using a creep model. By establishing constitutive equations between creep strain, stress and time, the intrinsic relationship between the three can be obtained to analyze the creep deformation law of rock. The common methods for establishing creep constitutive models are the empirical model method, element model method and fractional calculus method. The empirical model method has a low threshold for use and a simple model composition, but the parameters do not have specific physical meanings and have a limited range of applications26. The fractional calculus model proposes the Abel clay pot element and its constitutive equation based on the fractional calculus operator and theory, but this model is relatively complex and not suitable for subsequent simulation calculations. The element model method refers to establishing a rock creep model by combining different elements in series and parallel. This creep model can reflect the creep characteristics of different rocks, but can only reflect the linear part of rock creep and cannot describe the accelerated creep stage of rocks. So far, scholars have proposed various combination models to quantify the creep behavior of different types of rocks in tunnel engineering27,28. Note that the creep deformation of tunnel surrounding rock after plastic yield is nonlinear, the Burgers model(viscoelastic) and the Nishihara model(viscoelastic-plasticity) are inadequate to reflect the nonlinear deformation of tunnels. To compensate for the shortcomings of common linear components, it is proposed to replace the linear components in the creep element model with nonlinear components and obtain a new nonlinear rock creep constitutive model through the recombination of components. Therefore, it is highly worthwhile to establish a nonlinear creep model to describe the entire creep process. Moreover, few researchers have proposed creep models for the creep deformation of muddy siltstone in a shallow buried tunnel.

The study aims to construct a new nonlinear creep model to describe the time-dependent behavior of muddy siltstone and evaluate the long-term stability of metro cross-river tunnels. Taking the Changsha metro tunnel as an engineering background, uniaxial compression creep tests were conducted on the muddy siltstone in the surrounding rock layer of the tunnel. An improved nonlinear Burgers model was proposed based on experimental results to characterize the creep behavior of rocks throughout the entire process. Further numerical simulations were conducted on the long-term stability of the tunnel by combining creep models.

Project overview

This study focuses on the tunnel project of Changsha Metro Line 6, specifically the interval tunnel under the Liuyang River section, from Huaqiao Station to Longping Rice Museum Station. The platforms of Changsha Metro Line 6 are primarily of the underground island type, and the metro tunnel is in the form of double holes and double lines. The tunnel interval utilizes the shield excavation method. The left tunnel measures 1012.807 m in length, while the right tunnel measures 1007.285 m, making them both long tunnels. The spacing between the two lines in the interval is 13.2–15.2 m, and the minimum radius of the plane curve for both lines is 400 m. The excavation progress of the left line lags behind that of the right line. The tunnel interval features a ‘V’-shaped longitudinal section, resulting in a significant variation in depth, ranging from 9.2 m to 25.9 m. This interval passes through the Liuyang River, running perpendicular to the Bridge, as shown in Fig. 1(a).

Fig. 1
figure 1

Overview of the metro tunnel under the Liuyang River section. (a) Longitudinal profile of the Liuyang River interval tunnel (1:200). (b) Drilling and sampling.

Full size image

The primary surface water river in the tunnel area is the Liuyang River, which experiences a flood season from April to September and a main flood season from May to July each year. The water level of the Liuyang River is subject to seasonal influences. The project site is located in a subtropical monsoon climate with abundant precipitation, experiencing both rain and heat. The urban area has an average annual precipitation of 1394.6 mm. The section under Liuyang River has a minimum overburden of about 7.57 m. The surrounding rock is mainly composed of muddy siltstone with varying degrees of weathering, as shown in Fig. 1(b). This study focuses on the muddy siltstone and aims to investigate how the creep of the surrounding rock affects the stability of tunnels.

Uniaxial compression creep tests

Specimen preparation and experimental methods

The borehole samples of muddy siltstone were analyzed in-house, revealing that sericite was the most abundant mineral, followed by quartz. The samples were then cut and polished to produce indoor rock mechanics standard specimens with a diameter of 50 mm and a length of 100 mm, in accordance with the guidelines set by the International Society for Rock Mechanics (ISRM)29. To minimize the end effect of the specimen, the height-to-width ratio was set to 2.0, as suggested by ISRM30. To prevent weathering and alteration of the natural moisture content of the standard rock specimens when exposed to air, they were promptly sealed and stored in water-retaining film after cutting. A series of conventional physical and mechanical parameters were measured in the laboratory, and the results were averaged to obtain the mechanical parameter results shown in Table 1.

The uniaxial compression and creep test of muddy siltstone was conducted using the WHY-300/10 microcomputer-controlled pressure tester. The tester has a maximum load capacity of 300 kN with a relative error range of ± 1%, as shown in Fig. 2. To enhance the accuracy of the test, the displacement monitoring system utilized three sets of digital micrometers to collect the specimen’s creep deformation data.

Table 1 The main parameters of the muddy siltstone.
Full size table
Fig. 2
figure 2

Schematic diagram of the experimental system.

Full size image

This article presents research conducted through graded loading creep tests. To ensure that the results of the creep test are reliable, the first level of loading is set at 25% of the maximum uniaxial compressive strength, which corresponds to an initial stress loading level of 6 MPa. This approach accounts for the large dispersion of uniaxial compressive strength of the muddy siltstone. To reduce test duration, stress increments of 4 MPa are used in the early stages of testing. As the stress approaches the destructive strength of the specimen, the increment is reduced to 3 MPa. The stress is then increased gradually until creep damage of the muddy siltstone occurs. The phenomenon of rock creep in engineering can occur over several decades. However, experimental creep tests are subject to various constraints, and the applied load must be maintained for a long time during the test. To reduce test data errors, it is necessary to repeat the creep test. In this paper, the loading time for each stress level was determined to be 120 h after comprehensive consideration. The loading program of creep test is shown in Fig. 3. During the test, the loading rate of the stress growth section was set to 500 N/s, ignoring the effect of stress increase time on the test results, and the temperature and humidity of the test environment were kept as consistent as possible.

Fig. 3
figure 3

The loading program of creep test.

Full size image

Analysis of the creep test results

To reduce experimental errors, three specimens (MS1 ~ MS3) were carried out for repeated creep testing. Considering that there was no significant error in the repeatability creep test, specimen MS1 was selected for detailed analysis, Fig. 4(a) summarizes the instantaneous strain and creep strain at each creep loading level of specimen MS1. The creep of rock is characterized by nonlinearity and does not satisfy Boltzmann’s principle of linear iteration. Therefore, Chen’s loading method31 was used to process the creep test data of the muddy siltstone. With the loading pressure increasing, the instantaneous strain and creep strain increased. To gain the strain growth rate in the steady creep stage, this paper deemed that when the axial strain rate of the specimen is less than 0.001 mm/h, the creep deformation of the specimen has stabilized. With the loading pressure increasing, the steady-state growth rate was gradually increased, as shown in Fig. 4(b). A constant strain growth rate means the rock specimen had entered the steady creep stage. When there is no accelerated creep stage, the strain in the decaying creep stage increases with the increasing pressure, while the decaying creep strain in the failure pressure level decreases. The creep failure stress is close to the uniaxial compressive strength of muddy siltstone, the particle direction movement rate inside the rock is faster than before, and the time to enter the steady creep stage is shorter. The deformation caused by the closure of pores and cracks inside the rock specimens under loading pressure mainly occurs during the decaying creep stage.

Fig. 4
figure 4

Creep test results of muddy siltstone under different loading levels. (a) Strain growth of specimen MS1. (b) Steady-state growth rate of three specimens.

Full size image

The resulting creep curve of total creep process and each stress level is depicted in Fig. 5. Under each stress level, the deformation characteristic of muddy siltstone specimen was described as follows:

  1. (1)

    In the decaying creep stage, strain growth rate of specimen decreased rapidly, and strain increased rapidly.

  2. (2)

    In the steady creep stage, strain growth rate was kept at a constant value, and strain had two different curve expression styles. Under first and second level loading stress, the deformation was tended into a constant value, and strain growth rate could be seen as zero. However, under third to fifth level loading stress, the deformation was uniformly growing, and strain growth rate was a constant value greater than zero.

  3. (3)

    Under sixth level loading stress, the specimen entered the accelerated creep stage after decaying and steady creep stages, and the strain growth rate increased faster and faster until specimen destruction.

Fig. 5
figure 5

Uniaxial creep compression curves of muddy siltstone. (a) The total creep process of three specimens. (b) Each stress level of specimen MS1.

Full size image

To sum up, the creep characteristics of muddy siltstone are relatively consistent with those of soft rock, showing the decaying and steady creep stages under low loading stresses, while show all three creep stages under high loading stress. In the accelerated creep stage, the ratio of creep strain to total strain is relatively high, which is the maximum creep strain generated under all loading stress levels. The joints and cracks inside the rock begin to initiate and expand rapidly, and obvious cracks and failure trends appear on the surface of the specimen. The strain increases rapidly and the rock undergo creep failure in a short time. Therefore, the creep strain of the muddy siltstone specimen mainly occurs in the accelerated creep stage, and the ratio of creep strain at other stress levels to the total strain decreases with increasing stress. It is found that under high stress levels, creep accompanies the entire life cycle of the rock, and the accelerated creep characteristics of the rock become more and more obvious with the increase of stress.

Creep model of the muddy siltstone

Build of creep model

From the creep curves of muddy siltstone specimen, the creep characteristics can divide into three stages. In most cases, the generalized Kelvin, Burgers and Nishihara models can be used to depict the strain curves with multiple creep characteristics32,33,34. The creep strains in the Burgers model is expressed by follows:

$$\varepsilon \left( t \right)=\frac{\sigma }{{{E_1}}}+\frac{\sigma }{{{E_2}}}\left( {1 - {e^{ - \frac{{{E_2}}}{{{\eta _1}}}}}} \right)+\frac{\sigma }{{{\eta _2}}}t$$
(1)

where E1 is elastic modulus; E2 is viscoelastic modulus; η1 and η2 are the viscous coefficients of the decaying and steady creep stages, respectively; σ is the stress on the model.

The Burgers model can accurately describe the decaying of the rock creep process and the deformation characteristics of the steady creep stage. However, since the components used in the Burgers model are ideal linear components, it is difficult to describe the nonlinear deformation law of the rock in the accelerated creep stage. Corresponding to the subsequent simulation, MS1 was considered as a typical specimen to fit the test data using the Burgers models. The fitting curve and parameters value are as shown in Fig. 6. and Table 2.

Table 2 Burgers model fitting parameters of MS1 specimen.
Full size table
Fig. 6
figure 6

Fitting curves of MS1 specimen creep test data with Burgers model. (a) The first to fifth stress levels. (b) The sixth stress level.

Full size image

From Table 2, the R2 of stress level between 6 MPa and 22 MPa is larger, which means the reliability is higher. Therefore, the Burgers model can describe the decaying and steady creep stages. The MS1 specimen came into the accelerated creep stage in the sixth stress level, the deformation increased rapidly and presenting an upward concave curve. While the Burgers model was linear creep model, which could not accurately describe the nonlinear deformation in accelerated creep stage.

Considering the strain growth rate rapid growth and the deformation shows a nonlinear increase curve, the introduction of creep damage in the accelerated creep stage can better describe the creep process of muddy siltstone. Therefore, based on the Burgers model, a damage variable D was introduced and established a new nonlinear damage creep model as shown in Fig. 7. It is worth noting that, the nonlinear component only works when the loading level exceeds the stress threshold, a plastic component should be connected in parallel with the nonlinear viscous component. Although damage behavior also occurs during the decaying and steady creep stages of rock, the accumulated damage is small and negligible compared to the damage accumulation during the accelerated creep stage.

In summary, the damage factor variable D is introduced only in the accelerated creep stage of the rock as a way to quantitatively describe the damage evolution behavior in this stage. Without considering the damage behavior during decaying and steady creep stages, and ignoring the initial damage of the rock, the damage variable D is defined as follows:

$$D=1 - {e^{ - \alpha t}}$$
(2)

where α is the coefficient related to the properties of rock, α > 0; the damage variable 0 < D < 1.

The strain of the nonlinear viscoplastic model can be expressed by Eq. (3).

$$\varepsilon =\frac{{\sigma - {\sigma _S}}}{{\alpha {\eta _3}}}{e^{\alpha t}}$$
(3)

where σs is the yield limit, which serves as the stress threshold to determine the accelerated creep stage of the specimen; η3 is the viscous coefficients of the accelerated creep stages,

Fig. 7
figure 7

Diagram of nonlinear damage creep model.

Full size image

A series connection of the Burgers and a nonlinear viscoplastic model was used to construct the damage creep model of muddy siltstone, and the total strains are expressed by follows:

$$\left\{ \begin{gathered} \varepsilon =\frac{\sigma }{{{E_1}}}+\frac{\sigma }{{{E_2}}}\left( {1 - {e^{ - \frac{{{E_2}}}{{{\eta _1}}}t}}} \right){\text{ }}+\frac{\sigma }{{{\eta _2}}}t,{\text{ }}\sigma <{\sigma _S} \hfill \\ \varepsilon =\frac{\sigma }{{{E_1}}}+\frac{\sigma }{{{E_2}}}\left( {1 - {e^{ - \frac{{{E_2}}}{{{\eta _1}}}t}}} \right)+\frac{\sigma }{{{\eta _2}}}t+\frac{{\sigma - {\sigma _S}}}{{\alpha {\eta _3}}}{e^{\alpha t}},{\sigma _S} \leqslant \sigma \hfill \\ \end{gathered} \right.$$
(4)

If the loading stress is below the threshold, the model is change to the Burgers model, where the first and second stages of creep occur, and the creep strain rate is a constant positive value greater than or equal to zero. If the loading stress is greater than the threshold, the nonlinear viscoplastic model shares a portion of the stress load, and three stages of creep have occurred. The creep strain rate increases with the increase of stress level, and the specimen was broken during the accelerated creep stage. The established creep model was fitted and validated using the creep test results of the MS1 specimen, as shown in Fig. 8. and Table 3.

From Fig. 8, the fitting curve of the nonlinear damage Burgess model is in good agreement with the creep test results of the MS1 specimen. Nonlinear viscoplastic model can effectively describe the accelerated creep stage in creep deformation. Overall, the established nonlinear damage Burgers model can exactly describe the creep characteristics of muddy siltstone, and can be further used for subsequent numerical simulation calculations.

Table 3 Nonlinear damage creep model fitting parameters of MS1 specimen.
Full size table
Fig. 8
figure 8

Fitting curves of MS1 specimen creep test data with nonlinear damage creep model. (a) The first to fifth stress levels. (b) The sixth stress level.

Full size image

Creep equation for the 3D state

In actual geotechnical tunnel engineering, most the surrounding rock of tunnels is subjected to complex three-dimensional(3D) stress effects. Previous studies have shown that confining pressure can affect the creep failure process of rocks35,36,37. Therefore, converting the nonlinear damage Burgess model from one-dimensional to 3D stress state expressions is more conducive to accurately evaluating the long-term stability of metro tunnel engineering.

To more accurately analyze and reflect the creep curve characteristics of rocks under different confining pressures, it is necessary to first obtain the creep constitutive equation under 3D stress state. Assuming the following for rocks in a 3D state:

(1) Rock is a continuous, isotropic material, assuming the same damage in all directions.

(2) Rock damage only occurs during the accelerated creep stage, without considering the initial damage and other stages of damage.

(3) The damage time and creep time are consistent during the accelerated creep stage.

According to the superposition principle, the total strain of the damage creep constitutive model under 3D stress state is:

$${\varepsilon _{ij}}=\varepsilon _{b}^{{}}+\varepsilon _{d}^{{}}$$
(5)

where εij is total strain; εb and εd are the strain of Burgers model and a nonlinear viscoplastic model, respectively.

The creep model equations under a 3D stress state, the stress tensor σij at any point inside the rock under triaxial compression is composed of the spherical stress tensor σmδij and the deviatoric stress tensor Sij. Similarly, the strain tensor εij at any point inside the rock is composed of the spherical strain tensor εmδij and the deviatoric strain tensor eij.

$$\left\{ \begin{gathered} {\sigma _{ij}}={\sigma _m}{\delta _{ij}}+{S_{ij}} \hfill \\ {\varepsilon _{ij}}={\varepsilon _m}{\delta _{ij}}+{e_{ij}} \hfill \\ \end{gathered} \right.$$
(6)

where δij is the Kronecker tensor.

According to Hooke’s law, the stress-strain relationship under 3D stress state is as follows:

$$\left\{\begin{gathered}{\sigma_m}=3K{\varepsilon_m} \hfill \\ {S_{ij}}=2G{e_{ij}}\hfill\end{gathered}\right.$$
(7)

where K and G are the bulk and shear moduli of the Burgers model, respectively.

Through Eq. (4) to (7), the creep strain equation of the the Burgers model in 3D state can be derived as follows:

$$\varepsilon _{b} = \frac{{S_{{ij}} }}{{2G_{1} }} + \frac{{S_{{ij}} }}{{2\eta _{{2'}} }}t + \frac{{\sigma _{m} \delta _{{ij}} }}{{3K}} + \frac{{S_{{ij}} }}{{2G_{2} }}\left( {1 - e^{{ - \frac{{G_{2} }}{{\eta _{{1'}}^{{}} }}t}} } \right)$$
(8)

where G1 and G2 are the shear modulus in the Burgess model, corresponding to the elastic modulus in the one-dimensional constitutive equation; η1’ and η2’ are the viscosity parameters of the Maxwell and Kelvin models, respectively.

The 3D creep equation for a damage nonlinear viscoplastic model is:

$$\varepsilon _{d}^{{}} = \frac{{e^{{\alpha t}} }}{{\alpha \eta _{{3'}} }}\left\langle {\Phi \left[ {\frac{F}{{F_{0} }}} \right]} \right\rangle \frac{{\partial Q}}{{\partial \sigma _{{ij}} }}$$
(9)

where η3’ is the viscosity parameter of the damage nonlinear viscoplastic model; <Ф[F/F0] > is a switch function, as follows:

$$\left\langle {\Phi \left[ {\frac{F}{{{F_0}}}} \right]} \right\rangle =\left\{ \begin{gathered} 0{\text{ }},F<0 \hfill \\ \Phi \left[ {\frac{F}{{{F_0}}}} \right],F \geqslant 0 \hfill \\ \end{gathered} \right.$$
(10)

where F is the yield function of the rock; F0 is the initial value of the rock yield function, usually taken as F0 = 1; Ф[F/F0] is generally in the form of a power function, where Ф[F/F0] = (F/F0)k, where k is a specified constant and generally takes a value of 1; Q is the plastic deformation function, using the correlated flow rule, in which Q = F.

Assuming that under normal temperature conditions, ignoring the factor of σmδij, the key determining factor of rock creep is Sij. Therefore, the expression of the yield function F is:

$$F=\sqrt {{J_2}} - {{{\sigma _S}} \mathord{\left/ {\vphantom {{{\sigma _S}} {\sqrt 3 }}} \right. \kern-0pt} {\sqrt 3 }}$$
(11)

where σs is the yield stress value of the rock material; J2 is the second invariant of stress bias.

The axial deformation of the rock specimen is taken to describe the entire creep process of the rock, in which i = j = 1, and the first principal strain is obtained as ε11. According to the superposition principle, the creep constitutive equation of the nonlinear damage creep model under 3D stress state can be built as follows:

$$\left\{ \begin{gathered} \varepsilon _{{11}} = \frac{{\sigma _{1} - \sigma _{3} }}{{3G_{1} }} + \frac{{\sigma _{1} + 2\sigma _{3} }}{{9K}} + \frac{{\sigma _{1} - \sigma _{3} }}{{3G_{2} }}\left( {1 - e^{{ - \frac{{G_{2} }}{{\eta _{{_{1} '}}^{{}} }}t}} } \right){\text{ }} + \frac{{\sigma _{1} - \sigma _{3} }}{{3\eta _{{2'}}^{{}} }}t,\sigma _{1} - \sigma _{3} < \sigma _{S} \hfill \\ \varepsilon _{{11}} = \frac{{\sigma _{1} - \sigma _{3} }}{{3G_{1} }} + \frac{{\sigma _{1} + 2\sigma _{3} }}{{9K}} + \frac{{\sigma _{1} - \sigma _{3} }}{{3G_{2} }}\left( {1 - e^{{ - \frac{{G_{2} }}{{\eta _{{_{{1'}} }}^{{}} }}t}} } \right){\text{ }} + \frac{{\sigma _{1} - \sigma _{3} }}{{3\eta _{{2'}}^{{}} }}t + \frac{{e^{{\alpha t}} \cdot F}}{{\alpha \eta _{{3'}}^{{}} }}\frac{{\partial F}}{{\partial \sigma _{{ij}} }},\sigma _{S} \le \sigma _{1} - \sigma _{3} \hfill \\ \end{gathered} \right.$$
(12)

In conventional triaxial compression tests, σ1, σ2, and σ3 are the three principal stresses of the rock material, and in conventional triaxial compression tests, σ1 > σ2 = σ3. Therefore, the tensor relationship and the expression of J2 in the 3D stress state is expressed as.

$$\left\{ \begin{gathered} {\sigma _m}=\frac{{{\sigma _1}+{\sigma _2}+{\sigma _3}}}{3}=\frac{{{\sigma _1}+2{\sigma _3}}}{3} \hfill \\ {S_{11}}={\sigma _1} - {\sigma _m}=\frac{2}{3}\left( {{\sigma _1} - {\sigma _3}} \right) \hfill \\ \end{gathered} \right.$$
(13)
$$J_{2} = \frac{1}{2}S_{{ij}} S_{{ij}} = \frac{1}{6}\left[ {\left( {\sigma _{1} - \sigma _{2} } \right)^{2} + \left( {\sigma _{2} - \sigma _{3} } \right)^{2} + \left( {\sigma _{3} - \sigma _{1} } \right)^{2} } \right] = \frac{{\left( {\sigma _{1} - \sigma _{3} } \right)^{2} }}{3}$$
(14)

The simultaneous Eq. (11) and Eq. (13) can yield:

$$F=\frac{{\sqrt 3 }}{3}\left( {{\sigma _1} - {\sigma _3} - {\sigma _S}} \right)$$
(15)

After taking its derivative, it is shown as follows:

$$\frac{{\partial F}}{{\partial {\sigma _{11}}}}=\frac{{\sqrt 3 }}{3}$$
(16)

By combining the above formula, the creep equation of the established nonlinear damage creep model in the 3D stress state is:

$$\left\{ \begin{gathered} \varepsilon _{{11}} = \frac{{\sigma _{1} - \sigma _{3} }}{{3G_{1} }} + \frac{{\sigma _{1} + 2\sigma _{3} }}{{9K}} + \frac{{\sigma _{1} - \sigma _{3} }}{{3G_{2} }}\left( {1 - e^{{ - \frac{{G_{2} }}{{\eta _{{_{1} '}}^{{}} }}t}} } \right){\text{ }} + \frac{{\sigma _{1} - \sigma _{3} }}{{3\eta _{{2'}}^{{}} }}t,\sigma _{1} - \sigma _{3} < \sigma _{S} {\text{ }}and{\text{ }}F < 0 \hfill \\ \varepsilon _{{11}} = \frac{{\sigma _{1} - \sigma _{3} }}{{3G_{1} }} + \frac{{\sigma _{1} + 2\sigma _{3} }}{{9K}} + \frac{{\sigma _{1} - \sigma _{3} }}{{3G_{2} }}\left( {1 - e^{{ - \frac{{G_{2} }}{{\eta _{{_{{1'}} }}^{{}} }}t}} } \right){\text{ }} + \frac{{\sigma _{1} - \sigma _{3} }}{{3\eta _{{2'}}^{{}} }}t \hfill \\ {\text{ }} + \frac{{e^{{\alpha t}} \left( {\sigma _{1} - \sigma _{3} - \sigma _{S} } \right)}}{{3\alpha \eta _{{3'}}^{{}} }},\sigma _{S} \le \sigma _{1} - \sigma _{3} {\text{ }}and{\text{ }}F \ge 0 \hfill \\ \end{gathered} \right.$$
(17)

Generally, for one-dimensional(1D) creep models, they can be extended to 3D states based on the assumption of constant bulk modulus, and thus applied to 3D simulation research and engineering sites38. The assumption of constant bulk modulus refers to the assumption that the bulk modulus is a constant, indicating that viscoelastic deformation only occurs in shear deformation and no bulk deformation occurs. Bulk deformation is elastic deformation, and the bulk modulus during creep is a constant and consistent with elastic deformation. According to Table 3, the bulk modulus remains constant at 1.64GPa.

To avoid confusion caused by the conversion of elastic parameters from 1D tensile and compressive moduli to 3D shear moduli, two processing methods are generally adopted. One approach is to directly use the creep equation for experimental data fitting after hypothesis transformation, in order to ensure consistency of the hypothesis under 1D and 3D stress states and avoid model parameter transformation; The second is to still process the data according to the uniaxial creep test method, and then use model parameter transformation relationships (Eq. (18)) under various assumptions to obtain the creep equation parameters in the 3D state, as shown in Table 4.

$$\left\{ \begin{gathered} G_{1} = \frac{{E_{1} }}{{2(1 + \mu )}} \hfill \\ G_{2} = \frac{{E_{2} }}{3} \hfill \\ \eta _{{i'}} = \frac{{\eta _{i} }}{3},i = 1,2,3 \hfill \\ \end{gathered} \right.$$
(18)

where µ is the Poisson’s ratio, the value is 0.31.

Table 4 Parameter conversion results of muddy siltstone under the constant bulk modulus assumption.
Full size table

Numerical simulation of surrounding rock after tunnels excavated

Model construction and parameters

A numerical simulation was conducted to study the long-term stability of metro section tunnels, using the numerical simulation analysis software COMSOL Multiphysics (COMSOL). COMSOL is based on the finite element method, which simulates real physical phenomena by solving partial differential equations (single field) or partial differential equation systems (multi field). In the COMSOL, the built-in creep module does not include the rock creep equation described in this article. However, the unique partial differential equation modeling and solution in COMSOL namely PDE modeling, allows COMSOL to jointly and solve any field coupling equation, enabling secondary development of nonlinear Burgers models.

Table 5 Physical and mechanical parameters of each stratum within the model.
Full size table

A simplified model was established based on the cross-river tunnel section of Changsha Metro Line 6 (Fig. 1(a)), and the 3D view of the model is shown in Fig. 9(a). The diameter of the model tunnel is 6.2 m, the spacing between tunnels is 13.2 m, and the burial depth of the tunnel at the left embankment of Liuyang River is 25 m. The annual water level of Liuyang River is about 29 m. The horizontal (x) and vertical (y) lengths of the model are 250 m and 80 m, respectively. According to geological data, the model is divided into four strata, the parameter of strata is shown in Table 5. The four lateral boundaries of the model are constrained by roller supports, the bottom is constrained by fixed constraints, and the ground part of the top surface is set as a free surface.

Fig. 9
figure 9

Tunnel model. (a) Schematic diagram of model size and mesh division. (b) Schematic diagram of stratigraphic distribution.

Full size image

The gravity effect of the Liuyang River on the riverbed surface is applied in a uniformly distributed load manner, as shown in Fig. 9(b). Due to the small permeability coefficient of the inner layer in the model, the rate of seepage change is very slow. Therefore, it is believed that there is no seepage effect at the boundary far away from the tunnel. The two sides and bottom of the model are set as no flow boundaries, the riverbed side on the upper surface is a permeable boundary, and the ground side is a no flow interface. Set the water pressure boundary based on hydrogeological data, with a water head of 30 m on the ground side and 29 m on the riverbed side of the model.

The calculation steps of fluid structure coupling simulation process: (1) calculation of initial stress field and initial water stress field; (2) Tunnel excavation; (3) Lining support. At this stage, the deformation of the tunnel is mainly elastic-plastic deformation, and the Mohr-Coulomb model is chosen to describe it without considering the influence of rheology. To simulate the actual construction situation of metro tunnels, the right tunnel is excavated first, and then the left tunnel is excavated. Parameters can set to achieve step-by-step excavation of the tunnel and control the excavation progress difference between the left and right tunnels.

Above the mudstone sandstone formation where the tunnel is located is a round gravel layer, and the creep effect can be ignored. The influence of the upper layer creep on the tunnel surrounding rock is relatively small. Therefore, in this numerical simulation calculation, only the creep deformation occurring in the formation where the tunnel is located is considered.

After the tunnel support is completed, the customized creep calculation module is activated with the stable results of tunnel excavation as the initial value. The creep duration is set to 500 days (500 d), and the constructed damage creep constitutive model is used for creep fluid structure coupling analysis and calculation to obtain the creep deformation of the tunnel surrounding rock after 500 d.

To avoid the problem of convergence caused by frequently changing model parameters, according to the 3D model parameter conversion results of the muddy siltstone in Table 4, the average value of different stress loading intervals is taken as the parameter for the tunnel surrounding rock creep calculation. The stress range intervals are 0 to 12 MPa, 12 to 22 MPa, and above 22 MPa, respectively.

Analysis of simulation results

Select the ground side X = 100 m section for analysis, and set up 8 monitoring points on the twin tunnels.

Stress-time relationship of monitoring points during creep process

The curves of the first and third principal stresses over time at the monitoring points located at the section X = 100 m is presented in Figs. 10 and 11. The first principal stress and third principal stress of the twin tunnels at the section is both negative, indicating that the tunnel did not experience tensile stress during creep and was always subjected to compressive stress. Within 500 d creep period after excavation, the stress gradually increases with time increasing. The stress increases rapidly during the first 100 d period, the rate of stress growth slows down and the stress value tends to stabilize while after 100 d period.

Fig. 10
figure 10

Stress-time curve of monitoring points in the left tunnel at a 100 m section. (a) The first principal stress. (b) The third principal stress.

Full size image

By observing the stress values between monitoring points, it can be found that the stress at the arch waist of the twin tunnels is greater than that at the arch crown and inverted arch. The third principal stress value on the right tunnel is greater than that on the left tunnel, because the excavation process on the left tunnel has an impact on the surrounding rock stress of the already excavated right tunnel. Due to the mutual influence between the left and right tunnels, the rock mass at the center area is subjected to more concentrated principal stress, which is reflected in the stress values on inner side of twin tunnels are higher than those on outer side.

Fig. 11
figure 11

Stress-time curve of monitoring points in the right tunnel at a 100 m section. (a) The first principal stress. (b) The third principal stress.

Full size image

With the increasing time, the stress of the surrounding rock near the tunnel has undergone significant changes. During the process of stabilizing the strata and tunnel, the surrounding rock may continue to deform due to the creep of the rock layers. The continuous increase in surrounding rock stress is constantly adjusting until tending towards stability and uniformity. The strata further away from the tunnel are less affected by tunnel excavation, and the stress distribution of the surrounding rock is basically uniformly distributed according to the burial depth, showing a relatively stable state.

Deformation-time relationship of monitoring points during creep process

By tracking the creep deformation of monitoring points in the twin tunnels, the relationship curve between deformation and time is obtained, as shown in Figs. 12 and 13. The deformation rate of monitoring points gradually decreases over time and eventually tends to a fixed value, reflects the vertical displacement first rapidly decreases and then slowly decreases.

A negative vertical displacement indicates settlement at that point, while a positive value means uplift. The horizontal displacement of the arch crown and inverted arches of the tunnels is close to zero, and the horizontal displacement of the surrounding rock near the tunnel is distributed in a butterfly shape. Due to the creep process influence, the monitoring points at the arch waist of the tunnels have a gradually increasing horizontal displacement and approach towards the tunnel hole. Therefore, the deformation pointing towards the center of the circle as the positive value, the horizontal displacement in this case is defined as the convergence value.

From Fig. 12, the vertical displacement of monitoring points at the arch crown and arch waist remains negative, indicating settlement, with the settlement value at the arch crown being greater than that at the arch waist. Due to the creep characteristics of the rocks around the tunnel, the vicinity of the tunnel cavity was continuously squeezed by the deformation of the surrounding rock, resulting in an increasing deformation of the monitoring points and the tunnel shrinks inward. Since the mutual influence between the two tunnels, the settlement value of the arch waist monitoring point on the inner side is slightly greater than that on the outer side. After the tunnel excavation was completed, the inverted arch showed uplift, and the uplift value gradually decreased. The surrounding rock in the inverted arch area undergoes upward uplift deformation result from lateral pressure compression, the uplift value of the inverted arch is much smaller than the settlement value of the arch crown. However, the twin tunnels as a whole showed a tendency to sink.

Fig. 12
figure 12

Vertical displacement-time curve of monitoring points at a 100 m section. (a) The left tunnel. (b) The right tunnel.

Full size image

After 500 d creep process, the convergence value generated by the right tunnel excavated first is greater than that the left tunnel excavated later. The variation pattern of the arch waists is basically the same, while the maximum horizontal displacement of the twin tunnels is located at the inner side of arch waists. As the creep progresses, the distance between the tunnel does not shrink, but instead gradually moves away, as shown in Fig. 13. The overall value of horizontal deformation is smaller than vertical deformation, so the deformation of surrounding rock during creep process is mainly in the vertical direction. Both vertical and horizontal deformations occur at the arch waist of the tunnel, so more attention should be paid to the support and stability at the position during the maintenance and monitoring stages of the tunnel.

Fig. 13
figure 13

Convergence value-time curve of arch waist monitoring points in the twin tunnels a 100 m section.

Full size image

Comparison with engineering site monitoring data

To facilitate comparison with the simulation results, select ZDK40 + 737.684 section and YDK40 + 725.751 section for research. Convert the monitoring data into consecutive days as the horizontal axis, and the monitoring time inside the tunnel lasting a total of 109 d (from April 13th to July 31st ).

The tunnel deformation data was processed using the moving average method, and the processed monitoring data was imported into the data analysis software Origin. The exponential function was used to perform regression fitting analysis on the monitoring data, as shown in Fig. 14. The results indicate that the exponential function has a good fitting effect and can be used for predicting creep deformation in the later stage.

To verify the reliability of the simulation results, the on-site data of tunnel surrounding rock deformation after construction was compared with the simulation results, as shown in Fig. 15. The vertical displacement of the two curves follows a similar pattern, increasing continuously over time. The rapid increase in displacement within the first 20 d is the initial elastic-plastic deformation. Subsequently, the settlement deformation rate gradually decreases with time increasing, and the displacement growth slows down. The numerical simulation results show a smoother curve, while the actual monitoring results show fluctuations of different amplitudes. Numerical simulations can match the overall trend of settlement values with the measured results, and the two are in good agreement.

Fig. 14
figure 14

Monitoring section vertical displacement data and fitting curve. (a) Left tunnel arch crown. (b) Right tunnel arch crown. (c) Left tunnel inverted arch. (d) Right tunnel inverted arch.

Full size image

This simulation is based on idealized external conditions, with a constant water level and no consideration of the impact of surface vehicle loads on settlement. However, the tunnel monitoring environment is complex and influenced by multiple factors, resulting in some deviations between the two data. The specific manifestation is that the overall settlement trend is consistent, but the simulated data is smaller than the actual monitoring data.

The main surface water rivers in the tunnel section are Liuyang River and Guitang River. The main flood season is from May to July every year, and the water level of Liuyang River is greatly affected by seasonality. The monitoring data of the tunnel is from April 13th to July 31st, which happens to be during the flood season of Liuyang River, when the water level of Liuyang River is relatively high. Surface water is mainly formed by atmospheric precipitation, and greatly affected by rainfall. The water level fluctuates violently, and the highly permeable fine sand and round gravel strata are lost, resulting in local weak interlayers at the interface of the strata and exacerbating the subsidence deformation during tunnel operation. Further, the increase and fluctuation of water level exacerbate the water content of the strata where the tunnel is located, thereby promoting the creep deformation of the tunnel and being an important cause of settlement in the tunnel structure. The traffic load on the Liuyang River Bridge near the tunnel is transmitted to the riverbed through the bridge piers, and the ground vehicle load also have a certain impact on the settlement of the tunnel. Therefore, the above reasons collectively led to actual tunnel settlement being greater than the simulated calculation results.

Fig. 15
figure 15

Monitoring section vertical displacement on-site data and simulation result. (a) Left tunnel arch crown. (b) Right tunnel arch crown. (c) Left tunnel inverted arch. (d) Right tunnel inverted arch.

Full size image
Table 6 Vertical deformation of twin tunnels.
Full size table

With the passage of time, uneven settlement occurred on the left and right tunnel, and the settlement deformation on the right became more pronounced. Based on the Table 6, the maximum settlement position is at the arch crown of right tunnel. This is because the surrounding rock near the right tunnel continued to be disturbed by the left tunnel after the excavation was completed, making the internal structure of the surrounding rock more loose and more prone to creep deformation. During the tunnel operation period, when the arch crown is prone to significant deformation and cracking, it is necessary to maintain the maximum deformation within the deformation limit and strengthen monitoring measures during subsequent operation. In the maintenance and management process of long-term operation of metro tunnels, timely attention should be paid to the deformation of the tunnel arch crown, and perform secondary reinforcement as appropriate.

Conclusion

The creep characteristics of muddy siltstone were studied using uniaxial compression creep tests, and a nonlinear damage creep model was established based on the creep deformation curve, which can accurately describe the three creep stages of the grouting specimen. Further embedding the damage creep model in the simulation software COMSOL, the established creep-fluid-structure coupled tunnel model can be used to predict and analyze the long-term stability and safety of tunnels. The conclusion obtained is as follows.

(1) According to the principle of graded loading, uniaxial compression creep tests were conducted on the muddy siltstone specimen from the borehole. The specimen undergoes three stages of deformation throughout the creep process: decaying creep stage, steady creep stage, and accelerated creep stage. The total axial strain, instantaneous strain, and total creep of rocks all increase with the increasing stress level. As the loading stress increases, the steady-state creep rate of the specimen also increases.

(2) According to the nonlinear creep characteristics of muddy siltstone, the traditional Burgers model was improved and the creep constitutive equation of the unsteady damage Burgers model was built. By fitting the uniaxial creep test data, the model parameters were inverted and the variation law of stress level was analyzed. The result was verified that the model can accurately describe the entire creep process curve, especially the accelerated creep stage. Furthermore, based on the assumption of constant bulk modulus, the creep constitutive equation under 3D stress state was derived, and the conversion relationship between 1D and 3D state parameters was obtained.

(3) The constitutive equation of the constructed nonlinear damage creep model was embedded in COMSOL, and numerical simulation of the stability of surrounding rock under creep-fluid-structure coupling was carried out for the cross-river tunnel section. From the numerical results, the stress and deformation patterns of key monitoring points in the tunnel surrounding rock during the 500 d creep period were obtained. Comparing the on-site monitoring data with the simulation results, the overall deformation trend of the two is roughly consistent. Therefore, the established creep fluid structure coupling tunnel model can be used to predict and analyze the long-term stability and safety of the tunnel. Based on the vertical displacement analysis of different monitoring points, the most unfavorable monitoring point position is the arch top.