Introduction

Rheological behavior mainly occur in soft and loose rock masses and sediment-filled rock masses. Even for the rock mass with a high hardness, developed joints and fractures can also give rise to obvious rheological properties, which can significantly affect the long-term safety of the support system and the tunnel1,2,3. For anchorage support as one of the major means of tunnel support, how to keep tunnel support stability under the action of groundwater has become a major challenge.

Since the twentieth century, some scholars have started to study rock masses for deep engineering as an important topic4,5,6. Under high in-situ stress conditions, rock masses are in a state of tri-axial compression, and due to high bearing capacity, their deformation and failure exhibit certain time-varying effects and rheological properties7. A number of scholars8,9 have conducted a large number of tri-axial compression tests on surrounding rock masses, and concluded that the failed rock masses can experience brittle failure, followed by ductile failure, as the confining pressure increases. On this basis, the conditions and processes of the occurrence of creep phenomena were comprehensively described10,11,12,13,14. With the continuous improvement of test conditions, a large number of scholars have conducted studies on rheological properties of rock masses and structural planes15,16.

Although many existing studies have analyzed the rheological processes of deep rock masses from different perspectives, most of them haven’t considered the inherent time-varying effects of the rock masses themselves, or the effect of seepage on the support system, or even on the rheological process of bolts. As the depth of excavation in underground projects increases, the rheological phenomenon of the rock mass caused by high in-situ stresses is becoming more common. Many scholars have conducted theoretical and analytical studies on the time-dependent behaviors of rock-supporting systems from various perspectives17,18,19,20,21 simplified the action of a full-length bonded bolt on the surrounding rock as a volume force acting within the rock mass. Wang et al.22,23 treated the deeply-buried circular tunnel rock mass and supporting structure as independent rheological bodies, and established a coupled rheological model for the surrounding rock and end-anchored bolt supporting structure. Chu et al.2 proposed a coupled analytical solution that considers both rock rheology and tunnel development effects, and used five different viscoelastic models to obtain the convergence displacement characteristics of unlined tunnels with and without development effects.

As summarized above, the existing studies fall short of fully revealing the rheological law of tunnel anchoring support systems under seepage flow, in which the selection of the rheological models is too simplistic, with few attention paid to the impact of bolt pre-load on the rheological characteristics of anchoring support structures under seepage flow. Therefore, based on the seepage principle of groundwater and the load transfer mechanism of bolts and surrounding rock, we established rational rheological models for both bolts and surrounding rock to analyze the elastic and rheological properties of the circular tunnel anchoring support system under the influence of seepage flow, which reveals the influence of seepage flow on the surrounding rock and supporting parameters of the anchoring support system, and obtains a viscoelastic solution of the rheological model. Additionally, we conducted numerical simulation on the rheological model to analyze the consistence between the numerical and theoretical analytical solutions, and discuss the rationality and correctness of the theoretical analytical solution. The study results offer a dependable theoretical basis for the study of seepage flow in underground engineering.

Mechanical modelling of tunnel anchorage support system under seepage flow action

Establishment of mechanical model

We selected the distributed force model as the basis of the mechanical model of the anchoring support under seepage flow, which also considers the interaction between the bolt and the surrounding rock (the surrounding rock may deform due to the influence of the seepage flow force and the stress of primary rock. There is a force and reaction force relationship between the bolts and the surrounding rock: Under the action of seepage flow, the bolts deform, the axial force of the bolts changes accordingly, and thus the supporting effect for the surrounding rock is also influenced. This relationship between the bolts and surrounding rock is synergistic, which means that they are coupled. To aid theoretical analysis, we categorized the surrounding rock into two zones: anchored zone and non-anchored zone (original rock zone), as shown in Fig. 1. Hydrodynamic pressure, the pressure exerted by water in the pores on the surrounding rock in the form of seepage, causes groundwater seepage pressure.

Fig. 1
figure 1

Simplified model of tunnel anchoring support system under the effect of seepage flow.

Full size image

It is assumed that the excavation radius of the circular tunnel is r, the water head is ha, the pore water pressure is Pa; the radius of the boundary in the non-anchored area is R, the water head is h1; the seepage radius at the far end is R0, the water head is h0 and the pore water pressure is P0. The length of the bolt L = 4 m, the diameter φ = 20 mm, its average uniform force on the inner wall of the tunnel is T, the pre-stress is T0, and the average uniform force on the outer boundary of the anchored zone is (rT/R)24. Since the considered seepage flow is mainly in the radial direction and the influence of buoyancy is not considered for the time being, we adopted the principle of effective stress for theoretical derivation. All the following stress fields and displacement fields are effective fields.

In order to facilitate the theoretical derivation and analysis, we make the following basic assumptions:

(i) The tunnel is an infinite-length tunnel with equal cross-sections, which simplifies it into a plane strain issue; (ii) Tunnel excavation and bolt installation should be carried out simultaneously, regardless of the excavation effect; (iii) Both the tunnel and seepage flow are axis-symmetric and in accordance with Darcy’s law; (iv) The permeability coefficient of the surrounding rock does not change with time and stress; (v) The surrounding rock is a homogeneous, isotropic and porous continuum; (vi) There is no relative slip between the bolts and the surrounding rock; (vii) It is deeply buried, independent of gravity.

The anchored zone under seepage flow

The mechanical action of the mechanical model of the pre-anchored support system in the non-anchored zone is shown in Fig. 2.

Fig. 2
figure 2

Mechanical state, equilibrium state of the anchorage region and its further simplification in the mechanical model of the anchorage support system under seepage action.

Full size image

According to the mechanical model of the anchorage support system under seepage flow, the bolt force is simplified as a uniformly distributed force acting on the inner and outer boundaries of the surrounding rock in the anchored zone. Therefore, the mechanical problem of the anchored area can be simplified as a problem of a thick-wall cylinder under the uniform distribution pressure on the inner and outer boundaries of the seepage physical action. The state of the element in the anchored zone of the mechanical model, when reaching equilibrium, is shown in Fig. 2.

According to the static equilibrium condition of the anchored micro-element, the equilibrium differential equation of the anchored zone is established:

$$\frac{{\text{d}}}{{{\text{d}}\rho }}\sigma_{\rho 1} + \frac{{\sigma_{\rho 1} - \sigma_{\theta 1} }}{\rho } - \frac{{\gamma \left( {h_{1} - h_{a} } \right)}}{{\rho {\mkern 1mu} \ln \left( \alpha \right)}} = 0$$
(1)

Geometrical and physical equations are combined into:

$$\sigma_{\rho } = \frac{{E_{r} }}{{\left( {1 + \mu } \right)\left( {1 - 2{\mkern 1mu} \mu } \right)}}\left( {\left( {1 - \mu } \right)\frac{{\text{d}}}{{{\text{d}}\rho }}u_{\rho } + \frac{{\mu u_{\rho } }}{\rho }} \right)$$
(2)
$$\sigma_{\theta } = \frac{{E_{r} }}{{\left( {1 + \mu } \right)\left( {1 - 2{\mkern 1mu} \mu } \right)}}\left( {\frac{{\left( {1 - \mu } \right)u_{\rho } }}{\rho } + \mu {\mkern 1mu} \frac{{\text{d}}}{{{\text{d}}\rho }}u_{\rho } } \right)$$
(3)

The general solution for the displacement of the anchored zone can be obtained as follows by substituting Eqs. (2) and (3) into Eq. (1):

$$u_{\rho 1} = C_{1} {\mkern 1mu} \rho + \frac{{C_{2} }}{\rho } + \frac{{\gamma \rho {\mkern 1mu} \left( {1 + \mu } \right)\left( { - 1 + 2{\mkern 1mu} \mu } \right)\left( {h_{1} - h_{a} } \right)\left( {2{\mkern 1mu} \ln \left( \rho \right) - 1} \right)}}{{4\left( { - 1 + \mu } \right)E_{r} {\mkern 1mu} \ln \left( \alpha \right)}}$$
(4)

The general stress solutions of σρ1 and σθ1 can be obtained as follows by substituting Eq. (4) into Eqs. (2) and (3):

$$\sigma_{\rho 1} = \frac{{\left( { - 4E_{r} \left( { - 1 + \mu } \right)\left( {C_{1} \rho^{2} + 2C_{2} \mu - C_{2} } \right)\ln \left( \alpha \right) - 4\left( {h_{1} - h_{a} } \right)\left( {1 + \mu } \right)\left( { - 1/2 + \mu } \right)\rho^{2} \left( {\ln \left( \rho \right) - \mu + 1/2} \right)\gamma } \right)}}{{8\ln \left( \alpha \right)\left( {1 + \mu } \right)\left( { - 1/2 + \mu } \right)\left( { - 1 + \mu } \right)\rho^{2} }}$$
(5)
$$\sigma_{\theta 1} = {\mkern 1mu} \frac{{\left( { - 4{\mkern 1mu} E_{r} {\mkern 1mu} \left( { - 1 + \mu } \right)\left( {C_{1} {\mkern 1mu} \rho^{2} - 2{\mkern 1mu} C_{2} {\mkern 1mu} \mu + C_{2} } \right)\ln \left( \alpha \right) - 4{\mkern 1mu} \left( {\ln \left( \rho \right) + \mu - 1/2} \right)\left( {h_{1} - h_{a} } \right)\left( {1 + \mu } \right)\left( { - 1/2 + \mu } \right)\rho^{2} \gamma } \right)}}{{8\ln \left( \alpha \right)\left( {1 + \mu } \right)\left( { - 1/2 + \mu } \right)\left( { - 1 + \mu } \right)\rho^{2} }}$$
(6)

The general strain solutions of ερ1 and εθ1 can then be obtained:

$$\varepsilon_{\rho 1} = \frac{{E_{r} \left( { - 1 + \mu } \right)\left( {C_{1} \rho^{2} - C_{2} } \right)\ln \left( \alpha \right) + \left( {\ln \left( \rho \right) + 1/2} \right)\left( {h_{1} - h_{a} } \right)\left( {1 + \mu } \right)\left( { - 1/2 + \mu } \right)\rho^{2} \gamma }}{{\left( { - 1 + \mu } \right)E_{r} \ln \left( \alpha \right)\rho^{2} }}$$
(7)
$$\varepsilon_{\theta 1} = \frac{{E_{r} {\mkern 1mu} \left( { - 1 + \mu } \right)\left( {C_{1} \rho^{2} + C_{2} } \right)\ln \left( \alpha \right) + \left( {h_{1} - h_{a} } \right)\left( {1 + \mu } \right)\left( { - 1/2 + \mu } \right)\rho^{2} \left( {\ln \left( \rho \right) - 1/2} \right)\gamma }}{{\left( { - 1 + \mu } \right)E_{r} {\mkern 1mu} \ln \left( \alpha \right)\rho^{2} }}$$
(8)

where ρ is the radial coordinate, θ is the circumferential coordinate, σρ1 is the radial stress of the anchored zone, σθ1 is the tangential stress of the anchored zone, ερ1 is the radial strain of the anchored zone, εθ1 is the tangential strain of the anchored zone, μρ1 is the radial displacement of the anchored zone, Er is the Young’s modulus of the surrounding rock, μ is the Poisson ratio of the surrounding rock, α = R/r, C1 and C2 are the undetermined coefficients.

The non-anchored zone under seepage flow

The mechanical problem of the non-anchored zone can be regarded as the problem of the round hole bearing pressure in an infinite area under the action of the primary rock stress and the action of the seepage force. Since the stress forms of the anchored micro-elements in the anchored and non-anchored zones are the same, the differential equation can be obtained as follows according to the principle of stable continuous seepage flow between layers:

$$\frac{{\text{d}}}{{{\text{d}}\rho }}\sigma_{\rho 2} + \frac{{\sigma_{\rho 2} - \sigma_{\theta 2} }}{\rho } - \frac{{\gamma \left( {h_{0} - h_{1} } \right)}}{{\rho {\mkern 1mu} \ln \left( \beta \right)}} = 0$$
(9)

The general equation for the displacement of non-anchored zone can be obtained by solving Eq. (9) simultaneously with Eqs. (2) and (3):

$$u_{\rho 2} = C_{3} {\mkern 1mu} \rho + \frac{{C_{4} }}{\rho } + \frac{{\gamma \rho {\mkern 1mu} \left( {1 + \mu } \right)\left( { - 1 + 2{\mkern 1mu} \mu } \right)\left( {h_{0} - h_{1} } \right)\left( {2{\mkern 1mu} \ln \left( \rho \right) - 1} \right)}}{{4\left( { - 1 + \mu } \right)E_{r} {\mkern 1mu} \ln \left( \beta \right)}}$$
(10)

Substitute Eq. (10) into Eqs. (2) and (3) to obtain the general stress solutions of σρ2 and σθ2, which are substituted into Eqs. (8) and (9) to obtain the general strain solutions of ερ2 and εθ2:

$$\sigma_{\rho 2} = \frac{{\left( { - 4E_{r} \left( { - 1 + \mu } \right)\left( {C_{3} \rho^{2} + 2C_{4} \mu - C_{4} } \right)\ln \left( \beta \right) - 4\left( { - 1/2 + \mu } \right)\left( {1 + \mu } \right)\rho^{2} \left( {\ln \left( \rho \right) - \mu + 1/2} \right)\gamma \left( {h_{0} - h_{1} } \right)} \right)}}{{8\left( { - 1/2 + \mu } \right)\left( {1 + \mu } \right)\ln \left( \beta \right)\rho^{2} \left( { - 1 + \mu } \right)}}$$
(11)
$$\sigma_{\theta 2} = \frac{{\left( { - 4{\mkern 1mu} E_{r} {\mkern 1mu} \left( { - 1 + \mu } \right)\left( {C_{3} {\mkern 1mu} \rho^{2} - 2{\mkern 1mu} C_{4} {\mkern 1mu} \mu + C_{4} } \right)\ln \left( \beta \right) - 4{\mkern 1mu} \left( { - 1/2 + \mu } \right)\left( {1 + \mu } \right)\rho^{2} \left( {\ln \left( \rho \right) + \mu - 1/2} \right)\gamma \left( {h_{0} - h_{1} } \right)} \right)}}{{8\left( { - 1/2 + \mu } \right)\left( {1 + \mu } \right)\ln \left( \beta \right)\rho^{2} \left( { - 1 + \mu } \right)}}$$
(12)
$$\varepsilon_{\rho 2} = \frac{{\left( {E_{r} \left( { - 1 + \mu } \right)\left( {C_{3} \rho^{2} - C_{4} } \right)\ln \left( \beta \right) + \left( { - 1/2 + \mu } \right)\left( {1 + \mu } \right)\rho^{2} \left( {\ln \left( \rho \right) + 1/2} \right)\gamma \left( {h_{0} - h_{1} } \right)} \right)}}{{\left( { - 1 + \mu } \right)\ln \left( \beta \right)E_{r} \rho^{2} }}$$
(13)
$$\varepsilon_{\theta 2} = \frac{{\left( {E_{r} {\mkern 1mu} \left( { - 1 + \mu } \right)\left( {C_{3} {\mkern 1mu} \rho^{2} + C_{4} } \right)\ln \left( \beta \right) + \left( { - 1/2 + \mu } \right)\left( {1 + \mu } \right)\rho^{2} \gamma \left( {\ln \left( \rho \right) - 1/2} \right)\left( {h_{0} - h_{1} } \right)} \right)}}{{\left( { - 1 + \mu } \right)\ln \left( \beta \right)E_{r} {\mkern 1mu} \rho^{2} }}$$
(14)

where σρ2 is the radial stress of the non-anchored zone, σθ2 is the tangential stress of the non-anchored zone, ερ2 is the radial strain of the non-anchored zone, εθ2 is the tangential strain of the non-anchored zone, uρ2 is the radial displacement of the non-anchored zone, E and μ are the Young’s modulus and Poisson’s ratio of the surrounding rock, respectively, β = R0/R, C3 and C4 are the undetermined coefficients.

Bolt support under seepage flow

After the circular tunnel is supported with bolts, the actual displacement and strain of the surrounding rock under the action of seepage flow can be expressed as:

$$u_{\rho 1}^{\prime} = u_{\rho 1} - u_{\rho 0}$$
(15)
$$u_{\rho 2}^{\prime} = u_{\rho 2} - u_{\rho 0}$$
(16)
$$\varepsilon_{\rho 1}^{\prime} = \varepsilon_{\rho 1} - \varepsilon_{\rho 0}$$
(17)
$$\varepsilon_{\rho 2}^{\prime} = \varepsilon_{\rho 2} - \varepsilon_{\rho 0}$$
(18)

According to the hypothesis presented in this study, the deformation of the anchorage and surrounding rock occurs cooperatively. As a result, the strain ερ of the bolts in the axial direction is equal to that of the surrounding rock in the radial direction. Therefore, the axial force of the bolts can be expressed as:

$$T_{\rho } = E_{b} A_{b} \varepsilon_{\rho } - T_{0} = E_{b} A_{b} \left( {\varepsilon_{\rho 1}^{\prime} \left| {_{\rho = r} } \right. - \varepsilon_{\rho 0}^{\prime} \left| {_{\rho = r} } \right.} \right) - T_{0}$$
(19)

Then the average supporting force exerted by the bolts on the surrounding rock at the inner wall of the tunnel can be expressed as:

$$T = \frac{{T_{\rho } { + }T_{0} }}{{S_{\theta } S_{Z} }} = \frac{{E_{b} A_{b} }}{{S_{\theta } S_{Z} }}\left( {\varepsilon_{\rho 1}^{\prime} \left| {_{\rho = r} } \right. - \varepsilon_{\rho 0}^{\prime} \left| {_{\rho = r} } \right.} \right)$$
(20)

where Tρ is the axial force of the bolt at the radial coordinate of ρ, T0 is the pre-stress applied when the bolts are installed, Eb is the elastic modulus of the bolts, Ab is the cross-sectional area of the bolt body, uρ0 is the radial displacement of the primary rock before excavation, ερ0 is the radial strain of the primary rock before excavation, Sθ and SZ are the line spacing of the bolts in the radial and tangential directions of the tunnel, respectively. uρ1′, uρ2′, ερ1′ and ερ2′ represent the radial displacements in both anchored and non-anchored zones and the radial strains in both anchored and non-anchored zones.

Elastic analytical solution for the mechanical model of tunnel anchorage support system under seepage flow

The general equation for the elastic solutions of the anchored and non-anchored zones contains two unknowns. We determined the four boundary conditions based on the elastic theory and the mechanical model of the anchorage support system under the action of seepage flow. With these boundary conditions, we can determine the unknowns and obtain the elastic analytical solution. In the inner wall of the circular tunnel, the relationship between radial stress and bolt support force is a relationship between force and reaction force, thus the stress boundary condition is:

$$\left. {\sigma_{\rho 1} } \right|_{\rho = r} = - \frac{{T + T_{0} }}{{S_{\theta } S_{Z} }}$$
(21)

At the junction of the anchored and non-anchored zones, the stress boundary condition can be expressed using the following equilibrium equation:

$$\left. {\sigma_{\rho 1} } \right|_{{\rho = R^{ - } }} { + }\frac{r}{R}\left( {\frac{{T{ + }T_{0} }}{{S_{\theta } S_{Z} }}} \right) = \left. {\sigma_{\rho 2} } \right|_{{\rho = R^{ + } }}$$
(22)

As this study only considers the existence of hydrodynamic pressure, the stress boundary condition at the far end can be expressed as:

$$\sigma_{\rho 2} \left| {_{{\rho = R_{0} }} } \right. = \sigma_{0}$$
(23)

At the junction of the anchored and non-anchored zones, the displacement continuity conditions expressed as follows should be met:

$$\left. {u_{\rho 1}^{\prime} } \right|_{{\rho = R^{ - } }} = \left. {u_{\rho 2}^{\prime} } \right|_{{\rho = R^{ + } }}$$
(24)

The unknown parameters C1, C2, C3 and C4 can be determined by combining the above boundary conditions with the equations of elastic solutions of the anchored and non-anchored zones. These parameters can then be used in the general equations of stress, strain and displacement of both zones to obtain the elastic solution of the mechanical model of the tunnel anchorage support system under the action of seepage flow. The calculation details can be found in the Appendix SI.

Example analysis of elastic characteristics

After obtaining the elastic solutions of stress field and displacement field of the mechanical model, we analyzed the elastic characteristics of the anchoring support system under the action of seepage flow based on 11 calculation columns. The values of relevant parameter are listed in Table 1.

Table 1 Value of each parameter.
Full size table

Analysis of normal stress evolution characteristics of surrounding rock

Figure 3 shows the evolution curve of the relationship between normal stress and radial distance of surrounding rock with different seepage flow directions. The absolute value of the radial stress of the surrounding rock goes up with the increase of the radial distance and tends to keep constant; meanwhile, the absolute value of tangential stress decreases and also tends to keep constant eventually. It is important to note that the evolution curve of both radial stress and tangential stress is discontinuous: there is a sudden change in the stress at the anchorage end of the bolt, and the sudden change is the direction of stress absolute value reduction. The stress release of the surrounding rock in the anchored zone is limited by the bolting force, but the stress in the non-anchored zone is less affected by the bolting action, resulting in unequal stress on both sides of the junction between the anchored and the non-anchored zones.

Fig.3
figure 3

Evolution curve of surrounding rock normal stress with radial distance with different seepage flow directions.

Full size image

Additionally, Fig. 3 compares the normal stress variation of the surrounding rock under three seepage flow states. When there is no seepage flow, i.e. h0 = 0 m and ha = 0 m, the normal stress curves finally return to the original rock stress (σ0 = − 1 MPa) without deviation. In contrast, when seepage flow occurs both outside and inside the tunnel, i.e. h0 = 20 m, ha = 0 m, the normal stress curves finally shift downwards, indicating an increase in the absolute value of the stress; when seepage flow occurs from inside the tunnel to outside, i.e. h0 = 0 m, ha = 20 m, the normal stress curves finally shift upwards, indicating a decrease in the absolute value of the stress. This is because when the seepage is from inside the tunnel to outside, the seepage force points towards the inner tunnel, which is consistent with the direction of the primary rock stress, thereby increasing the stress on the surrounding rock. As a result, the absolute value of the normal stress tends to exceed that of the primary rock stress. However, when seepage is from inside the tunnel to outside, this direction is opposite to that of the primary rock stress, which contradicts some of the stress on the surrounding rock.

Figure 4 shows the evolution curves of the relationship between normal stress and radial distance of surrounding rock with different elastic moduli, Poisson’s ratios, elastic moduli of bolt and primary rock stresses. It can be seen from Fig. 4a that different values of elastic modulus of the surrounding rock do not affect the change trend of the normal stress of the surrounding rock with radial distance, but have a great influence on the absolute value of the normal stress in the anchored zone, i.e. the absolute value of the normal stress in the anchored zone decreases with the increase of the elastic modulus of the surrounding rock. When there is no seepage flow and the elastic modulus of the surrounding rock increases from 1 to 2 GPa, the absolute value of the radial stress decreases from 0.0321 to 0.0188 MPa at the junction of the anchored zone and non-anchored zone, a sudden decrease of 41%. Similarly, the absolute value of tangential stress decreases from 0.0080 to 0.0047 MPa, a sudden decrease of 41%. When the seepage flows from outside to inside the tunnel, the absolute value of radial stress decreases from 0.0349 to 0.0202 MPa, a sudden decrease of 42%; the absolute value of tangential stress decreases from 0.00873 to 0.00506 MPa, a sudden decrease of 42%. It is evident that the seepage flow can lead to a sudden decrease of the stress absolute value, and also enhance the influence of the increased elastic modulus on the sudden change rate at the junction of the anchored zone and non-anchored zone. Therefore, the seepage flow facilitates the influence of elastic modulus on the normal stress of the surrounding rock.

Fig. 4
figure 4

Evolution curves of the relationship between normal stress and radial distance of surrounding rock with different elastic moduli, Poisson’s ratios, elastic moduli of bolt and primary rock stresses.

Full size image

As shown in Fig. 4b–d, Poisson’s ratio of surrounding rock, elastic modulus of bolts and the value of primary rock stress don’t influence the trend of the normal stress changing with radial distance, but do have a significant influence on the absolute value of normal stress, especially on the sudden change rate of the stress at the junction. Besides, the variation of these three values has a greater impact on the normal stress of the anchored area, but a lesser impact on that of the non-anchored area. The primary rock stress has a relatively minor effect on the radial stress of the tunnel wall, while the tangential stress has a relatively significant effect on it. Therefore, when seepage flow occurs in the surrounding rock with lower primary rock stress, special attention should be paid to the influence of seepage force on the normal stress of the surrounding rock, and it is suggested that the supports be strengthened and observation be intensified to ensure the safe construction and operation.

Analysis of evolution characteristics of surrounding rock displacement

Figure 5 illustrates the evolution curve of the relationship between the displacement of the surrounding rock and radial distance in different seepage flow directions. As the radial distance increases, the radial displacement of the surrounding rock gradually decreases and eventually converges to a stable value. The displacement curve at the junction of the anchored zone and non-anchored zone is continuous, without any sudden change, which is consistent with the displacement boundary conditions mentioned above. When there is no seepage flow in the circular tunnel, the absolute value of the radial displacement of the surrounding rock converges to zero. When seepage flows from inside to outside the tunnel, the absolute value of the surrounding rock displacement will be further reduced. When seepage flows from the outside to the inside the tunnel, the absolute value of surrounding rock displacement will further increase. Its causes are as follows: When the seepage force is opposite to the original rock stress, the surrounding rock is hindered from deforming, causing a decrease in displacement absolute value, and displacement increases as the radial distance decreases; if the seepage force increases beyond a certain point, the displacement direction changes, and the surrounding rock deforms in the opposite direction, resulting in a positive displacement. When seepage flows from inside to outside tunnel, the seepage force shows a negative correlation with the displacement absolute value; conversely, when seepage occurs from outside to inside tunnel, the seepage force shows a positive correlation. Furthermore, in the tunnel wall, the displacement value is − 4.7 mm when there is no seepage, it is − 3.6 mm when the seepage flows from inside to outside tunnel, and it is − 5.7 mm when the seepage flows from outside to inside tunnel. Therefore, it is suggested that the supports be strengthened and observation be intensified in the radial displacement direction at the tunnel wall when seepage occurs from outside to inside tunnel.

Fig. 5
figure 5

Evolution curve of surrounding rock displacement with radial distance in different seepage directions.

Full size image

Figure 6 shows the evolution curve of the relationship between displacement and radial distance of surrounding rock with different values of elastic modulus, Poisson’s ratio, elastic modulus of bolts and primary rock stresses. From Fig. 6, it is evident that variations in these four parameters do not affect the trend of radial displacement changing with radial distance, but have a great influence on the absolute value of radial displacement. Besides, when the elastic modulus of surrounding rock increases, the decreasing rate of radial displacement drops with the increase of radial distance, especially in the anchored area. Therefore, when the value of surrounding rock elastic modulus is small, supports should be enhanced. As shown in Fig. 6b, when no seepage flow occurs, with the increase of Poisson’s ratio, the absolute value of wall displacement also increases from 4.65 to 5.03 mm, an increase of 0.38 mm; when seepage occurs, the absolute value of wall displacement increases from 5.67 to 6.09 mm, an increase of 0.42 mm. Therefore, seepage exacerbates the impact of increasing Poisson’s ratio on wall displacement. As shown in Fig. 6c, the effect of the elastic modulus of the bolts on the radial displacement of surrounding rock is only gently influenced by the seepage flow, which can be ignored. When the elastic modulus of bolts is increased, it has a significant impact on the supports of surrounding rock with large deformation, but has minimal impact on the supports after seepage. Additionally, it can also be seen from Fig. 6d that an increase in primary rock stress can mitigate the impact of seepage on distal radial displacement. For actual engineering projects, when seepage occurs in circular tunnels with relatively small primary rock stress, supports should be enhanced to prevent large displacements of the surrounding rock.

Fig. 6
figure 6

Evolution of the displacement of the surrounding rock with radial distance with the different parameters.

Full size image

Rheological modelling of tunnel anchorage support system under seepage flow

Rock rheological model is a crucial component of the rheological theory. A well-designed model can effectively express complex time-varying processes and rheological characteristics in a systematic and intuitive manner. It is essential for engineering project to select a model that conforms to the actual rheological characteristics and deformation forms of rock or rock mass to analyze engineering issues. Currently, the common method is to adopt rheological medium model and establish constitutive equation for calculus to analyse the rheological process of rock. This method can describe all the laws in the rheological process and is quite applicable to developing numerical software and engineering applications.

Selection of rheological model

Constitutive equation and rheological model of bolts

Compared to the surrounding rock, bolts have relatively simple material properties and mainly undergo axial deformation. This means that the axial stiffness of bolts is much higher than their transverse stiffness, making them an ideal one-dimensional viscoelastic material. The generalized one-dimensional viscoelastic constitutive equation of a bolt can be expressed as follows when the axial stress is σb and the axial strain is εb:

$$P_{b} \left( D \right)\sigma_{b} = Q_{b} \left( D \right)\varepsilon_{b} ,\;{\text{i}}.{\text{e}}{.}\;\sigma_{b} = \frac{{Q_{b} \left( D \right)}}{{P_{b} \left( D \right)}}\varepsilon_{b}$$
(25)
$$P_{b} (D) = \sum\limits_{K = 0}^{m} {p_{k} } \frac{{\partial^{k} }}{{\partial t^{k} }},\;Q_{b} (D) = \sum\limits_{K = 0}^{m} {q_{k} } \frac{{\partial^{k} }}{{\partial t^{k} }},\;D = \frac{\partial }{\partial t}$$
(26)

where D is the time differential operator, Pb(D) and Qb(D) are the operator function which also represent the viscoelastic parameters of the bolt, pk and qk are the functions that characterize the mechanical properties of the bolt.

One-dimensional Maxwell rheological model is selected to study and analyze the rheological properties of the bolt. Its one-dimensional constitutive equation is as follows:

$$\sigma { + }\frac{{\eta_{2} }}{{E_{0} }} = \eta_{2} \varepsilon$$
(27)

Then the operator function of Maxwell rheological model is:

$$P(D) = 1 + \frac{{\eta_{b} }}{{E_{{\text{b}}} }},Q_{b} (D) = \eta_{b} D$$
(28)

The forward Laplace transform of the operator function of the bolt is required for viscoelastic solutions that involve parameter transformation. Thus, the operator function of the Laplace space of the Maxwell rheological model is as follows:

$$\overline{{P_{b}^{\prime } }} (s) = 1{ + }\frac{{\eta_{b} }}{{G_{b} }}s,\overline{{Q_{b}^{\prime } }} (s) = \eta_{b} s,\overline{{P_{b}^{\prime \prime } }} (s) = 1,\overline{{Q_{b}^{\prime \prime } }} (s) = K$$
(29)

where Eb is the elastic modulus of the bolt, ηb is the viscosity coefficient of the bolt, s is the independent variable of the Laplace space, K is the bulk modulus of the bolt, and Gb is the shear elastic modulus of the bolt.

Constitutive equation and rheological model of surrounding rock

According to the elastic theory, the one-dimensional form of the elastic constitutive relation is σ = Erε, while the three-dimensional constitutive relation can be expressed as:

$$S_{ij} = 2G_{r} {\text{e}}_{ij} ,\sigma_{ij} = 3K\varepsilon_{ij}$$
(30)

where Sij is the stress deflection tensor, eij is the strain partial tensor, σij is the first invariant tensor of stress, εij is the first invariant tensor of strain, i is the row number of each component matrix, j is the column number of each component matrix.

The relationship between K, Gr, Er and μ is as follows:

$$E_{r} = \frac{{9G_{r} K}}{{3K + G_{r} }},\;\mu = \frac{{3K - 2G_{r} }}{{2(3K + G_{r} )}}$$
(31)

Compared with the one-dimensional constitutive relation, it can be obtained that the extension of the viscoelastic constitutive model of surrounding rock to the three-dimensional constitutive relation can be expressed as:

$$S_{ij} = 2\frac{{Q^{\prime \prime } (D)}}{{P^{\prime \prime } (D)}}{\text{e}}_{ij} ,\sigma_{ii} = 3\frac{{Q^{\prime \prime } (D)}}{{P^{\prime \prime } (D)}}\varepsilon_{ii}$$
(32)

In Eq. (32), Q'(D) and P'(D) correspond to Q(D) and P(D) mentioned above, but the viscoelastic modulus and coefficient in the operator functions Q(D) and P(D) need to be replaced by the corresponding shear viscoelastic modulus and shear viscoelastic coefficient in Q'(D) and P'(D). Q''(D) and P''(D) are operators representing the viscoelastic volume deformation of the material. If the volume deformation of the material is elastic, then Q'(D) = K, P'(D) = 1 through the comparison with the Eq. (30). From Eq. (30) and Eq. (31), we can obtain:

$$G_{r} = \frac{Q^{\prime}(D)}{{P^{\prime}(D)}},\;K = \frac{Q^{\prime\prime}(D)}{{P^{\prime\prime}(D)}}$$
(33)

Since the parameters of viscoelastic materials are similar to those of elastic materials, the variation between elastic parameters and three-dimensional viscoelastic parameters can be obtained as follows by substituting Eq. (33) into Eq. (31) :

$$E(D) = \frac{{9Q^{\prime } (D)Q^{\prime \prime } (D)}}{{3P^{\prime } (D)Q^{\prime \prime } ({\text{D}}) + P^{\prime \prime } (D)Q^{\prime}(D)}},\;\mu (D) = \frac{{3P^{\prime } (D)Q^{\prime \prime } D) - 2P^{\prime \prime } (D)Q^{\prime } (D)}}{{2(3P^{\prime } (D)Q^{\prime \prime } (D) + P^{\prime \prime } (D)Q^{\prime } (D))}}$$
(34)

Equation (34) is substituted into Eq. (35) for the forward Laplace transform, and the parametric expression of the Laplace space is obtained as follows:

$$L\left( s \right) = \int_{0}^{\infty } {f(t)e^{ - st} } dt$$
(35)
$$\overline{E}_{r} (s) = \frac{{9\overline{{Q^{\prime } }} (s)\overline{{Q^{\prime \prime } }} (s)}}{{3\overline{{P^{\prime } }} (s)\overline{{Q^{\prime \prime } }} (s) + \overline{{P^{\prime \prime } }} (s)\overline{{Q^{\prime } }} (s)}},\;\overline{\mu } (s) = \frac{{3\overline{{P^{\prime } }} (s)\overline{{Q^{\prime \prime } }} (s) - 2\overline{{P^{\prime \prime } }} (s)\overline{{Q^{\prime } }} (s)}}{{2(3\overline{{P^{\prime } }} (s)\overline{{Q^{\prime \prime } }} (s) + \overline{{P^{\prime } }} (s)\overline{{Q^{\prime } }} (s))}}$$
(36)

where \(\overline{{P^{\prime } }} (s)\), \(\overline{{Q^{\prime } }} (s)\), \(\overline{{P^{\prime \prime } }} (s)\) and \(\overline{{Q^{\prime \prime } }} (s)\) are the operator functions of the viscoelastic model of surrounding rock in Laplace space. Equation (36) is the transform expression of phase-space parameters in viscoelastic analysis.

The surrounding rock can be regarded as a three-dimensional viscoelastic material based on the spatial effect. We select the three-dimensional Burgers model to describe the rheological properties of surrounding rock, and its three-dimensional constitutive equation is:

$$S_{ij} + \left( {\frac{{\eta_{2} }}{{G_{0} }} + \frac{{\eta_{2} + \eta_{1} }}{{G_{1} }}} \right)\dot{S}_{ij} + \frac{{\eta_{1} \eta_{2} }}{{G_{0} G_{1} }}\ddot{S}_{ij} = 2\eta_{2} \dot{e}_{ij} + 2\frac{{\eta_{1} \eta_{2} }}{{G_{1} }}\ddot{e}_{ij} ,\sigma_{ii} = 3K\varepsilon_{ii}$$
(37)

Then the operator function of Burgers rheological model is:

$$P\left( D \right) = 1 + \left( {\frac{{\eta_{2} }}{{G_{0} }} + \frac{{\eta_{2} + \eta_{1} }}{{G_{1} }}} \right)D + \frac{{\eta_{1} \eta_{2} }}{{G_{0} G_{1} }}D^{2} ,Q\left( D \right) = 2\eta_{2} D + 2\frac{{\eta_{1} \eta_{2} }}{{G_{1} }}D^{2}$$
(38)

After Laplace transform, the operator function is:

$$\overline{P^{\prime}} (s) = 1 + \left( {\frac{{\eta_{2} }}{{G_{0} }} + \frac{{\eta_{1} \eta_{2} }}{{G_{1} }}} \right)s + \frac{{\eta_{1} \eta_{2} }}{{G_{0} G_{1} }}s^{2} ,\overline{{Q^{\prime } }} (s) = \eta_{2} s + \frac{{\eta_{1} \eta_{2} }}{{G_{1} }}s^{2} ,\overline{{P^{\prime \prime } }} (s) = 1,\overline{{Q^{\prime \prime } }} (s) = K$$
(39)

where η1 is the viscosity coefficient of the surrounding rock at the creep stage No. 1, η2 is the viscosity coefficient of the surrounding rock at the creep stage No. 2, G0 is the shear elastic modulus of the surrounding rock, and G1 is the shear viscoelastic modulus of the surrounding rock.

Based on the selected mechanical model of the tunnel anchorage support system and the rheological model of the bolt under the action of seepage flow, we obtain the distribution characteristics of the anchorage support system, as shown in Fig. 7. Due to the interaction between the surrounding rock and bolts in the anchored zone, Maxwell model and Burgers model coexist in this zone. However, in the non-anchored zone with primary rock, there is no bolt anchorage and thus only Burgers model existing.

Fig. 7
figure 7

Geometric distribution of rheological model of anchorage support system under the action of seepage flow.

Full size image

Solution of viscoelastic problems

Elastic solids and viscous fluids are two common materials used in engineering projects. Unlike elastic solids, viscous fluids do not have a regular shape. The deformation of viscous fluids is related to the duration of the external load. Only a small stress can make the viscous fluids flow, without the need for the external load exceeding a certain value. The flow deformation scale of the viscous fluid can increase with the loading time. Viscoelastic materials exhibit both elastic solid properties and viscous fluid deformation properties. Within the small deformation range of materials, the equilibrium equation, geometric equation, constitutive equation, and boundary condition of the linear elastic problem are basically the same as those of the viscoelastic problem, except that the viscoelastic solution needs to be completed in Laplace space. Therefore, the relationship between these two problems is called the elastic–viscoelastic correspondence principle. Therefore, for all viscoelastic problems, after solving the corresponding elastic solution, it is feasible to simply make forward Laplace transform of the elastic solution, replace the elastic parameters with viscoelastic parameters in the expression in Laplace space, then substitute it into the operator function (Laplace) of the rheological model selected for the viscoelastic problem, and perform inverse Laplace transform to obtain its viscoelastic solution. The process of parameter replacement in Laplace space can be expressed as:

$$E_{b} \Rightarrow \frac{{\overline{{Q_{b}^{\prime } }} (s)}}{{\overline{{P_{b}^{\prime } }} (s)}}$$
(40)
$$E_{{\text{r}}} \Rightarrow \overline{E}_{b} (s) = \frac{{9\overline{{Q^{\prime } }} (s)\overline{{Q^{\prime \prime } }} (s)}}{{3\overline{{P^{\prime } }} (s)\overline{{Q^{\prime \prime } }} (s) + \overline{{P^{\prime \prime } }} (s)\overline{Q^{\prime}} (s)}},\;\mu_{{\text{r}}} \Rightarrow \overline{\mu } (s) = \frac{{3\overline{{P^{\prime } }} (s)\overline{{Q^{\prime \prime } }} (s) - 2\overline{{P^{\prime \prime } }} (s)\overline{{Q^{\prime } }} (s)}}{{2(3\overline{P^{\prime}} (s)\overline{{Q^{\prime \prime } }} (s) + \overline{{P^{\prime } }} (s)\overline{{Q^{\prime } }} (s))}}$$
(41)

Numerical analysis of rheological properties

On the basis of obtaining the viscoelastic solutions of stress field and displacement field of the rheological model of the tunnel anchorage support system under the action seepage flow, we analyze the rheological properties using 15 calculation columns. The values of relevant parameters are shown in Table 2.

Table 2 Values of each parameter.
Full size table

Analysis of radial stress evolution characteristics of surrounding rock

  1. (1)

    Radial stress on tunnel wall

The evolution curve of the radial stress of the tunnel wall with time in different seepage directions is shown in Fig. 8. With the increase of duration, the absolute value of radial stress on the tunnel wall first increases and then decreases until it converges. The rheological process of radial stress on the tunnel wall is divided into three stages. Stage I features instantaneous increase of stress. Since the tunnel excavation and bolt reinforcement are carried out simultaneously, and time difference effects are not considered, the bolts need to bear the huge pressure brought by the surrounding rock in a short time, leading to a rapid increase of the absolute value of their stress to the peak value. Therefore, the stage I of the rheological process conforms to the model assumption. Stage II features slow release of stress. Since the stress shortly generated by the surrounding rock in stage I is limited by the bolt supports, it can only be slowly released over time. Stage III features stress convergence. In the process of stress release, the absolute value of the surrounding rock stress continues to decrease, and a new relationship between the surrounding rock and bolt supports gradually forms. When the surrounding rock and support structure reach a new balance state, the tunnel wall stress tends to be stable. It can be seen from Fig. 8 that the tunnel wall stress has converged at 2,000 h.

Fig. 8
figure 8

Evolution curve of radial stress on tunnel wall with time in different seepage directions.

Full size image

When the seepage flows from inside to outside tunnel, i.e. when h0 = 0 m and ha = 20 m, the absolute value of the radial stress on the tunnel wall decreases, the rate of stress increase in stage I decreases, the duration is shortened, and the peak absolute value of the stress becomes smaller, the stress release rate in stage II slows down, and the absolute value of the radial stress convergence in stage III also becomes smaller. When seepage flows from outside to inside tunnel, i.e. h0 = 20 m and ha = 0 m, the absolute value of radial stress on the tunnel wall increases, the rate of stress increase in Stage I increases, the duration increases, the peak value of the absolute value of stress increases, the stress release rate in Stage II increases, and the absolute value of stress convergence in Stage III also increases. The reason is that when the seepage flow direction is opposite to the stress direction of the primary rock, the force hindering the stress action of the surrounding rock will be generated, resulting in the reduction of the stress of the surrounding rock and the slowing down of the deformation. When the seepage flow direction is the same as the stress direction of the primary rock, the force that intensifies the stress action of the primary rock will be generated, and the stress of the surrounding rock will increase with the change of time, resulting in larger deformation.

Figure 9a and b show the relation evolution curve of the influence of different viscosity coefficients of surrounding rock on the radial stress on the tunnel wall under the action of seepage flow (h0 = 20 m, ha = 0 m). The viscosity coefficients of the surrounding rock in both stage I and stage II do not change the variation trend of the radial stress on the tunnel wall with time. However, they have a significant effect on the absolute value of stress. Furthermore, Fig. 9a and b illustrates that the viscosity coefficient of the surrounding rock in stage II has a significant impact on the subsequent stress on the tunnel wall, and the two are inversely correlated. The existence of seepage force intensifies the influence of η1 and η2 on the radial stress on tunnel wall.

Fig. 9
figure 9

Evolution curve of radial stress on tunnel wall with different viscosity coefficients of surrounding rock and bolt parameters.

Full size image

The evolution curve of radial stress on tunnel wall with time under the action of seepage flow is shown in Fig. 9c–f. It can be seen that changes of these parameters do not alter the changing trend of radial stress on the tunnel wall with time, but have a significant impact on the value and changing rate of radial stress on the tunnel wall. In Fig. 9c, when no seepage occurs, the initial absolute value of the radial stress on the tunnel wall increases with the increase of the prestressing force T0 of the bolt. The rate of stress increase in stage I increases, the peak value of the stress also goes up, the stress release in stage II speeds up, and the absolute value of stress convergence in stage III increases. When seepage occurs, the effect on the three stages of radial stress rheology is more significant. Therefore, there is a positive correlation between the absolute value of the radial stress on the tunnel wall and prestressing force. In Fig. 9d–f, when the length and diameter of the bolt increase and the spacing between rows decreases, both the rate of stress increase and peak stress value in stage I increase, the stress release in stage II accelerates, and the absolute value of stress convergence increases in stage III. However, it can be found that the longer the bolts, the smaller the amplitude of the change in the three stages of the rheological process. Therefore, the absolute value of the radial stress on the tunnel wall is positively correlated with the length and diameter of the bolt, and negatively correlated with the row spacing of the bolt.

  1. (2)

    Radial stress on surrounding rock.

When of t = 2000 h, the evolution curve of the radial stress on surrounding rock varying with the radial distance under the seepage flow action is shown in Fig. 10. It can be obtained from Fig. 10 that the absolute value of the radial stress on surrounding rock gradually decreases with the increase of the radial distance, and finally keeps stable. When no seepage occurs (h0 = 0 m, ha = 0 m), the absolute value of radial stress under rheological state (t = 2000 h) is greater than that under elastic state (t = 0 h). When seepage flows from inside to outside tunnel, the absolute value of radial stress in rheological state is smaller than that in elastic state. As the absolute value of radial stress in the non-anchored zone decreases greatly, the abrupt value of stress at the junction of the anchored zone and the non-anchored zone increases. The absolute value of radial stress in rheological state is greater than that in elastic state when seepage occurs from outside to inside tunnel.

Fig. 10
figure 10

Evolution curve of radial stress on surrounding rock with radial distance under seepage action when t = 2,000 h.

Full size image

Figure 11 illustrates the evolution curve of radial stress on the surrounding rock with radial distance for different bolt parameters under the seepage action when t = 2000 h. The bolt parameters do not affect the variation trend of radial stress on surrounding rock with time, but they significantly impact the absolute value of radial stress on the surrounding rock. When there is no seepage, the absolute value of radial stress increases with the increase of prestressing force in both the anchored and non-anchored zones. The increase of prestressing force has a greater impact on the anchored zone than on the non-anchored zone due to the effect of bolt supports. This results in a larger abrupt value of stress at the junction. Seepage flow from outside to inside tunnel causes an increase in the absolute value of the radial stress on the surrounding rock and a larger abrupt value of the junction, which is also reflected in the radial stress on the tunnel wall. Increasing the length and diameter of the bolts while decreasing the spacing between rows can result in an increase in the radial stress on the tunnel wall. Additionally, the abrupt stress value at the junction decreases due to the greater influence on the anchored zone compared to the non-anchored zone.

Fig. 11
figure 11

Evolution curve of radial stress on surrounding rock with radial distance for different bolt parameters at t = 2,000 h.

Full size image

Analysis of the evolution characteristics of surrounding rock displacement

  1. (1)

    Displacement of tunnel wall

The evolution curve of the relationship between tunnel wall displacement and time in different seepage directions is shown in Fig. 12. As duration increases, the absolute value of tunnel wall displacement remains stable for a short period of time, then gradually increases, and finally tends towards a stable value. The value and direction of the seepage force do not alter the changing trend of tunnel wall displacement with time, but have a significant impact on the absolute value of displacement. The rheological process of tunnel wall displacement also includes three stages. Stage I is the initial stable stage, stage II is the mid-term increasing stage, and stage III is the late convergence stage. During the initial stage of tunnel anchorage supports, as the stress on the surrounding rock is still in the early stage of rheological process, the tunnel wall displacement remains basically elastic under the action of the surrounding rock support system, with minimal rheological displacement. As time passes, the stress on the surrounding rock gradually increases, and in order to release pressure on the support system, the tunnel wall displacement begins to increase. When the surrounding rock support system reaches a new equilibrium state, the tunnel wall displacement converges. When seepage flows from inside to outside tunnel, the seepage has a relatively smaller impact on stage I of the rheological process of the tunnel wall displacement. In stage II, the increase of the displacement of the tunnel wall slows down. In stage III, the absolute value of the displacement convergence decreases. When seepage flows from outside to inside tunnel, it also has a relatively small impact on the stage I of the rheological process, with the increase of the tunnel wall displacement accelerating in stage II, and the absolute value of the displacement convergence increasing in stage III.

Fig. 12
figure 12

Evolution curve of wall displacement with time in different seepage flow directions.

Full size image

The relationship between the displacement of the tunnel wall and time under the action of seepage flow (h0 = 20 m, ha = 0 m) with different viscosity coefficients of the surrounding rock is shown in Fig. 13a and b. Both the viscosity coefficient η1 in stage I and the viscosity coefficient η2 in stage II of the surrounding rock do not change the trend of the tunnel wall displacement varying with time, but have a significant impact on the absolute value of the tunnel wall displacement. η1 only affects the value of tunnel wall displacement in the initial stable stage and the two are negatively correlated, while η2 has a significant impact on the tunnel wall displacement in the mid-term stage and the two are negatively correlated. In addition, it can be seen from Fig. 13a and b that the presence of seepage force intensifies the influence of η1 and η2 on the tunnel wall displacement.

Fig. 13
figure 13

Evolution curve of tunnel wall displacement with time for different viscosity coefficients of the surrounding rock and bolt parameters.

Full size image

The relationship between the tunnel wall displacement and time under the action of seepage flow with different bolt parameters is shown in Fig. 13c–f. It can be obtained that the change of bolt parameters does not alter the trend of the tunnel wall displacement with time, but has a significant impact on its absolute value. As shown in Fig. 13c, when no seepage occurs, with the increase of bolt prestressing force T0, the initial stable value of tunnel wall displacement remains basically unchanged, the rate of increase in stage II slows down, and the absolute value of displacement convergence in stage III decreases. When seepage occurs, it has a relatively larger impact on the last two stages of the displacement rheological process of the tunnel wall. That is, the larger the seepage force, the faster the displacement increases, and the larger the absolute value of displacement convergence. Therefore, the absolute value of radial stress on the tunnel wall is negatively correlated with prestressing force, but the magnitude of prestressing force change has a smaller impact on the tunnel wall displacement. Similarly, it can be seen from Fig. 13d–f that the absolute value of radial displacement of the tunnel wall is negatively correlated with the length and diameter of bolts, and positively correlated with the row spacing of the bolts.

  1. (2)

    Radial displacement of surrounding rock.

The relationship between displacement and time of surrounding rock under different seepage directions at t = 2000 h is shown in Fig. 14.

Fig. 14
figure 14

Evolution curve of surrounding rock displacement with radial distance in different seepage directions when t = 2,000 h.

Full size image

As the radial distance increases, the absolute value of the radial displacement of the surrounding rock gradually decreases and eventually tends towards a stable value of zero. The value and direction of the seepage force do not change the trend of displacement with time, but have a significant impact on the absolute value of the displacement. It can be seen from the Fig. 14 that when there is no seepage, the absolute value of the surrounding rock displacement is 0.33 m, and the stable value of displacement at the far end is 0.029 m. When seepage flows from inside to outside tunnel, the absolute value of radial displacement of the surrounding rock decreases, and the rate of change with the radial distance also decreases. The absolute value of displacement at the tunnel wall is 0.30 m, and the stable value of displacement at the far end is 0.025 m. When seepage flows from outside to inside tunnel, the absolute value of radial displacement of the surrounding rock increases, and the rate of change with the radial distance also increases. The absolute value of displacement at the tunnel wall is 0.36 m, and the stable value of displacement at the far end is 0.032 m. Therefore, when the direction of seepage is opposite to the direction of primary rock stress, it generates a force that resists the movement of the surrounding rock towards the inner tunnel. When the direction of seepage is the same as the direction of primary rock stress, it generates a force that intensifies the movement of the surrounding rock towards the inner tunnel.

The relationship between the surrounding rock displacement and the radial distance under the action of seepage flow with different bolt parameters when t = 2000 h is shown in Fig. 15. The absolute value of the surrounding rock displacement gradually decreases with the increase of the radial distance, and the curve eventually tends towards a stable value. The changes of the bolt parameters do not change the trend of the surrounding rock displacement with time, but have a significant impact on the absolute value of the surrounding rock displacement. The absolute value of the displacement of the surrounding rock is negatively correlated with the prestressing force, the length, and the diameter of the bolts, and positively correlated with the row spacing of the bolts.

Fig. 15
figure 15

Evolution curve of surrounding rock displacement with radial distance for different bolt parameters when t = 2,000 h.

Full size image

Axial force on bolts

The relationship between the axial force of the bolts and time in different seepage directions is shown in Fig. 16. As time passes, the axial force of the bolts shows a trend of first increasing, then decreasing, and finally converging. The value and directional changes of seepage force do not change the variation law of the axial force of the bolt with time, but have a significant impact on the value of the axial force on the bolts. In the mechanical model of the anchorage support system under the action of seepage flow, it is assumed that there is no relative slippage between the bolts and the surrounding rock, and the two are coupled. This means that the rheological process of the radial stress between the bolt and the surrounding rock wall is similar, thus the rheological process of the axial force on the bolts is also divided into three stages: the stage of instantaneous increase of axial force, the stage of axial force release, and the stage of axial force convergence. When seepage flows from outside to inside tunnel, in stage I of the axial force rheological process of the bolt, the axial force increase speeds up, with the duration shortened and peak axial force increased. In stage II, the axial force release speeds up. In stage III, the axial force convergence value increases. When seepage flows from inside to outside tunnel, in stage I, the axial force increase speeds up, with the duration increased and the peak axial force decreased. In stage II, the axial force release slows down. In stage III, the axial force convergence value decreases.

Fig. 16
figure 16

Evolution curve of bolt axial force with time in different seepage directions.

Full size image

Under the action seepage flow, the relationship between bolt axial force and time for different cohesion coefficients of surrounding rock is shown in Fig. 17a and b. Both η1 and η2 do not change the trend of the bolt axial force varying with time, but have a significant impact on the value of the axial force. As shown in Fig. 17a, as η1 increases, the initial axial force of the bolts decreases, the axial force increase rate in stage I of the rheological process decreases, and the axial force peak value decreases, but stage I takes longer; as time passes, it is clearly observed in Fig. 17 that η1 has little impact on the stages II and III of the rheological process, thus it can be concluded that η1 has a relatively large impact on the early rheological process for bolt axial force and the two are negatively correlated. Similarly, as shown in Fig. 17b, η2 has a relatively large impact on the late-stage stress on the tunnel wall and the two are negatively correlated, and the existence of seepage flow intensifies the impact of η1 and η2 on bolt axial force.

Fig. 17
figure 17

Evolution curve of bolt axial force varying with time for different viscosity coefficient of surrounding rock and bolt parameters.

Full size image

The evolution curve of bolt axial force varying with time for different bolt parameters is shown in Fig. 17c–f. It can be seen that the variation of bolt parameters doesn’t change the trend of bolt axial force varying with time, but has a significant impact on the value and change rate of bolt axial force. As shown in Fig. 17c, the axial force of bolt is positively correlated with the prestressing force. In Fig. 17d–f, the axial force of the bolt is positively correlated with the length and diameter of the bolt, and negatively correlated with the row spacing between the bolts. However, with the increase of length and the decrease of diameter, the influence is weakened.

Adaptive evaluation and variability analysis of rheological models for anchored support systems

Numerical model

A reasonable and accurate mechanical calculation model is the prerequisite for numerical simulations to ensure accurate results. To intuitively verify and analyze the characteristics and laws of the stress field, displacement field, and bolt stress field in the surrounding rock under the rheological model of tunnel anchoring and support systems considering seepage effects, the establishment of the numerical model in this paper follows several principles:

  1. (1)

    It is necessary to abstract and generalize various influencing factors in a realistic manner, ensuring the accuracy of the calculation results while considering other factors to the greatest extent possible.

  2. (2)

    Efforts should be made to ensure that the numerical model matches the theoretical model, fully reflecting the mechanical properties of each material.

  3. (3)

    To improve the precision of the model and ensure that it is not affected by boundary conditions during calculations, the model should be sufficiently large and include adequate boundaries, provided that Saint Venant’s Principle is satisfied.

  4. (4)

    The mesh in the primary research area should be refined, while the mesh size in secondary research areas or marginal areas can be appropriately increased to accelerate the calculation speed of the model.

A numerical simulation model of the rheological model of the tunnel anchorage support system under seepage was established using FLAC3D (Version 7.00)25, as shown in Fig. 18. Circular boundaries are not easy to be fixed and only one stress can be applied to one boundary, thus we used square boundaries instead of circular boundaries, which can effectively apply primary rock stress and facilitate the fixing of boundary conditions. The size of the numerical simulation model is 240 m × 240 m × 40 m, in which the tunnel is excavated in a circular shape with a radius of 4 m and a support radius of 4 m. According to field experience and support theories, the upper surface of the model is set as the stress boundary, the left, right, front, and back surfaces of the model are set as horizontal displacement constraints, and the upper and lower surfaces of the model are set as vertical displacement constraints. To simulate the deep stress environment, a uniform load of 1 MPa is applied on the upper and lower surfaces of the model. The initial horizontal stress in the left and right directions is set to 1 MPa in the model. The model has a total of 28,800 units and 31,691 nodes. Based on the theoretical model, the positive direction of the Y axis is set as the direction of tunnel excavation, that of the Z axis is set as the vertical upward direction, and that of the X axis is set as the left direction of the tunnel support section.

Fig. 18
figure 18

Numerical calculation model of rheology of tunnel anchoring support system under seepage.

Full size image

To simplify the calculation process, the model is treated as the same rock layer, and the seepage force is kept the same in each region, avoiding the complex transfer process caused by the differences between different soil or rock layers. Considering the differences between theoretical models and numerical models, some surrounding rock parameters are specially fine-tuned to ensure that the calculation results are more in line with the actual situation under reasonable conditions. The rock mechanical parameters and seepage parameters are shown in Tables 3, 4 and 5.

Table 3 Mechanical parameters of rock mass.
Full size table
Table 4 Mechanical parameters of bolts.
Full size table
Table 5 Seepage parameters.
Full size table

In order to compare with the theoretical analytical solution, two simulation comparison schemes were designed. Scheme I refers to the seepage flowing to outside tunnel, and Scheme II refers to the seepage flowing to inside tunnel. They are used to compare the stress field and displacement field of the numerical solution and theoretical solution, as shown in Fig. 19.

Fig. 19
figure 19

Two numerical simulation schemes with different percolation directions.

Full size image

Numerical calculation method

Based on theoretical models, in numerical simulations, the interaction between hydraulic and mechanics is considered and calculated simultaneously and interdependently. This approach allows for a more accurate representation of the actual behavior of support systems under seepage, including the dynamic changes in the stress field, displacement field of the surrounding rock, and the stress field of the bolts. In the numerical simulation process, the three-dimensional Burgers viscoelastic constitutive model is used to simulate the rheological properties of the tunnel surrounding rock, and the creep process is determined by the state of partial stress. In the rheological model of the tunnel anchorage support system under the action of seepage flow, the strain of the surrounding rock in the three-dimensional Burgers rheological model can be expressed as:

$$\varepsilon (t) = \left( {\frac{1}{{G_{0} }} + \frac{1}{{\eta_{2} }}t + \frac{1}{{G_{1} }}\left( {1 - e^{{\left( { - \frac{{G_{1} }}{{\eta_{t} }}} \right)}} } \right)} \right)\sigma_{0}$$
(42)

In numerical simulation, the default property of the CABLE element software is of linear elastic plasticity, without time-dependent properties. To verify the stress of the bolts in the theoretical model, a secondary program development shall be performed on the bolt element CABLE to enable it to exhibit some creep behaviors and is embedded in the simulation. When the Maxwell rheological model is selected as the constitutive model for the bolts, since its model consists of a Hooke solid and a Newtonian fluid in series, the value of the anchor axial force should be divided into two parts. One part is the stress increase caused by the transient behavior of the spring, and the other is the stress increase caused by the viscous flow of the Newtonian fluid. In the numerical simulation process, it is assumed that the two ends of the bolts are fixed on the surrounding rock, and their deformation is synchronized with the surrounding rock. Their strain and axial force are calculated based on the relative velocity v0(t) between two points on the surrounding rock. The strain and stress of the bolts can be expressed as a function of time:

$$\varepsilon \left( t \right) = \frac{1}{L}\int_{0}^{t} {v_{0} } \left( t \right)dt$$
(43)
$$\sigma \left( t \right) = E_{b} \cdot \varepsilon \left( t \right) + \eta_{b} \cdot \frac{d\varepsilon \left( t \right)}{{dt}} = \frac{{E_{b} }}{L} \cdot \int_{0}^{t} {v_{0} } \left( t \right)dt + \frac{{\eta_{b} v_{0} \left( t \right)}}{L}$$
(44)

where Eb is the linear elastic modulus of bolts; L is the length of bolt free section; ηb is the viscosity coefficient of bolts.

Calculation and analysis

Using the parameters shown in Tables 3, 4 and 5, we conducted a comparative analysis of analytical and numerical solutions for radial stress, bolt axial force, and displacement of surrounding rock in different seepage flow conditions, as shown in Fig. 20. Due to the simplification of seepage force in the theoretical solution of the model, where it is treated as a seepage volume force acting within the surrounding rock, equivalent to an internal force, during the simulation process, in order to more closely align with the theoretical solution, the FISH was utilized to assign pore water pressure to each node, thereby controlling the seepage force and hydraulic boundary conditions. When seepage flows from outside to inside tunnel, the pore water pressure at the far end is 2 × 105 Pa, and the pore water pressure at the tunnel wall is zero, i.e., it is a full-drainage state. When seepage flows from inside to outside tunnel, the pore water pressure at the far end is zero, and the pore water pressure at the tunnel wall is 2 × 105 Pa.

Fig. 20
figure 20

Comparison of analytical and numerical solutions (A-S: analytical solution; S–S: numerical solution).

Full size image

When t = 2000 h, i.e., in the rheological state, whether seepage flows from outside to inside or from inside to outside tunnel, the numerical simulation results are generally similar to the analytical solutions in terms of overall patterns, and the absolute values of the analytical solutions are slightly larger than those of the numerical solutions. In addition, the absolute values of radial stress on the tunnel wall are relatively larger. This is because in the mechanical model assumption, the effect of bolts on surrounding rock is simplified as a distributed force model, which is equivalent to applying lining support to the inner wall of tunnel. Therefore, theoretically, it strengthens the supporting effect of bolts. It is because of this strengthening effect that the stress released by surrounding rock during the rheological process is reduced, resulting in a larger absolute convergence value of radial stress in the anchored zone. In the numerical simulation process, the bolt only acts between two nodes, and its support range for surrounding rock is limited, resulting in a larger stress released by surrounding rock during the rheological process and a smaller absolute convergence value. In addition, there is no sudden change at the junction between the anchored and non-anchored zones. This is because discrete jumps cannot be realized in finite differential methods like discrete element method, and there is no stress jump phenomenon.

In addition, the difference in radial stress contrast during seepage flow from outside to inside tunnel is smaller than that during seepage flow from inside to outside tunnel. The is because the seepage boundary used in this numerical simulation is a square boundary, which can lead to an increase in seepage flow volume during axisymmetric circular seepage flow. This can also have a certain impact on the rheological process of radial stress.

In Fig. 20b, the variation laws of the analytical and numerical solutions of the bolt axial force with the same bolt parameters are the same, and their evolution curves are basically consistent, which proves the correctness of the analytical solution. The absolute value of the bolt axial force increases first, then releases, and finally converges as time passes. The absolute value of the bolt axial force of the numerical solution is slightly smaller than that of the analytical solution, and this difference is more obvious at the initial stage. Due to the differences in the action mode of the bolts in the numerical model and the theoretical model, the axial deformation of the bolts in the numerical simulation process is not large. This results in a smaller numerical value, which, however, is closer to the actual value.

Judging from the above comparative analysis, it is inevitable that there will be some errors between theoretical and numerical solutions, which is mainly attributed to the differences between theoretical and numerical models. Although the “distributed force” model of the bolt can evenly distribute the bolt force to the surrounding rock and well reflect the synergistic effect of the bolt and the surrounding rock, the bolt and the surrounding rock are bonded together by the anchoring agent, the anchorage has a certain range, and the length of the bolt does not represent the length of the free section of the bolt. Therefore, the theoretical model only considers the bolt axial force, without the consideration of the shear force existing in the anchorage section. In addition, pore water pressure is divided into dynamic water pressure and static water pressure, and the changes in the rheological process of tunnel wall displacement caused by static water pressure in the anchorage supporting system are not considered. Through comparative analysis between analytical and numerical solutions of tunnel wall displacement, as shown in Fig. 20c, it is found that the absolute value of simulated tunnel wall displacement is greater than the analytical result, and the numerical value is also greater than that of the analytical solution in the curve of the absolute value of surrounding rock displacement as a function of radial distance. In the numerical simulation, the bolt axial stress is less than the stress in the theoretical model, and the weakening of the supporting effect results in an increase in the absolute value of the surrounding rock displacement.

Conclusions

This study integrates theoretical analysis and numerical simulation methods to establish a rheological mechanical model for the tunnel anchoring support system under seepage effects. It reveals the time evolution laws of the stress field, displacement field of the surrounding rock, and the stress field of the bolts. The main conclusions obtained are as follows:

  1. (1)

    In the elastic state, the stress mutation value at the boundary between the anchored and non-anchored zones is related to the seepage state. When seepage occurs inwardly, the mutation value increases, and conversely, it decreases. Additionally, when the seepage direction points towards the tunnel interior, the final trend of the normal stress shifts towards an increase in the absolute value of the stress, resulting in increased surrounding rock displacement; conversely, it shifts towards a decrease, leading to reduced surrounding rock displacement.

  2. (2)

    The presence of seepage forces makes the influence of the elastic modulus, Poisson’s ratio of the surrounding rock, and the elastic modulus of the bolts on the stress of the surrounding rock more significant. When the initial in-situ stress increases to a certain extent, the seepage forces can be ignored. Moreover, when the seepage direction points towards the tunnel interior, it intensifies the influence of the initial in-situ stress, elastic modulus, and Poisson’s ratio of the surrounding rock on the tunnel wall displacement, while weakening the influence of Poisson’s ratio on the displacement at the far end.

  3. (3)

    The rheological process of radial stress on the tunnel wall, in the anchored and non-anchored zones, and the axial force of the bolts can be divided into three stages: instantaneous increase, slow release, and convergence. The absolute value of the normal stress of the surrounding rock and its rheological process are positively correlated with the prestressing force, length, and diameter of the bolts, and negatively correlated with the spacing between bolts. The stress mutation value at the boundary is positively correlated with the prestressing force and diameter of the bolts, and negatively correlated with the length and spacing of the bolts. η1 affects the axial force of the bolts in the first stage and is negatively correlated with it; η2 affects the axial force in the second and third stages and is also negatively correlated with it.

  4. (4)

    The rheological process of tunnel wall displacement can be divided into three stages: initial stability, mid-term increase, and late-stage convergence. The absolute value of surrounding rock displacement and its rheological process are negatively correlated with the prestressing force, length, and diameter of the bolts, and positively correlated with the spacing between bolts. When the seepage direction points towards the tunnel interior, the absolute value of surrounding rock displacement increases, and conversely, it decreases. The greater the seepage force, the larger η2 becomes, and the smaller the absolute value of the displacement.

  5. (5)

    The analytical and numerical solutions show good agreement. The numerical solutions for the radial stress, tangential stress of the surrounding rock, and axial force of the anchor bolts are all smaller than the analytical solutions, while the numerical solution for the absolute value of surrounding rock displacement is slightly larger than the analytical solution. In future research, further testing and optimising the theoretical models and simulation methods through engineering methods is an important development direction to improve the accuracy of theory and simulation.