Limit Equilibrium Models for Tunnel Face Stability in Composite Soft-Hard Strata

Abstract

The tunnel face stability in composite strata, commonly referred to as the soft upper and hard lower condition, is a critical challenge in tunnel construction. The soft–hard ratio (SA) strongly influences the limit support pressure as well as the failure mechanism experienced by a tunnel face. This study focused on the Xingye Tunnel project in the Xiangzhou District of Zhuhai City. By conducting numerical simulations, the impact of different SAs on the limit support pressure was investigated. Furthermore, a limit equilibrium model was established on the basis of the analysis of the results of numerical simulation. The findings were then compared and analyzed alongside those of relevant theoretical models. In the event of tunnel face instability of composite strata, the deformation tends to be concentrated mainly in the soft soil layer, with less noticeable deformation in the hard rock layer. The investigation of different SAs revealed a linear decrease in the limit support pressure ratio of the tunnel face in composite strata as SA decreases. The self-stability of the tunnel face was observed when SA ≤ 0.125. Moreover, the limit support pressure ratio predicted by the truncated log-spiral model (TLSM) exhibited a higher degree of agreement with the results of numerical simulation than those of other relevant models. The superiority of TLSM was mainly demonstrated in the range of SA = 0.25 to 1.0.

1. Introduction

During tunnel construction, composite strata with an uneven distribution of soft and hard layers are always present. Composite strata are characterized by various geological layers within an excavation section and along the extension direction, exhibiting wide variations in soil mechanics, engineering geology, hydrogeology, and other aspects. Among the composite strata encountered in practical engineering, the most common is the combination of soft upper soil layers and hard lower rock layers, often referred to as the soft upper–hard lower condition. Scholars [1] have established specific criteria to determine the physical and mechanical parameters of such composite strata. According to these criteria, the uniaxial compressive strength of the soft upper soil should be less than 1/10th of that of the hard lower rock. The unique characteristics of composite strata with the soft upper–hard lower configuration result in considerable variations in the self-stability capacity of the surrounding rock mass. These properties present formidable challenges in terms of ensuring the stability of the tunnel face. Therefore, determining a safe range of support pressure in composite strata is crucial for maintaining tunnel face stability and minimizing the impact on the surrounding environment [2].

The tunnel face stability in shield tunneling has been explored by numerous scholars through various research methods, such as experimental tests, numerical simulations, and theoretical analyses. Experimental tests have been commonly conducted to provide visualizations of failure mechanisms and to acquire monitoring results that are applicable in practical engineering, thus serving as valuable tools in studying tunnel face stability [3,4,5,6,7,8]. Liu et al. [3] investigated the stability of a tunnel face in dry sand through 1 g large-scale model tests. The internal movement of the soil was analyzed using particle image velocimetry, and the impact of movement speed on tunnel face stability was discussed. Chen et al. [5] conducted a series of 3D large-scale model tests to investigate tunnel face stability under different cover depths in dry sand. Moreover, with the advancement of computational technology, numerical simulation has become a key tool for studying large-scale and complex structures under different operating conditions [9,10,11,12,13,14]. Ren et al. [10] investigated the stability of tunnel faces reinforced with horizontal pregrouting using the finite element method (FEM). Paternesi et al. [14] used the FEM to analyze tunnel face stability under various conditions, including both reinforced and unreinforced scenarios. A strength reduction technique was employed to evaluate the safety of the tunnel faces.

In theoretical analyses, the determination of the limit support pressure often involves the use of two methods: the limit equilibrium and limit analysis methods. In the limit analysis method, plastic theory is used, and the shape of the slides in front of the tunnel face is assumed to estimate the limit support pressure [15,16,17,18,19,20,21,22]. Leca and Dormieux [15] proposed two active failure mechanisms and one passive failure mechanism for tunnel faces in cohesive soils, based on the principles of the Mohr–Coulomb yield criterion and the associated flow rule. These models were employed to conduct upper and lower bound analyses for the tunnel faces, which enabled the determination of the minimum and maximum limit support forces. Building upon the proposed failure mechanisms by Leca and Dormieux, subsequent researchers have made notable advancements. Soubra et al. [16] improved computational accuracy by employing multiple rigid truncated cones to describe the mechanisms of both active and passive failures in tunnel faces. Within the framework of kinematic design theory, Subrin and Wong [17] introduced a three-dimensional failure mechanism applicable to cohesive soil materials. Mollon et al. [18,19] investigated the active and passive failures of circular tunnel excavation faces using spatial discretization techniques and presented two rotational failure mechanisms for active and passive failures. In the limit equilibrium method, the limit support pressure is calculated by examining the static and moment balances of individual wedges within the sliding soil [23,24,25,26,27,28,29,30,31]. Horn [23] initially developed a limit equilibrium model, based on silo theory, to analyze the failure of the tunnel face. This model assumes that the failure zone consists of two parts, namely the wedge in front of the tunnel face and the prism above the wedge. Subsequently, several scholars have made a series of improvements and optimizations to Horn’s model, taking into account the specific characteristics of practical scenarios. Jancsecz and Steiner [24] considered the arching effect of the overlying soil, while Anagnostou and Kovári [25] studied the influence of seepage. Anagnostou [26] further enhanced the wedge failure model by incorporating the impact of horizontal stresses on the stability of the excavation surface using the strip method. Broere [27] analyzed the influence of arching effects on the vertical pressures at the top of the wedge. Chen et al. [28] proposed a failure model that takes into account the relationship between the height of the silo and the depth of the tunnel. The above study provided substantial contributions to the understanding of the failure mechanisms of tunnel faces under different working conditions. However, most of these studies were based on the assumption of a homogeneous soil layer in front of the tunnel face.

Few studies based on theoretical analysis have been conducted to investigate tunnel face stability in composite strata. Senent et al. [32] extended a rotational face collapse mechanism to calculate the limit support pressure of a tunnel face in layered or stratified ground. The findings suggested that the upper weak layers of the tunnel are susceptible to partial failure. Based on the upper bound theorem of kinematic limit analysis and the nonlinear Hoek–Brown yield criteria, Man et al. [33] investigated the influence of rock layer inclination and weak interlayers on tunnel face stability. Tu et al. [34] used discretization technology to improve the understanding of the collapse mechanism of a rotating rigid body based on limit analysis to evaluate tunnel face stability in inclined layered soil. Summarizing the above research findings, it is worth noting that there have been limited investigations conducted on tunnel face stability in horizontal composite layers. Moreover, there is currently a lack of a comprehensive limit equilibrium model capable of effectively predicting the limit support pressure under various ratios of soft and hard layers. Addressing this gap with further research could provide valuable insights for the field of tunneling engineering.

Tunnel face stability in soft upper–hard lower composite strata was investigated in this study. First, the impact of different ratios of soft and hard layers, denoted as the soft-hard ratio (SA), on the limit support pressure is analyzed using numerical simulation. Second, a limit equilibrium model was established to predict the limit support pressure of the tunnel face by analyzing the patterns of the variation in the displacement and shear strain contours. Finally, the model was compared and analyzed alongside numerical simulations and relevant theoretical models.

2. Engineering Background

The Xingye Expressway Tunnel Project in Zhuhai, Guangdong Province, China, was chosen as the study object, as depicted in Figure 1. The Xingye Expressway is located in Xiangzhou District. This road segment has a design speed of 60 km/h. The total length of the expressway is approximately 23.59 km. As one of the main north–south main vertical channels in the Nine Verticals and Five Horizontals road network, it plays an important role in connecting the eastern and western parts of Zhuhai City. It is also an important component of Zhuhai City’s integration with the Guangdong–Hong Kong–Macao Greater Bay Area, as well as a crucial link between the Hong Kong–Zhuhai–Macao Bridge and the Shen–Zhong Bridge. The Xingye Expressway will promote the coordinated development of Hong Kong, Macao, Zhuhai, Zhongshan, and Shenzhen, once completed, and is vital in strengthening the communication and exchange among these five cities.

The length of the shield tunneling section is 1739 m. The tunnel was constructed using the slurry balance shield tunneling method and has an outer diameter of 15.2 m and an inner diameter of 13.9 m, with the lining being 650 mm thick. The width of tunnel ring is 2.0 m. The tunnel has a minimum depth of 9.8 m and a maximum depth of 41.3 m. Furthermore, it has a minimum curvature radius of 599.5 m, and its cross-section resembles a V shape. The geological layers through which the tunnel is being excavated are highly complex. The upper layers consist of soft soil, including silty soil and fill, whereas the lower layers consist of granite with varying degrees of weathering; the geological layer distribution is shown in Figure 2. Throughout the excavation process, the distribution of soft and hard layers ahead of the tunnel face varies. In section 1-1, soft soil is predominant. In section 2-2, the proportion of soft and hard layers is relatively balanced. However, in section 3-3, the proportion of hard rock notably increases. This combination of soft and hard layers presents challenges in maintaining tunnel face stability, which, in turn, affects construction efficiency and safety.

3. Three-Dimensional Finite Element Modeling

3.1. Proportion of Soft–Hard Layers

Tunnel construction safety is affected by the uneven distribution of soft and hard layers in composite strata. This study aimed to investigate the stability of tunnel faces when encountering different combinations of soft and hard layers. The area method was used to define the different ratios of soft and hard layers [2], as illustrated in Figure 3. Additionally, we examined the influence of seven different SAs on the stability of the tunnel face (0, 0.125, 0.25, 0.5, 0.75, 0.875, and 1). Specifically, SA = 0 indicates that the area in front of the tunnel face consists entirely of hard rock, whereas SA = 1 signifies that the area in front of the tunnel face is entirely composed of soft soil.

3.2. Constitutive Model and Parameters

The stress–strain behavior of soil layers is described by the small strain-hardening model (HSS), which reasonably mimics the nonlinear characteristics and unloading behavior of the soil and has been proven to be applicable [35,36]. The values of the relevant parameters were derived from practical engineering reports based on the constant parameter value method proposed by [37] and others. The HSS model contained 13 parameters, as listed in

Table 1. Soil parameter values for HSS.

Parameter	Soft Soil	Hard Rock	Parameter	Soft Soil	Hard Rock
Effective friction angle, φ′ (°)	11	35	Reference shear modulus, G0ref (MPa)	20	240
Dilation angle, Ψ (°)	0.1	0.1	Reference stiffness stress, Pref (kPa)	100	100
Lateral pressure coefficient, K0	0.81	0.43	Poisson’s ratio for unloading/reloading, νur	0.2	0.2
Effective cohesion, c′ (kPa)	12	32	Unloading/reloading stiffness, Eurref (MPa)	7.2	150
Failure ratio, Rf	0.9	0.9	Power for stress-level dependency of stiffness, m	0.8	0.3
Secant stiffness, E50ref (MPa)	2.4	30	Shear strain corresponding to initial shear modulus of 70%, γ0.7	10−4	10−4
Tangent stiffness, Eoedref (MPa)	2.4	30

. In addition, to model the tunnel lining, the shell structural element was used. The key properties considered for the lining included an elastic modulus of 36.0 GPa, Poisson’s ratio of 0.2, a thickness of 25 cm, and a unit weight of 26 kN/m³.

3.3. Numerical Model and Simulation Process

In this study, numerical analyses were carried out employing the commercial software PLAXIS 3D. The utilization and validation of PLAXIS 3D by numerous researchers in investigating diverse challenges associated with tunnel faces underscore its reliability and credibility [10,14,38]. The model shown in Figure 4 was constructed to analyze the stability of the tunnel face in composite strata with varying proportions of soft and hard layers. The dimensions of the model, measured in length (X axis) by width (Y axis) by height (Z axis), were 50 m × 60 m × 70 m. The tunnel had a diameter, D, of 15.2 m and was buried at a depth of C = 15.2 m. The soil was divided into two layers, with the upper layer consisting of soft soil and the lower layer composed of hard rock. In this study, the tunnel cover depth (C = 15.2 m) was kept constant while modifying the thickness of the upper soft soil layer to explore different SAs. The boundary conditions of the model were defined as follows: the top face was considered a free boundary, the four side faces were constrained in the normal direction, and the bottom face was set as a fixed boundary. Seepage effects were not considered, and the groundwater level was established at −70 m. To accurately represent the tunnel face and the upper soft soil layer, mesh refinement was applied, resulting in a total of 32,898 mesh elements and 49,456 nodes.

Shield tunneling construction is a dynamic process. This study focused on the influence of different SAs in composite strata on tunnel face stability. The tunneling process was not considered in this study. The specific simulation steps were as follows [10]:

(i)

The numerical model was established, and the initial stresses were generated using the K₀ process;

(ii)

The excavation progressed in increments of a one-time advancing length, denoted as D, while incorporating the cooling of the soil elements;

(iii)

The support pressure on the tunnel face was set to be equal to the initial ground horizontal stress in the opposite direction;

(iv)

The support pressure gradually decreased while the displacement of the soil ahead of the tunnel face increased. The simulation continued until the support pressure showed minimal changes and a substantial increase in horizontal displacement of the soil was observed. This indicated that the failure of the face had occurred, and the calculation was terminated.

4. Analysis of Numerical Results

4.1. Limit Support Pressure

Figure 5 illustrates the correlation between the limit support pressure ratio and different SAs for the tunnel face. Figure 5a shows that as the support pressure ratio decreases within the range of SA = 0.25 to 1, the maximum horizontal displacement of the tunnel face gradually reduces. However, when the support pressure ratio exceeds the critical value, the horizontal displacement suddenly increases, indicating tunnel face instability. For SA values of 0 and 0.125, no significant change occurs in the maximum horizontal displacement of the tunnel face, even when the support pressure ratio is reduced to 0. This suggests the self-stability behavior of the tunnel face.

Figure 5b displays the relationship between the limit support pressure ratio and SA. The limit support pressure ratio linearly decreases as SA decreases. When the SA is equal to or less than 0.125, the limit support pressure ratio remains constant at 0. Furthermore, a linear fit with high accuracy (R² ≈ 0.994) can be observed for SA values ranging from 0.125 to 1. These findings emphasize the strong influence of the proportion of soft soil area on the tunnel face in terms of the limit support pressure ratio in composite strata.

4.2. Deformation Analysis

Figure 6 depicts the total displacement contours of composite strata at the limit state, with the interface between the soft and hard layers represented by a white dashed line. The legend is scaled from 0 to 1.1 m. The visualization highlights that the majority of soil displacement occurs within the soft soil, whereas the hard rock experiences relatively minimal displacement. As the SA progressively increases, the maximum displacement of the tunnel face shifts upward from its center located within the soft soil layer. For instance, when SA is 1.0, the largest displacement of the thrust surface is observed at 1/4 D (upper portion). Moreover, as the SA continues to increase, the disturbance originating from tunnel face instability propagates toward the surface ground, consequently amplifying ground deformation. These observed phenomena stem from the fact that, when maintaining constant values for tunnel diameter and cover depth, a larger SA indicates a larger proportion of soft soil within the tunnel face. This, in turn, increases the challenges associated with maintaining tunnel face stability and the radius of instability propagation, thus expanding the extent of disturbance within the soil layer.

Figure 7 depicts the shear strain contours of composite strata at the limit state, with the legend range set from 0 to 0.3. Similar phenomena can be observed, where the shear strain within the strata is predominantly concentrated in the soft soil layer, whereas the shear strain of the hard rock layer is less prominent. This observation aligns with the findings shown in Figure 6. In general, tunnel face instability is primarily induced by shear failure in the soil layer. The small strain-hardening model adheres to the Mohr–Coulomb failure criterion, which indicates that strata deformation sharply increases and results in substantial displacement and instability of the tunnel face when the stress in the soil layer exceeds the shear strength. Consequently, the contours presented in both Figure 6 and Figure 7 complement each other, providing a comprehensive understanding of the phenomena.

5. Analysis of Theoretical Model

5.1. Overview

In this study, the limit equilibrium method has been adopted to predict the limit support pressure ratio of composite strata with various SAs. The Horn’s classical wedge model is predominantly composed of a wedge and an overlying prism, as depicted in Figure 8a. However, in both the traditional and modified wedge models [23,24,25], the assumed planar failure surface is questioned due to experimental tests [5] and numerical simulations [21] indicating that the actual failure surface appears to be a logarithmic spiral curve in profile. Murayama et al. were the first to predict the limit support pressure of the tunnel face using a two-dimensional (2D) logarithmic spiral model, as shown in Figure 8b. Due to the fact that the problem of tunnel excavation is three-dimensional, many subsequent studies have expanded Murayama’s model to a three-dimensional log-spiral model [10,39].

Furthermore, the failure mechanism of tunnel face in soft upper–hard lower layers has been investigated by researchers through a combination of theoretical studies and experimental tests. The relevant theoretical models have been proposed, such as the log-spiral model [40] and partial-wedge model [41], to predict the limit support pressure ratio of composite strata with various SAs, as shown in Figure 9. In this study, numerical simulations were conducted to investigate the failure mechanism at the limit state under different ratios of soft and hard layers. The relevant theoretical models were optimized based on the phenomena shown in Figure 6 and Figure 7. An optimized model is shown in Figure 10.

Figure 6 and Figure 7 show that when the tunnel face fails, slip predominantly occurs within the soft soil layer. The slip surface is depicted by incorporating a log-spiral model, whereby the initiation of slip originates from the tunnel crown and terminates at the interface between the soft soil and hard rock. This optimized model is referred to as a truncated log-spiral limit equilibrium model (TLSM). In Figure 10, a collapse mechanism is shown, similar to classical wedge models, which involves the presence of a truncated log-spiral (TLS) slip surface in front of the tunnel face and an overlying prism. Furthermore, D represents the diameter of the tunnel, and C signifies the cover depth of the tunnel. Here are some basic assumptions for the TLSM:

(i).

The failure only occurs in areas of the upper soft soil layer [40], and the geometric shape of the slip surface is a logarithmic spiral.

(ii).

The upper soft soil layer is considered to be homogeneous and isotropic, adhering strictly to the Mohr–Coulomb failure criterion on each failure surface.

(iii).

The soil pressure generated by the overlying prism acts vertically on the top of TLS.

To facilitate the implementation of this straightforward model, a square with an equivalent area was used to approximate the circular cross-section of the tunnel face [39]. This simplification is acceptable as the pressure obtained from both the experimental tests and numerical simulations reasonably corresponds to the predicted limit support pressure [5,10].

From the force analysis of the longitudinal section of the TLS model, as shown in Figure 11, we established a coordinate system with point O as the origin. The Y and Z axes are parallel to the horizontal and vertical directions, respectively. The slip surface of the log-spiral model initiates at point a, terminates at point d, and is intercepted at point b. In Figure 6, the variables R_a, R_b, and R_d denote the lengths of lines Oa, Ob, and Od, respectively. L₁ and L₂ represent the horizontal distances from point O to the tunnel face and from the tunnel face to point a, respectively. The angles φ, α, and β represent the friction angle, the angle between line R_a and R_b, and the angle between line R_a and R_d, respectively. Furthermore, we assumed that R_a forms an angle, φ, with the horizontal line, resulting in the log-spiral being perpendicular to the horizontal line at point a [39].

The equation for the log-spiral model in a polar (R, x) coordinate system can be expressed as follows:

$$R(x)=Ra\cdot e^{x\tan \varphi }$$

(1)

The following equations can be derived based on the geometric relationships shown in Figure 11:

$$R_d\cdot \sin (\varphi +\beta )-Ra\cdot \sin \varphi =D$$

(2)

Furthermore, the geometric parameters R_a, R_d, L₁, and L₂ in Figure 9 can be expressed as follows:

$$R_a=\frac{D}{e^{\beta \cdot \tan \varphi }\sin (\varphi +\beta )-\sin \varphi }$$

(3)

$$R_d=\frac{D{e}^{\beta \cdot \tan \varphi }}{e^{\beta \cdot \tan \varphi }\sin (\varphi +\beta )-\sin \varphi }$$

(4)

$$L_1=\frac{D{e}^{\beta \cdot \tan \varphi }\cos (\varphi +\beta )}{e^{\beta \cdot \tan \varphi }\sin (\varphi +\beta )-\sin \varphi }$$

(5)

$$L_2=\frac{D\cos \varphi -D{e}^{\beta \cdot \tan \varphi }\cos (\varphi +\beta )}{e^{\beta \cdot \tan \varphi }\sin (\varphi +\beta )-\sin \varphi }$$

(6)

An important characteristic of the log-spiral model described by Equation (1) is that the tangent at the intersection of any radial line and the spiral forms an angle, φ, with the perpendicular direction [38].

5.2. Solving for Limit Support Pressure

The limit support pressure ratio of the tunnel face is mainly calculated based on the moment equilibrium at origin O (as shown in Figure 11), and the limit equilibrium equation is expressed as follows:

$$M_v+M_w-M_{NT}-2M_{TS}-M_p-M_b=0$$

(7)

where M_v represents the moment of the vertical earth pressure, σ_v; M_w represents the moment of the weight, W, in the log-spiral zone; M_NT represents the moment of the resistance to sliding on the log-spiral slip surface; M_TS represents the moment of the resistance to shear force, T_s, on the lateral sliding surface; M_p represents the moment of the support pressure, P, on the tunnel face; and M_b represents the moment of the resistance force on the bottom slip surface of the log-spiral zone.

(1)

Calculation of M_v

In scenarios where the thickness of the overburden soil is less than (or equal to) the diameter of tunnel and the ground exhibits relatively low stiffness, the soil arching effect is typically neglected, and the vertical earth pressure is determined using the full-overburden theory [38]. Hence, the vertical earth pressure, σ_v, exerted on the top of the tunnel can be described using the following equation:

$${\sigma }_v={\sigma }_s+\gamma C$$

(8)

where σ_s denotes the additional pressure of the surface; γ denotes soil density.

Moment, M_v, can be expressed as follows:

$$M_{\mathrm{v}}={\sigma }_{v}B{L}_2(R_a\cos \varphi -L_2/2)$$

(9)

$$B=\frac{\pi D}{4}$$

(10)

where B represents the equivalent width of the tunnel face.

(2)

Calculation of M_w

The forces acting on the ith unit obtained from the log-spiral zone are illustrated in Figure 12. The weight of the ith unit can be expressed as:

$$dW=[R(x)\cos (\varphi +x)-L_1]\gamma Bdz$$

(11)

where

$$dz=R(x)\sin \theta dx$$

(12)

$$\theta =\frac{\pi }{2}-x$$

(13)

and θ represents the angle between the tangent of log-spiral slip surface and the horizontal line.

The moment of soil weight of the ith unit with respect to origin O can be expressed as:

$$d{M}_w=\frac{R(x)\cos (\varphi +x)+L_1}{2}[R(x)\cos (\varphi +x)-L_1]\gamma Bdz$$

(14)

The calculation for the moment of the soil weight of the log-spiral zone with respect to the origin O is as follows:

$$M_w={\displaystyle {\int }_0^{\alpha }\frac{R(x)\cos (\varphi +x)+L_1}{2}[R(x)\cos (\varphi +x)-L_1]\gamma BR(x)\sin (\frac{\pi }{2}-x)}dx$$

(15)

(3)

Calculation of M_NT

To obtain the moment of the normal force and shear force on the log-spiral slip surface of the ith unit with respect to origin O, the following steps can be followed:

$$d{M}_{NT}=R(x)\cos \varphi dT-R(x)\sin \varphi dN$$

(16)

The relationship between the normal force and shear force of the ith unit is denoted as:

$$dT=\tan \varphi dN+BcR(x)dx$$

(17)

By combining Equations (16) and (17), the moment of the resistance to sliding on the log-spiral slip surface with respect to origin O can be calculated as follows:

$$M_{NT}={\displaystyle {\int }_0^{\alpha }cB{R}^2(x)}\cos \varphi dx$$

(18)

(4)

Calculation of M_TS

Assuming a linear distribution of vertical earth pressure along the Z axis on two vertical slip surfaces, the shear force of the ith unit can be expressed as follows:

$$d{T}_s=[K_0({\sigma }_v+\gamma z)\tan \varphi +c][R(x)\cos (\varphi +x)-L_2]dz$$

(19)

where

$$z=R(x)\sin (\varphi +x)-R_a\sin \varphi$$

(20)

$$K_0=1-\sin \varphi$$

(21)

The moment of the shear force on the ith unit vertical slip surface with respect to the origin O can be expressed as follows:

$${dM}_{TS}=[K_0({\sigma }_v+\gamma z)\tan \varphi +c][R(x)\cos (\varphi +x)-L_2]dz$$

(22)

Then, the moment of the resistance to shear force, T_s, on the lateral sliding surface with respect to origin O can be calculated as follows:

$$M_{TS}={\displaystyle {\int }_0^{\alpha }[K_0({\sigma }_v+\gamma z)\tan \varphi +c][R(x)\cos (\varphi +x)-L_2]R(x)\sin \theta }dx$$

(23)

(5)

Calculation of M_b

According to the vertical equilibrium condition of the TLS, as shown in Figure 9, the following equation can be obtained:

$$N_b=B{L}_2{\sigma }_v+W-(T+T_{\mathrm{s}})\sin \theta -N\cos \theta$$

(24)

The relationship between N_b and T_b is denoted as:

$$T_b=N_b\tan \varphi +cBl$$

(25)

where

$$l=R_b\cos (\varphi +\alpha )-L_2$$

(26)

The moment of the resistance force on the bottom slip surface of the log-spiral with respect to the origin O can be expressed as follows:

$$M_{\mathrm{b}}=\frac{R_b\cos (\varphi +\alpha )+L_2}{2}{T}_b-R_b\sin (\varphi +\alpha )N_b$$

(27)

(6)

Calculation of σ_p

According to the geometric relationship in Figure 9, the moments of P are calculated as follows:

$$M_p=[R_d\sin (\varphi +\alpha )-\frac{D}{2}]P$$

(28)

By substituting Equations (9), (15), (18), (23) and (27) into Equation (7), the limit support pressure, σ_p, can be calculated as follows:

$${\sigma }_{p,c}=\frac{M_v+M_w-M_{NT}-2M_{TS}-M_b}{[R_d\sin (\varphi +\alpha )-\frac{D}{2}]BD}$$

(29)

It is worth noting that within the model solving process, the independent variables include the mechanical parameters of the soil layer (cohesion, friction angle, and soil weight) and the basic parameters of the tunnel (tunnel burial depth and diameter). The remaining variables are dependent variables. Generally, the independent variables mentioned above have a profound influence on the final results attained [39].

5.3. Comparison and Analysis

The limit support pressure ratios obtained from the proposed model (TLSM), numerical simulation (NS), log-spiral model (LSM) [40], and partial-wedge model (PWM) [41] are shown in Figure 13. A general similarity can be observed in the trends of the changes in the four curves, with a gradual decrease in the limit support pressure ratio as SA reduces. However, substantial deviations can be noted between the predictions of the partial-wedge model and those of the numerical simulations. Notably, when SA ranges from 0.25 to 1.0, the limit support pressure ratio is consistently underestimated by the partial-wedge model, resulting in a serious risk of face instability. Similarly, as SA decreases within this range, the log-spiral model exhibits an increasing disparity from the numerical results, thereby indicating its relatively unreliable predictions. In contrast, the proposed TLSM demonstrates a close alignment of its predictions with the numerical simulations. Furthermore, the self-stability of the tunnel face occurs at SA ≤ 0.125.

As shown in

Table 2. Comparison between the theoretical and numerical results.

SA	0	0.125	0.25	0.5	0.75	0.875	1
NS	0	0	0.07	0.20	0.28	0.33	0.37
TLSM	0	0	0.08	0.22	0.28	0.31	0.39
LSM	0	0	0	0.14	0.26	0.31	0.39
PWM	0	0	0	0.07	0.15	0.18	0.23

, seven cases involving SA = 0, 0.125, 0.25, 0.5, 0.75, 0.875, and 1, with corresponding theoretical and numerical simulations were analyzed. When SA = 0 and 0.125, the results predicted by the three models are consistent with the numerical simulation. When SA ranges between 0.25 and 1.0, the results predicted by TLSM are closer to the numerical simulation compared to LSM and PWM. In addition, the difference between the predicted results of TLSM and the numerical simulation is controlled at around 10%. Through a comprehensive comparison of the different models, it was found that the superiority of the proposed model is mainly demonstrated in the range of SA = 0.25 to 1.0. The fitting of various theoretical models is depicted in Figure 14. In Figure 14a, the X axis corresponds to the limit support stress ratio predicted by TLSM for different SAs, while the Y axis represents the numerical simulation results. The detailed values can be found in Table 2. The coordinate meanings of Figure 14b,c are analogous to Figure 14a. The prediction outcomes of the three models were subjected to a proportional function (y = x) fitting. The TLSM exhibited an impressive fitting accuracy (R²) of 0.991, followed by the LSM with a fitting accuracy of 0.934, and the PWM with a comparatively lower fitting accuracy of 0.446. These results affirm that the proposed model (TLSM) outperforms both LSM and PWM in accurately predicting different SAs.

6. Conclusions

In this study, a series of three-dimensional numerical analyses was conducted to investigate the failure modes of composite strata with varying soft–hard ratios during tunnel face failure. Furthermore, the limit support pressure ratios of the tunnel face were determined by employing the limit equilibrium method, and the results were compared with those of relevant theoretical models. The following conclusions were derived:

(1)

The investigation of different SAs revealed a linear decrease in the limit support pressure ratio of the tunnel face in composite strata as SA decreases. Furthermore, the self-stability of the tunnel face was observed when SA was less than or equal to 0.125.

(2)

In the event of tunnel face instability, displacement predominantly occurred in the soft soil layer, whereas the deformation of the hard rock layer was relatively minimal. Additionally, an increase in the SA led to increased overall strata deformation and subsequently increased surface ground disturbance.

(3)

The comparative analysis of different theoretical models for predicting the limit support pressure ratio reveals that the TLSM exhibits predictions that closely match those obtained from numerical simulation, with a difference of approximately 10%.

(4)

The prediction results of the three models were subjected to a proportional function (y = x) fitting. The fitting accuracy (R²) for the three models (TLSM, LSM, and PWM) were 0.991, 0.934, and 0.446 respectively. This indicates that the proposed model (TLSM) outperforms both LSM and PWM in accurately predicting different SAs.