Research Article  Open Access
UpperBound Finite Element Adaptive Analysis of Plane Strain Heading in Soil with a Soft Upper Layer and Hard Lower Layer
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 upperbound finite element methods with a plasticdissipationbased mesh adaptive strategy (UBFEMPDMA). 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 wedgeshaped 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], modelbased 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 elasticplastic deformation process. The failuremechanismbased upperbound method (such as the rigidblock method) and the upperbound finite element method are the two main methods used to obtain the solutions of limit analysis. For circular tunnels, the failuremechanismbased upperbound method is usually used to find the upperbound solutions based on an assumed threedimensional 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 twodimensional failure mechanism, Huang et al. [27] obtained upperbound 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 lowerbound 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 upperbound finite element method with rigid translatory moving elements (UBFEMRTME) 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 upperbound finite element method in combination with a plasticdissipationbased mesh adaptive strategy (UBFEMPDMA). Upperbound 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 upperbound 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 twodimensional 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 (threequarters) of the tunnel face (case 3), (iv) the center (two quarters) of the tunnel face (case 4), (v) the top half (onequarter) 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. UBFEMPDMA
3.1. Establishment of a Linear Programming Model
Threenode 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 threenode 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.5D. The width of the model M + L_{2} is 5D, and the height of the model C + D + L_{1} is 3.5D. 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 pside 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 pside 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 xaxis, 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. PlasticDissipationBased Mesh Adaptive Strategy (PDMA)
To improve the accuracy of the upper bound solutions, a selfadaptive 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, threenode triangular elements, rather than sixnode triangular elements, are assumed to be refined. An indicator based on an evaluation of plastic dissipation of elements was presented to implement the selfadaptive 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 strainrate field with the lowest computational cost. The overall numerical procedure of the plasticdissipationbased mesh adaptive strategy is as follows.
The initial upperbound 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 upperbound solution can be obtained based on the refined mesh pattern. Based on the upperbound theorem, the value of the upperbound 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 upperbound solution can be expressed aswhere is the upperbound solution after the refinement, is the upperbound solution after the refinement, is the upperbound solution resulting from the initial mesh pattern, and is the gap of the upperbound 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 upperbound solutions by Yang et al. [29] obtained by using the UBFEMRTME, (ii) the upperbound solutions by Augarde et al. [26] obtained by using the upperbound finite element method (UBFEM), (iii) and the lowerbound solutions of Augarde et al. [26] obtained by using the lowerbound finite element method (LBFEM) are also included in Figure 2. Because the upperbound 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 upperbound solutions by Augarde et al. [26] and slightly greater than the lowerbound 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 upperbound 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 wedgeshaped 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 UBFEMRTME.
(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 wedgeshaped 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 xaxis 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 wedgeshaped 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 wedgeshaped 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 wedgeshaped 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 5D. The boundary conditions are consistent with those for the UBFEMPDMA. CPE4 (fournode 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 UBFEMPDMA. Figure 7 shows the plastic strain distribution of the numerical simulation method and the refined mesh of the UBFEMPDMA. As shown in Figure 7, the intensive element zones in Figure 7(b) show good agreement with the darkshaded areas in Figure 7(a). From the discussion above, the UBFEMPDMA 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 UBFEMPDMA 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 UBFEMPDMA. 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 wedgeshaped 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 wedgeshaped 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.
References
 G. Anagnostou and K. Kovári, “The face stability of slurryshielddriven tunnels,” Tunnelling and Underground Space Technology, vol. 9, no. 2, pp. 165–174, 1994. View at: Publisher Site  Google Scholar
 Z. Eisenstein and L. Samarasekera, “Stability of unsupported tunnels in clay,” Canadian Geotechnical Journal, vol. 29, no. 4, pp. 609–613, 1992. View at: Publisher Site  Google Scholar
 S. Jancsecz and W. Steiner, “Face support for a large mixshield in heterogeneous ground conditions,” Tunnelling’ 94, Springer, Boston, MA, USA, 1994. View at: Publisher Site  Google Scholar
 J. H. Atkinson and D. M. Potts, “Stability of a shallow circular tunnel in cohesionless soil,” Géotechnique, vol. 27, no. 2, pp. 203–215, 1977. View at: Publisher Site  Google Scholar
 E. H. Davis, M. J. Gunn, R. J. Mair, and H. N. Seneviratine, “The stability of shallow tunnels and underground openings in cohesive material,” Géotechnique, vol. 30, no. 4, pp. 397–416, 1980. View at: Publisher Site  Google Scholar
 R. J. Mair, R. N. Taylor, and A. Bracegirdle, “Subsurface settlement profiles above tunnels in clays,” Géotechnique, vol. 43, no. 2, pp. 315–320, 1993. View at: Publisher Site  Google Scholar
 A. N. Schofield, “Cambridge geotechnical centrifuge operations,” Géotechnique, vol. 30, no. 3, pp. 227–268, 1980. View at: Publisher Site  Google Scholar
 J. Ji, C. Zhang, Y. Gao, and J. Kodikara, “Effect of 2D spatial variability on slope reliability: a simplified FORM analysis,” Geoscience Frontiers, vol. 9, no. 6, pp. 1631–1638, 2018. View at: Publisher Site  Google Scholar
 Q. Lü, Z.P. Xiao, J. Ji, and J. Zheng, “Reliability based design optimization for a rock tunnel support system with multiple failure modes using response surface method,” Tunnelling and Underground Space Technology, vol. 70, pp. 1–10, 2017. View at: Publisher Site  Google Scholar
 U. Boonchai, Y. Kongkit, and K. Suraparb, “Threedimensional undrained tunnel face stability in clay with a linearly increasing shear strength with depth,” Computers and Geotechnics, vol. 88, pp. 146–151, 2017. View at: Publisher Site  Google Scholar
 S. H. Kim and F. Tonon, “Face stability and required support pressure for TBM driven tunnels with ideal face membrane—drained case,” Tunnelling and Underground Space Technology, vol. 25, no. 5, pp. 526–542, 2010. View at: Publisher Site  Google Scholar
 Y. Li, F. Emeriault, R. Kastner, and Z. X. Zhang, “Stability analysis of large slurry shielddriven tunnel in soft clay,” Tunnelling and Underground Space Technology, vol. 24, no. 4, pp. 472–481, 2009. View at: Publisher Site  Google Scholar
 M. J. Maynar and L. E. Rodríguez, “Discrete numerical model for analysis of earth pressure balance tunnel excavation,” Journal of Geotechnical and Geoenvironmental Engineering, vol. 131, no. 10, pp. 1234–1242, 2005. View at: Publisher Site  Google Scholar
 D. C. Drucker, H. J. Greenberg, and W. Prager, “The safety factor of an elastic plastic body in plane strain,” Journal of Applied MechanicsTransactions of the ASME, vol. 18, no. 4, pp. 371–378, 1951. View at: Google Scholar
 D. C. Drucker, W. Prager, and H. J. Greenberg, “Extended limit design theorems for continuous media,” Quarterly of Applied Mathematics, vol. 9, no. 4, pp. 381–389, 1952. View at: Publisher Site  Google Scholar
 M. Huang, Z. Tang, W. Zhou, and J. Yuan, “Upper bound solutions for face stability of circular tunnels in nonhomogeneous and anisotropic clays,” Computers and Geotechnics, vol. 98, pp. 189–196, 2018. View at: Publisher Site  Google Scholar
 A. Klar and B. Klein, “Energybased volume loss prediction for tunnel face advancement in clays,” Géotechnique, vol. 64, no. 10, pp. 776–786, 2014. View at: Publisher Site  Google Scholar
 E. Leca and L. Dormieux, “Upper and lower bound solutions for the face stability of shallow circular tunnels in frictional material,” Géotechnique, vol. 40, no. 4, pp. 581–606, 1990. View at: Publisher Site  Google Scholar
 I.M. Lee, S.W. Nam, and J.H. Ahn, “Effect of seepage forces on tunnel face stability,” Canadian Geotechnical Journal, vol. 40, no. 2, pp. 342–350, 2003. View at: Publisher Site  Google Scholar
 I.M. Lee, J.S. Lee, and S.W. Nam, “Effect of seepage force on tunnel face stability reinforced with multistep pipe grouting,” Tunnelling and Underground Space Technology, vol. 19, no. 6, pp. 551–565, 2004. View at: Publisher Site  Google Scholar
 G. Mollon, D. Dias, and A.H. Soubra, “Face stability analysis of circular tunnels driven by a pressurized shield,” Journal of Geotechnical and Geoenvironmental engineering, vol. 136, no. 1, pp. 215–229, 2010. View at: Publisher Site  Google Scholar
 G. Mollon, D. Dias, and A.H. Soubra, “Continuous velocity fields for collapse and blowout of a pressurized tunnel face in purely cohesive soil,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 37, no. 13, pp. 2061–2083, 2012. View at: Publisher Site  Google Scholar
 A. H. Soubra, “Kinematical approach to the face stability analysis of shallow circular tunnels,” in Proceedings of the 8th International Symposium on Plasticity, pp. 443–445, British Columbia, Canada, April 2002. View at: Google Scholar
 D. Subrin, D. Branque, N. Berthoz et al., “Kinematic 3D approaches to evaluate TBM face stability: comparison with experimental laboratory observations,” in Proceedings of 2th International Conference on Computational Methods in Tunneling, pp. 801–808, Ruhr University Bochum, Bochum, Germany, September 2009. View at: Google Scholar
 F. Zhang, Y. F. Gao, Y. X. Wu, and N. Zhang, “Upperbound solutions for face stability of circular tunnels in undrained clays,” Géotechnique, vol. 68, no. 1, pp. 76–85, 2018a. View at: Publisher Site  Google Scholar
 C. E. Augarde, A. V. Lyamin, and S. W. Sloan, “Stability of an undrained plane strain heading revisited,” Computers and Geotechnics, vol. 30, no. 5, pp. 419–430, 2003. View at: Publisher Site  Google Scholar
 M. Huang and C. Song, “Upperbound stability analysis of a plane strain heading in nonhomogeneous clay,” Tunnelling and Underground Space Technology, vol. 38, pp. 213–223, 2013. View at: Publisher Site  Google Scholar
 S. W. Sloan and A. Assadi, “Undrained stability of a plane strain heading,” Canadian Geotechnical Journal, vol. 31, no. 3, pp. 443–450, 1994. View at: Publisher Site  Google Scholar
 F. Yang, J. Zhang, L. Zhao, and J. Yang, “Upperbound finite element analysis of stability of tunnel face subjected to surcharge loading in cohesivefrictional soil,” KSCE Journal of Civil Engineering, vol. 20, no. 6, pp. 2270–2279, 2015. View at: Publisher Site  Google Scholar
 F. Yang, J. Zhang, J. Yang, L. Zhao, and X. Zheng, “Stability analysis of unlined elliptical tunnel using finite element upperbound method with rigid translatory moving elements,” Tunnelling and Underground Space Technology, vol. 50, pp. 13–22, 2015. View at: Publisher Site  Google Scholar
 S. W. Sloan and P. W. Kleeman, “Upper bound limit analysis using discontinuous velocity fields,” Computer Methods in Applied Mechanics and Engineering, vol. 127, no. 1–4, pp. 293–314, 1995. View at: Publisher Site  Google Scholar
 C. M. Martin, “Undrained collapse of a shallow planestrain trapdoor,” Géotechnique, vol. 59, no. 10, pp. 855–863, 2009. View at: Publisher Site  Google Scholar
 J. Zhang, T. Feng, J. Yang, F. Yang, and Y. Gao, “Upperbound stability analysis of dual unlined horseshoeshaped tunnels subjected to gravity,” Computers and Geotechnics, vol. 97, pp. 103–110, 2018. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2019 Jian Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.