Collapsed Shape of Shallow Unlined Tunnels Based on Functional Catastrophe Theory

This paper investigates the collapse mechanisms and possible collapsing block shapes of shallow unlined tunnels under conditions of plane strain. The analysis is performed following the framework from a branch of catastrophe theory, functional catastrophe theory. First, the basic principles of functional catastrophe theory are introduced. Then, an analytical solution for the shape curve of the collapsing block of a shallow unlined tunnel is derived using functional catastrophe theory based on the nonlinear HoekBrown failure criterion. The effects of the rock mass parameters of the proposed method on the shape and weight of the collapsing block are examined. Moreover, a critical cover depth expression to classify deep and shallow tunnels is proposed. The analytical results are consistent with those obtained by numerical simulation using the particle flow code, demonstrating the validity of the proposed analytical method.The obtained formulas can be used to predict the height and width of the collapsing block of a shallow unlined tunnel and to provide a direct estimate of the overburden on the tunnel lining. The obtained formulas can be easily used by tunnel engineers and researchers due to their simplicity.


Introduction
Shallow tunneling in soft ground inevitably results in ground movements.These movements may deform, rotate, distort, and irrevocably damage adjacent surface buildings in urban environments.In general, there are five different categories of tunnel collapses according to different factors as follows: daylight collapse, underground collapse, rock burst, inrush of water, and portal collapse [1].In the event of a daylight collapse (i.e., a tunnel collapse that reaches the ground surface), the ground is unraveled to the surface and the propagation of the failure to the surface can be extremely quick.Some examples of daylight collapse induced by shallow tunneling are shown in Figure 1 [2,3].Hence, tunnel stability is a main concern for the design and construction of shallow tunnels.The inherent uncertainties related to the natural surrounding rock may complicate the precise prediction of tunnel collapse.
Tunnel stability analyses have been investigated by several scholars using various analytical approaches.The limit equilibrium method is widely used in the theoretical analysis of tunnel stability.Based on the failure mechanism of a twodimensional log-spiral shaped slip plane, Murayama et al. [4] obtained the minimal support pressures that are needed to guarantee the tunnel face stability.Krause [5] calculated the minimal support pressures with a semicircular and spherical failure mechanism by using the limit equilibrium method.Horn [6] introduced a three-dimensional wedge model, which assumed a triangular slip wedge-shaped failure mechanism, to analyze the stability of shallow tunnels.Anagnostou and Kovári [7] subsequently used a similar model to Horn and analyzed the eroding influence of slurry infiltration on tunnel face stability in a homogeneous stratum.Broere [8] extended the wedge model to a layered stratum and performed a sensitivity analysis of different parameters to minimal support pressures.The limit analysis method, which is based on the upper and lower bound theorems of plasticity, is also an important method used to analyze tunnel stability.Davis et al. [9] developed failure mechanisms with assumed three-and four-variable plane strain conditions and obtained the upper-bound solutions of the limit support pressures.Sloan and Assadi [10] proposed a promising numerical approach (called finite element limit analysis) to analyze the undrained stability of shallow tunnels in the case of a plane strain circular tunnel in cohesive soil whose shear strength varied linearly with cover depth using linear programming techniques.Lyamin and Sloan [11] subsequently considered the stability of a plane strain circular tunnel in a cohesivefrictional soil using a more efficient nonlinear programming technique.More recently, several kinematical approaches based on continuous velocity fields in limit analysis theory were proposed.Osman et al. [12] developed upper-bound calculations of collapse loads in tunnels using a distributed shear deformation mechanism.They assumed that the soil deforms compatibly following a Gaussian distribution within the boundaries of this deformation mechanism and is rigid outside this mechanism.Klar et al. [13] suggested a new kinematical approach in limit analysis theory for the 2D and 3D stability analysis of circular tunnels in a purely cohesive soil.The velocity fields in the 2D and 3D stability analyses are based on the analytical works of Verruijt and Booker [14] and Sagaseta [15], respectively.Fraldi and Guarracino [16] and Yang and Huang [17] obtained the collapsed shape of a deep tunnel and shallow tunnel, respectively, based on the Hoek-Brown failure criterion by combining the upper bound of the limit analysis and the calculus of variations.
Because the deformation and failure of tunnel surrounding rock are characterized by nonlinear and discontinuous phenomena, catastrophe theory is suitable for investigating the behaviors of tunnel collapse.Catastrophe theory arises from the mathematical work of Thom [18], whose famous book, Structural Stability and Morphogenesis, was published in French in 1972 and in English in 1975.Catastrophe theory became widely used due to the efforts of Zeeman [19] and Arnold and Afraimovich [20].Recently, several researchers applied catastrophe theory to analyze stability problems in geology and geomechanics [21][22][23].Pan et al. [21] established a fold catastrophe model of a tunnel rock burst and reported that the occurrences of a rock burst are related to both the ratio of the elastic modulus to the descendent modulus of the rock mass and the rock mass crack growth degree.Based on fold catastrophe theory, Miao et al. [22] presented a model of a seepage flow system to analyze the dynamical behavior of water or gas flows in damaged rock, and the influences of different key parameters on the stability of seepage flow systems were obtained.Using catastrophe theory and the discontinuous deformation analysis (DDA) method, Xia et al. [23] studied the stability of tunnel surrounding rock and obtained the safety factors for tunnels, which can be adopted in real-life projects to guarantee the safety of the tunnels.
Although catastrophe theory has been proven to be a promising analysis tool for discontinuous phenomena, the current research on tunnel collapse is generally based on elementary catastrophe theory (ECT).In ECT, the potential function of the system is one of the eleven elementary functions defined by polynomial functions, whose independent variables are the state variable and control variable.Moreover,  and the corresponding potential function is expressed in a complex mathematical form, such as a functional of (), obtaining the catastrophic value   () of the state variable is challenging when using ECT.Functional catastrophe theory (FCT), which was proposed by Du [24] and in which the potential function is defined as a functional instead of an elementary function, provides a method to obtain the catastrophic value   () of the state variable.Du used FCT to successfully solve several complex physical and economic problems.The details of the solution procedures for ECT and FCT are described in Section 2, and the differences between the two methods are also presented.To investigate the shape curve () of the collapsing blocks in shallow unlined tunnels based on catastrophe theory, the potential function of the system is considered to be a functional based on ().Therefore, FCT is suitable for analyzing this problem.Using the FCT, Zhang et al. [25] studied the collapse mechanisms of deep tunnels and derived possible shapes of the collapsing block in deep tunnels.However, the boundary condition of ground surface has influence on the collapse mechanisms of shallow tunnels.The collapse of deep tunnel has not reached the ground surface (called underground collapse), whereas the collapse of shallow tunnels has always reached the ground surface (called daylight collapse).This paper investigates the collapse shape and mechanism of a shallow unlined tunnel under the conditions of plane strain in the FCT framework.First, the basic principles of FCT are introduced in Section 2.Then, an analytical solution for the shape curve of the collapsing block of a shallow unlined tunnel is derived based on the nonlinear Hoek-Brown failure criterion in Section 3. The effects of the rock mass parameters on the shape and weight of the collapsing block are studied.A critical cover depth expression to classify deep and shallow tunnels is proposed in Section 4. Finally, in Section 5, the analytical results are compared with those obtained by numerical simulation using the particle flow code.

Overview of Catastrophe Theory
2.1.Fundamental Theories for Catastrophe Theory.The research objects of catastrophe theory are the potential functions , whose independent variables are the state variable and control variable.Catastrophe theory analyzes the sudden changes of state variables resulting from a slow, smooth, and small change in one or more control variables.There are three fundamental theories behind catastrophe theory, namely, the implicit function theorem, Morse lemma, and Thorn splitting lemma [18].These theories state that the catastrophe of state variables would not occur (i.e., the system is stable) under the condition that the potential function is at the Morse critical point (or nonsingular), whereas the catastrophe of the state variables would occur (i.e., the system is unstable) under the condition that the potential function is at the non-Morse critical point (or singular).Therefore, it is important for both ECT and FCT to obtain the non-Morse critical point (or singular) of the potential function.

ECT.
In ECT, the potential function of the system is one of the eleven elementary functions [ = ( 1 ,  2 )] defined by the polynomial functions shown in Table 1.The non-Morse critical point ( 1 ,  2 ) of the potential function of the system is determined using the following [18]: where  and det() denote the gradient and determinant of the Hesse matrix of potential function ( 1 ,  2 ), respectively.

FCT.
If the potential function of the system is defined by a functional [ = [()]], as in ( 2), determining the non-Morse critical point   () of the potential function of the system becomes challenging.Consider where the primes indicate the derivatives of the functions with respect to their subscript coordinates; that is,   () = ()/.
To obtain the non-Morse critical point   () of the potential function of the system, the increment of functional should first be expanded in the forms of a two-variable Taylor series in small perturbations of .Here, we make the following substitutions [ = (),   =   ()] to facilitate the derivation: where ,  2 ,  3 , and  4  are the first-, second-, third-, and fourth-order variations of the functional [], respectively.
Importantly, the functional  synchronously reaches the catastrophic point when  reaches the catastrophic point.According to the Thom splitting lemma [18],  has a non-Morse critical point if the following equation is satisfied: According to ( 5), ( 6) can be easily obtained and are also the non-Morse critical point of functional [].Therefore, the necessary conditions of functional [] catastrophe are as shown in the following: The comparison of ECT and FCT is shown in Table 2.With a subsection integral, the two specific forms of

Solution Procedure for the Collapse Shape of a Shallow
Unlined Tunnel Based on FCT.First, we must obtain the potential function of the system [, (),   ()] to determine [, (),   ()].() is an unknown function that outlines the collapsing block of a shallow unlined tunnel.Second, (7) are used to calculate a group of differential equations related to ().By integrating these differential equations, we can obtain the shape curve [()] of the collapsing block of a shallow unlined tunnel.
Finally, the unknown constants in () can be determined by the boundary transversality conditions and geometric compatibility equations.

Catastrophe Analysis of the Stability of a Shallow Unlined Tunnel
For this study, to estimate the stability of a shallow unlined tunnel, we focus on determining the shape and dimensions of the collapsing block.This block can actually collapse from the roof of the tunnel.By applying the upper-bound theorem, Yang and Huang [17] proposed a mechanical model for shallow unlined tunnels and obtained the upper solution of the collapsing block shape based on the Hoek-Brown failure criterion.The current study analyzes the collapse of shallow unlined tunnels with FCT based on Yang's mechanical model (Figure 2).Several assumptions are made when solving the proposed problem via FCT.First, we consider only the gravity field and disregard the tectonic stress field for consistency with the relative reference (Yang and Huang [17]).The rock mass behavior is elastic-perfectly plastic.The yield surface is convex.Moreover, the plastic deformation rate can be derived from the yield function according to an associated flow rule.The geometry changes in the collapsing block can be regarded as insignificant at the onset of the collapse.Hence, the problem can be considered as plane strain.Given that the potential function along the failure surface is generated by normal and shear stresses, we investigate the rock mass using the Hoek-Brown failure criterion proposed by Hoek and Brown [26], which is expressed as follows: where   and   are the normal and shear stresses, respectively,  and  are the material parameters characterizing the rock mass, and   and   are the uniaxial compressive and tensile strengths of the rock mass, respectively.In Figure 2,  1 is the half-width of the upper edge of the collapsing block;  2 is the half-width of the lower edge of the collapsing block;  is the distance from the tunnel roof to the ground surface;  is the weight per unit volume of the rock mass;  is the thickness of the plastic detaching zone; () is a known function that describes the shape of a circular tunnel: In accordance with the formula presented by Yang and Huang [17], the potential function of the system is formulated by calculating the strain energy of the internal forces and the work of the applied loads on the detaching zone.
By following a purely geometrical line of reasoning and referencing Figure 2, the plastic strain components can be expressed as follows: At an impending collapse, the strain energy of the internal forces on the detaching zone   is and the work of the applied loads on the detaching zone   is In the detaching zone, the potential function induced by normal and shear stresses can be expressed as follows: where the function  is The key task in catastrophic state analysis is to obtain the specific expression of () using (7).By substituting ( 14) into (7), we obtain the explicit forms of the group of () differential equations as follows: By integrating (15), we compute the detaching curve () as follows: where  and ℎ are two integration constants.By substituting ( 17) into ( 16), we obtain the following: Given any value of , the result of (18) must always be zero.Thus, the value of  must be 0.5.Based on the value of , the form of () is generated as Equation ( 19) contains two unknown constants,  and ℎ, which can be determined by the boundary conditions.
First, the boundary condition of the ground surface is stress-free, which requires that the shear component of stress vanishes along this plane.According to Fraldi and Guarracino [16], the tangent line of the detaching curve  = () in  =  1 must be perpendicular to the axis of  =  1 (see Figure 3).Consider By substituting ( 20) into (19), the form of () is simplified as Second, Figure 3 shows that two geometric equations must be satisfied as follows: By substituting ( 21) into ( 22), we express the unknown constants in the expression of (), and the explicit form of the geometric equation is as follows: Finally, the boundary point  of the collapsing block moves along the tunnel boundary in Figure 3.The boundary conditions are derived from the study conducted by G. Q. Wang and T. Wang [27].To meet those conditions, (14) should satisfy By substituting ( 14) into ( 25), we obtain the explicit form of the transversality condition as follows: By simplifying (26), we can simplify the transversality condition into By substituting ( 27) into ( 24), we derive the implicit function equation of  2 as The value of  2 can then be calculated using (28).We can determine the value of  1 by substituting this value into (27).
We obtain the final form of the detaching curve () based on  1 and  2 .The shape of the failure surface can be illustrated according to the following equation: Equation (29) shows that the collapsing block is parabolic.Moreover, we can compute the overall weight of the collapsing block per unit length  as follows:

Influence of Rock Mass Parameters on the Shape and
Weight of the Collapsing Block.We can determine the shape curve of the collapsing block of a shallow unlined tunnel using (29).To investigate the influence of rock mass parameters on these shape curves, we refer to Figure 4, where different shapes of the collapsing block are plotted.The weight of the collapsing block under various parameters is plotted in Figure 5 based on (30).The weight of the collapsing block increases with increases in ,   ,   , and the / ratio (Figures 5(b), 5(c), and 5(d)).However, it decreases with increases in  (Figure 5(a)) when other parameters are fixed.The results directly estimate the overall burden on the tunnel lining and are in accordance with those observed in actual applications with regard to surrounding rock.Hence, our findings not only provide a reference to prevent instability and failure in shallow unlined tunnels but also guide reasonable construction.

Critical Cover Depth to Classify Deep and Shallow Tunnels.
Given the significant differences between shallow and deep tunnels in terms of their mechanical characteristics and failure mechanisms, establishing a critical cover depth is important because it can distinguish shallow and deep tunnels.
To express the critical cover depth, we must calculate the value of  2 by equating  1 with 0 in (27).Consider Based on this definition, the critical cover depth  cr can be derived from (28), which is expressed as follows: Equation (32) shows that the critical cover depth is not only determined by tunnel dimensions but is also associated with the mechanical parameters of the surrounding rock, including ,   ,   , and .Therefore, the tunnel can be considered shallow if the inequality in (33) is satisfied and deep if the inequality in (34) is satisfied: We plot the influences of the rock mass parameters on the critical cover depth of a shallow tunnel in Figure 6.The critical cover depth of a shallow tunnel increases with increases in the tunnel radius  and  (Figures 6(a

Comparison with Numerical Analysis
To validate the method proposed in this study, the analytical solution is compared to the numerical simulation with respect to the shape curve of the collapsing block based on the parameters of the surrounding rock.These results are examined based on the Hoek-Brown failure criterion.To investigate the stability of a shallow unlined tunnel, we used the two-dimensional distinct element program PFC2D [28], which simulates the mechanical behavior of a material by representing it as an assemblage of circular particles that can be bonded to one another.Such elastic properties can be input directly into a continuum model.In PFC2D, however, the mechanical behaviors of the assemblage are dominated by the microproperties of the particles and the bonds between them.These microproperties cannot be generated through laboratory tests.Thus, the relationship between the microand macroproperties of the rock mass should be determined prior to the simulation of the tunnel excavation.The microparameters input into the numerical simulation can be identified through biaxial tests before the problem is modeled [29].

Biaxial Tests to Determine the Microparameters of the Rock
Mass.Table 3 includes a list of rock mass parameters (  ,   , , , and ) employed for the validation.According to the analytical result of (29), the shape curve of the collapsing block of the tunnel is () = 2.22( − 1.03) 2.00 .In the following section, PFC2D is used to simulate that shape curve using the same rock mass parameters.As noted in the previous section, the input parameters of a rock mass are microparameters ( -bond ,  -bond ,   ,   ,   , and ) in the simulation.However, the macroproperties of the rock mass that correspond to the microscopic parameters of the rock mass are related to the Mohr-Coulomb criterion (, , and ) rather than to the Hoek-Brown criterion (  ,   , , , and ) in the biaxial tests.Table 3 depicts the transformation in the collapsing block of a standard rock mass characterized by the Hoek-Brown and equivalent Mohr-Coulomb parameters according to Hoek and Brown [30].Table 3 indicates that the target values of the cohesive force and of the friction angle of rock mass are 0.017 MPa and 13.34 ∘ , respectively, in the PFC2D biaxial tests.By referencing the selection criteria for microscopic parameters summarized by Zhou et al. [31] and Fang [32], we conducted a set of biaxial tests with the associated stress-strain behaviors under confining pressures of 0.004, 0.006, and 0.015 MPa. Figure 7 presents the model for the biaxial tests.Table 4 lists the input parameters for these tests.Figure 8 displays examples of the biaxial test results under different levels of confining pressures.Mohr's stress circles and the envelope curve are illustrated in parameter calibration, and the cohesive force and friction angle are 0.016 MPa and 12.42 ∘ , respectively, as shown in Figure 9.These values are highly similar to the target values listed in Table 4.

Numerical Model.
To ensure that the analysis results reflect the relevant characteristics on the ground, the numerical model must consider the boundary effects.The width and height of the numerical model are 56 and 30 m, respectively, and its cover depth is 5 m (Figure 10).The tunnel diameter is 10 m.The model soil is composed of 26,191 particles, with minimum and maximum radii of 0.050 and 0.075 m, respectively.After the particles settle, gravity is applied to   the model.The tunnel is excavated by removing the particles from within the tunnel area.

Results and Discussion.
To compare the numerical simulation and analytical results in terms of the shape and dimensions of the collapsing block, the fitting curve (the solid line in Figure 11) from the numerical simulation is generated.The collapsing block is rather asymmetric due to the limitation in the numerical simulation software PFC2D.Thus, we first analyze the right side of the collapsing block, followed by the symmetry to the left.Table 5 shows the detailed fitting results, which suggest that the percentage difference is quite small (the maximum is 8.74%).Hence, the proposed analytical solution of (29) can accurately predict the shape of the collapsing block.

Conclusions
The potential function based on FCT is presented in the form of a functional, whereas the potential function in ECT is expressed as an elementary function.The former solves complex problems more effectively than the latter.In this study, FCT is employed to investigate the collapse mechanisms and possible collapsing block shapes of shallow unlined tunnels for the first time.The following conclusions can be drawn.
(1) Based on the nonlinear Hoek-Brown failure criterion and FCT, we derive an analytical solution for the

Figure 1 :
Figure 1: Some examples of daylight collapse induced by shallow tunneling.

2 Figure 3 :
Figure 3: Boundary condition of the collapsing block.

𝐿 1 and 𝐿 2
of the collapsing block increase inversely with the values of  (Figure 4(a)) and  (Figure 4(b)).In contrast,  1 and  2 increase with the values of   and   (Figure 4(c)).High values of  and  and low values of   and   indicate that the collapsing blocks are small.

Figure 5 :
Figure 5: Weight of the collapsing block versus rock mass parameters.

Figure 6 :
Figure 6: Critical cover depth of the shallow unlined tunnel versus rock mass parameters.

Figure 7 :
Figure 7: Particle sample from a biaxial test.

Table 1 :
Potential functions used in ECT.

Table 2 :
Comparison of ECT and FCT.

Table 4 :
PFC parameters of the surrounding rock after calibration.

Table 5 :
Comparison of the numerical simulation and analytical results.Equation Results Gain coefficient Offset on  axis Power Adj. -square Result of curve fitting () =  ( − )  2.15 ( − 0.94) 1.93