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

: The small strain theory underestimates the self-bearing capacity of rock masses, especially for a soft rock tunnel under high geostress. To perform an efﬁcient and accurate calculation and provide a reference for the stiffness design of a tunnel, the ﬁnite strain solution for a circular tunnel in Mohr–Coulomb strain-softening rock masses with a non-associated ﬂow 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 ﬁnite strain solution was proposed to approximate the boundary condition at inﬁnity, the values of two adjacent boundaries, and the initial values on the excavation boundary. This numerical procedure was veriﬁed 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 − 8 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.


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 firstorder ordinary differential equations, for a spherical or circular opening in elasto-perfectlyplastic 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.

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 σ 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: After the excavation, the tunnel radius decreases from its initial value to the deformed radius . For a rock mass following elastic-strain-softening behaviour, the rock mass within a deformed radius reaches the residual state, the region outside a deformed radius stays in the elastic state and the rock mass between and is in the strain-softening state.

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).
where and 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: where and 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: where the shear modulus = /2/(1 + ), and are Young's modulus and Poisson's ratio, respectively; = 0 and = 0 in the elastic rock mass.

Softening Parameters, Yield Criterion, and Evolutional Law
The plastic shear strain is the most widely used softening index governing material parameters, which is defined as Equation (7). 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 r reaches the residual state, the region outside a deformed radius r p stays in the elastic state and the rock mass between r s and r p is in the strain-softening state.

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).
where σ r and σ θ 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: where ε θ and ε 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: where the shear modulus G = E/2/(1 + v), E and v are Young's modulus and Poisson's ratio, respectively; ε p r = 0 and ε p θ = 0 in the elastic rock mass.

Softening Parameters, Yield Criterion, and Evolutional Law
The plastic shear strain γ p is the most widely used softening index governing material parameters, which is defined as Equation (7).
Appl. Sci. 2022, 12, 2706 4 of 19 According to the non-associated MC potential function, the relation of radial and tangential plastic strains can be expressed as Equation (8).
where K(γ p ) = 1+sin ψ(γ p ) 1−sin ψ(γ p ) and ψ is the dilation angle. The failure of the rock mass is governed by the strain-softening Mohr-Coulomb yielding a function written as: ; c and ϕ 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 η p to the corresponding residual value η r as shown in Equation (10).
where η denotes any one of the strength parameters such as ψ, c, and ϕ etc.; γ p * is the critical plastic shear strain, the subscripts 'p' and 'r' of η denote the peak and residual values respectively.

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:

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: Deriving from the above equations, the stress and displacement solutions are expressed in the form of differential equations of the radial stress σ r , the plastic shear strain γ p , and the tangential stress σ θ .
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).
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).
From Equations (18) and (19), the second differential equation for the plastic shear strain regarding the initial radius is as follows: where k 1 = 1/(1 − sin ψ r ). 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).

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 dr/r is rewritten as follows: Substituting the tangential stress in Equation (7) and dr/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: 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.
dR . Hence, the differential equation for the plastic shear strain in the softening region is derived as Equation (26).
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).

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).
The displacement of the rock mass material point in the elastic region is calculated using Equation (30).

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 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 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 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.

Numerical Procedures for Calculating Limit Support Pressure
When the support pressure equals the softening limit value p * rp , the plastic shear strain on the tunnel inner boundary R 0 exactly equals the critical plastic shear strain γ p * . 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 * pe . 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 * rp , between p * rp and p * pe and greater than p * pe respectively. The flowcharts depicting the calculation procedures of p * rp and p * pe are shown in Figures 2 and 3 When the support pressure equals the softening limit value * , the plastic shear strain on the tunnel inner boundary exactly equals the critical plastic shear strain * . Comparatively, the plastic shear strain is equivalent to zero and the tangential stress and radial stress on the tunnel inner boundary just meet the yield criteria when the support pressure equals the elastic limit value * . 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 * , between * and * and greater than * respectively. The flowcharts depicting the calculation procedures of * and * are shown in Figure 2 and Figure 3, respectively.  In this paper, Tol p p , Tol( ), and are taken as to 10 , 10 , and 10 respectively, of which the rationality is verified in Section 5.
In Figure 2, the first values of and in the loop can be determined by the following method: (1) set ( ) = + + ( − 1) and calculate the plastic shear strain ( ) when = ( ) by solving Equations (17) and (20); (2) repeat the above calculation from = 1 until ( ) > 0 and ( + 1) < 0 and take ( ) and ( + 1) as and respectively. The value of and can be adjusted according to the accuracy and computational efficiency. For the stability of solutions, take and equal to 10 and 10 , respectively, in this paper.

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 is between * and , between * and , and smaller than * . The following sections elaborate the numerical procedures for the calculation of stress and displacement in the above three states.

Numerical Procedure for Elastic State
When the support pressure is greater than * and smaller than the geostress , 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.
Input the parameters of the tunnel and the excavated rock masses: 0 , 0 , , , r , p , r , p , r , p , p * .
Calculate r ∞ when = ∞ 0 .by solving Eqs. (28) and (29)   The error tolerances of α and γ R p p and the value of k ∞ have an important influence on the precision and accuracy of the calculation results. Ordinarily, greater k ∞ and smaller Tol γ R p p and Tol(α) lead to better precision and higher accuracy of the results. In this paper, Tol γ R p p , Tol(α), and k ∞ 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 p1 and R p2 in the loop can be determined by the following method: (1) set R p (k) = R 0 + δ 1 + (k − 1)δ 2 and calculate the plastic shear strain (17) and (20); (2) repeat the above calculation from k = 1 until γ R p p (k) > 0 and γ R p p (k + 1) < 0 and take R p (k) and R p (k + 1) as R p1 and R p2 respectively. The value of δ 1 and δ 2 can be adjusted according to the accuracy and computational efficiency. For the stability of solutions, take δ 1 and δ 2 equal to 10 −5 R 0 and 10 −4 R 0 , respectively, in this paper.

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 * pe and P 0 , between p * pe and P 0 , and smaller than p * pe . The following sections elaborate the numerical procedures for the calculation of stress and displacement in the above three states.

Numerical Procedure for Elastic State
When the support pressure is greater than p * pe 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  In Figure 4, the selecting method of the first values of and is like that of and in Section 3. At first, calculate when = .and ( ) = 0.1( − 1) by solving Equations (28) and (29) and repeat this calculation from = 1 until ( ) > and ( + 1) < . Then take ( ) and ( + 1) as the initial values of and in the loop respectively.

Numerical Procedure for Softening-Elastic State
There are the softening region and the elastic region in the rock masses when the support pressure is between * and * . 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 and is the same as that in Figure 3. In order to get the first values of and , calculate ( ) when ( ) = + ( − 1) and repeat this calculation from = 1 until the signs of ( ) and ( + 1) are contrary. Then take ( ) and ( + 1) as the initial values of and in the cycle computation. The value of and can be taken as 10 * and 10 * respectively, or smaller.
Input the parameters of the tunnel and the excavated rock masses: 0 , 0 , , , r , p , r , p , r , p , p * .

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 * rp and p * pe . 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 p1 and R p2 is the same as that in Figure 3. In order to get the first values of γ R 0 p1 and γ R 0 p2 , calculate α(k) when γ R 0 p (k) = δ 3 + (k − 1)δ 4 and repeat this calculation from k = 1 until the signs of α(k) and α(k + 1) are contrary. Then take γ R 0 p (k) and γ R 0 p (k + 1) as the initial values of γ R 0 p1 and γ R 0 p2 in the cycle computation. The value of δ 3 and δ 4 can be taken as 10 −6 γ p * and 10 −2 γ p * respectively, or smaller.

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 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). Appl. Sci. 2022, 12, x FOR PEER REVIEW 10 of 20

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 by the bisection method, compared with the procedure in Figure 5. Once the initial value of the plastic shear strain on the position 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)Page: 10 .
Input the parameters of the tunnel and the excavated rock masses: 0 , 0 , , , r , p , r , p , r , p , p * .
Take p 0 and r 0 as initial values of the differential equations of Eqs. (23) and (26) and calculate the stress and displacement in rock masses by Eqs. (23) (26) (27) (28) (29) and (30). In Figure 6., the searching process for the first values of R r1 and R r2 is like the above selecting process. Similarly, calculate iteratively the plastic shear strain γ R r p (k) when R = R r (k) = R 0 + δ 1 + (k − 1)δ 2 from k = 1 until γ R r p (k) > γ p * and γ R r p (k + 1) < γ p * and then take R r (k) and R r (k + 1) as the initial values.
The above five numerical procedures are integrated into the finite strain procedure. Appl. Sci. 2022, 12, x FOR PEER REVIEW 11 of 20 The above five numerical procedures are integrated into the finite strain procedure.

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 Input the parameters of the tunnel and the excavated rock masses: 0 , 0 , , , r , p , r , p , r , p , p * . Select 1 0 and 2 0 by test calculations.

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 p /R 0 , R 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). Table 1. Geometry and material property parameters.

MC Rock Value
Radius of opening, R 0 (m) 3 Initial stress, P 0 (MPa)  Figure 6' provided by Ref. [37]. The relations between the dimensionless displacement / , / , / and the dimensionless supporting pressure / 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 * and * are equal to 0.0945 MPa and 0.2006 MPa respectively. As shown in Figure 7 By the proposed procedure in this paper, the calculated limit support pressures p * rp and p * pe 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 p /R 0 , and R 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 ∞ is taken as 50 to approximate the infinite boundary, whereas k ∞ equals 10 5 in this section. As shown in Figure 8, the value of k ∞ obviously affects the computational accuracy of u 0 /R 0 , R p /R 0 and R r /R 0 . When k ∞ 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 p /R 0 and R r /R 0 gradually converge and the magnitude of the calculation error |α| decreases to about 10 −14 which reaches the highest accuracy when k ∞ = 10 7 . After iterating only 10 times, u 0 /R 0 , R p /R 0 and R r /R 0 become rather stable. When the tolerances of the calculation relative error Tol(α) 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. 10 . After iterating only 10 times, / , / and / become rather stable. When the tolerances of the calculation relative error Tol( ) are 10 , 10 , 10 , 10 , and 10 , 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.

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 Figures 2-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, , and , are smaller than the small strain solutions. That the gap between the two solutions decreases with the supporting pressure increasing, illus- 10 . After iterating only 10 times, / , / and / become rather stable. When the tolerances of the calculation relative error Tol( ) are 10 , 10 , 10 , 10 , and 10 , 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.

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 Figures 2-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, , and , are smaller than the small strain solutions. That the gap between the two solutions decreases with the supporting pressure increasing, illus-

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 Figures 2-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 p and R 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.
trates 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.

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: = 12MPa, = 0.95GPa, = 0.3, * = 0.035, = 0.48MPa, = 0.25MPa, = 26°, = 20°, = 16°, = 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 ( − 3 , + 3 ), where and 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.25 . According to the global sensitivity indices (Sobol indices) of each parameter in Figure 11, tunnel deformation is most sensitive to the support pressure compared to other factors except for the inexorable geostress . Therefore, improving the support stiffness is the most direct and effective means to control the tunnel deformation.

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: 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 (µ − 3σ, µ + 3σ), where µ and σ 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 [0, 0.5P 0 ]. 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. 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 [0, 0.5 ]. 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 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 significantly decreases with the stiffness increasing. When the combined stiffness is larger than 4.3 MPa/m, it can be guaranteed that the tunnel displacement 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.  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 [0, 0.5 ]. 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 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 significantly decreases with the stiffness increasing. When the combined stiffness is larger than 4.3 MPa/m, it can be guaranteed that the tunnel displacement 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. 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.

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 s and R 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 * rp and p * pe 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 ∞ and n, the calculation relative error |α| decreases to about 10 −14 . When k ∞ and n are greater than 1000 and 10 respectively for this numerical example, u 0 /R 0 , R p /R 0 and R 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 selfbearing 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. 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: where σ R p θ and σ R p r 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).