Advances in Civil Engineering
Volume 2021 (2021), Article ID 3591670, 14 pages
https://doi.org/10.1155/2021/3591670

Research Article

Limit Analysis of Collapse Mechanisms for Tunnel Roofs Subjected to Pore Water Pressure: A Numerical Approach

Jing-jing Liu,1,2 Tie-lin Chen,1,2 Chang-ling Xie,³ Jian-hua Tian,⁴ and Yu-xin Wei⁴

¹Key Laboratory of Urban Underground Engineering of Ministry of Education, Beijing Jiaotong University, Beijing 100044, China
²School of Civil Engineering, Beijing Jiaotong University, Beijing 100044, China
³Zhejiang Huazhan Institute of Engineering Research & Design Co., Ltd., Ningbo 315010, China
⁴Beijing Urban Rapid Transit Development Co., Ltd., Beijing 100027, China

Correspondence should be addressed to Tie-lin Chen; tlchen1@bjtu.edu.cn

Received 9 September 2019; Revised 10 August 2020; Accepted 4 May 2021; Published 27 May 2021

Academic Editor: Hui Yao

Copyright © 2021 Jing-jing Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The collapse mechanism of a circular unlined tunnel roof subjected to the pore water pressure under plane strain conditions is investigated in this article. First, the model of calculating the function expression of the detaching surface for the collapsing block is formed in the framework of the upper bound theorem of limit analysis and the extremum principle. The analytical solution of the pore water pressure around the tunnel in a two-dimensional steady seepage field is employed in the equations of the model. Then, the numerical approach based on the Runge–Kutta algorithm and traversal search method is proposed to solve the complex equations. The obtained expression of the detaching surface for the collapsing block provides the shape of the collapsing block and a theoretical basis for designing the support force for tunnels. The proposed limit analysis method and numerical approach are verified by comparing with existing theoretical solutions and the numerical simulation result, and they are suitable for deep, shallow tunnels and layered strata. Moreover, the effects of different parameters on the collapse mechanism are investigated, and qualitative results are provided.

1. Introduction

The collapse of a tunnel roof, which threatens the safety of tunnel construction or surface buildings, is a major challenge in tunnel engineering. The investigation of the collapse mechanism of tunnel roofs contributes to understanding the failure properties of surrounding rock masses and providing a theoretical basis for designing the support force for tunnels.

As one of the most efficient methods to analyze the collapse mechanism, the limit analysis method is more rigorous for deriving solutions because it considers the stress-strain relationship of soil mass in an idealized elastic-perfectly plastic manner and does not need to carry out the step-by-step elastic-plastic analysis [1–3]. Thus, the limit analysis method is widely used in the stability assessment of geotechnical engineering [4–6]. The typical approach of tunnel collapse analysis based on the upper bound theorem of the limit analysis method is that the predefined collapse mechanism is specified firstly, and then the critical support pressure of tunnels is obtained by choosing the variation of the dimension of the collapse mechanism [7–9]. The accuracy of an upper bound calculation depends, to a certain extent, on the closeness of the assumed mechanism to the real failure pattern [7]. Various collapse mechanisms were specified in some studies, which can be classified into translational rigid block mechanisms [10–12], rotational rigid block mechanisms [13–15], and continuous velocity fields [16–19]. Different from the existing upper bound analysis, Fraldi and Guarracino obtained the two-dimensional shape of the collapsing block for tunnel roofs directly by combining the upper bound theorem and the extremum principle, and the collapse load of tunnels is further obtained by integrating the area between the detaching surface of the collapsing block and the tunnel boundary, which provides the theoretical basis for designing the support force for tunnels [20]. Thus, the predefined collapse mechanism is not required in the calculation process for this method, and any shape of the detaching surface in the final result is possible [21]. Because this method is effective and simple to investigate the collapse mechanism of tunnel roofs, the application of the method is extended.

It is well known that groundwater is a key factor affecting tunnel roof stability. Thus, the study of the pore water pressure impacting the collapse mechanism of tunnel roofs employing the upper bound theorem of limit analysis draws a great deal of attention. By extending the study of Fraldi and Guarracino [20], Huang and Yang [22] developed a numerical solution for investigating the collapse mechanism in circular tunnel subjected to the pore water pressure, where the pore water pressure was expressed in terms of the pore pressure coefficient ru. It is a simplified result that cannot reflect the details of the pore water pressure around the tunnel. Based on this method, many researchers considered more conditions further, such as the progressive collapse mechanism [23], stratified rock [24, 25], nonassociative flow rule [26], and cover depths [27], while the pore water pressure was generally expressed by the pore pressure coefficient ru. In some other investigations, the pore water pressure distribution around the tunnel was obtained by numerical simulation method [28, 29], or their numerical simulation results were fitted to an analytical expression, which was used for stability analysis of tunnels [30–32].

In fact, the attention on estimating the analytical expression for the pore water pressure around a tunnel is growing, and many important works have been published [33–35]. The widespread approximative solution is ascribed to Goodman [36], which is based upon an approximative approach by Polubarinova-Kochina [37]. Especially, Kolymbas and Wagner [38] proposed an exact analytical solution for the distribution of hydraulic head in a steady seepage field for a tunnel with a circular cross-section. Thus, the analytical solution provides a new approach for analyzing the collapse mechanism of tunnels subjected to the pore water pressure.

This article focuses on the analysis of the collapse mechanism of a circular unlined tunnel roof subjected to pore water pressure. The study of Fraldi and Guarracino [20] is extended in this article, the upper bound theorem and extremum principle are combined to form the model of calculating the function expression of the detaching surface for the collapsing block, and the analytical solution of the pore water pressure distribution around the tunnel in two-dimensional steady seepage field is employed in the equations of the model. For the complexity of the equations, the numerical approach based on the Runge–Kutta algorithm and traversal search method is proposed to solve the equations. A comparison with existing theoretical solutions and the numerical simulation result is made to validate the proposed limit analysis method and numerical approach. The effects of parameters on the collapse mechanism are investigated. Finally, the applicability of the proposed limit analysis method and numerical approach under different cover depths or roof stratification conditions is discussed.

2. Limit Analysis of Tunnels Subjected to Pore Water Pressure

In this section, the model is formed to calculate the function expression of detaching surface of the collapsing block under the steady seepage conditions. The model contains the work equation and variational equation about the collapsing block.

The assumed collapse mechanism of tunnel roofs is that the collapsing block separated by the detaching zone moves downwards from stationary rock masses, and the detaching surface is described by an unknown function expression f (x). The plastic stress and strain rates in the detaching zone are derived from the yield criterion and associated flow rule and used to calculate the internal rate of dissipation. The external rate of work consists of the work rate of the total weight of collapsing block and the pore water pressure acting on the detaching surface. Equating the external rate of work to the internal rate of dissipation, the work equation containing the unknown function expression f (x) for the collapse mechanism is formed. In order to obtain the critical or least upper bound solution of the detaching surface, the extremum principle is employed to form the variational equation. By solving the work equation and variational equation, the function expression f (x) of the detaching surface for the collapsing block can be obtained.

The details of the work equation and variational equation of the mode are described in the following.

2.1. Work Equation Based on Upper Bound Theorem

To calculate the detaching surface of the collapsing block of the tunnel subjected to pore water pressure, the first condition is the work equation derived from the upper bound theorem in limit analysis.

In this case, an unlined tunnel with a circular cross-section in the steady seepage field is analyzed in the limit analysis. The following assumptions are introduced:

(1)The soil or rock is homogeneous and its permeability is isotropic(2)Circular cross-section tunnel with a constant hydraulic head(3)The water table is horizontal and remains constant(4)A state of steady flow

In this article, the case of the hydraulic head of ground surface  = 0 is taken as an example for analysis. The water table lies at the ground surface. Figure 1 shows the diagram of the hydraulic head around tunnel. It is noted that the datum is at the ground surface. Point A is at the ground surface, its pressure head and elevation head are zero, and its total head  = 0. Point B is in the tunnel, its pressure head is (pu is the pore water pressure, is unit weight of water), elevation head is –z (z > 0), and its total head ht = pu/ − z < 0. The head loss of flow line AB is  − ht = z − pu/.

Figure 1: Diagram of hydraulic head around tunnel.

As stated in the upper bound theorem, for any assumed collapse mechanism, if the external rate of work exceeds the internal rate of dissipation, it means the rock masses are unsteady and the external forces cannot be carried. If the external rate of work is less than the internal rate of dissipation, it means the rock mass is stable and no collapse happens. Equating the external rate of work to the internal rate of dissipation, the work equation for a collapse mechanism is formed, and it gives an unsafe upper bound solution, which means that collapse impends or has taken place. It is the judging criterion for the collapse of a tunnel. As shown in Figure 2(a), the collapse mechanism is that the collapsing block separated by the detaching zone moves downwards from stationary rock masses. An unknown function expression f (x) is used to describe the plastic detaching zone in the Cartesian reference frame x-y, which is named the detaching surface and can be calculated in the following. The detaching zone has finite thickness , and the stresses on it are shown in detail. The relative displacement and velocity of the entire collapsing block are vertical displacement and velocity in the admissible velocity field.

(a)

(b)

Figure 2: Collapse mechanism of rock mass for tunnels in steady seepage field. (a) Schematic diagram of collapse mechanism. (b) Geometric condition of detaching surface.

The Hoek–Brown criterion is adopted as the failure criterion of rock and soil mass in the limit analysis method proposed in this article because it can reflect the inherent nonlinear failure characteristics of rock and soil masses and belongs to the nonlinear strength criterion. The Hoek–Brown criterion is an empirical criterion and is put forward according to triaxial tests and scale field test values. Its Mohr envelope form expressed by normal and shear stresses is shown in equation (1) [39–41]:

where n represents the unit vector normal to the failure plane; is the normal effective stress; is the shear stress; A and B are the material constants; is the uniaxial compressive strength of the intact rock; is the tensile strength of the rock mass .

In limit analysis, the rock masses are assumed to be a perfectly plastic material and the plastic strain rates are derived from an associated flow rule corresponding to a yield function [2]. Employing the geometric condition of the detaching surface f (x) (Figure 2(b)), the plastic stress and strain rates are presented in equations (2) and (3) [20]:

where n represents the unit vector normal to the detaching surface (f (x)); and are normal and shear plastic strain rate, respectively; the prime represents the derivation of the function to its variable, i.e., .

Because the collapsing block is symmetrical with respect to the y-axis, half of the collapsing block is analyzed in the following. The internal rate of dissipation D along the detaching surface f (x) is

where S is the length of half detaching surface (f (x)); is the elementary length of the detaching surface (f (x)).

There are two choices relative to the external forces considered to act on the collapsing block [42–44]. For the first choice, a soil-water mass is considered, and the external force is a combination of the total weight of soil-water mass and the resultant of boundary water pressure; for the second choice, the soil skeleton only is considered, and the external force is a combination of the submerged weight and the seepage force. In the analysis of this article, the former is used and the external forces include the total weight of the collapsing block and pore water pressure acting on the boundary, as shown in Figure 3. Thus, the external rate of work consists of the work rate of the total weight of collapsing block and the work rate of the pore water pressure acting on the detaching surface . where L is the half of the span of the detaching surface (f (x)); γsat is the saturated unit weight of rock masses; V is the volume of the collapsing block per unit length, V = C (x) − f (x); C (x) is a known function describing the circular tunnel profile, ; pu (x, f (x)) is the pore water pressure acting on the detaching surface (f (x)); .

Figure 3: Diagram of external force on collapsing block.

Based on the above assumptions, the pore water pressure pu (x, y) in the two-dimensional seepage field for a tunnel is as follows [38]:

where is the unit weight of water; is the total head of ground surface; ht is the total head of the circumference of the tunnel;  − ht is the hydraulic head difference between the ground surface and the tunnel circumference; r is the circular tunnel radius; C is the depth that the vertical distance from the center of circular tunnel to the ground surface; , .

Let the functional Ω be equal to the external rate of work acting on the collapsing block minus the internal rate of dissipation along the detaching surface, and Ω is expressed as follows:

is a functional, which is expressed as follows:

The thickness of the plastic detaching zone (Figure 2(a)) ebb in the calculation [20]. It is worth noting that the internal rate of dissipation D is calculated to be negative in Cartesian coordinates x-y in Figure 2; thus, in equation (7), the minus is omitted in front of item D.

Therefore, the work equation of the upper bound theorem that the external rate of work is equal to the internal rate of dissipation expressed by the functional Ω is

that is,

2.2. Variational Equation Based on Extremum Principle

To calculate the detaching surface of the collapsing block of the tunnel subjected to the pore water pressure, the second condition is the variational equation derived from the extremum principle.

According to the extremum principle, in all the kinematically admissible velocity field, the actual velocity field makes the first variation of the functional (equation (7)) to be zero [45], which forms the variational equation (equation (11)), and it is the second condition for deriving the detaching surface.

Euler–Lagrange equation (equation (12)) transforms the variational equation (equation (11)) into the differential equation (equation (13)) [46]:

that is,

Thus, two conditions to calculate the detaching surface of the collapsing block of the tunnel subjected to pore water pressure derived from the upper bound theorem (equation (9)) and the extremum principle (equation (11)) are obtained as mentioned above, which are equations (10) and (13).

After the detaching surface of collapsing block f (x) is obtained through the numerical approach presented in the following, the weight of collapsing block per unit length, namely, the collapse load , can be calculated by integrating the area between the detaching surface of the collapsing block and the tunnel boundary (equation (14)), which can be used for designing the support force for tunnels.

3. Numerical Approach for Limit Analysis

The general procedures [20] for calculating the expression of the detaching surface f (x) through solving the variational equation based on the extremum principle and work equation derived from the upper bound theorem are as follows. Firstly, by integrating the second-order differential equation transformed from the variational equation, the expression f (x) with unknown parameters can be obtained. Secondly, the unknown parameters are determined by the geometric and boundary conditions and the work equation. As the pore water pressure is taken into account, the calculation of the second-order differential equation (equation (13)) is complicated, which cannot obtain the expression of detaching surface f (x) by integrating. Therefore, the numerical approach is proposed and solved by Matlab. The main steps of the process are as follows.

Step 1. Runge–Kutta algorithm [47, 48] is applied to solve the second-order differential equation (equation (13)). By assigning different initial values of variables, a group of surfaces is obtained, which are in the same shape. It is noted that the shape of the surface is determined by the second-order differential equation (equation (13)).

Step 2. Through calculating the functional Ω of each surface by equation (7), the surface satisfying the work equation (equation (9)) is the detaching surface of the actual collapse block.
The details of the numerical approach to calculate the expression of the detaching surface of the collapsing block for tunnels subjected to the pore water pressure are described in the following.

3.1. Runge–Kutta Algorithm Solving Differential Equation

Taking a kinematically admissible velocity field (one surface) as an example illustrates the numerical approach. Substituting the parameters of an engineering case into equations (13), (13) turns to be

where , , , , . Substituting the parameters of the engineering case into the parameters , t, p, q, m, then, all the parameters are known.

The nonlinear second-order differential equation (equation (15)) is transferred into differential equations.

where and .

For numerical calculation of equation (16), y₁ (x = 0) and y₂ (x = 0) should be assigned with initial values. () represents the ordinate of the apex of the surface y = f (x). It is noted that should not be assigned too small in the Cartesian coordinates (Figure 2(a)) so that the surface y = f (x) can intersect with the tunnel profile. Because of the symmetry of stress, , which deduces that and is taken in the numerical calculation process [20].

Runge–Kutta algorithm is used to solve the differential equations (equation (16)) to obtain a series of points in the Cartesian coordinates. The points are fitted by the power exponential model (equation (17)) referencing to the analytical solution of the detaching surface in Fraldi and Guarracino [20] to gain the expression of the surface [47].

Then, we obtain half of the surface in the first quadrant (x > 0, y > 0) in the Cartesian coordinates, the whole surface can be obtained by symmetry.

L is half of the span of the surface (f (x)), which can be solved by equations (equation (18)). The positive value of L is chosen to calculate the functional Ω.

The functional Ω of the surface f (x) is calculated through equation (7).

3.2. Traversal Search to Find Detaching Surface

Traversal search is used to find the detaching surface of the actual collapsing block. Based on the analysis process above, we can obtain a group of surfaces by assigning different initial values of ( is in constant and ). Through calculating the functional Ω (equation (7)) of each surface, the surface which satisfies the work equation (equation (9)) is the actual detaching surface. Actually, the functional Ω of the actual detaching surface is closest to zero in numerical analysis. It is noted that the range of is the range of traversal search, which is from the tunnel crown to the joint of the external envelope of the group of surfaces with the y-axis along the y-axis.

It should be noted that there are two special situations during the numerical analysis. If the functional of all the surfaces is always positive, it illustrates that the external rate of work is always larger than the internal rate of dissipation, which means the rock masses is unsteady and the external forces cannot be carried by the rock masses. If the functional of all the surfaces is always negative, it means the external rate of work is always less than the internal rate of dissipation, which means the rock mass is stable and no collapse happens. In summary, for the process of numerical calculation, the shape of the detaching surface is obtained by solving the variational equation (equation (11)), and the position of the detaching surface of the collapsing block is obtained by solving the work equation (equation (9)).

4. Numerical Results

4.1. Comparisons with Theoretical Solutions

To validate the numerical approach, taking the limit analysis of the collapse mechanism for tunnel roofs in the absence of the pore water pressure conducted by Fraldi and Guarracino [20] as an example, the example is analyzed by the proposed numerical approach, and the numerical results of expressions of detaching surfaces are compared with Fraldi and Guarracino’s theoretical solutions. When the numerical approach is used to analyze the case in the absence of the pore water pressure, in equation (7) is set to be zero.

The comparison of the height and width of the collapsing block between theoretical solutions [20] and numerical solutions with varied parameters (B, A, σt, and the unit weight of rock masses γ) are shown in Tables 1–4. Numerical solutions are generally consistent with Fraldi and Guarracino’s theoretical solutions [20]. The ratios of all the differences are almost lower than 0.5%, despite one value of 3.93%, as shown in Table 1. Therefore, the numerical approach based on the Runge–Kutta algorithm and traversal search presented in this article is validated to analyze the collapse mechanism of tunnels.

Table 1: Comparison of height H and width L of collapse block between theoretical solutions and numerical solutions with varied parameter B.

Variable
Theoretical solutions
Numerical solutions
Differences

B
Ht (m)
Lt (m)
Hn (m)
Ln (m)
ΔH = Ht − Hn (m)
ΔL = Lt − Ln (m)
ΔH/Ht (%)
ΔL/Lt (%)

0.90
8.4444
8.2800
8.4424
8.2800
0.0020
0.0000
0.0237
0.0000

0.80
9.0000
12.8149
8.9980
12.8150
0.0020
−0.0001
0.0222
−0.0008

0.70
9.7143
19.7567
9.7129
19.7568
0.0014
−0.0001
0.0144
−0.0005

0.60
10.6667
30.3075
10.7078
30.3078
−0.0411
−0.0003
−0.3853
−0.0010

0.50
12.0000
46.1880
12.4712
46.2126
−0.4712
−0.0246
−3.9267
−0.0533

Note: the parameters in the calculation are A = 2/3, σc = 1 × 10⁴ kPa, σt = σc/100 = 100 kPa, γ = 25.0 kN/m³.

Table 2: Comparison of height H and width L of collapse block between theoretical solutions and numerical solutions with varied parameter A.

Variable
Theoretical solutions
Numerical solutions
Differences

A
Ht (m)
Lt (m)
Hn (m)
Ln (m)
ΔH = Ht − Hn (m)
ΔL = Lt − Ln (m)
ΔH/Ht (%)
ΔL/Lt (%)

0.75
9.7143
22.2263
9.7132
22.2264
0.0011
−0.0001
0.0113
−0.0004

0.67
9.7143
19.7567
9.7130
19.7569
0.0013
−0.0002
0.0134
−0.0010

0.50
9.7143
14.8175
9.7120
14.8177
0.0023
−0.0002
0.0237
−0.0013

0.35
9.7143
10.3723
9.7116
10.3724
0.0027
−0.0001
0.0278
−0.0010

0.15
9.7143
4.4453
9.7122
4.4453
0.0021
0.0000
0.0216
0.0000

Note: the parameters in the calculation are B = 0.7, σc = 1 × 10⁴ kPa, σt = σc/100 = 100 kPa, γ = 25.0 kN/m³.

Table 3: Comparison of height H and width L of collapse block between theoretical solutions and numerical solutions with varied parameter σt/σc.

Variable
Theoretical solutions
Numerical solutions
Differences

σt/σc
Ht (m)
Lt (m)
Hn (m)
Ln (m)
ΔH = Ht − Hn (m)
ΔL = Lt − Ln (m)
ΔH/Ht (%)
ΔL/Lt (%)

1/50
19.4286
32.0949
19.4247
32.0949
0.0039
0.0000
0.0201
0.0000

1/75
12.9524
24.1641
12.9502
24.1643
0.0022
−0.0002
0.0170
−0.0008

1/100
9.7143
19.7567
9.7129
19.7568
0.0014
−0.0001
0.0144
−0.0005

1/200
4.8571
12.1617
4.8566
12.1617
0.0005
0.0000
0.0103
0.0000

1/300
3.2381
9.1565
3.2378
9.1565
0.0003
0.0000
0.0093
0.0000

Note: the parameters in the calculation are A = 2/3, B = 0.7, σc = 1 × 10⁴ kPa, γ = 25.0 kN/m³.

Table 4: Comparison of height H and width L of collapse block between theoretical solutions and numerical solutions with varied parameter γ.

Variable
Theoretical solutions
Numerical solutions
Differences

γ (kN/m³)
Ht (m)
Lt (m)
Hn (m)
Ln (m)
ΔH = Ht − Hn (m)
ΔL = Lt − Ln (m)
ΔH/Ht (%)
ΔL/Lt (%)

25.0
9.7143
19.7567
9.7129
19.7568
0.0014
−0.0001
0.0144
−0.0005

22.5
10.7937
21.9519
10.7922
21.9520
0.0015
−0.0001
0.0139
−0.0005

20.0
12.1429
24.6959
12.1412
24.6959
0.0017
0.0000
0.0140
0.0000

17.5
13.8776
28.2239
13.8757
28.2240
0.0019
−0.0001
0.0137
−0.0004

15.0
16.1905
32.9278
16.1883
32.9279
0.0022
−0.0001
0.0136
−0.0003

Note: the parameters in the calculation are A = 2/3, B = 0.7, σc = 1 × 10⁴ kPa, σt = σc/100 = 100 kPa.

4.2. Comparisons with Numerical Solutions

In this section, the program FLAC3D is used to simulate the process of the tunnel excavation under the steady seepage conditions to validate the proposed limit analysis method and numerical approach. Taking the condition of A = 0.65 in section 4.3 as an example, the collapsing block provided by the proposed numerical approach is compared with the numerical simulation results. The size of the numerical model is 100 m × 100 m × 1 m in the transversal, vertical, and longitudinal directions, and the circular tunnel radius r is 3.75 m and the depth C (C is the vertical distance from the center of the circular tunnel to the ground surface) is 50 m. A coordinate axis is defined with the origin at the center of the circular tunnel. Parameters of the Hoek–Brown failure criterion are shown in Table 5, and both the Hoek–Brown parameters (a, mb, s, σc, GSI) and the parameters of corresponding Mohr envelope form in equation (1) (A, B, σc) are presented [20, 39]. The saturated unit weight of rock masses γsat is 28.33 kN/m³. The lateral displacement boundaries are fixed in the normal direction and the displacement boundaries at the bottom are fixed in both the horizontal and vertical directions. The ground surface boundaries are free and permeable. The water table lies at the ground surface. In the process of numerical simulation, the tunnel is excavated after the initial stress reaches the equilibrium state.

Table 5: Parameters of Hoek–Brown failure criterion.

Equivalent Mohr envelope parameters in equation (1)
The generalized Hoek–Brown parameters in FLAC3D

A
B
σc (kPa)
a
mb
s
σc (kPa)
GSI

0.65
5/6
3 × 10⁴
0.5
4.5
0.02
3 × 10⁴
65

Figure 4 shows the contour of vertical displacement of the tunnel roof under the steady seepage conditions of the numerical simulation. It is shown that the shape of the contour lines is approximately consistent with the detaching surface of the collapsing block obtained by the proposed numerical approach. Therefore, the limit analysis method and numerical approach presented in this article are validated to analyze the collapse mechanism of tunnels.

Figure 4: Contour of vertical displacement of tunnel roof under steady seepage conditions of numerical simulation (m).

4.3. Effects of Parameters on Collapse Mechanisms

In this section, the effects of different parameters on the collapse mechanism for circular unlined tunnels under the steady seepage conditions are investigated by the numerical approach. The average quality rock mass is analyzed in the following limit analysis [41]. Detaching surfaces of collapsing blocks in tunnel roof strata varied with parameters are shown in Figures 5–11.

Figure 5: Effect of material parameter A of Hoek–Brown criterion on collapse mechanisms. The parameters are B = 5/6, σc = 3 × 10⁴ kPa, σt = σc/100 = 300 kPa,  − ht = 16.68 m, r = 3.75 m, C = 50 m, and γsat = 28.33 kN/m³.

Figure 6: Effect of material parameter B of Hoek–Brown criterion on collapse mechanisms. The parameters are A = 0.35, σc = 3 × 10⁴ kPa, σt = σc/100 = 300 kPa,  − ht = 16.68 m, r = 3.75 m, C = 50 m, and γsat = 28.33 kN/m³.

Figure 7: Effect of tensile strength σt and compressive strength σc on collapse mechanisms. The parameters are A = 0.35, B = 5/6,  − ht = 16.68 m, r = 3.75 m, C = 50 m, γsat = 28.33 kN/m³, and σt = σc/100.

Figure 8: Effect of saturated unit weight of rock masses γsat on collapse mechanisms. The parameters are A = 0.35, B = 5/6, σc = 3.1 × 10⁴ kPa, σt = σc/100 = 310 kPa,  − ht = 16.68 m, r = 3.75 m, and C = 50 m.

Figure 9: Effect of circular tunnel radius r on collapse mechanisms. The parameters are A = 0.35, B = 5/6, σc = 3 × 10⁴ kPa, σt = σc/100 = 300 kPa,  − ht = 16.68 m, C = 50 m, and γsat = 28.33 kN/m³.

Figure 10: Effect of hydraulic head difference ( − ht) on collapse mechanisms. The parameters are A = 0.35, B = 5/6, σc = 3.1 × 10⁴ kPa, σt = σc/100 = 310 kPa, r = 3.75 m, C = 50 m, γsat = 28.33 kN/m³, and  = 0.

(a)

(b)

(c)

(d)

Figure 11: Effect of depth of tunnel axis C on collapse mechanisms. The parameters are A = 0.35, B = 5/6, σc = 3 × 10⁴ kPa, σt = σc/100 = 300 kPa,  = 0, h (y = 50 m) = −16.68 m, r = 3.75 m, and γsat = 28.33 kN/m³.

Parameters A and B are material parameters of the Hoek–Brown criterion to describe the character of rock masses. The effects of parameters A and B on the collapse mechanism are contradictory (Figures 5 and 6). As parameter A increases, a wider collapsing block will form, while the height of the collapsing block is almost unchanged. Contrary to A, as parameter B increases, a narrower collapsing block is obtained and the height of the collapsing block decreases slightly. It is noted that parameter B is the power exponent of the Hoek–Brown criterion, which controls the curvature of the criterion in the σ-τ plane; when B = 1, the yield locus of Hoek–Brown criterion is straight and is coincident with Mohr–Coulomb criterion [20]. When B is closed to 1, the shape of the detaching surface is almost straight (Figure 6).

The effect of parameters, such as the compressive and tensile strength of rock masses σt and σc, the saturated unit weight of rock masses γsat, the hydraulic head difference between the ground surface and the tunnel circumference  − ht, and the depth of tunnel axis C, on the collapse mechanism can be explained by the upper bound theorem. The tensile strength of rock masses σt is one of the most important factors for the collapse mechanism of tunnels. As parameters σt and σc increase, a larger collapsing block is formed (Figure 7). According to the upper bound theorem, as σt and σc increase, the internal rate of dissipation D increases; thus, the external rates of work consisting of and increase accordingly. When the pore water pressure is constant, the collapsing block will enlarge to increase to balance D.

As parameter γsat increases, the dimension of the collapsing block decreases (Figure 8). According to the upper bound theorem, as γsat increases, the dimension of the collapsing block will decrease to maintain , when the internal rate of dissipation and the pore water pressure is stable.

Figure 9 shows the effect of the circular tunnel radius r on collapse mechanisms. As parameter r increases, the dimension of the collapsing block increases sharply, as shown in Table 6.

Table 6: Dimensions of collapsing block varied with tunnel radius r.

r
H (m)
L (m)
P (kN/m)

3.75
0.2960
0.4236
3.81

5.00
1.2822
1.6050
60.63

6.25
2.4070
2.9206
201.98

7.50
3.6527
4.4024
449.33

Note: H, L, and P represent the height, width, and weight per unit length of collapsing block, respectively. The parameters are A = 0.35, B = 5/6, σc = 3 × 10⁴ kPa, σt = σc/100 = 300 kPa,  − ht = 16.68 m, C = 50 m, andγsat = 28.33 kN/m³.

As parameter  − ht decreases, the dimension of the collapsing block reduces quickly, as shown in Figure 10. When the hydraulic head of ground surface is constant and the head difference  − ht decreases, the hydraulic head around tunnel ht increases (as ht < 0, |ht| decreases). Thus, the pore water pressure pu around tunnel increases and increases. According to the upper bound theorem, as the internal rate of dissipation is constant, the dimension of the collapsing block will reduce to decrease .

Figure 11 shows the effect of the depth of tunnel axis C on collapse mechanisms. In the two-dimensional steady seepage field, assuming the hydraulic head is constant, that is,  = 0 and h (y = 50 m) = −16.68 m, as parameter C decreases, the dimension of the collapsing block increases. When parameter C decrease, the head difference  − ht = |ht| ( = 0) decreases and |z| decreases faster, thus the pressure head decreases. Thus, the pore water pressure pu around the tunnel decreases and decreases. According to the upper bound theorem, as the internal rate of dissipation is constant, the dimension of the collapsing block will enlarge to increase .

5. Discussion about Application of Numerical Approach

As mentioned above, the limit analysis method and numerical approach for the collapse mechanism of deep tunnel roofs under the steady seepage conditions in a homogeneous stratum have been described. In fact, the proposed numerical approach based on the Runge–Kutta algorithm and traversal search method can solve complex equations effectively, which provides an effective technique for solving the limit analysis problem for collapse mechanisms of tunnel roofs under complex conditions. The application of the numerical approach in a shallow tunnel and layered strata will be discussed in the following.

5.1. Application in Shallow Tunnels

The schematic diagram of the collapse mechanism for shallow tunnels is shown in Figure 12. The limit analysis of the collapse mechanism for the shallow tunnel roof subjected to the pore water pressure can also be investigated by the proposed numerical approach. The solving procedures of the numerical approach for shallow tunnels are slightly different from that of deep tunnels, and the main steps of the process are as follows.

Figure 12: Schematic diagram of collapse mechanism for the shallow tunnel.

The work equation Ω = 0 and the variational equation are the conditions to calculate the detaching surface of the collapsing block. For a shallow tunnel subjected to the pore water pressure, height H of the collapsing block is known, while l is uncertain. Therefore, firstly the expression of the functional Ω should be set up according to the condition of a shallow tunnel subjected to the pore water pressure [49, 50]. Secondly, through solving variational equation by assigning , as initial values, a surface that passes through the origin of the coordinates in Figure 12 is obtained, which is expressed as f (x) = a (1) x ^ a (2). By assigning different initial values l, , , a group of surfaces is obtained, which are expressed as f (x) = a (1) (x − l) ^ a (2). Through calculating the functional Ω of each surface, the surface that satisfies Ω = 0 is the actual detaching surface.

5.2. Application in Layered Strata

The schematic diagram of the collapse mechanism of tunnels in layered strata is shown in Figure 13. To obtain the collapse mechanism of tunnels in layered strata through the proposed numerical approach, two assumptions are made: (1) the collapse process of the rock masses at the tunnel roof is progressive, and the collapse starts from the lower stratum to the upper stratum; (2) each stratum is connected with the other, but there is no cohesion between them.

Figure 13: Schematic diagram of collapse mechanism of the tunnel in layered strata.

Based on the assumptions, the collapse of the tunnel roof can be analyzed from the lower stratum to the upper stratum sequentially. If the collapsing block forms and l₁ > 0 in the lower stratum, the collapse mechanism of upper strata should be analyzed (it is called case 1), as shown in Figure 13; if the collapsing block is limited in the lower stratum, it is considered that the lower stratum is stable and upper strata are supported by lower strata (it is called case 2). Case 1 is investigated by the proposed numerical approach for the shallow tunnel; case 2 is investigated by the proposed numerical approach for the deep tunnel.

In addition, it is different from the existing analytical method [51], where the boundary of each stratum in the collapsing block is connected. When the collapse of the tunnel roof is analyzed by the proposed numerical approach from the lower stratum to the upper stratum sequentially, the boundary of collapsing block 1 and collapsing block 2 could be discrete, that is, xp ≤ xq.

Based on the discussion above, it can be summarized that the numerical approach can be used to analyze the cases under different cover depths, roof stratification, and pore water pressure conditions. Moreover, the proposed numerical approach extends the application of limit analysis of collapse mechanisms of tunnels.

6. Conclusions

The collapse mechanism of a circular unlined tunnel roof subjected to the pore water pressure is investigated. The analysis is performed in the framework of the upper bound theorem of limit analysis and the extremum principle. The work equation and the variational equation about the collapsing block are derived from the upper bound theorem and the extremum principle, and the analytical solution of the pore water pressure distribution around the tunnel in the two-dimensional steady seepage field is employed in the equations. The proposed numerical approach solves the complex equations effectively to obtain the function expression of the detaching surface of the collapsing block, where the Runge–Kutta algorithm is applied to solve the variational equation to obtain the possible surfaces, and the traversal search method is used to find the effective detaching surface of the collapsing block satisfied the work equation. The proposed limit analysis method and numerical approach are verified by comparison with existing theoretical solutions and the numerical simulation result. The proposed limit analysis method and numerical approach are suitable for deep, shallow tunnels and layered strata.

A simplified two-dimensional (2D) model is formed for analyzing collapse mechanisms of tunnel roofs based on the assumption of plane strain in this article. This is a strict approximate assumption, which will result in deviations. Therefore, the qualitative conclusions are provided in this article. Accurate results need to be drawn from the 3D model, while the complexity of 3D calculation increases remarkably, which will be the object of further studies. The effects of different parameters on the collapse mechanism are investigated, and the qualitative results are as follows:

(1)The soil mass parameters A and B of the Hoek–Brown criterion have a significant influence on the shape of the collapsing block. When parameter A increases or parameter B decreases, the width of the collapsing block increases correspondingly, while the height of the collapsing block is almost constant.(2)The effects of other parameters on the shape of the collapsing block are studied. It is found that the width and height of the collapsing block increase, when the compressive and tensile strength of soil mass, tunnel radius, and the hydraulic head difference between the ground surface and the tunnel circumference increase; the width and height of the collapsing block decrease, when the saturated unit weight of soil masses and depth of the tunnel axis increase.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

Acknowledgments

The authors gratefully acknowledge the financial support of the Major Projects of Scientific and Technological Research and Development Plan of China State Railway Group Co., Ltd. (Grant no. K2019G042).

References

D. C. Drucker, W. Prager, and H. J. Greenberg, “Extended limit design theorems for continuous media,” Quarterly of Applied Mathematics, vol. 9, no. 4, pp. 381–389, 1952.

W. F. Chen, Limit Analysis and Soil Plasticity, Elsevier Scientific Publishing Company, Amsterdam, Netherlands, 1975.

S. W. Sloan, “Geotechnical stability analysis,” Géotechnique, vol. 63, no. 7, pp. 531–571, 2013.

Z.-h. Qian, J.-f. Zou, J. Tian, and Q.-j. Pan, “Estimations of active and passive earth thrusts of non-homogeneous frictional soils using a discretisation technique,” Computers and Geotechnics, vol. 119, no. 3, Article ID 103366, 2020.

A. S. Osman, “Upper bound solutions for the shape factors of smooth rectangular footings on frictional materials,” Computers and Geotechnics, vol. 115, Article ID 103177, 2019.

H. W. Song and H. T. Wang, “Analytical approach for face stability estimate of tunnel with pipe roof reinforcment [J],” Advanced Materials Research, vol. 243-249, pp. 347–350, 2011.

J. H. Atkinson and D. M. Potts, “Stability of a shallow circular tunnel in cohesionless soil,” Geotechnique, vol. 27, no. 2, pp. 203–215, 1977.

E. H. Davis, M. J. Gunn, R. J. Mair et al., “The stability of shallow tunnels and underground openings in cohesive material,” Geotechnique, vol. 30, no. 4, pp. 397–416, 1980.

E. Leca and L. Dormieux, “Upper and lower bound solutions for the face stability of shallow circular tunnels in frictional material,” Geotechnique, vol. 40, no. 4, pp. 581–606, 1990.

G. Mollon, D. Dias, and A. H. Soubra, “Face stability analysis of circular tunnels driven by a pressurized shield,” Journal of Geotechnical and Geoenvironmental Engineering, vol. 136, no. 1, pp. 215–229, 2010.

C. P. Zhang, K. H. Han, and D. L. Zhang, “Face stability analysis of shallow circular tunnels in cohesive–frictional soils,” Tunnelling and Underground Space Technology, vol. 50, pp. 345–357, 2015.

K. H. Han, C. P. Zhang, and D. L. Zhang, “Upper-bound solutions for the face stability of a shield tunnel in multilayered cohesive-frictional soils,” Computers and Geotechnics, vol. 79, pp. 1–9, 2016.

G. Mollon, D. Dias, and A. H. Soubra, “Rotational failure mechanisms for the face stability analysis of tunnels driven by a pressurized shield,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 35, no. 12, pp. 1363–1388, 2011.

J. F. Zou, Z. H. Qian, X. H. Xiang et al., “Face stability of a tunnel excavated in saturated nonhomogeneous soils,” Tunnelling and Underground Space Technology, vol. 83, pp. 1–17, 2019.

W. Li and C. P. Zhang, “Face stability analysis for a shield tunnel in anisotropic sands,” International Journal of Geomechanics, vol. 20, no. 5, Article ID 04020043, 2020.

G. Mollon, D. Dias, and A. H. Soubra, “Continuous velocity fields for collapse and blowout of a pressurized tunnel face in purely cohesive soil,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 37, no. 13, pp. 2061–2083, 2013.

F. Zhang, Y. F. Gao, Y. X. Wu et al., “Upper-bound solutions for face stability of circular tunnels in undrained clays,” Géotechnique, vol. 68, no. 1, pp. 76–85, 2018.

C. P. Zhang, W. Li, W. J. Zhu et al., “Face stability analysis of a shallow horseshoe-shaped shield tunnel in clay with a linearly increasing shear strength with depth,” Tunnelling and Underground Space Technology, vol. 97, Article ID 103291, 2020.

W. Li, C. P. Zhang, W. J. Zhu et al., “Upper-bound solutions for the face stability of a non-circular NATM tunnel in clays with a linearly increasing undrained shear strength with depth,” Computers and Geotechnics, vol. 114, Article ID 103136, 2019.

M. Fraldi and F. Guarracino, “Limit analysis of collapse mechanisms in cavities and tunnels according to the Hoek–Brown failure criterion,” International Journal of Rock Mechanics and Mining Sciences, vol. 46, no. 4, pp. 665–673, 2009.

F. Huang, R. C. Ou, Z. L. Li et al., “Limit analysis for the face stability of a shallow-shield tunnel based on a variational approach to the blow-out failure mode,” International Journal of Geomechanics, vol. 18, no. 6, Article ID 4018038, 2018.

F. Huang and X. L. Yang, “Upper bound limit analysis of collapse shape for circular tunnel subjected to pore pressure based on the Hoek–Brown failure criterion,” Tunnelling and Underground Space Technology, vol. 26, no. 5, pp. 614–618, 2011.

C. B. Qin, S. C. Chian, and X. L. Yang, “2D and 3D limit analysis of progressive collapse mechanism for deep-buried tunnels under the condition of varying water table,” International Journal for Rock Mechanics and Mining Sciences, vol. 80, pp. 255–264, 2015.

C. B. Qin and S. C. Chian, “2D and 3D stability analysis of tunnel roof collapse in stratified rock: a kinematic approach,” International Journal of Rock Mechanics and Mining Sciences, vol. 100, pp. 269–277, 2017.

X. L. Yang and C. Yao, “Stability of tunnel roof in nonhomogeneous soils,” International Journal of Geomechanics, vol. 18, no. 3, Article ID 6018002, 2018.

R. Zhang and X. L. Yang, “Limit analysis of active and passive mechanisms of shallow tunnels in nonassociative soil with changing water table,” International Journal of Geomechanics, vol. 18, no. 7, Article ID 4018063, 2018.

H. T. Wang, L. G. Wang, and S. C. Li, “Roof collapse mechanisms for a shallow tunnel in two-layer rock strata incorporating the influence of groundwater,” Engineering Failure Analysis, vol. 98, pp. 215–227, 2019.

I. M. Lee and S. M. Nam, “Effect of tunnel advance rate on seepage forces acting on the underwater tunnel face,” Tunnelling and Underground Space Technology, vol. 19, no. 3, pp. 273–281, 2004.

Q. J. Pan and D. Dias, “The effect of pore water pressure on tunnel face stability,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 40, no. 15, pp. 2123–2136, 2016.

P. Perazzelli, T. Leone, and G. Anagnostou, “Tunnel face stability under seepage flow conditions,” Tunnelling and Underground Space Technology, vol. 43, pp. 459–469, 2014.

X. L. Lv, Y. C. Zhou, M. S. Huang et al., “Computation of the minimum limit support pressure for the shield tunnel face stability under seepage condition,” International Journal of Civil Engineering, vol. 15, pp. 849–863, 2017.

X. L. Lv and F. D. Li, ““Study on the stability of large cross-river shield tunnel face with seepage,” Progress in industrial and civil engineering II, PTS 1-4, vol. 405-408, pp. 1371–1374, 2013.

M. Huangfu, M. S. Wang, Z. S. Tan et al., “Analytical solutions for steady seepage into an underwater circular tunnel,” Tunnelling and Underground Space Technology, vol. 25, no. 4, pp. 391–396, 2010.

P. F. Li, F. Wang, Y. Y. Long et al., “Investigation of steady water inflow into a subsea grouted tunnel,” Tunnelling and Underground Space Technology, vol. 80, pp. 92–102, 2018.

J. F. Zou, A. Wei, and L. Liang, “Analytical solution for steady seepage and groundwater inflow into an underwater tunnel,” Geomechanics and Engineering, vol. 20, no. 3, pp. 267–273, 2020.

R. E. Goodman, D. G. Moye, A. Van Schalkwyk et al., “Ground water inflows during tunnel driving,” Bulletin of the Association of Engineering Geologists, vol. 2, pp. 35–56, 1965.

P. Y. Polubarinova-kochina, Theory of Ground Water Movement, Princeton University Press, Princeton, NJ, USA, 1962.

D. Kolymbas and P. Wagner, “Groundwater ingress to tunnels-The exact analytical solution,” Tunneling and Underground Space Technology, vol. 22, no. 1, pp. 23–27, 2007.

E. Hoek and E. T. Brown, “Empirical strength criterion for rock masses,” Journal of the Geotechnical Engineering Division-ASCE, vol. 106, no. 9, pp. 1013–1035, 1980.

E. Hoek and M. S. Diederichs, “Empirical estimation of rock mass modulus,” International Journal of Rock Mechanics and Mining Sciences, vol. 43, no. 2, pp. 203–215, 2006.

E. Hoek and E. T. Brown, “Practical estimate the rock mass strength,” International Journal of Rock Mechanics and Mining Sciences, vol. 34, no. 8, pp. 1165–1186, 1997.

D. W. Taylor, Fundamentals of Soil Mechanics, John Wiley and Sons, Inc., New York, NY, USA, 1948.

Z. Y. Chen, J. X. Zhou, and H. J. Wang, Soil Mechanics, Tsinghua University Press, Beijing, Chinese, 1994, in Chinese.

Z. Y. Chen, “Discussion on “Finite element calculation on circular slip of earth slope under seepage action,” Chinese Journal of Geotechnical Engineering, vol. 24, no. 3, pp. 394–396, 2002, in Chinese.

R. Hill, The Mathematical Theory of Plasticity, Oxford University Press, New York, NY, USA, 1950.

C. Lanczos, The Variational Principles of Mechanics, University of Toronto Press, Toronto, Canada, 1949.

C. B. Moler, Numerical Computing with Matlab Revised in 2013, Beihang University Press, Beijing, China, 2014.

S. C. Chapra, Applied Numerical Methods with Matlab for Engineers and Scientists, McGraw-Hill Education, New York, NY, USA, 2017.

C. P. Zhang and K. H. Han, “Collapsed shape of shallow unlined tunnels based on functional catastrophe theory,” Mathematical Problems in Engineering, vol. 2015, Article ID 681257, 2015.

M. Fraldi, R. Cavuoto, A. Cutolo et al., “Stability of tunnels according to depth and variability of rock mass parameters,” International Journal of Rock Mechanics and Mining Sciences, vol. 119, pp. 222–229, 2019.

H. T. Wang, P. Liu, C. Liu et al., “Three-dimensional upper bound limit analysis on the collapse of shallow soil tunnels considering roof stratification and pore water pressure,” Mathematical Problems in Engineering, vol. 2019, Article ID 8164702, 2019.