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

Authors: Xilin Lu; Haoran Wang; Maosong Huang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

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 cohesionc, surcharge load q, and soil gravity γ multiplied by their corresponding coefficients Nc, Nq, and Nγ, 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 condition was obtained. The total limit support pressure was shown to increase almost linearly with the water table.


---

## Body

## 1. 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–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. 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–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–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.


## 2. 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 isD=10 m and the tunnel depth is C=10 m; the finite element mesh constituted by 15-node triangular element employed in numerical simulation was shown in Figure 1. The elastic modulus and Poisson ratio are E=20 Mpa and ν=0.3. Mohr-Coulomb model was adopted to describe the constitutive relationship of soil in plastic stage, the cohesion c=2 kPa, the friction angle φ ranges from 5° to 45°, and the soil gravity γ=17 kN/m3. 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.Figure 1
The finite element mesh for the stability analysis of shield tunnel.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 byK0(q+γC). In case of the lateral earth pressure K0 has no influence on the support pressure at collapse [16], K0 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.Figure 2
The relationship between the support pressure and displacement at center-point of the tunnel face (C/D=1).The increments of displacement and plastic strain distributions at collapse state are shown in Figures3 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.The displacement increment around the tunnel face at collapse state.
(a)
C
/
D
=
0.5 (φ=15°, 30°, 45°)
(b)
C
/
D
=
1 (φ=15°, 30°, 45°)
(c)
C
/
D
=
2 (φ=15°, 30°, 45°)The equivalent plastic strain distribution around the tunnel face at collapse state.
(a)
C
/
D
=
0.5 (φ=15°, 30°, 45°)
(b)
C
/
D
=
1 (φ=15°, 30°, 45°)
(c)
C
/
D
=
2 (φ=15°, 30°, 45°)


## 3. 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 Section2, and referring to the Terzaghi failure mechanism for bearing capacity analysis, a failure mechanism which is composed of a shearing zone b and two rigid blocks a and c 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 O′OB which has an opening angle equal to 2φ is block a, and the axis of symmetry of the opening angle is vertical. The block c is isosceles triangle OAA′; the line AA′ has an angle of π/4+φ/2 with the horizontal direction. The shear zone b is a log-spiral curve with the center is point O. The geometric parameters of the failure mechanism are
(1)r0=D2sin⁡(π/4+φ/2)exp⁡[(π/4+φ/2)tan⁡(φ)],h2=r0(2tanφ),lB={0h-C≤02(h-C)tan⁡φh-C>0.
The comparison of the proposed failure mechanism and the numerical results is shown in Figure 4. The block a moves downward with a velocity va, the block c moves left with an angle of π/4+φ/2 with horizontal line, and the velocity in shear zone increases from va at OB line to vc at OA′ line.Figure 5
The proposed failure mechanism.


### 3.2. 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 a is
(2)PWa=12γva[r0h-lB(h-C)].The differential of the power of the weight of shearing zoneb is
(3)dPwb=va2γ∫0π/4+φ/2r02exp⁡⁡(3θtan⁡φ)cos⁡⁡θdθ.After integration, the power of the weight of shearing zoneb is
(4)Pwb=∫0π/4+φ/2dPwb=γvar022(1+9tan2φ)·(sin⁡(π4+φ2)exp⁡⁡[3(π4+φ2)tan⁡φ]mii+3tan⁡φ{cos⁡⁡(π4+φ2)miimmmmmm·exp⁡⁡[3(π4+φ2)tan⁡φ]-1}).The power of the weight of blockc is
(5)PWc=γ4D2tan⁡(π4-φ2)cos⁡⁡(π4+φ2)va·exp⁡⁡[(π4+φ2)tan⁡φ].The power of the surcharge load is(6)Pq=qlBva.The power of the support pressure on the tunnel face is(7)Pt=σtDvaexp⁡⁡[(π4+φ2)tan⁡φ]sin⁡(π4+φ2).The internal energy dissipation along lineOO′ and line O′B′ is
(8)Ea=(2h-lBcot⁡φ)cva.The internal energy dissipation along the lineAA′ is
(9)Ec=Dexp⁡[(π/4+φ/2)tan⁡φ]2sin⁡(π/4+φ/2)cva.The internal energy dissipation along lineBA′ which equals the dissipation power rate in block b is
(10)Eb=EBA′=12cvar0cot⁡φ{exp⁡⁡[(π2+φ)tan⁡φ]-1}.By equaling the power of external forces to the internal energy dissipation, we get(11)PWa+PWb+PWc+Pq-Pt=Ea+Eb+Ec+EBA′.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
(12)σt=cNc+qNq+γDNγ,
where the influence coefficients of cohesion Nc, surcharge load Nq, and soil gravity Nγ are
(13)Nc=-1Dexp⁡[(π/4+φ/2)tan⁡φ]sin⁡(π/4+φ/2)·〈(2h-lBcot⁡φ)π2iiiiii+Dexp⁡[(π/4+φ/2)tan⁡φ]cos⁡⁡φ2sin⁡(π/4+φ/2)iiiiii+r0cot⁡φ{exp⁡⁡[(π2+φ)tan⁡φ]-1}〉Nq=lBDexp⁡[(π/4+φ/2)tan⁡φ]sin⁡(π/4+φ/2)Nγ=1D2exp⁡[(π/4+ϕ/2)tan⁡ϕ]sin⁡(π/4+ϕ/2)·〈r0h-lB(h-C)D22iiiii+((3tan⁡ϕ{exp⁡⁡[3(π4+ϕ2)tan⁡ϕ]cos⁡⁡(π4+ϕ2)-1}mmmiiiii·{exp⁡⁡[3(π4+ϕ2)tan⁡ϕ]mmmmmim·cos⁡⁡(π4+ϕ2)-1}mmmiiiii+(3tan⁡ϕ{exp⁡⁡[3(π4+ϕ2)tan⁡ϕ]cos⁡⁡(π4+ϕ2)-1}exp⁡⁡[3(π4+ϕ2)tan⁡ϕ]sin⁡(π4+ϕ2)){exp⁡⁡[3(π4+ϕ2)tan⁡ϕ]cos⁡⁡(π4+ϕ2)-1}iiiiiiiiiii·(1+9tan2ϕ)-1)r02iiiii+D22tan⁡(π4-ϕ2)iiiii·cos⁡⁡(π4+ϕ2)exp⁡⁡[(π4+ϕ2)tan⁡ϕ]〉.As shown in Figure6, 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].The influence of the soil strength parameters on the limit support pressure (γ=17 kPa, D=10 m).
(a)
The influence of cohesion
(b)
The influence of friction angle (c=2kPa)


### 3.3. 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 Nc here is negative. As shown in Figure 7, Nc increases with the friction angle φ of soil and decreases with the depth-to-diameter ratio of shield tunnel. Vermeer et al. [16] stated that Nc could be calculated by cot⁡(φ) when the tunnel depth C is more than twice of diameter D or when friction angle φ>35°. The upper bound solutions are close to the solutions of Vermeer et al. [16] obtained from FEM when C/D=2 and φ>10°. When friction angle is less than a certain value (φ=10°), the Nc-φ curve changes, which indicates the intersection of the failure area on the ground surface. Comparatively, the formula Nc=cot⁡(φ) could not reflect this property and is only suitable under the condition of small soil friction angles. As shown in Figure 8, Nq is plotted as a function of the friction angle for common values of the depth ratio C/D 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 Nq 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 Nγ and friction angle φ of soil is shown in Figure 9. The value of Nγ 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 Nγ=1/(9tan⁡φ)-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.Figure 7
The relationship betweenNc and the friction angle of soil (D=10m).Figure 8
The relationship betweenNq and the friction angle of soil (D=10m).Figure 9
The relationship betweenNγ and the friction angle of soil (D=10m).


## 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 Section2, and referring to the Terzaghi failure mechanism for bearing capacity analysis, a failure mechanism which is composed of a shearing zone b and two rigid blocks a and c 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 O′OB which has an opening angle equal to 2φ is block a, and the axis of symmetry of the opening angle is vertical. The block c is isosceles triangle OAA′; the line AA′ has an angle of π/4+φ/2 with the horizontal direction. The shear zone b is a log-spiral curve with the center is point O. The geometric parameters of the failure mechanism are
(1)r0=D2sin⁡(π/4+φ/2)exp⁡[(π/4+φ/2)tan⁡(φ)],h2=r0(2tanφ),lB={0h-C≤02(h-C)tan⁡φh-C>0.
The comparison of the proposed failure mechanism and the numerical results is shown in Figure 4. The block a moves downward with a velocity va, the block c moves left with an angle of π/4+φ/2 with horizontal line, and the velocity in shear zone increases from va at OB line to vc at OA′ line.Figure 5
The proposed failure mechanism.


## 3.2. 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 a is
(2)PWa=12γva[r0h-lB(h-C)].The differential of the power of the weight of shearing zoneb is
(3)dPwb=va2γ∫0π/4+φ/2r02exp⁡⁡(3θtan⁡φ)cos⁡⁡θdθ.After integration, the power of the weight of shearing zoneb is
(4)Pwb=∫0π/4+φ/2dPwb=γvar022(1+9tan2φ)·(sin⁡(π4+φ2)exp⁡⁡[3(π4+φ2)tan⁡φ]mii+3tan⁡φ{cos⁡⁡(π4+φ2)miimmmmmm·exp⁡⁡[3(π4+φ2)tan⁡φ]-1}).The power of the weight of blockc is
(5)PWc=γ4D2tan⁡(π4-φ2)cos⁡⁡(π4+φ2)va·exp⁡⁡[(π4+φ2)tan⁡φ].The power of the surcharge load is(6)Pq=qlBva.The power of the support pressure on the tunnel face is(7)Pt=σtDvaexp⁡⁡[(π4+φ2)tan⁡φ]sin⁡(π4+φ2).The internal energy dissipation along lineOO′ and line O′B′ is
(8)Ea=(2h-lBcot⁡φ)cva.The internal energy dissipation along the lineAA′ is
(9)Ec=Dexp⁡[(π/4+φ/2)tan⁡φ]2sin⁡(π/4+φ/2)cva.The internal energy dissipation along lineBA′ which equals the dissipation power rate in block b is
(10)Eb=EBA′=12cvar0cot⁡φ{exp⁡⁡[(π2+φ)tan⁡φ]-1}.By equaling the power of external forces to the internal energy dissipation, we get(11)PWa+PWb+PWc+Pq-Pt=Ea+Eb+Ec+EBA′.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
(12)σt=cNc+qNq+γDNγ,
where the influence coefficients of cohesion Nc, surcharge load Nq, and soil gravity Nγ are
(13)Nc=-1Dexp⁡[(π/4+φ/2)tan⁡φ]sin⁡(π/4+φ/2)·〈(2h-lBcot⁡φ)π2iiiiii+Dexp⁡[(π/4+φ/2)tan⁡φ]cos⁡⁡φ2sin⁡(π/4+φ/2)iiiiii+r0cot⁡φ{exp⁡⁡[(π2+φ)tan⁡φ]-1}〉Nq=lBDexp⁡[(π/4+φ/2)tan⁡φ]sin⁡(π/4+φ/2)Nγ=1D2exp⁡[(π/4+ϕ/2)tan⁡ϕ]sin⁡(π/4+ϕ/2)·〈r0h-lB(h-C)D22iiiii+((3tan⁡ϕ{exp⁡⁡[3(π4+ϕ2)tan⁡ϕ]cos⁡⁡(π4+ϕ2)-1}mmmiiiii·{exp⁡⁡[3(π4+ϕ2)tan⁡ϕ]mmmmmim·cos⁡⁡(π4+ϕ2)-1}mmmiiiii+(3tan⁡ϕ{exp⁡⁡[3(π4+ϕ2)tan⁡ϕ]cos⁡⁡(π4+ϕ2)-1}exp⁡⁡[3(π4+ϕ2)tan⁡ϕ]sin⁡(π4+ϕ2)){exp⁡⁡[3(π4+ϕ2)tan⁡ϕ]cos⁡⁡(π4+ϕ2)-1}iiiiiiiiiii·(1+9tan2ϕ)-1)r02iiiii+D22tan⁡(π4-ϕ2)iiiii·cos⁡⁡(π4+ϕ2)exp⁡⁡[(π4+ϕ2)tan⁡ϕ]〉.As shown in Figure6, 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].The influence of the soil strength parameters on the limit support pressure (γ=17 kPa, D=10 m).
(a)
The influence of cohesion
(b)
The influence of friction angle (c=2kPa)


## 3.3. 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 Nc here is negative. As shown in Figure 7, Nc increases with the friction angle φ of soil and decreases with the depth-to-diameter ratio of shield tunnel. Vermeer et al. [16] stated that Nc could be calculated by cot⁡(φ) when the tunnel depth C is more than twice of diameter D or when friction angle φ>35°. The upper bound solutions are close to the solutions of Vermeer et al. [16] obtained from FEM when C/D=2 and φ>10°. When friction angle is less than a certain value (φ=10°), the Nc-φ curve changes, which indicates the intersection of the failure area on the ground surface. Comparatively, the formula Nc=cot⁡(φ) could not reflect this property and is only suitable under the condition of small soil friction angles. As shown in Figure 8, Nq is plotted as a function of the friction angle for common values of the depth ratio C/D 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 Nq 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 Nγ and friction angle φ of soil is shown in Figure 9. The value of Nγ 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 Nγ=1/(9tan⁡φ)-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.Figure 7
The relationship betweenNc and the friction angle of soil (D=10m).Figure 8
The relationship betweenNq and the friction angle of soil (D=10m).Figure 9
The relationship betweenNγ and the friction angle of soil (D=10m).


## 4. Influence of Seepage on the Stability of Tunnel Face

### 4.1. Seepage Analysis

When shield tunnel locates under water table line, the soil excavation often induces seepage and causes the failure of tunnel face. It is necessary to establish a calculation model to estimate the influence of seepage on the stability of tunnel face. The key point of the stability analysis under seepage condition is the calculation of the seepage force. By assuming the underground water seepage to follow the Darcy law, the seepage equation in steady state is(14)∂∂x(Kx∂H∂x)+∂∂y(Ky∂H∂y)=0,
where Kx, Ky are the seepage coefficients in x and y directions; in order to simplify the problem, these two coefficients are considered as the same; H 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 diameterD=10 m, tunnel depth-to-diameter C/D=2, water table-to-tunnel diameter H/D=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 A is a triangle, which is the same as block a in Figure 4, and the lower part B is blocks b and c. 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.Figure 10
The distribution of the pore water pressure around 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 forceF 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.Figure 11
The relationship between the average seepage force and the water table.In order to study the seepage force in detail, the horizontal and vertical components of the seepage forces in areaA (area a) and B (areas b and c) are studied separately. The total seepage force could be decomposed as
(15)F=(FAx+FBx)2+(FAy+FBy)2,
where FAx, FAy, FBx, and FBy are the horizontal and vertical components of seepage force in parts A and B.The ratios of each component of average seepage force over hydrostatic force are shown in Figure12; the biggest component FBx varies slightly with the water table, FAy and FBy are smaller, and FAx is the smallest.Figure 12
The relationship between the seepage force and the water table.


### 4.2. 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 forcePF is
(16)PF=PFx+PFy,
where PFx is the power rate of the horizontal component of seepage force; PFy is the power rate of the vertical component of seepage force. For the horizontal component velocity in part A is zero, so the power rates of the horizontal component of seepage force component vanishes; that is,
(17)PFAx=0.
The power rate of the vertical component of seepage force component is
(18)PFAy=12FAyva[r0h-lB(h-C)].
In shear zone b, the power rates of the seepage force in soil element are
(19)dPFbx=va2FBxr02exp⁡⁡(3θtan⁡φ)sin⁡θdθ,dPFby=va2FByr02exp⁡⁡(3θtan⁡φ)cos⁡⁡θdθ.After the integration of (21) and (22), we get
(20)PFbx=FBxvar022(1+9tan2φ)·{[-cos⁡⁡(π4+φ2)+3tan⁡φsin⁡(π4+φ2)]iiii·exp⁡⁡[3(π4+φ2)tan⁡φ]+1},(21)PFby=FByvar022(1+9tan2φ)·{[sin⁡(π4+φ2)+3tan⁡φcos⁡⁡(π4+φ2)]iiii·exp⁡⁡[3(π4+φ2)tan⁡φ]-3tan⁡φ}.
The power rate of the seepage force in block c is
(22)PFcx=FBx4D2tan⁡(π4-φ2)sin⁡(π4+φ2)va·exp⁡⁡[(π4+φ2)tan⁡φ],(23)PFcy=FBy4D2tan⁡(π4-φ2)cos⁡⁡(π4+φ2)va·exp⁡⁡[(π4+φ2)tan⁡φ].
By summing up, the power rates of the vertical and horizontal components of seepage forces are
(24)PFx=PFax+PFbx+PFcx,PFy=PFay+PFby+PFcy.
By noting that α=FAy/FBx, β=FBy/FBx and adding (24) into the left-hand side of (12), the limit support pressure is
(25)σt=cNc+qNq+γ′DNγ+FBxNf.Nf is the coefficient of seepage force and is
(26)NF=NFBx+αNFAy+βNFBy,
where
(27)NFBx=12D2exp⁡⁡[(π/4+ϕ/2)tan⁡ϕ]sin⁡(π/4+ϕ/2)·{r02·(([3tan⁡ϕsin⁡(π4+ϕ2)-cos⁡⁡(π4+ϕ2)]mmmimm·exp⁡⁡[3(π4+ϕ2)tan⁡ϕ]+1)mmmnm·[3tan⁡ϕsin⁡(π4+ϕ2)-cos⁡⁡(π4+ϕ2)](1+9tan2ϕ)-1)mm+D22tan⁡(π4-ϕ2)sin⁡(π4+ϕ2)mm·exp⁡⁡[(π4+ϕ2)tan⁡ϕ]},NFAy=r0h-lB(h-C)2D2exp⁡⁡[(π/4+ϕ/2)tan⁡ϕ]sin⁡(π/4+ϕ/2),NFBy=Nγ-NFAy.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/m3 and 19 kN/m3. For simplicity, the surcharge load is assumed to be zero. The calculated limit support pressures in dry sand case are 11.9 kPa (D=5 m) and 27.2 kPa (D=10 m). 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].Figure 13
The relationship between the total limit support pressure and the water table.


## 4.1. Seepage Analysis

When shield tunnel locates under water table line, the soil excavation often induces seepage and causes the failure of tunnel face. It is necessary to establish a calculation model to estimate the influence of seepage on the stability of tunnel face. The key point of the stability analysis under seepage condition is the calculation of the seepage force. By assuming the underground water seepage to follow the Darcy law, the seepage equation in steady state is(14)∂∂x(Kx∂H∂x)+∂∂y(Ky∂H∂y)=0,
where Kx, Ky are the seepage coefficients in x and y directions; in order to simplify the problem, these two coefficients are considered as the same; H 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 diameterD=10 m, tunnel depth-to-diameter C/D=2, water table-to-tunnel diameter H/D=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 A is a triangle, which is the same as block a in Figure 4, and the lower part B is blocks b and c. 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.Figure 10
The distribution of the pore water pressure around 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 forceF 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.Figure 11
The relationship between the average seepage force and the water table.In order to study the seepage force in detail, the horizontal and vertical components of the seepage forces in areaA (area a) and B (areas b and c) are studied separately. The total seepage force could be decomposed as
(15)F=(FAx+FBx)2+(FAy+FBy)2,
where FAx, FAy, FBx, and FBy are the horizontal and vertical components of seepage force in parts A and B.The ratios of each component of average seepage force over hydrostatic force are shown in Figure12; the biggest component FBx varies slightly with the water table, FAy and FBy are smaller, and FAx is the smallest.Figure 12
The relationship between the seepage force and the water table.


## 4.2. 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 forcePF is
(16)PF=PFx+PFy,
where PFx is the power rate of the horizontal component of seepage force; PFy is the power rate of the vertical component of seepage force. For the horizontal component velocity in part A is zero, so the power rates of the horizontal component of seepage force component vanishes; that is,
(17)PFAx=0.
The power rate of the vertical component of seepage force component is
(18)PFAy=12FAyva[r0h-lB(h-C)].
In shear zone b, the power rates of the seepage force in soil element are
(19)dPFbx=va2FBxr02exp⁡⁡(3θtan⁡φ)sin⁡θdθ,dPFby=va2FByr02exp⁡⁡(3θtan⁡φ)cos⁡⁡θdθ.After the integration of (21) and (22), we get
(20)PFbx=FBxvar022(1+9tan2φ)·{[-cos⁡⁡(π4+φ2)+3tan⁡φsin⁡(π4+φ2)]iiii·exp⁡⁡[3(π4+φ2)tan⁡φ]+1},(21)PFby=FByvar022(1+9tan2φ)·{[sin⁡(π4+φ2)+3tan⁡φcos⁡⁡(π4+φ2)]iiii·exp⁡⁡[3(π4+φ2)tan⁡φ]-3tan⁡φ}.
The power rate of the seepage force in block c is
(22)PFcx=FBx4D2tan⁡(π4-φ2)sin⁡(π4+φ2)va·exp⁡⁡[(π4+φ2)tan⁡φ],(23)PFcy=FBy4D2tan⁡(π4-φ2)cos⁡⁡(π4+φ2)va·exp⁡⁡[(π4+φ2)tan⁡φ].
By summing up, the power rates of the vertical and horizontal components of seepage forces are
(24)PFx=PFax+PFbx+PFcx,PFy=PFay+PFby+PFcy.
By noting that α=FAy/FBx, β=FBy/FBx and adding (24) into the left-hand side of (12), the limit support pressure is
(25)σt=cNc+qNq+γ′DNγ+FBxNf.Nf is the coefficient of seepage force and is
(26)NF=NFBx+αNFAy+βNFBy,
where
(27)NFBx=12D2exp⁡⁡[(π/4+ϕ/2)tan⁡ϕ]sin⁡(π/4+ϕ/2)·{r02·(([3tan⁡ϕsin⁡(π4+ϕ2)-cos⁡⁡(π4+ϕ2)]mmmimm·exp⁡⁡[3(π4+ϕ2)tan⁡ϕ]+1)mmmnm·[3tan⁡ϕsin⁡(π4+ϕ2)-cos⁡⁡(π4+ϕ2)](1+9tan2ϕ)-1)mm+D22tan⁡(π4-ϕ2)sin⁡(π4+ϕ2)mm·exp⁡⁡[(π4+ϕ2)tan⁡ϕ]},NFAy=r0h-lB(h-C)2D2exp⁡⁡[(π/4+ϕ/2)tan⁡ϕ]sin⁡(π/4+ϕ/2),NFBy=Nγ-NFAy.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/m3 and 19 kN/m3. For simplicity, the surcharge load is assumed to be zero. The calculated limit support pressures in dry sand case are 11.9 kPa (D=5 m) and 27.2 kPa (D=10 m). 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].Figure 13
The relationship between the total limit support pressure and the water table.


## 5. 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.



---

*Source: 727964-2014-10-09.xml*