Face stability analysis of circular tunnels in layered rock masses using the upper bound theorem

Jianhong Man,Mingliang Zhou,Dongming Zhang,Hongwei Huang,Jiayao Chen

Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education,Department of Geotechnical Engineering,Tongji University,Shanghai,200092,China

Keywords:Face stability Rock tunnel Layered rock masses Upper bound solution Hoek-Brown criterion

ABSTRACT An analysis of tunnel face stability generally assumes a single homogeneous rock mass.However,most rock tunnel projects are excavated in stratified rock masses.This paper presents a two-dimensional(2D)analytical model for estimating the face stability of a rock tunnel in the presence of rock mass stratification.The model uses the kinematical limit analysis approach combined with the block calculation technique.A virtual support force is applied to the tunnel face,and then solved using an optimization method based on the upper limit theorem of limit analysis and the nonlinear Hoek-Brown yield criterion.Several design charts are provided to analyze the effects of rock layer thickness on tunnel face stability,tunnel diameter,the arrangement sequence of weak and strong rock layers,and the variation in rock layer parameters at different positions.The results indicate that the thickness of the rock layer,tunnel diameter,and arrangement sequence of weak and strong rock layers significantly affect the tunnel face stability.Variations in the parameters of the lower layer of the tunnel face have a greater effect on tunnel stability than those of the upper layer.

1.Introduction

Layered rock masses are typical natural geological settings,accounting for approximately 66.7% of the land area.Hence,various engineering construction activities often encounter layered rock masses.A growing number of rock tunnel projects involve complex geological environments(Anagnostou et al.,2014;Chen et al.,2021;Zhao et al.,2021),most of which have tunnel faces excavated in layered rock masses (Zhang and Zhou,2017; Zhang et al.,2021a).Several studies have suggested that tunnel face failure,especially collapse,is often inextricably linked to the stratification of the layered rock masses(Arnáiz Ronda et al.,2003;Babendererde et al.,2006; Anagnostou and Zingg,2013; Zhou et al.,2021).Tunnel collapse accidents in China accounted for about 55.6% of construction accidents from 2006 to 2016,most of which were due to the instability of the tunnel face (Zhang et al.,2018).Construction experience indicates that most tunnel collapse accidents can be attributed to the instability of the tunnel face (Health and Executive,1996).When the tunnel face collapses,a large amount of rock mass flows into the tunnel,resulting in many casualties,property losses,damage to tunnel support structures(Wang et al.,2019; Huang et al.,2021),construction period delays,and other issues.Therefore,face stability analysis is of great importance for tunnel engineering in layered rock masses.

Over the past several decades,the stability of tunnel faces in homogenous rock mass/soil has been extensively analyzed by experimental tests,numerical simulations,and analytical solutions.Many scholars have studied the failure shape and mechanism of tunnel faces through the centrifuge model (Idinger et al.,2011;Chen et al.,2018;Weng et al.,2020)and 1-g physical model(Kirsch,2010; Chen et al.,2015; Liu et al.,2018).While experimental tests can accurately simulate the failure phenomenon of a tunnel face,it is difficult to test the stability of a tunnel face under complex geological conditions.Three types of numerical simulation method are widely used in the stability analysis of tunnel faces: the finite element method (FEM) (Vermeer et al.,2002; Li et al.,2009; Kim and Tonon,2010; Ukritchon et al.,2017),discrete element method (DEM) (Maynar and Rodriguez,2005; Chen et al.,2011;Jiang and Yin,2014; Zhang et al.,2021b),and FEM-DEM coupling method (Long and Tan,2020; Yin et al.,2020).These numericallybased methods can simulate complex stratum environment and changes in monitoring objects during the construction process.However,using any of these numerical methods to conduct the face stability analysis of a typical tunnel setting comes at a high computational cost,and takes a long time to complete the simulation.

In contrast,an analytical solution for face stability analysis can quickly produce design charts and can conveniently be used onsite.Analytical solutions typically have two types: the limit equilibrium method and the upper bound solution.The wedge model(Horn,1961),wedge-prism model(Perazzelli et al.,2014;Paternesi et al.,2017;Zou et al.,2019;Zhang et al.,2020a),and triangular base prism model(Oreste and Dias,2012)are common models using the limit equilibrium method,which do not consider the constitutive law of the material and have low computational accuracy.Therefore,the upper bound solution has been gradually adopted as an efficient method for tunnel face stability analysis,which has a variety of hypothesized face failure mechanisms,e.g.the multi-block failure mechanism (Chen,1975; Davis et al.,1980; Leca and Dormieux,1990; Subrin and Wong,2002; Mollon et al.,2010,2011a,b; Zou et al.,2019; Li et al.,2020; Zhou et al.,2020).Most current studies have focused on the face stability in homogeneous geomaterials.Only a few attempts have been made to account for the partial collapse of soil and the spatial variability of soil-based stratum (Pan and Dias,2018a,b; Zou and Qian,2018; Li and Yang,2020; Zhang et al.,2021c,d).To the authors’ knowledge,no analytical solutions have been proposed for the face stability analysis of tunnels in layered rock masses during excavation.Thus,this study proposes an upper bound-based analytical approach to conduct the face stability analysis of circular tunnels in layered rock masses.The Hoek-Brown yield criterion was adopted to better represent the rock mass failure behavior,and the equivalent Mohr-Coulomb model parameters were obtained for the analytical solution.The geological strength index,unconfined compressive strength,and rock material constants vary with the properties of each rock layer,resulting in segmented spiral failure surfaces.The proposed model is composed of multiple logarithmic spiral failure surfaces.

The rest of the paper is organized as follows.The nonlinear failure criterion of the rock masses and parameter equivalent method are first presented in Section 2.Section 3 shows the derivation details of the proposed analytical solution.In Section 4,the analytical approach is validated by comparison with the numerical calculation results in layered rock cases and other analytical solutions.In Section 5,the analytical approach is used to conduct a comprehensive stability analysis of circular tunnel faces in layered rock masses with different geometries and material properties.Finally,the major findings of this study and possible future studies are presented in Section 6.The calculation results obtained in this paper may be used as a reference for future tunnel excavation design.

2.Nonlinear failure criterion

2.1.Generalized Hoek-Brown criterion

Rock masses contain discontinuities of various sizes and orientations at different scales.Many geotechnical experiments have shown that almost all rock materials exhibit nonlinear characteristics at failure (Shen et al.,2012; Zhao et al.,2017).The wellaccepted Mohr-Coulomb criterion cannot adequately describe the nonlinear failure mechanism (Pan and Dias,2018b).Hoek and Brown (2019) proposed the generalized Hoek-Brown criterion,which can represent the nonlinear characteristics of rock failure and is extensively used for rocks with varying degrees of fracture.The generalized Hoek-Brown criterion reported by Hoek and Brown(2019) can be written as follows:

where σ1and σ3are respectively the maximum and minimum principal stresses at failure,and σcis the uniaxial compressive strength (UCS) of intact rock.The Hoek-Brown parameters (mb,s and n) are determined by the rock material constant mi,the disturbance coefficient d,and the geological strength index GSI by the following formulations:

2.2.Equivalent Mohr-Coulomb parameters

In the upper bound solution,the internal energy dissipation along the discontinuous surface is often calculated by the Mohr-Coulomb parameters of cohesion c and frictional angle φ.However,these two parameters are not used in the Hoek-Brown failure criterion.Hence,it is necessary to perform equivalent conversions between the Hoek-Brown parameters and the Mohr-Coulomb parameters.

There are two main ways to convert the Hoek-Brown parameters into the Mohr-Coulomb parameters (cohesion c and frictional angle φ).One of such methods is the tangent method,which can linearize the nonlinear failure criterion and obtain the equivalent parameters (Pan and Dias,2018a; Zhang et al.,2020b).The other is the direct equivalent method,which incorporates engineering experience,and has proven highly reliable when tested in practice(Wang et al.,2021).The direct equivalent method can also reflect the variation of parameters with the depth of the tunnel so that the computed result is more representative of the site condition (Hoek et al.,2002).Therefore,this study uses the direct equivalent method to derive the upper bound solution.

The equivalent cohesion c and frictional angle φ proposed by Hoek et al.(2002) are expressed as follows:

where σ3n=σ3max/σc,in which the parameter σ3maxis defined as follows:

where γ is the unit weight of the rock mass;H is the buried depth of the tunnel,being the distance between the top of the tunnel and the ground surface;and σcmis the rock mass strength,which can be written as

It is noteworthy that the equivalent method mentioned above does not consider the stratification of rock masses,especially in the solution process of Eq.(4).Layered rock masses have always been a complicated problem,thus certain simplifications are needed in the analytical solution.The simplified method is performed by assuming that the properties of all rock layers are the same when solving the equivalent parameters of a specific rock layer.

Hoek et al.(2002)converted the Hoek-Brown parameters to the Mohr-Coulomb parameters c and φ using curve fitting,which adopts the principle of balancing the areas above and below the Mohr-Coulomb plot.Therefore,the equivalent Mohr-Coulomb parameters inevitably have calculation differences.However,the exhibited calculation difference has been proven acceptable in a previous study (Wang et al.,2021).

3.The upper bound solution in layered rock masses

3.1.Analytical model description

The upper bound solution combines the kinematically admissible velocity field with the corresponding yield conditions,flow rules,and boundary conditions to study the stability of geotechnical structure.It is difficult to analyze the tunnel face stability due to the complex mechanical behavior of layered rock mass.However,Dong and Anagnostou (2013,2014) and Anagnostou et al.(2014)pointed out that in most cases,the complexity of mechanical behavior does not increase the complexity of engineering.From a practical point of view,the simplified model is still appropriate for engineering purposes.Consequently,this paper uses simplified continuum modeling to consider the stratification of rock mass in front of the tunnel face(Xu et al.,2017).The following assumptions were made to simplify the geological and mechanical complications of the problem:

(1) The stratum consists of three rock mass layers: the overburden layer,the upper tunnel face layer,and the lower tunnel face layer.The exposed tunnel face is assumed to lie within the upper and lower tunnel face layers.

(2) The excavated tunnel face within the upper and lower rock layers is assumed to have a perfect circular shape.

(3) The rock layer is assumed to be an ideal elastoplastic material that conforms to the associated flow rule.Each rock layer is homogeneous,and the dilation angle of each point along the sliding surface is equal to the equivalent frictional angle of the rock medium.

(4) The failure surface of each rock layer is continuous on the layer interface,and the associated flow rule also conforms to the layer interface.

(5) Each hypothesized sliding block is regarded as a rigid body,and the strain within each block is not considered so that energy dissipation only occurs on the failure surface.

To better analyze the influence of the stratification characteristics of rock masses on the stability of rock tunnel face,an improved failure model is proposed,with two cross-layers and one cover layer.Subrin and Wong (2002) proposed a doublelogarithmic spiral collapse model in a homogeneous rock mass configuration.Combining the double-logarithmic spiral with the approach of considering a multi-layer structure with progressive failure (Qin et al.,2017),this paper establishes a progressive collapse model of layered rock masses under two-dimensional(2D)conditions.The failure surfaces thus consist of multiple doublelogarithmic spirals,which can account for the failure of different rock layers according to their material properties(Fig.1).The failure surfaces of the proposed solution are composed of four logarithmic spirals.From a mathematical geometry perspective,each logarithmic spiral is determined by two angles with respect to the vertical direction (the starting and ending angles).However,comparing the upper and lower boundaries of the failure area,the angles corresponding to points A and C (relative to the vertical direction)can have two possible cases:θ2＜θ3and θ3＜θ2.Therefore,two collapse models are proposed,i.e.Cases 1 and 2,as shown in Fig.1.

3.2.Limit analysis and calculation method in layered rock mass

According to Fig.1,the diameter of the circular tunnel is denoted as D,the thickness of the overburden layer is denoted as H,the thickness of the upper layer of the tunnel face is denoted as h1,and the thickness of the lower layer of the tunnel face is denoted as h2.Based on the proposed failure models in Fig.1,there are four logarithmic spiral surface lines: BC,CD,DE,and AE.For the failure mechanism,each failure surface is assumed to rotate with a constant angular velocity ω around point O,where the lengths of OA,OB,OC,and OD are ra,rb,rcand rd,respectively.Meanwhile,the angles between OA,OB,OC,OD,and OE and the vertical direction are denoted as θ2,θ1,θ3,θ4,and θ5,respectively.Based on the associated flow rule of the upper bound theory of limit analysis,at any point on the logarithmic spiral failure lines,the frictional angle of each layer is the same as the angle between the velocity direction and the corresponding tangent direction of the failure lines.φtiand(i=0,1,2) are denoted as the frictional angle and the cohesion of the overburden layer,the upper layer of the tunnel face,and the lower layer of the tunnel face,respectively.As shown in Fig.2,the logarithmic spiral curves AE,BC,CD,and DE can be written as follows:

where θ is a variable denoting the rotational angle.According to the geometric relationship in Fig.1,ra,rb,rcand rdcan be expressed as the functions of θ2,θ1,θ3and θ4as follows:

Fig.1.Improved failure model composed of multiple double-logarithmic spirals:(a)Parameters of each rock strata;and(b,c)Longitudinal sections of the failure surfaces in Case 1(θ2 ＜θ3) and Case 2 (θ3 ＜θ2),respectively.

Fig.2.2D calculation models of (a) Case 1 (θ2 ＜θ3) and (b) Case 2 (θ3 ＜θ2).

The parameter θ5can also be established:

A rock tunnel face retains self-stability after excavation if the rock mass quality is good and the degree of jointing is not severe.However,when a tunnel is excavated in heavily jointed and “l(fā)owquality” rock masses,the original equilibrium state is broken.The rock mass in front of the tunnel face cannot be restored to a stable condition through stress redistribution,and tunnel collapse may occur.Hence,it is necessary to take pre-supporting measures for the tunnel face,such as advanced bolts,advanced small steel pipes,advanced small pipes,etc.This makes it necessary to determine in advance the magnitude of the force to be applied to the tunnel face.

In this study,the virtual support force determines whether the tunnel face is stable.We assume that the virtual support force σTis uniformly distributed at the tunnel face.Then,based on the upper bound theory of limit analysis,when the velocity boundary condition is satisfied,the work rate of external forces applied to the system is equal to or greater than the work rate of internal energy dissipation in the system (Chen,1975).For the calculation model shown in Fig.2,three aspects are included in this study: the work rate of gravity,the work rate of the virtual support force,and the work rate of internal energy dissipation along the failure lines.

Firstly,the work rate of the self-weight of the multiple doublelogarithmic spiral failure mechanisms is obtained from the four sliding blocks (EGD,AGD,BFC,and AFCD).We assume that the length of is rfand the angle between OF and the vertical direction is θ11.Meanwhile,the intersection of OC and AD(Case 1)or AF(Case 2)is point I,and the intersection of OD and the logarithmic spiral line AE is point G.The work rate of gravity of each part can then be calculated.

The work rate of the self-weight of regions EGD and AGD can be written as follows:

where v and Gdθ represent the tangent velocity at the center of gravity of the micro-element body and the gravity of the microelement body,respectively.

According to the geometric relationship in Fig.2,the auxiliary variables θ11and rfcan be expressed as follows:

The work rate of the self-weight of region BFC can be written as follows:

For Case 1,the region AFCD is composed of regions IDC and AFCI,for which the work rates of the self-weight are calculated as follows:

For Case 2,the region AFCD is composed of regions IFC and AICD,for which the work rates of the self-weight are calculated as follows:

Therefore,the work rate of the self-weight of the multiple double-logarithmic spiral failure mechanisms can be obtained:

The expressions of the parameters in the above equations are shown in Appendix A.It is worth noting that,when calculating the work rate of a specific region,the parameters used are only related to the rock layer where the calculation region is located.

Secondly,the work rate of the virtual support force on the tunnel face can be written as follows:

where vTrepresents the tangent velocity of the micro-element body around point O for line AB,and lTdθ represents the unit length of line AB.

Finally,the internal energy dissipation along the failure lines can be represented by the energy dissipation along the failure surfaces:

where vDrepresents the tangent velocity of the logarithmic spiral micro-element body around point O,and lDrepresents the arc length of the logarithmic spiral micro-element body.

Therefore,the total internal energy dissipation can be written as follows:

Equating the total work rate of external forces to the total internal energy dissipation rate,i.e,the following expression of virtual support force can be obtained:

The virtual support force σTcan be solved according to the upper bound theorem of limit analysis.As shown in the objective function of Eq.(18),σTis a function of the variables θ1,θ2,θ3,and θ4(other variables can be represented by these four variables),namely,σT=f(θ1,θ2,θ3,θ4).The function can then be transformed into a mathematical optimization problem,which can search the best upper bound solution by optimizing the objective function of Eq.(18) under the following constraints:

Several efficient nonlinear optimization methodologies have been developed in recent years,such as the sequence quadratic iterative algorithm,second-order cone programming,and genetic algorithm (Lysmer,1970; Zhang et al.,2015).The sequence quadratic iterative algorithm and second-order cone programming can efficiently solve the optimal solution,but when faced with complex nonlinear problems,it is easy for them to fall into a local optimal solution,thus missing the global optimal solution.However,the genetic algorithm uses a probabilistic mechanism for global search to effectively avoid such problems,and finds the global optimal solution in the sense of probability.The genetic algorithm is based on biological evolution,which has the advantages of good convergence,short computation time,and high robustness.Therefore,this study uses the genetic algorithm to find the optimal combination of θ1,θ2,θ3,and θ4for the global minimum solution of the objective function (virtual support force) using MATLAB software.The flow diagram for performing the stability analysis of the tunnel face in this paper is illustrated in Fig.3.

4.Validation of analytical model

Fig.3.Flowchart for optimization of the tunnel face stability analysis.

The proposed model is simplified to consider a homogeneous rock layer.To validate the proposed model,the results obtained using the same parameters of the proposed method are compared with the analytical results obtained by Senent et al.(2013) (see Table 1).The differences in computed face stability caused by the conversion of Hoek-Brown parameters to their equivalent Mohr-Coulomb parameters are shown in Table 1.Those results are obtained by the limit analysis and FLAC3D(fast Lagrangian analysis of continua in 3 dimensions) simulations performed by Senent et al.(2013).It is noteworthy that the simulation results of FLAC3Ddirectly use the Hoek-Brown parameters.There are differences between the analytical results in this paper and those given by Senent et al.(2013),which are acceptable according to the comparison conducted by Zhang et al.(2020b)and Senent et al.(2013).This also shows that it is feasible to use the parameter conversion in this study.

The face stability analysis in the layered rock mass is compared with the numerical analysis using the finite element software ABAQUS.Fig.4 illustrates the model geometry used for numerical simulation of the tunnel face stability,which only shows half the symmetrical model.The numerical model of the rock runnel has a diameter of 10 m and a buried depth of 30 m.The size of the numerical model is taken as 40 m in the x-direction,50 m in the ydirection,and 70 m in the z-direction.The model contains approximately 51,730 zones and 57,429 nodes.Hence,boundary effects can be avoided for these dimensions.The model’s boundary conditions are given by fixed displacements at the bottom and the lateral perimeter of the model.

A linear perfectly elastoplastic constitutive model based on the Mohr-Coulomb failure criterion (with an associated flow rule) is applied to the rock masses.The process of collapse simulation is implemented according to the sequence of in situ stress balancerock mass excavation-virtual support force-gradually decreasing the support force until the tunnel face collapses.The detailed steps were described by Vermeer et al.(2002),and the equivalent cohesion c and frictional angle φ were obtained by Hoek et al.(2002).The corresponding rock Young’s modulus can also be obtained through the alternative equation proposed by Hoek and Diederichs (2006).The tunnel concrete lining is simulated with shell structural elements with a Young’s modulus of 10 GPa,a Poisson’s ratio of 0.2,a density of 2500 kg/m3,and a thickness of 0.4 m.

In this section,the tunnel face of the calculation model is assumed to be in a two-layer rock mass,with the thickness of each layer being 5 m.The numerical analysis considers two types of rock masses: very poor rock mass (Case A) and relatively better rock mass (Case B).The material properties are listed in Table 2.The equivalent cohesion c,equivalent frictional angle φ,and the corresponding equivalent Young’s modulus proposed by Hoek et al.(2002) and Hoek and Diederichs (2006) are also shown in Table 2.The values are reserved to three decimal places to ensure the accuracy of the equivalent parameters as much as possible.

Table 1 Comparison of collapse pressure obtained in this work with analytical results gained by Senent et al.(2013).

In the numerical simulation,the minimum required virtual support force can be determined by the mutation position of the horizontal displacement (Vermeer et al.,2002; Li et al.,2009;Ibrahim et al.,2015).The results of the numerical simulation and the analytical analysis of this paper for these two cases are shown in Fig.5a.The results calculated by the proposed method are close to those of the numerical simulation.The geometry of the sliding surface was also captured in ABAQUS for Case B.We can see that the failure surface prediction of the upper bound solution is similar to that of the numerical model (see Fig.5b).Overall,the method in this paper has high computational efficiency,needing only 40-50 s,while the ABAQUS model needs several hours.Therefore,the calculation method is more suitable for construction applications.

Fig.4.Numerical model.

5.Design charts for face stability

Layered rock masses usually have different mechanical properties due to the existence of structural planes in the rock masses(Yang et al.,2018;Du et al.,2019;Wang et al.,2020),thus the tunnel face stability in layered rock is different from that in homogeneous rock.This section studies the influence of the thickness of the rock layer and the parameters of the rock layer at different positions on the stability of the tunnel face in the layered rock mass,and provides the corresponding design charts for engineering reference.

The basic parameters are as follows:the diameter of the tunnel D=10 m,the thickness of the overburden H=20 m,and σci=1 MPa,mi=5,GSIi=10,di=0,and γi=25 kN/m3(very poor rock mass),in which i=0,1,2 corresponds to the overburden,the upper layer of the tunnel face and the lower layer of the tunnel face,respectively.The detailed properties of each rock layer are shown in Table 3.The thickness of the upper layer of the tunnel faceh1and the lower layer of the tunnel face h2are determined by the computational cases.

Table 2 Parameters of the calculation model.

Table 3 Basic parameters of each rock layer.

5.1.Design charts with two different collapse models (Cases 1 and 2)

This section investigates the tunnel face stability of the two collapse models (Cases 1 and 2 in Fig.1) and determines the less favorable model geometry.The thickness of the upper and lower rock layers of the tunnel face is set to 5 m.The value ranges of the rock mass parameters are set as follows: m1=5-25,GSI1=5-25,and σc1=5-25 MPa.The normalized virtual support force is defined as σT/(γ1D).As shown in Fig.6,the stability of the tunnel face gradually increases as the quality of the rock mass increases.Hence,the influence of the rock mass quality on the tunnel face stability is similar for the two models.However,the failure model corresponding to Case 1 is prone to collapse compared with Case 2.Consequently,the subsequent analysis shall focus on the less stable model geometry (Case 1).

5.2.Design charts with different rock layer thicknesses

This section investigates the influence of rock layer thickness on the tunnel face stability,considering changes in the geological strength index GSI1,rock material constant m1,and unconfined compressive strength σc1.The thickness ratio is defined as α=h1/D to characterize the change in the thickness of the upper layer.As shown in Fig.7,the normalized virtual support force decreases nonlinearly as m1,GSI1,and σc1increase.In contrast,the nonlinear decrease in the virtual support force weakens as the thickness of the upper layer on the tunnel face decreases.Fig.7 also shows that the curves of the normalized virtual support force with different thicknesses of the upper layer have an intersection point.To the left side of the point,the greater the thickness of the upper layer,the more unstable the tunnel face,whereas the situation to the right side of this point is just the opposite.The average gradient parameter ?σT/?h1is then defined to characterize the influence of the upper layer thickness under different rock parameters on the stability of the tunnel face.As shown in Fig.8,the average gradient gradually decreases from positive to negative as the values of m1,GSI1and σc1increase.This further indicates that the higher the quality of the rock layer,the greater the thickness of the rock layer,and the more conducive to the stability of the tunnel face.In addition,the contrast between the three curves shows that change in the geological strength index GSI1,which characterizes the joint characteristics of the rock mass,has a more significant effect compared to the other two parameters on the stability of a tunnel face in a poor rock mass.However,for a better-quality rock mass,a change in unconfined compressive strength σc1,which characterizes the strength of the rock mass,has a more significant effect on the stability of a tunnel face compared to the other two parameters.

5.3.Design charts with different tunnel diameters

To characterize the influence of tunnel diameter changes on the stability of a tunnel face,Fig.9 presents the influence of different rock layer parameters changing with the thickness ratio on the stability of the tunnel face for three different diameters.Transverse comparison of the subfigures in Fig.9 shows that as the tunnel diameter increases,the stability of the tunnel face deteriorates,and the increase of thickness ratio α has a more significant influence on the stability of the tunnel face,regardless of m1,GSI1and σc1.Vertical comparison between the subfigures in Fig.9 shows that in the cases with a poor-quality upper layer rock mass(the black lines in Fig.9),the parameter GSI1is more sensitive to the change of thickness ratio.However,as the upper layer rock mass quality improves,the parameter σc1becomes the most sensitive to changes in the thickness ratio,and this influence becomes greater as the tunnel diameter increases.Therefore,it can be concluded that the stability of a tunnel face is significantly affected by the size of the tunnel diameter and the distribution of the rock layer thickness.

Fig.5.(a) Comparison of numerical simulation and analytical results for two cases; and (b) Comparison of failure mechanisms computed with the upper-bound solution and the numerical model (unit: m).

5.4.Design charts with different rock layer parameters

In this part,the thicknesses of the upper and lower rock layers of the tunnel face are set to 5 m to eliminate the influence of the thickness of the rock layers.The following analysis investigates the influence of the geological strength index GSI and the rock material constant m on the stability of the tunnel face when the other parameters are kept constant.The normalized virtual support force σT/(γiD)is plotted against the normalized UCS σci/(γiD).As shown in Fig.10,the normalized virtual support force declines nonlinearly as σci/(γiD) increases from 4 to 100,and the nonlinear behavior gradually weakens.The first column in Fig.10 shows that,as the parameter σci/(γiD) increases,the influence of changes in the parameter mion the virtual support force gradually weakens.However,for the geological strength index GSIi,the second column of Fig.10 indicates that the influence of changes in GSIion the virtual support force gradually increases as σci/(γiD) increases.More importantly,the comparison shows that the changes in the parameters of the lower rock layer have a greater influence on the tunnel face stability.Hence,the rock quality of the lower layer controls the stability of the tunnel face.

This section also investigates the influence of the rock layer arrangement on the tunnel face stability.It needs to be explained here that strength and weakness here are relative.Two types of layer arrangements are investigated: one with m1=5-25 and m2=15,and another with m2=5-25 and m1=15.As shown in Fig.11,to either the left or right of the intersection point,the stability of the tunnel face with the upper-weak and lower-strong arrangement is more stable than that of the lower-weak and upper-strong arrangement.Therefore,the arrangement of the weak and the strong rock layers has an important influence on the stability of the tunnel face.In practice,it is necessary to judge the stability of a tunnel face according to the arrangement of weak and strong layers.

Fig.12 shows the results obtained from 820 calculated examples.It was assumed that the tunnel face’s upper and lower rock layer parameters change simultaneously,mainly because the rock masses are cut into finite blocks by a weak plane (persistent joint plane).Note that selection of the parameters primarily corresponds to the poor-quality rock masses,where face instability problems are more likely to be found in engineering practice.Accordingly,σc/(γD)values vary from 4 to 100(equivalent to σcbetween 1 MPa and 25 MPa),typical of soft rocks to very soft rocks.GSI values are taken between 10 and 30,which are characteristic values of very poor-quality rock masses.As shown in Fig.12,the virtual support force gradually decreases to zero when the rock mass quality increases,which means that the tunnel face can maintain a selfstable state.However,comparison with Fig.10 shows that the stability of the tunnel face will be overestimated if the stratification of the rock mass and the differences in the properties of the upper and lower rock layers are not considered.Therefore,the calculation method proposed in this paper is more suitable for determining the stability state of a tunnel face in layered rock.

Fig.6.Influence of varying parameters on the tunnel face stability under two collapse models: (a) m1,(b) GSI1,and (c) σc1.

Fig.7.Influence of thickness on the virtual support force with variations in (a) m1,(b) GSI1,and (c) σc1.

6.Conclusions

Fig.8.The average gradient varying with rock parameters.

This study extends the upper bound limit analysis method to the face stability of a tunnel excavated in layered rock.The influence of rock layer characteristics on the stability of the tunnel face is considered,and a series of dimensionless parameter charts is presented using the upper bound limit analysis method.These charts are useful for evaluating the stability of a tunnel face in a layered rock mass.The results of this study are as follows:

(1) For layered rock masses with different layer thicknesses,the geological strength index GSI1,rock material constant m1,and UCS σc1significantly influence the stability of tunnel face characterized by the normalized virtual support force,which varies nonlinearly with increases in GSI1,m1and σc1.The nonlinear decreasing trend weakens with a decrease in the thickness of the upper layer of the tunnel face.As the rock mass quality improves,the key factor that controls the tunnel face stability changes from GSI1to σc1.

(2) The tunnel diameter and the distribution of the rock layer thickness significantly influence the stability of tunnel face.The larger the tunnel diameter,the greater the proportion of heavily fractured rock masses,the more unfavorable the stability of the tunnel face.The arrangement of the weak and strong rock layers has an important influence on the tunnel face stability.The research found that the tunnel face in a lower-weak and upper-strong composite layer fails more easily than that in an upper-weak and lower-strong composite layer.

Fig.9.Influence of parameter variations on the stability of tunnel face under different tunnel diameters.

(3) The influence of the parameters at different positions on the tunnel face stability was also investigated.As expected,the normalized virtual support force declines nonlinearly with an increase in σci/(γiD).The nonlinear decreasing trend weakens as the rock quality increases.In addition,the face stability is more sensitive to the variations in the lower layer parameters than that of the upper layer.As the parameter σci/(γiD) increases,the influence of increases in mion the tunnel face stability gradually weakens.However,the influence of increases in GSIion the face stability gradually increases with increases in σci/(γiD).

Fig.10.Influence of the rock material constants for a layered rock on the virtual support force: (a) m0; (b) GSI0; (c) m1; (d) GSI1; (e) m2; and (f) GSI2.

Fig.11.Influence of the arrangement of strong and weak rock layers on the tunnel face stability.

In summary,the method presented in this study can address stability in terms of the stratification characteristics of the tunnel face.Considering the cross-effects of rock parameters,the normalized virtual support force can be directly used to guide tunnel design.

It should be noted that the actual rock tunnel has a threedimensional (3D) arching effect,which inevitably leads to differences between the proposed 2D analytical solution and the actual tunnel face stability.As shown in Fig.5a and Table 1,the proposed 2D analytical solution overestimates the virtual support force of a tunnel face for extremely fractured rock masses,which may be associated with ignoring the 3D arching effect.Further studies should be conducted to derive an analytical solution that considers the 3D boundary condition with mechanical equivalence.The proposed solution in this article is an exploratory work to analyze tunnel face stability in stratified rock masses.In addition,due to the complexity of layered rock masses,several assumptions were adopted to derive the 2D analytical solution,such as ignoring the effect of stratification of the rock masses when the parameters are equivalent,assuming a progressive collapse model at failure,and so on.In-depth study shall be undertaken in the future to address those less realistic assumptions.

Fig.12.Stability of tunnel face under different parameters.

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 work described in this paper is supported by the Key Innovation Team Program of Innovation Talents Promotion Plan by MOST of China (Grant No.2016RA4059) and the Science and Technology Project of Yunnan Provincial Transportation Department (No.25 of 2018).

Appendix A.Supplementary data

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