Upper Bound Solution for the Face Stability of Shield Tunnel below the Water Table

By FE simulation with Mohr-Coulomb perfect elastoplasticity model, the relationship between the support pressure and displacement of the shield tunnel face was obtained. According to the plastic strain distribution at collapse state, an appropriate failure mechanism was proposed for upper bound limit analysis, and the formula to calculate the limit support pressure was deduced. The limit support pressure was rearranged to be the summation of soil cohesion c, surcharge load q, and soil gravity γ multiplied by their corresponding coefficientsN c ,N q , andN γ , and parametric studies were carried out on these coefficients. In order to consider the influence of seepage on the face stability, the pore water pressure distribution and the seepage force on the tunnel face were obtained by FE simulation. After adding the power of seepage force into the equation of the upper bound limit analysis, the total limit support pressure for stabilizing the tunnel face under seepage conditionwas obtained.The total limit support pressure was shown to increase almost linearly with the water table.


Introduction
The key issue during shield tunneling is to keep the stability of tunnel face, and this generally depends on the support pressure which was applied on the tunnel face after soil excavation. The pressure must be controlled at least no less than its limit value which corresponds to the active failure state of the tunnel face. So far, various kinds of methods have been proposed to study this problem. Experimental methods, including physical modeling and centrifuge modeling [1][2][3][4], were applied to study the failure mechanism of tunnel face and the limit support pressure. Some empirical method was proposed to evaluate the stability of the tunnel face; for example, Broms and Bennermark [5] proposed stability number to assess the stability of tunnel face in clay under undrained condition. Cornejo [6] gave out a formula to determine the stability number of the tunnel face in clay based on the limit equilibrium method. Based on the limit equilibrium of soil in front of the tunnel face, Murayama proposed a formula to calculate the limit support pressure, and the varying characteristic of the formula with soil strength was studied by Lu et al. [7]. Following Horn's theory, Jancsecz and Steiner [8] studied the face stability of shield tunnel under drained condition by considering the equilibrium of the sliding wedge at the tunnel face. By employing limit analysis, Davis et al. [9] studied the stability of tunnel face and presented lower and upper bound solutions of the limit support pressure to stabilize the tunnel face under undrained condition. Augarde et al. [10] employed finite element limit analysis to study the plane strain tunnel face under undrained condition. Leca and Dormieux [11] proposed a 3D failure mechanism to obtain the upper bound solution of limit support pressure under drained condition. After the utilization of multiblock failure mechanism, Soubra [12] and Mollon et al. [13] obtained better upper bound solutions which were closer to the experimental results of Chambon and Corté [4]. Mollon et al. [14] proposed a 2D multiblock limit analysis method to determine the critical collapse pressure of air pressurized shield tunnel, and Mollon et al. [15] further studied the face stability by probabilistic analysis. Numerical simulation method has been also adopted to study the stability of tunnel face. Vermeer et al. [16] studied the stability of the shield tunnel face and the influences of cohesion, surcharge load, and soil gravity on the limit support pressure by FEM.

Mathematical Problems in Engineering
Li et al. [17] studied the face stability of shield tunnel by FLAC 3D analysis. Discrete element method was also used for the face stability of shield tunnel [18][19][20][21]. Although the failure mechanism and limit support pressure of shield tunnel face could be obtained, the complicated calculation in numerical modeling makes it too difficult to be used in real engineering.
When shield tunnel locates under the water table line, the soil excavation often induces underground water seepage and apt to cause the collapse of shield tunnel face. The face stability analysis under seepage condition was studied by numerical simulation [22][23][24] or theoretical analysis which was based on the existing model. Anagnostou and Kovári [25] studied the influence of seepage on the stability of tunnel face based on the wedge model. de Buhan et al. [26] analyzed the face stability of shield tunnel by introducing the seepage force into the model of Leca and Dormieux [11]. Lee and Nam [27] considered the influence of seepage by superposing the results of seepage analysis on the mechanical analysis under drained condition. Lee et al. [28] compared the results obtained from theoretical analysis by taking seepage forces into account with the results of the coupled finite element analysis. Park et al. [29] studied the stability of pressurized shield tunnel by incorporating the results of Lee et al. [28] into the upper bound analysis of Leca and Dormieux [11]. The results from these works and recent study of Li et al. [30] showed that the underground water seepage played crucial role in the stability of tunnel face.
In this paper, the stability of shield tunnel face was studied by elastoplasticity FE simulation; the collapse mechanism and limit support pressure in active failure state were obtained. Based on the numerical results, a failure mechanism was proposed and a 2D upper-bound limit analysis model was established, and the formula for calculating the limit support pressure was also deduced. Following the Terzaghi superposition method which has been commonly used in bearing capacity analysis, the limit support pressure was rearranged as the summation of soil cohesion, surcharge load, and soil gravity multiplied by their corresponding coefficients, and the varying characteristics of these coefficients with the depth-to-diameter ratio of tunnel and the friction angle of soil were studied in detail. The influence of seepage on the stability of shield tunnel under water table was also studied. The pore water pressure distribution and seepage force on the shield tunnel face were obtained by FE numerical simulation. After the calculation of seepage force on the failure area of tunnel face, the proposed upper bound limit analysis model was extended into seepage condition.

Finite Element Modeling of Tunnel Face Stability
The relationship of deformation and support pressure of shield tunnel face was obtained by FE analysis with PLAXIS software, the constitutive model adopted is the widely used Mohr-Coulomb perfect elastoplasticity model. The tunnel diameter is = 10 m and the tunnel depth is = 10 m; the finite element mesh constituted by 15-node triangular element employed in numerical simulation was shown in  Mpa and ] = 0.3. Mohr-Coulomb model was adopted to describe the constitutive relationship of soil in plastic stage, the cohesion = 2 kPa, the friction angle ranges from 5 ∘ to 45 ∘ , and the soil gravity = 17 kN/m 3 . Considering that dilatancy angle has no influence on the limit support pressure [16], and in order to keep accordance with the assumption in upper bound limit analysis, the dilatancy angle is assumed to be equal to the friction angle.
The initial stress field induced by the soil weight and surcharge load was calculated. After the excavation of the soil during shield tunneling, the initial condition was recovered by applying lateral earth pressure with its value determined by 0 ( + ). In case of the lateral earth pressure 0 has no influence on the support pressure at collapse [16], 0 was set as 1 − sin( ). After the initial pressure was applied, the displacement of all nodes was set to zero. The pressure was reduced gradually from the initial value to obtain the curve of support pressure and displacement on the tunnel face. As shown in Figure 2, the support pressure decreases with the displacement, after attaining its critical value, it almost keeps constant even when the displacement keeps increasing, which indicates the collapse of the tunnel face. It is also shown in Figure 2 that the limit support pressure obviously decreases with the friction angle.
The increments of displacement and plastic strain distributions at collapse state are shown in Figures 3 and 4. The failure mode changes from global to local with the tunnel depth increases, and the deformation area around the tunnel face reduces with the increase of the friction angle of soil.

Upper Bound Limit Analysis of the Shield Tunnel Face Stability
3.1. Failure Mechanism. In order to analyze the stability of the tunnel face, an appropriate failure mechanism needs to be proposed. According to the plastic strain distribution at collapse state obtained from numerical modeling in Section 2, and referring to the Terzaghi failure mechanism for bearing capacity analysis, a failure mechanism which is composed of a shearing zone and two rigid blocks and was proposed. The proposed failure mechanism could reflect the transition from global failure mode to local mode with the increases of the tunnel depth which has been indicated by previous study [11,16]. As shown in Figure 5, the upper isosceles triangle which has an opening angle equal to 2 is block , and the axis of symmetry of the opening angle is vertical. The block is isosceles triangle ; the line has an angle of /4 + /2 with the horizontal direction. The shear zone is a log-spiral curve with the center is point . The geometric parameters of the failure mechanism are , Mathematical Problems in Engineering The comparison of the proposed failure mechanism and the numerical results is shown in Figure 4. The block moves downward with a velocity V , the block moves left with an angle of /4 + /2 with horizontal line, and the velocity in shear zone increases from V at line to V at line.

The Limit Support Pressure.
Based on the failure mechanisms proposed in previous section, the upper bound solution of the limit support pressure could be obtained by equaling the power of external force and plastic dissipation energy [31]. The power of the weight of block is The differential of the power of the weight of shearing zone is After integration, the power of the weight of shearing zone is The power of the weight of block is The power of the surcharge load is The power of the support pressure on the tunnel face is The internal energy dissipation along line and line is The internal energy dissipation along the line is The internal energy dissipation along line which equals the dissipation power rate in block is By equaling the power of external forces to the internal energy dissipation, we get By combining from (2) to (11), the formula for calculating the limit support pressure was obtained. Analogous to Terzaghi's superposition method which has been commonly used in bearing capacity analysis, the obtained limit support pressure could be rearranged to be Mathematical Problems in Engineering ⋅ cos ( 4 + 2 ) exp [( 4 + 2 ) tan ]⟩ .
As shown in Figure 6, the limit support pressure decreases with the friction angle ; it increases with the depth-to-diameter ratio of tunnel only when < 20 ∘ , while if > 20 ∘ , the variation of depth-to-diameter ratio shows almost no influence on the limit support pressure, and these results agree well with those of Vermeer et al. [16].

Parametric Studies on the Influence Coefficients.
The coefficients of cohesion, surcharge load, and soil gravity in (14) play important role in the calculation of limit support pressure and deserve further study. Different from the case of bearing capacity analysis, the coefficient of cohesion here is negative. As shown in Figure 7, increases with the friction angle of soil and decreases with the depthto-diameter ratio of shield tunnel. Vermeer et al. [16] stated that could be calculated by cot( ) when the tunnel depth is more than twice of diameter or when friction angle > 35 ∘ . The upper bound solutions are close to the solutions of Vermeer et al. [16] obtained from FEM when / = 2 and > 10 ∘ . When friction angle is less than a certain value ( = 10 ∘ ), the -curve changes, which indicates the intersection of the failure area on the ground surface. Comparatively, the formula = cot( ) could not reflect this property and is only suitable under the condition of small soil friction angles. As shown in Figure 8, is plotted as a function of the friction angle for common values of the depth ratio / equal to 0.5, 1, and 2. It is shown that the upper bound solutions decrease with the friction angle of soil and the depth-to-diameter ratio of tunnel. The value of turns to be zero when friction angle reaches a certain value, and this value decreases with the depth-to-diameter ratio of tunnel. The relationship between the upper bound solution of and friction angle of soil is shown in Figure 9. The value of decreases with , and it increases with the depth-to-diameter ratio of tunnel only when friction angle is small; otherwise it keeps constant. The upper bound solutions obtained in this paper agree well with the formula = 1/(9 tan ) − 0.05 proposed by Vermeer et al. [16] when the friction angle of soil is large. The figure also shows that the results in this paper are very close to the results of Mollon et al. [13], and it is applicable for wider range of the friction angle of soil. the seepage force. By assuming the underground water seepage to follow the Darcy law, the seepage equation in steady state is

Influence of Seepage on the Stability of Tunnel
where , are the seepage coefficients in and directions; in order to simplify the problem, these two coefficients are considered as the same; is the water table from the bottom of the tunnel.
FE simulation was employed to analyze the seepage characteristics of the tunnel face. In the simulation, two cases with tunnel diameter of 5 m and 10 m were studied; the water table varies from the top of the tunnel to three times of tunnel diameters. When tunnel diameter = 10 m, tunnel depth-to-diameter / = 2, water table-to-tunnel diameter / = 2, and the seepage coefficient is 0.3 × 10 −5 m/s; the pore water pressure distribution near the tunnel face obtained from numerical simulation is shown in Figure 10. The pore water pressure on tunnel face is zero, and the interval of pressure line is 1 m water level. The failure area was divided into two parts: the top part is a triangle, which is the same as block in Figure 4, and the lower part is blocks and . As shown in Figure 10, the pore water pressure distributes almost uniformly in the failure area of tunnel face, and it changes significantly near the top and bottom of the tunnel face.
From the pore water pressure distribution, the water head difference between the failure line and the tunnel face could be obtained, and then the seepage force on the failure area could be calculated; the detailed derivation could be found at Lee et al. [28]. The total seepage force was obtained and shown in Figure 11; the ratio of average seepage force over hydrostatic force keeps almost constant, and the value is slightly more than that of Lee et al. [28] obtained from 3D analysis.
In order to study the seepage force in detail, the horizontal and vertical components of the seepage forces in area (area ) and (areas and ) are studied separately. The total seepage force could be decomposed as where , , , and are the horizontal and vertical components of seepage force in parts and . The ratios of each component of average seepage force over hydrostatic force are shown in Figure 12; the biggest component varies slightly with the water table, and are smaller, and is the smallest.

Influence of Seepage on Limit Support
Pressure. By including the power rate of the seepage force in upper bound limit analysis, the influence of seepage on the limit support pressure of tunnel face could be studied. The power rate of the seepage force is where is the power rate of the horizontal component of seepage force; is the power rate of the vertical component of seepage force. For the horizontal component velocity in part is zero, so the power rates of the horizontal component of seepage force component vanishes; that is, The power rate of the vertical component of seepage force component is In shear zone , the power rates of the seepage force in soil element are Mathematical Problems in Engineering After the integration of (21) and (22) The power rate of the seepage force in block is By summing up, the power rates of the vertical and horizontal components of seepage forces are By noting that = / , = / and adding (24) into the left-hand side of (12), the limit support pressure is is the coefficient of seepage force and is The formula of (26) was validated when the tunnel diameters are 5 m and 10 m, and the tunnel depth is 20 m. The cohesion and friction angle of soil are 2 kPa and 30 ∘ ; the dry and saturated gravities of the soil are 17 kN/m 3 and 19 kN/m 3 . For simplicity, the surcharge load is assumed to be zero. The calculated limit support pressures in dry sand case are 11.9 kPa ( = 5m) and 27.2 kPa ( = 10m). By incorporating the seepage force in upper bound limit analysis, the upper bound solution of total limit support pressure under seepage condition can be calculated. As shown in Figure 13, the calculated total limit support pressure increases linearly with the water table and agrees very well with the results of Lee et al. [28].

Conclusions
The face stability of shield tunnel was studied by elastoplasticity FE simulation; the equivalent plastic strain distribution and limit support pressure at collapse state were obtained. Based on the numerical results, a failure mechanism was proposed to study the face stability of shield tunnel by upper bound limit analysis. The calculating formula of the limit support pressure was rearranged to be the summation of cohesion, surcharge load, and soil gravity multiplied by corresponding coefficients. Parametric analysis showed the coefficient of cohesion increases with friction angle of the soil and decreases with the depth-to-diameter ratio of tunnel. The coefficients of surcharge load and soil gravity decrease with the friction angle of soil. Both coefficients decrease with the tunnel depth-to-diameter ratio only when the friction angle is less than an appropriate value; otherwise they are independent of the depth-to-diameter ratio and the coefficient of surcharge load goes to zero. The seepage analysis was conducted by FE simulation; the pore water pressure distribution and seepage force on the tunnel face were obtained. By adding the power rate induced by seepage force, the proposed upper bound limit analysis was extended to seepage condition. The results showed that a large part of the limit support pressure was used to equilibrate the seepage force, and the total limit support pressure varied almost linearly with the water table.