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

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.


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 threedimensional 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][8][9] used to evaluate the stability problems has been replaced by the nonlinear criterion [10][11][12][13][14][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][17][18] to calculate the energy dissipation rate and the external work rate accurately.In particular on the basis of nonlinear yield criterion [19][20][21][22][23][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.

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 [ = ( 1 ,  2 )] defined by the polynomial functions shown in Table 1.
However, the total potential of the system always has a complex mathematical form, such as the state function (), 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 secondorder variations of the total potential.In a similar way, the critical value   () of the state variable is obtained when catastrophe occurs.As mentioned above, some researchers have investigated the shape curve () 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.

Limit Analysis on the Basis of Nonlinear
Failure Criterion where   stands for atmospheric pressure, {, , } 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 {, , } have clear physical significance.In particular,  is a scale parameter controlling the magnitude of shear strength;  is a shift parameter controlling the location of the envelope on the  axis;  NL ≡    is the tensile strength, with  representing a nondimensional tensile strength;  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   is positive, one has So the plastic strain rate can be written as follows: where  is a scalar parameter.According to Fraldi and Guarracino [19], the plastic strain rate components can be written in the form where  is the thickness of the plastic detaching zone.A dot denotes differentiation with respect to time and a prime with respect to ; that is, In order to enforce compatibility, from ( 3) and ( 4) it follows that And the normal component of stress can be written as 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, Ḋ , results: where ε  is normal plastic strain, γ  is shear plastic strain rates,  is the thickness of the plastic detaching zone, and u is the velocity of the collapsing block.

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: where   is the stress tensor and ε  is the strain rate in the kinematically admissible velocity field, respectively.  is the limit load exerted on the boundary surface. is the length of velocity discontinuity,   is the body strength which is caused by weight,  is the volume of the plastic zone, and V  is the velocity along the velocity discontinuity surface.

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 curves  = () and  = () 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.

4.1.
Assumptions for the Analysis.The failure shape of the collapsing block, as illustrated in Figure 1, 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, 5, is only determined by the depth from the center of tunnel to ground surface.Moreover the relationship of the  and  can be described as  =  ⋅ , where  stands for the distance between the tunnel centerline and the point of the trough inflexion,  is the depth of the center of tunnel, and  is a coefficient.By considering the shape of the collapsing block to be symmetrical with respect to the -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.

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 where  1 and  2 characterize the collapsing width of upper and lower blocks illustrated in Figure 2.  =   () and  =   () are the first derivative of  = () and  = (), respectively.

Calculation of External Work.
The work rate of failure block produced by weight can be calculated by integral process where   is the buoyant weight per unit volume of the rocks.  =  −   , in which  is the weight per unit volume of the rock, and   is the unit weight of water.The function () is the function describing the circular tunnel profile.() characterizes the shape of surface settlement 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: The distribution of excess pore pressure which is derived from the study of Saada et al. [38] can be written as where  is the pore water pressure at the considered point which can be obtained by a suitable method  =   ℎ,   stands for pore pressure coefficient, and   =   ℎ is the hydrostatic distribution for pore pressure; ℎ is the vertical distance between the roof of the tunnel and the top of the failure block.So − grad  can be defined as Therefore the work rate produced by seepage forces along the velocity discontinuity surface is 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 where  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 where   stands for the surcharge load put on the ground surface.

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 of () and () by constructing an objective function consisting of the external rate of work and the rate of the internal energy dissipation, 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: in which 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 By the variational calculation the explicit forms of the two Euler's equations for (21) can be obtained as 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 in which where ,   , ℎ, and ℎ  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 the -axis.It can be seen from Figure 2 that the shear stress on the ground surface equals zero at the point where its -height is 2.5: Furthermore, the explicit expressions of the function of velocity discontinuity surface should fulfill other boundary conditions, such as where  0 indicates the ratio of the upper height of roof collapse to the depth from the center of tunnel to ground surface .Substituting (24) into above expressions, they turn into For the purpose of keeping the whole curve look smooth, another boundary condition should be satisfied: Building on this result, the expressions of the function of velocity discontinuity surface can be obtained: in which On the basis of the expression of the profile of the circular tunnel, the piece of external work can be calculated by integrating () over the interval [0,  1 ]: Therefore, the objective function Λ can be calculated, which is Let Λ = 0; the equation can be written as For the purpose of getting the explicit forms of detaching curve profile consisting of () and (), the values of  1 ,  2 , and  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  can be analyzed.

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 function [ = ( 1 ,  2 )], as in (38), determining the non-Morse critical point   () of the potential function of the system becomes challenging.
in which the primes indicate the derivatives of the functions with respect to their subscript coordinates; namely,   () = ()/.In order to get the non-Morse critical point   () 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 .This work makes [  =   (),  = ()] facilitate the derivation: Importantly, the function  synchronously reaches the catastrophic point when  reaches the catastrophic point.According to the Thom splitting lemma [35],  has a non-Morse critical point if the following equation is satisfied: where  and det() denote the gradient and determinant of the Hesse matrix of potential function ( 1 ,  2 ), respectively.
According to (40), the conditions of function [] are illustrated below: The form of catastrophic conditions of function [] is illustrated in the following: During the process of the catastrophe analysis of the stability of a shallow tunnel, the specific expressions of () and () can be obtained by substituting (30) and ( 33) into (42); the results can be written as follows: Given any values of , (43) must be satisfied.Therefore, the value of  in (43) must be 0.5 synchronously.
The value of  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 { > 0, 1/2 ≤  ≤ 1,  ≥ 0}, because the radius of curvature of  NL () is less than the radius of the tangential Mohr circle [39] when  < 1/2.In other words, if  < 1/2, the  NL () 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  ≥ 1/2.This requirement is satisfied in applications of this criterion to real experimental data [37].

Numerical Results and Discussions
The explicit failure surfaces of circular tunnel profiles can be drawn according to the analytical solutions of velocity discontinuity surfaces () and () 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,  0 → 1, the whole collapse mode is saturated in solely dry rock layer.Meanwhile, when the water lever is close to the ground surface,  0 → 0, the total failure mechanism is situated below the water level.

Effects of Different Parameters on Shape of Roof Collapse.
In order to explore the influence of different parameters such as ,  0 , , and   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  = 0.25,  = 0.5,   = 100 kPa,   = 50 kPa,  = 5 m,   = 0.1 m, and  = 10 m.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  and .On the contrast, the table shows totally different results about the influence of parameters   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  0 .According to Osman et al. [31], the width of roof collapse extended to the ground surface is merely associated with the depth .It can be concluded that when  is fixed, the collapsing width stays a constant no matter how the other parameters vary.Moreover it can be concluded that the value of  2 tends to decrease with the increase of parameters  and  0 .However, the parameter  makes a totally different effect on the change of the value of  2 compared with parameters  and  0 .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.

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   , ,   ,  with respect to the depth , which are illustrated in Figure 3.According to the figure, the burial depth  has a positive correlation with the supporting force consistently.When the burial depth  is fixed, different parameters have different influences on the supporting pressure corresponding to   = 100 kPa,  = 0.5,  = 0.7,   = 0.1,  0 = 0.5, and  = 16 kN/m 3 .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  in NL criterion on supporting pressure present the opposite trend.
The figure shows that the value of maximum surface settlement   makes a bigger influence on the supporting force.Meanwhile the value of the supporting force varies more slightly with the change of the parameters  and   .The effect of the radius of the circular tunnels  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.

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

Figure 1 :
Figure 1: Potential falling blocks with consideration of the effects of surface settlement and pore water pressure.

Figure 2 :
Figure 2: Curved failure mechanism of shallow circular tunnels.

Table 2 :
The impending roof failure with regard to different parameters.