- Email: [email protected]
Continuous field based upper bound analysis for three-dimensional tunnel face stability in undrained clay
Computers and Geotechnics xxx (xxxx) xxx–xxx
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Technical Communication
Continuous field based upper bound analysis for three-dimensional tunnel face stability in undrained clay ⁎
Maosong Huanga, Sen Lia, Jian Yua, , Jorgin Qi Wen Tanb a b
Department of Geotechnical Engineering, Tongji University, Shanghai 20092, China Geotechnical and Tunneling, Land Transport Authority, 219428, Singapore
A R T I C L E I N F O
A B S T R A C T
Keywords: Undrained clay Tunnel face stability Upper-bound limit analysis Finite element Continuous field
A new numerical upper-bound method applicable to three-dimensional stability problems in undrained clay is proposed and applied in the analysis of tunnel face stability. Comparisons with the elasto-plastic finite element method (FEM) reveal that this method could provide sufficiently accurate upper-bound prediction for the tunnel face stability factors. For tunnels with large cover depth, the obtained results greatly improve the existing upperbound solutions. Their corresponding velocity fields are investigated, which show a similar evolving pattern with those of the FEM. A discussion is made to compare the various advantages of multi-block mechanisms and continuous mechanisms in tunnel face stability.
1. Introduction Underground tunnel excavation has become commonplace in infrastructure projects worldwide. An important design parameter in tunnel face stability is the support pressure that is applied at the tunnel face. Due to the practical significance of three-dimensional (3D) tunnel face stability, substantial research investigations have been conducted with various methods including the empirical method [2], lower bound limit analyses [3,8,13], upper-bound limit analyses [3,5,6,8–11,13] and elasto-plastic finite element (FE) analyses [14,16]. Broms and Bennermark [2] first proposed a hand calculation design expression for tunnel face stability in undrained clay by empirically relating the expected deformations with various factors such as support pressure, surcharge, soil unit weight, cover depth and tunnel diameter to determine whether the tunnel is sufficiently supported. The first theoretically rigorous analysis of 3D tunnel face stability was carried out by Davis et al. [3] using both upper and lower bound limit analyses. Leca and Dormieux [8] extended the work of [3] into fictional soil. Ukritchon et al. [14] and Zhou [16] further utilized the finite element method (FEM) to reevaluate the accuracy of the previously published results. As a widely used analysis method, the applications of the upperbound method in 3D tunnel face stability are further reviewed. The upper-bound method can be divided into three categories according to their underlying mechanisms. The first approach is to postulate mechanisms composed of either translational [3,8,10] or rotational rigid blocks [11] where the energy dissipation exists solely at the interfaces between adjacent blocks. The second approach is to develop continuous
⁎
failure mechanisms with no discontinuity, where the energy is solely dissipated by plastic deformation throughout the continuums [5,6,9]. The third approach is to construct failure mechanisms with both deformable zones and discontinuity surfaces [13]. The continuous velocity field of Klar et al. [6] was derived from both the elasticity solution corresponding to the original plastic problem and the concept of sinks and sources. The continuous mechanism of Mollon et al. [9] was obtained from the elasto-plastic finite difference method (FDM). The finite element limit analysis (FELA) was employed by Sloan [13] to investigate 3D tunnel face stability. However, the difference between the upper and lower bound for 3D tunnel problem is quite large when compared with the 2D FELA solutions [1]. Much effort has been put into constructing upper-bound solutions for 3D undrained tunnel stability. However, a significant gap still exists between the upper-bound results and the elasto-plastic FEM results for tunnels with large cover depth. Such a gap might be attributed to the difficulty in constructing complicated 3D upper-bound velocity fields. This paper attempts to propose a generally applicable 3D upper-bound method for undrained clay based on the total loading extended mobilisable strength design (T-EMSD) proposed by Huang et al. [4] and Yu et al. [15], which can acquire the upper-bound solution through elastic iterative analysis and thus could be conveniently replicated through an elastic displacement-based FEM package. The method will be applied to 3D undrained tunnel stability, and the obtained upper-bound solutions will be thoroughly compared with existing solutions.
Corresponding author. E-mail address: [email protected] (J. Yu).
http://dx.doi.org/10.1016/j.compgeo.2017.09.014 Received 27 June 2017; Received in revised form 30 August 2017; Accepted 19 September 2017 0266-352X/ © 2017 Elsevier Ltd. All rights reserved.
Please cite this article as: Huang, M., Computers and Geotechnics (2017), http://dx.doi.org/10.1016/j.compgeo.2017.09.014
Computers and Geotechnics xxx (xxxx) xxx–xxx
M. Huang et al.
2.1. Theory Shield and Drucker [12] demonstrated that the upper-bound theorem for Tresca material, conventionally adopted for undrained clay, can be expressed as
∫s Ti vi dS + ∫V
2. Method of analysis
∫V
The T-EMSD, developed from the EMSD, is a total loading method to obtain an upper-bound load-displacement solution of a plane-strain problem through iterative computation, rather than optimal computation. It has been verified in the study of the soil resistance to the lateral movement of a 2D circular pile. The T-EMSD is able to directly calculate an upper-bound lateral resistance under a given total pile displacement (without the need for an incremental loading process). Large enough displacements would enable the resistance to approach an upper-bound lateral limit load. This section will extend the T-EMSD theory to a generic 3D condition and illustrate its application in obtaining the upper-bound tunnel face stability factor.
∫V
Replace the old CPD with the current pressure difference s ; t Impose a new surchage s and support pressure t , the difference of which equals to half of CPD+NPD
Yes
uk
1 2
allowance value 1
Yes
Note: NPD=non-convergent pressure difference
NPD CPD NPD
allowance value 2
∫S
∫V
Fig. 2. Flow diagram for tunnel stability.
Perform elastic FE analysis with the equivalent shear modulus G d 1 d 3 , iteration number k=1
uk
fi vi dV (i = 1,2,3)
(2)
(3)
where dεmax = max(|dε1 |,|dε3 |,|dε3 |) is the increment of the largest absolute principal strain; cm (εs ) is the mobilised strength, which is a function n of engineering shear strain εs (= ∑1 |dε1−dε3 |, n is the number of the loading steps); for an elastic perfectly-plastic material cm (εs ) equals to Min (Gs εs,su ) , where Gs is the shear modulus of the soil. With εs accumulating during loading, cm (εs ) tends to su such that Ti∗ tends to an
Impose a surcharge s and support pressure t on the ground surface and tunnel face; record 0 as CPD and a sufficiently large value as NPD
NO
∫S Ti∗ vi dS + ∫V
Internal increment of work External increment of work 2cm (εs) dεmax dV = Ti∗ dui dS + fi dui dV (i = 1,2,3)
Perform elastic FE analysis with the updated equivalent shear modulus G k
k< maximum iteration number
̇ |dV = 2su |εmax
Klar and Osman [7] demonstrated that the upper-bound theorem of Eq. (1), expressed in the form of rate, could be rewritten as the energy conservation of the internal and external increment of work
Update the equivalent shear modulus at every integration point according to the stain field extracted from the previous FE analysis, k=k+1
NO
(1)
where Ti and fi are the vectors of surface traction and body force acting on the boundary S and V respectively; vi is the kinematically admissible ̇ | is the velocity field; su is the undrained shear strength of the soil; |εmax absolute largest principal component of the plastic stain rate εij̇ compatible with vi ; Δvt is the magnitude of the relative velocity change across the discontinuity surfaces A ; and Ti∗ is the upper-bound collapse load. For a continuous velocity field, Eq. (1) becomes
Fig. 1. Problem notation and example finite element layout of 3D undrained tunnel face stability.
Replace the old NPD with the current pressure difference s ; t Impose a new surcharge s and support pressure t , the difference of which equals to half of CPD+NPD
̇ |dV + ∫ Su Δvt dA ∫V 2Su |εmax A ∗ = ∫ Ti vi dS + ∫ fi vi dV (i = 1,2,3) S V
fi vi dV ⩽
NO
CPD=convergent pressure difference Yes Completed
2
Computers and Geotechnics xxx (xxxx) xxx–xxx
M. Huang et al.
14
ratio dε3 ranges from −1 to −0.5, η subsequently varies from 1 to dε1 1.077. The same results also can be obtained from the second case.
12
Thirdly, by defining ηMin Gs, |dε −u dε | 1 3 ∼ shear modulus G , Eq. (5) becomes:
(
( s- t ) /su
10 8
∫V
T-EMSD result C/D=0.5 C/D=1.0 C/D=2.0 C/D=3.0 C/D=4.0 C/D=5.0
6 4 2
40
80
120
160
200
uS2 center Ir /C Fig. 3. Load-displacement response of undrained tunnel face.
upper-bound solution. However, this relies on optimal computation. To break the dependency on optimal computation, Huang et al. [4] and Yu et al. [15] suggested transforming Eq. (3) into the form of the elastic virtual work equation, so that an upper-bound solution can be obtained by iterative computation. Firstly, at the initial incremental step of the EMSD method, εs equals to |dε1−dε3 |. Thus, Eq. (3) can be rewritten as
∫V
su ⎞ |dε1−dε3 |dεmax dV = 2Min ⎛Gs, − | dε dε3 | ⎠ 1 ⎝ ⎜
⎟
∫S Ti∗ dui dS + ∫V
(4)
su ⎞ (dε 2 + dε 2 + dε 2) dV = 2ηMin ⎛Gs, 1 2 3 |dε1−dε3 | ⎠ ⎝ ⎜
⎟
+
∫V
∫S Ti∗ dui dS
fi dui dV (i = 1,2,3)
(6)
At the end of the last section, a simple introduction has been given to illustrate the application of the T-EMSD in the problem with given displacement boundary conditions. More details can be found in Huang et al. [4] and Yu et al. [15]. However, tunnel face stability is a problem with mixed boundary conditions. The schematic diagram is shown in Fig. 1, where a rigid hollow cylinder tunnel of diameter D and cover depth C (measured from the crown to the ground surface) is modeled by fully fixed displacement boundaries. Following previous research, the tunnel is assumed to be ‘‘wish in place”, thereby neglecting the effect of tunneling process. Stress boundaries, uniform surcharge σs and support pressure σt , are applied over the ground surface S1 and on the tunnel face S2 respectively. Note that the compressive direction of σs and σt are defined as positive. Therefore, for the tunnel in undrained clay, Eq. (6) is simplified as:
Secondly, under the requirements of incompressibility (dε1 + dε2 + dε3 = 0 ), only two possible arrangements for the principal strain exist, namely dε1 ⩾ 0 ⩾ dε2 ⩾ dε3 (|dε1 | ⩾ |dε3 | ⩾ |dε2 |) and dε1 ⩾ dε2 ⩾ 0 ⩾ dε3 (|dε3 | ⩾ |dε1 | ⩾ |dε2 |). Taking the first case as an example, one can get |dε1−dε3 |dεmax = dε12−dε1 dε3 > 0 . Thus, Eq. (4) becomes
∫V
fi dui dV (i = 1,2,3)
2.2. Numerical implementation
fi dui dV (i
= 1,2,3)
∫S Ti∗ dui dS + ∫V
) as a nonlinear equivalent
which possesses the form of a typical elastic virtual work equation. In this study, η is taken as the upper value (1.077) for all strain states. This treatment ensures that the internal work is still an upper bound of the ∼ actual one, and that the equivalent shear modulus G is monotonic. As such, dui on the boundary S , the corresponding external force Ti∗ can be equivalently solved by nonlinear equilibrium iterations. As the derivation is based on the first incremental step of the EMSD approach, the load-displacement curve in this method is not computed by an incremental numerical procedure, but by a total displacement loading. More importantly, Huang et al. [4] and Yu et al. [15] have demonstrated that the external load Ti∗ would approach an upper-bound limit load when a sufficiently large dui is directly imposed on the boundary S . The velocity field for the upper-bound solution is proportional to the final displacement field.
0 0
∼ 2Gdεi dεi dV =
s
∫V (5)
∼ 2Gεi εi dV −
∫V
γdu v dV = σs
∫S
1
un dS−σt∗
∫S
2
un dS
(7)
where un is the normal displacement along boundary surfaces and u v is the vertical displacement; γ is the soil unit weight; S1 and S2 are the ground surface and tunnel face respectively. Further considering ∫S un dS = ∫S un dS = VL , due to incompressibility, where VL is the
dε 2 − dε dε
1 3 where η ⎛= dε 21+ dε 2 + dε 2 ⎞ is the function of soil strain state. According to ⎝ 1 2 3⎠ 1 − dε3 / dε1 . As the the incompressibility condition, η equals to 2
1
2[1 + dε3 / dε1 + (dε3 / dε1) ]
2
Fig. 4. Evolution of velocity field under various convergent pressure differences.
3
Computers and Geotechnics xxx (xxxx) xxx–xxx
M. Huang et al.
Upper-bound solutions with multi-block mechanism Davis et al. [3] Mollon et al. [10] Mollon et al. [11]
∫V
10
0
Elasto-plastic FE analysis Zhou [16] Ukritchon et al. [14] 0
1
2
3
4
5
Cover depth ratio, C/D
(a) Weightless clay with uniform strength
Tunnel face stability factor N s
60 Upper-bound solutions with continuous mechanism Zhang et al. [17] This study
50
Elasto-plastic FE analysis Zhou [16] Ukritchon et al. [14]
40 30 20 10
ϒ D/su0=0.0, C/D=3.0
0 0.00
0.25
0.50
0.75
r
Tunnel face stability factor Ns
40
in weightless clay
Upper-bound solutions with continuous mechanism Zhang et al. [17] This study Elasto-plastic FE analysis Zhou [16] Ukritchon et al. [14]
30 20 10 0 0.00
ϒ D/s u0=3.0, C/D=3.0 0.25
0.50
0.75
∫S
1
un dS
(8)
Gs su
does not affect the limit loads. Usually, a consistent rigidity index of 500 is taken in the whole model. If the effects of soil weight is considered, a geostatic step will be used in ABAQUS to establish an isotropic initial geostress field corresponding to a soil weight of γ . On the symmetrical plane, normal displacements are fixed and displacements in the other remaining directions are free. To satisfy the continuity requirement, displacements of the bottom and the other remaining vertical planes are fully fixed in all directions. Otherwise, the frictional energy dissipation due to the soil movements relative to external rigid soil on the boundaries would be ignored, thereby causing the present T-EMSD results to be no longer rigorous upper bounds. Yu et al. [15] demonstrated that upper-bound solutions in the TEMSD would decrease and tend to a stable value with the model domain increasing. Therefore, examinations on the boundary effect of the FE models used in this study have been carried out, and reveal that the current model domains are large enough to yield stable upper-bound solutions.
60 50
γdu v dV = (σs−σt )·
(I = ) only affects the development of load-displacement curves and
1.00
λD/s u0 (b) Effect of
∫V
is defined as the tunnel Following Davis et al. [3], the ratio face stability factor Ns . In terms of the T-EMSD approach, the pressure difference σs−σt would approach an upper-bound solution as VL increases. However, since the displacement patterns on the boundaries S1 and S2 are unknown, the problem has to be further converted to solve the maximum pressure difference to satisfy Eq. (8). The flow diagram in Fig. 2 illustrates the detailed calculation process, which consists of two parts: 1) the highlighted area exhibits the process of the equilibrium iteration under the given pressure difference; 2) the external loop is based on the method of bisection to determine the maximum pressure difference. The elastic FE analysis is implemented through the FEM program ABAQUS. Python, the built-in scripting language of ABAQUS, is used to establish external loops, internal equilibrium iteration loops and to extract the strain fields. During the equilibrium iteration, each loop ∼ involves updating the equivalent shear modulus G at every integration point using the strain field extracted from the previous loop and the elastic calculation with the newly updated shear modulus. The user subroutine USDFLD is used to update the material parameters in each loop. The equilibrium iteration will be stopped when the L2 norm of the displacement vector u between two neighboring loops is below the allowance value 1 of 10−4 or when the iteration number is larger than the prescribed number of 200. The external loop will be stopped when the relative difference between the convergent and non-convergent pressure difference is less than the allowance value 2 of 10−7. As an example, Fig. 1 illustrates the FE mesh for the 3D tunnel face with C / D = 1. Due to the existence of a plane of symmetry, only half of the domain is considered. A fine grid region close to the tunnel face and a far region with coarse grid are adopted to obtain accurate results with reasonable computation cost. To simulate incompressibility, a Poisson’s ratio of 0.499 is selected in the elastic FE analyses. In this study, the clay is assumed to be undrained with undrained shear strength su . Yu et al. [15] has proved that, in the T-EMSD, the soil rigidity index
15
5
∼ 2Gεi εi dV −
σs − σt su
Upper-bound solutions with continuous mechanism Klar et al. [6] Klar & Klein. [5] This study
20
Tunnel face stability factor N s
volume loss, Eq. (7) becomes
3. Results and discussions
1.00
λD/s u0
The inner loop (highlighted area in Fig. 2) carries out the equilibrium iteration under a given pressure difference σs−σt . Fig. 3 shows the load-displacement curves of an undrained tunnel face buried in weightless clay of uniform strength, where the y-axis and x-axis are the convergent pressure differences (CPD) σs−σt after the equilibrium iteration and their corresponding displacements uS 2,center at the center of the tunnel face. uS 2,center shows a rapid growth when the CPD approaches the upper-bound limit capacities. Fig. 4 further illustrates the changes
Fig. 5. Comparison of stability factor for undrained tunnel face under various geotechnical conditions.
4
Computers and Geotechnics xxx (xxxx) xxx–xxx
M. Huang et al.
Fig. 6. Velocity fields varying with cover depth in weightless clay with uniform strength.
of the T-EMSD velocity fields under different values of CPD for C / D = 1.0 . (Note that the velocity fields in Fig. 4 are proportional to the convergent displacement fields after the equilibrium iteration.) When the CPD is far less than the minimum non-convergent pressure difference (NPD), the velocity field spreads throughout the entire calculated region. With the increase of the CPD, the corresponding velocity field becomes more concentrated around the tunnel face. Eventually, an upper-bound velocity field is produced as the value of CPD approaches the NPD, which gives the upper-bound stability factor in Fig. 3. The stability factors for the tunnels in weightless clay with uniform strength are plotted in Fig. 5(a), together with previously published results. Compared with the other upper-bound solutions, the current results are closer to the trends of elasto-plastic FE analyses [14,16]. Especially for Ukritchon et al. [14], their difference are within 2%. In contrast, the results of the multi-block upper-bound methods are fairly close to the numerical solutions for C / D < 1.0 but gradually diverge from the FEM solutions for C / D > 1.0 . The T-EMSD method is an upper-bound analysis based on a continuous velocity field. Compared with the other continuous field based upper-bound solutions [5,6], this study presents the lowest prediction. The results of Mollon et al. [9] are not included, because the orthogonal strain tensor inconsistent with the adopted curvilinear velocity vectors
were used in energy calculation as pointed out by Klar and Klein [5]. The advantage of this method increases remarkably with an increasing C / D ratio, which is 41.07% lower than the result of [6] at C / D = 5.0 and 31.19% lower than the result of [5] at C / D = 4.0 . However, for C / D < 1.0 , this method fails in producing Ns smaller than those of the multi-block upper-bound method. This might be attributed to frictional energy dissipation being dominant for tunnel with small cover depth, such that fictional mechanisms become more suitable than purely continuous mechanisms. Fig. 6 illustrates the velocity fields obtained from the final displacement field, with C / D ranging from 1.0 to 5.0. For C / D ⩽ 2.0 , a pipe shaped mechanism similar to that proposed by [3] is mobilized; while for C / D ⩾ 3.0 , a hemisphere shaped mechanism with velocity distribution mainly concentrated around the tunnel face, which is similar to the FEM results [16], is produced. The velocity field shows a gradual and continuous evolution as the cover depth is increased. Consequently, with an increase in cover depth, the velocity field over the ground surface shows a reduction in magnitude as the field is spread within a larger ground surface. This means that there is a greater energy dissipation around the tunnel face as compared to the ground surface with an increasing cover depth. This phenomenon could not be captured using the multi-block mechanism. As such, for tunnels with large 5
Computers and Geotechnics xxx (xxxx) xxx–xxx
M. Huang et al.
Fig. 7. Effect of strength non-uniformity and soil weight on velocity fields.
(a) Weightless clay
(b) Non-weightless clay
factors, where an overburden stress factor γD is taken as 3. The insu0 troduction of soil weight would require larger support pressures σt to maintain the tunnel face stability, such that the current factors are much smaller than those from weightless clay (see Fig. 5(b)). A good agreement with the FE analyses [14,16] shows that the T-EMSD can provide sufficiently accurate upper-bound predictions for the tunnel face stability factors under complex geotechnical conditions. Fig. 7(b) shows that in the case of a non-weightless clay, the focal point of the velocity field shows a downwards shift when compared to that of a weightless clay. However, with the decrease in uniformity, the focal point shows a similar trend where the focal point is gradually shifted to the crown of the tunnel face. It should be noted that, in theory, the velocity fields shown in Figs. 6 and 7 should correspond to lower upper-bound stability factors, as η is taken as the upper value (1.077) for all soil strain states in this study. Therefore, better upper-bound results can be expected when more reasonable treatment to η is presented.
cover depth, the use of the multi-block upper-bound mechanism yields less accurate results as compared to the T-EMSD method. On the other hand, multi-block mechanism capable of producing large ratio of the velocity magnitude at ground surface over that at ground surface is more suitable for shallowly covered tunnels. This explains why the present method fails to yield more accurate results than the multi-block solutions for shallowly covered tunnel faces even though it is more successful in producing various mechanisms for different C / D values. Fig. 5(b) shows the application of the T-EMSD method in tunnel face stability in weightless clay with an undrained shear strength of su0 at ground surface and a strength gradient of λ , in which the cover depth ratio C / D of 3 is considered and the stability factor is defined as σs − σt . It su0 shows that the T-EMSD results are fairly close to those of the FE analyses [14,16], and much lower than the continuous field based upperbound solution of Zhang [17]. The corresponding velocity fields are illustrated in Fig. 7(a). Under a fully uniform strength condition
(
λD su0
)
= 0 , the velocity field tends to show a relatively uniform dis-
tribution about the tunnel face. However, as the uniformity decreases, there is a development of a concentrated velocity field about the face of the tunnel, with the focal point shifting gradually to the crown of the tunnel face. Fig. 5(c) further shows the effect of soil weight on the stability
4. Conclusions A numerical upper-bound method (T-EMSD), which could be implemented in displacement-based FEM packages, is extended to analyze 6
Computers and Geotechnics xxx (xxxx) xxx–xxx
M. Huang et al.
path approach. Chin J Geotech Eng 2015;37(3):400–9. (in Chinese). [5] Klar A, Klein B. Energy-based volume loss prediction for tunnel face advancement in clays. Géotechnique 2014;64(20):776–86. [6] 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. [7] Klar A, Osman AS. Load-displacement solutions for piles and shallow foundations based on deformation fields and energy conservation. Géotechnique 2008;58(7):581–90. [8] Leca E, Dormieux L. Upper and lower bound solutions for the face stability of shallow circular tunnels in frictional material. Géotechnique 1990;40(4):581–606. [9] Mollon G, Dias D, Soubra A. Continuous velocity fields for collapse and blowout of a pressurized tunnel face in purely cohesive soil. Int J Numer Anal Meth Geomech 2012;37(13):2061–83. [10] Mollon G, Dias D, Soubra A. Face stability analysis of circular tunnels driven by a pressurized shield. J Geotech Geoenviron Eng 2010;136(1):215–29. [11] Mollon G, Dias D, Soubra A. Rotational failure mechanisms for the face stability analysis of tunnels driven by a pressurized shield. Int J Numer Anal Meth Geomech 2011;35(12):1363–88. [12] Shield RT, Drucker DC. The application of limit analysis to punch indentation problems. J Appl Mech, ASME 1953;20:453–60. [13] Sloan SW. Geotechnical stability analysis. Géotechnique 2013;63(7):531–72. [14] 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. [15] Yu J, Huang M, Li S, et al. Load-displacement and upper-bound solutions of a loaded laterally pile in clay based on a total-displacement-loading EMSD method. Comput Geotech 2017;83:64–76. [16] Zhou W. Stability analysis of tunnel face in nonhomogeneous clays Master thesis. Shanghai: Tongji University; 2011. (In Chinese). [17] Zhang F, Gao YF, Wu YX, et al. Upper-bound solutions for face stability of circular tunnels in undrained clays. Géotechnique 2017:1–10.
3D undrained tunnel face stability with mixed boundary conditions. Comparisons with previous published results reveal that this method could provide sufficiently accurate upper-bound prediction for the tunnel face stability factors under complex geotechnical conditions. In addition, the velocity fields varying with cover depths are investigated, which have a similar evolving pattern with those of the FEM. A discussion is made to compare respective advantages of multiblock mechanisms and continuous mechanisms in tunnel face stability problem. Acknowledgements This work was financially supported by the National Key Research and Development Program (Grant No. 2016YFC0800202) and the National Fund for Distinguished Young Scholars of China (Grant No. 50825803). These supports are gratefully acknowledged. References [1] Augarde CE, Lyamin AV, Sloan SW. Stability of an undrained plane strain heading revisited. Comput Geotech 2003;30(5):419–30. [2] Broms BB, Bennermark H. Stability of clay in vertical openings. J Geotech Eng Div 1967;193:71–94. [3] Davis EH, Gunn MJ, Mair RJ, et al. The stability of shallow tunnels and underground openings in cohesive material. Géotechnique 1980;30(4):397–416. [4] Huang M, Yu J, Zhang C. P–y curves of laterally loaded piles in clay based on strain
7
Contents lists available at ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Technical Communication
Continuous field based upper bound analysis for three-dimensional tunnel face stability in undrained clay ⁎
Maosong Huanga, Sen Lia, Jian Yua, , Jorgin Qi Wen Tanb a b
Department of Geotechnical Engineering, Tongji University, Shanghai 20092, China Geotechnical and Tunneling, Land Transport Authority, 219428, Singapore
A R T I C L E I N F O
A B S T R A C T
Keywords: Undrained clay Tunnel face stability Upper-bound limit analysis Finite element Continuous field
A new numerical upper-bound method applicable to three-dimensional stability problems in undrained clay is proposed and applied in the analysis of tunnel face stability. Comparisons with the elasto-plastic finite element method (FEM) reveal that this method could provide sufficiently accurate upper-bound prediction for the tunnel face stability factors. For tunnels with large cover depth, the obtained results greatly improve the existing upperbound solutions. Their corresponding velocity fields are investigated, which show a similar evolving pattern with those of the FEM. A discussion is made to compare the various advantages of multi-block mechanisms and continuous mechanisms in tunnel face stability.
1. Introduction Underground tunnel excavation has become commonplace in infrastructure projects worldwide. An important design parameter in tunnel face stability is the support pressure that is applied at the tunnel face. Due to the practical significance of three-dimensional (3D) tunnel face stability, substantial research investigations have been conducted with various methods including the empirical method [2], lower bound limit analyses [3,8,13], upper-bound limit analyses [3,5,6,8–11,13] and elasto-plastic finite element (FE) analyses [14,16]. Broms and Bennermark [2] first proposed a hand calculation design expression for tunnel face stability in undrained clay by empirically relating the expected deformations with various factors such as support pressure, surcharge, soil unit weight, cover depth and tunnel diameter to determine whether the tunnel is sufficiently supported. The first theoretically rigorous analysis of 3D tunnel face stability was carried out by Davis et al. [3] using both upper and lower bound limit analyses. Leca and Dormieux [8] extended the work of [3] into fictional soil. Ukritchon et al. [14] and Zhou [16] further utilized the finite element method (FEM) to reevaluate the accuracy of the previously published results. As a widely used analysis method, the applications of the upperbound method in 3D tunnel face stability are further reviewed. The upper-bound method can be divided into three categories according to their underlying mechanisms. The first approach is to postulate mechanisms composed of either translational [3,8,10] or rotational rigid blocks [11] where the energy dissipation exists solely at the interfaces between adjacent blocks. The second approach is to develop continuous
⁎
failure mechanisms with no discontinuity, where the energy is solely dissipated by plastic deformation throughout the continuums [5,6,9]. The third approach is to construct failure mechanisms with both deformable zones and discontinuity surfaces [13]. The continuous velocity field of Klar et al. [6] was derived from both the elasticity solution corresponding to the original plastic problem and the concept of sinks and sources. The continuous mechanism of Mollon et al. [9] was obtained from the elasto-plastic finite difference method (FDM). The finite element limit analysis (FELA) was employed by Sloan [13] to investigate 3D tunnel face stability. However, the difference between the upper and lower bound for 3D tunnel problem is quite large when compared with the 2D FELA solutions [1]. Much effort has been put into constructing upper-bound solutions for 3D undrained tunnel stability. However, a significant gap still exists between the upper-bound results and the elasto-plastic FEM results for tunnels with large cover depth. Such a gap might be attributed to the difficulty in constructing complicated 3D upper-bound velocity fields. This paper attempts to propose a generally applicable 3D upper-bound method for undrained clay based on the total loading extended mobilisable strength design (T-EMSD) proposed by Huang et al. [4] and Yu et al. [15], which can acquire the upper-bound solution through elastic iterative analysis and thus could be conveniently replicated through an elastic displacement-based FEM package. The method will be applied to 3D undrained tunnel stability, and the obtained upper-bound solutions will be thoroughly compared with existing solutions.
Corresponding author. E-mail address: [email protected] (J. Yu).
http://dx.doi.org/10.1016/j.compgeo.2017.09.014 Received 27 June 2017; Received in revised form 30 August 2017; Accepted 19 September 2017 0266-352X/ © 2017 Elsevier Ltd. All rights reserved.
Please cite this article as: Huang, M., Computers and Geotechnics (2017), http://dx.doi.org/10.1016/j.compgeo.2017.09.014
Computers and Geotechnics xxx (xxxx) xxx–xxx
M. Huang et al.
2.1. Theory Shield and Drucker [12] demonstrated that the upper-bound theorem for Tresca material, conventionally adopted for undrained clay, can be expressed as
∫s Ti vi dS + ∫V
2. Method of analysis
∫V
The T-EMSD, developed from the EMSD, is a total loading method to obtain an upper-bound load-displacement solution of a plane-strain problem through iterative computation, rather than optimal computation. It has been verified in the study of the soil resistance to the lateral movement of a 2D circular pile. The T-EMSD is able to directly calculate an upper-bound lateral resistance under a given total pile displacement (without the need for an incremental loading process). Large enough displacements would enable the resistance to approach an upper-bound lateral limit load. This section will extend the T-EMSD theory to a generic 3D condition and illustrate its application in obtaining the upper-bound tunnel face stability factor.
∫V
Replace the old CPD with the current pressure difference s ; t Impose a new surchage s and support pressure t , the difference of which equals to half of CPD+NPD
Yes
uk
1 2
allowance value 1
Yes
Note: NPD=non-convergent pressure difference
NPD CPD NPD
allowance value 2
∫S
∫V
Fig. 2. Flow diagram for tunnel stability.
Perform elastic FE analysis with the equivalent shear modulus G d 1 d 3 , iteration number k=1
uk
fi vi dV (i = 1,2,3)
(2)
(3)
where dεmax = max(|dε1 |,|dε3 |,|dε3 |) is the increment of the largest absolute principal strain; cm (εs ) is the mobilised strength, which is a function n of engineering shear strain εs (= ∑1 |dε1−dε3 |, n is the number of the loading steps); for an elastic perfectly-plastic material cm (εs ) equals to Min (Gs εs,su ) , where Gs is the shear modulus of the soil. With εs accumulating during loading, cm (εs ) tends to su such that Ti∗ tends to an
Impose a surcharge s and support pressure t on the ground surface and tunnel face; record 0 as CPD and a sufficiently large value as NPD
NO
∫S Ti∗ vi dS + ∫V
Internal increment of work External increment of work 2cm (εs) dεmax dV = Ti∗ dui dS + fi dui dV (i = 1,2,3)
Perform elastic FE analysis with the updated equivalent shear modulus G k
k< maximum iteration number
̇ |dV = 2su |εmax
Klar and Osman [7] demonstrated that the upper-bound theorem of Eq. (1), expressed in the form of rate, could be rewritten as the energy conservation of the internal and external increment of work
Update the equivalent shear modulus at every integration point according to the stain field extracted from the previous FE analysis, k=k+1
NO
(1)
where Ti and fi are the vectors of surface traction and body force acting on the boundary S and V respectively; vi is the kinematically admissible ̇ | is the velocity field; su is the undrained shear strength of the soil; |εmax absolute largest principal component of the plastic stain rate εij̇ compatible with vi ; Δvt is the magnitude of the relative velocity change across the discontinuity surfaces A ; and Ti∗ is the upper-bound collapse load. For a continuous velocity field, Eq. (1) becomes
Fig. 1. Problem notation and example finite element layout of 3D undrained tunnel face stability.
Replace the old NPD with the current pressure difference s ; t Impose a new surcharge s and support pressure t , the difference of which equals to half of CPD+NPD
̇ |dV + ∫ Su Δvt dA ∫V 2Su |εmax A ∗ = ∫ Ti vi dS + ∫ fi vi dV (i = 1,2,3) S V
fi vi dV ⩽
NO
CPD=convergent pressure difference Yes Completed
2
Computers and Geotechnics xxx (xxxx) xxx–xxx
M. Huang et al.
14
ratio dε3 ranges from −1 to −0.5, η subsequently varies from 1 to dε1 1.077. The same results also can be obtained from the second case.
12
Thirdly, by defining ηMin Gs, |dε −u dε | 1 3 ∼ shear modulus G , Eq. (5) becomes:
(
( s- t ) /su
10 8
∫V
T-EMSD result C/D=0.5 C/D=1.0 C/D=2.0 C/D=3.0 C/D=4.0 C/D=5.0
6 4 2
40
80
120
160
200
uS2 center Ir /C Fig. 3. Load-displacement response of undrained tunnel face.
upper-bound solution. However, this relies on optimal computation. To break the dependency on optimal computation, Huang et al. [4] and Yu et al. [15] suggested transforming Eq. (3) into the form of the elastic virtual work equation, so that an upper-bound solution can be obtained by iterative computation. Firstly, at the initial incremental step of the EMSD method, εs equals to |dε1−dε3 |. Thus, Eq. (3) can be rewritten as
∫V
su ⎞ |dε1−dε3 |dεmax dV = 2Min ⎛Gs, − | dε dε3 | ⎠ 1 ⎝ ⎜
⎟
∫S Ti∗ dui dS + ∫V
(4)
su ⎞ (dε 2 + dε 2 + dε 2) dV = 2ηMin ⎛Gs, 1 2 3 |dε1−dε3 | ⎠ ⎝ ⎜
⎟
+
∫V
∫S Ti∗ dui dS
fi dui dV (i = 1,2,3)
(6)
At the end of the last section, a simple introduction has been given to illustrate the application of the T-EMSD in the problem with given displacement boundary conditions. More details can be found in Huang et al. [4] and Yu et al. [15]. However, tunnel face stability is a problem with mixed boundary conditions. The schematic diagram is shown in Fig. 1, where a rigid hollow cylinder tunnel of diameter D and cover depth C (measured from the crown to the ground surface) is modeled by fully fixed displacement boundaries. Following previous research, the tunnel is assumed to be ‘‘wish in place”, thereby neglecting the effect of tunneling process. Stress boundaries, uniform surcharge σs and support pressure σt , are applied over the ground surface S1 and on the tunnel face S2 respectively. Note that the compressive direction of σs and σt are defined as positive. Therefore, for the tunnel in undrained clay, Eq. (6) is simplified as:
Secondly, under the requirements of incompressibility (dε1 + dε2 + dε3 = 0 ), only two possible arrangements for the principal strain exist, namely dε1 ⩾ 0 ⩾ dε2 ⩾ dε3 (|dε1 | ⩾ |dε3 | ⩾ |dε2 |) and dε1 ⩾ dε2 ⩾ 0 ⩾ dε3 (|dε3 | ⩾ |dε1 | ⩾ |dε2 |). Taking the first case as an example, one can get |dε1−dε3 |dεmax = dε12−dε1 dε3 > 0 . Thus, Eq. (4) becomes
∫V
fi dui dV (i = 1,2,3)
2.2. Numerical implementation
fi dui dV (i
= 1,2,3)
∫S Ti∗ dui dS + ∫V
) as a nonlinear equivalent
which possesses the form of a typical elastic virtual work equation. In this study, η is taken as the upper value (1.077) for all strain states. This treatment ensures that the internal work is still an upper bound of the ∼ actual one, and that the equivalent shear modulus G is monotonic. As such, dui on the boundary S , the corresponding external force Ti∗ can be equivalently solved by nonlinear equilibrium iterations. As the derivation is based on the first incremental step of the EMSD approach, the load-displacement curve in this method is not computed by an incremental numerical procedure, but by a total displacement loading. More importantly, Huang et al. [4] and Yu et al. [15] have demonstrated that the external load Ti∗ would approach an upper-bound limit load when a sufficiently large dui is directly imposed on the boundary S . The velocity field for the upper-bound solution is proportional to the final displacement field.
0 0
∼ 2Gdεi dεi dV =
s
∫V (5)
∼ 2Gεi εi dV −
∫V
γdu v dV = σs
∫S
1
un dS−σt∗
∫S
2
un dS
(7)
where un is the normal displacement along boundary surfaces and u v is the vertical displacement; γ is the soil unit weight; S1 and S2 are the ground surface and tunnel face respectively. Further considering ∫S un dS = ∫S un dS = VL , due to incompressibility, where VL is the
dε 2 − dε dε
1 3 where η ⎛= dε 21+ dε 2 + dε 2 ⎞ is the function of soil strain state. According to ⎝ 1 2 3⎠ 1 − dε3 / dε1 . As the the incompressibility condition, η equals to 2
1
2[1 + dε3 / dε1 + (dε3 / dε1) ]
2
Fig. 4. Evolution of velocity field under various convergent pressure differences.
3
Computers and Geotechnics xxx (xxxx) xxx–xxx
M. Huang et al.
Upper-bound solutions with multi-block mechanism Davis et al. [3] Mollon et al. [10] Mollon et al. [11]
∫V
10
0
Elasto-plastic FE analysis Zhou [16] Ukritchon et al. [14] 0
1
2
3
4
5
Cover depth ratio, C/D
(a) Weightless clay with uniform strength
Tunnel face stability factor N s
60 Upper-bound solutions with continuous mechanism Zhang et al. [17] This study
50
Elasto-plastic FE analysis Zhou [16] Ukritchon et al. [14]
40 30 20 10
ϒ D/su0=0.0, C/D=3.0
0 0.00
0.25
0.50
0.75
r
Tunnel face stability factor Ns
40
in weightless clay
Upper-bound solutions with continuous mechanism Zhang et al. [17] This study Elasto-plastic FE analysis Zhou [16] Ukritchon et al. [14]
30 20 10 0 0.00
ϒ D/s u0=3.0, C/D=3.0 0.25
0.50
0.75
∫S
1
un dS
(8)
Gs su
does not affect the limit loads. Usually, a consistent rigidity index of 500 is taken in the whole model. If the effects of soil weight is considered, a geostatic step will be used in ABAQUS to establish an isotropic initial geostress field corresponding to a soil weight of γ . On the symmetrical plane, normal displacements are fixed and displacements in the other remaining directions are free. To satisfy the continuity requirement, displacements of the bottom and the other remaining vertical planes are fully fixed in all directions. Otherwise, the frictional energy dissipation due to the soil movements relative to external rigid soil on the boundaries would be ignored, thereby causing the present T-EMSD results to be no longer rigorous upper bounds. Yu et al. [15] demonstrated that upper-bound solutions in the TEMSD would decrease and tend to a stable value with the model domain increasing. Therefore, examinations on the boundary effect of the FE models used in this study have been carried out, and reveal that the current model domains are large enough to yield stable upper-bound solutions.
60 50
γdu v dV = (σs−σt )·
(I = ) only affects the development of load-displacement curves and
1.00
λD/s u0 (b) Effect of
∫V
is defined as the tunnel Following Davis et al. [3], the ratio face stability factor Ns . In terms of the T-EMSD approach, the pressure difference σs−σt would approach an upper-bound solution as VL increases. However, since the displacement patterns on the boundaries S1 and S2 are unknown, the problem has to be further converted to solve the maximum pressure difference to satisfy Eq. (8). The flow diagram in Fig. 2 illustrates the detailed calculation process, which consists of two parts: 1) the highlighted area exhibits the process of the equilibrium iteration under the given pressure difference; 2) the external loop is based on the method of bisection to determine the maximum pressure difference. The elastic FE analysis is implemented through the FEM program ABAQUS. Python, the built-in scripting language of ABAQUS, is used to establish external loops, internal equilibrium iteration loops and to extract the strain fields. During the equilibrium iteration, each loop ∼ involves updating the equivalent shear modulus G at every integration point using the strain field extracted from the previous loop and the elastic calculation with the newly updated shear modulus. The user subroutine USDFLD is used to update the material parameters in each loop. The equilibrium iteration will be stopped when the L2 norm of the displacement vector u between two neighboring loops is below the allowance value 1 of 10−4 or when the iteration number is larger than the prescribed number of 200. The external loop will be stopped when the relative difference between the convergent and non-convergent pressure difference is less than the allowance value 2 of 10−7. As an example, Fig. 1 illustrates the FE mesh for the 3D tunnel face with C / D = 1. Due to the existence of a plane of symmetry, only half of the domain is considered. A fine grid region close to the tunnel face and a far region with coarse grid are adopted to obtain accurate results with reasonable computation cost. To simulate incompressibility, a Poisson’s ratio of 0.499 is selected in the elastic FE analyses. In this study, the clay is assumed to be undrained with undrained shear strength su . Yu et al. [15] has proved that, in the T-EMSD, the soil rigidity index
15
5
∼ 2Gεi εi dV −
σs − σt su
Upper-bound solutions with continuous mechanism Klar et al. [6] Klar & Klein. [5] This study
20
Tunnel face stability factor N s
volume loss, Eq. (7) becomes
3. Results and discussions
1.00
λD/s u0
The inner loop (highlighted area in Fig. 2) carries out the equilibrium iteration under a given pressure difference σs−σt . Fig. 3 shows the load-displacement curves of an undrained tunnel face buried in weightless clay of uniform strength, where the y-axis and x-axis are the convergent pressure differences (CPD) σs−σt after the equilibrium iteration and their corresponding displacements uS 2,center at the center of the tunnel face. uS 2,center shows a rapid growth when the CPD approaches the upper-bound limit capacities. Fig. 4 further illustrates the changes
Fig. 5. Comparison of stability factor for undrained tunnel face under various geotechnical conditions.
4
Computers and Geotechnics xxx (xxxx) xxx–xxx
M. Huang et al.
Fig. 6. Velocity fields varying with cover depth in weightless clay with uniform strength.
of the T-EMSD velocity fields under different values of CPD for C / D = 1.0 . (Note that the velocity fields in Fig. 4 are proportional to the convergent displacement fields after the equilibrium iteration.) When the CPD is far less than the minimum non-convergent pressure difference (NPD), the velocity field spreads throughout the entire calculated region. With the increase of the CPD, the corresponding velocity field becomes more concentrated around the tunnel face. Eventually, an upper-bound velocity field is produced as the value of CPD approaches the NPD, which gives the upper-bound stability factor in Fig. 3. The stability factors for the tunnels in weightless clay with uniform strength are plotted in Fig. 5(a), together with previously published results. Compared with the other upper-bound solutions, the current results are closer to the trends of elasto-plastic FE analyses [14,16]. Especially for Ukritchon et al. [14], their difference are within 2%. In contrast, the results of the multi-block upper-bound methods are fairly close to the numerical solutions for C / D < 1.0 but gradually diverge from the FEM solutions for C / D > 1.0 . The T-EMSD method is an upper-bound analysis based on a continuous velocity field. Compared with the other continuous field based upper-bound solutions [5,6], this study presents the lowest prediction. The results of Mollon et al. [9] are not included, because the orthogonal strain tensor inconsistent with the adopted curvilinear velocity vectors
were used in energy calculation as pointed out by Klar and Klein [5]. The advantage of this method increases remarkably with an increasing C / D ratio, which is 41.07% lower than the result of [6] at C / D = 5.0 and 31.19% lower than the result of [5] at C / D = 4.0 . However, for C / D < 1.0 , this method fails in producing Ns smaller than those of the multi-block upper-bound method. This might be attributed to frictional energy dissipation being dominant for tunnel with small cover depth, such that fictional mechanisms become more suitable than purely continuous mechanisms. Fig. 6 illustrates the velocity fields obtained from the final displacement field, with C / D ranging from 1.0 to 5.0. For C / D ⩽ 2.0 , a pipe shaped mechanism similar to that proposed by [3] is mobilized; while for C / D ⩾ 3.0 , a hemisphere shaped mechanism with velocity distribution mainly concentrated around the tunnel face, which is similar to the FEM results [16], is produced. The velocity field shows a gradual and continuous evolution as the cover depth is increased. Consequently, with an increase in cover depth, the velocity field over the ground surface shows a reduction in magnitude as the field is spread within a larger ground surface. This means that there is a greater energy dissipation around the tunnel face as compared to the ground surface with an increasing cover depth. This phenomenon could not be captured using the multi-block mechanism. As such, for tunnels with large 5
Computers and Geotechnics xxx (xxxx) xxx–xxx
M. Huang et al.
Fig. 7. Effect of strength non-uniformity and soil weight on velocity fields.
(a) Weightless clay
(b) Non-weightless clay
factors, where an overburden stress factor γD is taken as 3. The insu0 troduction of soil weight would require larger support pressures σt to maintain the tunnel face stability, such that the current factors are much smaller than those from weightless clay (see Fig. 5(b)). A good agreement with the FE analyses [14,16] shows that the T-EMSD can provide sufficiently accurate upper-bound predictions for the tunnel face stability factors under complex geotechnical conditions. Fig. 7(b) shows that in the case of a non-weightless clay, the focal point of the velocity field shows a downwards shift when compared to that of a weightless clay. However, with the decrease in uniformity, the focal point shows a similar trend where the focal point is gradually shifted to the crown of the tunnel face. It should be noted that, in theory, the velocity fields shown in Figs. 6 and 7 should correspond to lower upper-bound stability factors, as η is taken as the upper value (1.077) for all soil strain states in this study. Therefore, better upper-bound results can be expected when more reasonable treatment to η is presented.
cover depth, the use of the multi-block upper-bound mechanism yields less accurate results as compared to the T-EMSD method. On the other hand, multi-block mechanism capable of producing large ratio of the velocity magnitude at ground surface over that at ground surface is more suitable for shallowly covered tunnels. This explains why the present method fails to yield more accurate results than the multi-block solutions for shallowly covered tunnel faces even though it is more successful in producing various mechanisms for different C / D values. Fig. 5(b) shows the application of the T-EMSD method in tunnel face stability in weightless clay with an undrained shear strength of su0 at ground surface and a strength gradient of λ , in which the cover depth ratio C / D of 3 is considered and the stability factor is defined as σs − σt . It su0 shows that the T-EMSD results are fairly close to those of the FE analyses [14,16], and much lower than the continuous field based upperbound solution of Zhang [17]. The corresponding velocity fields are illustrated in Fig. 7(a). Under a fully uniform strength condition
(
λD su0
)
= 0 , the velocity field tends to show a relatively uniform dis-
tribution about the tunnel face. However, as the uniformity decreases, there is a development of a concentrated velocity field about the face of the tunnel, with the focal point shifting gradually to the crown of the tunnel face. Fig. 5(c) further shows the effect of soil weight on the stability
4. Conclusions A numerical upper-bound method (T-EMSD), which could be implemented in displacement-based FEM packages, is extended to analyze 6
Computers and Geotechnics xxx (xxxx) xxx–xxx
M. Huang et al.
path approach. Chin J Geotech Eng 2015;37(3):400–9. (in Chinese). [5] Klar A, Klein B. Energy-based volume loss prediction for tunnel face advancement in clays. Géotechnique 2014;64(20):776–86. [6] 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. [7] Klar A, Osman AS. Load-displacement solutions for piles and shallow foundations based on deformation fields and energy conservation. Géotechnique 2008;58(7):581–90. [8] Leca E, Dormieux L. Upper and lower bound solutions for the face stability of shallow circular tunnels in frictional material. Géotechnique 1990;40(4):581–606. [9] Mollon G, Dias D, Soubra A. Continuous velocity fields for collapse and blowout of a pressurized tunnel face in purely cohesive soil. Int J Numer Anal Meth Geomech 2012;37(13):2061–83. [10] Mollon G, Dias D, Soubra A. Face stability analysis of circular tunnels driven by a pressurized shield. J Geotech Geoenviron Eng 2010;136(1):215–29. [11] Mollon G, Dias D, Soubra A. Rotational failure mechanisms for the face stability analysis of tunnels driven by a pressurized shield. Int J Numer Anal Meth Geomech 2011;35(12):1363–88. [12] Shield RT, Drucker DC. The application of limit analysis to punch indentation problems. J Appl Mech, ASME 1953;20:453–60. [13] Sloan SW. Geotechnical stability analysis. Géotechnique 2013;63(7):531–72. [14] 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. [15] Yu J, Huang M, Li S, et al. Load-displacement and upper-bound solutions of a loaded laterally pile in clay based on a total-displacement-loading EMSD method. Comput Geotech 2017;83:64–76. [16] Zhou W. Stability analysis of tunnel face in nonhomogeneous clays Master thesis. Shanghai: Tongji University; 2011. (In Chinese). [17] Zhang F, Gao YF, Wu YX, et al. Upper-bound solutions for face stability of circular tunnels in undrained clays. Géotechnique 2017:1–10.
3D undrained tunnel face stability with mixed boundary conditions. Comparisons with previous published results reveal that this method could provide sufficiently accurate upper-bound prediction for the tunnel face stability factors under complex geotechnical conditions. In addition, the velocity fields varying with cover depths are investigated, which have a similar evolving pattern with those of the FEM. A discussion is made to compare respective advantages of multiblock mechanisms and continuous mechanisms in tunnel face stability problem. Acknowledgements This work was financially supported by the National Key Research and Development Program (Grant No. 2016YFC0800202) and the National Fund for Distinguished Young Scholars of China (Grant No. 50825803). These supports are gratefully acknowledged. References [1] Augarde CE, Lyamin AV, Sloan SW. Stability of an undrained plane strain heading revisited. Comput Geotech 2003;30(5):419–30. [2] Broms BB, Bennermark H. Stability of clay in vertical openings. J Geotech Eng Div 1967;193:71–94. [3] Davis EH, Gunn MJ, Mair RJ, et al. The stability of shallow tunnels and underground openings in cohesive material. Géotechnique 1980;30(4):397–416. [4] Huang M, Yu J, Zhang C. P–y curves of laterally loaded piles in clay based on strain
7