sustainability

Article
The Face Stability Analysis of Shield Tunnels Subjected to
Seepage Based on the Variational Principle
Zhihui Zhou, Tonghua Ling, Fu Huang * and Min Zhang

School of Civil Engineering, Changsha University of Science and Technology, Changsha 410004, China
* Correspondence: [email protected]; Tel.: +86-139-7514-6439

Abstract: The stability of tunnel face remains a great challenge for tunnel engineers, especially when
excavation is in a complex stratum under the water table. This paper aims to investigate the stability
of soil in the front of a shield tunnel face induced by the presence of insufficient cabin pressure under
the water table, on the basis of a variational principle and an upper bound theorem. The analytical
expression of the rupture surface for the tunnel face is obtained, and the shapes of the rupture surface
are plotted. Comparisons are made to check the present approach against the solution provided by
numerical simulation techniques in order to show that the proposed method is valid. The parameter
analysis indicates that the groundwater seepage has a significant effect on the range of the rupture
surface for the shield tunnel face.

Keywords: collapse surface; seepage; variational principle; cabin pressure; Hoek–Brown failure criterion
Citation: Zhou, Z.; Ling, T.; Huang,
F.; Zhang, M. The Face Stability
Analysis of Shield Tunnels Subjected
to Seepage Based on the Variational 1. Introduction
Principle. Sustainability 2022, 14,
The shield tunneling method has been developed as a main construction method
16538. https://doi.org/10.3390/
su142416538
for subway construction owing to its safety for excavation and its low disturbance to the
surrounding environment. It is well known that stable soil in the front of the tunnel face,
Academic Editors: Zdravko Trivic, which is determined by the cabin pressure, is a key element in the constructive safety of a
Vesna Kosorić, Siu-Kit Lau, shield tunnel. The existing literature shows that insufficient cabin pressure will induce the
Abel Tablada, Miljana Horvat,
collapse of the tunnel face, whereas overlarge cabin pressure will lead to soil uplift in the
Milena Vukmirović,
front of the tunnel face. When a shield tunnel is being drilled in a water-rich stratum, the
Silvia Domingo-Irigoyen,
underground water table will be disturbed by the excavation, which induces underground
Marija Todorović, Jérôme H. Kaempf,
water seepage. However, numerous articles in the literature show that underground water
Kosa Golić and Ana Peric
seepage exerts a greatly adverse effect on the face stability. Moreover, the engineering
Received: 24 August 2022 accidents induced by instability of the tunnel face may result in serious casualties and
Accepted: 23 November 2022 property damage. Therefore, studying the stability of the tunnel face in shield tunnel
Published: 9 December 2022 excavation under the groundwater table is of primary significance.
Publisher’s Note: MDPI stays neutral
Currently, numerous scholars have employed various methods to investigate the
with regard to jurisdictional claims in stability of the tunnel face. Because a limit analysis theorem can be used to effectively
published maps and institutional affil- investigate the stability of the geotechnical engineering, many scholars have used this
iations. method to study this problem. Mollon et al. [1] established a three-dimensional (3D)
multiblock rupture mechanism to study the collapse and the blow-out failure modes of a
tunnel face by applying the upper bound method of limit analysis. Later, Mollon et al. [2]
used a spatial discretization technique to establish a 3D rupture mechanism, and they
Copyright: © 2022 by the authors. derived a corresponding theoretical solution of the supporting force on the basis of the
Licensee MDPI, Basel, Switzerland. upper bound approach. By modifying the rupture mechanism proposed by Mollon et al. [2],
This article is an open access article Senent et al. [3] used this improved mechanism in conjunction with the Hoek–Brown yield
distributed under the terms and criterion to study the stability of shield tunnels excavated in cracked rock. Moreover,
conditions of the Creative Commons
Yang et al. [4] studied the influence of the inhomogeneity and anisotropy of soil on the
Attribution (CC BY) license (https://
ultimate supporting force of a shield tunnel that was excavated in an anisotropic stratum.
creativecommons.org/licenses/by/
Consequently, Xu et al. [5] proposed a new rotational failure mechanism and used this
4.0/).

Sustainability 2022, 14, 16538. https://doi.org/10.3390/su142416538 https://www.mdpi.com/journal/sustainability

Sustainability 2022, 14, 16538 2 of 13

mechanism in conjunction with a nonlinear failure criterion to compute the theoretical

solution of the cabin pressure for tunnel face. Later, Zhang et al. [6] constructed three-
dimensional active and passive failure mechanisms of the tunnel face by investigating the
failure features of rock mass in front of the tunnel face. On the basis of these mechanisms,
Zhang et al. [6] obtained the safe range of supporting pressures of the shield tunnel face
and validated the effectiveness of the computing results by applying this method in an
actual project. To study the safety of a tunnel face when considering tensile strength cutoff,
Chen et al. [7] proposed two improved failure mechanisms in the tunnel face. Using these
new failure mechanisms, the limit supporting pressures of the tunnel face for collapse and
blow-out modes are calculated and the calculating results are compared with the existing
solutions.
Groundwater seepage is a one factor that affects tunnel stability when a tunnel is
excavated under the water table. Thus, the study of the influence of seepage on the safety
of a shield tunnel has drawn a great deal of attention, for some scholars. By improving on
an existing model, Perazzelli et al. [8] employed the limit equilibrium theory to investigate
the safety of the tunnel face under seepage flow conditions. To investigate the influence of
pore pressure on the stability of soil in front of a tunnel face, Pan and Dias [9] calculated the
distribution of pore pressure and obtained the theoretical solution for the cabin pressure
of the tunnel face. Subsequently, using a modified three-dimensional failure mechanism
as a base, Pan and Dias [10] studied the stability of soil in front of a tunnel face under
the water table by employing an upper bound theorem of limit analysis. Later, finding
that the study of the coupled flow deformation of a tunnel face is rare, Zou and Qian [11]
developed a theoretical method to calculate the supporting pressure of a tunnel face for a
shield tunnel drilled below the water table. Using the discretization technique as a base,
Yang and Zhong [12] proposed a new active failure mechanism to investigate the stability
of soil in front of a tunnel face when the tunnel is excavated below the water table. Using
this mechanism, Yang and Zhong [12] calculated the supporting pressure of the tunnel
while accounting for the pore pressure and validated the effectiveness of their method
by comparing their results with the existing solutions. Finding that most of the existing
studies have focused on the stability of tunnel faces excavated in saturated soil, Li et al. [13]
proposed a method to study the stability of a tunnel face by accounting for the unsaturated
seepage effect. Using this method, Li et al. [13] calculated the safety factor of the shield
tunnel face and assessed the stability of a tunnel face drilled in saturated stratum. Using
the computational fluid dynamics-discrete element method, Fu et al. [14] investigated the
seepage effect on the rupture mode of the underwater tunnel face. Their achievements show
that under seepage situations, the rupture region in the front of the tunnel face increases as
the water depth increases. Weng et al. [15] designed a centrifuge test to study the rupture
mode of tunnel faces caused by the lengthways slope angle of the tunnel and steady-state
seepage. Their model test results present the progressive rupture mode of the excavation
face caused by the two factors mentioned above.
The purpose of our study is to investigate the groundwater seepage effect on shield
tunnel face stability. To achieve this target, a new rupture mechanism that can be employed
to describe the failure feature of tunnel face induced by insufficient cabin pressure under
the water table is constructed. By regarding seepage force as an external force, the external
rate of work produced by seepage force is introduced into the virtual equation, and the
objective function of the rupture surface for the tunnel face is obtained. In the context of
a variational principle, the equation of the rupture surface for the tunnel face at the limit
state is derived, and likewise, the shapes of the rupture surface for various parameters are
drawn. The proposed method is validated by comparing it with the results that have been
provided by numerical simulation. The method can obtain the collapse region of the soil in
the front of the tunnel face, which provides a useful reference for a stability assessment of a
shield tunnel excavated under the water table.
Sustainability 2022, 14, 16538 3 of 13

2. Construction of a Failure Mechanism for a Tunnel Face, Accounting for

Groundwater Seepage
During the shield tunnel construction process, the pressure of the chamber of the
shield machine is a key element in maintaining the stability of the tunnel face. The soil in
front of the tunnel face may collapse when the soil chamber pressure is deficient to resist
the earth pressure produced by the soil in front of the tunnel face. To describe the collapse
characteristic of the tunnel face, an upper bound failure mechanism, which is composed of
an arbitrary curve z(y), is constructed, as shown in Figures 1 and 2. Curve z(y) extends
from the bottom to the top of the tunnel face, and a closed collapse surface forms in front of
the tunnel face. The soil within the region of curve z(y) may collapse in the tunnel because
of insufficient soil chamber pressure. Furthermore, D is the diameter of the shield tunnel,
and L is the depth of the tunnel. H is the space from the groundwater head to the shield
tunnel roof, and v is the velocity vector of the collapse block. It is assumed that the pressure
of the chamber provided by the shield machine is σT , which is evenly distributed over
the entire excavation face. For calculation simplicity, some assumptions about the failure
mechanism are presented, as follows:
Sustainability 2021, 13, x FOR PEER REVIEW
1. The stratum around the tunnel is homogeneous soil, and the stratum is regarded as
an ideal elastic-plastic body;
2. The collapse block is a rigid body, which means the deformation of the collapse block
can be ignored.

Water table elevation

y
H

Failure surface
D/2

v
Direction of tunneling Failure block
Sustainability 2021, 13, x FOR PEER REVIEW
z
D/2

Shield tunnel face

Figure 1. Diagram of soil collapse damage in front of excavation.

Figure 1. Diagram of soil collapse damage in front of excavation.

y
Water table elevation

y
H e
n α
L z Z=Z (y)
v

Z'
α
Failure block
σT β
1

D v

Z
σT

Shield tunnel face

Figure 2.
Figure 2. Mechanism
Mechanism ofof
soil collapse
soil damage
collapse in front
damage of excavation
in front in groundwater
of excavation seepage.seepage.
in groundwater
Sustainability 2022, 14, 16538 4 of 13

3. Basic Principles
3.1. Upper Bound Theorem
The upper bound theory states that for any kinematically admissible velocity field, the
internal energy dissipation power along the failure surface is not less than the power of the
external loads. The equation of the virtual work equation can be expressed as follows:
Z Z Z
V σij ε ij dV ≥ S Ti vi ds + V Xi vi dV (1)

where σij and ε ij are the stress and strain rate in the upper bound velocity field, respectively;
Ti is an overload on the border; Xi is the volume force; and vi is the velocity along the
velocity discontinuity surface.

3.2. Hoek–Brown Yield Criteria

To study the failure features of the rock mass, the Hoek–Brown yield criteria is applied,
which can be written as:  a
σ3
σ1 = σ3 + σci mb +s (2)
σci
where σ1 , σ3 , and σci are maximum main stress, minimum main stress, and the uniaxial
compressive strength of the rock mass, respectively; mb , s, and a are material parameters
that depend on geological indicators (GIs). The formula takes the following form:

GSI − 100
 
mb = mi exp (3)
28 − 14D

GSI − 100
 
s = exp (4)
9 − 3D
!
1 e−GSI/15 − e−20/3
a= + (5)
2 6

where mi is a constant related to geotechnical properties and D is the disturbance factor.

Because the internal dissipation rate of energy on the rupture surface is calculated by the
normal and tangential stress strain, the Hoek–Brown yield criteria expressed by normal
and shear stresses is used in this study, which can be written as:
B
σn − σt

τ = Aσc = E(σn − σt ) B (6)
σc

where σn is the normal stress; τ is the shear stress; A and B are material constants; and σc
and σt represent the uniaxial compressive strength and the tensile strength of the rock mass,
respectively.

3.3. The Analysis of Seepage Field

The influence of groundwater on the stability of geotechnical structures represents the
interaction of pore water between soil particles. Moreover, the formation of the pore water
pressure is related to the permeability of the soil. Therefore, when studying the stability
of geotechnical structures in the seepage field, the effect of seepage must be taken into
account. Because the seepage force is a volume force that acts on every point of the soil in
the seepage field, calculating the seepage force is difficult. In the two-dimensional seepage
model, the drag force of seepage water on soil particles in a unit area can be calculated by
using the following formula:

∂h ∂h
f z = − γw , f y = − γw (7)
∂z ∂y
Sustainability 2022, 14, 16538 5 of 13

where h is the total hydraulic head and γw is the water unit weight. By applying the Gauss
integral theory to the surface integral of seepage force along the rupture block in front of
the tunnel face, the expressions for the horizontal and vertical components of the seepage
force acting along the rupture block can be obtained:
 Z Z 
Fz ( β) = γw − sin β h∗ ds + h∗ ds (8)
S1 S2
 Z Z 
Fy ( β) = γw cos β h∗ ds + h∗ ds (9)
S1 S2

where S1 is the sliding surface, S2 is the tunnel face, and β is the angle between the sliding
block and the horizontal plane. h∗ is the average water head on the horizontal direction of
the sliding block, which can be calculated by using the following formula:

1
Z
h∗ = h(z, y)dy (10)
D
By combining Equations (8)–(10), the horizontal and vertical components of the seep-
age force can be expressed as:

Fz ( β) = γw − sinh β S1 h∗ds + S2h∗ ds

R R

Ry √  i (11)
= γw −0y − sin β H + D2 1 + z02 + H + D2 dy
0

Fy ( β) = γw coshβ S1 h∗ ds − S2 h∗ ds
R R

Ry  √  i (12)
= γw −0y cos β H + D2
0
1 + z 02 + H + D
2 dy

where D is the tunnel excavation diameter and H is the space between the groundwater
level and the tunnel roof.

4. Upper Bound Solution for the Collapse Surface for a Tunnel Face, Accounting for
Groundwater Seepage
As Chen [16] pointed out, the equation of virtual work in the kinematically admissible
velocity field comprises the external rate of work and the internal rate of dissipation.
Moreover, the external rate of work is produced by the external force of the presented
failure mechanism, which includes the soil gravity, the pressure of the chamber, and the
seepage force. Thus, the external rate of work generated by the external force can be
calculated. The rate of work done by soil gravity can be obtained by applying it to the
entire failure surface. Z y0
·
WG = γzv sin βdy (13)
− y0

where γ is the bulk density of rock/soil. Because the pressure of chamber σT is evenly
distributed over the entire tunnel face, the rate of work produced by the pressure of the
chamber can be expressed by:
· Z y0
WC = −σT cos βvdy (14)
− y0

As mentioned above, the effect of seepage is a significant factor that should be taken
into account in tunnel face stability analysis. Thus, using the expressions of horizontal and
vertical components of the seepage force, an expression for the rate of work produced by
seepage force can be written as:
· Z y0 Z y0  

D
WW = Fz cos β + Fy sin β vdy = γw H + (sin β − cos β)dy (15)
− y0 − y0 2
Sustainability 2022, 14, 16538 6 of 13

On the other hand, the computation of the dissipation rate in the presented rupture
mechanism is more complex than the calculation of the external rate of work. The soil is
assumed to obey an associated flow law. When the potential surface is in accordance with
the yield surface, the potential function is:

G = τ − E(σn − σt ) B (16)

Moreover, when the plastic flow occurs in the collapse surface, the strain increment
can be expressed as:  ·
 ε n = λ ∂G = −λEB(σn − σt ) B−1
 i ∂σn
(17)
 γ·n = ∂G = λ
∂σn
i
where λ is the plasticity multiplier. As mentioned above, any random point on the collapse
surface will obey the principle of small deformation, and the rate of deformation for
the point is approximately equal to the rate of strain. Given the geometric relationships
illustrated in Figure 3, the strain increment of any point on the collapse surface is:
 ·
 ε n = v sin(αi − β)
t
 i (18)
Sustainability 2021, 13, x FOR PEER REVIEW  γ·n = − v cos(αi − β)
i t

where αi is the angle between the tangent vector and the z-axis at any point on the collapse
surface. Moreover, tan αi = z0 (y), where t is the thickness of the plastic flow zone.

Figure 3. Collapse
Figure 3. Collapsesurface
surfaceof of a tunnel
a tunnel faceface for different
for different of D. of D .
valuesvalues

By combining Equations (16) and (17), the following equation can be obtained:
v
λ = − cos(αi − β) (19)
t
B
(
τi = E[ EB cot(αi − β)] 1−B
1 (20)
σi = σt + [ EB cot(αi − β)] 1−B
Because the energy dissipation occurs only at the velocity discontinuity surfaces,
the rate of energy dissipation for any point on the velocity discontinuity surfaces can be
calculated only by superimposing the rate of energy dissipation on the normal direction
and the tangential direction. Thus, the rate of energy dissipation for any point on the
rupture surface is:
Sustainability 2022, 14, 16538 7 of 13

· · ·
∆ Di = σi ε i + τi γn
i
v sin(αi − β) 1 1 B
o (21)
= t σt + cot(αi − β) 1−B [( EB) 1−B − E( EB) 1−B ]

For mathematical simplicity, the following equation is introduced:

1 B
E = ( EB) 1−B − E( EB) 1−B (22)

Merging part of Equation (21) with part of Equation (20), the rate of energy dissipation
for any point on the collapse surface is:
· v sin(αi − β) h 1
i
∆ Di = σt + E cot(αi − β) 1−B (23)
t
On the basis of the geometrical relation illustrated in Figure 2 and trigonometric
function, the following equations are determined.

z0 cos β − sin β 1 + z0 tan β

sin(αi − β) = , cot(αi − β) = 0 (24)
z − tan β
p
1 + z0 2
By applying Equation (22) to the entire collapse surface and combining Equation (23),
the total rate of energy dissipation along the collapse surface is:
1
v(z0 cos β−sin β) 1+z0 tan β
R y0  
1−B
WD = − y0
√ [− σt + E ]ds z0 − tanβ
d 1+ z 0 2
Ry   1 
0 tan β 1−B
(25)
= −0y0 (z0 cos β − sin β)[− σt + E 1z+0 −z tanβ vdy

On the basis of the kinematical approach, the relationship between the rate of external
work and the rate of energy dissipation is obtained:
· · ·
WG + WC + WW ≤ WD (26)

The aim of this study is to investigate the failure mode and failure range of the soil in
the front of the tunnel face at limit state. To obtain the analytical equations of the rupture
surface at the limit state, it is necessary to establish a target function that includes the
equation of the failure surface, by using the external rate of work and the rate of the energy
dissipation. By incorporating Equations (12)–(14) and (24) into Equation (25), the target
function is determined:
· · ·
 
J = WD − WG + WC + WW
 1
0 tan β 1−B (27)
Ry 
= −0y0 (z0 cos β − sin β)[− σt + E 1z+0 −z tanβ ]−γz sin β + σT cos β
 
−γw H + D2 (sin β − cos β)vdy

Because independent variable z0 in Equation (26) is a function, objective function J can

be regarded as a functional, which can be expressed as:
Z y0
F y, z(y), z0 (y) dy
 
J= (28)
− y0

where
 1 
1+z0 tan β
 
1− B
F [y, z, z0 ] = (z0 cos β − sin β) σt + E − γz sin β + σT cos β
z0 −tan β
  (29)
−γw H + D2 (sin β − cos β)
Sustainability 2022, 14, 16538 8 of 13

However, the upper bound solution, which derives from the relation between the
external rate of work and the rate of the energy dissipation, is not the real result. According
to Chen [16], the real result of the rupture surface can be obtained only when objective
function J has reached its extreme value. Furthermore, Equation (28) is a special case of the
Euler equation, which means that the extreme solution of the functional can be converted
into solving the solution of Equation (28) under fixed boundary conditions. Thus, the first
incorporation of Equation (28) is:

∂F
F − z0 = c1 (30)
∂z0
where
 1 B
1 + z0 tan β 1−B (1 + z0 tan β) 1−B

∂F E
= σt cos β + E cos β − (31)
∂z0 z0 − tan β (1 − B) cos β (z0 − tan β) 1−1 B

After some simplification, the following equation is derived:

1
1+z0 tan β
  h 0
i
1− B
E − sin β + (1− B)(1+zz0 tan β) cos β − σt sin β − γz sin β
z0 −tan β
  (32)
+σT cos β − γw H + D2 (sin β − cos β) = c1

On the basis of Equation (31), the expression of z0 is obtained. After calculating the
incorporated expressions, a set of curves that reaches extreme solutions is derived:

dz
Z
y= + c2 (33)
ϕ(z, c1 )
 
where c1 and c2 are integration constants. Using geometric boundary conditions z D2 =
 
z − D2 = 0, the values of c1 and c2 can be determined.
The value of parameter B is extremely important for solving Equation (31). According
to Hoek and Brown [17], the value of parameter B ranges from 0 to 1. However, Equa-
tion (31) is a linear differential equation only when B = 0.5, which can be solved analytically.
When B 6= 0.5, Equation (31) is a complex nonlinear partial differential equation, which
cannot be solved analytically. Therefore, in this study, the equation for a collapse surface
of a tunnel face at limit state is discussed only for the parameter B = 0.5. When parameter
B = 0.5, Equation (31) can be simplified as:

1+z0 tan β 2
  h i
2z0
E 0
z −tan β − sin β + − σt sin β − γz sin β
 (1+z0 tan β) cos β (34)
D
+σT cos β − γw H + 2 (sin β − cos β) = c1

A MATLAB [18] program is employed to solve this nonlinear partial differential

equation. Firstly, a form of Equation (33), which is expressed by separating variables z0 and
y, can be obtained by using a simplified calculation. Then, by substituting fixed boundary
conditions, the values of integration constants c1 and c2 are determined.
n Finally, a particular
o
solution among these solutions that satisfies boundary condition − D2 ≤ y ≤ D2 , z > 0 is
found. In summary, the equation for the collapse surface of the soil in front of the shield
tunnel face is obtained:
√
z(y) = {1 6E − 16 cos 2β + 16 2 sin2 β E [2 E + 2E cos4 β + cos2 β − cos4 β


(2σt + Dγw + 2Hγw ) + 2c1 cos2 β sin β + cos3 β sin β(2c2 γ − 2σT − Dγw
1
−2Hγw − 2γy) + cos β sin β(2 γy − 2c2 γ)]} 2 + (c2 γ − γy) (35)
(3 cos 3β sin β + sin 3β cos β − cos 3β sin 3β − cos β sin β )}/[γ(cos 3β
− cos β)2 ]
Sustainability 2022, 14, 16538 9 of 13

5. Effect of Various Parameters on the Shape of a Collapse Surface of a Tunnel Face

To study the effect of different parameters on the collapse range of the soil in the
front of the shield tunnel face, while accounting for groundwater seepage, the shapes
of the collapse surface for the soil in front of the shield tunnel face are plotted based on
Equation (34), with various parameters. These parameters, which are used in the parametric
Sustainability 2021, 13, x FOR PEER REVIEW
analysis, are presented as follows: γ = 7.5 kN/m3 , L = 10 m, H = 10 m, σT = 18 kN, σt = 0,
σc = 10~100 kPa, A = 0.4~0.7, and β = 33~66◦ . Moreover, the shapes of the collapse surface
Sustainability 2021, 13, x FOR PEER REVIEW
for the soil in the front of the shield tunnel face that correspond to the abovementioned
parameters are plotted in Figures 3–7.

Figure 4. Collapse surface of a tunnel face for different values of A .

Figure
Figure 4.
4. Collapse
Collapsesurface
surfaceofofa atunnel
tunnelface forfor
face different values
different A .A.
of of
values

Sustainability 2021, 13, x FOR PEER REVIEW

Figure
Figure5.5.Collapse
Collapsesurface
surfaceofofa atunnel face
tunnel forfor
face different values
different of 
values ofc σ. c .

Figure 5. Collapse surface of a tunnel face for different values of  c .

Figure 6.
Figure 6.Collapse
Collapsesurface of aoftunnel
surface faceface
a tunnel for different of  . of β.
valuesvalues
for different
 Figure 6. Collapse surface of a tunnel face for different values of  .
Sustainability 2022, 14, 16538 10 of 13

Figure
Figure7.7.Effect
Effectofofgroundwater
groundwaterseepage onon
seepage thethe
collapse surface
collapse of aof
surface tunnel face.face.
a tunnel

According to Figures 3–7, the ranges of the collapse surface increase with the increase
of parameters D, A, and β. Furthermore, the ranges of the collapse surface decrease with
the increase of parameter σc when other parameters are fixed.
To investigate the effect of underground water seepage on the collapse region of the soil
in the front of the shield tunnel face, the shapes of the rupture surface while accounting for
seepage and while not accounting for seepage are drawn in Figure 7. According to Figure 7,
the collapse range under the consideration of seepage effect is about 0.2 m larger than
the collapse range that doesn’t account for seepage when other parameters are the same.
Consequently, the effect of groundwater seepage on the stability of the soil in front of shield
tunnel face cannot be ignored when a shield tunnel is excavated in a groundwater-rich
stratum.

6. Comparison with Numerical Simulation Result

6.1. Numerical Model of Shield Tunnel That Accounts for Groundwater Seepage
To verify the validity of the theoretical result presented in this study, the theoretical
result is compared with the numerical solution provided by the numerical simulation
technique. A three-dimensional numerical model of a shield tunnel excavated under the
groundwater table is constructed with software FLAC3D [19]. The diameter of the tunnel
D is 6 m, the distance between the ground surface and the tunnel roof is 9 m, and the
groundwater table measured from the tunnel roof is also 9 m. The size of the 3D numerical
model is 100 m, 50 m, and 69.4 m in the transversal, longitudinal, and vertical directions,
respectively. The thickness of the tunnel lining is 0.3 m. Because the accuracy of the meshing
grid has a significant influence on the numerical calculation results, an inhomogeneous
grid is employed in the constructed model, which can be seen in Figure 8. To make the
Sustainability 2021, 13, x FOR PEER REVIEW
comparison under the same conditions, the rock mass is assumed as homogeneous.

Figure 8.
Figure 8. 3D
3D model
modelof
ofshield
shieldtunnel
tunnelexcavated
excavatedunder
underthethe
water table.
water table.
Sustainability 2022, 14, 16538 11 of 13

Because the rate of the energy dissipation of the collapse surface is computed by
superimposing the rate of energy dissipation on the normal direction and the tangential
direction, the H–B yield criteria expressed by normal and shear stresses is used in this study.
However, the H–B yield criteria embedded in the FLAC3D [19] program is represented
by the major and minor main stresses. Thus, to make the comparison under the same
circumstances, the parameters of H–B criteria in the two forms should be equivalently
converted. Using the method proposed by Hoek and Brown [16], the parameters of the
H–B failure criterion listed in Table 1 are equivalently converted into the parameters of the
H–B failure criterion listed in Table 2.

Table 1. Hoek–Brown parameters represented by normal and shear stresses.

Elastic Poisson’s Hoek–Brown Hoek–Brown

Weight Uniaxial Compressive Uniaxial Tensile
Modulus Ratio Parameter Parameter
kN/m3 Strength σci /kPa Strength σtm /kPa
G/MPa µ A B
18 80 0.3 150 0.5 0.2505 0.5

Sustainability 2021, 13, x FOR PEER REVIEW

Table 2. Hoek–Brown parameters represented by major and minor principal stresses.

Elastic Poisson’s Uniaxial Compressive Hoek–Brown Hoek–Brown Hoek–Brown

Weight
Modulus Ratio Strength Parameter Parameter Parameter
kN/m3
G/MPa µ σci /kPa a mb s
18 80 0.3 150 0.5 0.3155 0.0622

6.2. Comparison Analysis between Numerical Solution and Upper-Limit Solution

Using the parameters listed in Table 2 and the numerical model mentioned above,
the contours of the displacement for the tunnel face that reflect the shape of a collapse
surface are obtained, which can be seen in Figure 9. The collapse surface constituted by
the boundary of the displacement is in accordance with the rupture surface derived from
the theoretical analysis. Furthermore, the upper bound solution of the collapse surface
obtained from Equation (34) and Table 1 is superimposed in Figure 9. According to this
figure, the3Dtheoretical
Figure 8. resulttunnel
model of shield of the excavated
collapse surface is similar
under the to the collapse surface obtained
water table.
from the numerical simulation. The similarity between the two methods for the tunnel
collapse surface indicates that the theoretical methods proposed in this study are valid.

Figure 9.
Figure 9. Comparison
Comparison of
of aa collapse
collapse surface
surface between
between aa theoretical
theoreticalsolution
solutionand
andaanumerical
numericalsolution.
solution.
Sustainability 2022, 14, 16538 12 of 13

7. Conclusions
This paper aimed to study the rupture mode of a shield tunnel caused by insufficient
cabin pressure under the water table. A new collapse mechanism, which is composed
of an arbitrary curve, was constructed to illustrate the collapse characteristic of the soil
in the front of the tunnel face. Using this mechanism, an analytical expression of the
collapse surface of a shield tunnel face was obtained, according to a kinematical approach
while accounting for the effect of seepage. By comparing it with the numerical solution of
collapse failure surface obtained from the numerical simulation technique, the validity of
the proposed theoretical method was verified. Using the equation of the collapse surface for
a shield tunnel face, the shapes of the rupture surface were drawn for various parameters.
The parameter analysis indicated that the groundwater seepage had a significant effect
on the range of the rupture surface for the shield tunnel face. The range of the collapse
surface increased when the groundwater seepage was factored in. Moreover, the region of
the rupture surface increased with the increases in parameters D, A, and β and decreased
with the increase in σc . If a shield tunnel is excavated in groundwater-rich strata, the effect
of groundwater seepage on the stability of the tunnel face cannot be ignored.

Author Contributions: Conceptualization, F.H.; Formal analysis, Z.Z. and F.H.; Data curation, M.Z.;
Writing—original draft, Z.Z.; Project administration, Z.Z.; Funding acquisition, T.L. All authors have
read and agreed to the published version of the manuscript.
Funding: This study was supported by the National Natural Science Foundation of China (Grants
51878074, 52278395 and 52078061) and the Natural Science Foundation of Hunan Province, China
(Grant No. 2021JJ30714). Their financial support is greatly appreciated.
Data Availability Statement: The data presented in this study are available on request from the
corresponding author.
Conflicts of Interest: The authors declare no conflict of interest.

References
1. Mollon, G.; Dias, D.; Soubra, A.H. Probabilistic analysis and design of circular tunnels against face stability. Int. J. Geomech. 2009,
9, 237–249. [CrossRef]
2. Mollon, G.; Dias, D.; Soubra, A.H. Rotational failure mechanisms for the face stability analysis of tunnels driven by a pressurized
shield. Int. J. Numer. Anal. Methods Geomech. 2011, 35, 1363–1388. [CrossRef]
3. Senent, S.; Mollon, G.; Jimenez, R. Tunnel face stability in heavily fractured rock masses that follow the Hoek-Brown failure
criterion. Int. J. Rock Mech. Min. Sci. 2013, 60, 440–451. [CrossRef]
4. Yang, X.L.; Li, W.T.; Pan, Q.J. Influences of anisotropy and inhomogeneity on supporting pressure of tunnel face with kinematical
approach. J. Cent. South Univ. 2015, 22, 3536–3543. [CrossRef]
5. Xu, J.S.; Du, D.C.; Yang, Z.H. Upper bound analysis for deep tunnel face with joined failure mechanism of translation and rotation.
J. Cent. South Univ. 2015, 22, 4310–4317. [CrossRef]
6. Zhang, J.H.; Wang, W.J.; Zhang, D.B.; Zhang, B.; Meng, F. Safe range of retaining pressure for three-dimensional face of pressurized
tunnels based on limit analysis and reliability method. KSCE J. Civ. Eng. 2018, 22, 2625–2656. [CrossRef]
7. Chen, G.H.; Zou, J.F.; Liu, S.X. Stability Analysis of Pressurized 3D Tunnel Face with Tensile Strength Cutoff. Int. J. Geomech. 2021,
21, 04021226. [CrossRef]
8. Perazzelli, P.; Leone, T.; Anagnostou, G. Tunnel face stability under seepage flow conditions. Tunn. Undergr. Sp. Tech. 2014, 43,
459–469. [CrossRef]
9. Pan, Q.J.; Dias, D. The effect of pore water pressure on tunnel face stability. Int. J. Numer. Anal. Methods Geomech. 2016, 40,
2123–2136. [CrossRef]
10. Pan, Q.J.; Dias, D. Three-dimensional face stability of a tunnel in weak rock masses subjected to seepage forces. Tunn. Undergr. Sp.
Tech. 2018, 71, 555–566. [CrossRef]
11. Zou, J.F.; Qian, Z.H. Face-stability analysis of tunnels excavated below groundwater considering coupled flow deformation. Int. J.
Geomech. 2018, 18, 04018089. [CrossRef]
12. Yang, X.L.; Zhong, J.H. Stability analysis of tunnel face in nonlinear soil under seepage flow. KSCE J. Civ. Eng. 2019, 23, 4553–4563.
[CrossRef]
13. Li, Z.W.; Yang, X.L.; Li, T.Z. Face stability analysis of tunnels under steady unsaturated seepage conditions. Tunn. Undergr. Sp.
Tech. 2019, 93, 103095. [CrossRef]
14. Fu, Y.B.; Zeng, D.Q.; Xiong, H.; Li, X.H.; Chen, Y.L. Seepage effect on failure mechanisms of the underwater tunnel face
viaCFD–DEM coupling. Tunn. Undergr. Sp. Tech. 2022, 146, 104591.
Sustainability 2022, 14, 16538 13 of 13

15. Weng, X.L.; Sun, Y.F.; Yan, B.H.; Niu, H.S.; Lin, R.G.; Zhou, S.Q. Centrifuge testing and numerical modeling of tunnel face stability
considering longitudinal slope angle and steady state seepage in soft clay. Tunn. Undergr. Sp. Tech. 2020, 101, 103406. [CrossRef]
16. Chen, W.F. Limit Analysis and Soil Plasticity; Elsevier Science: Amsterdam, The Netherlands, 1975.
17. Hoek, E.; Brown, E.T. Practical estimate of rock mass strength. Int. J. Rock Mech. Min. Sci. 1997, 34, 1165–1186. [CrossRef]
18. MATLAB 9.1; Computer Software. MathWorks Inc.: Natick, MA, USA, 2016.
19. FLAC3D; Computer Software. Itasca Consulting Group: Minneapolis, MN, USA, 2012.