An Analytical Continuous Upper Bound Limit Analysis of Pore Water Effect on the Tail Stability of Underwater Shield Tunnels during Construction

An analytical continuous upper bound limit analysis is developed to analyse the effects of seepage on the transverse stability of underwater shield tunnels. The approach is based on an analytical continuous upper bound limit analysis method for cohesive-frictional soils. It employs the complex variables solution of the displacement field due to tunnel deformation and movement, and the analytical solution of the pore water pressure field for steady state seepage due to pore water influx at the tunnel perimeter. The most critical slip line position and the minimum required tunnel support pressure are determined by using a particle swarm optimization scheme for various generic situations. The method is verified via finite element simulation and comparison with the solution from using rigid block upper bound limit analysis. The parametric analysis revealed among other things that both the infimum of the necessary tunnel support pressure and the most critical plastic zone increase when the hydraulic head at the ground surface increases, but decrease when the tunnel influx increases due to the fact that pore water pressure at the tunnel perimeter decreases with the tunnel influx.


Introduction
Shield driven tunnels are common for under crossings of surface waters, such as lakes, rivers, and seas. Construction analysis of underwater shield tunnels for project safety and environmental stability should take into account of the influence of pore water pressure and inflow on the rational support pressures at both the face and the tail of the tunnel. The face support pressure is provided by the slurry or the muck in the work chamber of a slurry shield or an earth pressure balanced shield, respectively, and the tail support pressure by synchronized grouting into the annulus gap between excavation boundary and the outer perimeter of the segment lining.
Stability problems for shield tunneling have been studied via, in addition to laboratory experiments and filed tests, a number of different analytical and numerical approaches, which mainly include the (analytical) limit equilibrium method, the (analytical) rigid body upper bound limit analysis method, the (analytical) continuous upper bound limit analysis method, the (numerical) finite element upper bound limit analysis method, and the (numerical) discontinuity layout optimization method.
The continuous upper bound limit analysis method, which is sometimes also dubbed as mobilized strength design, is special mainly in two aspects. First, the method differs from the limit equilibrium method and the rigid body upper bound limit analysis method in that it may take into account of the continuities of both the ground material properties and the ground displacements, and that it does not need to prescribe any failure modes. Second, the method is, in contrast to the finite element upper bound limit analysis method and the discontinuity layout upper bound limit analysis method, analytical in the formulation, and simple and transparent in the calculation. There exist publications of the analytical continuous upper bound limit analysis method for shield tunnel stability problems, such as Osman et al. [1][2], Klar et al. [3], Mollon et al. [4], Klar and Klein [5], Xiang and Song [6], Zhang et al. [7], Huang et al. [8], Song and Xiang [9], Xiang and Song [10], however, none of them and other analytical continuous upper bound limit analysis developments in the existent literature has taken into consideration of groundwater seepage effect. In contrast, the limit equilibrium method and the rigid body upper bound limit analysis method have been widely developed with pore water effect, such as Anagnostou and Kovári [11], Lee et al. [12], Anagnostou [13], and Perazzelli [14], for the former, and Lee et al. [15], Park et al. [16], Huang et al. [17][18][19], Xie et al. [20], Pan and Dias [21], Liu et al. [22], and Sun et al. [23], for the latter.
In this paper, taking account of pore water effect, two types of formulations of the upper bound theorem are presented, and the necessary condition for the tail stability during construction of underwater shield tunnels is developed, in which the ground displacement velocity field is calculated by using a complex variables solution for tunnel convergence, oval deformation, and vertical translation, and the pore water pressure field is calculated by using an analytical solution for steady state seepage. The constraints due to the Mohr-Coulomb criterion and the associated flow rule, and the solution scheme for the analytical continuous upper bound limit analysis formulation are described. Verifications of the developed method are presented via finite element simulation and comparison with rigid block upper bound limit analysis. Parametric analysis with regard to the soil's friction angle, the tunnel cover to diameter ratio, the tunnel influx, and the surface water level are also presented.

Upper bound theorem with pore water effect
For any geotechnical stability problem with seepage effect in general, as described for slopes by Kim et al. [24], Michalowski [25], Miller & Hamilton [26], and Miller & Hamilton [27], the necessary condition for active stability of the ground can be stated via the upper bound theorem of limit analysis as ( ) or alternatively as ( ) where, the equal signs hold when at critical state;V is the assumed plastic volume of the physical domain; A is the natural surface of V ; S is the assumed plastic slip surface of V ; w p is the pore water pressure; pi f is the i -direction seepage force due to the pore water pressure gradient, The inequality in Equation (1) employs in its left hand side the plastic rates of work by the dry traction, the dry bulk body force, and the pore water related pressure and seepage force, and in its right hand side the sum of the rates of energy dissipation with the effective stress in the plastic volume and the effective traction on the plastic slip surface, E ′  ; thus the pore water related pressure and seepage force are treated as external forces.
On the other hand, the inequality in Equation (2) employs in its left hand side the plastic rates of work by the dry traction and the dry bulk body force only (without the pore water related pressure and seepage force), and in its right hand side the sum of the rates of energy dissipation with the total stress in the plastic volume and the total traction on the plastic slip surface, E  ; thus the pore water related pressure and seepage force are treated as internal forces.
These two inequalities in Equation (1) and in Equation (2) (3) is directly based on the two relationships, i.e. the relationship of the total stress, the effective stress, and the pore water pressure in the plastic volume, and the relationship of the total traction, the effective traction, and the pore water pressure on the plastic slip surface; Equation (4) can be derived by using the divergence theorem along with the geometric relationship between the plastic displacement velocity and the strain rate.
Using the following relationship of the pore water pressure, the elevation head, and the total head, Equation (4) can be written as where, w γ is the unit weight of pore water, y is the elevation head, w h is the total head, hi f is the i -direction seepage force due to the hydraulic It is obvious that the three terms in the right hand side of Equation (6) sequentially represents the plastic rates of work by the pore water pressure, the head gradient driven seepage force, and the buoyant force; note that,  , is the vertically upward plastic displacement velocity.

Necessary stability condition for underwater shield tunnels
Referring to Fig.1, let's consider the cross sectional tail stability problem of an underwater shield tunnel of diameter, D , and overburden thickness, C , in a homogeneous and isotropic cohesive-frictional soil with the dry bulk unit weight, Young's modulus, Poisson's ratio, effect cohesion, and effective angle of internal friction being denoted as γ , E , ν , c' , and φ ′ , respectively. According to Equation (1) or Equation (2), the necessary condition for active stability of the soil can be expressed as or alternatively as where, the equal signs hold when at critical state; px f and py f are the x direction and the y direction seepage forces due to the pore water pressure gradient, respectively; V is the volume (area times a unit length along the tunnel drive) of the assumed plastic zone; S is the assumed plastic slip surface (curve length times a unit length along the tunnel drive) the necessary tunnel support pressure (radially outward); for the assumed plastic kinematically admissible displacement velocity field,  y u is the distributed displacement velocity in the y direction in the plastic zone,  t u is the distributed radially inward displacement velocity of the tunnel circumference,  s u and  n u are the distributed tangential and normal displacement velocity discontinuities along the slip surface, respectively, and max γ is the distributed maximum shear strain rate. Now let's use Equation (9), assuming constant support pressure on the tunnel circumference, the necessary condition for stability can be derived as where, -σ t l stands for a lower bound of σ t and --max which corresponds to the minimum total rate of energy dissipation, m in E  .

Displacement velocity due to tunnel deformation and movement
Referring to Fig.2, we use the cross sectional plane strain composite volume loss model originated from Sagaseta [28] and expanded in Song and Xiang [9], which takes into account of convergence, oval deformation, and vertical translation of the tunnel; ε where θ stands for the counterclockwise angle of the outward radial direction with respect to the right springline of the tunnel (the x direction).
Take the following conformal mapping and use the relationships of θ with x and y , Equation e e e e (13) where u ρ  denotes the radially outward displacement According to the Kolosov-Muskhelishavili equation of complex variables in elasticity theory, the ground displacement velocities, x u  and x u  , due to the volume loss of the tunnel can be expressed as are the rates of the complex potentials, and can be written as Laurent series k= , 2 ,…) are the coefficient rates to be determined by using the boundary conditions, which include the free of stress condition at the ground surface, and the composite volume loss displacement velocity condition at the tunnel perimeter, as expressed by Equation (13)

Pore water pressure due to steady state tunnel influx
For steady state seepage, as described in Kolymbas and Wagner [29] and Xiang [30], with hydraulic heads given at the tunnel perimeter and the ground surface, complex variables method can be employed to solve the governing mass conservation equation to obtain the pore water pressure w w w  (16) and the tunnel influx per unit distance along the tunnel drive Where,   Referring to Sloan [31], take the following squared form of equation (18) (20), the associated flow rule and its approximation can be expressed as where λ  denotes the plastic rate multiplier, λ  k denotes the plastic rate multiplier rate for the where  n u and  s u are the normal and tangential components of the plastic displacement velocity, respectively; note that the normal and tangential components of the plastic displacement velocity for the outer side of the slip surface are presumed zero.

Solution of the continuous upper bound limit analysis problem
As can be seen in Equation (10), the task of determining the infimum of the necessary tunnel support pressure, --σ t l max , can be cast as finding the minimum total rate of energy dissipation, min  E . The analytical continuous upper bound limit analysis formulation here and its solution process are similar to Xiang and Song [10] with the difference in that the rates of work by the pore water pressure and the seepage force now need to be taken into consideration. The formulation uses the tunnel volume loss parameters, ε  u , α δ , v α , the starting position of the slip surface at the ground surface, x l , and the plastic rate multipliers, λ  k ( = k 1,2,…, p ), as the decision variables. The solution process involves a tracking scheme for identifying the most critical slip surface, and takes into account of, in an analytical and continuous manner, the flow rule constraint for the plastic zone, the flow rule constraint for the slip surface, the displacement velocity boundary condition constraint, the position constraint for the slip surface node at the ground surface, and the tunnel volume loss constraint. The minimization is conducted via particle swarm optimization. For further details of the tracking scheme, the set of constraints and the optimization formulation, the readers are referred to Xiang and Song [10].

Verification by numerical simulation
In order to verify the solutions from the analytical continuous upper bound limit analysis by numerical simulations, the finite element software MIDAS GTS/NX is utilized to simulate the behaviors of the surrounding cohesive-frictional soil under different tunnel support pressures. The finite element mesh is shown in Fig. 4. Note that the left hand and right hand side boundaries located ten times of the tunnel diameter away from the tunnel center are vertically aligned roller boundaries with a constant hydraulic head, the bottom boundary located six times of the tunnel diameter away from the tunnel center is a clamped boundary with a constant hydraulic head, and the ground surface is a stress-free boundary (considering cases without ground surface load) with a constant hydraulic head. The parameter values of the tunnel and the ground are given as, unless indicated otherwise, in Table 1. It can be seen in Fig.6a that when there is no tunnel support, plastic zone appears. It can be seen in Fig.6b that when there is tunnel support but the support pressure is smaller than the infimum of the necessary tunnel support pressure, plastic zone albeit much smaller still appears. As shown in Fig.6c, when there is tunnel support and the support pressure is equal to the infimum of the necessary tunnel support pressure, there is no plastic zone. Thus, the finite element simulations have verified the rationality of the infimum of the necessary support pressure obtained from the analytical continuous upper bound limit analysis with pore water effect for underwater shield tunnels. For the situation without tunnel support, as shown in Fig.7a and Fig.7b, and also in Fig.6a, the tunnel influx may cause the plastic zone more concentrated directly above and below the tunnel but narrower sideways, and the larger the tunnel influx, the smaller the plastic zone, which can be attributed to the fact that the larger the tunnel influx, the smaller the difference between the major and the minor principal effective stresses. The plastic zone is narrowly X shaped, as expected for a shallow tunnel and a soil obeying the Mohr-Coulomb rule, with the upper wings being somewhat longer than the lower wings, and symmetrical with respect to the tunnel in the horizontal direction. Note that although the wedge region contained in the upper dihedral wings of the plastic zone does not satisfy the Mohr-Coulomb rule, the involved displacements can still be very large.

Comparison with rigid block upper bound limit analysis solution
We employ the 6-variable failure mechanism shown in Fig.8 for a cross section of unit thickness of a tunnel and its surrounding ground. With symmetry in mind, let 0 V , 1 V , 2 V , and 3 V be the displacement velocities of the rigid blocks designated with their serial numbers as the subscripts; 0,1 V , 1,2 V , and 2,3 V be the inter-block relative velocities designated with their coupled serial numbers as the subscripts; and 1 α , 2 α , 3 α , 1 β , 2 β , and 3 β be the angular parameters to be optimized upon for minimization of the total rate of energy dissipation in the failure mechanism, which are subjected to the kinematic admissibility constraints shown in Fig.9. Note that where Ω stands for the area of a domain designated by the subscript, e.g.
where t σ stands for the uniformly distributed necessary tunnel support pressure.
The total rate of work by the pore water pressure on the slip surfaces is  Fig.8. Conceptual model of a 6-varible rigid block failure mechanism for underwater tunnels The total rate of energy dissipation on the slip surfaces in the failure mechanism is 4 3 Now, by employing the upper bound theorem of limit analysis, a lower bound of the necessary tunnel support pressure t σ can be expressed as ( ) 29) and the infimum of t σ is ( ) where min  S E is the minimum total rate of energy dissipation in the failure mechanism.  Table 2 and Table 3 are the comparisons between the results from the analytical continuous upper bound limit analysis and those from the 6-variable rigid block upper bound limit analysis, for the situation where the tunnel influx varies while the hydraulic head at the ground surface is held constant, and the situation where the hydraulic head at the ground surface varies while the tunnel influx remains constant, respectively. It can be seen that the relative differences between the solutions from the two different limit analysis methods are smaller than 8%.

Parametric analysis
Using the parameter values given in Table 1, unless indicated otherwise, a series of analytical continuous upper bound limit analysis calculations are carried out to examine the effects of various parameters on the solutions. It can be seen in Fig.10 that the infimum of the necessary tunnel support pressure, --max It can be seen in Fig.11 that the infimum of the necessary tunnel support pressure, --max t l σ , increases almost linearly with the tunnel cover-diameter ratio, / C D ,when both the hydraulic head at the ground surface, w H , and the tunnel in flux, Q , are held constant, and that such effect can be slightly more prominent if the effective angle of internal friction is smaller, e.g. 15 φ°′ = .
It can be seen in Fig.12 that the infimum of the necessary tunnel support pressure, --max t l σ , decreases linearly with the increase of the tunnel influx, Q , when both the tunnel cover-diameter ratio, / C D , and the hydraulic head at the ground surface, w H , are held constant, and that such effect is attributable to the fact that pore water pressure at the tunnel perimeter decreases with the increase of the tunnel influx, as illustrated in Fig.13 in terms of the relationship between the piezometer pressure at the tunnel perimeter, w t H γ , and the tunnel influx, Q .
As shown in Fig.14, the most critical plastic zone delineated by the most critical slip surface for a given soil, decreases with increase of the tunnel influx, Q , when both the tunnel cover-diameter ratio, / C D , and the hydraulic head at the ground surface, w H , are held constant. As shown in Fig.15, the most critical plastic zone for a given soil, increases with the increase of the hydraulic head at the ground surface, w H , when both the tunnel cover-diameter ratio, / C D , and the tunnel influx, Q , are held constant.
As shown in Fig.16, Fig.17, and Fig.18, the larger the angle of internal friction of the soil, the smaller the surface range and the area of the plastic zone (as well as the corresponding infimum of the necessary tunnel support pressure, not shown here), regardless of the tunnel cover-diameter ratio. It is also notable that when the angle of internal friction is small, e.g. 15 φ°′ = , the soil may behave more like, as far as the most critical state is concerned, a pure cohesive soil, especially so when tunnel cover-diameter ratio is not too large, e.g.

Conclusion
The upper bound theorem of limit analysis with pore water effect can be cast in two types of forms, one with the pore water related forces being treated as external forces, and the other with the pore water related forces being treated as internal forces; these two types of forms are different in appearance but identical in nature. An analytical formulation of continuous upper bound limit analysis is presented for studying the effects of seepage on the transverse stability of underwater shield tunnels. Both the displacement field due to tunnel deformation-movement and the pore water pressure field for steady state seepage due to pore water pressure relaxation at the tunnel perimeter are determined by using the complex variables method. The formulation is solved by using a particle swarm optimization scheme to obtain the most critical slip line position and the minimum required tunnel support pressure for a number of generic examples.
The formulation and its solution method is verified via finite element simulation and comparison with the solution from using rigid block upper bound limit analysis. The parametric analysis indicate that for underwater tunnels the infimum of the necessary tunnel support pressure and the most critical plastic zone are both smaller when the effective angle of internal friction of the soil is larger, and they both increase when the hydraulic head at the ground surface increases, but they both decrease when the tunnel influx increases due to the fact that pore water pressure at the tunnel perimeter decreases with the tunnel influx.