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

Authors: Jing-jing Liu; Tie-lin Chen; Chang-ling Xie; Jian-hua Tian; Yu-xin Wei
Journal: Advances in Civil Engineering
 (2021)

---

## 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.


---

## Body

## 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 expressionf (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 flowIn this article, the case of the hydraulic head of ground surfacehw = 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 hw = 0. Point B is in the tunnel, its pressure head is pu/γw (pu is the pore water pressure, γw is unit weight of water), elevation head is –z (z > 0), and its total head ht = pu/γw − z < 0. The head loss of flow line AB is hw − ht = z − pu/γw.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 Figure2(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 w, and the stresses on it are shown in detail. The relative displacement and velocity of the entire collapsing block are vertical displacement u and velocity u˙ in the admissible velocity field.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.
(a)(b)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]:(1)τn=Aσcσn′+σtσcB,A,B∈0,1,where n represents the unit vector normal to the failure plane; σn′ is the normal effective stress; τn is the shear stress; A and B are the material constants; σc is the uniaxial compressive strength of the intact rock; σt is the tensile strength of the rock mass σc>σt≥0.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]:(2)σn′=−σt+σcABf′x1/1−B,τn=Aσcσn′+σtσcB,(3)ε˙n=u˙w11+f′x2,γ˙n=−u˙wf′x1+f′x2,where n represents the unit vector normal to the detaching surface (f (x)); ε˙n and γ˙n are normal and shear plastic strain rate, respectively; the prime represents the derivation of the function to its variable, i.e., f′x=∂fx/∂x.Because the collapsing block is symmetrical with respect to they-axis, half of the collapsing block is analyzed in the following. The internal rate of dissipation D along the detaching surface f (x) is(4)D=∫0Sσn′ε˙n+τnγ˙ndS=∫0S−σt+1−B−1σcABf′x1/1−B⋅u˙w11+f′x2dS,where S is the length of half detaching surface (f (x)); dS=1+f′x2dx 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 W˙ consists of the work rate of the total weight of collapsing block W˙1 and the work rate of the pore water pressure acting on the detaching surface W˙2. 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, Cx=C−r2−x2; pu (x, f (x)) is the pore water pressure acting on the detaching surface (f (x)); cosn,y=1/1+f′x2.Figure 3
Diagram of external force on collapsing block.(5)W˙1=∫0LγsatCx−fxu˙dx,W˙2=∫0Spux,fxcosn,yu˙dS,Based on the above assumptions, the pore water pressurepu (x, y) in the two-dimensional seepage field for a tunnel is as follows [38]:(6)pu=pux,y=γwhw+hw−htlnR/r⋅lnrRx2+C−y−αr2r+αy−αC2+α2x2+y,where γw is the unit weight of water; hw is the total head of ground surface; ht is the total head of the circumference of the tunnel; hw − 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; α=C−C2−r2/r, R=r2/C−C2−r2.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:(7)Ω=D+W˙1+W˙2=∫0L−σt+1−B−1σcABf′x1/1−Bu˙wdx+∫0LγsatCx−fxu˙dx+∫0Lpux,fxu˙dx=∫0L−σt+1−B−1σcABf′x1/1−B+γsatCx−fx+pux,fxu˙dx=∫0LΛfx,f′x,xdx.Λfx,f′x,x is a functional, which is expressed as follows:(8)Λ=Λfx,f′x,x=−σt+1−B−1σcABf′x1/1−B+γsatCx−fx+pux,fx.The thickness of the plastic detaching zone (Figure2(a)) w 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(9)Ωfx,f′x,x=0,that is,(10)∫0L−σt+1−B−1σcABf′x1/1−B+γsatCx−fx+pux,fxu˙dx=0.


### 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.(11)δΩfx,f′x,x=0.Euler–Lagrange equation (equation (12)) transforms the variational equation (equation (11)) into the differential equation (equation (13)) [46]:(12)∂Λ∂fx−∂∂x∂Λ∂f′x=0,that is,(13)u˙1−B−1σcAB1/1−Bf′x2B−1/1−Bf″x+γwhw−htlnR/rfx−C2−r2x2+fx−C2−r22−fx+C2−r2x2+fx+C2−r22+γw−γsat=0.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 blockf (x) is obtained through the numerical approach presented in the following, the weight of collapsing block per unit length, namely, the collapse load P, 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.(14)P=2∫0LγsatCx−fxdx.


## 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 flowIn this article, the case of the hydraulic head of ground surfacehw = 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 hw = 0. Point B is in the tunnel, its pressure head is pu/γw (pu is the pore water pressure, γw is unit weight of water), elevation head is –z (z > 0), and its total head ht = pu/γw − z < 0. The head loss of flow line AB is hw − ht = z − pu/γw.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 Figure2(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 w, and the stresses on it are shown in detail. The relative displacement and velocity of the entire collapsing block are vertical displacement u and velocity u˙ in the admissible velocity field.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.
(a)(b)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]:(1)τn=Aσcσn′+σtσcB,A,B∈0,1,where n represents the unit vector normal to the failure plane; σn′ is the normal effective stress; τn is the shear stress; A and B are the material constants; σc is the uniaxial compressive strength of the intact rock; σt is the tensile strength of the rock mass σc>σt≥0.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]:(2)σn′=−σt+σcABf′x1/1−B,τn=Aσcσn′+σtσcB,(3)ε˙n=u˙w11+f′x2,γ˙n=−u˙wf′x1+f′x2,where n represents the unit vector normal to the detaching surface (f (x)); ε˙n and γ˙n are normal and shear plastic strain rate, respectively; the prime represents the derivation of the function to its variable, i.e., f′x=∂fx/∂x.Because the collapsing block is symmetrical with respect to they-axis, half of the collapsing block is analyzed in the following. The internal rate of dissipation D along the detaching surface f (x) is(4)D=∫0Sσn′ε˙n+τnγ˙ndS=∫0S−σt+1−B−1σcABf′x1/1−B⋅u˙w11+f′x2dS,where S is the length of half detaching surface (f (x)); dS=1+f′x2dx 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 W˙ consists of the work rate of the total weight of collapsing block W˙1 and the work rate of the pore water pressure acting on the detaching surface W˙2. 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, Cx=C−r2−x2; pu (x, f (x)) is the pore water pressure acting on the detaching surface (f (x)); cosn,y=1/1+f′x2.Figure 3
Diagram of external force on collapsing block.(5)W˙1=∫0LγsatCx−fxu˙dx,W˙2=∫0Spux,fxcosn,yu˙dS,Based on the above assumptions, the pore water pressurepu (x, y) in the two-dimensional seepage field for a tunnel is as follows [38]:(6)pu=pux,y=γwhw+hw−htlnR/r⋅lnrRx2+C−y−αr2r+αy−αC2+α2x2+y,where γw is the unit weight of water; hw is the total head of ground surface; ht is the total head of the circumference of the tunnel; hw − 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; α=C−C2−r2/r, R=r2/C−C2−r2.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:(7)Ω=D+W˙1+W˙2=∫0L−σt+1−B−1σcABf′x1/1−Bu˙wdx+∫0LγsatCx−fxu˙dx+∫0Lpux,fxu˙dx=∫0L−σt+1−B−1σcABf′x1/1−B+γsatCx−fx+pux,fxu˙dx=∫0LΛfx,f′x,xdx.Λfx,f′x,x is a functional, which is expressed as follows:(8)Λ=Λfx,f′x,x=−σt+1−B−1σcABf′x1/1−B+γsatCx−fx+pux,fx.The thickness of the plastic detaching zone (Figure2(a)) w 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(9)Ωfx,f′x,x=0,that is,(10)∫0L−σt+1−B−1σcABf′x1/1−B+γsatCx−fx+pux,fxu˙dx=0.


## 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.(11)δΩfx,f′x,x=0.Euler–Lagrange equation (equation (12)) transforms the variational equation (equation (11)) into the differential equation (equation (13)) [46]:(12)∂Λ∂fx−∂∂x∂Λ∂f′x=0,that is,(13)u˙1−B−1σcAB1/1−Bf′x2B−1/1−Bf″x+γwhw−htlnR/rfx−C2−r2x2+fx−C2−r22−fx+C2−r2x2+fx+C2−r22+γw−γsat=0.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 blockf (x) is obtained through the numerical approach presented in the following, the weight of collapsing block per unit length, namely, the collapse load P, 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.(14)P=2∫0LγsatCx−fxdx.


## 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(15)φf′xtf″x+pfx−mx2+fx−m2−fx+mx2+fx+m2+q=0,where φ=1−B−1σcAB1/1−B, t=2B−1/1−B, p=γwhw−ht/lnR/r, q=γw−γsat, m=C2−r2. 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.(16)dy1dx=y2,dy2dx=y2−tφ−py1−mx2+y1−m2−y1+mx2+y1+m2−q,where y1=fx and y2=f′x.For numerical calculation of equation (16), y1 (x = 0) and y2 (x = 0) should be assigned with initial values. y1|x=0 (0<y1|x=0<C−r) represents the ordinate of the apex of the surface y = f (x). It is noted that y1|x=0 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, τxyx=0=0, which deduces that y2|x=0=0 and y2|x=0=0.01 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].(17)fx=a1x^a2+a3.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 Ω.(18)fx=a1x^a2+a3,x2+y−C2−r2=0.The functional Ω of the surfacef (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 ofy1|x=0 (y2|x=0 is in constant and y2|x=0=0.01). 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 y1|x=0 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)).


## 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(15)φf′xtf″x+pfx−mx2+fx−m2−fx+mx2+fx+m2+q=0,where φ=1−B−1σcAB1/1−B, t=2B−1/1−B, p=γwhw−ht/lnR/r, q=γw−γsat, m=C2−r2. 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.(16)dy1dx=y2,dy2dx=y2−tφ−py1−mx2+y1−m2−y1+mx2+y1+m2−q,where y1=fx and y2=f′x.For numerical calculation of equation (16), y1 (x = 0) and y2 (x = 0) should be assigned with initial values. y1|x=0 (0<y1|x=0<C−r) represents the ordinate of the apex of the surface y = f (x). It is noted that y1|x=0 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, τxyx=0=0, which deduces that y2|x=0=0 and y2|x=0=0.01 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].(17)fx=a1x^a2+a3.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 Ω.(18)fx=a1x^a2+a3,x2+y−C2−r2=0.The functional Ω of the surfacef (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 ofy1|x=0 (y2|x=0 is in constant and y2|x=0=0.01). 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 y1|x=0 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, W˙2 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 heightH and width L of collapse block between theoretical solutions and numerical solutions with varied parameter B.
VariableTheoretical solutionsNumerical solutionsDifferencesBHt (m)Lt (m)Hn (m)Ln (m)ΔH = Ht − Hn (m)ΔL = Lt − Ln (m)ΔH/Ht (%)ΔL/Lt (%)0.908.44448.28008.44248.28000.00200.00000.02370.00000.809.000012.81498.998012.81500.0020−0.00010.0222−0.00080.709.714319.75679.712919.75680.0014−0.00010.0144−0.00050.6010.666730.307510.707830.3078−0.0411−0.0003−0.3853−0.00100.5012.000046.188012.471246.2126−0.4712−0.0246−3.9267−0.0533Note: the parameters in the calculation areA = 2/3, σc = 1 × 104 kPa, σt = σc/100 = 100 kPa, γ = 25.0 kN/m3.Table 2
Comparison of heightH and width L of collapse block between theoretical solutions and numerical solutions with varied parameter A.
VariableTheoretical solutionsNumerical solutionsDifferencesAHt (m)Lt (m)Hn (m)Ln (m)ΔH = Ht − Hn (m)ΔL = Lt − Ln (m)ΔH/Ht (%)ΔL/Lt (%)0.759.714322.22639.713222.22640.0011−0.00010.0113−0.00040.679.714319.75679.713019.75690.0013−0.00020.0134−0.00100.509.714314.81759.712014.81770.0023−0.00020.0237−0.00130.359.714310.37239.711610.37240.0027−0.00010.0278−0.00100.159.71434.44539.71224.44530.00210.00000.02160.0000Note: the parameters in the calculation areB = 0.7, σc = 1 × 104 kPa, σt = σc/100 = 100 kPa, γ = 25.0 kN/m3.Table 3
Comparison of heightH and width L of collapse block between theoretical solutions and numerical solutions with varied parameter σt/σc.
VariableTheoretical solutionsNumerical solutionsDifferencesσt/σcHt (m)Lt (m)Hn (m)Ln (m)ΔH = Ht − Hn (m)ΔL = Lt − Ln (m)ΔH/Ht (%)ΔL/Lt (%)1/5019.428632.094919.424732.09490.00390.00000.02010.00001/7512.952424.164112.950224.16430.0022−0.00020.0170−0.00081/1009.714319.75679.712919.75680.0014−0.00010.0144−0.00051/2004.857112.16174.856612.16170.00050.00000.01030.00001/3003.23819.15653.23789.15650.00030.00000.00930.0000Note: the parameters in the calculation areA = 2/3, B = 0.7, σc = 1 × 104 kPa, γ = 25.0 kN/m3.Table 4
Comparison of heightH and width L of collapse block between theoretical solutions and numerical solutions with varied parameter γ.
VariableTheoretical solutionsNumerical solutionsDifferencesγ (kN/m3)Ht (m)Lt (m)Hn (m)Ln (m)ΔH = Ht − Hn (m)ΔL = Lt − Ln (m)ΔH/Ht (%)ΔL/Lt (%)25.09.714319.75679.712919.75680.0014−0.00010.0144−0.000522.510.793721.951910.792221.95200.0015−0.00010.0139−0.000520.012.142924.695912.141224.69590.00170.00000.01400.000017.513.877628.223913.875728.22400.0019−0.00010.0137−0.000415.016.190532.927816.188332.92790.0022−0.00010.0136−0.0003Note: the parameters in the calculation areA = 2/3, B = 0.7, σc = 1 × 104 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/m3. 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 FLAC3DABσc (kPa)ambsσc (kPa)GSI0.655/63 × 1040.54.50.023 × 10465Figure4 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 parameterA of Hoek–Brown criterion on collapse mechanisms. The parameters are B = 5/6, σc = 3 × 104 kPa, σt = σc/100 = 300 kPa, hw − ht = 16.68 m, r = 3.75 m, C = 50 m, and γsat = 28.33 kN/m3.Figure 6
Effect of material parameterB of Hoek–Brown criterion on collapse mechanisms. The parameters are A = 0.35, σc = 3 × 104 kPa, σt = σc/100 = 300 kPa, hw − ht = 16.68 m, r = 3.75 m, C = 50 m, and γsat = 28.33 kN/m3.Figure 7
Effect of tensile strengthσt and compressive strength σc on collapse mechanisms. The parameters are A = 0.35, B = 5/6, hw − ht = 16.68 m, r = 3.75 m, C = 50 m, γsat = 28.33 kN/m3, 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 × 104 kPa, σt = σc/100 = 310 kPa, hw − 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 areA = 0.35, B = 5/6, σc = 3 × 104 kPa, σt = σc/100 = 300 kPa, hw − ht = 16.68 m, C = 50 m, and γsat = 28.33 kN/m3.Figure 10
Effect of hydraulic head difference (hw − ht) on collapse mechanisms. The parameters are A = 0.35, B = 5/6, σc = 3.1 × 104 kPa, σt = σc/100 = 310 kPa, r = 3.75 m, C = 50 m, γsat = 28.33 kN/m3, and hw = 0.Figure 11
Effect of depth of tunnel axisC on collapse mechanisms. The parameters are A = 0.35, B = 5/6, σc = 3 × 104 kPa, σt = σc/100 = 300 kPa, hw = 0, h (y = 50 m) = −16.68 m, r = 3.75 m, and γsat = 28.33 kN/m3.
(a)(b)(c)(d)ParametersA 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 hw − 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 W˙1 and W˙2 increase accordingly. When the pore water pressure is constant, the collapsing block will enlarge to increase W˙1 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 W˙1, when the internal rate of dissipation and the pore water pressure is stable.Figure9 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.
rH (m)L (m)P (kN/m)3.750.29600.42363.815.001.28221.605060.636.252.40702.9206201.987.503.65274.4024449.33Note: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 × 104 kPa, σt = σc/100 = 300 kPa, hw − ht = 16.68 m, C = 50 m, andγsat = 28.33 kN/m3.As parameterhw − ht decreases, the dimension of the collapsing block reduces quickly, as shown in Figure 10. When the hydraulic head of ground surface hw is constant and the head difference hw − ht decreases, the hydraulic head around tunnel ht increases (as ht < 0, |ht| decreases). Thus, the pore water pressure pu around tunnel increases and W˙2 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 W˙1.Figure11 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, hw = 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 hw − ht = |ht| (hw = 0) decreases and |z| decreases faster, thus the pressure head pu/γw decreases. Thus, the pore water pressure pu around the tunnel decreases and W˙2 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 W˙1.


## 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, W˙2 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 heightH and width L of collapse block between theoretical solutions and numerical solutions with varied parameter B.
VariableTheoretical solutionsNumerical solutionsDifferencesBHt (m)Lt (m)Hn (m)Ln (m)ΔH = Ht − Hn (m)ΔL = Lt − Ln (m)ΔH/Ht (%)ΔL/Lt (%)0.908.44448.28008.44248.28000.00200.00000.02370.00000.809.000012.81498.998012.81500.0020−0.00010.0222−0.00080.709.714319.75679.712919.75680.0014−0.00010.0144−0.00050.6010.666730.307510.707830.3078−0.0411−0.0003−0.3853−0.00100.5012.000046.188012.471246.2126−0.4712−0.0246−3.9267−0.0533Note: the parameters in the calculation areA = 2/3, σc = 1 × 104 kPa, σt = σc/100 = 100 kPa, γ = 25.0 kN/m3.Table 2
Comparison of heightH and width L of collapse block between theoretical solutions and numerical solutions with varied parameter A.
VariableTheoretical solutionsNumerical solutionsDifferencesAHt (m)Lt (m)Hn (m)Ln (m)ΔH = Ht − Hn (m)ΔL = Lt − Ln (m)ΔH/Ht (%)ΔL/Lt (%)0.759.714322.22639.713222.22640.0011−0.00010.0113−0.00040.679.714319.75679.713019.75690.0013−0.00020.0134−0.00100.509.714314.81759.712014.81770.0023−0.00020.0237−0.00130.359.714310.37239.711610.37240.0027−0.00010.0278−0.00100.159.71434.44539.71224.44530.00210.00000.02160.0000Note: the parameters in the calculation areB = 0.7, σc = 1 × 104 kPa, σt = σc/100 = 100 kPa, γ = 25.0 kN/m3.Table 3
Comparison of heightH and width L of collapse block between theoretical solutions and numerical solutions with varied parameter σt/σc.
VariableTheoretical solutionsNumerical solutionsDifferencesσt/σcHt (m)Lt (m)Hn (m)Ln (m)ΔH = Ht − Hn (m)ΔL = Lt − Ln (m)ΔH/Ht (%)ΔL/Lt (%)1/5019.428632.094919.424732.09490.00390.00000.02010.00001/7512.952424.164112.950224.16430.0022−0.00020.0170−0.00081/1009.714319.75679.712919.75680.0014−0.00010.0144−0.00051/2004.857112.16174.856612.16170.00050.00000.01030.00001/3003.23819.15653.23789.15650.00030.00000.00930.0000Note: the parameters in the calculation areA = 2/3, B = 0.7, σc = 1 × 104 kPa, γ = 25.0 kN/m3.Table 4
Comparison of heightH and width L of collapse block between theoretical solutions and numerical solutions with varied parameter γ.
VariableTheoretical solutionsNumerical solutionsDifferencesγ (kN/m3)Ht (m)Lt (m)Hn (m)Ln (m)ΔH = Ht − Hn (m)ΔL = Lt − Ln (m)ΔH/Ht (%)ΔL/Lt (%)25.09.714319.75679.712919.75680.0014−0.00010.0144−0.000522.510.793721.951910.792221.95200.0015−0.00010.0139−0.000520.012.142924.695912.141224.69590.00170.00000.01400.000017.513.877628.223913.875728.22400.0019−0.00010.0137−0.000415.016.190532.927816.188332.92790.0022−0.00010.0136−0.0003Note: the parameters in the calculation areA = 2/3, B = 0.7, σc = 1 × 104 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/m3. 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 FLAC3DABσc (kPa)ambsσc (kPa)GSI0.655/63 × 1040.54.50.023 × 10465Figure4 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 parameterA of Hoek–Brown criterion on collapse mechanisms. The parameters are B = 5/6, σc = 3 × 104 kPa, σt = σc/100 = 300 kPa, hw − ht = 16.68 m, r = 3.75 m, C = 50 m, and γsat = 28.33 kN/m3.Figure 6
Effect of material parameterB of Hoek–Brown criterion on collapse mechanisms. The parameters are A = 0.35, σc = 3 × 104 kPa, σt = σc/100 = 300 kPa, hw − ht = 16.68 m, r = 3.75 m, C = 50 m, and γsat = 28.33 kN/m3.Figure 7
Effect of tensile strengthσt and compressive strength σc on collapse mechanisms. The parameters are A = 0.35, B = 5/6, hw − ht = 16.68 m, r = 3.75 m, C = 50 m, γsat = 28.33 kN/m3, 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 × 104 kPa, σt = σc/100 = 310 kPa, hw − 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 areA = 0.35, B = 5/6, σc = 3 × 104 kPa, σt = σc/100 = 300 kPa, hw − ht = 16.68 m, C = 50 m, and γsat = 28.33 kN/m3.Figure 10
Effect of hydraulic head difference (hw − ht) on collapse mechanisms. The parameters are A = 0.35, B = 5/6, σc = 3.1 × 104 kPa, σt = σc/100 = 310 kPa, r = 3.75 m, C = 50 m, γsat = 28.33 kN/m3, and hw = 0.Figure 11
Effect of depth of tunnel axisC on collapse mechanisms. The parameters are A = 0.35, B = 5/6, σc = 3 × 104 kPa, σt = σc/100 = 300 kPa, hw = 0, h (y = 50 m) = −16.68 m, r = 3.75 m, and γsat = 28.33 kN/m3.
(a)(b)(c)(d)ParametersA 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 hw − 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 W˙1 and W˙2 increase accordingly. When the pore water pressure is constant, the collapsing block will enlarge to increase W˙1 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 W˙1, when the internal rate of dissipation and the pore water pressure is stable.Figure9 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.
rH (m)L (m)P (kN/m)3.750.29600.42363.815.001.28221.605060.636.252.40702.9206201.987.503.65274.4024449.33Note: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 × 104 kPa, σt = σc/100 = 300 kPa, hw − ht = 16.68 m, C = 50 m, andγsat = 28.33 kN/m3.As parameterhw − ht decreases, the dimension of the collapsing block reduces quickly, as shown in Figure 10. When the hydraulic head of ground surface hw is constant and the head difference hw − ht decreases, the hydraulic head around tunnel ht increases (as ht < 0, |ht| decreases). Thus, the pore water pressure pu around tunnel increases and W˙2 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 W˙1.Figure11 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, hw = 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 hw − ht = |ht| (hw = 0) decreases and |z| decreases faster, thus the pressure head pu/γw decreases. Thus, the pore water pressure pu around the tunnel decreases and W˙2 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 W˙1.


## 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 Figure12. 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δΩ=0 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 y1x=0=fx=0=0, y2x=0=f′x=0=0.01 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, y1x=l=fx=l=0, y2x=l=f′x=l=0.01, 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 Figure13. 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 andl1 > 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.


## 5.1. Application in Shallow Tunnels

The schematic diagram of the collapse mechanism for shallow tunnels is shown in Figure12. 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δΩ=0 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 y1x=0=fx=0=0, y2x=0=f′x=0=0.01 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, y1x=l=fx=l=0, y2x=l=f′x=l=0.01, 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 Figure13. 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 andl1 > 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 parametersA 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.



---

*Source: 3591670-2021-05-27.xml*