Face stability of shield tunnels considering a kinematically admissible velocity field of soil arching

Wei Li, Chengping Zhang,*, Dingli Zhang, Zijian Ye, Zhibiao Tan

a Key Laboratory for Urban Underground Engineering of the Ministry of Education, Beijing Jiaotong University, Beijing,100044, China

b School of Civil Engineering, Beijing Jiaotong University, Beijing,100044, China

Keywords:Tunnel face stability Velocity field Failure pattern Improved failure mechanism Critical face pressure

ABSTRACT Existing mechanism of simulating soil movement at tunnel face is generally based on the translational or rotational velocity field, which is, to some extent, different from the real soil movement in the arching zone.Numerical simulations are carried out first to investigate the characteristics of the velocity distribution at tunnel face and above tunnel vault.Then a new kinematically admissible velocity field is proposed to improve the description of the soil movement according to the results of the numerical simulation.Based on the proposed velocity field,an improved failure mechanism is constructed adopting the spatial discretization technique,which takes into account soil arching effect and plastic deformation within soil mass.Finally, the critical face pressure and the proposed mechanism are compared with the results of the numerical simulation, existing analytical studies and experimental tests to verify the accuracy and improvement of the presented method.The proposed mechanism can serve as an alternative approach for the face stability analysis.

1.Introduction

Tunnel face stability problems have been receiving increasing attention.For engineering safety, it is crucial to control the face support pressure during the tunneling process for preventing tunnel face collapse and the subsidence of ground surface.To avoid face collapse and maintain face stability, the required minimal support pressure should be determined.It is thus desirable to establish a practical face failure mechanism and obtain a precise solution to the critical face pressure.

The stability of tunnel face has been investigated by many researchers.The existing approaches include the analytical method,experimental tests and numerical simulation.Both the numerical simulations (Augarde et al., 2003; De Buhan et al.,1999; Vermeer et al., 2002; Chen et al., 2011; Ukritchon and Keawsawasvong,2017; Ukritchon et al., 2017a, b; Keawsawasvong and Ukritchon,2019; Du et al., 2020; Shiau and Al-Asadi, 2020) and experimental tests (Broms and Bennermark,1967; Mair,1969; Schofield,1980; Atkinson and Potts, 1977; Chambon and Corte, 1994;Takano et al., 2006; Kirsch, 2010; Idinger et al., 2011; Chen et al.,2013) have been extensively performed to visualize the failure pattern of a tunnel face.The advantage of the numerical simulation lies in its good repeatability (Huang et al., 2018a; Ukritchon and Keawsawasvong, 2019a, b, c), while the experimental tests are good at capturing the characteristic of the face failure.As shown in Fig.1, when the tunnel face collapses, a soil arching and a shear band are formed above the tunnel face and ahead of the tunnel face,respectively,which could be adequately simulated with logarithmic spirals.

The limit analysis and limit equilibrium are two effective analytical methods to examine the face stability.The limit equilibrium method usually assumes the calculation model to deduce critical face pressures and consider different influence factors(Anagnostou and Kovári,1996;Anagnostou,2012;Anagnostou and Perazzelli, 2013, 2015; Perazzelli et al., 2014).The limit analysis consists of the upper and lower bound theorems, which result in the upper and lower bound solutions of critical loads by constructing kinematic and static mechanisms, respectively.A simple and intuitive kinematic mechanism was the rigid block failure mechanism based on the translational velocity field, which was first introduced by Davis et al.(1980) and Leca and Dormieux(1990), as shown in Fig.2a.Then, the failure mechanisms of thetranslational multi-block (Mollon et al., 2009) and the rotational rigid block(Subrin and Wong,2002)were proposed to improve the velocity field and allow freer development of the shear failure plane,as shown in Fig.2b and c.However,the failure mechanisms mentioned above suffered from the fact that only a portion of the tunnel face collapses but the remaining region was undisturbed.To overcome this shortcoming, an innovative spatial discretization technique was proposed by Mollon et al.(2010, 2011) to enhance the knowledge about the failure mechanism.Based on these failure mechanisms, researches on face stability influence factors and improvement methods were widely conducted(Zhang et al.,2015;Han et al., 2016a, b; Pan and Dias, 2016a, b, 2017; Zou and Qian,2018; Ding et al., 2018, 2019; Li et al., 2018, 2019a, b; Zou et al.,2019a, b; Li and Zhang, 2020).Moreover, a continuous velocity field was proposed by Mollon et al.(2013)for purely cohesive soils,which significantly improved the existing upper solutions.The proposed velocity field closely simulates the movement of the purely cohesive soils and provided better results for the related geotechnical problems (Osman et al., 2006; Klar et al., 2007; Klar and Klein, 2014; Huang et al., 2018b; Zhang et al., 2018a, b; Li et al., 2019c; Zhang et al., 2020).

Fig.1.Three-dimensional (3D) failure zone (Takano et al., 2006).

It can be found from the extensive literature review that the accuracy of the upper bound solution is closely related to the adopted kinematically admissible mechanism.The accuracy of the obtained critical face pressure and the correspondence of the failure mechanism to the real situation can be improved with increasing sensibility and credibility of the velocity field.However,most reported studies neglected or simplified the soil movement in the arching zone.Fig.2 shows that the simplified translational or rotational rigid movement used in the analytical method leads to a cone or‘horn’shape of the failure mechanism.Therefore,this paper intends to mitigate this problem by proposing a new kinematically admissible velocity field for the soil archinging zone.

Inspired by the advantages of a more promising and reasonable kinematically admissible velocity field,this paper aims at exploring a more realistic velocity field for the tunnel face in frictional soils.Numerical simulations will be carried out first to investigate the characteristics of the velocity distribution at the tunnel face and above the tunnel vault.Then, a kinematically admissible velocity field will be defined and deduced in detail to better correspond to the real soil movement.Finally,comparisons among the presented study, the numerical simulation and other existing studies will be performed to verify the accuracy and improvement of the proposed mechanism.

2.Numerical simulations

2.1.Numerical model

To investigate the face stability of shield tunnels, FLAC3D is adopted to perform the numerical simulation.The failure or unstable state of the system is determined by the unbalanced force ratio, which is defined as the ratio of the average unbalanced mechanical force to the average applied mechanical force for all the grid points at each calculation step.This unbalanced force ratio will decrease with the calculation step.A steady stable state will be achieved when the unbalanced force ratio is under a prescribed tolerance value(10-5is usually set as the default value).But when the unbalanced force ratio tends to attain a quasi-constant value that is larger than the prescribed tolerance value,it means that an unstable state or a steady plastic flow of the system is reached.In this case,the infinitely increasing displacements and failure of the numerical model will occur.

A 3D numerical model with the diameter to depth ratios C/D varying from 0.5 to 3 is built.Fig.3 presents a half model of a shield tunnel with a diameter of D = 10 m and C/D = 2.The model contains 90,675 zones and the dimensions are 3D × 4D × (2D + C) in the X, Y and Z directions to counteract the boundary effect.As for the boundary conditions, the top of the model (representing ground surface) is set free,the four sides are constrained horizontally and the bottom face is fixed.

Four types of soils are adopted and shown in Table 1.Soils are assigned with the Mohr-Coulomb failure criterion, which corresponds to a shear failure condition of a frictional soil.The Young’smodulus E = 20 MPa is used (Mollon et al., 2013; Ukritchon et al.,2017b; Huang et al., 2018b).The empirical equation K = 1 - sinφ for the earth pressure coefficient at rest is adopted to calculate the value of K, where φ is the internal friction angle of the soil.The value of the Poisson’s ratio υ is determined according to the relation υ = K/(1 + K) between the earth pressure coefficient and the Poisson’s ratio.Moreover, the shell structural element with a thickness of 0.35 m is adopted to model the tunnel lining,which is assigned with the Young’s modulus E = 35 GPa and the Poisson’s ratio υ=0.25.During the simulation process,a length of 10 m of the tunnel is excavated.After the excavation, the zone behind tunnel face is installed with lining,simulated by shell elements,instantly.

Table 1Soil conditions of the tunnel face.

Fig.2.Existing failure mechanisms: (a) Two-block, (b) Multi-block, and (c) ‘Horn’ shape.v and φ represent the velocity vector and internal friction angle, respectively.

Fig.3.Numerical model of a circular shield tunnel.

Without a priori assumption on the critical state of the tunnel face, the determination of the critical collapse pressure is proceeded by the stress-controlled method.First, the tunnel face is applied with a uniformly distributed face pressure to maintain its force balance.Second,the displacement of tunnel face is monitored with the decrease of the support pressure in each calculation cycle.Finally, when the tunnel face collapses or the plastic flow occurs after a tiny decrement of the support pressure,the face pressure at this moment is noted and considered as the critical face pressure.

2.2.Results of the numerical simulations

2.2.1.Critical face pressure

Fig.4 shows the face pressure ratio n versus the horizontal displacement of tunnel face with different relative buried depth ratios C/D and different types of soils.The face pressure ratio is defined as follows:

where n defines the ratio of the face pressure to the horizontal earth stress,is the face pressure,and γ is the unit weight of soils.

It is shown that the horizontal face displacement increases rapidly with the decrease of the face pressure ratio.The critical value of n is determined when the slope of the curve approximates 0.The dotted lines in Fig.4 represent the critical value of n.It is found that a higher critical support pressure ratio is obtained for asmaller C/D.Besides,differences between Fig.4a and b shows that a smaller critical face pressure ratio is needed for soils with a higher shear strength.It is correspondingly expected that when the shear strength of soils is higher, a smaller face pressure is required.

2.2.2.Velocity distribution at the tunnel face

As aforementioned, the translational and rotational velocity fields are two commonly used mechanisms for tunnel face problem.In order to make further investigation on failure mechanism of tunnel face,this section will perform several numerical simulations to study the velocity distribution at the tunnel face.The velocity field of the numerical simulations will be compared with the translational and rotational velocity fields usually adopted in the existing analytical approaches.The obtained velocity inclination at the tunnel face is defined in Fig.5.The velocity vectors of monitoring points distributed on the whole tunnel face are collected and recorded in the numerical simulation.The velocity inclination α represents the angle between the velocity vector of the monitoring point and the opposite excavation direction (negative Y-coordinate), which can be calculated by the equation α = arccos[-vy/(vx2+vy2+vz2)0.5].

The velocity inclinations at the tunnel face in clays and sands for different depth to diameter ratios C/D are provided in Figs.6 and 7.The velocity inclinations at the tunnel face are similar in clays and sands, both increase with the Z-coordinate.The rotational velocity field is basically consistent with the velocity inclination obtained from the numerical simulations.It is shown that the rotational velocity field well estimates the velocity distribution at the tunnel face.But the velocity inclination of the translational velocity field keeps constant with the Z-coordinate,which oversimplifies the soil movement at the tunnel face.Besides,it is found that the points of the numerical simulations on both ends of curves and the points far away from the vertical symmetric plane of the tunnel face(X ≥0.4D) deviate from the curves estimated by the rotational velocity field.This is because the velocity direction near the tunnel face edge will be inevitably affected by the face boundary in the numerical simulation,but the rotational velocity field assumes that the soils at the tunnel face move together as a rigid block.It is considered that the rotational velocity field still accurately simulates the soil movement of most parts of the tunnel face and the rotational velocity field is sufficiently reliable to represent the velocity at the tunnel face.

2.2.3.Velocity distribution at the tunnel vault

This section will present the velocity distribution at the tunnel vault obtained from the numerical simulation.The comparisons of the velocity inclination provided by the numerical simulations,the translational and rotational velocity fields are performed.The velocity inclination at the tunnel vault is defined in Fig.8.As shown in Fig.8, the velocity vectors of the monitoring points distributed on the whole horizontal plane Z1are collected and recorded in the numerical simulation.The plane Z1has a length of D and a width of D, which is located at the tunnel vault.Similarly, the velocity inclination of the monitoring point of the tunnel vault is alsodefined as the angle between the velocity vector and the opposite excavation direction.

Fig.4.Curves of the face pressure ratio versus the horizontal displacement of tunnel face:(a)c=7 kPa and φ=17?,(b)c=10 kPa and φ=25?,(c)c=0 kPa and φ=20?,and(d)c = 0 kPa and φ = 40?.

Fig.5.Schematic diagram of the velocity inclination at the tunnel face in the numerical simulation.

Figs.9 and 10 show the velocity inclination at the tunnel vault obtained from the numerical simulations for clays and sands.It is shown that the velocity inclinations at the tunnel vault in the numerical simulation generally decrease with the Y-coordinate,which is almost the same in both clays and sands.However, the velocity inclination estimated by the rotational velocity field increases with the Y-coordinate, which is totally different from the results of the numerical simulation.This difference reveals that the rotational velocity field fails to represent the soil movement with arching effect.The constant velocity inclination of the translational velocity field also falls short of representing the variation characteristic of the velocity at the tunnel vault.Besides,it is shown that the velocity inclinations of the grid points far away from the vertical symmetric plane (X > 0.5D) have a different variation trend compared with those of the other points.This situation is caused by the soil arching effect.When X/D>0.5,the grid points are beyond the soil arching and located at the region of the undisturbed soils, and thus the velocity inclinations of these undisturbed points are different from others.

2.3.Comparisons between the numerical simulations and the experimental tests

To verify the obtained numerical results, the failure patterns of three typical experimental tests (Chambon and Corte, 1994;Kirsch, 2010; Idinger et al., 2011) are adopted to compare with the failure zone obtained from FLAC3D.It is shown in Fig.11 that the failure patterns of the experimental tests are basicallyconsistent with that of the numerical simulation, especially for the shear failure band obtained from the experimental test by Idinger et al.(2011).Compared with the failure patterns obtained from the experimental tests, the numerical simulation provides a more conservative estimate of the failure range of the tunnel face.Thus, although the finite difference code FLAC cannot predict the soil behavior perfectly, the numerical simulation of FLAC3D is still a reliable tool to study the stability problem of tunnel faces.

Fig.6.Velocity inclination at the tunnel face in sands: (a) C/D = 0.5, (b) C/D = 1, (c) C/D = 2, and (d) C/D = 3.

Furthermore, Table 2 presents the comparisons of the normalized critical face pressure provided by the numerical simulation in this paper and other experimental tests.The numerical method and experimental test are inherently different in the determination of the critical state of the tunnel face failure.But it is shown that the normalized critical face pressures obtained from this study approximate the existing experimental results.Only some numerical results are slightly of larger values than the experimental results.This comparison also implies that the numerical simulation can be used to investigate the tunnel face stability problem and the numerical results are relatively effective and reasonable.

Table 2Comparisons of the normalized critical face pressure.

Table 3Dimensionless parameters Nγ, Nc and Ns.

3.Construction of the improved failure mechanism

As mentioned above, the rotational and translational velocity fields are two commonly used kinematically admissible velocity fields for the face stability analysis in frictional soils.Compared with the translational velocity field, the rotational velocity field can simulate the soil movement at the tunnel face well.However,for the region above the tunnel vault, both the rotational and translational velocity fields are different from that of the numerical simulation, which fail to reflect the characteristic of the velocity distribution.Thus, this section aims at describing a new kinematically admissible velocity field for the soil arching zone based on the results of the numerical simulations to improve the representation of the soil movement with arching effect.An improved failure mechanism is newly constructed based on this proposed velocity field using the spatial discretization technique.In the proposed failure mechanism, a perfectly plastic soil material is assumed.The Mohr-Coulomb failure criterion and the associated flow rule are adopted.Finally, based on the upper bound method of the limit analysis theorem, the critical face pressure is derived by equating the rate of external work to the rate of energy dissipation.

This paper adopts two different methods to construct the upper and lower failure mechanisms, respectively.The whole failure mechanism includes two parts, as shown in Fig.12.Zone I represents the area ahead of the tunnel face(lower part)and the Zone II represents the region above the tunnel vault (upper part).Both zones are constructed by the spatial discretization technique.The lower part of the failure mechanism is constructed based on the rotational velocity field, while the upper part adopts a new kinematically admissible velocity field to generate the 3D failure surface.

3.1.Generation of the failure mechanism in Zone I

3.1.1.Principle of the point generation in Section 1 of Zone I

As shown in Fig.13, there are 2n discretized points on the circular tunnel face in Section 1,denoted by Ajand(1 ≤j ≤n).Therotation center O of the rotational velocity field is located in the vertical symmetric plane.Each radial plane passing through the point O is named as Πj,where 1 ≤j ≤n.A1andare the two first points to generate the mechanism of Section 1.These two points and point O constitute the radial plane Π1,and the points A1andare renamed as Pi,1and Pi+1,1, respectively.A new point Pi,2in the radial plane Π2is generated from the points Pi,1and Pi+1,1of Π1.More generally, a new point Pi,j+1in the radial plane Πj+1is generated from the points Pi,jand Pi+1,jof Πj,as shown in Fig.14.The facet Fi,jconsists of the points Pi,j+1, Pi,jand Pi+1,j, and all the triangular facets constitute the 3D failure surface.The point generation process should obey the normality condition and the quasiuniformly distribution.The normality condition means that the angle between the outside normal vectorand the velocity vector is π/2+φ.The quasi-uniformly distribution requires that all the new generated points Pi,j+1are uniformly distributed in the plane Πj+1,i.e.θi,j+1=(θi,j+θi+1,j)/2,where Pi,j+1=(Cj+1,ri,j+1,θi,j+1).In Section 1, there are j-1 new generated points and two existing points.This process of point generation continues until the last plane Πnof Section 1.

Fig.7.Velocity inclination at the tunnel face in clays: (a) C/D = 0.5, (b) C/D = 1, (c) C/D = 2, and (d) C/D = 3.

Fig.8.Schematic diagram of the velocity inclination at the tunnel vault in the numerical simulation.

3.1.2.Principle of the point generation in Section 2 of Zone I

After finishing the point generation in Section 1 of Zone I,a total of n +1 points will be generated in the last radial plane Πn.These n+1 points will generate n+1 new points in the radial plane Πn+1of Section 2 of Zone I.Following the same normality condition and the quasi-uniformly distribution, the point generation will be successively conducted until the final point F, which represents the closure of the failure mechanism.

The spatial discretization technique is to overcome the shortcoming that only a portion of the circular tunnel face collapses,while the remaining area is at rest.Thus the 3D failure surface generated by the spatial discretization technique could not be described by a simple geometrical shape but a discretized surface.The accuracy of the discretized failure mechanism will highly depend on the fineness of the discretization parameters.More details on the mathematical formulation for the point generation ofSections 1 and 2 of Zone I can be referred to in Mollon et al.(2010,2011).

Fig.9.Velocity inclination at the tunnel vault in sands: (a) C/D = 0.5, (b) C/D = 1, (c) C/D = 2, and (d) C/D = 3.

3.2.Generation of the failure mechanism in Zone II

The generation of the 3D failure surface of Zone II begins with the intersection plane between Zones I and II.This section will propose a new kinematically admissible velocity field for the soil arching zone.The spatial discretization technique will also be adopted to generate the failure mechanism based on the proposed velocity field.

3.2.1.Proposed velocity field for Zone II

This section will present the new proposed velocity field above the tunnel vault.Zone II is divided into infinite vertical rigid blocks,each two adjacent rigid blocks are separated by a vertical X-Z plane.The velocity field of every rigid block is assumed as shown in Fig.15.The velocity is assumed to be independent of both the X- and Zcoordinates but vary with the Y-coordinate.The formula of the velocity in Zone II is defined as follows:

where vU(y)represents the velocity of Zones II;vUcis the velocity at the point Cs1, which is assumed to be vertically downward; and θv(y) is the angle between the velocity vector and the negative direction of the Z-coordinate.In addition, when θv(y) rotates counterclockwise, θv(y) is defined as positive, otherwise negative.The verification of normality condition satisfaction of the proposed velocity field is provided in Appendix A.

According to the assumption on vUc,the condition of θv(y)at the point Cs1can be obtained:

where ycis the Y-coordinate of the point Cs1.

As shown in Fig.16, the velocity field should not only obey the normality condition in Zone II,but also be kinematically compatible with the rotational velocity field in Zone I.Thus the following conditions are required:

where vL(y)is the velocity in Zone I along the intersection plane of Zones I and II;ω represents the angular velocity in Zone I;vLU(y)is the relative velocity between vL(y)in Zone I and vU(y)in Zone II;yOand zOare the coordinates of the rotational center in the vertical symmetric plane;and α(y)is the angle between the velocity vector vL(y) and the direction of Y-coordinate.

Combining Eqs.(2)-(4) can yield two relations between the rotational velocity field of Zone I and the proposed velocity field of Zone II.The first relation of vU(y)is built based on any point in the intersection plane of Zones I and II.Any point in the intersectionplane of Zones I and II must follow the triangular relationship as shown in Fig.16.Particularly,when y=yc,the second relation in Eq.(5) can be established.

Fig.10.Velocity inclination at the tunnel vault in clays: (a) C/D = 0.5, (b) C/D = 1, (c) C/D = 2, and (d) C/D = 3.

where vLcis the velocity of Zone II at the point Cs1, and αcis the corresponding angle between the direction of Y-coordinate and the velocity vL(y) at the point Cs1.The formulae of vLcand αcare expressed as follows:

According to Eq.(5), the expression of θv(y) can be given by

where

By substituting Eqs.(5)-(7) into Eq.(2), and with some simplification, the expression of vUcis represented as follows:

Thus Eq.(2) can be rewritten as follows:

Eqs.(2)-(9) above can ensure that the normality condition is considered in both Zones I and II.

As shown in Fig.17,the cross-section of the failure mechanism in Zone II is bounded by two curves emerging from points A and C,respectively(the tunnel vault and the right end of the intersection plane between Zones I and II).The two curves are defined as z = f1(y) and z = f2(y), respectively, and they will intersect at the final point Csnh.According to the normality condition,the following geometric conditions of the failure surface can be easily obtained:

Fig.11.Comparisons of the failure pattern between the numerical simulation and the experimental tests: (a) Comparison with the test by Chambon and Corte (1994), (b) Comparison with the test by Kirsch (2010), and (c) Comparison with the test by Idinger et al.(2011).

By the integration of Eq.(10), the two curves of z = f1(y) and z = f2(y) can be obtained as follows:

where C1and C2are two integration constants, which can be obtained by substituting two boundary conditions into Eq.(11).The boundary conditions are as follows:

where yAand yCrepresent the Y-coordinates of the points A and C,respectively.

By equating f1(y) to f2(y), the coordinates of the point Csnhare obtained:

Fig.12.Proposed 3D failure mechanism for a circular shield tunnel.

where(0,yCs,zCs)are the coordinates of the point Csnh.Similarly,the coordinates of the point Cs1are (0, yCs, 0).Zone II is uniformly divided into nhparts in the vertical direction by several horizontal planes,in which nh=(zCs/δH)and δHis the distance between each two adjacent horizontal planes.

Both the velocity fields in Zones I and II are defined by the coordinates of the point O (yOand zO).The shape of the proposed failure mechanism varies with yOand zO, and the critical failure pattern can be obtained by the optimization of yOand zO.

3.2.2.Principle of the point generation in Zone II

As shown in Fig.18, all the points Psj,1discretize the contour of the intersection plane between Zones I and II.The determination on Psj,1is performed by testing all the facets Fi,j.When bothand Pi,j+1are beyond the tunnel vault, then the point Pi,j+1will be deleted.When bothand Pi,j+1are below the tunnel vault, then the point Pi,j+1will be preserved.But whenis below the tunnel vault while Pi,j+1is beyond it, the intersection point between the linePi,j+1and the plane at the level of the tunnel vault will replace Pi,j+1using a linear interpolation.

The obtained discretized points Pi,jin the intersection plane between Zones I and II are all renamed as Psi,1,where 1 ≤i ≤nsand nsis the number of the obtained discretized points.The serial number i of the point Psi,1is ranked counterclockwise as shown in Fig.18.The angular parameter θsi,1in the plane Πs1is calculated as follows:

Fig.13.Spatial discretization technique for the generation of the mechanism in Zone I.

Fig.14.Details of the point generation in Zone I.

Fig.15.Proposed velocity field in Zone II.

Each two points Psi,1and Psi+1,1of plane Πs1can generate a new point Psi,2of plane Πs2.More generally, Psi,jand Psi+1,jof plane Πsjcan generate Psi,j+1of plane Πsj+1.Similarly,the point generation in Zone II should also respect the following three conditions:

(1) The angle between the normal vectorand the velocity vector of facet Fsi,jis equal to π/2 + φ.

(2) Psi,j+1belongs to the plane Πsj+1.

(3) The quasi-uniform distribution requires θsi,j+1=(θsi,j+θsi+1,j)/2.

In the plane Πs1(the intersection plane between Zones I and II.These nsdiscretized points will generate nsnew points in the next plane Πs2.This generation process will continue until the end of Zone II (the point Csnh).

3.2.3.Mathematical formulation for the point generation in Zone II

This section will illustrate the mathematical formulation for the point generation in detail,as shown in Fig.19.Psi,j,Psi+1,jand Psi,j+1are given by

where

Moreover, the normal vector to the facet Fsi,jand the velocity vector of the facet Fsi,jare defined as follows:

According to the velocity field in Zone II,the expressions of the unit velocity vector are as follows:

In order to obey the normality condition, the following conditions must be satisfied:

By substituting Eqs.(16)-(18) into Eq.(19), the following conditions are obtained:

By solving Eq.(20), the unit normal vectoris obtained as follows:

where the negative or positive sign of znsi,jis to ensurepointing outside.This condition can be satisfied by)>0.The intermediate parameters in Eq.(21) are given as follows:

Fig.16.Velocity relationship between Zones I and II.

Fig.17.Derivation of the failure mechanism in Zone II.

The vector from the point Csj+1to the point Psi,j+1is given by

where rsi,j+1is the norm of the vectorThe coordinates of the pointis the unit vector from the point Csj+1to the point Psi,j+1and can be expressed as follows:

By substituting Eq.(24) into Eq.(25), rsi,j+1is expressed as follows:

Thus the coordinates of the point Psi,j+1is obtained:

The point generation will successively proceed until the Zcoordinate of the newly generated point is greater than zCs.Besides, when the point Csnhis beyond the ground surface, the proposed mechanism outcrops.The intersection plane between the ground surface and Zone II needs to be determined.To determine this intersection plane, all the facets Fsi,jshould be tested.When bothand Psi,j+1are beyond the ground surface,Psi,j+1will be deleted.When both P′si,jand Psi,j+1are below the ground surface, Psi,j+1will be preserved.Whenis below the ground surface while Psi,j+1is beyond it, the intersection point between the linePsi,j+1and the ground surface will replace Psi,j+1using a linear interpolation.

3.3.Work equation

The computation of the proposed model includes the calculations on the rate of the energy dissipation and the total work rates of all the external forces of the failure mechanism.To illustrate the procedure of the computation of the proposed model,a flowchart is provided in Fig.20.First,the conditions of the soil property and the tunnel size should be given.Second, based on the associated flow rule, the failure mechanism of the tunnel face should be constructed to approximate the reality.Then, both the rate of the energy dissipation and the total work rates of all the external forces of the failure mechanism need to be calculated.Subsequently, the critical face pressure can be determined according to the limit analysis theorem.Thereafter,the obtained face pressure and failure pattern of the tunnel face should be verified by comparing with the existing studies, the numerical simulation, the model test and the field monitoring.At last,the constructed failure mechanism can be proven to be valid and reasonable if the obtained results are close to that of other studies.

According to Chen(1975),the following condition is required in the limit analysis method:

Fig.18.Spatial discretization technique for the generation of the mechanism in Zone II.

Fig.19.Details of the point generation in Zone II.

where ˙WDrepresents the rate of energy dissipation; ˙Wsand ˙Wdare the rates of energy dissipation within the soil mass and along the failure surface, respectively; and ˙WErepresents the total rate of work done by external forces, which includes the face pressure σt,the surcharge σsand the soil unit weight γ.

The rate of work is computed by the summation of the rate of work of each elementary area or volume.The calculation method of the elementary area and volume can be referred to Mollon et al.(2010, 2011).

The rate of work of the face pressure σtis calculated as follows:

where Stis the tunnel face; Stjis the elementary area of St; and Rtjand βtjare the radius and angular parameters of the barycenter of the elementary facet Stj, respectively.

The rate of work of the surcharge is calculated as follows:

Fig.20.Flowchart of the computation of the analytical model based on the limit analysis theorem.

where Ssuis the intersection plane between Zone II and the ground surface, Ssujis the elementary area of Ssu, and θvjis the anglebetween the velocity vUjof Ssujand the negative direction of the Zcoordinate.

The rate of work of the soil weight includes two parts: ˙WγI of Zone I and ˙WγIIof Zone II:

where V and Vsare the volumes of Zones I and II,respectively;Vi,jand Vsi,jare the corresponding elementary volumes of Zones I and II,respectively;Ri,jand βjrepresent the polar coordinates of Vi,j;vUi,jis the velocity of the elementary block Vsi,j;and θvi,jis the angle between the velocity vUi,jand the negative direction of the Z-coordinate.where S and Ssare the failure surfaces in Zones I and II,respectively;Si,jand Ssi,jare the corresponding elementary surfaces of S and Ss,respectively;Ri,jrepresents the radial coordinate of the elementary surface Si,j; ˙εmaxis the maximal principal strain rate; ˙εmaxi,jis the maximal principal strain rate of the elementary block Vsi,j;and vLUjis the relative velocity on the elementary surface SIIIjof the intersection plane SIIIbetween Zones I and II.The formulae for the strain rate tensor are provided in Appendix A.

By substituting Eqs.(29)-(32) into Eq.(28) and equating the total rates of external force to the total rates of energy dissipation in Eq.(28),the critical face pressure can be obtained after some simplifications:

where Nγ,Ncand Nsrepresent the dimensionless parameters of soil unit weight γ,soil cohesion c and surcharge σs,respectively.Nγ,Ncand Nsare given as follows:

The rate of energy dissipation also includes three parts:the rate of energy dissipation ˙WDIof Zone I, the rate of energy dissipation ˙WDIIof Zone II,and ˙WDIIIin the intersection plane between Zones I and II.In Zone I, a rotational velocity field is considered, thus only the soil plastic deformation ˙WdIalong the failure surface leads to the energy dissipation.In Zone II, the proposed velocity field is continuous along the Y-coordinate, thus the energy dissipation includes ˙WdIIoccurring along the failure surface and ˙WsIIundergoing within the soil mass.Besides, the energy dissipation along the velocity discontinuity plane between Zones I and II is represented by ˙WdIII.The rate of energy dissipation is provided as follows:

The critical face pressure is determined by maximizing Eq.(33)with different coordinates (xO, yO).The constrained optimization procedure is performed with the software MATLAB.The discretization parameters are selected as follows: n = 200, δβ= 0.1?and δH=0.01 m,which are the compromises between accuracy and time-cost.The maximization formula is given as follows:

Fig.21.Comparisons of the critical face pressure versus C/D for sands and clays:(a)c=0 kPa and φ=20?,(b)c=0 kPa and φ=40?,(c)c=7 kPa and φ=17?,and(d)c=10 kPa and φ = 25?.

Fig.22.Variations of the critical face pressure versus (a) internal friction angle φ and (b) cohesion c.

where these constrained conditions are set to ensure the normality condition in both Zones I and II.

4.Results and discussion

This section presents and discusses the results of the critical face pressure and failure pattern obtained from the proposed mechanism.To validate its accuracy, the proposed mechanism is compared with the numerical simulation and other studies.All the results are calculated based on the conditions of D = 10 m,γ = 18 kN/m3and σs= 0 kPa.

4.1.Comparison of the face pressure

Fig.21a and b presents the comparisons of the critical face pressure for sands of φ=20?and 40?,respectively.It is shown that the critical face pressure increases with C/D when C/D = 0-1 and then remains constant when C/D>1.This is because that the whole failure pattern is below the ground surface when C/D > 1, and the increase of the depth to diameter ratio will not influence the value of critical face pressure.In terms of the upper bound solutions,the proposed mechanism corresponds better with the numerical simulation and obviously improves the solutions obtained from the rotational and translational rigid block mechanisms adopting thespatial discretization technique (S.D.T.) presented by Mollon et al.(2010, 2011).The improvements are 6.1% and 5.8% for φ = 20?and 40?compared with Mollon et al.(2011), respectively.Compared with Mollon et al.(2010), the improvements can reach 12.4% and 9.7% for φ = 20?and 40?, respectively.With respect to other kinematic approaches of the limit analysis (Leca and Dormieux, 1990; Subrin and Wong, 2002; Mollon et al., 2009),over 15% of the improvements are achieved.Besides, it is shown that the limit equilibrium methods generally provide a solution of greater value for the critical face pressure especially for Anagnostou and Kovári (1996) and Anagnostou (2012).Since the kinematic approach of the limit analysis and the limit equilibrium method are considered in two different ways,it is suggested that both methods should be taken into account to avoid over-conservative or critical estimate of the face stability.

Fig.23.Variations of the dimensionless parameters (a) Nγ, (b) Nc, and (c) Ns.

Fig.21c and d presents the comparisons of the critical face pressure for clays of c = 7 kPa, φ =17?and c = 10 kPa, φ = 25?,respectively.The results of the proposed mechanism are basically the same as that of the numerical simulation.The proposed mechanism improves the solutions obtained from Mollon et al.(2011) by 7.4% and 5.6% for c = 7 kPa, φ = 17?and c = 10 kPa,φ = 25?, respectively.Compared with the results provided by Mollon et al.(2010),the improvements can attain 18.6%and 32.5%for c = 7 kPa, φ = 17?and c = 10 kPa, φ = 25?, respectively.Moreover, it is shown that the solutions provided by the limit equilibrium method (Anagnostou and Kovári, 1996; Anagnostou,2012; Zhang et al., 2015) generate larger values compared with the limit analysis (Leca and Dormieux, 1990; Mollon et al., 2009,2010, 2011) in the case of sands.It is implied that the influence of cohesion is more pronounced in the kinematic approach of the limit analysis.

Fig.22a and b presents the variations of the critical face pressure versus the internal friction angle φ and cohesion c,respectively.On the one hand, it can be seen that a nonlinear decreasing function between the critical face pressure σtand φ is observed.The improvement of the proposed mechanism with respect to the other kinematic mechanisms decreases with φ.This is because that the proposed mechanism improves the existing studies by providing anew velocity field above the tunnel vault.But the region of the failure mechanism above the tunnel vault decreases with φ, and thus the improvement decreases accordingly.On the other hand,the critical face pressure σtfollows a linear function relationship of c.But the slopes of the lines provided by the limit equilibrium method are different from those provided by the limit analysis.It is suggested that the cohesion c has different effects on the face pressure in different methods.Compared with the limit analysis method, the influence of the cohesion c is reduced in the limit equilibrium method (Anagnostou and Kovári, 1996; Anagnostou,2012; Zhang et al.,2015).

Fig.24.Comparisons of the velocity inclination above tunnel vault in Y-Z cross-section: (a) X/D = 0, (b) X/D = 0.08, (c) X/D = 0.24, (d) X/D = 0.39, and (e) X/D = 0.55.

4.2.Design chart for the critical face pressure

According to Eq.(33),the critical face pressure can be calculated based on three dimensionless parameters Nγ, Ncand Ns.These parameters are almost independent of cohesion c but dependent on internal friction angle φ and relative depth ratio C/D.Fig.23 provides the variations of Nγ, Ncand Nsversus C/D with different φ values.It is shown that both Nγ and Ncincrease with C/D but decrease with φ.Nγ and Ncwill remain constant when C/D > 1.Moreover, Nsdecreases with both C/D and φ.The coefficient Nsrepresents the influence of the surcharge σson the tunnel face,thusNscertainly becomes zero when the whole failure mechanism is below the ground surface and when C/D reaches a certain value.

Considering the use in the practical tunneling engineering, Nγ,Ncand Nsare provided in Table 3.This table can provide a quick calculation on the critical face pressure for different C/D and φ values.

4.3.Comparison of the velocity distribution

To validate the credibility of the proposed velocity field above the tunnel vault, the comparison of the velocity inclination above the tunnel vault in Y-Z cross-section obtained from the numerical simulation and the proposed velocity field for c=0 kPa and φ=20?is shown in Fig.24.The proposed failure mechanism above the tunnel vault extends to approximately 0.7D height above the tunnel vault(Z-direction),0.4D ahead of tunnel face(Y-direction)and 0.5D in the horizontal direction(X-direction).Five vertical cross-sections(X/D = 0-0.55) with a range of Z/D = 0-1 and Y/D = 0-0.45 are selected to calculate the velocity inclination in the numerical simulation.It is shown that the velocity inclination of the proposed velocity field is basically consistent with that of the numerical simulation.The curve of the proposed velocity field is approximately the average value of the velocity inclination of Z/D = 0-1.The proposed velocity field can be considered as a reasonable hypothetical velocity field for the region above the tunnel vault in terms of the theoretical analysis.

Furthermore,the velocity vectors at the level of the tunnel vault provided by the proposed velocity field, the numerical simulation,the existing rotational and translational velocity fields are compared in Fig.25.It is shown that the distributions of the existing rotational and translational velocity fields are totally different from that of the numerical simulation,while the proposed velocity field corresponds well to that of the numerical simulation.But it needs to be pointed out that the proposed velocity field provides a greater velocity magnitude on the far end at the level of the tunnel vault compared with the numerical simulation.It is caused by the requirement of the velocity compatibility between Zones I and II, and the velocity is assumed to keep constant in the vertical direction.Consequently, the proposed velocity field assumes a more severe situation of soil movement above the tunnel vault, which will provide a more conservative estimate of the face stability.

Fig.25.Comparison of the velocity vector at the tunnel vault.

4.4.Comparison of the failure mechanism

Fig.26a and b presents the 3D failure mechanism for sands and clays, respectively.Three perspectives are presented to show the failure pattern.The failure mechanism in sands of c = 0 kPa and φ = 20?extends to about 0.7D height above the tunnel vault (Zdirection), 0.4D ahead of tunnel face (Y-direction) and approximately 0.5D in the horizontal direction(X-direction).With respect to the failure mechanism in clays of c = 10 kPa and φ = 25?, a stronger soil condition leads to a smaller failure region.The failure region is about 0.5D height above the tunnel vault, 0.4D ahead of tunnel face and approximately 0.5D in the horizontal direction for clays of c = 10 kPa and φ = 25?.The proposed failure mechanism presents a curved arch shape above the tunnel vault rather than a cone or horn shape in the existing rotational and translational rigid block mechanism.Besides, it is clear that the presented failure mechanisms are similar to the 3D images of the failure zone, as shown in Fig.1.

Fig.27 shows the comparisons among the proposed failure mechanism, the displacement contours of the numerical simulation and the exiting rotational and translational failure mechanisms.Zhang et al.(2015)defined the failure boundary according to the displacement contour.Based on this criterion, the proposed failure mechanism and the existing rotational and translational failure mechanisms all correspond well to the numerical simulation for the region ahead of the tunnel face.However, for the region above the tunnel vault, the shape, position and height of the soil arching of the rotational and translational failure mechanisms are different from the displacement contour.It is clear that the soilarching of the proposed failure mechanism closely resembles that of the numerical simulation.

Fig.26.Layout of the proposed failure mechanism for sands and clays: (a) c = 0 kPa,φ = 20?; and (b) c = 10 kPa, φ = 25?.

Fig.27.Comparisons of the failure pattern: (a) c = 0 kPa, φ = 20?; and (b) c = 10 kPa, φ = 25?.

Fig.28.Comparisons between the proposed mechanism and the incremental displacements.

Fig.29.Comparisons between the proposed mechanism and the contours of shear strain.

Fig.30.Influence of the parameters of the Mohr-Coulomb yield criterion on the critical face pressure of the tunnel face.

Fig.28 shows the comparisons between the proposed mechanism and the incremental displacements obtained from the experimental test by Kirsch(2010).The incremental displacements can present both the moving soil body and the undisturbed soil region.The failure pattern of the tunnel face is thus easily determined from the visualization of the incremental displacement fields.It is well known that the failure pattern in the experimental test is closely related to the moment when the experimenter captures the failed soils.In this study,two cases of the failure patterns in the experimental test are compared with the proposed mechanism.Fig.28a shows the incremental displacement for the piston simulating the tunnel face that advances from 1 mm to 1.25 mm.In this case, the proposed mechanism fully describes the failure pattern provided by the experimental test, which extends to only 0.3D above the tunnel vault and 0.27D ahead of the tunnel face.But when the piston advances from 2.25 mm to 2.5 mm as shown in Fig.28b, the failure pattern provided by the experimental test propagates to about 0.6D above the tunnel vault and 0.3D ahead of the tunnel face, which is comparable to the proposed mechanism with a range of 0.4D high and 0.3D wide.In general, the proposedfailure mechanism agrees well with the incremental displacement fields.

Fig.31.Influences of the parameters of the Mohr-Coulomb yield criterion on the failure pattern of the tunnel face: (a)Influence of Cohesion, and (b) Influence of internal friction angle.

Fig.29a and b presents the comparisons between the proposed mechanism and the contours of shear strain provided by the experimental test of Idinger et al.(2011)for piston displacements of 1.5 mm and 2.5 mm,respectively.The contours of shear strain can accurately describe the shear bands of soils.It is shown that the proposed mechanism complies with the shear bands of the experimental test in both cases.Both the proposed mechanism and the experimental failure pattern are approximately 0.3D ahead of the tunnel face and 0.4D above the tunnel vault.It is implied that the shear strain contour obtained from the experimental test is suitable for determining the boundary of the failed soils.

4.5.Influence of the Mohr-Coulomb parameters

Fig.30 shows the influence of the parameters of the Mohr-Coulomb yield criterion on the critical face pressure of the tunnel.It is shown that the dimensionless critical face pressure decreases linearly with the dimensionless parameter of cohesion, and the variation gradient decreases with the internal friction angle.It is suggested that the cohesion of the Mohr-Coulomb yield criterion has a greater influence on the critical face pressure with a smaller internal friction angle.Fig.30 also implies that the critical face pressure of the tunnel face can be easily obtained from the interpolation calculation according to the linear relationship between the critical face pressure and cohesion, which provides a convenient approach for the practicing engineers to estimate the stability of the tunnel face with the given conditions.

Fig.31 shows the influence of the parameters of the Mohr-Coulomb yield criterion on the failure pattern of the tunnel face.Compared with the internal friction angle,the cohesion has a subtle influence on the failure pattern of the tunnel face as the failure zone decreases slightly with the cohesion.But the internal friction angle obviously affects both the range and position of the failure pattern of the tunnel face.With the decrease of the internal friction angle,the failure pattern of the tunnel face tends to expand in both the longitudinal and vertical directions.For the tunnel of C/D=1 in the soils havingγ=18 kN/m3and c=0 kPa,thecollapse of the tunnel facewill influence the ground surface when the internal friction angle decreases to 15?.Thus,it is considered that the internal friction angle of the Mohr-Coulomb yield criterion plays a more vital role in the performance of the proposed analytical model compared with the cohesion.Because the internal friction angle is closely related to the associated flow rule adopted in the limit analysis theorem, which inherently changes the failure pattern of the tunnel face and influences the upper solution of the critical face pressure.

5.Conclusions

Both numerical simulation and theoretical analysis were carried out to assess the face stability in frictional soils.A series of numerical simulations for different C/D and soil conditions was performed to investigate the velocity distribution of the tunnel face.According to the quantitative results of the numerical simulations,a kinematically admissible velocity field for the soil arching zone was proposed to construct a new failure mechanism based on the spatial discretization technique.Both the critical face pressure and the failure pattern were compared with the results of the numerical simulations and the existing studies to verify its accuracy.The main conclusions are given as follows:

(1) The results of the numerical simulations showed that the existing rotational velocity field can reasonably simulate the soil movement at the tunnel face.But both the rotational and translational velocity fields could not represent the velocity distribution above the tunnel vault.

(2) The proposed mechanism obviously improved the critical face pressure estimation of the existing analytical studies.The critical face pressure provided by the proposed mechanism corresponded well to that of the numerical simulation.

(3) The proposed mechanism presented a curved arch shape for the soil arching zone rather than a cone or‘horn’shape in the existing analytical methods.The proposed mechanism adequately reflect the failure patterns obtained from the numerical simulation and the experimental tests.

(4) The dimensionless parameters Nγ, Ncand Nsfor different φ and C/D values were provided to calculate the critical face pressure.Nγ,Ncand Nsall decreased with φ.C/D would only impact Nγ, Ncand Nswhen the failure mechanism outcropped at the ground surface.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors acknowledge the financial support provided by the National Natural Science Foundation of China (Grant No.51978042).

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2021.10.006.