Face stability analysis of shallow circular tunnels driven by a pressurized shield in purely cohesive soils under undrained conditions

Computers and Geotechnics 107 (2019) 110–127

Contents lists available at ScienceDirect

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

Research Paper

Face stability analysis of shallow circular tunnels driven by a pressurized shield in purely cohesive soils under undrained conditions

T

Wantao Dinga,b, Keqi Liub, Peihe Shib, Mingjiang Lib, Minglei Houb a b

School of Qilu Transportation, Shandong University, Jinan, China Research Center of Geotechnical and Structural Engineering, Shandong University, Jinan, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Slip line Limit analysis Face stability Collapse mechanism Purely cohesive soils

Based on the slip line and limit analysis theories, this paper presents a new 3D collapse mechanism to analyze the face stability of a circular shield-driven tunnel in purely cohesive soil. The new collapse mechanism is a 4variable multizone mechanism that depends on cover-to-diameter ratios (C/D). The proposed mechanism suggests that failure zones at the tunnel face are primarily formed by shear failure. A comparison of this approach with numerical methods and other kinematic approaches shows that the proposed failure mechanism significantly improves the existing upper-bound solutions for the face stability of circular tunnels in purely cohesive soils.

1. Introduction

homogeneous cohesive soils, and this solution was in good agreement with the experimental results by Kimura et al. [3]. Based on the upper and lower-bound theorems of plasticity, Davis et al. [5] proposed design values for the stability of a shallow heading in soft clay under undrained conditions. When conducting an upper-bound analysis using the finite element method (FEM), predetermination of the failure mechanism is not required. Thus, Sloan et al. [6] examined the undrained stability of a shallow heading under plane strain loading using two numerical techniques that employ finite elements in conjunction with the limit theorem of classical plasticity. Augarde et al. [7] investigated the stability of a heading under plane strain undrained conditions by means of the finite element limit analysis method in the framework of classical plasticity theory. Combining the results of a series of centrifuge model tests with corresponding numerical simulations, Lee et al. [8] studied the face stability of a tunnel and the influence of soil arching effects in soft clay soils and determined the boundaries of the positive and negative arching zones. Recently, several kinematic methods for limit analysis that consider a continuous velocity field were proposed by various authors. For example, Klar et al. [9] derived 2D and 3D upper-bound solutions to analyze the face stability of circular tunnels in purely cohesive soils by means of an admissible continuous velocity field, which is achieved using elasticity theory. Mollon et al. [10] suggested two new continuous velocity fields for analyzing the face stability of a tunnel driven by pressurized shields, and these velocity fields more closely represented the actual failures observed in undrained clay. Wu et al. [11] discussed the workface stability of shield tunnels, considering the horizontal arching effect in front of the tunnel face. Hu et al. [12] presented a novel wedge model for the analysis of face stability, taking into consideration the observed failure modes, the character differences in layered soil and the arching effect of soils. Lu

In shield tunnel projects, face stability is always an important issue when tunnel heading. Face stability analysis aims to ensure the safety of the project with regard to soil collapse in front of the tunnel face. Thus, the minimum supporting force (air, slurry, or earth) required for the tunnel face needs to be determined to prevent collapse. As a tunnel is excavated, new rings of the lining behind the shield are installed rapidly, and the collapse of the tunnel face is often sudden due to a sudden loss of the supporting force. Hence, it is reasonable to describe the strength of soils using the undrained shear strength cu . This paper focuses on a case in which the support pressures σT applied on the tunnel face are uniform and the soils are purely cohesive. Additionally, the soils around the tunnel are assumed to be homogeneous, and the undrained shear strength cu is assumed to be constant with depth. Face stability during the excavation of soft clay ground tunnels has been investigated by many contributors to the relevant body of literature. Based on field observation data collected during tests, the socalled stability ratio, N , was proposed by Broms et al. [1], and it is used to analyze the stability of the tunnel face. The stability ratio is defined as N = [σS − σT + γ ·(C + D /2)]/ cu , where σS is the possible surcharge loading acting on the ground surface, σT is the uniform pressure applied on the tunnel face, γ is the soil unit weight, C is the cover depth of the tunnel, D is the diameter of the circular tunnel, and cu is the undrained soil cohesion. Schofield [2] carried out a centrifuge model test and presented a general stability criterion that depended on the tunnel cover. Considering the results of the centrifuge test, Kimura et al. [3] reported that the critical value of N depends on the tunnel depth when N is between 5 and 10. An analytic solution of N was obtained by Ellstein [4] using the limit equilibrium analytical method for

https://doi.org/10.1016/j.compgeo.2018.11.025 Received 1 April 2018; Received in revised form 22 November 2018; Accepted 28 November 2018 Available online 07 December 2018 0266-352X/ © 2018 Elsevier Ltd. All rights reserved.

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

Nomenclature C D C/D γ cu φ N σS

σT qV PT PS Pγ Pq Pe PV H

cover depth of the tunnel diameter of the tunnel relative depth of the tunnel soil unit weight undrained soil cohesion internal friction angle stability ratio possible surcharge loading

uniform support pressure or the critical collapse pressure equivalent uniformly distributed load pressure power of the critical collapse pressure power of the possible surcharge loading power of the soil weight of three zones (I, II and III) power of the uniform distributed load pressure power of the external forces power of internal energy dissipation possible collapse thickness

by Davis et al. [5] show that the failure zone on the tunnel face in cohesive materials covers only part of the circular tunnel face. Leca et al. [14] suggested that the shape of the failure area in frictional materials is an inscribed ellipse with a major semi-axis length equal to half of the diameter of the circular tunnel face. The mechanisms presented by Mollon et al. suggest that the shape of the failure zone on the tunnel face is an entirely circular tunnel face. According to the results of the centrifuge model tests conducted by Idinger et al. [26], Zhang et al. [22] proposed that the failure plane of the tunnel face in cohesivefrictional soils is formed due to shear failure and that its shape is similar to a semi-ellipse with a major semi-axis length equal to the diameter of the circular tunnel face. Thus far, the above studies have shown that the failure zone on the tunnel face or in front of the tunnel face is not clear for cohesive or frictional materials. The purpose of this study is to construct an admissible continuous velocity field using the slip line and limit analysis theories, to analyze the critical collapse face pressure and to find more precise failure zone boundaries in undrained clay. This study focused on purely cohesive soils (e.g. saturated undrained clay), in which the internal frictional angle φ is equal to zero. The results are compared with existing upperbound solutions and illustrated in the form of practical equations that facilitate their application in practice. This study should provide a theoretical basis for the more accurate calculation of the support force acting on the shield tunnel excavation face and the more accurate prediction of the surface influence area on the top of the excavation face for clay stratum to facilitate better tunnel design.

et al. [13] proposed 3D numerical and analytical solutions of limit support pressure at a shield tunnel face by using the elasto-plasticity FEM. This paper aims to analyze the face stability of a tunnel in purely cohesive soil within the framework of a kinematic limit analysis method. At present, one of the most frequently used methods based on limit analysis for analyzing the face stability of a tunnel is the rigid block failure mechanism method. For example, Davis et al. [5] reported an upper-bound solution for a plane strain heading. Subsequently, Leca et al. [14] proposed failure mechanisms based on the assumption of rigid conical blocks and achieved three lower and upper-bound solutions of the critical support pressure for frictional materials. Combining this work with the results of real projects in cohesive and frictional soils, Mollon et al. [15–17] improved these mechanisms and obtained more accurate solutions. More recently, based on numerical simulations and experimental studies, several analysis models of collapse face pressure have been proposed. Zou et al. [18,19] developed a new numerical model by combining the innovative mesh-dividing optimization technology with the 3D rotational failure mechanism to analyze the face stability of tunnels. Chen et al. [20] established a 3D discrete element method (DEM) model to calculate the critical support pressure of shallow shield tunnels. Later, combining this work with an experimental study, Chen et al. [21] carried out a numerical simulation of face stability by means of a 3D finite difference method (FDM) model and proposed a two-stage failure pattern. Considering a series of tunnel diameter-to-depth ratios and corresponding soil properties, Zhang et al. [22] conducted a face stability analysis of a circular tunnel and proposed a new 3D failure criterion that consists of four truncated cones on which a distributed force is imposed. The existing studies show that the rigid block failure mechanisms are a simple and intuitive approach for achieving the critical support pressure of shield tunnels and that they are often assumed to be translational or rotational. The rigid block failure mechanisms are adopted to analyze the face stability of shield tunnels in cohesive and frictional materials. In addition, the shapes of these blocks adhere to the normality condition. Another commonly used method based on limit analysis for analyzing the face stability of a tunnel is to use the continuous velocity field. For example, Mollon et al. [10] suggested that the rigid block failure mechanisms were quite different from the actual velocity field and proposed two new continuous velocity fields for analyzing the collapse pressure of a tunnel face in purely cohesive soil. Later, Zhang et al. [23] proposed lowest upperbound solutions for undrained clays by carrying out an optimization procedure based on the continuous velocity field proposed by Mollon et al. [10]. Huang et al. [24] proposed a new 3D upper-bound method to analyze the undrained clay based on the total loading extended mobilizable strength design(T-EMSD). Considering the linear increase in undrained shear strength with depth, Ukritchon et al. [25] analyzed 3D undrained tunnel face stability in clay. Although these new continuous velocity fields can significantly improve the best existing bounds for the collapse pressure, they are still unable to describe a reasonable failure zone in front of the tunnel face based on these mechanisms. In particular, there are ongoing disputes about the failure area on the tunnel face. For example, the failure mechanisms proposed

2. Slip line and limit analysis theories 2.1. Slip line theory The Coulomb criterion is widely used as the yield condition of soils. In this paper, the Coulomb criterion is simplified as the Tresca criterion because the internal friction angle, φ , of a purely cohesive soil is zero. Taking into account the equations of equilibrium and the yield condition, a group of plastic equilibrium solutions can be obtained for the yield zone. Considering the stress boundary conditions, these solutions can be used to analyze the stress distribution of the yield zone in front of the tunnel face. Therefore, assuming that the stress characteristics are the same as the velocity characteristics for purely cohesive soils, the slip line method can be used to construct the slip line field. In the slip line method, the constitutive equation of the soil is omitted, while the equilibrium and yield conditions are included. Thus, to avoid addressing the plastic equilibrium issue, the theory of limit analysis is used to analyze the collapse face pressures of the tunnel. 2.2. Upper-bound method of limit analysis The theorems on which limit analysis is based are established by assuming that soil is a perfectly plastic material whose flow rule is associated. A considerable amount of test evidence shows that this assumption is feasible for purely cohesive soils. Within the framework of this assumption, the critical collapse pressure of a tunnel face can be 111

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

⌢ ⌢ curved surface FNE and the curved surface FME , as shown in Fig. 5(c). The plane MFNOM is formed by cutting the cone [as shown in the dashed line in Fig. 5(c)] with a vertical plane at point F. The circular plane is formed by cutting the cone with a horizontal plane at point F. The plane MENOM is formed by cutting the cone with a sloping plane (at point E) at an angle of δ + α − π/2 (equal to π/4 ) from the horizontal plane. The curved surfaces are formed by cutting the cone with a vertical plane (at point F), a horizontal plane (at point F) and a sloping plane (at point E). Using the Cartesian coordinate system shown in ⌢ Fig. 5(c), the hyperbola MFN can be expressed as

obtained by using the rigorous limit analysis approach (or the upperbound theorems) for purely cohesive soils. Application of the bound theorems involves the calculation of an admissible stress field, in which case the external loads cannot cause collapse, and then the selection of a mechanism of collapse and an appropriate work rate calculation, in which case the external loads must cause collapse. This paper focuses on analyzing the face stability of a tunnel in purely cohesive soil using the kinematical method of limit analysis. For a purely cohesive soil, the internal friction angle φ is equal to zero. Since the internal friction angle φ is equal to the dilatancy angle due to the associated flow rule, there is no volume change with plastic deformation, and the normality condition is followed. Thus, it is assumed that the actual sliding surface is described by both a velocity characteristic and a stress characteristic.

(z − 3D /4)2 /(D /4)2 − y 2 /(D /4)2 = 1

(1)

where 0 ⩽ z ⩽ D /2 . ⌢ ⌢ The curved surfaces FNE and FME are expressed as

(z − 3D /4)2 = (x + D /4)2 + y 2

3. Admissible continuous velocity field

(2)

where − D /2 ⩽ x ⩽ 0 and − 2 D/2 ⩽ y ⩽ 2 D/2 . The second zone, zone II, is a transitional zone [27]. Due to the constraint of the soil in zone I, zone III and the elastic area in front of zone II, the principal stress axes of the soil in zone II rotate. Hence, zone II is formed by rotating the points on the plane MENOM 90 degrees clockwise around axis MN. The curved surface MBNEM is expressed in a Cartesian coordinate system by the equation

Fig. 1 shows that the Coulomb yield criterion can be simplified to the Tresca yield criterion because the internal friction angle φ is zero for a purely cohesive soil. To satisfy the normality condition and obey the yield criterion, the feasible slip line field needs to be analyzed when establishing the admissible continuous velocity field. For purely cohesive soils, when the material is in the plastic state, four orthogonal shear planes (planes χ , η , ξ and ζ )may exist at each point, as shown in Fig. 2. The possible shear planes are connected to form the so-called slip line network. For purely cohesive soils, it is assumed that the stress characteristics are the same as the velocity characteristics. This paper aims to find the accurate failure zone (failure mechanism) in front of the tunnel face and obtain the critical collapse pressure when tunnel heading, as shown in Fig. 3. Thus, the stress state of the soil in front of the tunnel face needs to be analyzed. To accurately classify the different stress state zones in the soil, that is, the components of the failure mechanism, the following assumptions are made; and the assumed failure mechanism and its components are shown in Figs. 4 and 5:

y 2 = D2 /2 −

2 D x 2 + z 2 /2

(3)

where − 2 D /2 ⩽ x ⩽ 0 and − D/2 ⩽ z ⩽ D/2 , as shown in Fig. 5(b). Meanwhile, zone II obeys the normality condition. The first zone, zone I, is also in the Rankine zone. Since the tunnel face can be regarded as a movable surface when the support force acting on the tunnel face is not enough to maintain the stability of the tunnel face, the direction of the principal stress axes of the soil in zone I remains unchanged [28]. Thus, zone I can be formed by projecting the points on plane MBNOM to the z-coordinate plane through a rotation of 45 degrees in the positive direction of the z-axis; the curved surface MANBM is expressed in the Cartesian coordinate system by the equation

y 2 = D 2 /2 + D (x + z )/2

(a) The friction on the tunnel face caused by the rotation of the cutter head of the shield machine is ignored. (b) The effect of Terzaghi arching in the upper soil on the lining at the rear of the tunnel face is ignored (C / D ⩽ 0.5) for purely cohesive soils. (c) σS / cu + γC / cu ⩾ 2 is satisfied.

(4)

where − D /2 ⩽ x ⩽ 0 and - D ⩽ z ⩽ 0 , as shown in Fig. 5(a). The admissible continuous velocity field with an axially symmetric surface is constructed using the slip line and limit analysis theories, as shown in Fig. 6. This admissible continuous velocity field respects the normality condition and obeys Tresca’s yield criterion.

Fig. 4(a) shows that the failure mechanism consists of four zones. The shape of the failure mechanism can be optimized with respect to the four variable angles, as shown in Fig. 4(b). In the present failure mechanism, since a purely cohesive soil satisfies the Tresca yield criterion and the internal friction angle φ is zero, angles α , β , and θ are equal to π/4 , and angle δ is π/2. Since the geometry of a circular tunnel and the distribution of the possible surcharge loading on the ground surface are symmetrical along the longitudinal axis of the tunnel in the half-space, all failure zones of the present failure mechanism are also symmetrical. Due to 3D equilibrium conditions around the tunnel face, the orientation of the principal stress axes in the soil varies among the failure zones. Considering the yield criterion and the normality condition (e.g., associate flow rule), each failure zone in front of the tunnel face can be determined as follows, as shown in Fig. 5. The fourth zone, zone IV, is in the Rankine zone, a possible failure zone that is affected by the vertical soil arching effect at the top of zone III. The shape of zone IV is a circular truncated cone [see Fig. 5(d)]. The third zone, zone III, is also in the Rankine zone. Based on the second assumption mentioned above, the orientation of the principal stress axes in the soil in zone III remains unchanged. Zone III is formed by three planes and two curved surfaces, e.g., the plane MFNOM, the plane MENOM, the plane created by the circle with diameter EF, the

4. Limit collapse thickness criterion Fig. 4(a) shows that two possible fracture mechanism shapes exist under different values of C / D . For values of C / D greater than 0.5, the

Fig. 1. Tresca’s representation of stress and a yield criterion. 112

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

Fig. 2. Possible slip planes at one point of the Tresca yield criterion.

Fig. 3. 3D graph of tunnel heading.

circular truncated cone (JKEF) with height H may be the fourth possible fracture zone influenced by the vertical soil arching effect, as shown in Fig. 4(a). This fracture model is in good agreement with the experimental phenomena obtained by Schofield [2]. In Fig. 4(a), three critical factors, namely, the weight of the soil in the circular truncated cone (JKEF), the uniform surcharge loading σS on the ground surface, and the undrained soil cohesion cu , determine the stability of the circular truncated cone (JKEF). Thus, to calculate the critical collapse pressure, that is, to determine the limit collapse thickness (C = H + D/2 ), the stability of the circular truncated cone needs to be evaluated. Assuming that the fracture reaches the ground surface, the soil in the circular truncated cone should be loose, and the equivalent pressure acting on the circle with diameter EF can be obtained based on Terzaghi’s soil pressure relation. Fig. 7 shows the modified theory of Terzaghi’s soil pressure relation. Fig. 7 shows that the equivalent uniformly distributed load pressure qV acting on the circle with diameter EF can be given by

qV = γH − 2cu ln(1 + 4

H ) + σS D

Fig. 4. Assumed failure mechanism: (a) all failure zones; (b) failure zones in the plane of symmetry of the tunnel.

(5)

4.1. Analysis of the limit collapse thickness

where H is the possible collapse thickness. D In Fig. 4(a), H = C − 2 , and Eq. (5) then becomes

qV = γD (C / D − 1/2) − 2cu ln(4C / D − 1) + σS

4.1.1. No uniform surcharge loading σS on the ground surface In this case, Eq. (6) reduces to (6)

qV = γD (C / D − 1/2) − 2cu ln(4C / D − 1)

(7)

When qV ⩽ 0 , that is, γD / cu ⩽ 2 ln(4C / D − 1)/(C / D − 1/2) , the soil in the circular truncated cone (JKEF) is stable. Therefore, the range of 113

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

Fig. 5. Components of the assumed failure mechanism: (a) zone I of the failure zones; (b) zone II of the failure zones; (c) zone III of the failure zones; and (d) zone IV of the failure zones.

unstable, and the equivalent uniformly distributed load pressure, qV , must be considered when calculating the critical collapse pressure. Therefore, the limit collapse thickness is equal to C (i.e., H + D/2 ).

the failure mechanism includes zones I, II and III, and the limit collapse thickness is equal to 0.5D ; otherwise, the limit collapse thickness is equal to C (i.e., H + D/2 ). Fig. 8 more intuitively indicates the limit collapse thickness. Fig. 8 shows that, for a certain value of C / D , the soil in the circular truncated cone (JKEF) is stable when the value of γD / cu is located to the left of the curve or along the curve. Thus, the limit collapse thickness is equal to 0.5D . Additionally, when the value of γD / cu is located to the right of the curve, the soil in the circular truncated cone (JKEF) is

4.1.2. Uniform surcharge loading, σS , on the ground surface In this case, for that is, qV ⩽ 0 , γD (C /D − 1/2) − 2cu ln(4C/ D − 1) + σS ⩽ 0 , the soil in the circular truncated cone (JKEF) is stable. Therefore, the range of the failure mechanism includes zones I, II and III, and the limit collapse thickness is equal to 0.5D . Otherwise, the limit collapse thickness is equal to C (i.e., H + D/2 ). To more intuitively determine the limit collapse thickness, a chart describing the relationships among C / D , γD / cu and σS / cu is used, as shown in Fig. 9. The soil in the circular truncated cone (JKEF) is stable when the point with coordinates (C / D , γD / cu , σS / cu ) is located to the left of the curved surface in Fig. 9 or the curved surface described by Eq. (6) and when the limit collapse thickness is equal to 0.5D . Otherwise, the limit collapse thickness is equal to C (i.e., H + D/2 ). 5. Critical collapse pressure The critical collapse pressure of the face of a circular tunnel constructed by a pressurized shield is the main issue in a purely cohesive soil. Based on the kinematical approach of LA (limit analysis) and the slip line theory, a 3D multizone failure mechanism is proposed. As mentioned in the previous sections, the shape of the fracture mechanism depends on the values of C / D . Therefore, for a failure mechanism shape, the critical collapse pressure of the face of a circular

Fig. 5. (continued) 114

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

Fig. 5. (continued)

is obtained. The external force consists of the weight of three zones, the possible uniform surcharge loading on the ground surface and the critical collapse pressure acting on the tunnel face. The internal force mainly refers to the undrained cohesion. The power of the external forces and the internal energy dissipation associated with this fracture mechanism are detailed in the Appendix A. Therefore, the work equation equates the power of the external

tunnel under different values of C / D is given as follows. 5.1 C/D. equal to 0.5 Fig. 10 shows the fracture mechanism for values of C / D equal to 0.5. This fracture mechanism consists of three zones. To calculate the critical collapse pressure, the power of the internal and external forces 115

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

collapse pressure becomes

σT = 0.9706γD − 4.9577cu + σS

(14)

When σT = 0, Eq. (14) becomes

[σS + 0.9706γD]/ cu = 4.9577

(15)

5.2 C/D. greater than 0.5 The fracture mechanism for values of C / D greater than 0.5 is shown in Fig. 4(a). This fracture mechanism consists of four zones. To calculate the critical collapse pressure, the power of the internal and external forces is obtained. The external force consists of the weight of four zones, the possible uniform surcharge loading on the ground surface and the critical collapse pressure acting on the tunnel face. Again, the internal force mainly refers to the undrained cohesion, and the power of the external forces and the internal energy dissipation associated with this fracture mechanism are detailed in the Appendix A. Thus, the work equation equates the power of the external forces to the power of internal energy dissipation, and it is given by Eq. (8). In this case, Eq. (8) can be expressed by

Fig. 6. Admissible continuous velocity field with an axially symmetric surface.

forces Pe to the power of internal energy dissipation PV and it is given as follows:

Pe = PV

Pcu, I + Pcu, II + Pcu, III + Pcu, MNF = Pγ , I + Pγ , II + Pγ , III + Pq + PT

(8)

where Pq is the power of the uniform distributed load pressure, qV . After some simplifications, the tunnel collapse pressure can be expressed as

(9)

σT = γDNγ − cu Nc + σS NS

Eq. (8) can be rewritten as

Pcu, I + Pcu, II + Pcu, III + Pcu, MNF = Pγ , I + Pγ , II + Pγ , III + PS + PT

(10)

where Nγ , Nc , and NS are nondimensional coefficients that represent the effects of soil weight, cohesion, and surcharge loading, respectively. The expressions of coefficients Nγ , Nc , and NS are given as follows:

Nγ = 0.9706

(11)

Nc = 4.9577

(12)

NS = 1.0

(13)

(17)

where Nγ , Nc , and NS are nondimensional coefficients that represent the effects of soil weight, cohesion, and surcharge loading, respectively. In this case, the expressions of coefficients Nγ , Nc , and NS are given as follows:

where Pcu, I , Pcu, II , Pcu, III and Pcu, MNF are the contributions of zone I, zone II, zone III, and plane MNF , respectively; Pγ , I , Pγ , II andPγ , III are the powers of the soil weight of zones I, II, and III, respectively; PS is the power of the possible uniform surcharge loading σS ; and PT is the power of the critical collapse pressure σT . After some simplifications, the tunnel collapse pressure can be expressed as

σT = γDNγ − cu Nc + σS NS

(16)

Nγ = Nγ′ + Nγ″ = 0.9706 + (C / D − 1/2)

(18)

Nc = Nc′ + Nc″ = 4.9577 + 2 ln(4C / D − 1)

(19)

NS = 1

(20)

where Nγ′ = 0.9706 is the effect of the soil weight of zones I, II and III; Nγ″ = (C / D − 1/2) is the possible effect of the soil weight of zone IV; Nc′ = 4.9577 is the effect of cohesion along the surfaces in zones I, II and III; Nc″ = 2 ln(4C / D − 1) is the possible effect of cohesion of the soil in zone IV; and NS = 1.0 is the effect of a possible uniform surcharge loading on the ground surface. Notably, whether the influences of Nγ″, Nc″ and NS need to be considered in Eqs. (18), (19) and (20) depends on dimensionless parameters, such as γD / cu , σS / cu and C / D . Thus, the calculation of the

Substituting Eqs. (11), (12) and (13) into Eq. (10), the critical

Fig. 7. Modified theory of Terzaghi’s soil pressure relation. 116

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

In this case, the equivalent uniform distribution load, qV , is less than or equal to zero, and Eq. (20) becomes

NS = 0

(23)

Substituting Eqs. (21), (22) and (23) into Eq. (17), the critical collapse pressure is given as follows:

σT = 0.9706γD − 4.9577cu

(24)

When σT = 0, Eq. (24) is given by

γD / cu = 5.108

(25)

(b) γD / cu (C / D − 1/2) + σS / cu − 2 ln(4C / D − 1) > 0 In this case, Nγ″, Nc″ and NS must be considered. Substituting Eqs. (18), (19) and (20) into Eq. (17), the critical collapse pressure is given as follows:

σT = [0.9706 + (C / D − 1/2)] γD − [4.9577 + 2 ln(4C / D − 1)] cu + σS

Fig. 8. Design chart for the relationship between γD / cu and C / D .

(26) When σT = 0, Eq. (26) is given by

[σS + γD (C / D + 0.4706)]/ cu = 4.9577 + 2 ln(4C / D − 1)

(27)

6. Analytical comparisons For purely cohesive soils, it is sufficient to analyze the face stability of a circular shield-driven tunnel using six variables, i.e., σT , σS , C , D , cu and γ . Dimensional analysis is often used as a tool for addressing multiparameter problems that arise in practical projects by reducing groups of variables, e.g., σT / cu , σS / cu , C / D and γD / cu . Since undrained behavior is assumed for purely cohesive soils, σT / cu and σS / cu can be replaced by (σS − σT )/ cu . For purely cohesive soils, the method commonly used to evaluate the stability of a circular tunnel face is based on the so-called stability ratio, N , which is defined by Broms et al. [1]. Hence, to validate the results of the present analysis model, the comparisons of the stability ratio, N , between the present model and the existing models, including analytical, numerical and experimental solutions, are analyzed. Because several researchers believe that the contribution of the soil weight to the stability ratio, N , of (C / D + 1/2) is not strict enough, analytical comparisons are made under different values of γD / cu , as discussed below.

Fig. 9. Design chart for the relationships among C / D , γD / cu and σS / cu .

6.1. γD / cu = 0 In this case, the stability ratio, N , becomes

N = (σS − σT )/ cu

(28)

In this study, the stability ratio can be expressed as

N = Nc = 4.9577 + 2 ln(4C / D − 1)

Thus, a comparison of the stability ratio N between the present model and the existing solutions is shown in Fig. 11. Fig. 11 shows that the upper-bound solutions proposed by the present model are between the upper-bound solutions and the lower-bound solutions proposed by Davis et al. [5]. Therefore, the upper-bound solutions proposed by the present model are better than those proposed by Davis et al. [5]. Fig. 11 shows that the upper-bound solutions proposed by the present model are closer to the lower-bound solutions (thick cylinder) proposed by Davis et al. [5] for values of C / D < 0.86, whereas for values of C / D > 0.86, the upper-bound solutions proposed by the present model are closer to the lower-bound solutions (thick sphere) proposed by Davis et al. [5] The maximum difference between the upper-bound solutions proposed by the present model and the lower-bound solutions proposed by Davis et al. [5] is 2.76 for C / D = 1.0 . For C / D > 1.0 , the difference between the upper-bound solutions proposed by the present model and the lower-bound solutions

Fig. 10. Fracture mechanism forC / D = 0.5.

critical collapse pressure for the different values of C / D according to Eq. (6) is discussed below. (a) γD / cu (C / D − 1/2) + σS / cu − 2 ln(4C / D − 1) ⩽ 0 Under this condition, Nγ″ and Nc″ are equal to zero; thus, Eqs. (18) and (19) are reduced to

Nγ = 0.9706

(21)

Nc = 4.9577

(22)

(29)

117

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

values of γD / cu on the stability ratio N is insignificant; hence, the stability ratio N of the present model can be approximated by Eq. (29). For the M2 model proposed by Mollon et al. [10], the maximum influence degree δ is 6.98% when C / D ⩽ 5.0 . Thus, the effect of γD / cu needs to be taken into account when calculating the stability ratio N of the M2 model proposed by Mollon et al. [10]. Fig. 13 shows the analytical comparisons of the stability ratio N between the present model and the existing solutions. Fig. 13 shows that the upper-bound solutions proposed by the present model are lower than those proposed by other models for values of C / D from 0.5 to 3.0 when γD / cu > 0 . The upper-bound solutions proposed by the present model are closer to the analytical solutions, experimental solutions and lower-bound solutions than those proposed by other models. This result illustrates that the upper-bound solutions proposed by the present model are better than those proposed by other models. Additionally, the upper-bound solutions proposed by the present model provide a safer estimate of the tunnel pressure required to maintain stability than those proposed by other models. It can be seen from Fig. 13 that the stability ratioN of the present model increases with increasing values of γD / cu for a certain value of C / D , whereas the stability ratio N of the M2 model proposed by Mollon et al. [10] decreases with increasing values of γD / cu , corresponding to a certain value of C / D . For the actual project in purely cohesive soils with the same parameters cu , C and D , it is obvious that the greater the values of γ are, the more unsafe the actual project. That is, the stability ratio N should be large. Thus, the variation in stability ratio N of the M2 model proposed by Mollon et al. [10], with the values of γD / cu corresponding to a certain value of C / D , may conflict with the actual situation. Meanwhile, Fig. 13 indicates that the effect of the stability ratio of the present model on the values of γD / cu corresponding to a certain value of C / D is insignificant. The reason for this phenomenon is that the three-dimensional equilibrium conditions around the tunnel face can improve the stability of the tunnel face and weaken the effect of the unit weight γ of soils on the stability of the tunnel face. Thus, the stability ratio N defined by Broms et al. [1], is reasonable and safe, neglecting the effect of γD / cu . The initial shape of the failure zone in front of the tunnel face is difficult to observe using experimental methods. Hence, the critical collapse pressure obtained by the centrifuge model test should be greater than that required in an actual project. In this study, the present failure mechanism is constructed based on the slip line and limit analysis theories, and it satisfies the yield criterion and the normality condition (associate flow rule). Thus, this mechanism may accurately explain the progressing failure of the soil in front of the tunnel face. The upper-bound solutions proposed by the present model significantly improve the best existing upper-bound solutions for critical collapse pressure. Fig. 14 illustrates comparisons of the critical coefficients (i.e., Nc , Nγ and NS ) between the existing models and the present model. In Fig. 14,

Fig. 11. Comparisons of N between the present model and the existing solutions

proposed by Davis et al. [5] decreases with increasing C / D . When C / D = 5.0 , this difference is 1.26. Thus, for the present model, the greater the values of C / D are, the safer the tunnel face. Meanwhile, the present model should be closer to the failure mechanism in the actual project. 6.2. γD / cu > 0 For the more general case of a cohesive soil, the traditional method used to assess the critical collapse pressure uses the following equation:

σT = γDNγ − cu Nc + σS NS

(30)

Since the critical coefficient, NS , is equal to 1.0 for purely cohesive soils under undrained conditions, Eq. (30) becomes

σT − σS = γDNγ − cu Nc

(31)

Eq. (31) can be transformed into

[(σS − σT ) + γDNγ ]/ cu = Nc

(32)

Adding γD (C / D + 1/2−Nγ )/ cu to both sides of Eq. (32), it becomes

[(σS − σT ) + γD (C / D + 1/2)]/ cu = Nc + γD [(C / D + 1/2) − Nγ ]/ cu (33) After simplifying Eq. (33), the stability ratio, N , can be expressed as

N = Nc + γD [(C / D + 1/2) − Nγ ]/ cu

(34)

Eq. (34) shows that the stability ratio N is dependent on γD / cu when considering the effect of the self-weight of the soil. To analyze the effect of the self-weight of the soil, the influence degree δ of the self-weight of the soil on the stability ratio N is defined as

δ = {1 + γD [(C / D + 1/2) − Nγ ]/(cu Nc )} × 100%

(35)

Table 1 The stability ratio, N , of the present model in this study.

In this study, the stability ratio N and the influence degree δ corresponding to different values of C / D , can be obtained according to Eqs. (34) and (35), as shown in Tables 1 and 2. Based on the research data of the M2 model proposed by Mollon et al. [10], the stability ratio N and the influence degree δ calculated by Eqs. (34) and (35) are shown in Tables 3 and 4 at certain values of C / D . Table 1 shows that the stability ratio N of the present model increases with increasing C / D and γD / cu . However, Table 3 shows that the stability ratio N of the M2 model increases with increasing C / D , whereas it decreases with increasing γD / cu . Tables 2 and 4 show that the influence degree δ of the self-weight of the soil on the stability ratio N decreases with increasing C / D , whereas it increases with increasing γD / cu . For the present model, the maximum influence degree δ is only 2.97% when C / D ⩽ 5.0 . Meanwhile, Fig. 12 shows that the effect of the

γD/cu N

118

C/D

0.0

1.0

2.0

3.0

4.0

5.0

0.5 0.6 0.8 1.0 1.3 1.5 1.6 2.0 2.5 3.0

4.96 5.63 6.53 7.15 7.83 8.18 8.33 8.85 9.35 9.75

4.99 5.66 6.56 7.18 7.86 8.21 8.36 8.88 9.38 9.78

5.02 5.69 6.59 7.21 7.89 8.24 8.39 8.91 9.41 9.81

5.05 5.72 6.62 7.24 7.92 8.26 8.42 8.94 9.44 9.84

5.08 5.75 6.65 7.27 7.95 8.29 8.45 8.97 9.47 9.87

5.10 5.78 6.68 7.30 7.97 8.32 8.48 9.00 9.50 9.90

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

Table 2 The influence degree of the weight of soil on the stability ratio of the present model in this study. γD/cu δ/% C/D

1.0

2.0

3.0

4.0

5.0

0.5 0.6 0.8 1.0 1.3 1.5 1.6 2.0 2.5 3.0

0.60 0.53 0.45 0.42 0.38 0.36 0.36 0.34 0.32 0.30

1.19 1.04 0.90 0.82 0.75 0.72 0.71 0.66 0.63 0.60

1.78 1.57 1.35 1.23 1.13 1.08 1.06 1.00 0.94 0.90

2.37 2.09 1.80 1.64 1.50 1.44 1.41 1.33 1.26 1.21

2.97 2.61 2.25 2.05 1.88 1.80 1.76 1.66 1.57 1.51

Table 3 The stability ratio N of the M2 model proposed by Mollon et al. [10]. Fig. 12. Curved surface of the relationships among C / D , γD / cu and N .

γD/cu N C/D

0.0

1.0

2.0

3.0

4.0

5.0

0.6 0.8 1.0 1.3 1.6 2.0 2.5 3.0

6.45 7.19 7.87 8.81 9.64 10.64 11.73 12.68

6.36 7.1 7.78 8.71 9.54 10.53 11.61 12.56

6.27 7.01 7.69 8.61 9.44 10.42 11.49 12.44

6.18 6.92 7.6 8.51 9.34 10.31 11.37 12.32

6.09 6.83 7.51 8.41 9.24 10.2 11.25 12.2

6 6.74 7.42 8.31 9.14 10.09 11.13 12.08

Table 4 The influence degree δ of the M2 model proposed by Mollon [10]. γD/cu δ/% C/D

1.0

2.0

3.0

4.0

5.0

0.6 0.8 1.0 1.3 1.6 2.0 2.5 3.0

1.40 1.25 1.14 1.14 1.04 1.03 1.02 0.95

2.79 2.50 2.29 2.27 2.07 2.07 2.05 1.89

4.19 3.76 3.43 3.41 3.11 3.10 3.07 2.84

5.58 5.01 4.57 4.54 4.15 4.14 4.09 3.79

6.98 6.26 5.72 5.68 5.19 5.17 5.12 4.73

Fig. 13. Analytical comparisons of N between the present model and the existing solutions (γD/cu > 0).

stability ratioN is independent of γD / cu when a kinematically permissible mechanism is used and that a safe upper-bound solution to assess the critical collapse pressure required to maintain a stable tunnel face is obtained.

the critical coefficient NS (i.e., the effect of surcharge loading) of each model is equal to 1.0, which shows that the effect of surcharge loading on the critical collapse pressure does not vary with C / D . Since these methods for cohesive materials involve no volume change, the decrease in area of the tunnel must equal the area of ground loss at the surface. Fig. 14(a) shows that the order, from high to low, of the critical coefficients corresponding to a certain value of C / D is Nc , Nγ and NS when 0.5 ⩽ C / D ⩽ 3.0 , indicating that the effect of cohesion on the critical collapse pressure is the most remarkable, followed by the effect of the soil weight; the weakest effect is due to surcharge loading. Fig. 14(b) shows that the critical coefficient Nγ of all models linearly increases with C / D . The critical coefficient Nγ of these models, excluding the present model and the M2 model proposed by Mollon et al.[10], is fixed (C / D + 1/2 ). It can be seen from Fig. 14(b) that the increase in the critical coefficient Nγ with respect to Nγ = C / D + 1/2 when using the M2 model proposed by Mollon et al. [10] is smaller than 6.98% when C / D ⩾ 0.5. In addition, the decrease in the critical coefficient Nγ with respect to Nγ = C / D + 1/2 when using the present model is smaller than 2.97% when C / D ⩾ 0.5. Hence, it is reasonable to assume that the

7. Case analysis In this case, the geometrical and material parameters of the circular tunnel example in purely cohesive soil include the following: diameterD = 10 m, soil unit weight γ = 18 kN/m3, and undrained soil cohesioncu = 20 kPa or 30 kPa. Hence, the two corresponding dimensionless parameters γD / cu are 9.0 and 6.0, respectively. Meanwhile, the uniform surcharge loading σS on the ground surface is not considered in this case. To calculate the critical collapse pressure, a suitable fracture mechanism that depends on C / D and the corresponding equations need to be determined. ForC / D = 0.5, the fracture mechanism of the two cu values is shown in Fig. 10 (σS =0). The critical collapse pressures σT calculated using Eq. (14) withσS = 0 are 75.55 kPa and 25.98 kPa. When C / D > 0.5, whether the soil in zone IV is stable needs to be assessed using Eq. (7) for certain values of γD / cu . When γD / cu is 9.0, the soil in zone IV is unstable, and the equivalent uniform distribution load qV is greater than zero; the fracture mechanism is shown in Fig. 4(a) (whereσS = 0). The critical collapse pressure σT can 119

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

zone in front of the tunnel face has a clear physical meaning and satisfies the abovementioned conditions. Hence, the present model should significantly improve the best existing upper-bound solutions for the collapse pressure. Additionally, the present model can effectively explain the failure mechanism of the shield-driven tunnel face. Comparisons of the sizes of the failure zones of various models in a half-space with various values of C / D are shown in Fig. 16. Fig. 16 shows that the upper boundary of the failure zones of the present model appears above, behind and in front of the tunnel face, whereas those obtained by Davis et al. [5] and Mollon et al. [10] is only in front of the tunnel face. Additionally, Fig. 16 shows that those of Mollon et al. [10] suggest that the entire circular tunnel face is failing, and those of Davis et al. [5] suggest that part of an ellipse inscribed in the circular area is failing. However, in the present model, the area of failure at the circular tunnel face consists of two parts: the upper half of the circular area, and the area above the intersection of the lower half of the circular tunnel face and plane MANOM (as shown in Fig. 22). The shape of the failure zones at the tunnel face obtained by the present model agrees well with the results of the centrifuge model tests performed by Idinger et al. [26] and the results of the numerical simulations proposed by Zhang et al. [22]. These results illustrate that the present model is closer to the failure mode observed in actual engineering than that of other models. Because it considers the yield criterion and the normality condition (e.g., associate flow rule), the present model can effectively explain the progressive failure of the shield-driven tunnel face. Meanwhile, more accurate failure zone

(a) Nc-N -Ns

(b) N -Ns Fig. 14. Comparison of Nc, Nγ and Ns between the existing 3D models and the present model: (a) Nc-Nγ-Ns; (b) Nγ-Ns.

be calculated using Eq. (26) withσS = 0. WhenγD / cu = 6.0, the soil in zone IV is stable for values of C / D from 0.5 to 0.6834, and the equivalent uniform distribution load qV is equal to zero. The fracture mechanism is shown in Fig. 10 (σS =0). The critical collapse pressure σT calculated using Eq. (14) withσS = 0 is 25.98 kPa [see Fig. 15(b)]. ForC / D > 0.6834, the soil in zone IV is unstable forγD / cu = 6.0, and the equivalent uniform distribution load qV is greater than zero. The fracture mechanism is shown in Fig. 4(a) (σS =0). Therefore, the critical collapse pressure σT should be calculated using Eq. (26) with σS = 0. For this case, the solutions of the critical collapse pressure obtained by Davis et al. [5], Klar et al. [9] and Mollon et al. [10,15–17] and the present approach are given in Fig. 15(a) and (b). In addition to the numerical solution, note that the other results are based on the kinematical approach of LA. Fig. 15 shows that the solutions of the M2 model and the present model are the two solutions that fit the numerical model best. In Fig. 15(a), when γD / cu is 9.0, the values of σT provided by M2 are very satisfying for C/D ⩽ 1.5, whereas the values of σT provided by the present model are very satisfying for C/D > 1.5. Fig. 15(b) shows that the values of σT provided by the present model are very satisfying for values of C/D from 0.5 to 3.0 when γD / cu = 6.0. Since the failure shape in front of the tunnel face from the M2 model was roughly constructed as a simple geometrical object that envelopes the plastically deformed region of the soil mass observed in numerical simulations, the M2 model is inaccurate and includes many uncertainties. In the present model, the yield criterion, the normality condition (e.g., associate flow rule), and the vertical soil arching effect are considered. Each failure

(a) cu = 20kPa

(b) cu = 30kPa Fig. 15. Comparisons of collapse pressures from various models. 120

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

zones in front of the tunnel face in the present model are closer to those observed in actual projects than those of other models. Moreover, we provided comparisons of the results from the present model and other approaches. The main conclusions are as follows: 1. The present failure mechanism is a 4-variable multizone mechanism that consists of four zones, e.g., zones I, II, and III and a possible zone IV. Since the geometry of the tunnel and the distribution of the possible surcharge loading on the ground surface are symmetric along the longitudinal axis of the circular tunnel, the shapes of these four zones are also symmetric. While zones I, III, and IV are in the Rankine zone, zone II is a transitional zone. Whether zone IV is a failure zone needs to be determined by assessing the stability of the soil in zone IV using Terzaghi’s soil pressure relation. 2. The proposed failure mechanism indicates that the critical values of the stability ratio N vary with the burial depth (i.e., C / D ) of the tunnel and the weight (i.e., γ D/ cu ) of the soil. However, Figs. 12 and 13 show that the variation in the critical values of the stability ratio, N , caused by the weight (i.e., γ D/ cu ) of soil is insignificant. Thus, the tunnel stability estimate is reasonable when using the stability ratio N . The comparisons of the collapse pressures between the present model, numerical simulations and other theoretical approaches show that a good agreement was attained with the results of the numerical simulation, especially for values ofC / D > 1.5. 3. The proposed failure mechanism suggests that the failure zones at the tunnel face consist of the upper half of the circular tunnel face and an inscribed parabola on the lower half of the tunnel face instead of the entire circular tunnel face or an inscribed ellipse on the tunnel face, as shown in Fig. 22. This result agrees well with the results of Idinger et al. [26] and Zhang et al. [22]. Moreover, this result demonstrates that the failure zones at the tunnel face are primarily formed by shear failure. 4. Comparisons of the critical values of Nγ , Nc , and NS between the present model and other approaches are given in this paper. The present model significantly improves the critical values of Nc compared with other limit analysis mechanisms. The critical values of Nc were improved by using the proposed failure mechanism, showing good agreement with the results of the analytical [4] and experimental [3] approaches. Although the critical values of Nγ are improved, the difference is insignificant. Moreover, Fig. 14 indicates that the critical values of NS from all approaches are independent of the values of C/D. When 0.5 ⩽ C / D ⩽ 3.0 , the contribution of the critical values of Nc to the collapse pressure is more significant than those of the other two critical values. 5. Equations to express the external envelope surface of the failure zones are proposed, allowing the precise boundary of the failure zones to be plotted. These analytical equations may also be used to calculate the critical collapse pressure required to maintain the stability of a tunnel face. Since the yield criterion, the normality condition and the vertical soil arching effect are considered in the proposed failure mechanism, the failure zones obtained using the present model may be very similar to those in actual projects.

(a) C/D =0.5

(b) C/D =1.0 Fig. 16. Comparisons of the failure zones of various models in a half-space: C/ D = 0.5; (b) C/D = 1.0.

volumes can be obtained using the present model. Hence, the proposed failure mechanism significantly improves the existing upper-bound solutions for the face stability of circular tunnels in purely cohesive soils. 8. Discussion The shape of the failure zone at the tunnel face has been described from various points of view. For cohesive or frictional material, it has been reported in the literature [5,14] that the shape of the failure zone at the tunnel face is an ellipse inscribed on the circular tunnel face. This failure zone has a major semi-axis length equal to D /2, where D is the diameter of the circular tunnel. Subsequently, based on the results observed in their experimental test [29] and numerical simulation [30] work, Mollon et al. [10,15–17] proposed that the shape of the failure zone at the tunnel face was the entire circular tunnel face. However, increasingly more evidence from centrifuge model tests [26,31], and numerical simulations [22,32,33] shows that the envelope of the tunnel face failure should be bulb-shaped, indicating that the failure plane at the tunnel face is primarily formed by shear failure. In this paper, the shape of the failure zone at the tunnel face is in agreement with the latter point of view. 9. Conclusion A new 3D multizone translational failure mechanism based on the slip line theory and kinematical approach of LA was proposed in this study to improve the existing solutions of the critical collapse pressure of a circular shield-driven tunnel. Compared with the existing limit analysis mechanisms, the present model represents a significant advantage in terms of the accurate determination of the critical collapse pressure. Since this model takes into account the yield criterion, the normality condition and the vertical soil arching effect, the failure

Acknowledgements The authors acknowledge the financial support provided by the National Natural Science Foundation of China (Grant No. 41572275) and the Natural Science Foundation of Shandong Province (Grant No. ZR2012EEM006).

121

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

Appendix A. Derivation of the external power and dissipation power A.1 Geometric properties This study considered two types of fracture mechanisms corresponding to different values of C / D . For C / D equal to 0.5, the first fracture mechanism consists of three zones, i.e., zones I, II and III (see Fig. 10). For C / D greater than 0.5, the second fracture mechanism involves four zones, i.e., zones I, II, III and possibly IV [see Fig. 4(a)]. The geometric properties of three zones (i.e., zones I, II and III) of the two fracture mechanisms are the same. Because of the contribution to the collapse loads of the surcharge loading on the ground surface of the top of the possible zone IV, the soil weight and the undrained soil cohesion of the possible zone IV can be equivalent to the uniformly distributed load pressure, qV , acting on the circle with diameter EF based on Terzaghi’s soil pressure relation (see Fig. 7). The geometric properties of the possible zone IV can be disregarded. Thus, the geometric quantities of zones I, II and III need to be determined to deduce the power of the external force and the internal dissipation energy. The geometric properties include the volumes of the zones (i.e., zones I, II and III), the area of the discontinuous surfaces (i.e., the curved surfaces of zones I, II and III; and the plane MFN) and the area of failure on the tunnel face, which are calculated as follows. A.2 Volumes of zones A.1.1 Volume of Zone I Zone I is shown in Fig. 5(a). The curved surface of zone I can be expressed as Eq. (4) in the Cartesian coordinate system. The space domain Ω1, which is enclosed by the curved surface MANBM, and the planes MAN and MBN in zone I, consist of the space domains Ω11 and Ω12 .

Ω11 = {(z , x , y )| − D /2 ⩽ z ⩽ 0, z ⩽ x ⩽ 0, −

D2 D (x + z ) + ⩽y⩽ 2 2

Ω12 = {(z , x , y )| − D ⩽ z ⩽ −D /2, −D − z ⩽ x ⩽ 0, −

D2 D (x + z ) + } 2 2

D2 D (x + z ) + ⩽y⩽ 2 2

D2 D (x + z ) + } 2 2

(36)

(37)

where Ω1 = Ω11 ∪ Ω12 . Thus, the volume of zone I is given by

V1 =

∭ dV = ∭ dV + ∭ dV = Ω1

Ω11

Ω12

2 2 3 D 15

(38)

A.1.2 Volume of Zone II Zone II is shown in Fig. 5(b). The curved surface of zone II can be expressed as Eq. (3) in the Cartesian coordinate system. Since the planes xoy and xoz are the symmetric surface of zone II, one quarter of the volume of zone II is calculated, and the local cylindrical coordinate system ( ρ , θ , y ) is used, as shown in Fig. 17. The space domain Ω2 occupied by one quarter of the volume of the zone II can be expressed as

Ω2 = {(ρ , θ , y )| 0 ⩽ ρ ⩽

2 D /2, 0 ⩽ θ ⩽ π /4, 0 ⩽ y ⩽

D 2 /2 −

2 Dρ /2 }

(39)

Therefore, the volume of zone II is given by

V2 = 4

∭ dV = 4 ∭ ρdρdθdy = Ω2

Ω2

2π 3 D 15

(40)

A.1.3 Volume of Zone III Zone III is shown in Fig. 5(c). The curved surface of zone III can be expressed as Eq. (2) in the Cartesian coordinate system. Thus, the space

Fig. 17. Local cylindrical coordinate system in zone II. 122

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

domain Ω3 can be defined as

Ω3 = \{ (z , x , y )| 0 ⩽ z ⩽ D /2, −z ⩽ x ⩽ 0, -

(z − 3D /4)2 − (x + D /4)2 y ⩽

(z − 3D /4)2 − (x + D /4)2 \}

(41)

Therefore, the volume of zone III is given by

V3 =

∭ dxdzdy = 0.0972D3 (42)

Ω3

A.3 Areas of the discontinuous surfaces A.3.1 Area of the curved surface of zone I The curved surface of zone I is shown in Fig. 18 (gray shaded area). The curved surface in the Cartesian coordinate system can be expressed as Eq. (4). Thus, the area of the curved surface S1 of zone I is given by

S1 =

∬

1 + D /16(D + x + z ) dA = 0.6323D 2 (43)

Γ1

where Γ1 is the projection domain of the curved surface in the coordinate plane xoz . A.3.2 Area of the curved surface of zone II The curved surface of zone II is shown in Fig. 19 (gray shaded area). The curved surface in the Cartesian coordinate system can be expressed as Eq. (3). Given the symmetry of the curved surface, the area of the curved surface S2 of zone II is calculated as follows:

S2 = 4

∬ Γ2

1+

4D (4 − 2 2 ) π 2 dA = D 3 D − 2ρ

where Γ2 = {(ρ , θ)| 0 ⩽ ρ ⩽

(44)

2 D /2, 0 ⩽ θ ⩽ π /4} .

A.3.3 Area of the curved surface of zone III The curved surface of zone III is shown in Fig. 20 (gray shaded area). The curved surface in the Cartesian coordinate system can be expressed as Eq. (2). Given the symmetry of the curved surface, the area of the curved surface S3 of zone III is calculated as follows:

S3 = 2

∬ Γ3

2a2 32 − 3 2 π 2 dA = D a2 − b2 48

(45)

where Γ3 = {(a, b)| − b − D /2 ⩽ a ⩽ D /4, −D /4 ⩽ b ⩽ D /4} , a = z − 3D /4 , and b = x + D/4 . A.3.4 Area of the plane MFN The plane MFN is shown in Fig. 21. The hyperbola MFN in the Cartesian coordinate system can be expressed by Eq. (1). Thus, the area of the plane MFN is given by

SMNF = [6 2 − ln(3 + 2 2 )] D 2 /16

(46)

Fig. 18. Curved surface of zone I. 123

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

Fig. 19. The whole curved surface and 1/4 of the curved surface of zone II.

A.4 Area of failure on the tunnel face The area of failure on the tunnel face is shown in Fig. 22. The parabola MAN in the Cartesian coordinate system can be expressed by

y 2 = D2 /2 + Dz /2

(47)

Thus, the area of failure on the tunnel face is given by

A = (1/3 + π /8) D 2

(48)

A.5 Velocity field Based on the slip line theory and the kinematical approach of limit analysis, a kinematically admissible velocity field with an axially symmetric surface is constructed, as shown in Fig. 23. Since the failure mechanism is assumed to be translational, the magnitudes of the velocity of the soils in the failure zones are the same. The soil in failure zone IV is assumed to follow Terzaghi’s soil pressure relation. Thus, the magnitude of the kinematic velocity in the failure zones can be expressed as follows:

v1 = v 2 = v 3 = v

(49)

The magnitudes of the kinematic velocity at points B, E and F are given by (50)

v1 = vB = vE = vF = v The magnitude of the kinematic velocity of the tunnel face can be written as

vT =

2 v1/2 =

(51)

2 v /2

The dissipation of plastic energy occurs only along discontinuous surfaces, and the dissipation energy per unit area dPv / d Σ is (52)

dPv / d Σ = cu v∙n

where n is the unit vector normal to the discontinuous surface at the point where dPv / d Σ is computed. Since the angles between the velocity v and the discontinuity surfaces S1, S2 and S3 are zero, Eq. (52) can be written as (53)

dPv / d Σ = cu v Since the angle between the velocity v and the discontinuity surface MNF is π/4 , Eq. (52) can be written as

dPv / d Σ =

(54)

2 cu v /2

Fig. 20. Curved surface of zone III. 124

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

Fig. 21. Plane MFN.

Fig. 22. Area of failure on the tunnel face.

A.6 Power of the external forces A.6.1 C / D equal to 0.5 In this case, the power of the external forces Pe consists of three components

Pe = PT + PS + Pγ

(55)

where PT is the power of the critical collapse pressure σT ; PS is the power of the possible uniform surcharge loading σS ; and Pγ is the power of the soil weight of three zones, i.e., zones I, II and III.

Fig. 23. Kinematically admissible velocity field with an axially symmetric surface. 125

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

PT = −

∬Σ

T

σT vT d Σ = −σT vT A

(56)

Substituting expressions (48) and (51) into (56), we find:

PT = − 2 (1/3 + π /8) σT D 2v /2

(57)

The power of a possible uniform surcharge loading, σS , on the ground surface can be expressed as follows:

PS =

∬Σ

s

σS vns dΣ

(58)

vns

is the downward normal velocity of the surface and Σs is the deforming area of the surface. where In this study, since undrained behavior is assumed, the soil deforms at a constant volume [6] and

∬Σ

s

vns dA =

∬Σ

T

vT dA

(59)

Substituting Eq. (59) into (58), the power of a possible uniform surcharge loading σS on the ground surface becomes

∬Σ

PS = σS

T

vT dA= 2 (1/3 + π /8) σS D 2v /2

(60)

The power of the soil weight of the three zones can be expressed as follows:

Pγ = Pγ , I + Pγ , II + Pγ , III =

∭V

1

γvn1 dV +

∭V

2

γvn2 dV +

∭V

3

γvn3 dV

(61)

where vn1, vn2 and vn3 are the downward normal velocities of the centroid of zones I, II and III, respectively, as shown in Fig. 5(a), (b) and (c); and V1, V2 and V3are the volumes of zones I, II and III, respectively. The magnitudes of the downward normal velocity vn1, vn2 and vn3 are given as follows:

vn1 =

vn2

2 v /2

(62)

=v

vn3 =

(63)

2 v /2

(64)

Substituting Eqs. (38), (40), (42), (62), (63) and (64) into (61), Eq. (61) can be rewritten:

Pγ =

2 2π 3 γD3v + γD v + 9.7218 × 10−2 2 γD3v /2 = 0.4983γD3v 15 15

(65)

Substituting Eqs. (57), (60) and (65) into (55), Eq. (55) can then be expressed as

Pe = 0.5134(σS − σT ) D 2v + 0.4983γD3v

(66)

A.6.2 C / D greater than 0.5 In this case, the fracture mechanism consists of four zones, as shown in Fig. 4(a). In the fracture mechanism, zone IV is a possible fracture zone. The stability of the fourth zone depends on three factors: the weight of zone IV, the uniform surcharge loadings on the ground surface and the undrained cohesion of the soil in zone IV. To obtain the power of these factors of zone IV, Terzaghi’s soil pressure relation is used to calculate the uniform distributed load pressure qV , which is caused by these factors. The pressure acts on a circle with diameter EF, as shown in Fig. 7. In addition, the uniform distributed load pressure qV can be expressed as Eq. (6). Thus, the power of the external force Pe in the fracture mechanism should be given by

Pe = PT + Pq + Pγ

(67)

where PT and PS are the same as those in Eq. (55) and Pq is the power of the uniform distributed load pressure qV . The method for determining Pq is the same as the method for determining PS in Eq. (60).

Pq = qV

∬Σ

T

vT dA= 2 (1/3 + π /8) qV D 2v /2

(68)

Substituting Eq. (6) into Eq. (68), we obtain

Pq =

2 (1/3 + π /8)[γD (

C 1 4C − ) − 2cu ln( − 1) + σs ] D 2v /2 D 2 D

(69)

Substituting Eqs. (57), (65) and (69) into (67), Eq. (67) can then be written as

Pe = 0.5134[σS + γD (C /D − 1/2) − 2cu ln(4C / D − 1) − σT ] D 2v + 0.4983γD3v

(70)

A.7 Dissipation power Whether or not the values of C / D are equal to zero or greater than zero, the dissipation energy can be written as

Pv = Pcu, I + Pcu, II + Pcu, III + Pcu, MNF

(71)

where Pcu, I , Pcu, II , Pcu, III and Pcu, MNF are the contributions of zone I, zone II, zone III, and plane MNF , respectively.

Pcu, I = cu vS1

(72)

Pcu, II = cu vS2

(73) 126

Computers and Geotechnics 107 (2019) 110–127

W. Ding et al.

Pcu, III = cu vS3

(74)

Pcu, MNF =

(75)

2 cu vSMNF /2

Eqs. (43), (44), (45) and (46) allow one to write (72), (73), (74) and (75) in the form:

Pcu, I = 0.6323cu vD 2

(76)

Pcu, II =

(4 − 2 2 ) π cu vD 2 3

(77)

Pcu, III =

32 − 3 2 π cu vD 2 48

(78)

Pcu, MNF =

[12 −

2 ln(3 + 2 2 )] cu vD 2 32

(79)

Substituting Eqs. (76), (77), (78) and (79) for Pcu, I , Pcu, II , Pcu, III and Pcu, MNF in Eq. (71) leads to

Pv = 2.5453cu vD 2

(80)

2011;35(12):1263–388. [18] Zou JF, Ze HQ. Face stability analysis of tunnels excavated below groundwater considering coupled flow deformation. Int J Geomech, ASCE 2018;18(8). https:// doi.org/10.1061/(ASCE)GM.1943-5622.0001199. [19] Zou JF, Chen KC, Pan QJ. An improved numerical approach in surrounding rock incorporating rockbolt effectiveness and seepage force. Acta Geotech 2018;13(3):707–27. [20] 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. [21] Chen RP, Li J, Kong LG, Tang LJ. Experimental study on face instability of shield tunnel in sand. Tunn Undergr Space Technol 2013;33(1):12–21. [22] Zhang CP, Han KH, Zhang DL. Face stability analysis of shallow circular tunnels in cohesive-frictional soils. Tunn Undergr Space Technol 2015;50:345–57. [23] 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. [24] Huang MS, Li S, Yu J, Tan Wen JQ. Continuous field based upper bound analysis for three-dimensional tunnel face stability in undrained clay. Comput Geotech 2018;94:207–13. [25] Boonchai K, Kongkit Y, Suraparb K. Three-dimensional undrained tunnel face stability in clay with a linearly increasing shear strength with depth. Comput Geotech 2017;88:146–51. [26] Idinger G, Aklik P, Wu W, Borja RI. Centrifuge model test on the face stability of shallow tunnel. Acta Geotech 2011;6(2):105–17. [27] Chen WF. Limit analysis and soil plasticity. Elsevier Press; 1975. [28] Chevalier B, Takano D, Otani J. Comparison of X-ray CT and discrete element method in the evaluation of tunnel face failure, GeoX2010: Advances in Computed Tomography for Geomaterials. John Wiley & Sons Inc; 2013. p. 397–405. [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. Proc Geocongr, Proc ASCE, Atlanta 2006;187:1–6. [30] Eisentein AR, Ezzeldine O. The role of face pressure for shieldswith positive ground control. Tunneling and Ground conditions. Balkema: Rotterdam 1994:557–71. [31] Chambon P, Corte JF. Shallow tunnels in cohesionless soil: stability of tunnel face. J Geotech Eng 1994;120(7):1148–65. [32] 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:187–95. [33] Wang J, He C, Wang C, Chen ZQ, Tang R. Face stability analysis of EPB shield tunnel in sand. Chinese J Geotech Eng 2018;40(1):177–84.

References [1] Broms BB, Bennermark H. Stability of clay at vertical openings. J Soil Mech Found Div 1967;SM1(93):71–94. [2] Schofield AN. Cambridge geotechnical centrifuge operations. Geotechnique 1980;30(3):227–68. [3] Kimura T, Mair RJ. Centrifugal testing of model tunnels in clay. Proceedings of the 10th international conference of soil mechanics and foundation engineering, Stockholm, Rotterdam: Balkema, V1. 1981. p. 319–22. [4] Ellstein AR. Heading failure of lined tunnels in soft soils. Tunnels and Tunnelling 1986;18(6):51–4. [5] 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. [6] Sloan SW, Assadi A. Undrained stability of a plane strain heading. Can Geotech J 1994:443–50. [7] Augarde CE, Lyamin AV, Sloan SW. Stability of an undrained plane strain heading revisited. Comput Geotech 2003;30(5):419–30. [8] Lee CJ, Wu BR, Chen HT, Chiang KH. Tunneling stability and arching effects during tunneling in soft clayey soil. Tunn Undergr Space Technol 2006;21(2):119–32. [9] Klar A, Osman AS, Bolton M. 2D and 3D upper bound solutions for tunnel excavation using ‘elastic’ flow fields. Int J Numer Anal Meth Geomech 2007;31(12):1367–74. [10] 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 Anal Meth Geomech 2013;37(13):2061–83. [11] Wu J, Liao SM, Shi ZH. Workface stability of shield tunnel considering arching effect. J Tongji Univ (Nat Sci) 2015;43(2):213–20. [12] Hu XY, Zhang ZX. Calculation model of shield tunneling stability based on real face failure models. J Shanghai Jiao Tong Univ 2013;47(9). 1469–1456. [13] Lu XL, Li FD, Huang MS, Wan JL. Three-dimensional numerical and analytical solutions of limit support pressure at shield tunnel face. J Tongji Univ (Nat Sci) 2012:1469–73. [14] 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. [15] Mollon G, Dias D, Soubra AH. Probabilistic analysis and design of circular tunnels against face stability. Int J Geomech, ASCE 2009;9(6):237–49. [16] Mollon G, Dias D, Soubra AH. Face stability analysis of circular tunnels driven by a pressurized shield. J Geotech Geoenviron Eng, ASCE 2010;136(1):215–29. [17] Mollon G, Dias D, Soubra AH. Rotational failure mechanisms for the face stability analysis of tunnels driven by pressurized shields. Int J Numer Anal Meth Geomech

127