Assessment of the Stability of an Unlined Rectangular Tunnel with an Overload on the Ground Surface

City tunnels are often constructed at shallow depths, and tunnel failure may be initiated by overloads resulting from surrounding buildings, structures, heavy-haul trailers, and other installations. Although several works have been reported on tunnel stability, stability numbers have mainly been obtained for cases with fully cohesive soils. Moreover, little information has been presented about the influence of overloads on the failure patterns for unlined rectangular tunnels. -is paper uses upper-bound finite element methods to assess the stability of an unlined rectangular tunnel in cohesive-frictional soils with an overload acting on the ground surface. A complete set of dimensionless parameters covering the tunnel size and shallow tunnel depth and MohrCoulomb material parameters are determined to obtain the dimensionless overload. In addition, failure modes that are similar to slip line fields are acquired. A failure mechanism that may cause base heave is proposed in this paper to improve the accuracy of the results. -ese failure patterns are more complex for cases with larger dimensionless depth, larger internal friction angle, and smaller dimensionless unit weight. Compared with the rigid-block mechanisms from the upper-bound rigid-block analysis method, these computed failure mechanisms are better suited for rectangular tunnel stability analysis.


Introduction
Due to the rapidly increasing demand for urban construction and the improvement in social welfare, tunnels have been widely constructed in underground engineering applications, such as subways, electric utilities, and city infrastructure. ese tunnels are often constructed at shallow depths, and tunnel failures may be initiated by overloads resulting from surrounding buildings, structures, heavyhaul trailers, and other installations. To avoid economic losses and adverse social impacts resulting from tunnel construction accidents, civil engineers must accurately assess tunnel stability.
Various methods, such as finite element limit analysis methods [1][2][3][4][5][6][7][8][9][10][11], assumed-failure-mechanism-based methods [12][13][14][15][16][17][18][19][20][21][22][23][24], boundary-element-based methods [25,26], element-free-based methods [27], and laboratory and centrifuge test methods [28][29][30][31], have been used to evaluate the stability of underground excavations. Due to its simplicity and effectiveness, the finite element limit analysis method, which can handle complicated geometries and complex loads without assuming any failure mechanism, is one of the most effective stability evaluation methods. Because circular tunnels can provide good stability and are convenient to construct, they are commonly used in practical engineering applications. Circular tunnel stability has been widely analyzed [32][33][34][35][36]. Compared to circular tunnels, rectangular tunnels have a higher degree of utilization of underground space and a lower volume of earthwork excavation; thus, a rectangular tunnel is often a valid option due to the development of advanced tunneling machines. ese tunnels are highly different in terms of the stability with changes in the tunnel span; moreover, rectangular tunnels exhibit different failure mechanisms with circular tunnels. Several studies have been performed to determine the stability of single or dual square tunnels. For undrained soils, Assadi and Sloan [37] and Wilson et al. [38,39] analyzed the failure patterns of a square tunnel using the upper-bound finite element method (UBFEM), the lowerbound finite element method (LBFEM), and the semianalytical upper-bound rigid-block (UBRB) method. For cohesive-frictional soils, Yamamoto et al. [40,41] determined square tunnel stability with an overload on the ground surface using UBFEM and LBFEM; in these studies, they determined the ultimate overloads and presented an empirical equation for approximating the ultimate overload. At present, little information is available in the literature for investigating rectangular tunnel stability. Based on the UBRB, Yang and Yang [42] performed a limited study on rectangular tunnels in cohesive-frictional soils using a series of rigid blocks, through which they derived an expression for the ultimate support pressure along the tunnel. Later, Abbo et al. [43] and Wilson et al. [44] determined the influence of the tunnel span on the wide rectangular tunnel stability under undrained conditions and proposed an empirical expression for computing stability numbers using both the UBFEM and the LBFEM.
Although several works have been reported on rectangular tunnel stability, stability numbers have mainly been obtained for cases with fully cohesive soils. Moreover, little information has been presented about the effect of overloads on the failure patterns of an unlined rectangular tunnel.
is paper aims to determine an unlined rectangular tunnel stability with an overload on the ground surface for cohesive-frictional soils using a modified version of the UBFEM with rigid triangular elements considering translational freedom. e optimal upper-bound solutions can be obtained by adaptively adjusting the locations and directions of velocity discontinuities during optimization. A set of the dimensionless parameters covering the shallow cover the span ratios of the tunnels and depth ratios of the tunnels and the Mohr-Coulomb material parameters are comprehensively investigated. e stability of an unlined rectangular tunnel is investigated through a dimensionless overload σ s /c, and failure patterns are presented in the form of slip lines that are similar to those assumed in the work of Yang and Yang [42]. Moreover, a new failure pattern that has not been presented in the literature is proposed in this study through artificial amendments of the mesh patterns. Compared with the rigid-block mechanisms from the UBRB in the work of Yamamoto et al. [40], these computed failure mechanisms are better suited for rectangular tunnel stability assessment, especially for high frictional angles. Figure 1 shows the cross section of the rectangular tunnel. It is assumed that the rectangular tunnel length is sufficiently long relative to the tunnel dimension, so that the idealized problem of a shallow rectangular tunnel with an overload on the ground surface is simplified as a plane strain model. e soil mass, with uniform internal friction angle (φ), unit weight (c), and cohesion (c), is established as a Mohr-Coulomb material. e rectangular tunnel (under a cover depth H) has a span B and a height D. In this paper, it is assumed that D � 6 m, and the values of B and H can be obtained through the defined ratios of B/D and H/D. No support force is applied along the tunnel geometry. e overload is presented by applying a continuous load (σ s ) at the ground surface. e horizontal velocity is assumed to be zero to model a rough interface condition between the soil and the load. e collapse is initiated through the action of overload. In this paper, the collapse failure is driven by the surcharge loading, and the computed ultimate surcharge loading is the maximum external load, which can be borne by the tunnel. When the value of the surcharge loading increases to the ultimate value, the tunnel is in a state of critical instability. In addition, in order to facilitate the analysis, dimensionless treatment is carried out. e unlined rectangular tunnel stability can be described using a dimensionless overload (σ s /c), which is a function of the dimensionless unit weight cD/c, φ, dimensionless span ratio B/D, and dimensionless cover depth ratio H/D, as shown in the following expression:

UBFEM-RTME
e UBFEM with rigid triangular elements [37,45,46], considered as an extension of the UBRB, can abandon the complex and long derivation procedures of velocity relationships and geometrical relationships in establishing programming models. For this method, the locations and directions of the velocity discontinuities are two main reasons influencing the computed results. To improve the accuracy, Milani and Lourenco [47] and Hambleton and Sloan [48] optimized the layout of the velocity discontinuities based on multiple successive perturbations, and the second-order cone programming was used to solve each perturbation step. Later, through the direct establishment of a nonlinear programming model, Yang et al. [49] proposed a version of the UBFEM (UBFEM-RTME) whose model is discretized with rigid triangular translator moving elements. Compared to the UBRB, the UBFEM-RTME does not require the construction of an appropriate failure pattern to obtain the solution; when using the UBFEM-RTME, the complex and long derivation procedures of velocity relationships and geometrical relationships can be abandoned owing to the introduction of finite element technology. Moreover, the failure pattern obtained from the UBFEM-RTME is more reasonable than that obtained from the UBRB. ree-node rigid triangular elements with velocity discontinuities are applied to discretize the computational domain where no plastic deformation occurs within the elements. e node coordinates (x i , y i ) and element velocities (u i , v i ) are defined as decision variables to be determined during optimization, and each node corresponds to a particular element. A kinematically admissible velocity field, which is compatible with the velocity boundary conditions along velocity discontinuities and at the boundaries, can be obtained during the solution procedure. Figure 1(b) shows typical upper-bound rigid meshes for a model with B/D � 0.75 and H/D � 1.5. Due to the symmetries of the tunnel, the boundary, and the load, only the right part of the whole region is presented. Figure 1(b) shows that the model is discretized with 461 rigid elements and 669 velocity 2 Advances in Civil Engineering discontinuities. Mesh refinement is performed to ensure that the mesh density gradually increases toward the tunnel periphery. e horizontal velocity component is set to zero along the boundaries P 1 P 2 and P 4 P 6 (i.e., symmetric boundaries), and both the vertical and horizontal velocities are set to zero along the boundaries P 6 P 7 and P 7 P 8 . e chosen domain should be such that further increases in the values L 1 and L 2 (defined in Figure 1) do not influence the results. Note that the mesh presented in Figure 1(b) is used only to acquire an initial solution. To obtain a better upper-bound solution, the element form is optimized through small artificial corrections in the optimization process according to the computed results.
Although the discrete method in this study is similar to that in the work of Yang et al. [49], due to the differences in the external load types and tunnel profiles, a new nonlinear programming model and a new upper-bound rigid mesh type needed to be established in this paper. e ultimate overload is acquired by equating the power dissipated along the velocity discontinuities with the power exerted by the overload and soil weight. e nonlinear programming model is expressed in the following form: Subject to

Advances in Civil Engineering
where P d,i � c · ξ ″ is the power dissipation along the i th velocity discontinuity, P e,i � −A i · c · v i represents the gravity power exerted by the i th element, ξ i ′ and ξ i ″ are the auxiliary parameters that can be expressed with nodal coordinates (x i , y i ) and element velocities (u i , v i ), and A i represents the area of the i th element.
Equation (2a) represents the associate flow rule along the velocity discontinuity, in which n d defines the sum of the velocity discontinuity. Equation (2b) represents the geometry constraints of each element, in which n e defines the sum of elements. Equation (2c) is the displacement boundary condition along the ground surface, wherein the external power expended by the overload becomes the ultimate overload (σ s ); in this equation, n s defines the sum of elements along P 1 P 8 , l i represents the length of the i th element along P 1 P 8 , and v i defines the vertical velocity of the i th element. Equations (2d)-(2k) define the boundary constraints along P 1 P 8 , P 8 P 7 , P 7 P 6 P 6 P 4 , P 2 P 1 , P 2 P 3 , P 4 P 5 , and P 3 P 5 , respectively. In these cases, the nodal coordinates (x i , y i ) and element velocities (u i , v i ) define the nodes and elements at the geometric boundaries, respectively, and parameters (n g1 , n g2 , n g3 , n g4 , n g5 , n g6 , n g7 , n g8 ) and (n v2 , n v3 , n v4 , n v5 ) are the sum of corresponding nodes and elements, respectively. e detailed explanation of the variables mentioned above is identical to those in the work of Yang et al. [8].

Verification.
When the critical overload (σ s ) is acquired through the UBFEM-RTME, the dimensionless overload is conveniently computed using the expression σ s /c. As shown in Table 1, the square tunnel (B/D � 1) stability is determined in terms of the stability tables. A positive value of σ s /c indicates that tunnel failure occurs when a compressive normal stress acts at the ground surface. A negative value indicates that the tunnel has poor stability and that a tensile normal stress needs to be applied at the ground surface to achieve tunnel stability. For comparative analysis, (i) the results from the UBRB and (ii) the average of the lowerbound results and upper-bound results from the UBFEM and the LBFEM for a square tunnel in the work of Yamamoto et al. [40] are listed in Table 1. Because the UBFEM is applied to obtain the ultimate load at the critical state, a smaller value of the overload corresponds to a higher precision solution. Table 1 shows that the UBFEM-RTME shows a relatively good ability to evaluate the stability numbers compared with the UBRB, especially for cases with a higher φ and larger H/D (when φ and H/D are large, there is severe localized plastic flow deformation, and the UBRB has difficulty searching for the best slip line pattern with small numbers of rigid blocks). Although the computed upper-bound results are slightly greater than the average of the upper-bound results and lower-bound results of limit analysis in the work of Yamamoto et al. [40], the errors between these results are generally less than 2%.
When the distance between the centers of dual tunnels is over a critical value, the failure mechanisms for dual tunnels become similar to those for a single tunnel. In this paper, with a critical center-to-center space, (i) the averages of the lower-bound results and upper-bound results for dual square tunnels from the UBFEM and LBFEM [41] are listed in Table 1. In addition, the upper-bound solutions of dual circular tunnels from the UBFEM-RTME [50] are also included. e present dimensionless overloads have the same variation trend as that in the work of Yamamoto et al. [41]. As the interface condition between the soil and the loadings is smooth in the work of Yamamoto et al. [41], the average value is slightly smaller than that of the present analysis. Table 1 also shows that the dimensionless overload of the square tunnel is smaller than that of the circular tunnel, especially for cases with large φ and large cD/c. e dimensionless overload for a square tunnel decreases approximately in the ranges of (i) 13 , and (iii) 15-44% for φ � 10°a nd 41-55% for φ � 30°(H/D � 5) compared with the cases with a circular geometry. Given the same tunnel height, the square tunnel provides a lower stability than the circular tunnel, and this conclusion has been confirmed through quantitative calculations, as shown in Table 1.

Variation in the Stability Numbers for Rectangular
Tunnels. As shown in Table 2, a series of dimensionless overloads σ s /c of a rectangular tunnel for (i) cD/c varying from 0 to 2.5, (ii) B/D varying from 0.5 to 1.5, (iii) φ varying from 10°to 35°, and (iv) H/D varying from 1 to 5 were Table 1: Comparisons between the computed results and those reported in the literature.

H/D φ(°)
Yamamoto et al. [40] Yamamoto et al. [ Note. e results of Yang et al. [23] are obtained with a circular geometry, and the interface condition between the soil and the load is smooth in the work of Yamamoto et al. [41]. Advances in Civil Engineering obtained using the UBFEM-RTME. To more explicitly describe the rectangular tunnel stability, σ s /c was also presented in the forms of upper-bound solution figures for some conditions, as shown in Figures 2 and 3. As shown in Figure 6 Advances in Civil Engineering increasing the tunnel depth for soil with large φ (φ ≥30°), especially for small cD/c. However, for soil with small φ (φ ≤10°), the improvement in the tunnel stability is limited, and the stability number is found to decrease even for cD/c ≥1.5. ese figures also indicate that the value of B/D has a significant influence on σ s /c. σ s /c decreases in the range of approximately (i) 60-148% for H/D � 1 and 31-1227% for H/D � 5 (φ �10°) and (ii) 67-85% for H/D � 1 and 65-88% for H/D � 5 (φ � 30°) as B/D increases from 0.5 to 1.5. Figure 4 shows the optimal failure mechanisms of square tunnels using the UBFEM-RTME (solid line). ese failure modes are acquired by removing invalid velocity discontinuities, where the velocities of adjacent rigid elements are zero or there is no relative motion between adjacent rigid elements. e failure mode consists of a whole shear-gliding wedge block adjacent to the soil surface and a shear-gliding wedge-shaped field comprising a series of slip lines. e major slip lines of the sheargliding wedge-shaped field start at the top and bottom corners of the rectangular tunnel wall and evolve toward the soil surface. It is noticeable that the slip lines intersect around the middle part of the tunnel, and they do not extend to the tunnel wall or the roof. ese mechanisms agree well with the power dissipation presented in the work of Yamamoto et al. [40].

Failure Mechanisms of Rectangular Tunnels.
For comparison purposes, the failure mechanism of a circular tunnel with a critical center-to-center space (dotted line) in the work of Yang et al. [23] for φ � 10°, cD/c � 1, and H/D � 3 is also included in Figure 4(a). Note that although this failure mechanism also consists of a shear-gliding wedge-shaped field and a whole shear-gliding wedge zone, the difference in the failure mechanisms between the square tunnels and circular tunnels is easy to identify. e main differences are (i) the failure pattern around the tunnel, (ii) the slip line pattern within the wedge-shaped zone, and (iii) the starting positions of the slip lines between a rectangular tunnel and a circular tunnel. Figures 4(b) and 4(c) also show the rigid-block mechanisms of square tunnels (dotted line) in the work of Yamamoto et al. [40]. As shown in Figures 4(b) and 4(c), due to the introduction of finite element technology, the final  Advances in Civil Engineering failure patterns from the present analysis are more refined than those from the UBRB, and the computed upper-bound solutions are smaller (higher precision solutions). e high errors in the stability numbers from the UBRB, as shown in Table 1, are explained by the fact that simple rigid block failure modes are less accurate for these cases. Figure 5 illustrates the optimal failure modes of rectangular tunnels for cD/c � 1, H/D � 3, and φ � 10°with various B/D. As the tunnel span increases, the shear-gliding wedge-shaped field extends upward, and the size of the whole shear-gliding wedge block decreases accordingly. Although the pattern of the slip lines shows no significant change, the slope of the major slip line starting at the bottom corner of tunnel wall increases as the tunnel span increases, and the stability number decreases obviously when B/D increases. e maximum horizontal position of the failure mode (i) beneath the surface increases from 1.64D to 1.69D and (ii) that at the ground surface increases from 1.30D to 1.41D when B/D increases from 0.5 to 1.5. Figures 6(a)-6(c) illustrate the failure mechanisms of rectangular tunnels for various cD/c, φ, and H/D, respectively. Only the contours of the failure modes are included in these figures. Figure 6(a) indicates that the pattern of the failure mechanism shows small changes with respect to cD/c, and the slip lines around the tunnel cannot be recognized on the scale of the plot. As cD/c decreases, the maximum horizontal portion of the failure mechanism increases, and the influence of a failure from the right base of the wall is more noticeable. Figure 6(b) shows that, owing to the assumption of the associated flow rule along the velocity discontinuities, the whole shear-gliding wedge block extends upward to the soil surface with an angle of φ between the slip line and the vertical line, and the failure zone at the ground surface decreases from 1.47D to 1.01D with increasing φ. Moreover, when φ increases, the shear-gliding wedgeshaped field extends upward, and the whole shear-gliding wedge zone decreases. Figure 6(c) indicates that the sizes of the shear-gliding wedge-shaped field and the whole sheargliding wedge increase obviously in both the horizontal and vertical directions when H/D increases from 2 to 4. Compared with shallow tunnels, the effect of the failure from the right base of the tunnel is more remarkable for H/D � 4.
Computations show that, for some cases, the errors between the computed results and those in the work of

Advances in Civil Engineering
Yamamoto et al. [40] are relatively high with the failure mechanism mentioned above (mechanism 1). us, another failure pattern (mechanism 2) is presented during the process of artificial amendments of the mesh patterns. Figure 7 shows the failure patterns of square tunnels (B/D � 1) with a large tunnel depth (H/D � 5) using both mechanism 1 and mechanism 2. As shown in Figure 7,   the computed stability number increases by more than 20% over those obtained using mechanism 2 for some cases. In this paper, the computed results are obtained with mechanism 2. Note that although mechanism 2 is used to determine the rectangular tunnel stability, the failure pattern can degenerate into mechanism 1 by removing invalid velocity discontinuities. In this paper, two kinds of failure mechanism are both used to compute the upper-bound solutions for the cases listed in Table 2, and the better result (the smaller one) is selected as the optimal upper-bound solution. For larger friction angles, deeper tunnels, and smaller values of cD/c, these failure patterns are more complex than those for the other cases. As the solutions are obtained using a nonlinear programming model, mechanism 1 and mechanism 2 may be local optimal solutions, and other solutions may exist for some cases. However, the computed results in this paper are still upper-bound solutions.

Conclusions
e stability of an unlined rectangular tunnel in cohesivefrictional soil affected by overload was assessed using a version of the UBFEM with rigid triangular translator moving elements. e results are presented in the form of dimensionless overloads and failure modes. e dimensionless overload increases with increasing internal friction angle φ, whereas the dimensionless overload decreases with increasing dimensionless span B/D and dimensionless unit weight cD/c. Note that σ s /c decreases in the range of approximately (i) 60-148% for H/D � 1 and 31-1227% for H/D � 5 (φ �10°) and (ii) 67-85% for H/D � 1 and 65-88% for H/D � 5 (φ � 30°) as B/D increases from 0.5 to 1.5. e dimensionless overload for the square tunnel decreases by (i) 59% with H/D � 1 and (ii) 55% with H/D � 5 compared with the cases with a circular geometry. e failure modes are mostly composed of a shear-gliding wedge-shaped field around the tunnel and a whole sheargliding wedge block adjacent to the surface. For shallow tunnels, the major slip lines start at the top and bottom corners of the wall, and these slip lines intersect around the middle part of the tunnel. Note that although the blocks around the roof and wall do not substantially change, the influence of a failure from the zone at the floor of the tunnel is more conspicuous for larger H/D, larger φ, and smaller cD/c. ese failure patterns are more complex than those for other cases, and the construction quality of the tunnel floor should be better controlled to avoid tunnel heave.
In this paper, the Mohr-Coulomb criterion is used to characterize the strength of soils and the associated flow rule is assumed. For the cases with Mohr-Coulomb materials under drained conditions, nonassociated flow rule can be introduced to analyze the tunnel stability to avoid excess volume expansion. However, for soils with various cohesion (c) and friction angle (φ), the reduction degrees of the shear strengths are also different. To be close to the practical situation, more laboratory tests need to be done. Further studies combining with Mohr-Coulomb materials considering nonassociated flow rule will be studied in future work.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
ere are no conflicts of interest regarding the publication of this paper.