Numerical analyses for three-dimensional face stability of circular tunnels in purely cohesive soils with linearly increasing strength

The tunnel face stability in purely cohesive soils with linearly increasing strength was investigated using three-dimensional finite element limit analysis (FELA). Both the collapse (active failure) and blow-out (passive failure) of the tunnel face were considered in the analysis. The rigorous upper bound (UB) and lower bound (LB) solutions of the load factor were calculated with a wide range of ground conditions to cover a broad scope of practical application. The results showed that the whole surface of the face is at failure in the collapse case; while in the blow-out case, there exists a gradual evolution process from partial failure to global failure within the tunnel face with increasing buried depth. Later, based on 960 finite element limit analysis results, a series of practical equations were proposed for tunnel face stability analysis in purely cohesive soils. These equations can be employed to quickly calculate the UB and LB solutions of the limit support pressure and the stability number of a tunnel face in both the collapse and blow-out cases. Finally, the calculation results from these equations were compared with those from previous studies in detail. The comparisons showed that the proposed equations make an improvement over existing methods and can be used as an efficient tool in practical engineering.

The tunnel face stability in purely cohesive soils with linearly increasing strength was investigated using three-dimensional finite element limit analysis (FELA).Both the collapse (active failure) and blow-out (passive failure) of the tunnel face were considered in the analysis.The rigorous upper bound (UB) and lower bound (LB) solutions of the load factor were calculated with a wide range of ground conditions to cover a broad scope of practical application.The results showed that the whole surface of the face is at failure in the collapse case; while in the blow-out case, there exists a gradual evolution process from partial failure to global failure within the tunnel face with increasing buried depth.Later, based on 960 finite element limit analysis results, a series of practical equations were proposed for tunnel face stability analysis in purely cohesive soils.These equations can be employed to quickly calculate the UB and LB solutions of the limit support pressure and the stability number of a tunnel face in both the collapse and blow-out cases.Finally, the calculation results from these equations were compared with those from previous studies in detail.The comparisons showed that the proposed equations make an improvement over existing methods and can be used as an efficient tool in practical engineering.In shield tunnelling, one of the most crucial issues is to prevent the tunnel face from failure.Specifically, the collapse (i.e., active failure) of the tunnel face occurs when the support pressure is not sufficient to prevent the movement of the soils toward the tunnel.The blow-out (i.e., passive failure) of the tunnel face occurs when the support pressure is high enough to push the soils toward the ground surface 1 .Therefore, the determination of reasonable support pressure is a crucial problem in practical engineering.To solve this problem, a series of research have been conducted to investigate the tunnel face stability in frictional, frictional-cohesive, and purely cohesive soils, which can be mainly divided into three aspects: experimental, numerical, and theoretical methods.

Abbreviations
In frictional and frictional-cohesive soils, several centrifuge model tests [2][3][4][5][6] , and 1 g model tests [7][8][9] have been performed to study the limit support pressure and failure mechanism of the tunnel face.The test results revealed the failure mechanism and provided a reference for the verification of numerical and theoretical methods.Over the last decades, numerical simulations have become a common tool.Many numerical models have been established to assess the tunnel face stability using the finite element method (FEM) 10,11 , finite difference method (FDM) 12,13 , discrete element method (DEM) 14,15 and finite element limit analysis (FELA) 16,17 .The failure mechanism of the tunnel face was investigated in detail, and a series of equations to assess the tunnel face stability were proposed.However, the numerical methods were time-consuming, therefore, many analytical models have been developed for practitioners to efficiently analyze the tunnel face stability.Horn 18 proposed the wedge-silo model, and the first practical application of this model was conducted by Jancsecz and Steiner 19 based on the limit equilibrium method.Later, this model was improved to be more accurate [20][21][22] and modified to consider the effects of the heterogeneity of the ground 23 , reinforcements 6,24,25 , unsupported span 26 and seepage flow 27 .However, even the wedge-silo model has been widely employed owing to its simplicity, the necessary assumptions (especially the simplified failure mechanism and the lateral pressure coefficient) make this model not rigorous enough.Another commonly used method in tunnel face stability analysis is the limit analysis method 28 (including upper and lower bounds theorems), which is more rigorous than the limit equilibrium method.Based on the upper bound theorem of limit analysis, a series of translational rigid-block models [29][30][31][32] and rotational rigid-block models 33,34 have been established.The results from these analytical models, especially those established in the last decade, agree well with those from experimental tests and numerical simulations, which indicates that the failure mechanism employ in these models is reasonable.
In purely cohesive soils, experimental tests [35][36][37][38][39] also have played essential roles in tunnel face stability study.The results showed that compared with the failure mechanism of the tunnel face in frictional and frictionalcohesive soils, the failure mechanism in purely cohesive soils do not have apparent localized shear band, and the soil moves more like a 'flow' rather than a rigid block.This conclusion also indicated that the rigid-block models that already exist for frictional and frictional-cohesive soils maybe not suitable for purely cohesive soil.In terms of numerical simulations, a series of numerical models have been established to study this problem [40][41][42][43] .It is worth noting that the equations proposed by Ukritchon et al. 41 based on FEM results provide the optimal solution so far.However, in FEM, the limit state is obtained by the so-called stress-controlled approach instead of rigorous analysis.Besides, Ukritchon et al. 41 only considered the collapse case and neglected the blow-out case.In terms of theoretical methods, Broms and Bennermark 38 defined the stability number N based on field observation data collected during tests, and N has been widely employed in practical applications.Davis et al. 44 first proposed upper and lower bound solutions for tunnel face stability analysis in undrained clay.Ellstein 45 obtained an analytic solution employing the limit equilibrium method, and the calculation results were in close agreement with the experimental results reported by Kimura and Mair 39 .Klar et al. 46 established two new upper bound solutions for 2D and 3D tunnel face stability analysis based on a continuous deformation of the soil.Later, Mollon et al. 1 suggested two continuous velocity fields for the collapse and blow-out of a pressurized tunnel face in purely cohesive soils, and these velocity fields significantly improve the best existing upper bound solutions.More recently, this failure mechanism was extended for tunnel face stability analysis in nonhomogeneous undrained clay 47,48 , horseshoe-shaped shield tunnels 49,50 , and volume loss prediction for tunnel face advancement 51 .However, there are still some non-ignorable differences between the results calculated by the most advanced analytical models and those from numerical simulations 1,48,50 , which indicates that the failure mechanism employed in these analytical models need to be further improved.Besides, the partial failure mechanism in the blow-out case reported by several authors 12,52,53 was neglected in these analytical models.
To solve these issues, this paper aims at proposing a quick and accurate method for tunnel face stability analysis in purely cohesive soils in both the collapse and blow-out cases."Problem definition and method of analysis" describes the problem of the tunnel face stability and introduces the FELA."Results and discussions" plots and discusses the calculation results of the load factor and failure mechanism."New design equations" proposes a series of equations for tunnel face stability analysis in purely cohesive soils."Comparisons with previous studies" compares the calculation results of the proposed equations with those from the existing methods.

Problem definition and method of analysis
The problem definition of the present study is shown in Fig. 1, and only half of a tunnel is considered due to the problem symmetry.The tunnel has a diameter D and a cover depth C in purely cohesive soils.The soil is modeled as an elastic-perfectly plastic Tresca material.The undrained cohesion c u of soil is linearly increasing with depth y, i.e., c u = c u0 + ρy, where c u0 is the undrained cohesion of soil at the ground surface, and ρ is the variation gradient of undrained cohesion.Besides, the tunnel face is applied by a uniform support pressure σ t , and the ground surface is loaded by a uniform surcharge σ s .
Most of the available 3D numerical simulations such as finite element method (FEM), finite difference method (FDM) and discrete element method (DEM) are considered very time consuming for simulating the failure of a 3D tunnel face.The new development of 3D FELA technique is a powerful and efficient tool in calculating the rigorous upper bound (UB) and lower bound (LB) solutions of stability problems in geotechnical engineering.
This numerical limit analysis method assumes a perfectly plastic material with associated flow rule, and combines the classical plasticity theorems, finite element method and mathematical optimization method 54 .In addition, the FELA does not need to make any assumptions about the failure mechanism, which allows FELA to get a more rigorous result than analytical limit analysis methods.In this work, the advanced three-dimensional commercial FELA software OptumG3 55 was employed to calculate the rigorous UB and LB solutions of the limit collapse and blow-out pressures σ t of a tunnel face.Figure 2 shows an example of the mesh and boundary conditions for the 3D stability analysis of a circular tunnel face in C/D = 2.The adaptive meshing method improve the solution accuracy as the mesh density is the greatest in zones with significant plastic shear strains, which is achieved by gradually optimizing the mesh size in zones with significant plastic shear strains.So, an automatically adaptive meshing method was employed herein, and four iterations of adaptive meshing with the initial number of elements increasing from 3000 to 10,000 were used in all analyses based on Shiau and Al-Asadi's work 56 .In addition, using the load multiplier method in FELA, the rigorous UB and LB solutions of the limit pressure can be easily calculated.

Load factor (σ s − σ t )/c u0
For tunnel face stability analysis in purely cohesive soils, there are seven calculation parameters: cover depth C, tunnel diameter D, support pressure σ t , surface surcharge σ s , unit weight of soil γ, undrained cohesion of soil at ground surface c u0 , and variation gradient of undrained cohesion ρ.According to previous studies, four dimensionless parameters were often used in undrained stability analysis 41,42 : cover depth ratio C/D, overburden stress factor γD/c u0 , linear strength gradient ratio ρD/c u0 , and load factor (σ s − σ t )/c u0 .The latter parameter is considered as a dependent parameter, and the other three parameters are regarded as independent parameters.
Nine hundred and sixty limit analyses were conducted to investigate this problem, in which C/D ranges from 0.25 to 5, γD/c u0 ranges from 0 to 10, and ρD/c u0 ranges from 0 to 1.These values were selected to cover a wide range of practical engineering conditions.Figure 3 and Table A1 shows the UB and LB solutions of www.nature.com/scientificreports/with different C/D, γD/c u0 , and ρD/c u0 .It can be observed that, in the collapse case, the value of (σ s − σ t )/c u0 linearly decreases with increasing γD/c u0 and linearly increases with increasing ρD/c u0 .However, in the blow-out case, the value of (σ s − σ t )/c u0 linearly decreases with increasing γD/c u0 and ρD/c u0 .This because under undrained conditions, the strength of soil is mainly determined by the undrained cohesion c u , while the friction angle of soil has no effect on the strength.And these findings align well with existing models 41,42 .Besides, the relationship between (σ s − σ t )/c u0 and C/D are nonlinear in both the collapse and blow-out cases.

Failure mechanism
Figure 4a-l show the shear dissipation contours in limit states of tunnel face collapse and blow-out for C/D = 0.25-3.0,which were output from FELA.The soil parameters γD/c u0 = 8 and ρD/c u0 = 0 were selected due to the limited space of the publication.Note that, in each figure, the left side is the lateral view of the tunnel, and the right side is the front view of the tunnel face.Apparently, the failure mechanisms in the collapse and blowout cases have both similarities and differences.
In the collapse case, the failure zone in front of the tunnel face resembles a chimney-like mechanism, which complies with the existing numerical and experimental results 37,50 .In addition, the failure mechanism stretches out from the tunnel invert and then propagates to the ground surface, and the whole surface of the face is at In the blow-out case, the failure zone in front of the tunnel face is also a chimney-like mechanism.However, it is worth noting that for C/D < 2.0, the failure mechanism starts from a position above the tunnel invent instead of the tunnel invert, and only the upper part of the tunnel face is concerned by failure.This failure mechanism can be called the partial failure mechanism.Besides, the global failure mechanism can be observed in C/D > 2.0.This phenomenon indicates that in the blow-out case, there exists a gradual evolution process from partial failure www.nature.com/scientificreports/ to global failure within the tunnel face with increasing C/D.This partial failure phenomenon has been reported by many authors 12,52,53 .Still, the existing analytical models for tunnel face stability analysis in purely cohesive soils do not consider this mechanism, which may lead to inaccurate calculation results.

New design equations
In practical engineering, it is necessary to evaluate the stability of the tunnel face quickly and accurately.However, numerical simulation usually requires a high-performance computer to achieve high computational efficiency, while analytical models are not rigorous enough owing to some essential assumptions (especially the failure mechanism).For the sake of application, two practical equations are presented based on the calculation results listed in Table A1 using the nonlinear regression analysis, as shown in Eqs.(1a) and (1b).
where coefficients a 1 -a 6 for different conditions are listed in Table 1.
As shown in Fig. 5a-d, the results calculated by Eqs.(1a) and (1b) are compared with those from FELA.As expected, these two methods are in excellent agreement with a coefficient of determination (R 2 ) of 99.96% (in the collapse case) and 99.99% (in the blow-out case).
Generally, the stability of a tunnel face is evaluated by the limit support pressure.Therefore, the equations that can be employed to calculate the UB and LB solutions of limit collapse and blow-out pressures are presented by adjusting Eqs.(1a) and (1b), as follows: where N c0 , N cρ , and N γ are the stability factors that respectively represent the effect of undrained cohesion c u0 , variation gradient of undrained cohesion ρ and unit weight of soil γ on tunnel face stability: In addition, a traditional method to empirically assess the stability of a tunnel face in purely cohesive soils is based on the stability number N that proposed by Broms and Bennermark 38 , as shown in Eq. ( 3).Substituting Eqs.(2a) and (2b) into Eq.( 3) results in two equations of stability number N:

Comparisons with previous studies
To verify the present equations, the calculation results from proposed equations are compared with those from previous studies.These include three terms: stability number N, stability factors N c0 , N cρ , N γ , and limit support pressure σ t .

Comparisons of stability number N
The ground is assumed to be homogeneous in most previous studies, therefore, ρ equals 0 in this section for the sake of comparison with existing methods.Since the stability factor N γ in Eq. (2e) is not equal to (C/D + 0.5), the value of stability number N dependents on γD/c u0 .To consider the effect of γD/c u0 , three values of γD/c u0 = 0, 5, and 10 are considered in Fig. 6.Note that the equations proposed by Ukritchon et al. 41 were derived from the results of FEM, and the analytical models 1,31,33,44,46,57 were established based on the upper bound theory of the limit analysis method.Figure 6 shows that the effect of γD/c u0 on stability number N is not apparent, and the equations established in this study make a significant improvement over existing methods.This conclusion can be explained as that the analytical models need to presuppose a failure mechanism, but the failure mechanism employed in these methods maybe not the optimal one; in FEA, the limit support pressure is obtained by the  7a, the results of N γ calculated by all methods are in close agreement.Then, it can be observed in Fig. 7b,c that the analytical models 1,47,50 significantly overestimate N c0 and N cρ , and the equations proposed by Ukritchon et al. 41 provide slightly larger estimates on N c0 and N cρ .

Comparisons of limit support pressure σ t
Figures 8 and 9 present the comparison in terms of the limit collapse and blow-out pressures for a tunnel of D = 10 m in homogeneous (ρD/c u0 = 0) and nonhomogeneous (ρD/c u0 > 0) soils with a unit weight γ = 18 kN/ m 3 .It is evident that the present study provides the best result, because the limit collapse pressure calculated by Eq. (2a) is the largest of all methods, while the limit blow-out pressure calculated by Eq. (2b) is the smallest of all methods 1 .This because the partial failure mechanism is not considered in these methods.To be specific, the limit collapse pressure calculated by Mollon et al. 1 agrees well with those from Eq. (2a), but the agreement decreases in the blow-out case.Besides, the inhomogeneity of soil is neglected in this model, which may limit its application.Meanwhile, other analytical models significantly underestimate the limit collapse pressure and overestimate the limit blow-out pressure, which indicates that the failure mechanism employed in these models can be further improved.As shown in Figs.8b and 9b, the inhomogeneity of soil has a significant influence on tunnel face stability.The limit collapse pressure decreases with the increases of ρD/c u0 , while the limit blow-out pressure increases with increasing ρD/c u0 .Therefore, the inhomogeneity of soil cannot be ignored in practical engineering.Also, in the collapse case, the results obtained by the equations proposed by Ukritchon et al. 41 slightly underestimate the limit support pressure.However, in the blow-out case, the difference between Eq. (2b) and Ukritchon et al. 's equations are visible.Additionally, as shown in Fig. 9, because the partial failure mechanism is considered in this paper, the results of the limit blow-out pressure is smaller than those of other models.So, if this partial failure mechanism is not considered, the limit blow-out pressure of the tunnel face will be overestimated.

Application
A proposed Yakang highway tunnel has a diameter (D) of 10 m and a cover depth (C) of 5 m.The soil is found to be purely cohesive soil with a unit weight (γ) of 18 kN/m 3 , undrained cohesion at ground (c u0 ) of 20 kPa, variation gradient of undrained cohesion (ρ) of 0.4.Determine the limit support pressure (σ t ) when the surcharge pressure (σ s ) is zero.According to Eqs. (2a)-(2e) and Table 1, the following results can be obtained, as shown in Table 2.
The comparisons of limit pressures between this study and previous studies are listed in Table 3.It is evident that the present study provides the best result, because the limit collapse pressure calculated by this method is the largest of all methods, while the limit blow-out pressure calculated by this method is the smallest of all methods.

Conclusions
To study the face stability of circular tunnels in purely cohesive soils with linearly increasing strength, numerical simulations were performed employing the 3D FELA.A wide range of ground conditions was selected to cover a broad scope of practical application, and the rigorous UB and LB solutions of load factor (σ s − σ t )/c u0 were calculated.Using the nonlinear regression analysis method, a series of equations of the limit support pressure www.nature.com/scientificreports/σ t and the stability number N in both the collapse and blow-out cases were proposed based on 960 numerical results.The following conclusions are drawn: 1.A nonlinear relationship between (σ s − σ t )/c u0 and C/D was found in both the collapse and blow-out cases.Moreover, in the collapse case, the value of (σ s − σ t )/c u0 linearly decreases with increasing γD/c u0 and linearly increase with increasing ρD/c u0 .In the blow-out case, the value of (σ s − σ t )/c u0 linearly decreases with increasing γD/c u0 and ρD/c u0 .2. In both the collapse and blow-out cases, the failure zone in front of the tunnel face is chimney-like.In the collapse case, the failure mechanism is global and stretches out from the tunnel invert.Also, the entire surface of the face is at failure.In the blow-out case, the failure mechanism is partial in C/D < 2.0, which starts from a position above the tunnel invent instead of the tunnel invert, and the failure zone includes only the upper part of the tunnel face; when C/D > 2, the failure mechanism is also global.3. The comparisons showed that the proposed equations make an improvement over existing methods, and these equations provide a quick and accurate estimation on tunnel face stability in purely cohesive soils.4.However, our results are applicable to C/D = 0.25-5, which can already cover most practical projects.However, in an extreme larger burial depth condition, our results may be invalid, and more calculations need to be performed further.In the future research, it is suggested to further establish a more widely applicable

Figure 1 .Figure 2 .
Figure 1.Problem definition of the 3D tunnel face stability in purely cohesive soils.

Figure 4 .
Figure 4. Shear dissipation contours in limit states of tunnel face collapse and blow-out.

Figure 6 .
Figure 6.Comparisons of stability number N between this study and previous studies.

Figure 7 .Figure 8 .
Figure 7. Comparisons of the stability factors between this study and previous studies: (a) N γ , (b) N c0 and (c) N cρ .

Table 1 .
Values of coefficients a 1 -a 6 .

Table 2 .
Calculation results.www.nature.com/scientificreports/theoretical analytical model for tunnel stability analysis in purely cohesive soils based on the research results of this paper.

Table 3 .
Comparisons of limit pressures between this study and previous studies.