#### Abstract

This paper investigates the stability of a rectangular tunnel face affected by surcharge loading in soil with a soft upper layer and hard lower layer using upper-bound finite element methods with a plastic-dissipation-based mesh adaptive strategy (UBFEM-PDMA). Seven different positions for the soil interface are selected to study this problem. The upper bounds on the ultimate surcharge loads *σ*_{s} are presented in terms of dimensionless stability charts. The *σ*_{s} increases with tunnel depth, and it increases when the position of the soil interface moves up along the tunnel face. The failure mechanism primarily involves a wedge-shaped zone around the tunnel face and two slip lines originating from the top and bottom of the tunnel face, and it is mainly influenced by three factors, i.e., the position of the soil interface, the soil properties, and the tunnel depth. In contrast to the failure mechanism for uniform soil, multiple slip lines exist in the tunnel face in soil with a soft upper layer and hard lower layer. The results compare reasonably well with those in the literature and those from the numerical method.

#### 1. Introduction

For cities located on soft ground adjacent to rivers and coastlines in China, city tunnels are mostly constructed at shallow depths with large fluctuations in the soil interface. Compared with uniform formations, the collapse mechanisms of the tunnel face in these formations are complex, and instability in the tunnel face can easily occur during the construction process. Research on the effects of surcharge loading on the stability and failure mechanism of a tunnel face in heterogeneous formations is important, and these issues must be addressed in tunnel engineering.

The limit equilibrium method [1–3], model-based experiments [4–7], reliability method [8, 9], and numerical simulations [10–13] are the most common methods used to analyze the stability of the underground engineering. Recently, limit analysis, based on the plastic theorems proposed by Drucker et al. [14, 15], has also been widely used to study tunnel face stability. Compared with the common methods, this method has a strict theory, and it can provide a convenient way of obtaining bounds on parameters describing the stability, as well as the failure mechanisms, without considering a complicated elastic-plastic deformation process. The failure-mechanism-based upper-bound method (such as the rigid-block method) and the upper-bound finite element method are the two main methods used to obtain the solutions of limit analysis. For circular tunnels, the failure-mechanism-based upper-bound method is usually used to find the upper-bound solutions based on an assumed three-dimensional failure mechanism [16–25]. For rectangular tunnels whose spans are sufficiently large relative to their heights, the problem of underground excavation stability can be simplified as a plane strain analysis model, and some research results are obtained using these two methods [26–28]. Based on an assumed two-dimensional failure mechanism, Huang et al. [27] obtained upper-bound stability solutions of a plane strain tunnel heading in nonhomogeneous clay and proposed an improved simple failure mechanism. Sloan and Assadi [28] and Augarde et al. [26] applied the finite element limit analysis method to study the undrained stability of a shallow heading under conditions of plane strain surcharge loading, and rigorous upper- and lower-bound solutions were obtained to generally bracket the true solution. Later, Yang et al. [29] extended the research to cohesive frictional soils, and the failure mechanisms are clearly discussed using the upper-bound finite element method with rigid translatory moving elements (UBFEM-RTME) presented by Yang et al. [30].

Although significant information on the stability of a rectangular tunnel face subjected to surcharge loading has been reported, the solutions are obtained mainly for uniform soils (undrained clay or cohesive frictional soil) or soils whose strength increases linearly with depth. Discussions on tunnel face stability in soils with large differences in soil properties are seldom available. Moreover, limited research has been undertaken to study the effect of surcharge loads on the characteristics of failure mechanisms for plane strain heading. Therefore, this paper aims to analyze the rectangular tunnel face stability in soil with a soft upper layer and hard lower layer using an upper-bound finite element method in combination with a plastic-dissipation-based mesh adaptive strategy (UBFEM-PDMA). Upper-bound solutions with high accuracy and the failure mechanism exhibiting the evolutionary characteristics of slip lines can be obtained by automatically refining the failure regions with high strain rates. On the one hand, the upper-bound solutions for uniform formation in Augarde et al. [26] and Yang et al. [29] are improved. On the other hand, the effects of different positions of the soil interface, soil properties, and different tunnel depths on critical surcharge loading and the failure mechanism are investigated, and the differences in the results between soil with a soft upper layer and hard lower layer and uniform soil are discussed. The results are compared with those in the literature and those from the numerical method.

#### 2. Problem Description

Figure 1 shows the analytical model of a rectangular tunnel face in soil with a soft upper layer and hard lower layer. It is assumed that the span of the rectangular tunnel is large enough relative to the tunnel height, and thus, the problem can be simplified as a two-dimensional plane strain analysis problem. Figure 1(a) shows that the rectangular tunnel has a height *D* under a depth of cover *C*. The soil mass is established as a Mohr–Coulomb material with a drained internal friction angle (*ϕ*), cohesion (*c*), and unit weight (*γ*). Continuous loading (*σ*_{s}) is applied at the ground surface to act as the surcharge loading, and the tunnel face is free to move. As the lining ground interface condition is complex, based on the operations in Augarde et al. [26] and Yang et al. [29], an infinitely rigid lining is placed along the heading to simplify this problem. This assumption seems reasonable because this paper aims to analyze the tunnel stability rather than determine the deformation. The dotted lines 1∼7 indicate that the soil interface is located at (i) the bottom of the model (case 1), (ii) the bottom of the tunnel face (case 2), (iii) the bottom half (three-quarters) of the tunnel face (case 3), (iv) the center (two quarters) of the tunnel face (case 4), (v) the top half (one-quarter) of the tunnel face (case 5), (vi) the top of the tunnel face (case 6), and (vii) the top of the model (case 7). The soil above the soil interface is assumed to be soft soil, and the soil below the soil interface is assumed to be hard soil. Therefore, case 1 corresponds to uniform formation of the soft soil, and case 7 corresponds to uniform formation of the hard soil. To simplify the analysis, typical soil parameters, rather than specific soil parameters, are selected for the soft soil and hard soil. The strength parameters of the soil interface are the same as those of the soft soil. The failure of the tunnel face is assumed to be driven by the action of gravity and surcharge loading, and a critical surcharge loading *σ*_{s} can be used to describe the stability of the tunnel face. The value of *σ*_{s} is influenced by the position of the soil interface, the soil strength, and the tunnel depth.

**(a)**

**(b)**

#### 3. UBFEM-PDMA

##### 3.1. Establishment of a Linear Programming Model

Three-node triangular plastic deformation elements with a constant strain rate are applied to discretize the model. These elements possess the characteristic of translation without rotating freedom, and the velocity variables changed in a linear fashion within the element. Because the incompressibility of three-node triangular plastic deformation elements is overestimated, discontinuous velocity is permitted to occur among elements to improve the accuracy of the results. Figure 1(b) shows the typical element pattern (511 elements) of a rectangular tunnel face with *C/D* = 1 and the soil interface located at the center of the tunnel face (case 4). Similar mesh patterns are applied to other cases. The length of the heading *M* is assumed to be 1.5*D*. The width of the model *M* + *L*_{2} is 5*D*, and the height of the model *C* + *D* + *L*_{1} is 3.5*D*. The velocity constraints (*u* *=* 0*, v* *=* 0) are imposed along the boundaries DE, AF, AB, BC, EG, and FH. No velocity constraints need to be set on the tunnel face. The constraint created by surcharge loading is set as .

Considering the difficulties in obtaining the initial values of decision variables and solving a nonlinear programming model efficiently, a regular *p*-side polygon, which is external to the yield surface, is employed to linearize the Mohr–Coulomb yield criterion. This handing method can ensure that the result always remains a rigorous upper bound on the exact solution, which is also explained by Sloan and Kleeman [31]. The objective function can be obtained by minimizing the internal power dissipation (dissipated within elements and along velocity discontinuities) minus the rate of work expended by external forces (gravity). The linear programming model takes the following form:

Subject to

The dissipated power of the velocity discontinuity can be obtained through the expression . defines the tangential velocity jump for the velocity discontinuity, and is the length of the velocity discontinuity. The external power exerted by the gravity of the element can be formulated as . *A*_{i} is the area of the element, and is the vertical velocity of the element. The dissipated power of the element can be represented as , and is a nonnegative plastic multiplier rate.

The constraints listed in equations (2a), (2c), (2d), and (2i) are generated from the satisfaction of the plastic flow in continuum and along velocity discontinuities, respectively. are the normal and shear strain rates, and are the normal and shear stresses. is the function for the side of the *p*-side yield polygon, and is the amount of elements. are the normal velocity jumps of the endpoints of the velocity discontinuity, and are the tangential velocity jumps of the endpoints of the velocity discontinuity. represent the nonnegative auxiliary variables of the endpoints, which are introduced to eliminate the nonlinear expressions. are the horizontal velocities of the nodes, and are the vertical velocities of the nodes. is internal friction angle of the soft soil and is internal friction angle of the hard soil. defines the angle of the discontinuity relative to the *x*-axis, and is the account of velocity discontinuities. Equations (2j)–(2p) define the boundary constraint of velocity imposed along DC, AB, BC, DE, FA, EG, and FH. (*ng*_{1}, *ng*_{2}, *ng*_{3}, *ng*_{4}, *ng*_{5}, *ng*_{6}) are the account of corresponding nodes.

##### 3.2. Plastic-Dissipation-Based Mesh Adaptive Strategy (PDMA)

To improve the accuracy of the upper bound solutions, a self-adaptive refining algorithm, which is similar to those reported by Martin [32] and Zhang et al. [33] is introduced. In contrast to Martin [32] and Zhang et al. [33], in this paper, three-node triangular elements, rather than six-node triangular elements, are assumed to be refined. An indicator based on an evaluation of plastic dissipation of elements was presented to implement the self-adaptive refining algorithm, and the plastic collapse domain (with high strain rates) can be captured automatically in the optimizing process. The present mesh adaptive strategy can help achieve a goal that each element does not make an excessive contribution to the internal power dissipation, which means that the elements located at failure regions can be refined. Due to local refinement, the plastic failure problems can be solved using a strain-rate field with the lowest computational cost. The overall numerical procedure of the plastic-dissipation-based mesh adaptive strategy is as follows.

The initial upper-bound solution of the analysis model is first computed using a MOSEK solver based on the initial mesh pattern, and the indicator (internal power dissipation) *P*_{p} for all the elements is obtained through the following expression:where defines the nonnegative plastic multiplier rate of the element and *A*_{i} defines the area of the element.

Then, the elements are ordered from largest to smallest based on indicators of all elements, and the quantity and corresponding number of refinement elements can be identified using the following expression of *η*:where defines the account of marked elements and defines the account of elements of the model.

Through this operation, the elements with high strain rates can be selected automatically based on computed requirements. A suitable value of *η* is obtained through trial calculation, and 0.2∼0.6 are often preferred. Note that if the value of *η* is too small, a large number of loop iterations are needed, and some plastic zones cannot be refined when the elements in local regions make an excessive contribution to the internal power dissipation. If the value of *η* is too large, those elements with low power dissipation may be refined, and the solving efficiency decreases. The algorithm of vertex bisection is applied for mesh refinement. For each triangular element that needs to be refined, two new triangular elements are generated by connecting the midpoint of the longest edge to the peak or vertex that is opposite of the longest edge.

For each adaptive step, the information regarding nodes, edges, and elements is updated with that of the new refined mesh, and the newest upper-bound solution can be obtained based on the refined mesh pattern. Based on the upper-bound theorem, the value of the upper-bound solution decreases when the iteration increases. When the gap between the result of critical surcharge loading in the current iteration and that in the last iteration is sufficiently small, the iteration is suspended, and the latest result can be taken as the final optimal critical surcharge loading. The gap of the upper-bound solution can be expressed aswhere is the upper-bound solution after the refinement, is the upper-bound solution after the refinement, is the upper-bound solution resulting from the initial mesh pattern, and is the gap of the upper-bound solution after the refinement.

#### 4. Results and Comparisons

##### 4.1. Comparison of the Results with Those in the Literature for Cases with Uniform Soil

The computations were first carried out to obtain the dimensionless critical surcharge loading *σ*_{s}*/c* for cases with uniform soil by varying (i) *ϕ* between 0° and 25°, (ii) *γ*_{D}/*c* between 0 and 3, and (iii) *C/D* between 1 and 4, as shown in Figure 2. In this paper, the proper value of *η* changes from 0.25 to 0.35. During computation, it was determined that the value of *p* can be set as 48 to maintain the economy and accuracy of the calculations. Note that a positive value of *σ*_{s}*/c* shows that a compressive normal stress should be applied at the ground surface to initiate tunnel failure, whereas a negative value implies that, theoretically, a tensile normal stress has to be applied to achieve tunnel stability. For comparison purposes, (i) the upper-bound solutions by Yang et al. [29] obtained by using the UBFEM-RTME, (ii) the upper-bound solutions by Augarde et al. [26] obtained by using the upper-bound finite element method (UBFEM), (iii) and the lower-bound solutions of Augarde et al. [26] obtained by using the lower-bound finite element method (LBFEM) are also included in Figure 2. Because the upper-bound method is used to obtain the ultimate loading at the critical state, in this paper, a smaller surcharge loading results in a more accurate solution. As seen in Figure 2(a), the computed *σ*_{s}*/c* is marginally lower than the upper-bound solutions by Augarde et al. [26] and slightly greater than the lower-bound solutions by Augarde et al. [26]. In addition, Figures 2(b)–2(d) show that (i) for cases with small *C/D*, *ϕ*, and *γ*_{D}*/c*, the upper-bound solutions from the present method are in good agreement with those from Yang et al. [29], and (ii) for other cases, the computed results are found lower than those from Yang et al. [29] (when the value of *ϕ* is large, the plastic deformation of the local zone is intensive, and it is difficult to search the best slip line pattern with small numbers of elements even though node coordinates of elements can be adjusted during the solving process; thus, the slip pattern of the plastic collapse zone with high strain rates from Yang et al. [29] is less reasonable than that from the presented method). These comparisons indirectly indicate that the obtained solutions can obtain better results than those in the literature.

**(a)**

**(b)**

**(c)**

**(d)**

##### 4.2. Variation in the Surcharge Loading

Notably, because many types of soils exist, this paper cannot discuss all the scenarios. The typical soil parameters for three soils (two kinds of soft soils and one kind of hard soil) in cities located on soft ground adjacent to rivers and coastlines in China and seven different positions of the soil interface are selected to analyze this problem. Much can still be learned from the discussion of variation rules of surcharge loading and evolution characteristics of failure mechanisms for these situations. The main parameters of soft 1 are *γ* = 19 kN/m^{3}, *ϕ* = 12°, and *c* = 40 kPa, and the main parameters of soft 2 are *γ* = 19 kN/m^{3}, *ϕ* = 25°, and *c* = 20 kPa. The main parameters of the hard soil are as follows: *γ* = 23 kN/m^{3}, *ϕ* = 32°, and *c* = 80 kPa.

The computed surcharge loads *σ*_{s} in seven cases for two kinds of soft soils are all listed in Tables 1 and 2. As seen in the tables, when the soil interface moves from the bottom of the model to the bottom of the tunnel face, the position of the soil interface has less impact on *σ*_{s}. When the position of the soil interface moves up along the tunnel face, the value of *σ*_{s} increases significantly, especially for the cases with a large tunnel depth. The *σ*_{s} increases by approximately (a) (i) 612% (*C/D* = 1), (ii) 891% (*C/D* = 2), (iii) 1384% (*C/D* = 3), and (iv) 3071% (*C/D* = 4) for soft 1 and (b) (i)1956% (*C/D* = 1), (ii) 2318% (*C/D* = 2), (iii) 2386% (*C/D* = 3), and (iv) 2313% (*C/D* = 4) for soft 2 when the soil interface moves from the bottom to the top of the tunnel face. When the soil interface is located at the top of the model, *σ*_{s} is much larger than that for the other cases. Based on the above discussion, when the soil interface is below the bottom of the tunnel face, the value of *σ*_{s} is close to that for uniform soft soil (case 1). When the soil interface changes along the tunnel face, the difference in *σ*_{s} between soil with a soft upper layer and hard lower layer and uniform soil is large. Thus, in tunnel engineering, if the difference in soil properties is ignored by engineers, the difference between the evaluation result and the actual value is large.

##### 4.3. Discussion of the Failure Mechanism

Figure 3 shows the failure mechanisms of a rectangular tunnel face with uniform soil for soft 1, soft 2, and hard soil. As seen in Figure 3(a), the collapse mechanism for soft 1 primarily involves three slip zones: (i) zone 1: the slip line originates from the top of the tunnel, and it curves toward the ground surface; (ii) zone 2: two slip lines originate from the top and the bottom of the tunnel, and a wedge-shaped zone exists around the top of the tunnel; and (iii) zone 3: the slip line originates from the bottom of the tunnel, and it curves toward the ground surface. As shown in Figure 3(b), for soft 2, although the failure mechanism consists of three slip zones, plastic deformation occurs mainly in zone 2. When the soil mass is hard soil, as shown in Figure 3(c), the collapse mechanism primarily involves zone 2, and the failure zone is mainly located around the tunnel face. These failure mechanisms show good agreement with those in [29] using the UBFEM-RTME.

**(a)**

**(b)**

**(c)**

For the cases of soil with a soft upper layer and hard lower layer, the failure mechanisms of the tunnel face for soft 1 and soft 2 are displayed in Figures 4 and 5. A solid arrow is used to mark the position of the soil interface. Comparing Figure 3 with Figures 4 and 5, it can be found that for both soft 1 and soft 2, no significant change occurs in the failure mechanism when the soil interface moves from the bottom of the model to the bottom of the tunnel face. This finding also explains why the value of *σ*_{s} for case 2 is close to that for uniform soft soil (case 1).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

Figures 4(b)–4(d) show that the failure mechanism consists of zone 1, zone 2, and zone 3; however, the failure mechanism differs from those for uniform deformation. (i) For zone 2, the two slip lines of the wedge-shaped zone originate from the top of the tunnel face and the soil interface, and the area of this zone decreases. (ii) For zone 3, the angles of the slip line relative to the *x*-axis increases when the soil interface moves up along the tunnel face. When the soil interface is located at the top of the tunnel face, as shown in Figure 4(e), the failure mechanism mainly consists of zone 2, and the two slip lines for the wedge-shaped zone originate from the top and the bottom of the tunnel face.

Figures 5(b)–5(c) show that the failure mechanism mainly consists of zone 2 and zone 3. The area of the wedge-shaped zone in zone 2 decreases when the soil interface moves up along the tunnel face. Differing from the failure mechanism for soft 1, the slip line in zone 3 for soft 2 originates from the bottom of the tunnel face, and it curves toward the soil interface. When the soil interface continues to move up, as shown in Figures 5(d)–5(e), the failure mechanism mainly consists of zone 2, and the slip lines of the wedge-shaped zone originate from the top and bottom of the tunnel face.

Figure 6 shows the failure mechanism for soft 1 with a large tunnel depth. Because the failure mechanism for soft 2 mainly consists of zone 2 (no significant changes occur in the failure mechanism for various C values), this paper does not cover those failure mechanisms again. Comparing Figure 4 with Figure 6, it can be found that when C increases, (a) for case 2, (i) the area of the failure zone increases and (ii) the curvature of the slip line in zone 3 decreases; (b) for case 4, the plastic strain rates in zone 2 increase and the plastic strain rates in zone 1 and zone 3 decrease; and (c) for case 6, zone 2 becomes the main failure zone.

**(a)**

**(b)**

**(c)**

##### 4.4. Comparison of the Results with Those from the Numerical Method

To verify the correctness of the presented method in analyzing the tunnel face stability in soil with a soft upper layer and hard lower layer, the ABAQUS software is also applied to compute the critical surcharge loading and failure mechanism for soft 1 with case 4 and *C/D* = 1 by using 2D model. The length and width of the model are all set to 5*D*. The boundary conditions are consistent with those for the UBFEM-PDMA. CPE4 (four-node rectangular elements) is applied to discretize the computational domain. The Mohr–Coulomb model is adopted to simulate the soil behavior, and the parameters for soil are listed in Table 3.

The overall numerical procedure for analyzing tunnel face stability is described as follows:(1)Establish an analysis model with CPE4, and obtain the initial geostress field(2)Remove the elements located in the excavated region, and apply a small vertical uniform surcharge loading at the ground surface(3)Solve the model and monitor the displacement of elements located at the tunnel face(4)Increase the surcharge loading and repeat steps1 to 3(5)If the element displacement increases significantly when the surcharge loading increases slightly, then the load is considered to be the critical surcharge loading

Through the calculation, the computed critical surcharge loading is determined to be 362 kPa, and it is greater than that (353.79 kPa) resulting from the UBFEM-PDMA. Figure 7 shows the plastic strain distribution of the numerical simulation method and the refined mesh of the UBFEM-PDMA. As shown in Figure 7, the intensive element zones in Figure 7(b) show good agreement with the dark-shaded areas in Figure 7(a). From the discussion above, the UBFEM-PDMA is concluded to be an effective tool to analyze the tunnel face stability in both soil with a soft upper layer and hard lower layer and uniform soil, and the precision of calculation of the UBFEM-PDMA is better than that of the numerical simulation method. In addition, the solving process of the critical load using the numerical simulation method is tedious, and the accurate critical state is difficult to obtain. However, if the object of the study is the determination of deformation, rather than failure conditions, the numerical simulation method would be appropriate for application.

**(a)**

**(b)**

#### 5. Conclusions

The stability of a plane strain rectangular tunnel face in soil with a soft upper layer and hard lower layer and in uniform soil subjected to surcharge loading was investigated using the UBFEM-PDMA. The results of these analyses have been presented in the form of dimensionless stability charts and failure mechanisms.(1)The critical surcharge loading *σ*_{s} increases when C increases. The value of *σ*_{s} is close to that for uniform soft soil when the soil interface is below the bottom of the tunnel face, and it increases significantly when the position of the soil interface moves up along the tunnel face. The *σ*_{s} increases by approximately 3071% (*C/D* = 4) when the soft soil is soft 1, and it increases by approximately 2386% (*C/D* = 3) when the soft soil is soft 2.(2)When the soil interface is below the bottom of the tunnel face, the collapse mechanism primarily involves a wedge-shaped zone around the tunnel face and two main slip lines originating from the top and bottom of the tunnel face. When the soil interface moves along the tunnel face, the different soil interface positions have a large effect on the characteristics of the wedge-shaped zone, and multiple slip lines exist in the tunnel face, which is different from the failure mechanism for uniform soil. When *C* increases, the area of the failure zone increases and the failure pattern has no significant change.

#### Data Availability

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

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (grant nos. 51808193 and 51708564), the Fundamental Research Funds for the Central Universities (grant nos. 2018B57014, 2018B42414, and 18lgpy31), the China Postdoctoral Science Foundation (grant no. 2018M633223), the Guangzhou Science & Technology Program of China (grant no. 201804010107), and the Jiangsu Province Postdoctoral Research Funding Scheme (grant no. 2018K044B). The authors are grateful for this support.