The traditional Nishihara model (Nishihara, 1957; Zhou et al., 2010) is composed of a Hooke body, Kelvin body, and Bingham body, as shown in Figure 1.

The viscosity coefficient η 2 in the traditional Nishihara model is replaced by the nonlinear viscoplastic coefficient η ( t ) related to time, so the model becomes an unsteady Nishihara model. This study divides the model into two parts: the generalized Kelvin model and the nonlinear Bingham model. The viscoelastic and viscoplastic regions of the surrounding rock are analyzed. The improved Nishihara model is shown in Figure 2:where η 1 , η 2 , and η ( t ) represent the coefficients of viscosity. E 1 and E 2 indicate the elastic modulus. σ s stands for the ultimate frictional resistance of St. Venant’s body.

In this study, the calculation model assumes that the surrounding rock is a continuous isotropic rock mass mechanics model and does not consider the action of groundwater, and when the lateral pressure coefficient is λ(t), the viscoelastic plastic problem of stress and deformation of surrounding rock of a circular tunnel is assumed, and the original rock stress field is axisymmetric along the vertical direction. The viscoelastic–plastic surrounding rock mechanical model is shown in Figure 3. The radius of excavated tunnels is R₁ , R₁-R₂ , and R₂-R₃ , representing ranges of the crack zone and plastic zone, while the elastic deformation zone should be from R₃ to infinity in theory.

In addition, many scholars have made theoretical assumptions and derivations based on experiments for the research on the expression of the nonlinear viscosity coefficient (Xiong et al., 2010; Yan et al., 2010). However, the results do not have universal adaptability. According to the nonlinear viscosity coefficient characteristics, this experiment is firstly expressed as Eq. 1: (whether it is a reasonable hypothesis, rationality will be verified in the sixth part of this study.) η ( t ) = t A + e − t B (1) where A and B are constants related to the surrounding rock properties, which can be obtained from experimental data.

Based on the derivation method of the tunnel lining displacement in the Zhao et al. (2016), the theory and the relationship between “spatial effects” were improved by combining LDP curves, in which Hoek’s formula (Nishihara, 1957) was also introduced to conduct a comprehensive analysis and theoretical derivation of “space-time effects.”

The generalized Kelvin body analyzes the viscoelastic deformation. Considering that the shape changes only when the geometry undergoes elastoplastic deformation, it can be calculated because the volume deformation is almost negligible. According to the elastic mechanics, in the cylindrical coordinate system ( r , θ , z ) , due to the axial symmetry of the tunnel calculation model, all stress components are independent of the coordinate θ , and the stress-strain is only a function of r, u. Therefore, the relationship between strain and displacement along the ( r , θ ) plane is as follows: ε r = ∂ u ∂ r , ε θ = u r (2) ε θ + ε r = u r + ∂ u ∂ r = 0 (3) then, u = A ( t ) r (4) where A ( t ) is a function of time, and substituting Eq. 4 into Eq. 2 yields: ε r = ∂ u ∂ r = − A ( t ) r , ε θ = u r = A ( t ) r 2 (5)

The generalized Kelvin constitutive equation is transformed into a stress bias form, and the boundary condition, σ r = σ θ = p , σ ¯ = p when r → ∞ , is considered, and the axisymmetric condition can be expressed as follows: ∂ ∂ t ( σ r − p ) + E 1 + E 2 η 1 ( σ r − p ) = − E 1 ∂ A ( t ) r 2 ∂ t − E 1 E 2 A ( t ) r 2 η 1 (6) where σ r , σ θ , σ ¯ —radial, circumferential, average stress;

ε r , ε θ —radial strain and circumferential strain;

For the periphery of the tunnel, the “equivalent initial geo-stress” p ( x ) corresponding to the displacement was introduced and substituted into Eq. 6 in the form of p ( x ) = p 0 λ ( t ) . Considering the boundary condition of r = r 0 , σ = 0 , the following relationship can be obtained: p 0 d λ ( t ) d t + E 1 + E 2 η 1 p 0 λ ( t ) = E 1 d A ( t ) r 0 2 d t − E 1 E 2 A ( t ) r 0 2 η 1 (7) where p 0 is geo-stress.

To facilitate the study of the “space effects” of the working surface, we first consider the effects of elastic media. The “space effect” of the excavation surface is directly reflected in the radial displacement of the tunnel along the tunneling direction. This distribution curve is the longitudinal deformation profile curve, which is regarded as the LDP curve. At present, many scholars have made in-depth studies on LDP curves within the scope of “space effect.” Panet (1995) summarized the following approximate expressions by applying the finite-element analysis method of elasticity theory and combining the relationship between the measured radial displacement around the tunnel and the distance between the driving faces: u m u m ∞ = 0.25 + 0.75 [ 1 - ( 0.75 0.75 + x / R ) 2 ] (8) where R is the radius of the tunnel cave and x is the distance between the calculation section and the excavation face.

Hoek (1999) analyzed and fitted the field monitoring data of the Mingtam underground cavern project and summarized the following empirical expressions by using the related theory analysis method: u m u m ∞ = [ 1 + e ( − x 1.1 R ) ] − 1.7 (9)

For the applicability of the above two formulas, Sun (2007) thinks that the Panet method overestimates the radial displacement of the tunnel, making it is easy to underestimate the supporting load in practical engineering and inducing the design unsafety. The result of the Hoek method agrees with the numerical analysis in the tunnel vault, waist, sidewall, and other parts, which shows that Hoek’s formula has stronger applicability and can better describe the “space effect” of the tunnel excavation surface. The radial displacement described by Hoek’s formula varies with the excavation surface, as shown in Figure 4.

Therefore, this study transformed Hoek’s empirical formula:

Let χ = 10 x / ln [ 1 + e − x 1.1 R [ ( 1 + e − x 1.1 R ) 1.7 − 1 ] 10 ] , substitution Eq. 9:

The relationship between the radial displacement u m measured around the tunnel and the cross-section distance x measured at a distance from the working face is as follows: u m = u m ∞ ( 1 + e − x χ ) (10)

There is a relationship between u m ∞ and constants χ as shown in Figure 5 below:

Since the radial free displacement u 0 of surrounding rock had been generated before the excavation entry surface reached the cross-section of the survey point, the total displacement should be calculated as follows: u = u m ∞ + u 0 = ( u ∞ − u 0 ) ( 1 + e − x χ ) + u 0 (11)

If the equivalent initial geo-stress is p ( x ) , the relationship between displacement and equivalent initial geo-stress can be expressed as follows: u = 1 K g p ( x ) (12)

Therefore, the above displacement should be proportional to the initial geo-stress, and Eq. 12 should be rewritten as follows: u = 1 K g p 0 λ ( x ) , p 0 λ ( x ) = p 0 ( x ) (13) where T is a constant relating to the characteristics of surrounding rock and the radius of the tunnel; Y is dimensionless.

Because u ∞ = p 0 K g , then u = u ∞ λ ( x ) , after substituting Eq. 13 into Eq. 11: λ ( x ) = ( 1 − u 0 u ∞ ) ( 1 + e − x χ ) + u 0 u ∞ (14)

Let λ 0 = u 0 u ∞ , as shown by Eq. 9, when x / R = 0 , then λ 0 = 0.3078 . Then Eq. 14 can be simplified as: λ ( x ) = ( 1 − λ 0 ) ( 1 + e − x χ ) + λ 0 (15)

According to Hoek’s empirical formula, the value of λ 0 is 0.3078. According to the symmetry and continuity of Figure 5, we can see that: λ ( x ) = λ 0 e − x χ 1 , χ 1 = λ 0 1 − λ 0 χ (16)

Considering the influence of rock rheology, which has been neglected before, the Hoek formula is used to comprehensively analyze the “space effect” caused by the working face’s excavation and the “time effect” caused by the viscous effect of the surrounding rock. Assuming that the tunnel is driven forward at a continuous and uniform speed, and the driving speed is υ , the simultaneous Eqs 15, 16 show that: λ e ( t ) = { 1 + ( 1 − λ 0 ) e − t T a , t > 0 λ 0 e − t T a 1 , t ≤ 0 T a = χ υ , T a 1 = χ 1 υ , (17)

When t > 0 , the above λ e ( t ) expression is introduced into Eq. 7 to solve differential equations. A ( t ) = − E 2 p 0 r 0 2 T a λ 0 + E 2 p 0 r 0 2 T a − p 0 r 0 2 λ 0 η 1 + p 0 r 0 2 η 1 E 1 ( T a E 2 − η 1 ) e − t T a − p 0 r 0 2 T a λ 0 + p 0 r 0 2 T a E 2 T a − η 1 e − t T a + p 0 r 0 2 E 1 + p 0 r 0 2 E 2 + C e − E 2 η 1 t (18)

From the above formula, it can be seen that the displacement around the tunnel is u e = A ( t ) r ; by taking into account r = r 0 , Eq. 18 can be transformed into: u e = A ( t ) r 0 = − E 2 p 0 r 0 T a λ 0 + E 2 p 0 r 0 T a − p 0 r 0 λ 0 η 1 + p 0 r 0 η 1 E 1 ( T a E 2 − η 1 ) e − t T a − p 0 r 0 T a λ 0 + p 0 r 0 T a E 2 T a − η 1 e − t T a + p 0 r 0 E 1 + p 0 r 0 E 2 + C r 0 e − E 2 η 1 t (19)

According to Eq. 19, when t → + ∞ , the limit displacement around the tunnel is u ∞ e = ( E 1 + E 2 ) p 0 r 0 E 1 E 2 by calculating the upper limit. During the tunnel’s excavation, the surrounding rock in front of the working face is affected by the construction disturbance, promoting the redistribution of the stress in the surrounding rock, thus causing the initial displacement u 0 of the surrounding rock ahead.

To study the initial displacement, we should first discuss the case when t ≤ 0 .

Equation 17 shows that when t ≤ 0 , λ e ( t ) = λ 0 e − t T a 1 ; when t → − ∞ , A ( t ) = 0 , so we can obtain: A ( t ) = u e ( 1 − λ 0 ) t T a λ 0 (20)

Because u = [ η 1 + ( E 1 + E 2 ) T a 1 ] E 1 ( E 2 T a 1 + η 1 ) p 0 λ 0 r 0 2 , when t = 0, there is: u 0 = u r 0 = [ η 1 + ( E 1 + E 2 ) T a 1 ] E 1 ( E 2 T a 1 + η 1 ) p 0 λ 0 r 0 (21)

Equation 21 shows that excavation speed directly affects the initial displacement u 0 . If the tunneling speed is plodding, the viscous stress and strain of the surrounding rock will have enough time to be fully released. So, it can be concluded that as υ → 0 , T a 1 = χ υ → ∞ . Then transform Eq. 21, then u 0 = η T a 1 + E 1 + E 2 E 1 E 2 + E 1 η T a 1 p 0 λ 0 r 0 , When T a 1 → ∞ , u 0 → E 1 + E 2 E 1 E 2 p 0 λ 0 r 0 .

On the contrary, if the tunneling face is speedy, it is instantaneously annoying from infinity to the study’s monitoring section. The viscous stress-strain of the surrounding rock does not have sufficient time to release completely, and only the elastic deformation is emancipated. At this point, let us say: υ → ∞ , then T a 1 = χ υ → 0 , so there is u 0 → p 0 λ 0 r 0 E 1 . According to the above analysis, the initial displacement u 0 of the surrounding rock expressed by Eq. 20 is between the two extremes, namely: p 0 λ 0 r 0 E 1 < u 0 < E 1 + E 2 E 1 E 2 p 0 λ 0 r 0 (22)

Let E ∞ = E 1 E 2 E 1 + E 2 , the upper form is transformed into: p 0 λ 0 r 0 E 1 < u 0 < p 0 λ 0 r 0 E ∞ (23)

Therefore, the viscoelastic deformation of the tunnel is S e = 2 ( u e − u 0 ) .

In the process of tunnel excavation, besides the instantaneous elastic deformation released from the working face, the secondary stress generated by the stress redistribution of surrounding rock exceeds the yielding stress of the rock mass at that point, which results in the viscoplastic state that the deformation of surrounding rock increases with time. The “time effect” of rock rheology is related to viscoplastic flow in surrounding rock mass (Sun, 2007). This kind of viscous flow usually occurs in the weak surrounding rock with high stress, especially for the surrounding rock with well-developed joints and fissures and dissolution. It is also found that the nonlinear Bingham model is suitable for the rheological study of this kind of surrounding rock.

Considering that only when geo-materials undergo plastic deformation and volume deformation, the change of shape can be neglected, it is assumed that the volume of the plastic zone is constant, that is, the volume strain is zero. Then there are: ε θ + ε r = u p r + d u p d r = 0 (24) thus, u p = A p ( t ) r (25)

So the Bingham viscoplastic constitutive equation can be rewritten into the form of stress deviation. Furthermore, the boundary conditions are considered when r → ∞ , σ r = σ θ = p , σ ¯ = p and the axisymmetric state can be expressed as follows:

when σ < σ s , { σ r − σ ¯ = − E ( ε r − ε ) σ θ − σ ¯ = − E ( ε θ − ε ) (26)

Because there is no creep and stress relaxation in this case, we only analyze the second case

When σ > σ s , { ∂ ( ε r − ε ¯ ) ∂ t = ∂ ( σ r − σ ¯ ) E ∂ t + ( σ r − σ ¯ ) − σ s η ( t ) ∂ ( ε θ − ε ¯ ) ∂ t = ∂ ( σ θ − σ ¯ ) E ∂ t + ( σ θ − σ ¯ ) − σ s η ( t ) (27)

When r = r 0 , the σ r = 0 , σ ¯ = p , ε r = − A p ( t ) r 0 2 , ε θ = A p ( t ) r 0 2 and axisymmetric conditions are substituted into Eq. 27, so p is calculated in the form of p = p ( x ) = p 0 λ p ( t ) : d A p ( t ) r 0 2 d t = p 0 d λ p ( t ) E d t + p 0 λ p ( t ) - σ s η ( t ) (28)

When t < 0 , the tunnel has not been excavated to the monitoring section, and the viscoplastic deformation of the surrounding rock of the section is minimal, which can be almost ignored. Therefore, it is only necessary to study the deformation after excavation. when t > 0 , λ ( t ) = 1 + ( 1 − λ 0 ) e − t T a is substituted into Eq. 28 to get: d A p ( t ) r 0 2 d t + p 0 ( 1 − λ 0 ) e − t T a E T a − p 0 [ 1 + ( 1 − λ 0 ) e − t T a ] + σ s η ( t ) = 0 (29)

The above Eq. 29 is solved to obtain the expression A p ( t ) as follows: A p ( t ) = p 0 T a λ 0 r 0 2 e − t T a η ( t ) − p 0 λ 0 r 0 2 e − t T a E − p 0 T a r 0 2 e − t T a η ( t ) + p 0 r 0 2 e − t T a E + r 0 2 ( p 0 + σ s ) η ( t ) t (30) η ( t ) = t A + e − t B is substituted into Eq. 30 can be obtained: u p = p 0 T a λ 0 r 0 e − t T a t ( A + e − t B ) − p 0 λ 0 r 0 e − t T a E − p 0 T a r 0 e − t T a t ( A + e − t B ) + p 0 r 0 e − t T a E + r 0 ( p 0 + σ s ) ( A + e − t B ) (31)

When t → + ∞ , u ∞ p = A p ( t ) / r 0 = r 0 ( p 0 + σ s ) η ( t ) t = A r 0 ( p 0 + σ s ) .

Therefore, the viscoplastic deformation of the tunnel is S p = 2 u p .

The viscoelastic theoretical solution of surrounding rock deformation mentioned above can be considered as the stress redistribution of surrounding rock after tunnel excavation is completed instantaneously with the elastic or elastic-plastic wave’s propagation velocity. The viscoelastic convergence expression of surrounding rock is S e = 2 ( u e + u 0 ) , and the viscoplastic convergence expression is S p = 2 u p . Therefore, the analytical formula of viscoelastic–viscoplastic convergence of surrounding rock is developed as follows: { S = 2 ( u e + u p + u 0 ) u e = − E 2 p 0 r 0 T a λ 0 + E 2 p 0 r 0 T a − p 0 r 0 λ 0 η 1 + p 0 r 0 η 1 E 1 ( T a E 2 − η 1 ) e − t T a − p 0 r 0 T a λ 0 + p 0 r 0 T a E 2 T a − η 1 e − t T a + ( E 2 + E 1 ) p 0 r 0 E 1 E 2 + C r 0 e − E 2 η 1 t u p = p 0 T a λ 0 r 0 e − t T a t ( A + e − t B ) − p 0 λ 0 r 0 e − t T a E − p 0 T a r 0 e − t T a t ( A + e − t B ) + p 0 r 0 e − t T a E + r 0 ( p 0 + σ s ) ( A + e − t B ) u ∞ e = ( E 1 + E 2 ) p 0 r 0 E 1 E 2 , u ∞ p = A r 0 ( p 0 + σ s ) , p 0 λ 0 r 0 E 1 < u 0 < p 0 λ 0 r 0 E ∞ (32) where σ s is the yield stress; E ∞ = E 1 E 2 E 1 + E 2 ; η 1 is the viscosity coefficients; A and B are constants relating to the properties of the surrounding rock, and C is a constant.

In Eq. 32, excavation speed υ cannot approach infinity in tunnel engineering; the elastic deformation cannot be completed instantaneously in actual engineering, and the initial deformation u 0 can hardly be measured. Therefore, in order to make the calculated value more in line with the field monitored data, the above expression is simplified as S s = 2 ( u e + u p ) for calculation.

The Dingxi tunnel of the Wujing Expressway in Hunan Province is located in Huaihua City. The tunnel area is in a low mountainous landform with large terrain fluctuations, where the midline elevation is 510–648 m and the maximum height difference is about 148 m, as shown in Figure 6. The starting pile number of the left tunnel is ZK66 + 940-ZK68 + 066, with a length of 1,126 m, while the right tunnel is K66 + 926-K67 + 999 with a length of 1,073 m.

In this case, three typical sections of ZK67 + 220 (Grade V), ZK67 + 500 (Grade VI), and ZK67 + 900 (Grade V) on the left line of the Dingxi tunnel are selected for research. The detailed geological evaluation is as follows:

(1) ZK67 + 220 section: the surrounding rock is the intensely weathered sandy slate, with a small amount of moderately weathered sandy slate. The joints and fissures of the rock mass are well developed, and the rock mass is broken.

Because the hydraulic fracturing method is simple and easy to operate in the testing process without needing to drill rock core samples or using exact electronic instruments, the testing depth can reach more than 5,000 m. Therefore, this study applies it to the initial in situ stress test of tunnel-surrounding rock.

The principle of the hydraulic fracturing method is to make the excellent wall break by increasing the water pressure based on the designed measurement depth of the hydraulic pump, then to determine the pressure and fracture direction of the characteristic point in the fracturing process, and calculate the initial stress state of the rock mass at the measurement point. Hydraulic fracturing is a two-dimensional testing method that can determine the maximum principal stress’s magnitude and direction and the minimum principal stress perpendicular from the borehole plane. Figure 8 illustrates hydraulic fracturing test device.

According to the geological engineering report of the Dingxi tunnel, the in situ stress of three representative sections, ZK67 + 220, ZK67 + 500, and ZK67 + 900 on the left line, was measured with the hydraulic fracturing method. The self-weight stress value was calculated according to the designed thickness of overlying strata ( σ v = ρ g h , rock bulk density ρ = 2.7 g / c m 3 ). The test results are shown in Table 2.

Because the surrounding rock in this study area is intensely weathered sandy slate, it is suitable for the viscoelastic–viscoplastic rheological model. The mechanical parameters of the surrounding rock of the Dingxi tunnel are shown in Table 3. The model parameters in Table 3 were derived from the creep data under the corresponding stress by the formula E 1 = σ 0 ε c ( ∞ ) , E 2 = 1 ε c ( ∞ ) + ε c i r ( ∞ ) σ − 1 E 1 and η = E 1 t ln ( 1 − E 1 ε c ( t ) σ 0 ) .

Because the design of the Dingxi tunnel is a three-center circular section, in order to ensure the accuracy of a theoretical calculation, it is necessary to convert the circular radius of the tunnel used in theoretical calculation into the equivalent circle radius using Eq. 33. r 0 = ( b 2 ) 2 + H 2 2 H (33)

The equivalent radius of the Dingxi tunnel section is r 0 = 7.03  m, which was calculated by substituting b = 12.68 m and H = 10.08 m for the upper formula.

The vault settlement of ZK67 + 220, ZK67 + 500, and ZK67 + 900 sections was calculated by the analytic equation group S s = 2 ( u e + u p ) and compared with the measured results.

The simultaneous Eq. 32 was used to calculate all the parameters into the equation group S s = 2 ( u e + u p ) , and the displacement expressions S s and ultimate displacement expressions S ∞ of surrounding rock deformation can be obtained as follows: S s = 2 [ − E 2 p 0 r 0 T a λ 0 + E 2 p 0 r 0 T a − p 0 r 0 λ 0 η 1 + p 0 r 0 η 1 E 1 ( T a E 2 − η 1 ) e − t T a − p 0 r 0 T a λ 0 + p 0 r 0 T a E 2 T a − η 1 e − t T a + ( E 2 + E 1 ) p 0 r 0 E 1 E 2 + C r 0 e − E 2 η 1 t + p 0 T a λ 0 r 0 e − t T a t ( A + e − t B ) − p 0 λ 0 r 0 e − t T a E − p 0 T a r 0 e − t T a t ( A + e − t B ) + p 0 r 0 e − t T a E + r 0 ( p 0 + σ s ) ( A + e − t B ) ] (34) S ∞ = 2 ( u ∞ e + u ∞ p ) = 2 ( E 1 + E 2 ) p 0 r 0 E 1 E 2 + 2 A r 0 ( p 0 + σ s ) (35)