Upper-bound solutions for the face stability of a non-circular NATM tunnel in clays with a linearly increasing undrained shear strength with depth

Computers and Geotechnics 114 (2019) 103136

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Upper-bound solutions for the face stability of a non-circular NATM tunnel in clays with a linearly increasing undrained shear strength with depth

T

Wei Li, Chengping Zhang , Wenjun Zhu, Dingli Zhang ⁎

Key Laboratory of Urban Underground Engineering of the Education Ministry, Beijing Jiaotong University Beijing 100044, China School of Civil Engineering, Beijing Jiaotong University, Beijing 100044, China

ARTICLE INFO

ABSTRACT

Keywords: NATM tunnel Face stability Safety factor Limit analysis Continuous velocity field

The face stability of the NATM tunnels (the New Austrian Tunneling Method) in undrained clays with a linearly increasing strength is numerically and analytically investigated by defining two types of the safety factors F in this paper. Based on the kinematic approach of the limit analysis, a continuous velocity field for an open-face tunnel is proposed. Moreover, the calculated solutions are compared with the results obtained from the numerical simulation and other existing literatures to validate its effectivity and superiority. The stability charts for different influence factors are also provided. Finally, an improved velocity field is presented, which provides a closer result to the real situation of the face collapse.

1. Introduction Due to the rapid development of urban subways, the stability analysis of the tunneling is widely studied during the last decades. When dealing with the urban shield tunnels, it is crucial to make sure the support pressure uniformly distributed on the tunnel face. Otherwise, the shield tunnel face becomes unstable or even collapses, which might finally lead to the subsidence of the surface ground. As for the open-face excavation of the New Austrian Tunneling Method (NATM), no supporting pressure is applied on the tunnel face. A more appropriate and convenient evaluation criterion for the face stability is required for the NATM tunnel face. Thus the stability analysis of the tunnel face is essential for the practical NATM tunneling projects. The face stability has been investigated by various methods including the analytical analysis, the experimental methods and the numerical simulation. Many researchers adopted the upper bound theorem of the limit analysis to study the stability issue of the geotechnical engineering. Based on the normal flow rule, the stress-strain relationship in the upper bound method of the limit analysis is considered in the limit analysis (Chen [1]; Drucker et al. [2]). So both the rigid-block mechanism (the cone or the logarithmic spiral surface) and the continuous deformation mechanism can be applied in the limit analysis of geotechnical problems and provide a rigorous upper bound solution. The existing researches adopted the upper bound theorem of the limit analysis to study the stability issue of the tunnel face by proposing a translational multi-rigid block mechanism (Davis et al. [3]; Leca and

⁎

Dormieux [4]; Soubra [5]; Mollon et al. [6,7]; Zhang et al. [8]; Han et al [9,10])or rotational rigid block failure mechanism (Subrin and Wong [11]; Mollon et al. [12]) for the frictional soils and a continuous velocity field (Osman and Mair [13]; Mollon et al. [14]; Klar et al. [15]; Klar and Klein [16]; Zhang et al. [17,18]; Huang et al. [19]) for the purely cohesive soils. It is worth noting that the failure mechanism of the tunnel face within the frame of the continuous velocity field significantly improves the upper solution of face stability in clays, because the motion of rigid blocks was not involved by the face collapse in purely cohesive soils but rather a continuous deformation of the soil. The rigid motion and the discontinuity surface (shear band) are only suitable for the frictional soils. In order to visualize the collapse pattern, both the numerical simulation (Vermeer et al. [20]; Chen et al. [21]; Ukritchon et al. [22–24]) and the experimental tests (Broms and Bennermark [25]; Mair [26]; Schofield [27]; Chambon and Corté [28]; Takano et al. [29]; Idinger et al. [30]; Chen et al. [31]) were performed to investigate the stability problem. On one hand, the numerical simulation can study the stability issue effectively by building an advanced numerical model or combining the 3D finite element analysis with the upper bound or lower bound theorem of the limit analysis (Ukritchon et al. [32–35]; Huang et al. [36]). On the other hand, the experimental tests are performed to capture the failure characteristic of the tunnel face. In summary, for the frictional soils, there will be two obvious slip bands starting from the tunnel invert and crown, respectively and then propagating towards the ground surface. Whereas for the purely cohesive soils, there is no

Corresponding author at: Key Laboratory of Urban Underground Engineering of the Education Ministry, Beijing Jiaotong University Beijing 100044, China. E-mail address: [email protected] (C. Zhang).

https://doi.org/10.1016/j.compgeo.2019.103136 Received 10 December 2018; Received in revised form 4 April 2019; Accepted 12 June 2019 Available online 15 June 2019 0266-352X/ © 2019 Elsevier Ltd. All rights reserved.

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Nomenclature

distance between the origin O and the point E Rf Rmax(θ, β)distance between the point E and the boundary of the cross-section Πβ L1 position of the maximal tangential velocity L1, r2, r3 distances of three centres O1, O2, O3 to the point E vβ axial component of velocity vr radial component of velocity vθ orthoradial component of velocity vm(β) maximum value of the axial velocity strain rate tensor W rate of work of the soil mass Ws rate of work of the surcharge WD rate of the energy dissipation F safety factor safety factor of unit area F* Nγ, Ns, ND non-dimensional coefficients of soil unit weight, surface surcharge and cohesion

R1, R2, R3 three radii of the horseshoe-shaped tunnel face O1, O2, O3 three centres of the horseshoe-shaped tunnel face α1, α2, α3 ranges of three sections of circular arc W equivalent width of the tunnel face D equivalent height of the tunnel face C/D buried depth to diameter ratio γ soil unit weight cu undrained strength of the soils E Young’s modulus v Poisson’s ratio ρ variation gradient of soil undrained strength σs surface surcharge Πβ cross-section of the failure mechanism Ri distance between the tunnel crown and the maximal tangential velocity point E localized shear band and the soil mass moves more like a ‘flow’ rather than a rigid block, as shown in Fig. 1. Most of the research on the face stability is focused on the shield tunneling and determining the face supporting pressures, but the attention to the stability of the open-face excavation for an NATM tunnel (Khezri et al. [37]; Pan and Dias [38]) is relatively limited. Moreover, An NATM tunnel usually has an arbitrary cross-section. So it will be very interesting to investigate the influence of the cross-section shape on the stability of the tunnel face. The existing studies (Khezri et al. [37]; Pan and Dias [38]) on the NATM tunnels involved the translational or rotational rigid movement, which was only suitable for the frictional soils and couldn’t be applied for the continuous deformation of the clays. This paper will build a continuous velocity field for a noncircular NATM tunnel in clays with a linearly increasing undrained strength with depth, which can take into account of an arbitrary face shape and the full face failure. This proposed continuous velocity field presents a more realistic failure mechanism for the NATM tunnels, which describes the continuous deformation of the soil mass between the tunnel face and the ground surface more previsely compared with the rigid-block failure mechanism (Pan and Dias [38]). With the idea of open-face excavation, this paper aims at adopting the upper bound method of the limit analysis to investigate the stability of NATM tunnels with arbitrary cross-sections in purely cohesive soils (φ = 0). A continuous velocity field will be defined in the 3D failure mechanism of the NATM tunnels. The safety factor is adopted to assess the stability of the tunnel face of an NATM tunnel. A numerical simulation is carried out at first without any priori assumption regarding the shape of the failure mechanisms. The numerical results will be helpful for the definition of the continuous velocity field of an NATM tunnel with an arbitrary cross section. Based on the results of numerical simulations, a 3D failure mechanism combined with a continuous velocity field is proposed for the NATM tunnels in the framework of the kinematic approach of the limit analysis. And then the upper bound solutions obtained from the proposed continuous velocity field in this paper will be compared with the results of numerical simulations and other existing studies to validate its reasonability and superiority. Moreover, the influences of both the geometric and geologic parameters on the face stability will be investigated, and the charts of the dimensionless coefficients Nγ and ND versus C/D are provided. Finally, an improved velocity field of the NATM tunnels is proposed to provide a better upper bound solution especially for the deep buried tunnels.

a forward scheme for nonlinear problem, which does not require iteration and could deal with the large deformation problem well. Thus, FLAC3D is adopted in this paper to investigate the stability of the tunnel face and outline the velocity field of the failure mechanism. Various shapes of the tunnel face are adopted in the construction of the NATM tunnels, such as the circular, multi-circular and non-circular (including elliptical, rectangular and horseshoe-shaped cross-sections) tunnel faces. The most widespread ones, however, are multi-circular and oblate horseshoe-shaped tunnels. This is because the horseshoeshaped NATM tunnels can significantly improve the space utilization and have a better stress state for the tunnel face compared with the common circular or rectangular tunnels. There are many kinds of horseshoe-shaped tunnel profiles including different numbers of the component circular sections and different shapes of the tunnel invert and side wall. In general, the most economical section and other technical aspects of the NATN tunnels are not fixed, which mainly depend on the specific conditions of the practical engineering project. The choice of the tunnel profiles lies in the compromise among the operating requirements of the tunnel, minimizing stresses in the lining and the costs for the excavation and supporting. For the sake of a generalized explanation, a typical geometric shape of the urban tunnel design is adopted herein, as is shown in Fig. 2. More details on the the tunnel design can be referred to the practical cases of the urban tunnel construction [39–42]. The following relations refer to the geometrical properties of the horse-shaped tunnel. With the initial geometric parameters R1, R2, R3 and α1, α2, the following parameters are obtained:

2. Numerical simulations with FLAC3D The 3D finite difference code FLAC3D uses explicit Lagrangian calculation scheme and a mixed discretization zoning technique. This is

Fig. 1. Continuous deformation of soils (Schofield [27]). 2

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 2. Typical geometric shape of the urban tunnel design.

d1 =

R1) 2 + R22

(R2

2 cos

2· R2 (R2

enough to avoid the influence of the boundary effects. To study a purely cohesive soil, an elastic perfectly plastic constitutive model based on Tresca failure criterion (which is identical to a Mohr–Coulomb criterion with φ = 0) is assigned to soils, which means that a purely clayey soil is in an undrained condition and the soil plastic deformation takes place without any volume change. The Young’s modulus E and Poisson’s ratio ν of the clay are taken as 500cu and 0.495, respectively (Huang et al. [19]; Ukritchon et al [24]). It is worth noting that a higher value of the E/cu ratio is selected because it can increase the computation speed with no significant impact on the stability of the tunnel face. Moreover, the tunnel lining is modelled by the shell structural element with the Yong’s modulus of 33.5 GPa, Poisson’s ratio of 0.2 and the thickness of 0.35 m. The full face excavation are performed in the FLAC3D. The excavation process is simulated using a simplified single-step excavation scheme that assumes that part of the tunnel (10 m in length) is excavated instantaneously. And the lining shells are installed simultaneously. The safety factors provided by the FLAC3D are calculated by the automatic search combined with the strength reduction technique. Fig. 4 shows the velocity vector and contour of displacement of the NATM tunnels for C/D = 1.0 and 2.0. It can be seen that the failure mechanism involves a continuous plastic deformation of purely cohesive soils and the maximal velocity flow line locates at a distance below

R1)

d2 = R3 d1/sin 3·sin( 1 + 3 + ) 3 = arcsin[d1·sin( 1 + )/ R3 ] =

3 2

1

arccos{[(R2

Area =

2 1 R1

+ [ 3 R32

+[ (R3

R1 ) 2 + d12 2 2 R2

(R2

d2 ) R3·sin

R22]/[2d1 (R2 R1 ) R2 ·sin 3]

R1 )]}

2]

(1)

where α3, β, d1 and d2 are shown in Fig. 2. D and W are the equivalent height and width of a horseshoe-shaped NATM tunnel. Generally, the specific shape of the horseshoe-shaped tunnel can be determined with these geometrical parameters. And then the numerical model of an NATM tunnel is built. The detailed geometric parameters of the tunnel face are calculated and summarized in Table 1. These parameters of the horseshoe-shaped tunnel are not specific for any practical tunneling project. They are presented herein to provide a direct and generalized explanation for a typical urban tunnel. Since the face stability is closely related with the dimension of the tunnel face, the values of all the geometric parameters shown in Table 1 are adopted to obtain a horseshoe-shaped face with equivalent height D and width L of approximately 9.5 m and 10 m, respectively, which corresponds to an equivalent size of the circular face with a diameter D of 10 m. Thus the solutions obtained from this paper can be compared with the results of the existing researches [3,6,7,12,14] on the circular tunnel face with a diameter of 10 m. Considering the symmetry of the tunnel on the vertical plane, a 3D numerical model is built based on half of a tunnel cut lengthwise along the central axis. The nodes at the bottom of mesh are fully fixed for the boundary conditions. Nodes on the vertical symmetric plane are only fixed in the X direction but to assure the vertical movement. Fig. 3 presents a full face excavation numerical model of an NATM tunnel for C/D = 1.0 with a discretization of 58,410 zones (‘zone’ is the FLAC3D terminology for each discretized element) and a total of 62,186 nodes. The dimensions of the 3D numerical model are 20 × 30 × 50 m in the transversal, vertical, and longitudinal directions, which are large

Table 1 The geometric parameters of the tunnel face.

3

Geometric parameters

Values

R1 (m) R2 (m) R3 (m) α1 (°) α2 (°) α3 (°) β (°) d1 (m) d2 (m)

5.00 7.00 10.00 90.00 30.00 24.00 40.75 5.36 4.36

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

the center O1 of tunnel face. Zhang et al. [8] gave a simple and feasible criterion to outline the boundary trip of the failure zone at collapse in displacement contour, which suggested that the contour with a sudden increasing gradient could be defined as the boundary strip of the failure zone. However, this assessing criterion is proposed based on the rigidblocks failure mechanism. For the failure mechanism with a continuous deformation, the outermost layer of the displacement contour is more suitable to be considered as the failure zone. Moreover, it can be seen from Fig. 4 that the velocity field follows a similar parabola distribution in the tunnel face, and the velocity vector varies continuously in magnitudes and direction with the positions in the failure mechanism, which also indicates that the ground movement is a continuous deformation of soils. Next section will be devoted to the description of the continuous velocity field of an NATM tunnel based on these numerical results. The contours of the maximum shear strain increment for C/D = 1.0 and 2.0 obtained from FLAC3D are shown in Fig. 5. When C/D is equal to 1.0, the contour of the maximum shear strain plots the failure pattern of the tunnel face well. However, when C/D is greater than or equal to 2.0, it can be seen that the maximum shear strain mainly concentrates in front of the tunnel face. It is considered that the contour of the

Fig. 3. The numerical model of an NATM tunnel.

Fig. 4. Displacement contour and velocity vector of the NATM tunnels: (a) C/D = 1.0; (b) C/D = 2.0. 4

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Mollon et al. [14] had solved the problems by closed-form analytical formulas and explicit finite-difference discretization scheme, respectively. For the irregular profiles of the NATM tunnels, a rigorous analytical solution concerning the complexity of the horseshoe-shaped tunnel geometry and the velocity distribution is derived in this section. Based on the numerical results obtained from the FLAC3D in the previous section, this section will build a continuous velocity field in the face stability analysis of an NATM tunnel in the undrained clay and calculate the highest upper bound of the safety factor. Compared with the rigid-block failure mechanism of the frictional soils (Pan and Dias [38]), this continuous velocity field presents a more realistic failure mechanism for the non-circular NATM tunnels in clays, which represents the continuous movement of the clayey soils well and has no shear band within the failure mechanism. Besides, only the failure mechanism of the face collapse is considered in this paper because the situation of the face blowout is fairly rare in the practical engineering of the NATM tunnels. And a typical geometrical shape of the urban tunnel design is adopted herein for the engineering reference. 3.2. Failure mechanism 3.2.1. Geometry As shown in Fig. 6, a tunnel with equivalent height D and width L is excavated under a cover depth of C in a purely cohesive soil with a linearly increasing undrained strength of cu = cu0 + ρ(−y). A surcharge σs is applied on the ground surface and no pressure are supported on the tunnel face. The failure mechanism is generated by rotating the horseshoe-shaped tunnel face of increasing area (shaded zone in Fig. 6) about the rotation center O on the symmetric vertical plane y-z. The whole 3D failure region is defined by the orthogonal curvilinear coordinate (r, θ, β). In an arbitrary rotational plane Πβ, the radial distance Rmax(θ, β) represents the radius from the maximal axial velocity point Eβ to the boundary of the rotational plane Πβ, which increases with the rotation angle β. The formula of Rmax(θ, β) is expressed as follows:

Rmax ( , ) = Rmax ( , 0) ×

R( ) Ri

(2)

where Rmax(θ, 0) is the radius of the tunnel face. R(β) is the corresponding length of Rmax(θ, β) of the rotational plane Πβ when θ = 0. And R(β) is assumed to linearly increase from Ri at the tunnel face to Rf at the ground surface as follows:

Fig. 5. Contour of the maximum shear strain increment of the NATM tunnels: (a) C/D = 1.0; (b) C/D = 2.0.

R ( ) = Ri + (Rf Ri )/( /2) Ri = R1 + L1 Rf = R1 + L1 + C

maximum shear strain increment is more suitable for plotting the range of the failure pattern of the frictional soils, which have an obvious shear band in the soil ground. But for the clayey soils, the outermost layer of the displacement contour corresponds better to the failure zone, where the soil deformation is continuous like a flow and no obvious shear band occurs. Similar situations are also found in other numerical simulations of the tunnel face [19,22,24,43].

(3)

where L1 is the distance between O1 and E at the tunnel face. And the radius Rmax(θ, 0) of the tunnel face can be calculated as follows:

L1 cos Rmax ( , 0) = r2 cos( +

3. Limit analysis of the NATM tunnels

0)

+

L12cos2

+

r22cos2 (

r3 cos( ) +

r32cos2

L12 + R12 0 +

0)

r22

r32 + R32

+

R22 2

< 1

1

<

2

< (4)

3.1. Problem statement

where r2 and r3 are the distances from O2 and O3 to center E, respectively. θ0 is the angle between the EO1 and EO2. θ1 is the angle between EO1 and EP1, and θ2 is the angle between EO1 and EP2. The geometric relation of the horseshoe-shaped tunnel face is shown in Fig. 7. Substituting Eq. (4) into Eq. (2) yields, after some simplifications:

According to the research mentioned above, it could be found that there are many differences between the horseshoe-shaped tunnels and the circular tunnels in terms of the velocity field, the 3D collapse region and the assessment on the stability of the tunnel face. Therefore, in order to estimate the stability of the NATM tunnels with an arbitrary cross-section, it is necessary to build a more realistic continuous velocity field in the framework of the kinematic theorem of limit analysis. For the regular circular face of shield tunnels, Klar and Klein [16] and 5

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 6. 3D failure mechanism for an NATM tunnel.

velocity of the failure mechanism in Fig. 6 is divided into three components: the axial component vβ perpendicular to the rotational plane Πβ (velocity profile) and pointing towards the tunnel face, which is assumed to follow a parabolic distribution in the plane Πβ; the radial component vr located on the plane Πβ and pointing towards center Eβ; and the orthoradial component vθ located on the plane Πβ and orthogonal to the radial component, which is assumed to be zero. According to the results of FLAC3D, the maximal axial velocity point Eβ of every velocity profile Πβ is assumed to be located at a distance of R(β) below the vertex of the plane. Thus the formulas of axial component vβ and the orthoradial component vθ of the velocity vector are defined as follows: r2

v (r , , ) = vm ( ) × 1

2 ( , ) Rmax

(6)

v (r , , ) = 0

where vm(β) is a normalization function used to ensure that the velocity flux through any plane Πβ is constant and independent of angle β based on the principle of mass conservation, which is given by:

1

vm ( ) =

2

Fig. 7. Calculation of Rmax of an NATM tunnel.

Rmax ( , ) =

Ri + 2 (Rf

Rmax ( , ) 0

0

× r2 cos( +

0)

+

L12cos2

+

r22cos2 ( +

r3 cos( ) +

r32cos2

r2

rdrd

2 ( , ) Rmax

(7)

Accomplishing the integration of vm(β) in Eq. (7), the formula of vβ(r, θ, β) is expressed as follows:

Ri )

Ri L1 cos

1

L12 + R12 0 0)

r22 + R22

r32 + R32

2

< 1

<

<

r2

f( ) × 1 A

v (r , , ) =

1

2 R max (

(8)

, )

where f(β) is a function of β. And A is a constant concerned with the geometry of the tunnel face.

2

(5)

f( )=

Any point in the space can be expressed in terms of the angle β (angle between the plane Πβ and the plane of the tunnel face) and the polar coordinates (r, θ) in the plane Πβ, which are defined with respect to center Eβ of the plane Πβ.

2R 2 i

[ Ri + 2 (Rf

1

1

A = 2{ +

2 1

[L1 cos

[r2 cos( +

+

3.2.2. Velocity field The proposed 3D failure mechanism presents the plastic deformation space of the soil ground. Due to the fact that the face collapse in purely cohesive involves a continuous deformation of the soils, the

0

2

Ri)]2

0)

+ +

[r3 cos( ) +

and

L12cos2 r22cos2 ( + r32cos2

L12 + R12 ]2 d 0)

r22 + R22 ]2 d

r32 + R32 ]2 d

}

According to the normality condition of the limit analysis, the plastic deformation of a purely cohesive soil subjected to the undrained conditions occurs without any change of volume. Thus the normality 6

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

condition can be implemented in the equation of strain rate as follows:

div ( ) =

rr

+

+

is the rate of the energy dissipation within the soil mass undergoing a plastic deformation. In the stability analysis of an NATM tunnel, the external forces includes the gravity of the soil unit weight γ and the surface surcharge σs. The rate of work of the soil mass W is given by:

(9)

=0

The strain rate tensor in the orthogonal curvilinear coordinates (r, θ, β) can be expressed as follows: rr

=

r

W =

r

(10)

r

rr

=

r

=

1 r cos )

2(Rf

1 2r

(

v cos =

=

1 2r (Rf

r cos )

=

Rf

vr

r

1 r cos

vr r

=

+ 1 r

v

v r

+r

(v

r

vr

+

v

v r sin

(

v sin

)

r cos )

vr + (Rf r

r cos )

r

+ (Rf

r cos ) +

v

v

Ws =

)

2r cos ) vr + r

v

K( , ) =

2f ( ) Rmax ( , )

Rmax ( , )

=0

4 ( , ) ARmax

df ( ) d

(11)

We = W + Ws

As

WD = 2 + (Rf

(15)

v (r , ,

=

2

) dAs =2

Rmax ( , = ) 2

s

0

0

v (r , ,

V

2cu·max(| i|) dV = 2

Rmax ( , )

2

0

0

0

r cos ) cos ]·2 max(| i|) dr·rd ·(Rf

cu (y ) = cu0 + ( - y ) and M ( ) =

=

2

) dr· (16)

[cu0

r cos ) d

(17)

where max(| i|) denotes the maximal absolute value of the principal strain rate i , and cu is the undrained shear strength of the clay, which increases linearly with the ground depth y

1 df ( ) · . A d

(18)

where cu0 is the undrained shear strength of the clay at the ground surface and ρ is the gradient of the undrained shear strength varying with the ground depth. Besides, soils are degenerated into homogeneous case when ρ = 0.

3.2.3. Calculation on the rate of work Based on the upper bound theorem of the limit analysis, a necessary stability condition for the tunnel face is considered Chen [1]:

WD

r cos ) d

where vβ(r, θ, β = π/2) and As represent the velocity and area of the soil mass on the ground surface, respectively. The rate of the energy dissipation WD in the failure mechanism can be calculated based on the reference Chen [1]:

(12) 2 ( , ) Rmax

+ vr cos cos

rd

K ( , ) r 3 2M ( ) r 4(Rf r cos )

where

+ vr cos cos ) dr· rd ·

The rate of work Ws of a possible uniform surcharge σs acting on the ground surface is given by:

Substituting Eq. (8) into Eq. (11), the radial component vr is obtained:

vr (r , , ) =

(v sin

(14)

dV = dr·rd ·(Rf

Thus the radial component vr located on the plane Πβ can be derived according to the equilibrium equation of the strain rate tensor according to Eqs. (9) and (10):

r (Rf

0

r cos ) d

v Y = v sin v

)

vr cos

0

where VY is the vertical component of the velocity for a given soil element and dV is the volume of an elementary soil mass. And the formulas of VY and dV can be expressed as follows:

+ (Rf v

Rmax ( , )

0

(Rf

where

v Y dV 2

=2

r

r

V

3.3. Definition of the safety factor F

(13)

The safety factor is a crucial and effective criterion to assess the face

where We represents the total rate of work of the external forces and WD

Fig. 8. Comparisons of SRM-F with the results of FLAC3D: (a) γD/cu0 = 2.0; (b) γD/cu0 = 5.0. 7

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 9. Comparisons of ratio-defined F with the results of FLAC3D: (a) γD/cu0 = 2.0; (b) γD/cu0 = 5.0.

Fig. 11. Comparison of the safety factor F with the results of the circular face (γD/cu0 = 5).

Fig. 10. Comparison of the safety factor F for C/D = 2.0 with the results of other existing researches.

continuously reduced until the following relations are obtained:

stability of the NATM tunnels. Two kinds of safety factors are introduced in this paper to estimate the stability of tunnel face based on the upper bound method of the limit analysis. The first effective technique is the strength reduction method (SRM), which has been widely adopted by many researchers to investigate the stability of slope (Lo [44]; Chen et al. [45]; Su and Liao [46]; Al-Karni and Al-Shamrani [47]; Farzaneh and Askari [48]) and tunnel face (Khezri et al. [37]; Pan and Dias [38]). The second method is to define the safety factor as the ratio of the load-resisting capacity to the effect of externally destructive forces (Anthoine [49]; Michalowski and Drescher [50]; Han et al. [51,52]), which is also a convenient way to provide the stability estimation on tunnel face.

| t| = | D ·N +

s ·Ns

cu0/ F · ND|

(19)

0.1

where Nγ, Ns and ND are the non-dimensional coefficients that represent influences of soil unit weight γ, surface surcharge σs and cohesion cu0, respectively, which can be determined by mathematical integral once the geometric parameters of the NATM tunnels including the face shape, the tunnel size and the buried depth are given. The formulas are given as follows:

N =

2 2 D 0 0

Ns =

3.3.1. The strength reduction method (SRM) The SRM combined with the upper bound theorem of the limit analysis is adopted to determine the safety factor of an NATM tunnel. First of all, the initial values of the safety factor F is provided with 0.1, which gives a high value of the undrained shear strength cu and the tunnel face will be in a stable condition. And then the safety factor F is incrementally increased by 0.01 and the strength of the clays is

ND =

2 02 0

Rmax ( , ) (v sin + vr cos cos ) dr·rd ·(Rf 0 R ( , = 0) 2 0 0 max v (r , , = 0) dr·rd

2 0 2 0

Rmax ( , ) [1 + 0

2 0

Rmax ( , = ) 2 v 0

(r ,

, =

Rmax ( , = 0) v (r , , 0

/ cu0 (Rf

2

r cos ) d

) dr·rd

= 0) dr·rd

r cos ) cos ]·2 max(| i |) dr·rd ·(Rf

Rmax ( , = 0) v (r , , 0

r cos ) d

= 0) dr·rd

When Eq. (19) is satisfied, a highest upper bound of the safety factor F is obtained which means that the excavation face of an NATM tunnel reaches a critical failure state and the tunnel face will collapse with 8

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 12. Influences of the geometric parameters of the face shape on the safety factor F*: (a) R1; (b) R2; (c) R3; (d) α2.

only a small increment of safety factor F.

4. Results and discussions This section will compare the safety factors calculated by the proposed kinematic approach using continuous velocity field with the results obtained by FLAC3D and other existing studies to validate the reasonability and superiority of the proposed failure mechanism of the NATM tunnels. The safety factor is calculated based on a horseshoeshape NATM tunnel (shown in Table 1) buried in a clay with a unit weight γ = 10 kN/m and for two values of cu0 = 0.2γD and 0.5γD (corresponding to about 18 kPa and 46 kPa, respectively). The surcharge σs of the ground surface is assumed to be equal to zero. All the geologic and geometric parameters are transformed into non-dimensional factors to present the stability charts by the normalization method for better comparisons of stability evaluation on tunnel face and practical references for similar engineering projects.

3.3.2. The ratio of the load-resisting capacity to the externally destructive effect In this case, the safety factor F is defined as the ratio of the rate of energy dissipation WD to the rate of work of external forces We as follows:

F =

WD We

(20)

The larger rate of energy dissipation represents better load-resisting capacity of the failure system. The safety factor F can be achieved by numerically integration of the rate of energy dissipation and the rate of work of external forces with the help of the optimization tool implemented in the Matlab. It is worth noting that the value of this ratiodefined safety factor F represents the conditions of the tunnel face as follows:

< 1 unstable F = 1 limit state > 1 stable

4.1. Comparisons with the numerical simulation Figs. 8 and 9 present the comparisons of the safety factor F deduced by the numerical simulation and the analytical solutions, which are defined by the SRM and the ratio of the rate of energy dissipation to the work rate of external forces, respectively. In Figs. 8 and 9, it is shown

(21) 9

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 13. Geometric shape of the top heading.

Vermeer et al. [20] proposed a simple formula of the fitting curve based on the numerical calculations and the strength reduction technique to calculate the safety factor for an open face tunneling as follows:

F =

0.9 tan + 18cu / D 2 + 3(d / D)6 tan / F

(22)

where cu, φ, and γ are cohesion, frictional angle and unit weight of soils. d is the unlined distance of tunnels and D is the equivalent diameter of an NATM tunnels. Pan and Dias [38] adopted the limit analysis using the spatial discretization technique (rotational rigid block mechanism) in combination with the strength reduction method to calculate the safety factor F of an NATM tunnel of height D and width L. t

Fig. 14. Comparison of the safety factor F with different excavation methods.

RV sin

=

cF cos F (

F

RS

R cos )

=0

cF = c /F = arctan(tan / F )

(23)

where R is the rotational radius of the discrete element and β is the rotational angle. V and S are the unit volume and area of the discrete element, respectively. Fig. 10 compares the ratio-defined safety factor F and the value defined by the SRM with the results obtained from Vermeer et al. [20], Pan and Dias [38] and FLAC3D. It can be seen that both ratio-defined F and SRM methods generally lead to a same assessment on the safety factor F. In addition, the proposed mechanism of the continuous velocity field provides the most critical solution compared with other results, which is quite close to the results deduced by FLAC3D. The dimensionless parameter γD/cu0 has obvious influence on the safety factor F, but the safety factor F is prone to constant when γD/cu0 ≥ 7, which implies that the influence of dimensionless geologic parameter γD/cu0 is negligible when its value is relatively large.

that two methods provide a basically consistent value of safety factor F. Moreover, it can be found that the analytical assumption of the position of maximal axial velocity has obvious influence on the stability of the face of the NATM tunnels. The critical value of L1 varies with the buried depth ratio C/D. For the shallow tunnels (C/D = 0.5), L1/d2 = 0.8 leads to a critical value of F. But when the buried depth ratio is C/D ≥ 3.0, the critical value of L1/d2 is equal to 0.2. Considering the computation time (a larger value of L1/d2 leads to a larger gradient in the radial direction on the bottom of the cross-section Πβ, which significantly increases the computation time) and the critical value of safety factor F, a compromising value of L1/d2 between the computation time and the accuracy of the analytical solution is thus taken to be equal to 0.6, which gives a solution of F accurate enough. On the other hand, the differences between the most critical analytical solutions and the results of FLAC3D increase from 5% to 18% with the buried depth ratio C/D as shown in Figs. 8 and 9. When the buried depth ratio C/D ≥ 2.0, the difference will exceed 12%, which shows some shortcomings happen to the proposed failure mechanisms for deep-buried NATM tunnels.

4.2. Comparisons with the results of circular cross-section For purely cohesive soils, a traditional method used to estimate the 10

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 15. Stability charts of the coefficients Nγ and ND: (a) Nγ; (b) ND; (c) F.

stability of the tunnel face is to calculate the load factor N, which is a dimensionless parameter N (stability number) introduced by Davis et al. [3] as follows:

N =

s

+ H cu0

t

(24)

where H = C + D/2. σs is the surcharge of ground surface and σt is the critical supporting pressure applied on the tunnel face. By using Eqs. (19) and (24), the following formulas can be obtained:

N = H / D = C / D + 0.5 ND = N Ns = 1

(25)

And the critical safety factor F is calculated as follows:

F = Fig. 16. The coefficient ND for the nonhomogeneous soils.

cu 0· N D· N +

= s

(

D N cu0

C D

+ 0.5 +

cu0· Nc + D

(

)

) s

(26)

Eq. (26) can be used in the calculation of the safety factor F of the tunnel face both in the non-homogeneous and homogeneous soils. The load factor N (N = ND) can be divided into two parts. Nc and Nρ are the 11

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 17. The safety factor of the NATM tunnels in the nonhomogeneous soils: (a) γD/cu0 = 2.0; (b) γD/cu0 = 5.0.

By using Eq. (26), Fig. 11 presents the comparisons of the safety factor F of the NATM tunnels with the results of other existing research of the circular tunnel face. It is shown that the horseshoe-shaped crosssection obviously improves the stability of the tunnel face compared with that of the circular tunnel. The maximum improvement on the face stability by the horseshoe-shaped tunnels can achieve 24% when C/D is equal to 0.5. Besides, the continuous velocity field builds a more realistic and critical failure mechanism than the rigid blocks, thus a minimum of the safety factor F is obtained by the M2 (Mollon et al. [14]) among all the results of the studies on the circular tunnel face. This fact also shows that the safety factors of the horseshoe-shaped tunnels are smaller than some results of the circular tunnel fact, and this phenomenon is more significant when the buried depth C/D is larger than 2.0, which can lead to a difference of approximately 5–14% between the horseshoe-shaped tunnels and the circular tunnels when C/D is from 2.0 to 3.0.

Table 2 Non-dimensional coefficients for the practical design of the face stability. C/D

ρD/cu0

Nγ

Nc

Nρ

ND(=Nc + ρDNρ/cu0)

0.5

0 0.25 0.50 0.75 1.0

1.08

7.37

5.85

7.37 8.83 10.29 11.75 13.22

1.0

0 0.25 0.50 0.75 1.0

1.58

9.42

10.68

9.42 12.09 14.76 17.42 20.09

2.0

0 0.25 0.50 0.75 1.0

2.58

12.45

23.16

12.45 18.24 24.03 29.82 35.61

3.0

0 0.25 0.50 0.75 1.0

3.58

14.67

38.41

14.67 24.27 33.87 43.47 53.07

0 0.25 0.50 0.75 1.0

4.58

0 0.25 0.50 0.75 1.0

5.58

4.0

5.0

16.42

17.86

55.77

74.85

4.3. Effects of the geometric parameters For the analytical solution, the stability of the tunnel face is closely related with the area of the tunnel face. To investigate the influence of the geometric parameters on the face of the NATM tunnels, a safety factor F* of a unit area of the tunnel face is defined as follows:

16.42 30.36 44.30 58.24 72.18

F* =

N =

2 2 D 0 0

Rmax ( , ) 2 max(| i |) dr·rd ·(Rf r cos ) d 0 R ( , = 0) 2 0 0 max v (r , , = 0) dr·rd

2 02 0

Rmax ( , ) (Rf r cos ) cos ·2 max(| i |) dr·rd ·(Rf 0 R ( , = 0) 2 0 0 max v (r , , = 0) dr·rd

(28)

where Area represents the area of the tunnel face. Fig. 12 shows the variations of the safety factor F* with four geometric parameters. The lengths of R2 and R3 influence on the stability of the tunnel face slightly, but the R1 and α2 influence the stability of the tunnel face obviously. When the length of R1 increases, the tunnel vault becomes flatter leading to a worse stress condition, and the area of the tunnel face will greatly enlarge accordingly, which thus lowers the stability of the tunnel face. With regard to the geometric parameter α2, it mainly controls the range of the invert and side wall of the NATM tunnels. It can be seen from Fig. 12(d) that the stability of the tunnel face gets decreased with the angle α2, which implies that the increase of range of tunnel invert is beneficial to the stability of the NATM tunnels but the enlargement of the side wall structure doesn’t improve the stability of the tunnel face. Therefore, the design of the parameter α2 should be a compromise among the operating requirements of the tunnels, the stability control of the tunnel face and the costs for the

17.86 36.58 55.29 74.00 92.72

non-dimensional coefficients that represent the influences of the homogeneous soil strength and the soil inhomogeneity, respectively. The formulas of Nc and Nρ are expressed as follows:

Nc =

F Area

r cos ) d

(27) 12

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 18. Comparisons of Nγ, Nc and Nρ between the solutions of this paper and design equation proposed by Ukritchon et al. [24]: (a) Nγ; (b) Nc; (c) Nρ.

excavation. Furthermore, let the angle α2 decrease to 0 and then the excavation method becomes the common benching tunneling method in the practical engineering of the NATM tunnels as shown in Fig. 13. By using the proposed analytical method, a comparison of the face stability between the top heading excavation and the full face excavation is conducted, as shown in Fig. 14. The safety factor F of the top heading excavation (the benching tunneling method) is obviously greater than that of the full face excavation with a maximum difference of 30%, which means that the stability of the tunnel face excavated by the benching tunneling method can be theoretically improved by up to 30% compared with the full face excavation.

linear fitting curve, which is corresponding to Eq. (23). Moreover, substituting the values of Nγ and ND into Eq. (24), the corresponding safety factor F with different parameters γD/cu0 is obtained as shown in Fig. 15(c). For soils with the undrained strength increasing linearly with depth (ρ ≠ 0), the coefficient ND with different ρD/cu0 is shown in Fig. 16. It is obviously expected that the coefficient ND increases with ρD/cu0, which illustrates that the inhomogeneous soils have greater influence on the stability of tunnel face than the homogeneous soils. Since the parameter ρD/cu0 has no effect on the coefficient Nγ, the safety factor F with different ρD/cu0 can be calculated by Eq. (24). Fig. 17 plots the safety factor F against the ratio of C/D. It can be found that the inhomogeneity of soils ground has significant influence on the stability of the tunnel face. When the parameter is ρD/cu0 ≥ 2.0, the safety factor F increases with the ratio of C/D, which shows an inverse situation compared with Fig. 15(c). From the view of the upper bound solutions, it is suggested that the neglect of the inhomogeneity of soils will lead a conservative evaluation on the stability of the tunnel face. Thus it is necessary to take the inhomogeneity of soils into account in the analytical method considering the economic benefit of the practical engineering. Besides, for the practical use in the geotechnical engineering, the dimensionless coefficients Nγ and ND versus C/D are provided as shown

4.4. Effects of the geologic parameters This section will present the stability charts of the influences of the geologic parameters on the horseshoe-shape NATM tunnel. The face geometry shown in Fig. 2 is adopted herein. The values of Nγ and ND for different γD/cu0 are plotted as a function of the buried depth ratio C/D, as shown in Fig. 15(a) and (b). It is shown that the parameter γD/cu0 has negligible influence on Nγ and ND, but the coefficients of Nγ and ND depend on the ratio of C/D. As shown in Fig. 15(a), the variation of the coefficient Nγ with C/D conforms to a 13

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 19. Layout of the derived continuous velocity field for an NATM tunnel: (a) C/D = 2.0, L1/d2 = 0, γD/cu0 = 2; (b) C/D = 1.0, L1/d2 = 0.6, γD/cu0 = 2.

in Table 2. Fig. 18 shows the comparisons between the solutions obtained from this paper and design equation proposed by Ukritchon et al. [24]. It can be seen that the variation rules of Nγ, Nc and Nρ with C/D are similar with that of the design equation by Ukritchon et al. [24]. Non-linear increasing functions between Nc, Nρ and C/D are observed, whereas a linearly increasing relationship between Nγ and C/D is found. Besides, Nγ obtained from this paper shows a good agreement with the design equation. And Nc and Nρ are comparable to the design equation when C/D ≤ 2.0. But it is worth noting that Nc and Nρ of the present study are

obviously greater than the results of the design equation when C/ D = 2–5, which thus leads to a greater value of F according to Eq. (24). The differences between the present study and the design equation are approximately 25–31% for Nc and Nρ. Two reasons are considered to cause this difference. The first reason is the difference between the kinematical approach of the upper bound method and the finite element method of the numerical simulation, which is also seen in the comparison between the design equation by Ukritchon et al. [24] and the upper bound solution by Mollon et al. [14]. The second reason is the difference between the horseshoe-shaped tunnels and the circular 14

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 20. Comparison of the 3D failure mechanisms between the horseshoe-shaped tunnel and the circular tunnel: (a) Horseshoe-shaped tunnel face; (b) Circular tunnel face.

tunnels. As shown in Fig. 11, the horseshoe-shaped tunnels improve the face stability of the circular tunnels. Thus Nc and Nρ become greater accordingly for a horseshoe-shaped NATM tunnel.

Fig. 19 plots two sketches of the failure mechanism for the NATM tunnels with different buried depths (C/D = 1.0 and 2.0) and different positions of maximal (L1/D = 0 and 0.3). The distributions of the axial (vβ) and radial (vr) components are also shown in Fig. 19. The axial component vβ of velocity follows a parabola distribution in any crosssection passing through origin O, which is rather consistent with the analytical assumption in the previous section. It is also worth noting that the derived radial velocity is zero at the point of maximal axial velocity and outside the envelope of the failure mechanism based on the normality condition of the limit analysis. The envelopes of vβ and vr well depict the continuous plastic deformation of soils influenced by the face collapse of the NATM tunnels. Moreover, the ranges of the failure zone and the surface subsidence increase with the buried depth ratio C/D.

velocity field became the most realistic when L1/D was equal to 0.2 times of the tunnel radius D/2. For a horseshoe-shaped tunnel, L1/d2 is taken to be 0.6 according to the results mentioned above. The failure mechanism of the horseshoe-shaped tunnel is similar with that of the circular tunnel. But compared with the common circular tunnel face, the failure mechanism of the horseshoe-shape tunnel face extends to a larger region both in transverse and longitudinal range of the soils. Fig. 21 shows the comparison of the failure mechanism obtained from the analytical approach and numerical simulation. It can be found that the failure pattern of the analytical method propagates from the tunnel face towards the ground surface, which corresponds well with that of the numerical simulation when C/D is equal to 1.0. But when C/D is equal or greater than 2.0, the proposed failure mechanism will extend to a larger range than the results of the numerical simulation and this unconformity will be more significant with a larger buried depth ratio C/D. It is because the failure zone will not infinitely extend its region in the soils ground with the buried depth C due to the soil arch effect. Thus, to solve this problem, an improved velocity field will be presented and compared with the results of the numerical simulation in the next section.

4.6. Comparisons of the failure mechanism

5. Improvement of the velocity field for NATM tunnels

Fig. 20 shows the comparison of the failure mechanisms of the horseshoe-shaped tunnel and circular tunnel with the same equivalent diameter D. The failure mechanism of the tunnel face mainly depends on the face shape and the position L1 of the maximal axial velocity point. For a circular tunnel, Mollon et al. [14] recommended that the

In the previous section, it was found that the proposed analytical solutions of the safety factor F for C/D ≥ 2.0 have a difference of more than 10% compared with the results of FLAC3D. And the failure zone for a deepburied NATM tunnel is different with that obtained from the numerical simulation. This is because that the range of the proposed continuous

4.5. Profile of the failure mechanism

15

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 21. Comparison of the failure mechanisms between the analytical method and the numerical simulation: (a) C/D = 1.0; (b) C/D = 2.0; (c) C/D = 3.0; (d) C/ D = 4.0.

Fig. 22. Improved velocity field for an NATM tunnel.

16

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 24. Critical L0 versus C/D.

the ground surface with a constant area of velocity profile, which is similar with the soil arch model presented by Terzaghi [53]. The 3D improved failure mechanism is presented in Fig. 23. The pressure arch theory was usually adopted in the limit equilibrium method and those studies [54–57] mainly focused on its stress distribution in the failure mechanism, however, for the kinematic approach of the limit analysis, the velocity distribution within the zone of soil arch was the key issue. Therefore, by assuming the velocity profile within the soil arch and the vertical boundary of the arch zone, the improved velocity field will be adopted in the limit analysis. The Zone I of the improved mechanism is the same as that proposed in the previous section. But it is worth pointing out that C should be replaced by a height of L0 in Eq. (3). Moreover, the work rates of the external forces and the energy dissipation of the Zone II can be simply obtained by substituting β = π/2 into the calculations of Zone I as follows: C L0

W =2

Ws = 2 C L0

WD = 2

0

0

0

s

0

0

0

Rmax ( , = 2 )

v (r , ,

= 2 ) dr·rd ·dy

Rmax ( , = 2 )

v (r , ,

= 2 ) dr· rd

Rmax ( , = 2 )

[cu0 + (C

L0

y ) cos ]·2 max(| i|) dr·

rd · dy (29) The corresponding strain rate tensor in the cylindrical coordinates (r, θ, y) of the Zone II can be expressed as follows:

Fig. 23. 3D improved failure mechanisms for the deep-buried NATM tunnels: (a) C/D = 2.0; (b) C/D = 4.0.

rr

=

velocity field is related with the buried depth C/D of the tunnels, and the failure zone will infinitely extend the region in the soil ground with C/D. Therefore, in order to solve this problem, this section will improve the continuous velocity field for the deep-buried NATM tunnels by posing the limit on the boundary of the continuous velocity field.

r

ry

y

yy

r

y

yr

(30)

where rr r

5.1. Improvement of the failure mechanism

= ry

Fig. 22 plots the improved velocity field of the NATM tunnels, which involves two zones of the velocity distribution: Zone I is the same as the proposed continuous velocity field, but it is assumed that the axial component of the velocity turns vertical and the range of the velocity profile keeps constant at a height of L0 above the tunnel crown rather than the buried depth C; Zone II propagates vertically towards

(

1 2r

vr

=

1 r

= y

=

1 2r

17

vr y

(v

r

(r

yy

v r

+r

(

1 2

vr r

=

=

v

+

vy

+

v

v y

+ vy y

r

)

) vy

)

)

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 25. The improvement of the safety factor F: (a) γD/cu0 = 5; (b) γD/cu0 = 2.

Fig. 26. Comparison between the improved velocity field and the continuous velocity field for the deep-buried NATM tunnels: (a) C/D = 3.0; (b) C/D = 4.0.

5.2. Results and discussion

collapse especially for the deep-buried NATM tunnels. As shown in Fig. 27, the maximum shear strain mainly concentrates in front of the tunnel face for the deep-buried NATM tunnels, which is similar with that of Fig. 5. It is considered that the contour of the maximum shear strain increment is not suitable for plotting the range of the failure pattern of the clayey soils, which has no obvious shear band when the soils achieve the critical state. Thus for the purely cohesive soils with a continuous flow characteristic, the outermost layer of the displacement contour shown in Fig. 26 represents the failure zone of the tunnel face better.

The safety factor F can be calculated by substituting Eq. (29) into Eqs. (19) and (20). To determine the critical upper bound solutions of F, a maximization procedure will be performed with the help of the optimization tool implemented in the software Matlab with respect to the parameter L0. Fig. 24 shows the corresponding values of L0 when the safety factor F reaches a critical value. Fig. 25 presents the comparisons of the safety factor F. It is shown that the improved velocity field obviously provides a more critical upper bound solution, which is closer to the results of FLAC3D and the corresponding maximum difference is no more than 8%. In general, by imposing the limit on the boundary of the continuous velocity field, the new proposed failure mechanism obviously improves the shortcoming for the deep-buried failure mechanism. Fig. 26 presents the comparison of the improved velocity field with the continuous velocity field. Compared with the continuous velocity field, the improved velocity field restricts the velocity distribution in a narrower region, which is mainly concentrated in the vicinity of the tunnel face and the zone directly above the tunnel face. And this concentrated failure pattern is closer to the realistic situation of the face

6. Conclusions Based on the kinematic method of the limit analysis, this paper proposed a continuous velocity field to assess the face stability of the NATM tunnels in undrained clays by introducing two kinds of safety factors F. The calculated upper bound solutions were compared with the results of 3D numerical simulations for a typical horseshoe-shaped tunnel face and the solutions obtained from other research programs. The influences of the geometric parameters on the stability of the tunnel face are investigated and the stability charts for different C/D, γD/cu0 18

Computers and Geotechnics 114 (2019) 103136

W. Li, et al.

Fig. 27. Comparison among the improved velocity field, the continuous velocity field and the contour of the maximum shear strain increment: (a) C/D = 3.0; (b) C/ D = 4.0.

and ρD/cu0 are provided. Besides, an improved velocity field is proposed for the NATM tunnels due to the shortcoming that a big difference occurs between the analytical solutions and numerical results. The main conclusions are drawn as follows:

References [1] Chen WF. Limit analysis and soil plasticity. Amsterdam: Elsevier Science; 1975. [2] Drucker DC, Prager W, Greenberg HJ. Extended limit design theorems for continuous media. Q J Mech Appl Math 1952;9(4):381–9. https://doi.org/10.1090/ qam/45573. [3] Davis EH, Gunn MJ, Mair RJ, Seneviratne HN. The stability of shallow tunnels and underground openings in cohesive material. Geotechnique 1980;30(4):397–416. https://doi.org/10.1680/geot.1980.30.4.397. [4] Leca E, Dormieux L. Upper and lower bound solutions for the face stability of shallow circular tunnels in frictional material. Geotechnique 1990;40(4):581–606. https://doi.org/10.1680/geot.1990.40.4.581. [5] Soubra AH, Dias D, Emeriault F. Three-dimensional face stability analysis of circular tunnels by a kinematical approach. GeoCongress 2008. United States, New Orleans. doi: 10.1061/40972(311)112. [6] Mollon G, Dias D, Soubra AH. Probabilistic analysis and design of circular tunnels against face stability. Int J Geomech 2009;9(6):237–49. https://doi.org/10.1061/ 41022(336)45. [7] Mollon G, Dias D, Soubra AH. Face stability analysis of circular tunnels driven by a pressurized shield. J Geotech Geoenviron 2010;136(1):215–29. https://doi.org/10. 1061/(ASCE)GT.1943-5606.0000194. [8] Zhang CP, Han KH, Zhang DL. Face stability analysis of shallow circular tunnels in cohesive-frictional soils. Tunn Undergr Sp Tech 2015;50:345–57. https://doi.org/ 10.1016/j.tust.2015.08.007. [9] Han KH, Zhang CP, Zhang DL. Upper-bound solutions for the face stability of a shield tunnel in multilayered cohesive-frictional soils. Comput Geotech 2016;79:1–9. https://doi.org/10.1016/j.compgeo.2016.05.018. [10] Han KH, Zhang CP, Li W, Guo CX. Face stability analysis of shield tunnels in homogeneous soil overlaid by multilayered cohesive-frictional soils. Math Probl Eng 2016;2016:9. https://doi.org/10.1155/2016/1378274. Article ID 1378274. [11] Subrin D, Wong H. Tunnel face stability in frictional material: a new 3D failure mechanism. C.R. Mecanique 2002;330:513–519 [in French]. [12] Mollon G, Dias D, Soubra AH. Rotational failure mechanisms for the face stability analysis of tunnels driven by pressurized shields. Int J Numer Analyt Met Geomech 2011;35(12):1363–88. https://doi.org/10.1002/nag.962. [13] Osman AS, Mair RJ, Bolton MD. On the kinematics of 2D tunnel collapse in undrained clay. Geotechnique 2006;56(9):585–95. https://doi.org/10.1680/geot. 2006.56.9.585. [14] Mollon G, Dias D, Soubra AH. Continuous velocity fields for collapse and blowout of a pressurized tunnel face in purely cohesive soil. Int J Numer Analyt Met Geomech 2013;37(13):2061–83. https://doi.org/10.1002/nag.2121. [15] Klar A, Osman AS, Bolton M. 2D and 3D upper bound solutions for tunnel excavation using ‘elastic’ flow fields. Int J Numer Analyt Methods Geomech 2007;31(12):1367–74. https://doi.org/10.1002/nag.597. [16] Klar A, Klein B. Energy-based volume loss prediction for tunnel face advancement in clays. Géotechnique 2014;64(10):776–86. https://doi.org/10.1680/geot.14.P.024. [17] Zhang F, Gao YF, Wu YX, Zhang N. Upper-bound solutions for face stability of circular tunnels in undrained clays. Géotechnique 2018;68(1):76–85. https://doi. org/10.1680/jgeot.16.T.028. [18] Zhang F, Gao YF, Wu YX, Wang ZX. Face stability analysis of large-diameter slurry shield-driven tunnels with linearly increasing undrained strength. Tunn Undergr Sp Tech 2018;78:178–87. https://doi.org/10.1016/j.tust.2018.04.018. [19] Huang MS, Tang Z, Zhou WX, Yuan JY. Upper bound solutions for face stability of circular tunnels in non-homogeneous and anisotropic clays. Comput Geotech 2018;98:189–96. https://doi.org/10.1016/j.compgeo.2018.02.015.

(1) Compared with other existing studies and FLAC3D, the proposed mechanism of the continuous velocity field provides the most critical solution, which is closest to the results of FLAC3D and have a most realistic deformation mechanism. Compared with the results of other studies of the circular cross-section, the horseshoe-shaped cross-section significantly improved the stability of the tunnel face, which implied that the investigation on the critical velocity field of the NATM tunnel face was necessary and worthwhile. (2) In terms of the geometry of the tunnel face, R1 and α2 had a great influence the stability of the tunnel face. A flatter tunnel vault and a smaller range of tunnel invert would decrease the stability of the tunnel face by up to 30%. For the geologic parameters, γD/cu0 has a negligible influence on Nγ and ND, but the coefficients of Nγ and ND depend on the buried depth ratio C/D. And the coefficient Nγ with C/D conforms to a linear fitting curve. Moreover, the inhomogeneity ρD/cu0 of soils ground has a significant influence on the stability and the neglect of the inhomogeneity of soils will lead a conservative assessment on the face stability. (3) By imposing the limit on the boundary of the continuous velocity field, an improved failure mechanism is proposed based on the soil arch theory, and it obviously improves the shortcoming for the deep-buried tunnels with a maximum difference of no more than 8% compared with the results of numerical simulation. (4) This paper adopts the safety factor to assess the face stability of a horseshoe-shaped NATM tunnel. Some certain limitation still happens to the safety factor, because it depends on the specific values of the tunnel problems including the soil undrained strength (cu), the loading condition (γ and σs) and the problem dimension (C and D). Thus the safety factor could be determined only if the tunnel conditions of the practical engineering are given. Acknowledgments The authors acknowledge the financial support provided by Beijing Municipal Natural Science Foundation of China (Grant No. 8172037) and the National Natural Science Foundation of China (Grant No. 51378002). 19

Computers and Geotechnics 114 (2019) 103136

W. Li, et al. [20] Vermeer PA, Ruse N, Marcher T. Tunnel heading stability in drained ground. Felsbau 2002;20(6):8–18. [21] Chen RP, Tang LJ, Ling DS, Chen YM. Face stability analysis of shallow shield tunnels in dry sandy ground using the discrete element method. Comput Geotech 2011;38(2):187–95. https://doi.org/10.1016/j.compgeo.2010.11.003. [22] Ukritchon B, Keawsawasvong S, Yingchaloenkitkhajorn K. Undrained face stability of tunnels in Bangkok subsoils. J Geotech Eng 2017;11(3):262–77. https://doi.org/ 10.1080/19386362.2016.1214773. [23] Ukritchon B, Keawsawasvong S. Design equations for undrained stability of opening in underground walls. Tunn Undergr Sp Tech 2017;70:214–20. https://doi.org/10. 1016/j.tust.2017.08.004. [24] Ukritchon B, Yingchaloenkitkhajorn K, Keawsawasvong S. Three-dimensional undrained tunnel face stability in clay with a linearly increasing shear strength with depth. Comput Geotech 2017;88:146–51. https://doi.org/10.1016/j.compgeo. 2017.03.013. [25] Broms BB, Bennermark H. Stability of clay at vertical openings. J Soil Mech Found Eng 1967;193(SM1):71–94. [26] Mair RJ. Centrifugal modelling of tunnel construction in soft clay PhD thesis University of Cambridge; 1969. [27] Schofield AN. Cambridge geotechnical centrifuge operations. Geotechnique 1980;30(3):227–68. https://doi.org/10.1680/geot.1980.30.3.227. [28] Chambon P, Corté JF. Shallow tunnels in cohesionless soil: stability of tunnel face. J Geotech Eng 1994;120(7):1148–65. https://doi.org/10.1061/(ASCE)07339410(1994) 120:7(1148). [29] Takano D, Otani J, Nagatani H, Mukunoki T. Application of X-ray CT on boundary value problems in geotechnical engineering-research on tunnel face failure. In: Proc Geo ASCE, Atlanta; 2006. [30] Idinger G, Aklik P, Wu W, Borja RI. Centrifuge model test on the face stability of shallow tunnel. Acta Geotech 2011;6:105–17. https://doi.org/10.1007/s11440011-0139-2. [31] Chen RP, Li J, Kong LG, Tang LJ. Experimental study on face instability of shield tunnel in sand. Tunn Undergr Sp Technol 2013;33:12–21. https://doi.org/10.1016/ j.tust.2012.08.001. [32] Ukritchon B, Keawsawasvong S. Stability of retained soils behind underground walls with an opening using lower bound limit analysis and second-order cone programming. Geotech Geol Eng 2018:1–17. https://doi.org/10.1007/s10706-0180710-9. [33] Ukritchon B, Keawsawasvong S. Lower bound limit analysis of an anisotropic undrained strength criterion using second-order cone programming. Int J Nmer Anal Methods Geomech 2018;42:1016–33. https://doi.org/10.1002/nag.2781. [34] Ukritchon B, Keawsawasvong S. Stability of unlined square tunnels in Hoek-Brown rock masses based on lower bound analysis. Comput Geotech 2019;105:249–64. https://doi.org/10.1016/j.compgeo.2018.10.006. [35] Ukritchon B, Keawsawasvong S. Lower bound stability analysis of plane strain headings in Hoek-Brown rock masses. Tunn Undergr Sp Technol 2019;84:99–112. https://doi.org/10.1016/j.tust.2018.11.002. [36] Huang MS, Li S, Yu J, et al. Continuous field based upper bound analysis for threedimensional tunnel face stability in undrained clay. Comput Geotech 2018;94:207–13. https://doi.org/10.1016/j.compgeo.2017.09.014. [37] Khezri N, Mohamad H, HajiHassani M, Fatahi B. The stability of shallow circular tunnels in soil considering variations in cohesion with depth. Tunn Undergr Space Technol 2015;49:230–40. https://doi.org/10.1016/j.tust.2015.04.014. [38] Pan QJ, Dias D. Upper-bound analysis on the face stability of a non-circular tunnel.

[39] [40] [41] [42]

[43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57]

20

Tunn Undergr Sp Technol 2017;62:96–102. https://doi.org/10.1016/j.tust.2016. 11.010. Fang Q, Zhang DL, Li QQ, Wong LNY. Effects of twin tunnels construction beneath existing shield-driven twin tunnels. Tunn Undergr Space Technol 2015;45:128–37. https://doi.org/10.1016/j.tust.2014.10.001. Fang Q, Zhang DL, Wong LNY. Shallow tunnelling method (STM) for subway station construction in soft ground. Tunn Undergr Space Technol 2012;29:10–30. https:// doi.org/10.1016/j.tust.2011.12.007. Li XG, Zhang CP, Yuan DJ. An in-tunnel jacking above tunnel protection methodology for excavating a tunnel under a tunnel in service. Tunn Undergr Space Technol 2013;34:22–37. https://doi.org/10.1016/j.tust.2012.10.004. Zhang DL, Fang Q, Hou YJ, Li PF, Wong LNY. Protection of buildings against damages as a result of adjacent large-span tunneling in shallowly buried soft ground. J Geotech Geoenviron Eng 2013;139(6):903–13. https://doi.org/10.1061/(ASCE)GT. 1943-5606.0000823. Li Y, Emeriault F, Kastner R, Zhang ZX. Stability analysis of large slurry shielddriven tunnel in soft clay. Tunn Undergr Space Technol 2009;24:472–81. https:// doi.org/10.1016/j.tust.2008.10.007. Lo KY. Stability of slopes in anisotropic soils. J Soil Mech Found Div 1965;91(4):85–106. https://doi.org/10.5194/acp-12-11679-2012. Chen WF, Snitbhan N, Fang HY. Stability of slopes in anisotropic, nonhomogeneous soils. Can Geotech J 1975;12(1):146–52. https://doi.org/10.1139/t75-014. Su SF, Liao HJ. Effect of strength anisotropy on undrained slope stability in clay. Geotechnique 1999;49(2):215–30. https://doi.org/10.1680/geot.1999.49.2.215. Al-Karni AA, Al-Shamrani MA. Study of the effect of soil anisotropy on slope stability using method of slices. Comput Geotech 2000;26(2):83–103. https://doi.org/ 10.1016/S0266-352X(99)00046-4. Farzaneh O, Askari F. Three-dimensional analysis of nonhomogeneous slopes. J Geotech Geoenviron Eng 2003;2(137):137–45. https://doi.org/10.1061/(ASCE) 1090-0241(2003) 129:2(137). Anthoine A. Mixed modelling of reinforced soils within the framework of the yield design theory. Comput Geotech 1989;7(1):67–82. https://doi.org/10.1016/0266352X(89)90007-4. Michalowski RL, Drescher A. Three-dimensional stability of slopes and excavations. Géotechnique 2009;59(10):839–50. https://doi.org/10.1680/geot.8.P.136. Han CY, Wang JH, Xia XH, et al. Limit analysis for local and overall stability of a slurry trench in cohesive soil. Int J Geomech 2012;15(5):06014026. https://doi. org/10.1061/(ASCE)GM.1943-5622.0000268. Han CY, Chen JJ, Wang JH, et al. 2D and 3D stability analysis of slurry trench in frictional/cohesive soil. J Zhejiang Univ-SCI A 2013;14(2):94–100. https://doi.org/ 10.1631/jzus.A1200257. Terzaghi K. Theoretical soil mechanics. New Jersey: John Wiley & Sons; 1943. Anagnostou G, Kovari K. Face stability conditions with earth-pressure-balanced shields. Tunn Undergr Space Technol 1996;11(2):165–73. https://doi.org/10. 1016/0886-7798(96)00017-X. Anagnostou G. The contribution of horizontal arching to tunnel face stability. Geotechnik 2012;35(Heft 1):34–44. https://doi.org/10.1002/gete.201100024. Anagnostou G, Perazzelli P. The contribution of horizontal arching to tunnel face stability. Geotechnik 2013;36(Heft 1):40–50. https://doi.org/10.1002/gete. 201200014. Chen RP, Tang LJ, Yin XS, Chen YM, Bian XC. An improved 3D wedge-prism model for the face stability analysis of the shield tunnel in cohesionless soils. Acta Geotech 2015;10:683–92. https://doi.org/10.1007/s11440-014-0304-5.