A New Numerical Finite Strain Procedure for a Circular Tunnel Excavated in Strain-Softening Rock Masses and Its Engineering Application

Abstract

The small strain theory underestimates the self-bearing capacity of rock masses, especially for a soft rock tunnel under high geostress. To perform an efficient and accurate calculation and provide a reference for the stiffness design of a tunnel, the finite strain solution for a circular tunnel in Mohr–Coulomb strain-softening rock masses with a non-associated flow rule was derived as three sets of differential equations under the Lagrangian coordinate, which are in the residue region, the softening region, and the elastic region, respectively. Based on the bisection method, an iteration procedure for solving the finite strain solution was proposed to approximate the boundary condition at infinity, the values of two adjacent boundaries, and the initial values on the excavation boundary. This numerical procedure was verified by comparing with self-similar solutions, recursive solutions, and FLAC simulation results. In the calculation example, the relative error on boundaries can be decreased to less than 10⁻⁸ after only 10 times iteration and the time for each calculation is less than 15 s. Applying this procedure on the sensibility analysis and stiffness reliability design for the Zhongyi tunnel, a support stiffness of 4.3 MPa/m is recommended to guarantee a tunnel displacement lower than 0.5 m.

1. Introduction

With the rapid construction of transportation infrastructure, more and more tunnels have been built in soft rock under high geostress. Large deformation occurs frequently in such tunnel engineering and even causes roof collapse, steel ribs’ twisting and breaking, shotcrete cracking, lining cracking, the invasion of primary lining into the space of secondary lining, and other disasters [1]. For instance, the maximum deformation of the Muzhailing Tunnel located upon the Lanzhou–Haikou National Expressway exceeds 2000 mm, far exceeding the deformation allowance [2]. Similar problems are encountered in the Muzhailing Railway Tunnel [3], Zhegu Mountain Tunnel [4], Huangjiazhai Tunnel [5], Minxian Tunnel [6] and have also been predicted to appear during the construction of 43 soft rock tunnels with high geostress in the Sichuan–Tibet railway [7]. The large deformation hazard raises significant challenges for the design of the tunnel support system.

In recent years, scholars have been conducting considerable research on the mechanical calculation of tunnels excavated in strain-softening rock masses. Alonso et al. [8] and Park [9] derived a similarity solution for the circle tunnel in strain-softening rock mass, which it still needed to be solved by a numerical method. Guan et al. [10] compared a simplified theoretical method with a rigorous theoretical method for tunnel excavation in the rock masses exhibiting strain-softening behavior. Park et al. [11], Lee et al. [12], Wang et al. [13], Zhang et al. [14], Han et al. [15], Cui et al. [16], Zou et al. [17], Wang et al. [18], Zhang et al. [19], Xia et al. [20], and GhorbaniHadi et al. [21] proposed a series of new semi-analytical solutions and numerical procedures for calculating the stress and displacement around a circular tunnel excavated in strain-softening Mohr–Coulomb, generalized Hoek–Brown, or three-dimensional nonlinear Hoek–Brown rock mass. Although these methods were applied well on the mechanical analysis and the tunnel design [16,22,23,24,25,26,27], they are based on the small deformation hypothesis and not applicable to analyze a tunnel with a high risk of large deformation.

For the large strain problem, Durban [28], Yu and Rowe [29], and Vrakas [30,31] derived finite strain closed-form solutions of circular opening issues in the elastic and the elasto-perfectly-plastic continuous medium based on a geometric equation in Lagrangian coordinates [32]. Park [33] proposed a large strain similarity solution, expressed as first-order ordinary differential equations, for a spherical or circular opening in elasto-perfectly-plastic rock masses. Mo et al. [34] described the confinement-convergence responses for deep tunnels in linearly elastic, non-associated MC, and brittle Hoek–Brown rock masses using large-strain solutions of unified spherical and cylindrical cavity contraction. Zhang et al. [35] derived a large strain differential equations for elasto-brittle-plastic rock masses and implemented a numerical procedure.

For a tunnel in strain-softening rock masses, Guan et al. [36] analyzed a highly deformed circular tunnel in the Lagrangian coordinate by dividing the plastic zone into $n$ concentric annuli and using analytical solutions for the elastic-brittle-plastic rock mass within each annulus. Zhang et al. [37] derived the differential equations in the original coordinate by the stress equilibrium equation and deformation compatibility equation and proposed one numerical procedure which updated the variables of (i + 1)th annulus in the softening region by the fifth-order Runge–Kutta method according to the data of the ith annulus. Then Zhang et al. [38] converted differential equations into integral equations and implemented a recursive procedure named ‘softening annulus loop’. In a similar method, Xu et al. [39] developed a numerical solution for strain-softening rock masses obeying the GZZ strength criterion. The recursion method in [36,37,38,39] inevitably leads to an accumulative error which cannot be evaluated quantitatively and adjusted automatically. Increasing the number of concentric annuli can reduce the accumulative error only by an unknown extent and results in a loss of computation efficiency.

This paper aims at an accurate and efficient procedure for numerical finite strain solutions of circular tunnels excavated in strain-softening rock masses to provide more practical and applicable design references for soft rock tunnels with high geostress. In the Lagrangian coordinate, the solutions in the elastic, softening and residue region of a circular tunnel are derived in the form of the differential equations. A new procedure for a numerical finite strain solution is implemented based on the iteration and bisection method. The procedure’s accuracy and efficiency are verified by comparing it with previous studies. This study performs a global sensitivity analysis of tunnel displacement using this procedure and recommends the supporting stiffness for the design of the Zhongyi tunnel based on the stochastic analysis results.

2. Problem Statement

2.1. Model Description

The calculation model is shown in Figure 1. It depicts an axisymmetric plane strain issue of a circular tunnel excavated in an infinite, continuous, homogeneous, and isotropic rock mass, subjected to uniform hydrostatic stress $P_0$ at infinity and uniform support pressure ${\sigma }_0$ exerted on the wall. The original radius $R$ and the deformed radius $r$ represent the initial position and deformed(current) position of a material point, respectively. Consequently, the radial displacement of the material point $u$ can be defined as:

(1)$u=R-r$

After the excavation, the tunnel radius decreases from its initial value $R_0$ to the deformed radius $r_0$. For a rock mass following elastic-strain-softening behaviour, the rock mass within a deformed radius $r_{\mathrm{r}}$ reaches the residual state, the region outside a deformed radius $r_{\mathrm{p}}$ stays in the elastic state and the rock mass between $r_{\mathrm{s}}$ and $r_{\mathrm{p}}$ is in the strain-softening state.

2.2. Governing Equations

For the finite strain model, the stress equilibrium equation of the rock masses in the deformed state can be expressed as the differential equation, Equation (2).

(2)$\frac{\mathrm{d}{\sigma }_{\mathrm{r}}}{\mathrm{d}r}+\frac{{\sigma }_{\mathrm{r}}-{\sigma }_{\mathsf{\theta }}}{r}=0$

where ${\sigma }_{\mathrm{r}}$ and ${\sigma }_{\mathsf{\theta }}$ denote the radial stress and tangential stress, respectively.

The logarithmic function can accurately describe the geometric relationship between radial displacement and strains. In the current polar coordinates, the geometric equation can be expressed as follows:

(3)${\epsilon }_{\mathsf{\theta }}=\ln \left(\frac{u}{r}+1\right)$

(4)${\epsilon }_{\mathrm{r}}=\ln \left(\frac{du}{dr}+1\right)$

where ${\epsilon }_{\mathsf{\theta }}$ and ${\epsilon }_{\mathrm{r}}$ denote the radial and tangential strains, respectively.

The total strain induced by excavation can be decomposed into the elastic part and the plastic part as follows:

(5)${\epsilon }_{\mathrm{r}}=\frac{1}{2G}\left[\left(1-v\right)\left({\sigma }_{\mathrm{r}}-P_0\right)-v\left({\sigma }_{\mathsf{\theta }}-P_0\right)\right]+{\epsilon }_{\mathrm{r}}^{\mathrm{p}}$

(6)${\epsilon }_{\mathsf{\theta }}=\frac{1}{2G}\left[\left(1-v\right)\left({\sigma }_{\mathsf{\theta }}-P_0\right)-v\left({\sigma }_{\mathrm{r}}-P_0\right)\right]+{\epsilon }_{\mathsf{\theta }}^{\mathrm{p}}$

where the shear modulus $G=E/2/\left(1+v\right)$, $E$ and $v$ are Young’s modulus and Poisson’s ratio, respectively; ${\epsilon }_{\mathrm{r}}^{\mathrm{p}}=0$ and ${\epsilon }_{\mathsf{\theta }}^{\mathrm{p}}=0$ in the elastic rock mass.

2.3. Softening Parameters, Yield Criterion, and Evolutional Law

The plastic shear strain ${\gamma }^{\mathrm{p}}$ is the most widely used softening index governing material parameters, which is defined as Equation (7).

(7)${\gamma }^{\mathrm{p}}={\epsilon }_{\mathsf{\theta }}^{\mathrm{p}}-{\epsilon }_{\mathrm{r}}^{\mathrm{p}}$

According to the non-associated MC potential function, the relation of radial and tangential plastic strains can be expressed as Equation (8).

(8)${\epsilon }_{\mathrm{r}}^{\mathrm{p}}+K\left({\gamma }^{\mathrm{p}}\right){\epsilon }_{\mathsf{\theta }}^{\mathrm{p}}=0$

where $K\left({\gamma }^{\mathrm{p}}\right)=\frac{1+\sin \psi \left({\gamma }^{\mathrm{p}}\right)}{1-\sin \psi \left({\gamma }^{\mathrm{p}}\right)}$ and $\psi$ is the dilation angle.

The failure of the rock mass is governed by the strain-softening Mohr–Coulomb yielding a function written as:

(9)${\sigma }_{\mathsf{\theta }}=k\left({\gamma }^{\mathrm{p}}\right){\sigma }_{\mathrm{r}}+{\sigma }_{\mathrm{c}}\left({\gamma }^{\mathrm{p}}\right)$

where $k\left({\gamma }^{\mathrm{p}}\right)=\frac{1+\sin \phi \left({\gamma }^{\mathrm{p}}\right)}{1-\sin \phi \left({\gamma }^{\mathrm{p}}\right)}$, ${\sigma }_{\mathrm{c}}\left({\gamma }^{\mathrm{p}}\right)=\frac{2c\left({\gamma }^{\mathrm{p}}\right)\cos \phi \left({\gamma }^{\mathrm{p}}\right)}{1-\sin \phi \left({\gamma }^{\mathrm{p}}\right)}$; $c$ and $\phi$ are the cohesive strength and the friction angle respectively.

The evolution of material property parameters can be described by the piecewise linear function [38]. Any strength parameter of the rock mass is linearly decreased from the peak value ${\eta }_{\mathrm{p}}$ to the corresponding residual value ${\eta }_{\mathrm{r}}$ as shown in Equation (10).

(10)$\eta \left({\gamma }^{\mathrm{p}}\right)=\{\begin{array}{cc}{\eta }_{\mathrm{p}}& {\gamma }^{\mathrm{p}}\le 0\\ {\eta }_{\mathrm{p}}-\left({\eta }_{\mathrm{p}}-{\eta }_{\mathrm{r}}\right){\gamma }^{\mathrm{p}}/{\gamma }^{\mathrm{p}\ast }& 0<{\gamma }^{\mathrm{p}}<{\gamma }^{\mathrm{p}\ast }\\ {\eta }_{\mathrm{r}}& {\gamma }^{\mathrm{p}}\ge {\gamma }^{\mathrm{p}\ast }\end{array}$

where $\eta$ denotes any one of the strength parameters such as $\psi$, $c,$ and $\phi$ etc.; ${\gamma }^{\mathrm{p}\ast }$ is the critical plastic shear strain, the subscripts ‘p’ and ‘r’ of $\eta$ denote the peak and residual values respectively.

2.4. Boundary Conditions

In the original or deformed state, the excavation disturbance for the rock mass at infinity is infinitely small. Hence, the stress boundary conditions of the rock mass can be expressed in the original coordinate as follows:

(11)$\{\begin{array}{cc}{\sigma }_{\mathrm{r}}={\sigma }_{\mathsf{\theta }}=P_0& R\to \infty \\ {\sigma }_{\mathrm{r}}=p& R=R_0\end{array}$

3. Solutions for Stress and Displacement of Rock Mass

Using Equations (1), (3), and (4), the differential relation between the original coordinate and the deformed is expressed as:

(12)$\frac{\mathrm{d}r}{r}=\frac{\mathrm{d}R}{R}{e}^{{\epsilon }_{\mathsf{\theta }}-{\epsilon }_{\mathrm{r}}}$

Deriving from the above equations, the stress and displacement solutions are expressed in the form of differential equations of the radial stress ${\sigma }_{\mathrm{r}}$, the plastic shear strain ${\gamma }^{\mathrm{p}}$, and the tangential stress ${\sigma }_{\mathsf{\theta }}$.

3.1. Solution in the Residue Region

In the residue region, Equations (7)–(9) are simplified as Equations (13)–(15).

(13)${\epsilon }_{\mathsf{\theta }}^{\mathrm{p}}=\frac{\left(1-\sin {\psi }_{\mathrm{r}}\right){\gamma }^{\mathrm{p}}}{2}$

(14)${\epsilon }_{\mathrm{r}}^{\mathrm{p}}+\frac{1+\sin {\psi }_{\mathrm{r}}}{1-\sin {\psi }_{\mathrm{r}}}{\epsilon }_{\mathsf{\theta }}^{\mathrm{p}}=0$

(15)${\sigma }_{\mathsf{\theta }}=\frac{1-\sin {\phi }_{\mathrm{r}}}{1+\sin {\phi }_{\mathrm{r}}}{\sigma }_{\mathrm{r}}+\frac{2c_{\mathrm{r}}\cos {\phi }_{\mathrm{r}}}{1-\sin {\phi }_{\mathrm{r}}}$

Using Equations (5), (6), and (9), the differential formula $\mathrm{d}r/r$ can be rewritten as:

(16)$\frac{\mathrm{d}r}{r}=\frac{\mathrm{d}R}{R}\exp \left[\frac{1}{2G}\left(\frac{-2\sin {\phi }_{\mathrm{r}}}{1+\sin {\phi }_{\mathrm{r}}}{\sigma }_{\mathrm{r}}+\frac{2c_{\mathrm{r}}\cos {\phi }_{\mathrm{r}}}{1-\sin {\phi }_{\mathrm{r}}}\right)+{\gamma }^{\mathrm{p}}\right]$

By substituting radial stress and $\mathrm{d}r/r$ in Equation (16) into Equation (2), the first differential equation of the stress equilibrium equation can be obtained, expressed as Equation (17).

(17)$\frac{\mathrm{d}{\sigma }_{\mathrm{r}}}{\mathrm{d}R}=\frac{1}{R}\left(\frac{-2\sin {\phi }_{\mathrm{r}}}{1+\sin {\phi }_{\mathrm{r}}}{\sigma }_{\mathrm{r}}+\frac{2c_{\mathrm{r}}\cos {\phi }_{\mathrm{r}}}{1-\sin {\phi }_{\mathrm{r}}}\right)\exp \left[\frac{1}{2G}\left(\frac{-2\sin {\phi }_{\mathrm{r}}}{1+\sin {\phi }_{\mathrm{r}}}{\sigma }_{\mathrm{r}}+\frac{2c_{\mathrm{r}}\cos {\phi }_{\mathrm{r}}}{1-\sin {\phi }_{\mathrm{r}}}\right)+{\gamma }^{\mathrm{p}}\right]$

The differential form of the tangential strain can be derived as Equation (18) by differentiating the tangential strain regarding the initial radius in Equation (3) and substituting it into Equation (16).

(18)$\frac{\mathrm{d}{\epsilon }_{\mathsf{\theta }}}{\mathrm{d}R}=\frac{1}{R}-\frac{\mathrm{d}r}{r\mathrm{d}R}=\frac{1}{R}\left\{1-\exp \left[\frac{1}{2G}\left(\frac{-2\sin {\phi }_{\mathrm{r}}}{1+\sin {\phi }_{\mathrm{r}}}{\sigma }_r+\frac{2c_{\mathrm{r}}\cos {\phi }_{\mathrm{r}}}{1-\sin {\phi }_{\mathrm{r}}}\right)+{\gamma }^{\mathrm{p}}\right]\right\}$

Seeking the derivative of each term in Equation (6) regarding the initial radius $r$, and another partial differential of the tangential strain can be expressed as Equation (19).

(19)$\frac{\mathrm{d}{\epsilon }_{\mathsf{\theta }}}{\mathrm{d}R}=\frac{1}{2G}\left[\left(1-v\right)\frac{1-\sin {\phi }_{\mathrm{r}}}{1+\sin {\phi }_{\mathrm{r}}}-v\right]\frac{\mathrm{d}{\sigma }_{\mathrm{r}}}{\mathrm{d}R}+\frac{1-\sin {\psi }_{\mathrm{r}}}{2}\frac{\mathrm{d}{\gamma }^{\mathrm{p}}}{\mathrm{d}R}$

From Equations (18) and (19), the second differential equation for the plastic shear strain regarding the initial radius is as follows:

(20)$\frac{\mathrm{d}{\gamma }^{\mathrm{p}}}{\mathrm{d}R}=\frac{k_1}{G}\left(v-\left(1-v\right)\frac{1}{k}\right)\frac{\mathrm{d}{\sigma }_{\mathrm{r}}}{\mathrm{d}R}+\frac{2k_1}{R}\left\{1-\exp \left[\frac{1}{2G}\left(\frac{-2\sin {\phi }_{\mathrm{r}}}{1+\sin {\phi }_{\mathrm{r}}}{\sigma }_{\mathrm{r}}+{\sigma }_{\mathrm{c}}\right)+{\gamma }^{\mathrm{p}}\right]\right\}$

where $k_1=1/\left(1-\sin {\psi }_{\mathrm{r}}\right)$.

The tangential stress can be easily obtained by substituting the radial stress into Equation (7) and the displacement of rock mass material point can be calculated using Equations (1), (3), (6), and (13).

(21)$u=R\left\{1-\exp \left[\frac{-\left(1-v\right)\left({\sigma }_{\mathrm{r}}-P_0\right)+v\left({\sigma }_{\mathsf{\theta }}-P_0\right)}{2G}-\frac{\left(1-\sin {\psi }_{\mathrm{r}}\right){\gamma }^{\mathrm{p}}}{2}\right]\right\}$

3.2. Solution in the Softening Region

The derivation method of two differential equations for the radial stress and the plastic shear strain in the softening region is like that in the residue region. Because the rock mass strength properties vary with the plastic shear strain, the differential equations become more complex.

Using Equations (5)–(7) and (9), the differential formula $\mathrm{d}r/r$ is rewritten as follows:

(22)$\frac{\mathrm{d}r}{r}=\frac{\mathrm{d}R}{R}\exp \left\{\frac{1}{2G}\left[\left(k\left({\gamma }^{\mathrm{p}}\right)-1\right){\sigma }_{\mathrm{r}}+{\sigma }_{\mathrm{c}}\left({\gamma }^{\mathrm{p}}\right)\right]+{\gamma }^{\mathrm{p}}\right\}$

Substituting the tangential stress in Equation (7) and $\mathrm{d}r/r$ in Equation (22) into the stress equilibrium equation, the first differential equation of the radial stress regarding the original radius is obtained as follows:

(23)$\frac{\mathrm{d}{\sigma }_{\mathrm{r}}}{\mathrm{d}R}=\frac{\left(k\left({\gamma }^{\mathrm{p}}\right)-1\right){\sigma }_{\mathrm{r}}+{\sigma }_{\mathrm{c}}\left({\gamma }^{\mathrm{p}}\right)}{R}\exp \left\{\frac{1}{2G}\left[\left(k_{\mathrm{p}}\left({\gamma }^{\mathrm{p}}\right)-1\right){\sigma }_{\mathrm{r}}+{\sigma }_{\mathrm{c}}\left({\gamma }^{\mathrm{p}}\right)\right]+{\gamma }^{\mathrm{p}}\right\}$

Like Equations (18) and (19), the differential form of the tangential strain to the original radius can be written as Equations (24) and (25) respectively.

(24)$\frac{\mathrm{d}{\epsilon }_{\mathsf{\theta }}}{\mathrm{d}R}=\frac{1}{R}-\frac{1}{R}\exp \left[\frac{\left(k_{\mathrm{p}}\left({\gamma }^{\mathrm{p}}\right)-1\right){\sigma }_{\mathrm{r}}+{\sigma }_{\mathrm{c}}\left({\gamma }^{\mathrm{p}}\right)}{2G}+{\gamma }^{\mathrm{p}}\right]$

(25)$\frac{\mathrm{d}{\epsilon }_{\mathsf{\theta }}}{\mathrm{d}R}=\frac{1}{2G}\left[\left(1-v\right)\frac{\mathrm{d}{\sigma }_{\mathsf{\theta }}}{\mathrm{d}R}-v\frac{\mathrm{d}{\sigma }_{\mathrm{r}}}{\mathrm{d}R}\right]+\frac{\mathrm{d}{\epsilon }_{\mathsf{\theta }}^{\mathrm{p}}}{\mathrm{d}R}$

where $\frac{\mathrm{d}{\sigma }_{\mathsf{\theta }}}{\mathrm{d}R}=\left[{\sigma }_{\mathrm{r}}\frac{\partial k\left({\gamma }^{\mathrm{p}}\right)}{\partial {\gamma }^{\mathrm{p}}}+\frac{\partial {\sigma }_{\mathrm{c}}\left({\gamma }^{\mathrm{p}}\right)}{\partial {\gamma }^{\mathrm{p}}}\right]\frac{\mathrm{d}{\gamma }^{\mathrm{p}}}{\mathrm{d}R}+k\left({\gamma }^{\mathrm{p}}\right)\frac{\mathrm{d}{\sigma }_{\mathrm{r}}}{\mathrm{d}R}$ and $\frac{\mathrm{d}{\epsilon }_{\mathsf{\theta }}^{\mathrm{p}}}{\mathrm{d}R}=\frac{\partial \left[\frac{{\gamma }^{\mathrm{p}}}{1+K\left({\gamma }^{\mathrm{p}}\right)}\right]}{\partial {\gamma }^{\mathrm{p}}}\frac{\mathrm{d}{\gamma }^{\mathrm{p}}}{\mathrm{d}R}$.

Hence, the differential equation for the plastic shear strain in the softening region is derived as Equation (26).

(26)$\frac{\mathrm{d}{\gamma }^{\mathrm{p}}}{\mathrm{d}R}=\frac{\frac{1}{R}-\frac{1}{R}exp\left[\frac{\left(k\left({\gamma }^{\mathrm{p}}\right)-1\right){\sigma }_{\mathrm{r}}+{\sigma }_{\mathrm{c}}\left({\gamma }^{\mathrm{p}}\right)}{2G}+{\gamma }^{\mathrm{p}}\right]-\frac{\left(1-v\right)k\left({\gamma }^{\mathrm{p}}\right)-v}{2G}\frac{\mathrm{d}{\sigma }_{\mathrm{r}}}{\mathrm{d}R}}{\frac{1-v}{2G}\left[{\sigma }_{\mathrm{r}}\frac{\partial k\left({\gamma }^{\mathrm{p}}\right)}{\partial {\gamma }^{\mathrm{p}}}+\frac{\partial {\sigma }_{\mathrm{c}}\left({\gamma }^{\mathrm{p}}\right)}{\partial {\gamma }^{\mathrm{p}}}\right]+\frac{\partial \left[\frac{{\gamma }^{\mathrm{p}}}{1+K\left({\gamma }^{\mathrm{p}}\right)}\right]}{\partial {\gamma }^{\mathrm{p}}}}$

After solving Equations (23) and (26), the tangential stress can be obtained easily by substituting the radial stress into Equation (7) and the displacement of the rock mass material point in the softening region can be calculated by Equation (27).

(27)$u=R\left\{1-\exp \left[\frac{\left(vk\left({\gamma }^p\right)-\left(1-v\right)\right){\sigma }_{\mathrm{r}}+\left(1-2v\right)P_0+v{\sigma }_{\mathrm{c}}\left({\gamma }^{\mathrm{p}}\right)}{2G}-\frac{\left(1-\sin \psi \left({\gamma }^{\mathrm{p}}\right)\right){\gamma }^{\mathrm{p}}}{2}\right]\right\}$

3.3. Solution in the Elastic Region

In the elastic region, there is no plasticity. Hence, the radial stress and tangential stress are taken as basic variables. The stress equilibrium equation and the compatibility equation can be easily rewritten into differential equations of the radial stress and the tangential stress, as Equations (28) and (29).

(28)$\frac{\mathrm{d}{\sigma }_{\mathrm{r}}}{\mathrm{d}R}+\frac{{\sigma }_{\mathrm{r}}-{\sigma }_{\mathsf{\theta }}}{R}\exp \left(\frac{{\sigma }_{\mathsf{\theta }}-{\sigma }_{\mathrm{r}}}{2G}\right)=0$

(29)$\frac{\mathrm{d}{\sigma }_{\mathsf{\theta }}}{\mathrm{d}R}=\frac{1}{1-v}\left\{\frac{2G}{R}\left[1-\exp \left(\frac{{\sigma }_{\mathsf{\theta }}-{\sigma }_{\mathrm{r}}}{2G}\right)\right]+v\frac{{\sigma }_{\mathsf{\theta }}-{\sigma }_{\mathrm{r}}}{R}\exp \left(\frac{{\sigma }_{\mathsf{\theta }}-{\sigma }_{\mathrm{r}}}{2G}\right)\right\}$

The displacement of the rock mass material point in the elastic region is calculated using Equation (30).

(30)$u=R\left\{1-\exp \left[\frac{-\left(1-v\right)\left({\sigma }_{\mathrm{r}}-P_0\right)+v\left({\sigma }_{\mathsf{\theta }}-P_0\right)}{2G}\right]\right\}$

4. Numerical Implementation Procedure

According to Equations (28)–(30), the stress and displacement of the rock mass in the elastic region are determined by the initial values of the radial stress and tangential stress on the boundary between the elastic region and the plastic region ($R_{\mathrm{p}}$). In the softening region, the tangential stress can be calculated by the radial stress and the plastic shear strain by the yielding function. Accordingly, the radial stress and tangential stress on $R_{\mathrm{p}}$ can be obtained by solving Equations (23) and (26), once the initial values of the radial stress and the plastic shear strain on the boundary between the softening region and the residue region ($R_{\mathrm{r}}$) are known. Similarly, the initial values of the radial stress and the plastic shear strain on the tunnel inner boundary ($R_0$) determine the stress, the plastic shear strain, and the displacement of any point in the residue region according to Equations (17), (20), and (21). Hence, there are strict one-to-one correspondences between the rock mass stress $P_0$ at infinity and the plastic shear strain which decrease monotonously in the residue and softening regions with $R$ increasing when the support pressure is given.

4.1. Numerical Procedures for Calculating Limit Support Pressure

When the support pressure equals the softening limit value $p_{\mathrm{rp}}^{\ast }$, the plastic shear strain on the tunnel inner boundary $R_0$ exactly equals the critical plastic shear strain ${\gamma }^{\mathrm{p}\ast }$. Comparatively, the plastic shear strain is equivalent to zero and the tangential stress and radial stress on the tunnel inner boundary $R_0$ just meet the yield criteria when the support pressure equals the elastic limit value $p_{\mathrm{pe}}^{\ast }$. There are three states of the surrounding rock, the residue-softening-elastic state, the softening-elastic state, and the elastic state, with the support pressure being less than $p_{\mathrm{rp}}^{\ast }$, between $p_{\mathrm{rp}}^{\ast }$ and $p_{\mathrm{pe}}^{\ast }$ and greater than $p_{\mathrm{pe}}^{\ast }$ respectively. The flowcharts depicting the calculation procedures of $p_{\mathrm{rp}}^{\ast }$ and $p_{\mathrm{pe}}^{\ast }$ are shown in Figure 2 and Figure 3, respectively.

The error tolerances of $\alpha$ and ${\gamma }_{\mathrm{p}}^{R_{\mathrm{p}}}$ and the value of $k_{\infty }$ have an important influence on the precision and accuracy of the calculation results. Ordinarily, greater $k_{\infty }$ and smaller $\mathrm{Tol}\left({\gamma }_{\mathrm{p}}^{R_{\mathrm{p}}}\right)$ and $\mathrm{Tol}\left(\alpha \right)$ lead to better precision and higher accuracy of the results. In this paper, $\mathrm{Tol}\left({\gamma }_{\mathrm{p}}^{R_{\mathrm{p}}}\right)$, $\mathrm{Tol}\left(\alpha \right)$, and $k_{\infty }$ are taken as to ${10}^{-11}$, ${10}^{-8}$, and ${10}^5$ respectively, of which the rationality is verified in Section 5.

In Figure 2, the first values of $R_{\mathrm{p}1}$ and $R_{\mathrm{p}2}$ in the loop can be determined by the following method: (1) set $R_{\mathrm{p}}\left(k\right)=R_0+{\delta }_1+\left(k-1\right){\delta }_2$ and calculate the plastic shear strain ${\gamma }_{\mathrm{p}}^{R_{\mathrm{p}}}\left(k\right)$ when $R=R_{\mathrm{p}}\left(k\right)$ by solving Equations (17) and (20); (2) repeat the above calculation from $k=1$ until ${\gamma }_{\mathrm{p}}^{R_{\mathrm{p}}}\left(k\right)>0$ and ${\gamma }_{\mathrm{p}}^{R_{\mathrm{p}}}\left(k+1\right)<0$ and take $R_{\mathrm{p}}\left(k\right)$ and $R_{\mathrm{p}}\left(k+1\right)$ as $R_{\mathrm{p}1}$ and $R_{\mathrm{p}2}$ respectively. The value of ${\delta }_1$ and ${\delta }_2$ can be adjusted according to the accuracy and computational efficiency. For the stability of solutions, take ${\delta }_1$ and ${\delta }_2$ equal to ${10}^{-5}{R}_0$ and ${10}^{-4}{R}_0$, respectively, in this paper.

4.2. Numerical Procedures for Three States

On the condition that the support pressure and other tunnel parameters are given, the stress and displacement calculation can be divided into three scenarios when the support pressure $p$ is between $p_{\mathrm{pe}}^{\ast }$ and $P_0$, between $p_{\mathrm{pe}}^{\ast }$ and $P_0$, and smaller than $p_{\mathrm{pe}}^{\ast }$. The following sections elaborate the numerical procedures for the calculation of stress and displacement in the above three states.

4.2.1. Numerical Procedure for Elastic State

When the support pressure is greater than $p_{\mathrm{pe}}^{\ast }$ and smaller than the geostress $P_0$, the rock mass stays elastic and the stress and displacement of the rock masses can be obtained by the procedure coded as the flowchart in Figure 4.

In Figure 4, the selecting method of the first values of ${\sigma }_{\mathsf{\theta }1}^{R_0}$ and ${\sigma }_{\mathsf{\theta }2}^{R_0}$ is like that of $R_{\mathrm{p}1}$ and $R_{\mathrm{p}2}$ in Section 3. At first, calculate ${\sigma }_r^{R_{\infty }}$ when $R=k_{\infty }{R}_0$.and ${\sigma }_{\mathsf{\theta }}^{R_0}\left(k\right)=0.1\left(k-1\right)P_0$ by solving Equations (28) and (29) and repeat this calculation from $k=1$ until ${\sigma }_{\mathrm{r}}^{R_{\infty }}\left(k\right)>P_0$ and ${\sigma }_{\mathrm{r}}^{R_{\infty }}\left(k+1\right)<P_0$. Then take ${\sigma }_{\mathrm{r}}^{R_{\infty }}\left(k\right)$ and ${\sigma }_{\mathrm{r}}^{R_{\infty }}\left(k+1\right)$ as the initial values of ${\sigma }_{\mathsf{\theta }1}^{R_0}$ and ${\sigma }_{\mathsf{\theta }2}^{R_0}$ in the loop respectively.

4.2.2. Numerical Procedure for Softening-Elastic State

There are the softening region and the elastic region in the rock masses when the support pressure $p$ is between $p_{\mathrm{rp}}^{\ast }$ and $p_{\mathrm{pe}}^{\ast }$. The calculation for stress and displacement is implemented as the flowchart in Figure 5.

In Figure 5, the selecting method of the first values of $R_{\mathrm{p}1}$ and $R_{\mathrm{p}2}$ is the same as that in Figure 3. In order to get the first values of ${\gamma }_{\mathrm{p}1}^{R_0}$ and ${\gamma }_{\mathrm{p}2}^{R_0}$, calculate $\alpha \left(k\right)$ when ${\gamma }_{\mathrm{p}}^{R_0}\left(k\right)={\delta }_3+\left(k-1\right){\delta }_4$ and repeat this calculation from $k=1$ until the signs of $\alpha \left(k\right)$ and $\alpha \left(k+1\right)$ are contrary. Then take ${\gamma }_{\mathrm{p}}^{R_0}\left(k\right)$ and ${\gamma }_{\mathrm{p}}^{R_0}\left(k+1\right)$ as the initial values of ${\gamma }_{\mathrm{p}1}^{R_0}$ and ${\gamma }_{\mathrm{p}2}^{R_0}$ in the cycle computation. The value of ${\delta }_3$ and ${\delta }_4$ can be taken as ${10}^{-6}{\gamma }^{\mathrm{p}\ast }$ and ${10}^{-2}{\gamma }^{\mathrm{p}\ast }$ respectively, or smaller.

4.2.3. Numerical Procedure for Elastic-Softening-Residue State

The numerical procedure for the elastic-softening-residue state shown in Figure 6 adds the searching procedure for the radius of the residue region $R_{\mathrm{r}}$ by the bisection method, compared with the procedure in Figure 5. Once the initial value of the plastic shear strain on the position $R_0$ and the radii of the residue region and the plastic region are approximated adequately accurately by the method of bisection while the other initial value of the radial stress equals the support pressure, the stress and displacement of any point in the rock masses can be calculated by solving Equations (17), (20), (21), (23), and (26)–(30).

In Figure 6., the searching process for the first values of $R_{\mathrm{r}1}$ and $R_{\mathrm{r}2}$ is like the above selecting process. Similarly, calculate iteratively the plastic shear strain ${\gamma }_{\mathrm{p}}^{R_{\mathrm{r}}}\left(k\right)$ when $R=R_{\mathrm{r}}\left(k\right)=R_0+{\delta }_1+\left(k-1\right){\delta }_2$ from $k=1$ until ${\gamma }_{\mathrm{p}}^{R_{\mathrm{r}}}\left(k\right)>{\gamma }^{\mathrm{p}\ast }$ and ${\gamma }_{\mathrm{p}}^{R_{\mathrm{r}}}\left(k+1\right)<{\gamma }^{\mathrm{p}\ast }$ and then take $R_{\mathrm{r}}\left(k\right)$ and $R_{\mathrm{r}}\left(k+1\right)$ as the initial values.

The above five numerical procedures are integrated into the finite strain procedure.

5. Verification, Comparison, and Discussion

5.1. Verification for Finite Strain Procedure

For the verification, this study compares the proposed solution with the numerical solutions in [37,38]. The parameters of the MC rock mass listed in Table 1 are the same as those in ‘Figure 6’ provided by Ref. [37]. The relations between the dimensionless displacement $u_0/R_0$, $R_{\mathrm{p}}/R_0$, $R_{\mathrm{r}}/R_0$ and the dimensionless supporting pressure $p/P_0$ are visualized by the scattered diagrams in Figure 7 with two sets of data calculated by the procedure in this paper and former studies. Ref. [37] presents the results solved by the recursive program (Zhang’s solution) and obtained by the finite finite-difference method (FDM).

By the proposed procedure in this paper, the calculated limit support pressures $p_{\mathrm{rp}}^{\ast }$ and $p_{\mathrm{pe}}^{\ast }$ are equal to 0.0945 MPa and 0.2006 MPa respectively. As shown in Figure 7, the GRC and the evolutions of the softening region radius and residual region radius calculated by the proposed procedure are virtually identical with those of Zhang’s solution, which verifies the accuracy of the proposed numerical finite strain procedure. Without the supporting pressure, $u_0/R_0$, $R_{\mathrm{p}}/R_0$, and $R_{\mathrm{r}}/R_0$ are 0.3196, 1.884, and 2.260 respectively, slightly different with 0.3235, 1.87 and 2.25 in Zhang’s study. This difference is caused mainly by the approximant treatment of the infinity.

In Zhang’s study, $k_{\infty }$ is taken as 50 to approximate the infinite boundary, whereas $k_{\infty }$ equals ${10}^5$ in this section. As shown in Figure 8, the value of $k_{\infty }$ obviously affects the computational accuracy of $u_0/R_0$, $R_{\mathrm{p}}/R_0$ and $R_{\mathrm{r}}/R_0$. When $k_{\infty }$ reaches 1000, the calculation results stabilize. In addition, the iteration times $n$ have a decisive influence on the accuracy as well. As shown in Figure 9, with the number of iterations increasing, the values of $u_0/R_0$, $R_{\mathrm{p}}/R_0$ and $R_{\mathrm{r}}/R_0$ gradually converge and the magnitude of the calculation error $\left|\alpha \right|$ decreases to about ${10}^{-14}$ which reaches the highest accuracy when $k_{\infty }={10}^7$. After iterating only 10 times, $u_0/R_0$, $R_{\mathrm{p}}/R_0$ and $R_{\mathrm{r}}/R_0$ become rather stable. When the tolerances of the calculation relative error $\mathrm{Tol}\left(\alpha \right)$ are ${10}^{-6}$, ${10}^{-8}$, ${10}^{-10}$, ${10}^{-12}$, and ${10}^{-14}$, the time for one whole calculation is 0.89 s, 1.18 s, 1.72 s, 1.98 s, and 2.82 s in this numerical example. It shows that for the proposed procedure the significant improvement of the precision does not affect the calculation efficiency excessively.

5.2. Comparison for Small Strain Procedure

Likewise, this iteration method can be applied to solve the excavating problem in strain-softening rock masses under the small deformation hypothesis. The differential equations and the calculation formula for the displacement in the elastic region, the softening region, and the residue region are derived and listed in Appendix A. Following the procedures in Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6, the numerical procedure for circular tunnels excavated in strain-softening rock masses is programmed based on the small deformation hypothesis.

In Figure 10, solutions obtained by finite strain procedure and small strain procedure are compared with the authoritative Alonso’s self-similar solution [8], which is based on the small deformation hypothesis. The small strain solution is consistent with Alonso’s result, verifying the correctness and strong applicability of this iteration method. The finite strain solutions, $u_0$, $R_{\mathrm{p}}$ and $R_{\mathrm{r}}$, are smaller than the small strain solutions. That the gap between the two solutions decreases with the supporting pressure increasing, illustrates that the small strain hypothesis underestimates the self-bearing capacity of the surrounding rock and the finite deformation theory is more appropriate for the analysis of large deformation tunnels.

6. Application for Tunnel Design

During the construction of the Zhongyi tunnel in the Yunnan province of China, over 44% of the tunnel encountered large deformation and the maximum deformation was larger than 600 mm. The main reason for this problem was the mismatch between the support stiffness and the weak surrounding rock conditions [37,40]. This section adopts the finite strain procedure to calculate the needed stiffness of the whole supporting system. The three-cantered arch section is equivalent to a circular tunnel with a radius of 4.25 m. The in-situ stress and rock mass parameters are as follows: $P_0=12\mathrm{MPa},E=0.95\mathrm{GPa}$, $v=0.3,{\gamma }^{p\ast }=0.035,c_{\mathrm{p}}=0.48\mathrm{MPa},c_{\mathrm{r}}=0.25\mathrm{MPa},{\phi }_{\mathrm{p}}={26}^{°},{\phi }_{\mathrm{r}}={20}^{°},{\psi }_{\mathrm{p}}={16}^{°},{\psi }_{\mathrm{r}}={12}^{°}$ [37].

First, the finite strain procedure is integrated into one MATLAB toolbox called SAFE (Sensitivity Analysis For Everybody) [41] to perform the sensitivity analysis of tunnel deformation. The boundary of rock mass parameters is $\left(\mu -3\sigma , \mu +3\sigma \right)$, where $\mu$ and $\sigma$ denote the mean value and the standard deviation of each parameter. The distribution type of the rock mass parameters is normal distribution [42]. It is supposed that the coefficient of variations (COVs) of the in-situ stress and of other parameters are 0.1 and 0.05, respectively. The support pressure is sampled randomly and uniformly between 0 and $0.25P_0$. According to the global sensitivity indices (Sobol indices) of each parameter in Figure 11, tunnel deformation is most sensitive to the support pressure $p$ compared to other factors except for the inexorable geostress $P_0$. Therefore, improving the support stiffness is the most direct and effective means to control the tunnel deformation.

Furthermore, the stochastic analysis is performed by the Monte Carlo method to optimize the stiffness of the support structure for the Zhongyi tunnel. All distribution parameters are the same as the above but the range of the support pressure is $\left[0, 0.5P_0\right]$. The sampling number is 20,000 and the time for each calculation is between 0.9 s and 15 s in a computer with Intel Core I 7-6700 CPU. The tunnel displacement and support pressure are scattered in Figure 12a and the upper and lower envelope curves are depicted. As the support pressure increases, the variation range of $u_0$ becomes smaller. When the deformation limit is 0.5 m, the range of the support pressure is between 0.72 MPa and 1.63 MPa.

In the specific design process for the tunnel, calculating the needed stiffness is more practical than predicting the tunnel displacement. Hence, the combined stiffness of the primary support and the secondary lining is obtained by taking the displacement release coefficient as 0.25 before installing the support structure [43]. As shown in Figure 12b, the discrete degree of the displacement $u_0$ significantly decreases with the stiffness $k$ increasing. When the combined stiffness is larger than 4.3 MPa/m, it can be guaranteed that the tunnel displacement $u_0$ will not exceed 0.5 m, smaller than the limit deformation. As for the detail design, references [44,45,46,47,48,49,50,51,52,53] can be used to calculate the support stiffness, optimize the support parameters, analyze the tunnel stability and perform the whole stiffness design for the tunnel-support system.

7. Conclusions

This paper proposes a new numerical finite strain procedure for a circular tunnel in MC strain-softening rock masses with a non-associated flow rule. First, the explicit form of the differential relation between the original coordinate and the deformed coordinate is established in the Lagrangian coordinate. Using the governing equations, yield criterion, evolutional law, the solutions in the residue region, and the softening region are derived as two sets of differential equations of radial stress and plastic shear strain and the solutions in the elastic region are expressed as the differential equations of radial stress and tangential stress. By directly solving the differential equations using the Runge–Kutta method, the problem of solving for the finite strain solution transforms into the problem of searching for two adjacent boundaries of three regions ($R_{\mathrm{s}}$ and $R_{\mathrm{p}}$) and the initial values on the excavation boundary ($R_0$).

Then the bisection method is adopted to approximate the true value based on the characteristics that the dependent variables (radial stress, plastic shear strain and tangential stress) in each region increase or decrease monotonously with the deformed radius. Two numerical procedures of calculating $p_{\mathrm{rp}}^{\ast }$ and $p_{pe}^{\ast }$ and three procedures for the solutions in the residue-softening-elastic state, the softening-elastic state, and the elastic state respectively are programmed and integrated into the finite strain numerical procedure. The proposed finite strain procedure is validated by comparing with recursive solutions and FLAC simulation results. By increasing $k_{\infty }$ and $n$, the calculation relative error $\left|\alpha \right|$ decreases to about ${10}^{-14}$. When $k_{\infty }$ and $n$ are greater than 1000 and 10 respectively for this numerical example, $u_0/R_0$, $R_{\mathrm{p}}/R_0$ and $R_{\mathrm{r}}/R_0$ become rather stable. The proposed procedure directly solves three sets of differential equations and uses the bisection method to approximate the boundary conditions, realizing quantitative control for the calculation error. More meaningfully, the highly efficient computing performance of the proposed procedure creates conditions for its engineering application in the global sensitivity analysis and the reliability design of the Zhongyi tunnel, which requires large sample data.

Moreover, this paper develops a numerical procedure for circular tunnels excavated in strain-softening rock masses based on the small deformation hypothesis still using the iteration method. The correctness of this iteration method was verified again by comparing with a self-similar solution. However, this solution underestimates the self-bearing capacity of the surrounding rock and is inappropriate to be used to analyze a large deformation tunnel.

In the end, the finite strain procedure was adopted for the sensitivity analysis and the reliability design of Zhongyi tunnel. The sensitivity indices of tunnel deformation to support pressure and geostress are much higher than others. This study recommends 4.3 MPa/m as the designed support stiffness to guarantee a tunnel displacement lower than 0.5 m.

Author Contributions

Conceptualization, W.C. and D.Z.; methodology, W.C. and Q.F.; software, W.C.; validation, W.C.; formal analysis, W.C.; investigation, W.C., X.C. and T.X.; writing—original draft preparation, W.C.; writing—review and editing, D.Z. and Q.F.; visualization, W.C., X.C. and T.X.; supervision, D.Z.; project administration, D.Z. and Q.F.; funding acquisition, D.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported financially by the National Key Research and Development Program of China (Grant No. 2021YFB2300803 and Grant No. 2017YFC0805401) and National Natural Science of China (Grant No. 51738002).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within this article.

Acknowledgments

The authors are grateful to all editors’ and reviewers’ work which improved the results and the presentation of this paper. Further, the authors thank Francesca Pianosi for making the MATLAB toolbox SAFE [41] publicly available.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

MC	Mohr–Coulomb
DF	Dimensional factor
DA	Divergence angles

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. Tunnel excavation calculation model based on finite strain theory.

Figure 2. Calculation flowchart of [Forumla omitted. See PDF.].

Figure 3. Calculation flowchart of [Forumla omitted. See PDF.].

Figure 4. Calculation flowchart for elastic state.

Figure 5. Calculation flowchart for softening-elastic state.

Figure 6. Calculation flowchart for elastic-softening-residue state.

Figure 7. Comparison between the results of this paper and Zhang’s: (a) Ground response curve (GRC); (b) radii of residue region and softening region.

Figure 8. Influence of [Forumla omitted. See PDF.]: (a) [Forumla omitted. See PDF.] varying with [Forumla omitted. See PDF.]; (b) plastic radii varying with [Forumla omitted. See PDF.].

Figure 9. Influence of iteration times: (a) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] varying with [Forumla omitted. See PDF.]; (b) plastic radii and [Forumla omitted. See PDF.] varying with [Forumla omitted. See PDF.].

Figure 10. Comparison of finite strain solution, small strain solution, and Alonso’s self-similar solution: (a) [Forumla omitted. See PDF.] varying with [Forumla omitted. See PDF.]; (b) plastic radii varying with [Forumla omitted. See PDF.].

Figure 11. Sensitivity index of tunnel deformation.

Figure 12. Stochastic sampling results: (a) tunnel displacement and support pressure; (b) tunnel displacement and support stiffness.

Table 1

Geometry and material property parameters.

MC Rock	Value
Radius of opening, $R_0$ (m)	3
Initial stress, $P_0$ (MPa)	1
Young’s modulus, $E$ (MPa)	30
Poisson ratio, $v$ (-)	0.3
Critical plastic shear strain, ${\gamma }^{\mathrm{p}\ast }$ (-)	0.15
$c_{\mathrm{p}}$ (MPa)	0.2
$c_{\mathrm{r}}$ (MPa)	0.02
${\phi }_{\mathrm{p}}$ (Deg)	40
${\phi }_{\mathrm{r}}$ (Deg)	20
${\psi }_{\mathrm{p}}$ (Deg)	10
${\psi }_{\mathrm{r}}$ (Deg)	5

Appendix A

With the small deformation hypothesis, the formulas for stress and displacement of rock mass in three regions are derived and listed as follows. All formula symbols are consistent with the above formulas.

(1) Solutions in residue region

The solution of the radial stress and the tangential stress in residue region can be expressed as Equations (A1) and (A2). (A1)$${\sigma }_{\mathrm{r}}=\frac{{\sigma }_0\left(k_{\mathrm{r}}-1\right)+{\sigma }_{\mathrm{c}}}{k_{\mathrm{r}}-1}{\left(\frac{r}{R_0}\right)}^{k_{\mathrm{r}}-1}-\frac{{\sigma }_{\mathrm{c}}}{k_{\mathrm{r}}-1}$$ (A2)$${\sigma }_{\mathsf{\theta }}=\frac{k_{\mathrm{r}}\left[{\sigma }_0\left(k_{\mathrm{r}}-1\right)+{\sigma }_{\mathrm{c}}\right]}{k_{\mathrm{r}}-1}{\left(\frac{r}{R_0}\right)}^{k_{\mathrm{r}}-1}-\frac{{\sigma }_{\mathrm{c}}}{k_{\mathrm{r}}-1}$$ where $k_r=\frac{1+\sin {\phi }_r}{1-\sin {\phi }_r}$ and ${\sigma }_c\left({\gamma }^p\right)=\frac{2c_r\cos {\phi }_r}{1-\sin {\phi }_r}.$ The solution of the displacement in the residue region is expressed as Equation (A3). (A3)$$u=C\frac{r^{k_{\mathrm{r}}}-r^{-K}{R}_0{}^{k_{\mathrm{r}}+K}}{K+k_{\mathrm{r}}}+D\frac{r-r^{-K}{R}_0{}^{K+1}}{K+1}+r^{-K}{R}_0{}^K{u}_0$$ where $C=\frac{1}{2G}\frac{{\sigma }_0\left(k_{\mathrm{r}}-1\right)+{\sigma }_c}{\left(k_{\mathrm{r}}-1\right)R_0{}^{k_{\mathrm{r}}-1}}\left[1-v-Kv+k_{\mathrm{r}}\left(K-Kv-v\right)\right]$, $D=-\frac{\left(K+1\right)\left(1-2v\right)}{2G}\left(\frac{{\sigma }_{\mathrm{c}}}{k_{\mathrm{r}}-1}+P_0\right)$ and $K=\frac{1+\sin {\psi }_{\mathrm{r}}}{1-\sin {\psi }_{\mathrm{r}}}$.

(2) Solutions in softening region

The differential equations of the radial stress and of the plastic shear strain with respect to the original radius are derived as Equations (A4) and (A5). (A4)$$\frac{\mathrm{d}{\sigma }_{\mathrm{r}}}{\mathrm{d}r}=\frac{1}{r}\left[\left(k_{\mathrm{p}}\left({\gamma }^{\mathrm{p}}\right)-1\right){\sigma }_{\mathrm{r}}+{\sigma }_{\mathrm{c}}\left({\gamma }^{\mathrm{p}}\right)\right]$$ (A5)$$\frac{\mathrm{d}{\gamma }^{\mathrm{p}}}{\mathrm{d}r}=\frac{-\frac{1}{2G}\left[\left(1-v\right)k_{\mathrm{p}}\left({\gamma }^{\mathrm{p}}\right)-v\right]\frac{\mathrm{d}{\sigma }_{\mathrm{r}}}{\mathrm{d}r}-\frac{1}{r}\frac{1}{2G}\left\{\left[k_{\mathrm{p}}\left({\gamma }^{\mathrm{p}}\right)-1\right]{\sigma }_{\mathrm{r}}+{\sigma }_{\mathrm{c}}\left({\gamma }^{\mathrm{p}}\right)\right\}-\frac{{\gamma }^{\mathrm{p}}}{r}}{\frac{1}{2G}\left[\left(1-v\right){\sigma }_{\mathrm{r}}\frac{\partial k_{\mathrm{p}}\left({\gamma }^{\mathrm{p}}\right)}{\partial {\gamma }^{\mathrm{p}}}+\left(1-v\right)\frac{\partial {\sigma }_{\mathrm{c}}\left({\gamma }^{\mathrm{p}}\right)}{\partial {\gamma }^{\mathrm{p}}}\right]+\frac{\partial \left[\frac{{\gamma }^{\mathrm{p}}}{1+K\left({\gamma }^{\mathrm{p}}\right)}\right]}{\partial {\gamma }^{\mathrm{p}}}}$$

The displacement of rock mass material points in the softening region can be calculated by Equation (A6). (A6)$$u=\frac{r}{2G}\left[\left(1-v\right)\left({\sigma }_{\mathsf{\theta }}-P_0\right)-v\left({\sigma }_{\mathrm{r}}-P_0\right)\right]+r\frac{{\gamma }^{\mathrm{p}}}{1+K\left({\gamma }^{\mathrm{p}}\right)}$$

(3) Solutions in elastic region

The solution of the radial stress and the tangential stress in the elastic region can be expressed as Equations (A7) and (A8). (A7)$${\sigma }_{\mathrm{r}}=\frac{{\sigma }_{\mathsf{\theta }}^{R_{\mathrm{p}}}+{\sigma }_{\mathrm{r}}^{R_{\mathrm{p}}}}{2}-\frac{{\sigma }_{\mathsf{\theta }}^{R_{\mathrm{p}}}-{\sigma }_{\mathrm{r}}^{R_{\mathrm{p}}}}{2}\frac{R_{\mathrm{p}}{}^2}{r^2}$$ (A8)$${\sigma }_{\mathsf{\theta }}=\frac{{\sigma }_{\mathsf{\theta }}^{R_{\mathrm{p}}}+{\sigma }_{\mathrm{r}}^{R_{\mathrm{p}}}}{2}+\frac{{\sigma }_{\mathsf{\theta }}^{R_{\mathrm{p}}}-{\sigma }_{\mathrm{r}}^{R_{\mathrm{p}}}}{2}\frac{R_{\mathrm{p}}{}^2}{r^2}$$ where ${\sigma }_{\mathsf{\theta }}^{R_{\mathrm{p}}}$ and ${\sigma }_{\mathrm{r}}^{R_{\mathrm{p}}}$ are the radial stress and tangential stress on the boundary between the elastic region and the plastic region ($R_p$).

The displacement of rock mass in the elastic region can be calculated by Equation (A9). (A9)$$u=\frac{r}{2G}\left[\left(1-v\right)\left({\sigma }_{\mathsf{\theta }}-P_0\right)-v\left({\sigma }_{\mathrm{r}}-P_0\right)\right]$$

References

1. Zhang, D.L.; Sun, Z.Y.; Fang, Q. Scientific problems and research proposals for Sichuan-Tibet railway tunnel construction. Undergr. Space; 2021; in press [DOI: https://dx.doi.org/10.1016/j.undsp.2021.10.002]

2. Liu, W.W.; Chen, J.X.; Chen, L.J.; Luo, Y.B.; Shi, Z.; Wu, Y.F. Nonlinear deformation behaviors and a new approach for the classification and prediction of large deformation in tunnel construction stage: A case study. Eur. J. Environ. Civ. Eng.; 2020; 2, pp. 1-29. [DOI: https://dx.doi.org/10.1080/19648189.2020.1744482]

3. Chen, Z.M.; Zhao, D.A.; Yu, Y.Y. The interaction mechanism of surrounding rock and supporting structure in high geostress extrusion fault. Adv. Mat. Res.; 2011; 243–249, pp. 3588-3598. [DOI: https://dx.doi.org/10.4028/www.scientific.net/AMR.243-249.3588]

4. Meng, L.B.; Li, T.B.; Jiang, Y.P.; Wang, R.; Li, Y.R. Characteristics and mechanisms of large deformation in the Zhegu mountain tunnel on the Sichuan-Tibet highway. Tunn. Undergr. Space Technol.; 2013; 37, pp. 157-164. [DOI: https://dx.doi.org/10.1016/j.tust.2013.03.009]

5. Bian, K.; Liu, J.; Liu, Z.P.; Liu, S.G.; Ai, F.; Zheng, X.Q.; Ni, S.H.; Zhang, W. Mechanisms of large deformation in soft rock tunnels: A case study of Huangjiazhai Tunnel. Bull. Eng. Geol. Environ.; 2019; 78, pp. 431-444. [DOI: https://dx.doi.org/10.1007/s10064-017-1155-8]

6. Wang, F.N.; Guo, Z.B.; Qiao, X.B.; Fan, J.Y.; Li, W.; Mi, M.; Tao, Z.G.; He, M.C. Large deformation mechanism of thin-layered carbonaceous slate and energy coupling support technology of NPR anchor cable in Minxian Tunnel: A case study. Tunn. Undergr. Space Technol.; 2021; 117, 104151. [DOI: https://dx.doi.org/10.1016/j.tust.2021.104151]

7. Xu, Z.X.; Zhang, L.G.; Jiang, L.W.; Wang, K.; Zhang, G.Z.; Feng, T.; Wang, D.; Song, Z.; Yi, X.J.; Wang, Z.W. et al. Engineering geological environment and main engineering geological problems of Ya’an-Linzhi section of the SiChuan-Tibet railway. Adv. Eng. Mater.; 2021; 53, pp. 29-42. (In Chinese) [DOI: https://dx.doi.org/10.15961/j.jsuese.202100224]

8. Alonso, E.; Alejano, L.R.; Varas, F.; Fdez-Manin, G.; Carranza-Torres, C. Ground response curves for rock masses exhibiting strain-softening behaviour. Int. J. Numer. Anal. Methods Geomech.; 2003; 27, pp. 1153-1185. [DOI: https://dx.doi.org/10.1002/nag.315]

9. Park, K.H. Similarity solution for a spherical or circular opening in elastic-strain softening rock mass. Int. J. Rock Mech. Min. Sci.; 2014; 71, pp. 151-159. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2014.07.003]

10. Guan, Z.; Jiang, Y.; Tanabasi, Y. Ground reaction analyses in conventional tunnelling excavations. Tunn. Undergr. Space Technol.; 2007; 22, pp. 230-237. [DOI: https://dx.doi.org/10.1016/j.tust.2006.06.004]

11. Park, K.H.; Tontavanich, B.; Lee, J.G. A simple procedure for ground response of circular tunnel in elastic-strain softening rock masses. Tunn. Undergr. Space Technol.; 2008; 23, pp. 151-159. [DOI: https://dx.doi.org/10.1016/j.tust.2007.03.002]

12. Lee, Y.; Pietruszczak, S. A new numerical procedure for elasto-plastic analysis of a circular opening excavated in a strain-softening rock mass. Tunn. Undergr. Space Technol.; 2008; 23, pp. 588-599. [DOI: https://dx.doi.org/10.1016/j.tust.2007.11.002]

13. Wang, S.L.; Yin, X.T.; Tang, H.; Ge, X.R. A new approach for analyzing circular tunnel in strain softening rock masses. Int. J. Rock Mech. Min. Sci.; 2010; 47, pp. 170-178. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2009.02.011]

14. Zhang, Q.; Jiang, B.S.; Wang, S.L.; Ge, X.R.; Zhang, H.Q. Elasto-plastic analysis of a circular opening in strain-softening rock mass. Int. J. Rock Mech. Min. Sci.; 2012; 50, pp. 38-46. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2011.11.011]

15. Han, J.X.; Li, S.C.; Li, S.C.; Yang, W.M. A procedure of strain-softening model for elasto-plastic analysis of a circular opening considering elasto-plastic coupling. Tunn. Undergr. Space Technol.; 2013; 37, pp. 128-134. [DOI: https://dx.doi.org/10.1016/j.tust.2013.04.001]

16. Cui, L.; Zheng, J.J.; Zhang, R.J.; Lai, H.J. A numerical procedure for the fictitious support pressure in the application of the convergence–confinement method for circular tunnel design. Int. J. Rock Mech. Min. Sci.; 2015; 78, pp. 336-349. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2015.07.001]

17. Zou, J.F.; Li, C.; Wang, F. A new procedure for ground response curve (GRC) in strainsoftening surrounding rock. Comput. Geotech.; 2017; 89, pp. 81-91. [DOI: https://dx.doi.org/10.1016/j.compgeo.2017.04.009]

18. Wang, F.Y.; Qian, D. Difference solution for a circular tunnel excavated in strainsoftening rock mass considering decayed confinement. Tunn. Undergr. Space Technol.; 2018; 82, pp. 66-81. [DOI: https://dx.doi.org/10.1016/j.tust.2018.08.001]

19. Zhang, Q.; Zhang, C.H.; Jiang, B.S.; Li, N.; Wang, Y.C. Elastoplastic coupling solution of circular openings in strain-softening rock mass considering pressure-dependent effect. Int. J. Geomech.; 2018; 18, 04017132. [DOI: https://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0001043]

20. Xia, C.C.; Xu, C.; Liu, Y.P.; Han, C.L. Elastoplastic solution of deep buried tunnel considering strain-softening characteristics based on GZZ strength criterion. Chin. J. Rock Mech. Eng.; 2018; 37, pp. 2468-2477. (In Chinese) [DOI: https://dx.doi.org/10.13722/j.cnki.jrme.2018.0611]

21. GhorbaniHadi, A.; Hasanzadehshooiili, H. A comprehensive solution for the calculation of ground reaction curve in the crown and sidewalls of circular tunnels in the elastic-plastic-EDZ rock mass considering strain softening. Tunn. Undergr. Space Technol.; 2019; 84, pp. 413-431. [DOI: https://dx.doi.org/10.1016/j.tust.2018.11.045]

22. Alejano, L.R.; Rodriguez-Dono, A.; Alonso, E.; Fdez-Manín, G. Ground reaction curves for tunnels excavated in different quality rock masses showing several types of post-failure behaviour. Tunn. Undergr. Space Technol.; 2009; 24, pp. 689-705. [DOI: https://dx.doi.org/10.1016/j.tust.2009.07.004]

23. Alejano, L.R.; Rodriguez-Dono, A.; Alonso, E.; Fdez-Manín, G. Application of the convergence-confinement method to tunnels in rock masses exhibiting Hoek–Brown strain-softening behaviour. Int. J. Rock Mech. Min. Sci.; 2010; 47, pp. 150-160. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2009.07.008]

24. González-Cao, J.; Varas, F.; Bastante, F.G.; Alejano, L.R. Ground reaction curves for circular excavations in non-homogeneous, axisymmetric strain-softening rock masses. J. Rock Mech. Geotech. Eng.; 2013; 5, pp. 431-442. [DOI: https://dx.doi.org/10.1016/j.jrmge.2013.08.001]

25. Ketabian, E.; Molladavoodi, H. Practical ground response curve considering post-peak rock mass behavior. Eur. J. Environ. Civ. Eng.; 2015; 21, pp. 1-23. [DOI: https://dx.doi.org/10.1080/19648189.2015.1090928]

26. Zareifard, M.R. Ground reaction curve for deep circular tunnels in strain-softening Mohr-Coulomb rock masses considering the damaged zone. Int. J. Geomech.; 2020; 20, 10. [DOI: https://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0001822]

27. Zareifard, M.R. Ground response curve of deep circular tunnel in rock mass exhibiting Hoek–Brown strain-softening behaviour considering the dead weight loading. Eur. J. Environ. Civ. Eng.; 2021; 25, pp. 2509-2539. [DOI: https://dx.doi.org/10.1080/19648189.2019.1632745]

28. Durban, D. A finite strain axially-symmetric solution for elastic tubes. Int. J. Solids Struct.; 1988; 24, pp. 675-682. [DOI: https://dx.doi.org/10.1016/0020-7683(88)90016-9]

29. Yu, H.S.; Rowe, R.K. Plasticity solutions for soil behaviour around contracting cavities and tunnels. Int. J. Numer. Anal. Methods Geomech.; 1999; 23, pp. 1245-1279. [DOI: https://dx.doi.org/10.1002/(SICI)1096-9853(199910)23:12<1245::AID-NAG30>3.0.CO;2-W]

30. Vrakas, A.; Anagnostou, G. A finite strain closed-form solution for the elastoplastic ground response curve in tunnelling. Int. J. Numer. Anal. Methods Geomech.; 2014; 38, pp. 1131-1148. [DOI: https://dx.doi.org/10.1002/nag.2250]

31. Vrakas, A. A finite strain solution for the elastoplastic ground response curve in tunnelling: Rocks with non-linear failure envelopes. Int. J. Numer. Anal. Methods Geomech.; 2017; 41, pp. 1077-1090. [DOI: https://dx.doi.org/10.1002/nag.2665]

32. Chadwick, P. The quasi-static expansion of a spherical cavity in metals and ideal soils. Q. J. Mech. Appl Math.; 1959; 12, pp. 52-71. [DOI: https://dx.doi.org/10.1093/qjmam/12.1.52]

33. Park, K.H. Large strain similarity solution for a spherical or circular opening excavated in elastic-perfectly plastic media. Int. J. Numer. Anal. Methods Geomech.; 2015; 39, pp. 724-737. [DOI: https://dx.doi.org/10.1002/nag.2333]

34. Mo, P.Q.; Fang, Y.; Yu, H.S. Benchmark solutions of large-strain cavity contraction for deep tunnel convergence in geomaterials. J. Rock Mech. Geotech. Eng.; 2020; 12, pp. 596-607. [DOI: https://dx.doi.org/10.1016/j.jrmge.2019.07.015]

35. Zhang, Q.; Li, C.; Jiang, B.S.; Yu, L.Y. Elastoplastic analysis of circular openings in elasto-brittle-plastic rock mass based on logarithmic strain. Math. Probl. Eng.; 2017; 2017, 7503912. [DOI: https://dx.doi.org/10.1155/2017/7503912]

36. Guan, K.; Zhu, W.; Wei, J.; Liu, X.; Niu, L.; Wang, X. A finite strain numerical procedure for a circular tunnel in strain-softening rock mass with large deformation. Int. J. Rock Mech. Min. Sci.; 2018; 112, pp. 266-280. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2018.10.016]

37. Zhang, Q.; Wang, H.; Jiang, Y.; Lu, M.; Jiang, B. A numerical large strain solution for circular tunnels excavated in strain-softening rock masses. Comput. Geotech.; 2019; 114, 103142. [DOI: https://dx.doi.org/10.1016/j.compgeo.2019.103142]

38. Zhang, Q.; Wang, X.F.; Jiang, B.S.; Liu, R.C.; Li, G.M. A finite strain solution for strain-softening rock mass around circular roadways. Tunn. Undergr. Space Technol.; 2021; 111, 103873. [DOI: https://dx.doi.org/10.1016/j.tust.2021.103873]

39. Xu, C.; Xia, C.C. A new large strain approach for predicting tunnel deformation in strain-softening rock mass based on the generalized Zhang-Zhu strength criterion. Int. J. Rock Mech. Min. Sci.; 2021; 143, 104786. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2021.104786]

40. Cui, G.Y.; Qi, J.S.; Wang, D.Y. Research on large deformation control technology of tunnels in squeezing rock and its application. Sci. Prog.; 2020; 103, 003685042092316. [DOI: https://dx.doi.org/10.1177/0036850420923167]

41. Pianosi, F.; Sarrazin, F.; Wagener, T. A Matlab toolbox for Global Sensitivity Analysis. Environ. Model. Softw.; 2015; 70, pp. 80-85. [DOI: https://dx.doi.org/10.1016/j.envsoft.2015.04.009]

42. Bjureland, W. Reliability-Based Design of Rock Tunnel Support. Ph.D. Thesis; KTH Royal Institute of Technology: Stockholm, Sweden, 2020; Available online: https://www.diva-portal.org/smash/record.jsf?pid=diva2%3A1426755&dswid=6322 (accessed on 4 March 2022).

43. Zhang, Y.J.; Su, K.; Zhu, H.Z.; Qian, Z.D.; Wu, H.G. Installation time of an initial support for tunnel excavation upon the safety factors of surrounding rock. Appl. Sci.; 2020; 10, 5653. [DOI: https://dx.doi.org/10.3390/app10165653]

44. Wong, L.N.Y.; Zhang, D.L.; Fang, Q. Mechanical analysis of circular tunnels supported by steel sets embedded in primary linings. Tunn. Undergr. Space Technol.; 2013; 37, pp. 80-88. [DOI: https://dx.doi.org/10.1016/j.tust.2013.03.011]

45. Zhang, D.H.; Liu, S.H.; Ren, S.Q. Research on selection of steel and steel grid for tunnel support in soft rock with high geostress. Chin. J. Rock Mech. Eng.; 2014; 33, pp. 2258-2266. (In Chinese) [DOI: https://dx.doi.org/10.13722/j.cnki.jrme.2014.11.010]

46. Yu, F.C. Study on Mechanical Characteristics of the System of Support and Surrounding Rock. Ph.D. Thesis; Beijing Jiaotong University: Beijing, China, 2017; Available online: https://kns.cnki.net/KCMS/detail/detail.aspx?dbname=CDFDLAST2017&filename=1017064234.nh (accessed on 4 March 2022). (In Chinese)

47. Zhang, D.L.; Sun, Z.Y.; Fang, Q. Tunnel support structure system and its synergistic effect. Chin. J. Theor. Appl. Mech.; 2019; 51, pp. 577-593. (In Chinese) [DOI: https://dx.doi.org/10.6052/0459-1879-18-322]

48. Lei, S.Y.; Zhao, W. Study on the mechanism of circumferential yielding support for soft rock tunnel with large deformation. Rock Soil Mech.; 2020; 41, pp. 1039-1047. [DOI: https://dx.doi.org/10.16285/j.rsm.2019.5545]

49. Zhang, D.L.; Fang, H.C.; Chen, L.P.; Sun, Z.Y. Stiffness design theory for tunnel-support system. Chin. J. Rock Mech. Eng.; 2021; 40, pp. 649-662. (In Chinese) [DOI: https://dx.doi.org/10.13722/j.cnki.jrme.2020.0822]

50. Fang, Q.; Wang, G.; Yu, F.C.; Du, J.M. Analytical algorithm for longitudinal deformation profile of a deep tunnel. J. Rock Mech. Geotech. Eng.; 2021; 13, pp. 845-854. [DOI: https://dx.doi.org/10.1016/j.jrmge.2021.01.012]

51. Li, P.F.; Wei, Y.J.; Zhang, M.J.; Huang, Q.F.; Wang, F. Influence of non-associated flow rule on passive face instability for shallow shield tunnels. Tunn. Undergr. Space Technol.; 2022; 119, 104202. [DOI: https://dx.doi.org/10.1016/j.tust.2021.104202]

52. Sun, Z.Y.; Zhang, D.L.; Fang, Q. Analysis of the interaction between tunnel support and surrounding rock considering pre-reinforcement. Tunn. Undergr. Space Technol.; 2022; 115, 104074. [DOI: https://dx.doi.org/10.1016/j.tust.2021.104074]

53. Zhou, M.Z.; Fang, Q.; Peng, C. A mortar segment-to-segment contact method for stabilized total-Lagrangian smoothed particle hydrodynamics. Appl. Math. Model.; 2022; in press [DOI: https://dx.doi.org/10.1016/j.apm.2022.02.022]