Analyses on face stability of shallow tunnel considering different constitutive models

Based on the finite element limit analysis method, the stability of the face in case of active failure under three constitutive models, the Mohr-Coulomb model (MC), the modified Cambridge model (MCC) and the Drucker-Prager model (DP), were analyzed. The ultimate support pressure of the face and the influence of factors such as different burial depth ratios (C/D), cohesion (c) and friction angle (φ) in the MC model are also discussed. The results show that the safety factor obtained by the MCC model under the same support pressure is always smaller than that of the MC model, and the difference is the largest when there is no support pressure. As the support pressure increases, it will gradually approach the MC model. When the support pressure is small, the safety factor obtained by the DP model is larger than the MC model, but when the support pressure is large, it is smaller than the MC model, and the final difference tends to be stable. It is necessary to select an appropriate constitutive model according to different rock masses in practical engineering. The self-stabilizing performance of the face is not affected by C/D, and the ultimate support pressure will increase with the increase of C/D, decrease linearly with the increase of cohesion, and decrease with the increase of friction angle. When the friction angle is small, the ultimate support pressure is greatly affected by C/D, and when the friction angle is large, it is hardly affected by C/D.


Introduction
The development of underground space in cities has become a trend with the rapid growth of urban rail transportation in China. The current shield construction technology is one of the leading engineering methods for urban underground tunnel construction (Sui et al., 2021) the shield excavation process, the support on the tunnel face is the guarantee to maintain the stability of the tunnel. When the support pressure is too small, the tunnel will collapse, and when the support pressure is too enormous the ground surface will bulge (Zamora Hernández et al., 2019). Therefore, how to solve the ultimate support pressure more accurately has become a sensitive research topic for most scholars.
In the early days, the stability of the tunnel face after tunnel excavation was studied by empirical formulas and simple qualitative or quantitative analysis (Broms and Bennermark, 1967). However, these methods were relatively crude in the calculation. Subsequently, the limit equilibrium method and the limit analysis method were rapidly developed for geotechnical applications (He et al., 2019a;He et al., 2019b;He et al., 2021;He et al., 2022) and were also widely used in the stability analysis of tunnel faces and could be used to find out the ultimate support pressure of tunnel faces (Chen, 1975;Davis et al., 1980;Comejo, 1989;Leca and Dormieus, 1990;Mollon et al., 2010;Han et al., 2016;Lu et al., 2017). However, these two methods have many assumptions that make it difficult to solve complex problems. The emergence of finite element methods have brought new solvers to geotechnical issues (Sloan and Assadi, 1991;Wilson et al., 2011;Zheng et al., 2012;He et al., 2019c;Shiau Al-Asadi, 2020), and with the development of computer and numerical computation technology, various numerical simulation software has become more and more widely used, and not only finite element numerical software such as ABAQUS, ANSYS, FLAC 3D and OPTUM G2 have been used for stability analysis of tunnel faces (Do et al., 2014;Li and Li, 2019;Zamora Hernández et al., 2019;Huang et al., 2020;Zheng et al., 2022), but many discrete element numerical software is also maturing and being used in the study of geotechnical properties (Bai et al., 2022).
The Mohr-Coulomb constitutive model (MC model) is mainly used in the above tunnel face stability analysis (Ji et al., 2021;Feng et al., 2022), which is popular among scholars and engineers because of its simple and practical features, and as research progresses, more and more Constitutive models are proposed and gradually used in the stability analysis of tunnel faces (Bai et al., 2019;Bai et al., 2021). The Mohr-Coulomb model considers that the compression modulus and rebound modulus of the soil are the same, and there are some cases where the Mohr-Coulomb model is not applicable, and it is necessary to select a suitable constitutive model according to different tunneling conditions. Zhao et al. (2015) used the HSS model to simulate the shield supported mechanized excavation of the Western Scheldt tunnel in the Netherlands. Keawsawasvong and Ukritchon, (2020) used the Hoek-Brown model to develop a new design equation for stability analyses of shallow unlined circular tunnels in rock masses. Liu et al. (2020) used the DP model to present a semianalytical solution as well as finite element numerical simulations for the tunnel excavation problem. Fang et al. (2022) used a number of constitutive models to account for initial soil anisotropy and non-coaxial plasticity, which were confirmed in a field investigation through the Tsinghua Park Tunnel of the Beijing-Zhang High Speed Railway in China. Therefore, this paper discusses the stability of the tunnel face under active failure calculated with three constitutive models and solves for the ultimate support pressure, and discusses the effects of parameters such as burial depth ratio C/D and cohesion c, friction angle φ. The MC model is very widely used in geotechnical analysis nowadays. It is not only a relatively simple model, but more importantly is that all of its parameters have direct physical meaning and can be measured by conventional tests such as direct shear and triaxial tests. The principal stresses σ 1 and σ 3 are used to represent the yield surface (as shown in Figure 1) as a function of In order to avoid explicit calculation of principal stresses, Nayak and Zienkiewicz (Nayak and Zienkiewicz, 1972) proposed to use the following variables based on the stress tensor principle: 1 2 s x 2 + s y 2 + s z 2 + τ xy 2 + τ xz 2 + τ yz 2 (4) Frontiers in Materials frontiersin.org 02 J 3 s x s y s z + 2τ xy τ xz τ yz + s x τ 2 yz + s y τ 2 xz + s z τ 2 xy (6) The MC model is controlled by the elastic modulus E and Poisson's ratio υ for elastic deformation, and the cohesion c, friction angle φ and shear expansion angle ψ for plastic deformation.

Modified cam-clay model (MCC)
The Cambridge model was constructed by Roscoe and Burland (Roscoe and Burland, 1968); in 1968 and was later modified and refined for the analysis of elastic-plastic deformation of clay soils (both consolidated and superconsolidated). A slightly extended version of this model was implemented according to the scheme proposed by Krabbenhoft (Krabbenhoft and Lyamin, 2012), and the yield surface (as shown in Figure 2) functioned as: Where the superscripts of the variables indicate the effective stress and p c is the prior consolidation pressure. The MCC model is controlled by the parameters of compression index λ, resilience index κ, initial pore ratio e 0 , cohesion c and friction angle φ, all of which can be obtained experimentally.   Frontiers in Materials frontiersin.org 03

Drucker-Prager model (DP)
Drucker-Prager modified the Mises criterion to propose a new yielding criterion for the damage analysis of rocks with the yield surface shown in Figure 3. The projection in the π-plane is the inner tangent circle of the MC yield surface projection shown in Figure 4, and the yield function is expressed as Where the parameters M and k are constants related to the soil friction angle φ and cohesion c respectively. To ensure that the soil material parameters match the numerical calculations and to consider the associated flow rule, the DP model and MC model parameters can be equivalently transformed by the following equation.
3 Ultimate support pressure analysis at the tunnel face In this paper, the stability of the tunnel face is analyzed by three methods: finite element limit analysis, finite difference method and limit analysis upper limit method, and the ultimate support pressure on the tunnel face is obtained.
Based on the OPTUM G2 (Finite Element Limit Analysis) and FLAC (Finite Difference Method) platforms. To simplify the analysis, the median plane of the 3D tunnel is intercepted, and a two-dimensional analysis is carried out considering the plane strain problem, as shown in Figure 5. The tunnel diameter D is 8 m. In order to ignore the influence of boundary conditions as far as possible, the length of the tunnel extension below L 1 = 8 m, the longitudinal length of the tunnel L 2 = 16 m, and the length of the model extension in the horizontal direction L 3 = 16 m. The entire bottom boundary of the model is fully constrained, and the left and right boundaries of the model and the upper and lower boundaries inside the tunnel are constrained normally to each other.
Meanwhile, this paper adopts the multi-block damage model of the tunnel face proposed by Yang Feng (Yang et al., 2010), as shown in Figure 6; the red line in the figure indicates the damaged surface, based on the theory of limit analysis upper limit method, construct the constraint and objective function and use the nonlinear solver fmincon in MATLAB to program the ultimate support pressure of the tunnel face can be found, the specific solution process and calculation formulae are shown in the literature (Nayak and Zienkiewicz, 1972).
Where the soils have a weight of 20 kN/m 3 and a cohesive force of 10 kPa, calculations were carried out using the Mohr-Coulomb yielding constitutive model. The variation of the ultimate support pressure with friction angle at the tunnel face for two working conditions (C/D = 1, C/D = 2) was explored using the three methods mentioned above, and the results are shown in Figures 7, 8.  During the analysis, the simulation of the working conditions at different C/D was controlled only by changing the thickness C of the overburden layer, while other values were kept constant. It can be seen from Figures 7, 8 that the ultimate support pressure at the tunnel face decreases as the friction angle increases, and the larger the C/D, the greater the ultimate support pressure required. All three calculation methods show the same variation pattern, and the results are close to each other, with the limit analysis upper limit method being more dangerous, FLAC (finite difference method) being more conservative, and OPTUM G2 (finite element limit analysis method) being in between and the tunnel face damage pattern (shown in Figure 9) is consistent with the damage pattern used in the limit analysis upper limit method (shown in Figure 6) hence this method is used consistently in the subsequent discussions of the constitutive model analysis.

Analysis of different constitution models 4.1 Model building
In order to analyze the effect of different constitutive models on the ultimate support pressure at the tunnel face, the three constitutive models mentioned in Section 1 were used for the analysis: i) the ideal elastoplastic constitutive Mohr-Coulomb model (MC); ii) the Modified Cambridge model (MCC), which reflects the elastic-plastic damage of soft soils well; and iii) the Drucker-Prager model (DP), which reflects the elastic-plastic damage of rocks. The calculation method was selected from the finite element limit analysis method used in the previous section.
Refer to Figure 5 for a geometric model of the tunnel with geometric parameters taken from the literature , where the diameter of the tunnel D is 9 m, the extension length below L 1 = 9 m, the longitudinal length L 2 = 18 m, and the horizontal extension length L 3 = 18 m.

Soil parameters
The soil parameters in the tunnel model were taken from the literature  and the values of the soil When C/D = 1 in MC model, the ultimate support pressure varies with friction.

FIGURE 8
When C/D = 2 in MC model, the ultimate support pressure varies with friction.

Analysis of results
Based on the finite element limit analysis and strength reduction method, this paper analyses the stability of the tunnel face under three different constitutive models and discusses the effect of varying burial depth ratios C/D on the ultimate support pressure.

Analysis of the results of different constitution models
Considering two working conditions with C/D is 1 and 1.5, calculate the variation of safety factor with support pressure at the tunnel face under three constitutive models (see Figure 10).
It can be seen from Figure 10 that the factor of safety increases with the increase in support pressure regardless of the working conditions or the constitutive model (only active damage is considered). For the same support pressure, the safety factor obtained with the MC model is always slightly higher than that of the MCC model; when the safety factor is less than 1, the safety factor obtained with the DP model is greater than that of the MC and MCC models for the same support pressure, but when the required safety factor is greater than 1, the safety factor calculated with the DP model is less than that of the MC model (this is consistent with the fact that the yield surface of the DP model is tangent to the yield surface of the MC model). The ultimate support pressures (safety factor equal to 1) obtained with the MC and DP models are almost identical, while the ultimate support pressures obtained with the MCC model are relatively large.
To reflect the influence of different constitutive models on the stability of the tunnel face more intuitively, the percentage difference in the safety factor under different constitutive models (compared with the MC model) is discussed, and the burial depth ratio C/D is taken as 1 for calculation, and the results are shown in Figure 11.
The positive sign in Figure 11 indicates that the result is larger than that calculated by the MC model, and the negative sign indicates that it is smaller. It can be seen that the safety

FIGURE 10
The safety factor varies with support pressure under different constitutive conditions. (A) C/D = 1, (B) C/D = 1.5.
Frontiers in Materials frontiersin.org 06 coefficient calculated by the MCC model is smaller than that of the MC model under the same support pressure, and the difference is the largest when the support pressure is not considered, which is 7.8% smaller than that of the MC. The coefficient of safety at the intersection of the DP model and the MC model is 1, indicating that the ultimate support pressure under the MC model and the DP model is the same. When the support pressure is small, the coefficient of safety obtained under the DP model is larger than that under the MC model, and the difference is the largest when the support pressure is not considered, with a value of 17.5%. The difference tends to stabilize and gradually converges to 6.6% in this paper.

Analysis of C/D results for different burial depth ratios
In the discussion of the previous section, it was found that the burial depth ratio C/D has a greater effect on the ultimate support pressure, so the effect of the C/D ratio on the stability of the tunnel face and the ultimate support pressure under active failure was further investigated, and the results are shown in Figure 12.
When the support pressure is zero, the coefficient of safety obtained from different C/D is almost the same, indicating that the self-stabilizing performance of the tunnel face is almost independent of C/D. As the support pressure increases, the smaller the C/D the greater the coefficient of safety, and the greater the coefficient of safety is influenced by the support pressure (the slope of the curve in the graph). And the greater the C/D, the greater the ultimate support pressure obtained and the more dangerous the tunnel face is.

FIGURE 11
Percentage of safety factor gap when C/D = 1. The MC model was selected in order to analyze the influence of the strength parameters of the soil on the ultimate support pressure at the tunnel face further. The tunnel geometry model constructed in Section 3.1 was used to explore the influence of the friction angle c and cohesion φ on the ultimate support pressure at the tunnel face in the MC model based on the finite element limit analysis method, and the influence of four different C/D working conditions was also analyzed, and the calculation results are shown in Figure 13.
It can be seen that the ultimate support pressure at the tunnel face decreases with the increase of cohesion, and presents a linear relationship, the larger the C/D, the greater the ultimate support pressure, and the ultimate support pressure-viscosity curve under different C/D conditions is approximately parallel; the ultimate support pressure also decreases with the increase of friction angle, presenting a non-linear relationship, and finally tends to a stable value. When the friction angle is small, it is influenced by C/D, however, when the friction angle is large, it is almost independent of C/D, and the ultimate support pressure under different C/D is close to the same.

Conclusion
This paper adopts the finite element limit analysis method to establish a two-dimensional tunnel face model, analyses the stability of the tunnel face under three different constitutive models, and solves for the ultimate support pressure, taking into account the effects of parameters such as burial depth ratio, cohesion and friction angle in the MC model. And leading to the following conclusions: (1) Compared to the MC model, the safety factor obtained by the MCC model for the same support pressure is always smaller than that of the MC model, with a maximum difference of 7.8% when there is no support pressure, but gradually converges to that of the MC model as the support pressure increases. This shows that the MCC model is more conservative than the MC model. (2) Compared to the MC model, the DP model gives a greater safety factor than the MC model when the support pressure is small and 17.5% greater than the MC model when there is no support pressure; however, when the support pressure is large, the safety factor is less than the MC model and the final difference leveled off at 6.6%. (3) The larger the C/D, the greater the required ultimate support pressure of the tunnel face. When the support pressure is not considered, the safety coefficients under different C/D calculations are the same, which means that the self-stabilizing performance of the tunnel face is almost independent of C/D. When the support pressure is greater, the smaller the C/D is the greater the influence. (4) The ultimate support pressure at the tunnel face decreases linearly with the increase of cohesion, and the curves of the ultimate support pressure with cohesion are approximately parallel to each other at different burial depth ratios C/D. The larger the friction angle is, the higher the friction angle is, the ultimate support pressure at the tunnel face decreases non-linearly with the friction angle and finally tends to a stable value. When the friction angle is small, the ultimate support pressure is more influenced by C/D, and when the friction angle is large, it is less influenced by C/D. At this time, the ultimate support pressure under different C/D is almost the same. Frontiers in Materials frontiersin.org 08

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/supplementary material.

Conflicts of interest
Authors LC was employed by the company Road Tunnel Branch of Sichuan Highway Bridge Construction Group Co., Ltd.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher's note All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.