# Functional Catastrophe Analysis of Collapse Mechanism for Shallow Tunnels with Considering Settlement

Authors: Rui Zhang; Hai-Bo Xiao; Wen-Tao Li
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Limit analysis is a practical and meaningful method to predict the stability of geomechanical properties. This work investigates the pore water effect on new collapse mechanisms and possible collapsing block shapes of shallow tunnels with considering the effects of surface settlement. The analysis is performed within the framework of upper bound theorem. Furthermore, the NL nonlinear failure criterion is used to examine the influence of different factors on the collapsing shape and the minimum supporting pressure in shallow tunnels. Analytical solutions derived by functional catastrophe theory for the two different shape curves which describe the distinct characteristics of falling blocks up and down the water level are obtained by virtual work equations under the variational principle. By considering that the mechanical properties of soil are not affected by the presence of underground water, the strength parameters in NL failure criterion can be taken to be the same under and above the water table. According to the numerical results in this work, the influences on the size of collapsing block different parameters have are presented in the tables and the upper bounds on the loads required to resist collapse are derived and illustrated in the form of supporting forces graphs that account for the variation of the embedded depth and other factors.


---

## Body

## 1. Introduction

With the development of cities, it is necessary to use a large amount of the underground area for the construction of transportation infrastructures. The ground movements are caused inevitably during the construction of the shallow tunneling in soft ground. Due to the adverse impact of the presence of underground water, as the occurrence of collapse of the tunnel exerts great threats to people’s lives and engineering loss [1], it is a great significance to predict the stability of the shallow tunnels in the engineering.To estimate the face tunnel stability problem is a hot topic in the tunnel engineering. The shallow tunnels are often chosen in the urban subway projects. Among a large number of methods which are suitable for solving the stability problem of shallow tunnels, limit equilibrium method, numerical simulation, and the limit analysis method are widely used. Owing to the ignorance of the relationship between stress and strain, the equilibrium method only considering stress balance has great defect. Due to the limitations of the limit equilibrium method and other methods, some scholars adopted the limit analysis approach to predict the stability and failure modes of the face and crown of tunnels, which shows an extreme simplicity and great effectiveness. For instance, the rigorous bounds of supporting pressure were obtained by Sloan and Assadi [2] with the help of limit analysis theory and finite element technique. In the 1970s the upper bound theorem was proposed by Chen [3]. Then this theorem had great importance in the field of geotechnical engineering because of its great validity in dealing with the stability problems in underground structures. With the extensive use of it, it was improved by many scholars. By introducing a linear multiblock collapse mechanism in 1980, Davis et al. [4] obtained the upper bound solutions of the stability coefficient. In 1994 the accuracy of supporting pressure derived under the three-dimensional collapse mechanism with the centrifugal model test was studied by Chambon and Corté [5]. Yang and Long [6] studied the influences of two angles on the 3D slope stability by limit analysis method and obtained the critical stability number. The upper bound theorem will be illustrated at length in the following.With the development of the limit analysis method, the linear criterion [7–9] used to evaluate the stability problems has been replaced by the nonlinear criterion [10–15] in terms of the nonlinear mechanical characteristics of geotechnical material in tunnel project. During the process of applying nonlinear failure criterion, a generalized tangential methodology was suggested [16–18] to calculate the energy dissipation rate and the external work rate accurately. In particular on the basis of nonlinear yield criterion [19–24] many researchers put forward the analytical solution for characterizing the collapsing shape. The curved failure mechanism was constructed according to the theory that the energy dissipation rate and external work rate were calculated along the velocity discontinuity [25]. The equation about it will be set in the following.To choose appropriate failure criteria is one of the key steps in predicting the stability of the tunnel roofs. Numerous failure criteria have been used for rock failure analysis, but there is no common agreement of which failure criterion to select. The Mohr-Coulomb failure criterion was recommended for stability analysis because of the more realistic results compared with the different forms of Drucker-Prager [26]. The Modified-Lade failure criterion was developed by Ewy [27]. Furthermore the Hoek-Brown failure criterion [28] has been widely used in geotechnical analysis and still has been improved by some scholars. This work uses a new nonlinear strength function proposed by Baker [29] to analyze the stability of tunnel roof. The NL failure criterion will be presented at length in the following.Owing to a big difference in the tunnel roof collapse mechanism between shallow tunnels and deep tunnels, a new curved failure mechanism should be proposed to reflect the identities more appropriately. Yang and Wang [30] put forward a new failure mechanism of shallow tunnels by considering the surface settlement. On the basis of previous work in which a kinematic plastic solution with ground movements is derived by Osman et al. [31], a compatible displacement field to calculate the stability problem of shallow twin tunnels excavated in the soft layer is constructed by Osman [32].The catastrophe theory has been proved to be effective in analyzing the stability problems in geology and geomechanics. Many scholars have applied this theory in prediction of the stability in the engineering. A fold catastrophe model of a tunnel rock burst was established by Pan et al. [33] to predict the occurrences of a rock burst. Ren et al. [34] studied a cusp catastrophe model to analyze the potential damage mode of the surrounding rock in the tunnels.According to the introduction of the previous works which focus on the predictions of the stability of the tunnel buried in shallow soil layer, this paper establishes a failure mechanism of shallow tunnels with regard to surface settlement. Referring to the NL failure criterion, research on failure is conducted due to the presence of varying water table in the limit analysis. Moreover the functional catastrophe theory is employed to investigate the mechanisms of tunnel roof collapse. The analytical solution of the collapsing block shape curve is obtained and the effects of different parameters on the collapsing block shape are also discussed in this work. Furthermore the upper bound supporting pressure is obtained to ensure the safety of the roof of the shallow tunnels during the constructions.


## 2. Overview of Catastrophe Theory

The catastrophe theory was put forward by Thom in 1972 [35]. As a branch of nonlinear theory, catastrophe theory has been developed rapidly and applied widely. From the main starting points which are the basis of the structural stability, singularity theory, and bifurcation theory; the mathematical model is set for the purpose of the discussions of the general laws of the jump change of the state in dynamic system. The catastrophe theory is a good method to illustrate the mechanics of the change from continuous gradient to stable state in nonlinear system. Due to the complexity and diversity of the constitutive relation of the engineering soils and the uncertainty of the boundaries of the elastic and plastic during the process of the tunnel excavation, it is difficult to describe this process with the method of classical mathematical method. The quantitative state of the system can be predicted with a small number of control variables even without considering the constitutive equation and the mechanical properties of the soil mass system, which is the advantage of the catastrophe theory compared with other theories. In order to explain the phenomenon of various mutations, Thom put forward seven different kinds of models; these models are generally based on elementary catastrophe theory (ECT). In ECT, the potential function of the system is one of the seven elementary functions V=fx1,x2 defined by the polynomial functions shown in Table 1.Table 1
Potential functions used in elementary catastrophe theory.
Name Potential function State variable Fold x 1 3 + t 1 x 1 x 1 Cusp x 1 4 + t 1 x 1 2 + t 2 x 1 x 1 Swallowtail x 1 5 + t 1 x 1 3 + t 2 x 1 2 + t 3 x 1 x 1 Butterfly x 1 6 + t 1 x 1 4 + t 2 x 1 3 + t 3 x 1 2 + t 4 x 1 x 1 Elliptic umbilic x 1 3 - x 1 x 2 2 + t 1 x 2 2 + t 2 x 1 + t 3 x 2 x 1, x2 Hyperbolic umbilic x 1 3 + x 1 x 2 2 + t 1 x 2 2 + t 2 x 1 + t 3 x 2 x 1, x2 Parabolic umbilic x 1 4 + x 1 x 2 2 + t 1 x 2 2 + t 2 x 2 2 + t 3 x 1 + t 4 x 2 x 1, x2However, the total potential of the system always has a complex mathematical form, such as the state functionf(x), rather than the elementary functions. Based on the ECT, Du [36] proposed the functional catastrophe theory (FCT), of which the total potential is employed in the form of a functional instead of an elementary function. Some practical problems were successfully solved in the FCT framework. In the FCT, it is required to obtain the degenerate critical point of the total potential by analyzing the first- and second-order variations of the total potential. In a similar way, the critical value fc(x) of the state variable is obtained when catastrophe occurs. As mentioned above, some researchers have investigated the shape curve f(x) of the collapsing blocks in tunnels in the FCT framework. In this study, the detaching zone of a collapsing block is the studied system. Referring to the NL strength law, the FCT is used to explore the mechanisms of collapse in shallow-buried tunnel with regard to surface settlement. The details of the solution procedures for FCT are illustrated in the following at length.


## 3. Limit Analysis on the Basis of Nonlinear Failure Criterion

### 3.1. NL Strength Law

Various strength functions (Mohr envelopes) have been proposed to represent nonlinear strength functions for soils. Many publications utilized a simple power law relation of the form. Based on the work of Baker [29], the present work employs a slight generalization of this nonlinear relation expressed in the form(1)τn=PaAσPa+Tn,where Pa stands for atmospheric pressure, A,n,T are nondimensional strength parameters, and strength function is defined by (1).The nonlinear strength criterion for rock masses [37] associated with a Mohr envelope is convenient for limiting equilibrium applications (e.g., slope stability). An additional advantage of this form is that the parameters A,n,T have clear physical significance. In particular, A is a scale parameter controlling the magnitude of shear strength; T is a shift parameter controlling the location of the envelope on the σ axis; tNL≡PaT is the tensile strength, with T representing a nondimensional tensile strength; n controls the curvature of the envelope.By assuming the plastic potential,ξ, to be coincident with the Mohr envelope and considering without any loss of generality that τn is positive, one has(2)ξ=τn-PaAσPa+Tn.So the plastic strain rate can be written as follows:(3)ε˙n=λ∂ξ∂σ=-λnAσPa+Tn-1,γ˙n=λ∂ξ∂τ=λ,where λ is a scalar parameter. According to Fraldi and Guarracino [19], the plastic strain rate components can be written in the form(4)ε˙n=u˙w1+f′x2-1/2,γ˙n=-u˙wf′x1+f′x2-1/2,where w is the thickness of the plastic detaching zone. A dot denotes differentiation with respect to time and a prime with respect to x; that is,(5)u˙=∂u∂t,f′x=∂fx∂x.In order to enforce compatibility, from (3) and (4) it follows that(6)λ=-u˙wf′x1+f′x2-1/2.And the normal component of stress can be written as(7)σn=-T+nAf′x1/1-nPa.So, by virtue of the Greenberg minimum principle, the effective collapse mechanism can be found by minimizing the total dissipation; the dissipation density of the internal forces on the detaching surface,D˙i, results:(8)D˙i=σnε˙n+τnγ˙n=Pa·-T+1-n-1nAf′x1/1-nw1+f′x2u˙,where ε˙n is normal plastic strain, γ˙n is shear plastic strain rates, w is the thickness of the plastic detaching zone, and u˙ is the velocity of the collapsing block.


### 3.2. Upper Bound Theorem

The theory of the upper bound has been widely used to the predictions of the stability of the tunnels. The upper bound theorem of limit analysis can be depicted as follows: when the velocity boundary conditions and consistency conditions for strain and velocity are satisfied by the maneuvering-allowable velocity field which is built, the actual loads should be no more than the values of the calculated loads which is derived from the equation constituted by equating the external rate of work and the rates of the internal energy dissipation.According to Chen [3], the upper bound theorem with seepage forces effect can be written as follows:(9)∫Vσij·ε˙ijdV≥∫STi·vidS+∫VXi·vidV+∫V-grad⁡u·vidV,where σij is the stress tensor and ε˙ij is the strain rate in the kinematically admissible velocity field, respectively. Ti is the limit load exerted on the boundary surface. S is the length of velocity discontinuity, Xi is the body strength which is caused by weight, V is the volume of the plastic zone, and vi is the velocity along the velocity discontinuity surface.


## 3.1. NL Strength Law

Various strength functions (Mohr envelopes) have been proposed to represent nonlinear strength functions for soils. Many publications utilized a simple power law relation of the form. Based on the work of Baker [29], the present work employs a slight generalization of this nonlinear relation expressed in the form(1)τn=PaAσPa+Tn,where Pa stands for atmospheric pressure, A,n,T are nondimensional strength parameters, and strength function is defined by (1).The nonlinear strength criterion for rock masses [37] associated with a Mohr envelope is convenient for limiting equilibrium applications (e.g., slope stability). An additional advantage of this form is that the parameters A,n,T have clear physical significance. In particular, A is a scale parameter controlling the magnitude of shear strength; T is a shift parameter controlling the location of the envelope on the σ axis; tNL≡PaT is the tensile strength, with T representing a nondimensional tensile strength; n controls the curvature of the envelope.By assuming the plastic potential,ξ, to be coincident with the Mohr envelope and considering without any loss of generality that τn is positive, one has(2)ξ=τn-PaAσPa+Tn.So the plastic strain rate can be written as follows:(3)ε˙n=λ∂ξ∂σ=-λnAσPa+Tn-1,γ˙n=λ∂ξ∂τ=λ,where λ is a scalar parameter. According to Fraldi and Guarracino [19], the plastic strain rate components can be written in the form(4)ε˙n=u˙w1+f′x2-1/2,γ˙n=-u˙wf′x1+f′x2-1/2,where w is the thickness of the plastic detaching zone. A dot denotes differentiation with respect to time and a prime with respect to x; that is,(5)u˙=∂u∂t,f′x=∂fx∂x.In order to enforce compatibility, from (3) and (4) it follows that(6)λ=-u˙wf′x1+f′x2-1/2.And the normal component of stress can be written as(7)σn=-T+nAf′x1/1-nPa.So, by virtue of the Greenberg minimum principle, the effective collapse mechanism can be found by minimizing the total dissipation; the dissipation density of the internal forces on the detaching surface,D˙i, results:(8)D˙i=σnε˙n+τnγ˙n=Pa·-T+1-n-1nAf′x1/1-nw1+f′x2u˙,where ε˙n is normal plastic strain, γ˙n is shear plastic strain rates, w is the thickness of the plastic detaching zone, and u˙ is the velocity of the collapsing block.


## 3.2. Upper Bound Theorem

The theory of the upper bound has been widely used to the predictions of the stability of the tunnels. The upper bound theorem of limit analysis can be depicted as follows: when the velocity boundary conditions and consistency conditions for strain and velocity are satisfied by the maneuvering-allowable velocity field which is built, the actual loads should be no more than the values of the calculated loads which is derived from the equation constituted by equating the external rate of work and the rates of the internal energy dissipation.According to Chen [3], the upper bound theorem with seepage forces effect can be written as follows:(9)∫Vσij·ε˙ijdV≥∫STi·vidS+∫VXi·vidV+∫V-grad⁡u·vidV,where σij is the stress tensor and ε˙ij is the strain rate in the kinematically admissible velocity field, respectively. Ti is the limit load exerted on the boundary surface. S is the length of velocity discontinuity, Xi is the body strength which is caused by weight, V is the volume of the plastic zone, and vi is the velocity along the velocity discontinuity surface.


## 4. Upper Bound Analysis of Collapse

Many scholars have studied the failure modes of the deep tunnel roof with arbitrary cross sections in layered rocks. By considering the fact that a large number of shallow-buried tunnels are also constructed in the layered stratums, this work investigates the failure mechanism of shallow tunnels under the condition of varying water table considering the surface settlement. The upper bound solutions derived from this work have general applicability and can be used more widely. Due to the deformation continuity of the failure shape of the collapsing mass, the boundary conditions along the detaching lines should be satisfied to ensure the geometry continuity. In order to get the upper bound solutions to describe the shape of the collapsing block, the first work is to calculate the internal energy dissipation rate produced by the shear stress and normal stress along the two different detaching lines. Furthermore the objective function consisting of the internal and external work should be constructed. Lastly two failure shape curvesy=f(x) and y=g(x) can be obtained by the variational approach to reflect the distinct characteristics of falling blocks up and down the water level, as shown in Figure 1.Figure 1
Potential falling blocks with consideration of the effects of surface settlement and pore water pressure.


### 4.1. Assumptions for the Analysis

The failure shape of the collapsing block, as illustrated in Figure1, should be obtained with conditions. Given the behavior of the layered rocks which are under the condition of varying water table should be ideally plastic and follow an associated flow rule; that is, the stress points of the ideally plastic soils should not exceed the yield surface. The width of collapsing block at the ground level, 5i, is only determined by the depth from the center of tunnel to ground surface. Moreover the relationship of the i and H can be described as i=k·H, where i stands for the distance between the tunnel centerline and the point of the trough inflexion, H is the depth of the center of tunnel, and k is a coefficient. By considering the shape of the collapsing block to be symmetrical with respect to the y-axis, the profile of the detaching curves should be smooth and continuous. The collapsing block can be seen as a rigid block without considering the arch effect of shallow circle tunnels.


### 4.2. Computation of Internal Energy Dissipation Rate

From the assumptions above, without considering the geometric difference in different soil layers, the parameters of the NL strength law are assumed to be the same. Due to the presence of velocity detaching line existing in the soil layer of tunnel roofs, the impending failure would slide in a limit state along with the velocity discontinuous surfaces. During the process of the impending collapse, the dissipation densities of the internal forces on the detaching surface are(10)WD=∫L1L2D˙i·w1+f′x2dx+∫L22.5iD˙i·w1+g′x2dx=∫L1L2-T+1-n-1nAf′x1/1-nPa·u˙dx+∫L22.5i-T+1-n-1nAg′x1/1-nPa·u˙dx,where L1 and L2 characterize the collapsing width of upper and lower blocks illustrated in Figure 2. y=f′(x) and y=g′(x) are the first derivative of y=f(x) and y=g(x), respectively.Figure 2
Curved failure mechanism of shallow circular tunnels.


### 4.3. Calculation of External Work

The work rate of failure block produced by weight can be calculated by integral process(11)We=ρ′∫0L1cx·u˙dx+ρ′∫L1L2fx·u˙dx+ρ∫L22.5igx·u˙dx-ρ∫02.5iSx·u˙dx+u˙·ρwL2fL2,where ρ′ is the buoyant weight per unit volume of the rocks. ρ′=ρ-ρw, in which ρ is the weight per unit volume of the rock, and ρw is the unit weight of water. The function c(x) is the function describing the circular tunnel profile. S(x) characterizes the shape of surface settlement(12)cx=H+R-R2-x2.According to the research of Osman [32], the surface settlement is defined by a Gaussian distribution curve. The area of surface settlement trough can be obtained by calculation:(13)∫-2.5i2.5iSxdx=2πiSm.The distribution of excess pore pressure which is derived from the study of Saada et al. [38] can be written as(14)u=p-pw=p-ρwh,where p is the pore water pressure at the considered point which can be obtained by a suitable method p=ruρh, ru stands for pore pressure coefficient, and pw=ρwh is the hydrostatic distribution for pore pressure; h is the vertical distance between the roof of the tunnel and the top of the failure block. So -grad⁡u can be defined as(15)-grad⁡u=ρw-ruρ.Therefore the work rate produced by seepage forces along the velocity discontinuity surface is(16)Wu=u˙∫0L1ρw-ruρcx-fL2dx+u˙∫L1L2ρw-ruρfx-fL2dx.Due to the fact that the tunnel is buried in the shallow strata, the supporting structure is unavoidable for the requirement of safety and stability. Therefore, the work rate of supporting pressure in the shallow circular tunnel is(17)Wq=Rqu˙arcsin⁡L1Rcos⁡π,where q is the supporting pressure exerting on the circumference of tunnel lining.Meanwhile, because of the external force always exerted on the underground structure, the work rate of extra force which puts on the ground surface cannot be ignored. The expressions can be written as(18)WS=2.5iσSu˙,where σS stands for the surcharge load put on the ground surface.


### 4.4. Analytic Solution for Characterizing the Collapsing Shape and Supporting Pressure

In order to describe the shape and extension of the failure collapsing block, it is essential to obtain the explicit expressions ofg(x) and f(x) by constructing an objective function consisting of the external rate of work and the rate of the internal energy dissipation,(19)Λ=WD-We-Wq-Wu-WS.So, the effective collapse mechanism can be obtained by minimizing the objective functionΛ according to the kinematic theorem of limit analysis.Then substituting (10), (11), (16), (17), and (18) into (19), the expression of objective function is given:(20)Λ=∫L1L2ψ1fx,f′x,xu˙dx+∫L22.5iψ2gx,g′x,xu˙dx-ρ′∫0L1cxu˙dx-ρw-ruρ∫0L1cx-fL2·u˙dx+ρ∫02.5iSxu˙dx-2.5iσSu˙+Rqu˙arcsin⁡L1R,in which(21)ψ1fx,f′x,x=-T+1-n-1nAf′x1/1-nPa+ρru-1fx,ψ2gx,g′x,x=-T+1-n-1nAg′x1/1-nPa-ρgx.In order to get the extremum of the objective functionΛ, it is necessary to search for the extremum values of objective functions ψ1 and ψ2. In other words, the analytical solutions of collapse dimension could be obtained by seeking for the minimum values of objective functions ψ1 and ψ2 with the method of variational principle in the realm of plasticity theory. As a consequence the expressions of ψ1 and ψ2 should be turned into two Euler’s equations through the variational method. The expressions of variational equations of ψ1 and ψ2 on stationary conditions can be written as(22)∂ψ1∂fx-∂∂x∂ψ1∂f′x=0,∂ψ2∂gx-∂∂x∂ψ2∂g′x=0.By the variational calculation the explicit forms of the two Euler’s equations for (21) can be obtained as(23)ru-1ρ-PanA1/1-nn-1-1f′x2n-1/1-nf′′x=0,-ρ-PanA1/1-nn-1-1g′x2n-1/1-ng′′x=0.Obviously, (23) are two nonlinear second-order homogeneous differential equations. By the integral calculation process the expressions of velocity discontinuity surface of lower and upper soil layers are(24)fx=k1x+θ1-ruρ1/n+h,gx=k2x+θ′ρ1/n+h′in which(25)k1=A-1/n1-ruρPa1-n/n,k2=A-1/nρPa1-n/n,where θ,θ′,h, and h′ stand for the integration constants coefficients determined by mechanical and geometric boundary conditions, respectively.From what is mentioned above, the detaching curves are supposed to be symmetric with respect to they-axis. It can be seen from Figure 2 that the shear stress on the ground surface equals zero at the point where its x-height is 2.5i:(26)τxyx=2.5i,y=0=0.Furthermore, the explicit expressions of the function of velocity discontinuity surface should fulfill other boundary conditions, such as(27)gx=2.5i=0,(28)fx=L2=gx=L2=n0H,(29)cx=L1=fx=L1,where n0 indicates the ratio of the upper height of roof collapse to the depth from the center of tunnel to ground surface H.Substituting (24) into above expressions, they turn into(30)gx=k2x-2.5i1/n,(31)h=-k1L2+θ1-ruρ1/n+n0H.For the purpose of keeping the whole curve look smooth, another boundary condition should be satisfied:(32)f′x=L2=g′x=L2.Building on this result, the expressions of the function of velocity discontinuity surface can be obtained:(33)fx=k1x-Z1/n-k1L2-Z1/n+n0H,in which(34)Z=2.5i-L2ru1-ru.On the basis of the expression of the profile of the circular tunnel, the piece of external work can be calculated by integratingc(x) over the interval 0,L1: (35)∫0L1cxdx=H+RL1-L12R2-L12-R22arcsin⁡L1R.Therefore, the objective functionΛ can be calculated, which is(36)Λ=u˙ru-1ρHn0-k1ru-1ρL2-Z1/nL2-L1-TPa2.5i-L1+nn+1Pa1-n-1Ak11/1-n+k1ru-1ρ·L2-Zn+1/n-L1-Zn+1/n+nn+1k2ρ-Pa1-n-1Ak21/1-n·L2-2.5in+1/n-ruρL2Hn0+ru-1ρH+RL1-L12R2-L12-R22arcsin⁡L1R+Rqarcsin⁡L1R-2.5iσS+2πiSmρ2.Let Λ=0; the equation can be written as(37)q=-1R·arcsin⁡L1/Rru-1ρHn0-k1ru-1ρL2-Z1/nL2-L1-TPa2.5i-L1+nn+1Pa1-n-1Ak11/1-n+k1ru-1ρ·L2-Zn+1/n-L1-Zn+1/n+nn+1k2ρ-Pa1-n-1Ak21/1-n·L2-2.5in+1/n-ruρL2Hn0+ru-1ρH+RL1-L12R2-L12-R22arcsin⁡L1R-2.5iσS+2πiSmρ2.For the purpose of getting the explicit forms of detaching curve profile consisting off(x) and g(x), the values of L1,L2, and q must be obtained by combining and solving (28), (29), and (37). Then it is not difficult to draw the shape of failure surface by (30) and (33). Meanwhile influences of different factors on the upper bound solution q can be analyzed.


## 4.1. Assumptions for the Analysis

The failure shape of the collapsing block, as illustrated in Figure1, should be obtained with conditions. Given the behavior of the layered rocks which are under the condition of varying water table should be ideally plastic and follow an associated flow rule; that is, the stress points of the ideally plastic soils should not exceed the yield surface. The width of collapsing block at the ground level, 5i, is only determined by the depth from the center of tunnel to ground surface. Moreover the relationship of the i and H can be described as i=k·H, where i stands for the distance between the tunnel centerline and the point of the trough inflexion, H is the depth of the center of tunnel, and k is a coefficient. By considering the shape of the collapsing block to be symmetrical with respect to the y-axis, the profile of the detaching curves should be smooth and continuous. The collapsing block can be seen as a rigid block without considering the arch effect of shallow circle tunnels.


## 4.2. Computation of Internal Energy Dissipation Rate

From the assumptions above, without considering the geometric difference in different soil layers, the parameters of the NL strength law are assumed to be the same. Due to the presence of velocity detaching line existing in the soil layer of tunnel roofs, the impending failure would slide in a limit state along with the velocity discontinuous surfaces. During the process of the impending collapse, the dissipation densities of the internal forces on the detaching surface are(10)WD=∫L1L2D˙i·w1+f′x2dx+∫L22.5iD˙i·w1+g′x2dx=∫L1L2-T+1-n-1nAf′x1/1-nPa·u˙dx+∫L22.5i-T+1-n-1nAg′x1/1-nPa·u˙dx,where L1 and L2 characterize the collapsing width of upper and lower blocks illustrated in Figure 2. y=f′(x) and y=g′(x) are the first derivative of y=f(x) and y=g(x), respectively.Figure 2
Curved failure mechanism of shallow circular tunnels.


## 4.3. Calculation of External Work

The work rate of failure block produced by weight can be calculated by integral process(11)We=ρ′∫0L1cx·u˙dx+ρ′∫L1L2fx·u˙dx+ρ∫L22.5igx·u˙dx-ρ∫02.5iSx·u˙dx+u˙·ρwL2fL2,where ρ′ is the buoyant weight per unit volume of the rocks. ρ′=ρ-ρw, in which ρ is the weight per unit volume of the rock, and ρw is the unit weight of water. The function c(x) is the function describing the circular tunnel profile. S(x) characterizes the shape of surface settlement(12)cx=H+R-R2-x2.According to the research of Osman [32], the surface settlement is defined by a Gaussian distribution curve. The area of surface settlement trough can be obtained by calculation:(13)∫-2.5i2.5iSxdx=2πiSm.The distribution of excess pore pressure which is derived from the study of Saada et al. [38] can be written as(14)u=p-pw=p-ρwh,where p is the pore water pressure at the considered point which can be obtained by a suitable method p=ruρh, ru stands for pore pressure coefficient, and pw=ρwh is the hydrostatic distribution for pore pressure; h is the vertical distance between the roof of the tunnel and the top of the failure block. So -grad⁡u can be defined as(15)-grad⁡u=ρw-ruρ.Therefore the work rate produced by seepage forces along the velocity discontinuity surface is(16)Wu=u˙∫0L1ρw-ruρcx-fL2dx+u˙∫L1L2ρw-ruρfx-fL2dx.Due to the fact that the tunnel is buried in the shallow strata, the supporting structure is unavoidable for the requirement of safety and stability. Therefore, the work rate of supporting pressure in the shallow circular tunnel is(17)Wq=Rqu˙arcsin⁡L1Rcos⁡π,where q is the supporting pressure exerting on the circumference of tunnel lining.Meanwhile, because of the external force always exerted on the underground structure, the work rate of extra force which puts on the ground surface cannot be ignored. The expressions can be written as(18)WS=2.5iσSu˙,where σS stands for the surcharge load put on the ground surface.


## 4.4. Analytic Solution for Characterizing the Collapsing Shape and Supporting Pressure

In order to describe the shape and extension of the failure collapsing block, it is essential to obtain the explicit expressions ofg(x) and f(x) by constructing an objective function consisting of the external rate of work and the rate of the internal energy dissipation,(19)Λ=WD-We-Wq-Wu-WS.So, the effective collapse mechanism can be obtained by minimizing the objective functionΛ according to the kinematic theorem of limit analysis.Then substituting (10), (11), (16), (17), and (18) into (19), the expression of objective function is given:(20)Λ=∫L1L2ψ1fx,f′x,xu˙dx+∫L22.5iψ2gx,g′x,xu˙dx-ρ′∫0L1cxu˙dx-ρw-ruρ∫0L1cx-fL2·u˙dx+ρ∫02.5iSxu˙dx-2.5iσSu˙+Rqu˙arcsin⁡L1R,in which(21)ψ1fx,f′x,x=-T+1-n-1nAf′x1/1-nPa+ρru-1fx,ψ2gx,g′x,x=-T+1-n-1nAg′x1/1-nPa-ρgx.In order to get the extremum of the objective functionΛ, it is necessary to search for the extremum values of objective functions ψ1 and ψ2. In other words, the analytical solutions of collapse dimension could be obtained by seeking for the minimum values of objective functions ψ1 and ψ2 with the method of variational principle in the realm of plasticity theory. As a consequence the expressions of ψ1 and ψ2 should be turned into two Euler’s equations through the variational method. The expressions of variational equations of ψ1 and ψ2 on stationary conditions can be written as(22)∂ψ1∂fx-∂∂x∂ψ1∂f′x=0,∂ψ2∂gx-∂∂x∂ψ2∂g′x=0.By the variational calculation the explicit forms of the two Euler’s equations for (21) can be obtained as(23)ru-1ρ-PanA1/1-nn-1-1f′x2n-1/1-nf′′x=0,-ρ-PanA1/1-nn-1-1g′x2n-1/1-ng′′x=0.Obviously, (23) are two nonlinear second-order homogeneous differential equations. By the integral calculation process the expressions of velocity discontinuity surface of lower and upper soil layers are(24)fx=k1x+θ1-ruρ1/n+h,gx=k2x+θ′ρ1/n+h′in which(25)k1=A-1/n1-ruρPa1-n/n,k2=A-1/nρPa1-n/n,where θ,θ′,h, and h′ stand for the integration constants coefficients determined by mechanical and geometric boundary conditions, respectively.From what is mentioned above, the detaching curves are supposed to be symmetric with respect to they-axis. It can be seen from Figure 2 that the shear stress on the ground surface equals zero at the point where its x-height is 2.5i:(26)τxyx=2.5i,y=0=0.Furthermore, the explicit expressions of the function of velocity discontinuity surface should fulfill other boundary conditions, such as(27)gx=2.5i=0,(28)fx=L2=gx=L2=n0H,(29)cx=L1=fx=L1,where n0 indicates the ratio of the upper height of roof collapse to the depth from the center of tunnel to ground surface H.Substituting (24) into above expressions, they turn into(30)gx=k2x-2.5i1/n,(31)h=-k1L2+θ1-ruρ1/n+n0H.For the purpose of keeping the whole curve look smooth, another boundary condition should be satisfied:(32)f′x=L2=g′x=L2.Building on this result, the expressions of the function of velocity discontinuity surface can be obtained:(33)fx=k1x-Z1/n-k1L2-Z1/n+n0H,in which(34)Z=2.5i-L2ru1-ru.On the basis of the expression of the profile of the circular tunnel, the piece of external work can be calculated by integratingc(x) over the interval 0,L1: (35)∫0L1cxdx=H+RL1-L12R2-L12-R22arcsin⁡L1R.Therefore, the objective functionΛ can be calculated, which is(36)Λ=u˙ru-1ρHn0-k1ru-1ρL2-Z1/nL2-L1-TPa2.5i-L1+nn+1Pa1-n-1Ak11/1-n+k1ru-1ρ·L2-Zn+1/n-L1-Zn+1/n+nn+1k2ρ-Pa1-n-1Ak21/1-n·L2-2.5in+1/n-ruρL2Hn0+ru-1ρH+RL1-L12R2-L12-R22arcsin⁡L1R+Rqarcsin⁡L1R-2.5iσS+2πiSmρ2.Let Λ=0; the equation can be written as(37)q=-1R·arcsin⁡L1/Rru-1ρHn0-k1ru-1ρL2-Z1/nL2-L1-TPa2.5i-L1+nn+1Pa1-n-1Ak11/1-n+k1ru-1ρ·L2-Zn+1/n-L1-Zn+1/n+nn+1k2ρ-Pa1-n-1Ak21/1-n·L2-2.5in+1/n-ruρL2Hn0+ru-1ρH+RL1-L12R2-L12-R22arcsin⁡L1R-2.5iσS+2πiSmρ2.For the purpose of getting the explicit forms of detaching curve profile consisting off(x) and g(x), the values of L1,L2, and q must be obtained by combining and solving (28), (29), and (37). Then it is not difficult to draw the shape of failure surface by (30) and (33). Meanwhile influences of different factors on the upper bound solution q can be analyzed.


## 5. Catastrophe Analysis of the Stability of Shallow Tunnel

As mentioned in the previous section, if the potential function of the system is defined by a functionV=fx1,x2, as in (38), determining the non-Morse critical point fc(x) of the potential function of the system becomes challenging. Consider(38)Jfx=∫x1x2Ffx,f′x,xdxin which the primes indicate the derivatives of the functions with respect to their subscript coordinates; namely, f′(x)=∂f(x)/∂x.In order to get the non-Morse critical pointfc(x) of the potential function of the system, it is necessary to expand the increment of function to the forms of a two-variable Taylor series in small perturbations of y. This work makes y′=f′(x),y=fx facilitate the derivation:(39)ΔJ=Jy+δy-Jy=∫x1x2Fy+δy,y′+δy′,x-Fy,y′,xdx=∫x1x2∂F∂yδy+∂F∂y′δy′+12!∂2F∂y2δy2+2∂2F∂y∂y′δy∂y′+∂2F∂y′2δy′2+13!∂3F∂y3δy3+3∂3F∂y2∂y′δy2∂y′+3∂3F∂y∂y′2δy∂y′2+∂3F∂y′3δy′3+⋯dx.Importantly, the functionJ synchronously reaches the catastrophic point when F reaches the catastrophic point. According to the Thom splitting lemma [35], F has a non-Morse critical point if the following equation is satisfied: (40)Df=0,det⁡Hf=0,where Df and det⁡(Hf) denote the gradient and determinant of the Hesse matrix of potential function fx1,x2, respectively.According to (40), the conditions of function Jy are illustrated below:(41)δJ=∫x1x2∂F∂yδy+∂F∂y′δy′dx=∫x1x2∂F∂y∂F∂y′δyδy′dx=0,δJ=∫x1x2∂2F∂y2δy2+2∂2F∂y∂y′δy∂y′+∂2F∂y′2δy′2dx=∫x1x2δyδy′∂2F∂y2∂2F∂y∂y′∂2F∂y∂y′∂2F∂y′2δyδy′dx=0.The form of catastrophic conditions of functionJy is illustrated in the following:(42)∂2Λ∂fx2-2∂∂x∂2Λ∂fx∂f′x+∂2∂x2∂2Λ∂f′x2=0.During the process of the catastrophe analysis of the stability of a shallow tunnel, the specific expressions ofg(x) and f(x) can be obtained by substituting (30) and (33) into (42); the results can be written as follows:(43)A1/nρ2n-1/nPa1-n/n2n-1x+ρ-1θ′-1/n=0,A1/nρ2n-1/nPa1-n/n2n-1x+1-ruρ-1θ-1/n=0.Given any values ofx, (43) must be satisfied. Therefore, the value of n in (43) must be 0.5 synchronously.The value ofn in this work is consistent with the general requirements specified in nonlinear Mohr envelopes general considerations. According to Baker [29], it is necessary to impose the restrictions A>0,1/2≤n≤1,T≥0, because the radius of curvature of SNLσ is less than the radius of the tangential Mohr circle [39] when n<1/2. In other words, if n<1/2, the SNLσ intersects twice with the same Mohr circle, thus contradicting the basic definition of Mohr envelopes. The original form of the HB empirical criterion was derived without using the notion of Mohr envelopes, and the extensive literature dealing with this criterion does not include the requirement n≥1/2. This requirement is satisfied in applications of this criterion to real experimental data [37].


## 6. Numerical Results and Discussions

The explicit failure surfaces of circular tunnel profiles can be drawn according to the analytical solutions of velocity discontinuity surfacesf(x) and g(x) which are shown as (30) and (33). According to the failure mechanism proposed by Qin et al. [40], the strength parameters are considered to be the same with ignoring the pore water effect made by the underground water. If the water level is under the line of the center of the tunnel, n0→1, the whole collapse mode is saturated in solely dry rock layer. Meanwhile, when the water lever is close to the ground surface, n0→0, the total failure mechanism is situated below the water level.


### 6.1. Effects of Different Parameters on Shape of Roof Collapse

In order to explore the influence of different parameters such asA,n0,ρ, and ru on the shape of roof collapse, the difference of failure mechanism of falling blocks up and down the water level should be considered. At the same time, other factors should be fixed. Table 2 shows the failure mechanism of tunnel roof corresponding to T=0.25, n=0.5, Pa=100kPa, σS=50kPa, R=5m, Sm=0.1m, and H=10m. According to Table 2, it can be obviously seen that the width of collapsing block around the circumference of tunnel lining tends to decrease with the increase of parameters A and ρ. On the contrast, the table shows totally different results about the influence of parameters ru on the width of collapsing block around the circumference of tunnel lining. Furthermore the size of the failure collapsing block varies more slightly with the change of parameter n0. According to Osman et al. [31], the width of roof collapse extended to the ground surface is merely associated with the depth H. It can be concluded that when H is fixed, the collapsing width stays a constant no matter how the other parameters vary. Moreover it can be concluded that the value of L2 tends to decrease with the increase of parameters A and n0. However, the parameter ρ makes a totally different effect on the change of the value of L2 compared with parameters A and n0. From the perspective of engineering, the position of the water table will contribute to controlling the size of collapse block. The shallower the underground water lever is buried, the larger the scope of the failure block will be during the process of excavation. In consequence, supporting system for shallow ground water table should be intensified to avoid collapse and more measurements should be taken to protect the soil saturated in the underground water.Table 2
The impending roof failure with regard to different parameters.
A n n 0 r u P a ρ T ⋅ P a L 1 L 2 k P a kN/m3 k P a m m 0.9 0.5 0.2 0.1 100 16 2.5 1.0092 7.4688 0.8 0.5 0.2 0.1 100 16 2.5 2.1009 8.0279 0.7 0.5 0.2 0.1 100 16 2.5 3.1228 8.5869 0.6 0.5 0.2 0.1 100 16 2.5 4.0889 9.1459 0.7 0.5 0.1 0.1 100 16 2.5 3.0639 9.7330 0.7 0.5 0.2 0.1 100 16 2.5 3.1228 8.5869 0.7 0.5 0.3 0.1 100 16 2.5 3.1596 7.7074 0.7 0.5 0.4 0.1 100 16 2.5 3.1852 6.9660 0.7 0.5 0.6 0.1 100 16 2.5 3.2173 5.7223 0.7 0.5 0.8 0.1 100 16 2.5 3.2341 4.6738 0.7 0.5 0.2 0.1 100 16 2.5 3.1228 8.5869 0.7 0.5 0.2 0.2 100 16 2.5 2.9906 8.5869 0.7 0.5 0.2 0.3 100 16 2.5 2.8433 8.5869 0.7 0.5 0.2 0.4 100 16 2.5 2.6775 8.5869 0.7 0.5 0.2 0.1 100 14 2.5 2.6359 8.3167 0.7 0.5 0.2 0.1 100 16 2.5 3.1228 8.5869 0.7 0.5 0.2 0.1 100 18 2.5 3.5168 8.8107 0.7 0.5 0.2 0.1 100 20 2.5 3.8430 9.0000


### 6.2. Effects of Different Parameters on Supporting Pressure

To endure the surcharge load, soil weight, and other external forces more efficiently, it is meaningful to obtain the minimal upper bound supporting pressure exerted on the tunnel lining within the framework of upper bound theorem with considering the diverse impacts which different parameters have. Parametric analysis is conducted, such asσS,R,Sm,T with respect to the depth H, which are illustrated in Figure 3.Figure 3
Effects of different parameters on supporting force of shallow tunnels.
(a)
Effects ofSm on supporting pressure (b)
Effects ofσS on supporting pressure (c)
Effects ofT on supporting pressure (d)
Effects ofR on supporting pressureAccording to the figure, the burial depthH has a positive correlation with the supporting force consistently. When the burial depth H is fixed, different parameters have different influences on the supporting pressure corresponding to Pa=100kPa, n=0.5, A=0.7, ru=0.1, n0=0.5, and ρ=16kN/m3. It can be concluded that the supporting force tends to decrease with the increase of the radius of the circular tunnels and the value of maximum surface settlement, while the influences of surcharge load and the parameter T in NL criterion on supporting pressure present the opposite trend.The figure shows that the value of maximum surface settlementSm makes a bigger influence on the supporting force. Meanwhile the value of the supporting force varies more slightly with the change of the parameters T and σS. The effect of the radius of the circular tunnels R becomes bigger as the tunnels buries shallower. From the perspective of engineering, the shallow-buried circular tunnels with heavier ground load and smaller radius need larger supporting forces and supporting system of shallower tunnels should be intensified to keep the stability.


## 6.1. Effects of Different Parameters on Shape of Roof Collapse

In order to explore the influence of different parameters such asA,n0,ρ, and ru on the shape of roof collapse, the difference of failure mechanism of falling blocks up and down the water level should be considered. At the same time, other factors should be fixed. Table 2 shows the failure mechanism of tunnel roof corresponding to T=0.25, n=0.5, Pa=100kPa, σS=50kPa, R=5m, Sm=0.1m, and H=10m. According to Table 2, it can be obviously seen that the width of collapsing block around the circumference of tunnel lining tends to decrease with the increase of parameters A and ρ. On the contrast, the table shows totally different results about the influence of parameters ru on the width of collapsing block around the circumference of tunnel lining. Furthermore the size of the failure collapsing block varies more slightly with the change of parameter n0. According to Osman et al. [31], the width of roof collapse extended to the ground surface is merely associated with the depth H. It can be concluded that when H is fixed, the collapsing width stays a constant no matter how the other parameters vary. Moreover it can be concluded that the value of L2 tends to decrease with the increase of parameters A and n0. However, the parameter ρ makes a totally different effect on the change of the value of L2 compared with parameters A and n0. From the perspective of engineering, the position of the water table will contribute to controlling the size of collapse block. The shallower the underground water lever is buried, the larger the scope of the failure block will be during the process of excavation. In consequence, supporting system for shallow ground water table should be intensified to avoid collapse and more measurements should be taken to protect the soil saturated in the underground water.Table 2
The impending roof failure with regard to different parameters.
A n n 0 r u P a ρ T ⋅ P a L 1 L 2 k P a kN/m3 k P a m m 0.9 0.5 0.2 0.1 100 16 2.5 1.0092 7.4688 0.8 0.5 0.2 0.1 100 16 2.5 2.1009 8.0279 0.7 0.5 0.2 0.1 100 16 2.5 3.1228 8.5869 0.6 0.5 0.2 0.1 100 16 2.5 4.0889 9.1459 0.7 0.5 0.1 0.1 100 16 2.5 3.0639 9.7330 0.7 0.5 0.2 0.1 100 16 2.5 3.1228 8.5869 0.7 0.5 0.3 0.1 100 16 2.5 3.1596 7.7074 0.7 0.5 0.4 0.1 100 16 2.5 3.1852 6.9660 0.7 0.5 0.6 0.1 100 16 2.5 3.2173 5.7223 0.7 0.5 0.8 0.1 100 16 2.5 3.2341 4.6738 0.7 0.5 0.2 0.1 100 16 2.5 3.1228 8.5869 0.7 0.5 0.2 0.2 100 16 2.5 2.9906 8.5869 0.7 0.5 0.2 0.3 100 16 2.5 2.8433 8.5869 0.7 0.5 0.2 0.4 100 16 2.5 2.6775 8.5869 0.7 0.5 0.2 0.1 100 14 2.5 2.6359 8.3167 0.7 0.5 0.2 0.1 100 16 2.5 3.1228 8.5869 0.7 0.5 0.2 0.1 100 18 2.5 3.5168 8.8107 0.7 0.5 0.2 0.1 100 20 2.5 3.8430 9.0000


## 6.2. Effects of Different Parameters on Supporting Pressure

To endure the surcharge load, soil weight, and other external forces more efficiently, it is meaningful to obtain the minimal upper bound supporting pressure exerted on the tunnel lining within the framework of upper bound theorem with considering the diverse impacts which different parameters have. Parametric analysis is conducted, such asσS,R,Sm,T with respect to the depth H, which are illustrated in Figure 3.Figure 3
Effects of different parameters on supporting force of shallow tunnels.
(a)
Effects ofSm on supporting pressure (b)
Effects ofσS on supporting pressure (c)
Effects ofT on supporting pressure (d)
Effects ofR on supporting pressureAccording to the figure, the burial depthH has a positive correlation with the supporting force consistently. When the burial depth H is fixed, different parameters have different influences on the supporting pressure corresponding to Pa=100kPa, n=0.5, A=0.7, ru=0.1, n0=0.5, and ρ=16kN/m3. It can be concluded that the supporting force tends to decrease with the increase of the radius of the circular tunnels and the value of maximum surface settlement, while the influences of surcharge load and the parameter T in NL criterion on supporting pressure present the opposite trend.The figure shows that the value of maximum surface settlementSm makes a bigger influence on the supporting force. Meanwhile the value of the supporting force varies more slightly with the change of the parameters T and σS. The effect of the radius of the circular tunnels R becomes bigger as the tunnels buries shallower. From the perspective of engineering, the shallow-buried circular tunnels with heavier ground load and smaller radius need larger supporting forces and supporting system of shallower tunnels should be intensified to keep the stability.


## 7. Summary and Conclusions

On the basis of previous work which has focused the efforts on the collapse mechanism in deep tunnels, a new curved failure mechanism of circular shallow tunnel considering the joint effects of surface settlement and seepage forces is proposed to estimate the stability of tunnel roof under a limit state. Furthermore this work proved the validity of the NL strength function in predicting the stability of the tunnel roof by comparing with HB failure criterion. With NL failure criterion and FCT theory the numerical solution for the shape of collapse mechanism is obtained by setting an objective function consisting of energy dissipation rate and the external work rate. Some conclusions are drawn from above:(1)
With the increase of the shear strength controlling coefficientA and unit weight ρ in NL failure criterion, the value of the width of collapsing block around the circumference of tunnel lining starts to decrease. Meanwhile pore pressure coefficient ru shows contrary law on it. Moreover, the effects of the position of the water table on the size of the failure collapsing block in shallow-buried tunnel tend to be small.(2)
The supporting forces tend to decrease with the increase of the radius of the circular tunnels and the surface settlement while the surcharge load and the parameterT both have a positive correlation with the scope of supporting forces. Furthermore, the size of the failure collapsing block and the supporting forces both increase as the tunnels bury deeper. Particularly the upper width of collapsing block is proportional to the depth H. However the effect of the value of radius of circular tunnels R gets smaller as the tunnel buries deeper.



---

*Source: 4820716-2016-08-03.xml*