Stability of Unlined Elliptical Tunnels in Rock Masses

•A state of art approach to evaluate the elliptical tunnel stability of Hoek–Brown rock mass.•Rigorous upper bound and lower bound solutions of elliptical tunnel stability are derived using advanced finite element limit analysis.•Comprehensive design tables and equations are proposed for stability evaluation.


Introduction
The application of elliptical or quasi-rectangular tunnel sections has received much attention in recent years due to the need to reduce the volume of soils excavated as well as to better utilize tunnel space. The spaces at the top and the bottom of circular tunnels are indeed unneeded as the trains are straight-edged in shape. Elliptical (oval-shaped) tunnels are, therefore, preferred, and a new "composite circular" shield tunneling machine was first developed with an adjustable cutter to bore the elliptical tunnel in the Fukutoshin Line, a new subway line in Tokyo in 2008 (Akagi 2004).
In recent years, a new quasi-rectangular earth pressure balance (EPB) shield machine has recently been developed for Ningbo railway transit in China (Liu et al. 2021). The machine consists of 4 arcs and the lining ring is divided into 11 pieces. As stated by Liu et al. (2021), the filling and diffusion process of synchronous grouting in the shield tail void is significantly different from a circular shield and the quality of the tail grouting has a dominant effect on the overall stability during the construction. Indeed, when comparing with other standard tunnel sections, elliptical tunnels are more challenging to construct and the stability evaluation is of paramount importance to ensure the safety of such a unique tunnel shape (Amberg 1983;Hochmuth et al. 1987;Wone et al. 2003;Miura 2003;Miura et al. 2003).
Very few studies were reported in relation to the stability study of elliptical tunnels. Only recently in soil stability, Yang et al. (2015Yang et al. ( , 2016 studied the stability of unlined elliptical tunnel using the numerical upper bound method. An extended work of the stability solutions of elliptical tunnels considering internal normal stress along the tunnel periphery was later presented by Zhang et al. (2017) and Bhattacharya and Dutta (2021) for cohesionless and cohesive-frictional soils. For elliptical tunnel excavations in rock masses, to the best of our knowledge, there is none. To fill the research gap, this paper proposes a stability study for elliptical tunnel stability in rock mass using advanced upper and lower bound analyses using finite elements and nonlinear programming techniques.
For rock stability study, the Hoek-Brown (HB) model (Hoek and Brown 1980;Hoek et al. 2002) is a well-recognized failure criterion that includes the nonlinearity of the minor principal (compressive) stress. The formula of the HB failure criterion is in the form of a power-law relationship between the major and minor principal stresses (i.e., σ 1 and σ 3 ). Taking tensile normal stresses as positive, Eq. (1) describes the HB failure criteria in a mathematical form: where σ ci is the uniaxial compressive strength of intact rock mass and the parameters m b , s, and a are expressed in Eqs. (2)-(4):
(2)-(4), the geological strength index (GSI) has typical values from 10 to 100 (extremely poor rock mass to a perfectly intact rock mass). DF represents the degree of disturbance and it has typical values from 0 (undisturbed in-situ rock masses) to 1 (extremely disturbing in-situ rock masses). The parameter m i a material constant that is related to the frictional strength of an intact rock mass and has typical values from 5 to 30. Noting that these empirical parameters have been widely adopted by the rock mechanics community, in spite that geological observations are not always in line with the semi-empirical criterion.
Several researchers have recently studied the stability of tunnels in rock masses obeying the HB failure criterion. The rock stability of unlined circular was investigated by Zhang et al. (2019) and Keawsawasvong and Ukritchon (2020), while unlined square and rectangular tunnels in rock mass were studied by Xiao et al. (2019Xiao et al. ( , 2021 and Ukritchon and Keawsawasvong (2019a). Further, the rock stability of plane strain heading of tunnels was investigated by Ukritchon and Keawsawasvong (2019b). As stated before, the stability solutions of unlined elliptical tunnels in HB rock mass have never been presented in the literature. It is, therefore, the aim of this paper to produce useful stability charts and equations for geotechnical practitioners to estimate the stability of shallow unlined elliptical tunnels in rock masses obeying the Hoek-Brown failure criterion. Figure 1 shows the problem definition of an unlined elliptical tunnel in a rock mass. The tunnel has a cover depth of C. The shape of the tunnel is an ellipse with a horizontal dimension of B and a vertical dimension of D. The standard equation of an ellipse with center (0,0) and major axis parallel to the x-axis is presented in the following equation:

Statement of the Problem
where B = 2a and D = 2b. Equation (5) was used to generate the tunnel geometry throughout the parametric study in the paper. Figure 2a-c presents three typical model geometry so generated for a depth ratio C/D = 2 and different values of B/D = 0.5, 1, 2, respectively. The rock mass has a unit weight of γ, and at the surface area of rock masses, a uniform surcharge pressure σ s is applied over the area. It is hypothesized that no disturbance in the surrounding rock mass during tunnel excavation and the undisturbed  ( a,0) in-situ condition is applicable to the HB model with the disturbance factor DF be set to zero. This leads to the six design parameters in this study (i.e., B, D, σ ci , GSI, m i , γ). Using the dimensionless output parameter (σ s /σ ci ), where σ s is the uniform surcharge at collapse, Eq. (6) defines the stability factor (σ s /σ ci ) that can be written as a function of five dimensionless parameters: where B/D is the width ratio and C/D denotes the cover depth ratio. σ ci /γD represents the normalized uniaxial compressive strength ratio. σ s /σ ci is the stability factor. The parametric ranges considered in the current stability study are: (1). the cover depth ratio of tunnels is C/D = 1-5; (2). the width ratio is B/D = 0.5-2; (3). the yield parameter m i for the frictional strength of intact rock mass is m i = 5-30; (4). The geological strength index is GSI = 40-100, and the dimensionless uniaxial compressive strength ratio σ ci /γD is set to be 100-∞. Note that σ ci /γD = ∞ corresponds to a rock mass with extremely high strength (either σ ci is very large or γ is very small). The practical ranges of dimensionless parameters were chosen by following the previous research in the stability of other tunnel shapes in rock masses (e.g., Keawsawasvong and Ukritchon 2020; Ukritchon and Keawsawasvong 2019a, b;Xiao et al. 2019Xiao et al. , 2021Zhang et al. 2019).

Limitations
The rock masses in this study are defined to obey the Hoek-Brown (HB) failure criterion version 2002 (Hoek et al. 2002), where the considered parameters and their ranges include the yield parameter m i = 5-30, the geological strength index GSI = 40-100, and the dimensionless uniaxial compressive strength ratio σ ci /γD = 100-∞. The study assumes an excellent quality-controlled blasting or excavation by tunnel boring machine (TBM) during the tunnel construction. In this aspect, the disturbance factor (DF) of the confined rock mass surrounding a tunnel can be considered as very low or close to zero (DF = 0). More information on the guidelines for estimating the proper number of the disturbance factor (DF) can be found in Hoek et al. (2002). An example showing the effect of the disturbance factor on the tunnel stability can be found in Xiao et al. (2021). In addition, the cover depth ratio of tunnels is limited to the cases of C/D = 1-5 and the width ratio is limited to the cases of B/D = 0.5-2 in this study. It should be also noted that, when the uniaxial compressive strength is very small (or almost zero), it corresponding to highly fractured rock masses and the HB failure criterion does present certain limitation in the present stability analysis.

FELA Modelling
The proposed computational limit analysis is based on the theoretical work presented in Sloan (2013). It employs plastic bounding theorems, finite element discretization, and nonlinear programming. The associated upper and lower bound theorems (UB and LB) follow a perfectly plastic material with an associated flow rule. It is expected that a true stability solution can be bracketed from above (UB) and below (LB). A new computer program OptumG2 (OptumCE 2020) has been widely used to solve several geotechnical stability problems (e.g., Shiau and Smith 2006;Shiau et al. 2006a, b, 2021, Shiau and Al-Asadi 20182020a, b, c, 2021Ukritchon and Keawsawasvong 2017;Keawsawasvong and Ukritchon 2020). The program is, therefore, selected in this paper to compute the active collapse of unlined elliptical tunnels in rock masses using the latest powerful adaptivity meshing technique.
Using the model domains under plane strain conditions (see Fig. 2), three typical FELA meshes are shown in Fig. 3. The symmetrical plane on the left-hand side has a boundary condition that allows vertical movement only. This boundary condition is also applied to the right-hand side boundary of the domain. At the bottom plane, the boundary condition is set to be no movements in both vertical and horizontal directions. It is important to ensure that no error occurred in the computed bound solutions due to the problem size. Therefore, the sizes of the domain are carefully chosen to be so large that there is no intersection of the plastic shear zone at the far side and the bottom boundaries. The nodes on the periphery of the tunnel are free to move and there is no internal pressure presented inside the tunnel. The mesh adaptivity feature proposed by Ciria et al. (2008) is employed in this study. Using this powerful feature, a large number of elements can be automatically generated in zones where plastic shear strains are high. Throughout the study of the paper, five iterations of mesh adaptivity and 10,000 elements for the final mesh were set for all analyses.
The compressive uniform surcharge σ s is applied on the rock surface and it is the main objective function to be optimized at the active collapse state. The upper bound solution of the problem is obtained by solving the optimization problem that minimizes the surcharge (σ s ) (i.e., the objective function) and is subjected to the kinematically admissible velocity constraints. For the calculations of a lower bound solution, it is achieved by solving the optimization problem that maximizes the surcharge pressure (σ s ) (i.e., the objective function), while subjecting to the statically admissible stress constraints including equilibrium equations within the elements and along stress discontinuities, stress boundary conditions, and nowhere violation of the HB failure criterion. Once the surcharge pressure (σ s ) is computed for each parametric analysis, Eq. (6) is then used to calculate the stability factor (σ s /σ ci ).

Comparison of Results
Since the differences between UB and LB solutions obtained in this study are within an acceptable limit of 5%, all numerical results are hereafter presented as the average solutions (Ave) of the upper bound (UB) and lower bound (LB). In spite of this confidence in numerical results, Fig. 4 still shows a comparison of stability factor (σ s /σ ci ) between the present study and that of Keawsawasvong and Ukritchon (2020). The selected cases for the comparison are: B/D = 1; m i = 20; σ ci /γD = 100; C/D = 1, 2, 3, 4, 5; GSI = 40, 60, 80, 100. It is positive to see the excellent agreement between the two solutions, further giving great confidence before  producing the comprehensive design charts, tables, and equations in the next section.

Stability Charts and Tables
A total of 1,200 computed solutions (Ave) of the stability factor (σ s /σ ci ) are presented in Tables 1, 2, 3, 4, 5 for B/D = (0.5, 0.75, 1, 1.333, and 2). Numerical results of the stability factor are also presented graphically in Figs. 5,6,7,8,9 to study the individual effect of the considered parameters (e.g., B/D, C/D, σ ci /γD, GSI, and m i ). Figure 5a and b presents the effect of GSI on the stability factor (σ s /σ ci ) for various values of B/D. The data of the plots are for C/D = 3, σ ci /γD = 1000, and m i = 5 and 30, respectively. Numerical results have shown an exponential relationship between GSI and σ s /σ ci , where an increase of GSI results in a nonlinear increase of σ s /σ ci for all B/D. This observation is compatible with the exponential function of the HB failure criterion in Eqs. (2)-(4).
The effect of m i on σ s /σ ci is illustrated in Fig. 6a and b for the cases of GSI = 40 and 100, respectively. The plots are for (C/D = 3, σ ci /γD = 1000) and for five different width ratios B/D = (0.5, 0.75, 1, 1.333, and 2). Numerical results have shown a linear relationship between σ s /σ ci and m i . Figure 7a and b presents the effect of σ ci /γD on σ s /σ ci for m i = 5 and 30, respectively. The plot is for the case of C/D = 3 and GSI = 80, and it is concluded that σ ci /γD has little to no effect on the stability factor σ s /σ ci . This is not dissimilar to the undrained stability problems in soils, in which the normalized strength ratio has a negligible effect on the overall stability of rock tunnels too (Shiau and Al-Asadi 2018). The effect of cover depth ratio C/D on σ s /σ ci is shown in Fig. 8a and b, respectively, for m i = 5 and 30. A nonlinear relationship between C/D and σ s /σ ci exists. A larger C/D value results in greater stability (σ s /σ ci ). Finally, Fig. 9a and b depicts the influence of B/D on σ s /σ ci for m i = 5 and 30, respectively. The plots are for five different values of C/D = 1-5. Numerical results have shown a nonlinear decrease in the stability factor σ s /σ ci with the increasing B/D. A larger B/D ratio results in a smaller value of σ s /σ ci , which is in line with common engineering judgment.   Figure 10 compares the failure mechanisms of three B/D = 0.5, 1, and 2. In general, the failure zone of unlined elliptical tunnels is in the form of an elliptical shape beginning from the rock surface to the base of the tunnel. It is also interesting to see that, as B/D increases from 0.5 to 2, the ending point of the failure surface inside the tunnel has moved from the base of the tunnel to somewhere near the mid-height of the tunnel. Noting that the absolute values of shear dissipation are not important in such a perfect plasticity model, they are not presented in the figure. The depth ratio effects (C/D = 1, 2, 4, 5) on the overall failure mechanisms of the three various width ratios (B/D = 0.5, 1, and 2) are shown in Figs. 11, 12, 13, respectively. The plots are for σ ci /γD = 1000, GSI = 80, and m i = 20. A larger C/D results in a greater spreading of the failure zone around the tunnel, given the current study under the surcharge effects. For all cases presented, the failure zones extend to the base of tunnels when C/D ≥ 2.  A quick comparison of the stability factors σ s /σ ci for the different shapes of tunnels is shown in Fig. 14. Two published results are considered in the comparison, and they are for the square tunnel (Ukritchon and Keawsawasvong 2019a) and the plane strain heading (Ukritchon and Keawsawasvong 2019b). Together with our elliptical tunnels (B/D = 0.5 and 2) and circular tunnels (B/D = 1), numerical results have shown that the stability factors σ s /σ ci of plane strain heading is the largest, and it is then followed by the elliptical tunnel with B/D = 0.5, the circular tunnel B/D = 1, the square tunnel, and the elliptical tunnel with B/D = 2. This comparison figure is useful and is of great value to design practitioners in making engineering decisions.

Stability Equations
Using the average stability solutions (Ave) of upper and lower bounds and a curve fitting method, several design equations are developed to estimate the stability factor σ s / σ ci of shallow unlined horseshoe tunnels in Hoek-Brown rock masses. The mathematical form of approximate expressions based on Keawsawasvong and Ukritchon (2020) are expressed in the following equations: where (a i , b i , c i , d i , e i , f i , and g i ) are constant coefficients for the equations. The least-square method proposed by Sauer (2014) was used to determine the optimal values of the constant coefficients. Sauer's method can be used to minimize the sum of squares of the deviation (i.e., the error) in the stability factor (σ s /σ ci ) between the computed Ave solutions and the approximate solutions (i.e., the predictions). The comprehensive constants so obtained for the design Eqs. (7a)-(7d) are shown in Table 6 for width ratios B/D = 0.5, 0.75, 1, 1.33, and 2. The coefficient of determination (R 2 ) for each B/D equation is approximately 99.98%, indicating a very high precision of the developed constants and equations, which can be further used with great confidence by practitioners to estimate the stability factor σ s /σ ci of shallow unlined elliptical tunnels in rock masses.

Example
The practical use of the proposed equations is best explained using an example. The selected example of an elliptical tunnel has a horizontal dimension B = 6 m, a vertical dimension D = 3 m, a cover depth C = 3 m. The rock is found to have GSI = 50, m i = 17, σ ci = 63 MPa, and γ = 22 kN/m 3 . Determine the maximum surcharge pressure (σ s ) allowed before the tunnel reaches a collapse state.

Conclusions
To the best of the authors' knowledge, there was no published stability solution for unlined elliptical tunnels in Hoek-Brown rock masses. With the advanced adaptive meshing technique, finite elements, and nonlinear programming, this paper has successfully studied the stability of unlined elliptical tunnel in Hoek-Brown rock mass under the effect of surcharge pressure. Both the upper and lower bound limit analyses were used to solve for the stability solutions of a wide range of geometrical and Hoek-Brown material parameters. Using the average bound solutions, new design equations for computing the stability factors were developed using a least-square method. The proposed design equations are accurate with the coefficient of determination R 2 = 99.98%. This paper provides information that can be used as a reference at the preliminary design stage.