Stability of Spherical Cavity in Hoek–Brown Rock Mass

A state of art approach to evaluate the stability of Hoek–Brown rock mass with a spherical cavity. Rigorous upper bound and lower bound solutions of stability factors are solved using advanced finite element limit analysis. Comprehensive design charts, tables and equations are presented for stability evaluation. A state of art approach to evaluate the stability of Hoek–Brown rock mass with a spherical cavity. Rigorous upper bound and lower bound solutions of stability factors are solved using advanced finite element limit analysis. Comprehensive design charts, tables and equations are presented for stability evaluation.


Introduction
A karst cavity is a type of underground cavities in rock masses that can be formed from the dissolution of soluble rocks such as limestone, dolomite, and gypsum (Huang et al. 2017). Collapse of a karst cavity is a serious threat to human lives as well as an economic or environmental disaster. Therefore, obtaining information about karst cavity sizes and locations as well as the stability assessment of a spherical opening is of paramount importance to assure the safety of the existence of underground cavities.
The failure patterns of underground cavities in soils were investigated by Craig (1990) and Abdulla (1995) using centrifuge model tests as well as numerical and analytical techniques (e.g., Drumm et al. 1990; Abdulla and Goodings 1996;Tharp 1999;Vaziri et al. 2001;Augarde et al. 2003;Keawsawasvong and Ukritchon 2019;Shiau and Al-Asadi 2020;Shiau et al. 2021;Keawsawasvong 2021). The finite element limit analysis (FELA) is a powerful numerical technique that has been widely used by many researchers to determine the upper bound (UB) or the lower bound (LB) solutions of plastic collapse load of various stability problems based on the plastic bound theorems (Sloan 2013; Drucker et al. 1952). By using FELA, exact collapse loads can be accurately bracketed from UB and LB solutions. For the stability analysis of underground cavities, several researchers employed FELA to numerically solve the solutions to this problem (Augarde et al. 2003;Smith 2006, Keawsawasvong andUkritchon 2019;Shiau et al. 2016a, b;Keawsawasvong 2021). However, their solutions are limited to the cases of underground cavities in soils obeying the Mohr-Coulomb failure criterion. It is well known that rock masses are discontinuous materials with joints and fractures. The curved shape of the true failure envelope of rock masses cannot be replaced by a linear expression of the Mohr-Coulomb failure criterion. Consequently, it is seldom used to predict the stability of underground cavities in rock masses, in spite that the Mohr-Coulomb failure criterion constitutes a good approximation of the true failure envelope for some weak rock masses (Hoek et al. 2002).
The Hoek-Brown (HB) failure criterion (Hoek et al. 2002) has been widely accepted and adopted by many geotechnical engineers to predict the failure of intact rocks. It is an empirical failure criterion and the basic idea of the Hoek-Brown criterion was to start with the properties of intact rock, and then the factors to reduce those properties 1 3 are added to the equation because of the existence of joints in the rock (Hoek et al. 2002). By adopting the HB failure criterion into FELA, many researchers have investigated the stability of tunnels and underground openings in HB rock masses under plane strain conditions (e.g., Fraldi and Guarracino 2009;Ukritchon and Keawsawasvong 2019a, b;Keawsawasvong and Ukritchon 2020;Rahaman and Kumar 2020;Xiao et al. 2018Xiao et al. , 2019Xiao et al. , 2021Wu et al. 2020;Zhang et al. 2019).
Currently, there is no solution or stability criterion published for stability assessment of cavities in HB rock masses under axisymmetric conditions in the literature. This technical note aims to develop and propose a stability criterion for stability assessment of cavities in HB rock masses based on UB and LB solutions obtained from axisymmetric FELA. The considered parameters include the cover-depth ratio of cavities and the HB material parameters that have significant influences on the normalized collapse pressure applied at the rock surface above a cavity. Nonlinear regression analysis is employed to develop a closed-form approximate equation of this problem. This developed equation is valuable for engineers in practice to estimate the stability of underground cavities in rock masses. Figure 1 shows the problem definition of a spherical cavity in a rock mass. The cavity has a diameter (D) and a cover depth (C). A uniform surcharge (σ s ) at the collapse is applied over the surface area. The stability problem is investigated under 2D axisymmetric conditions. The work stated herein used the Hoek-Brown (HB) failure criterion (Hoek et al. 2002) to investigate the failure of a cavity in a rock mass. The Hoek-Brown (HB) parameters for a rock mass include σ ci , GSI, and m i , and a unit weight of γ.

Problem Statement and FELA Modelling
The expression of the HB failure criterion is in the form of a power-law relationship between the effective major and minor principal stresses (σ 1 and σ 3 ) as shown in Eq. (1).
It is to be noted that the compression negative sign convention applies to Eq. 1 and σ ci denotes the uniaxial compressive strength of intact rock masses. Other parameters such as m b , s, and α are expressed in Eqs. (2)-(4).
In the HB failure criterion, GSI is the geological strength index describing the quality of an in-situ rock mass. A GSI of 10 represents an extremely poor rock mass whilst 100 is used for intact rock. m i is the parameter used to describe the frictional strength of the intact rock mass. Noting that DF is the disturbance factor that has the range of 0-1, an undisturbed in-situ rock mass with DF = 0 is studied throughout the paper.
Using the concept of dimensionless ratios for practical design purposes, the stability solutions are determined through the use of five dimensionless variables as shown in Eq. (5).
where σ s /σ ci is the normalized collapse surcharge; C/D is the cover depth ratio; σ ci /γD is the normalized uniaxial compressive strength.
(1) The computer program OptumG2 (FELA, OptumCE 2019) is employed to perform the numerical analyses of the upper bound (UB) and lower bound (LB) finite element limit analysis (FELA). The FELA is based on the plastic bound theorems for a perfectly plastic material with an associated flow rule in conjunction with the finite element discretization and the mathematical optimization (Sloan 2013). The results from FELA include the UB and LB solutions that can bracket the true limit load from above and below.
In UB FELA, the rock mass is discretized by using sixnoded quadratic triangular elements to describe the overall velocity fields. In LB FELA, the rock mass is modelled using three-noded triangular elements to describe the linear stress field. The UB and LB solutions of this problem are computed by solving the optimization problem that minimizes (for the UB method) or maximizes (for LB method) the active surface pressure (σ s ) i.e., the collapse pressure at the ground surface. The mesh adaptivity technique (e.g., Ciria et al. 2008) is a powerful feature for improving LB and UB solutions. By activating this feature, more elements are added to the sensitive regions with large shear strain gradients at any iteration step, aiming to bridge the differences between UB and LB solutions. Five iterations of mesh adaptivity were used for all UB and LB simulations in the study, with 5000-10,000 elements in all analyses. It is interesting to note that the current technique reveals the location of a possible failure mechanism at the final stage of mesh adaptivity. Figure 2 shows a typical domain for the analysis of a spherical cavity. The left boundary is the plane of axisymmetry where only vertical movements can take place. The same condition applies to the right boundary. Nevertheless, velocities are fixed in both vertical and horizontal directions at the bottom boundary. At the rock surface, there is a uniform surcharge σ s applied over the surface area. The size of the domain is chosen to be large enough to avoid any interferences due to boundary effects. The current model does not allow internal pressure inside the cavity.

Results and Discussions
The chosen ranges of dimensionless parameters for the study are for σ ci /γD = 100-∞, GSI = 40 − 100, m i = 5-30 and C/D = 1-5. A total of 320 computed LB and UB solutions are obtained. and 6 for demonstrating the effects of C/D, GSI, m i , and σ ci / γD respectively. Figure 3a-d show that σ s /σ ci increases as the cover depth C/D increases. This trend is for σ ci /γD = 100 and all GSI values presented. In general, the larger the value of m i is, the greater the increase rate. The effect of GSI on σ s /σ ci is shown in Fig. 4, where σ ci /γD = 100 and m i = 5, 10, 20, and 30. The figure shows a nonlinear relationship between GSI and σ s /σ ci . An increase in GSI results in an exponential increase in σ s / ci , as reflected in the mathematic equations of HB failure criterion in Eqs. (2)-(4).
The effect of m i on σ s /σ ci is illustrated in Fig. 5a-d for four different values of GSI = 40, 60, 80, and 100. The study is for σ ci /γD = ∞ and five cover depth ratios C/D = 1-5. In general, σ s /σ ci increases linearly as m i increases for all C/D. Note that as C/D increases, the gradient of the linear line also increases. Figure 6 shows the effect of σ ci /γD on σ s /σ ci . This is for the cases of m i = 20. It is clear that the effect σ ci /γD on σ s /σ ci is insignificant for all C/D, given all the horizontal lines in the figure. Although it is physically impossible to have a weightless rock, the obtained results in Figs. 5 and 6 for σ ci /γD = ∞ are simply representatives of a very "strong" rock (i.e., very large strength ratio). Numerically, we can either put a very large value of σ ci or a very small γ to achieve the numerical results, which are needed to develop the stability criterion in the next section. Figure 7 presents upper bound adaptive meshes of four various cover depth ratios C/D. The chosen rock mass is for GSI = 80, m i = 10, and σ ci /γD = 100. It is important to note that the automatic adaptive meshing technique utilizes shear power dissipation as the control variable for re-meshing estimation. The number of elements in areas with very high shear power dissipation is automatically increased through successive iterations using the adaptive technique. Although the final adaptive mesh so produced resembles a so-called failure mechanism, it is common to use the contour plots of shear power dissipation to depict possible failure mechanisms of a soil structure, as it provides a good indicator of the intensity of non-zero plastic strains. Technically speaking, the actual values of the contour are not important in a perfectly plastic material model using FELA, and therefore the contour bars for

Stability Criterion
Design tables and charts of σ s /σ ci have been presented in the earlier sections. Nevertheless, in most cases, the limit pressure σ s is greater than the unconfined compressive strength of the intact rock σ ci . To obtain solutions with small values of σ s /σ ci , an approximate expression for calculating collapse pressures at the rock surface above a cavity is developed by using a curve fitting method. The average values of UB and LB are employed to determine an appropriate mathematical expression. The proposed stability criterion is presented in Eq. (6). Both N c and N γ are to be determined using Eqs. (7)-(11) with known values of C/D, GSI, m i , and σ ci /γD and g i are constant coefficients that were determined by performing a least square method (Sauer 2014). The optimum values of these constant coefficients are shown in Table 2. The value of R 2 of the proposed new stability criterion is about 99.98%, meaning that the approximation from Eq. (6) fits the FELA results very well. Figure 9 shows that N γ is a function of C/D only and the relationship between the two is linear. Note that N γ is negative for C/D < 2.5. This means that for cavities at relatively small depths, the self-weight of the ground has a positive effect on the limit load (i.e., σ s increases with γ). Intuitively, this is not right. Nevertheless, by taking another look at Eqs. (7, 9, and 10) would have explained this. Noting that the values of (F 1 , F 2 , and N c ) become small when C/D < 2.5, it may result in a total decrease of σ s (see Eq. 6). On the other note, Fig. 10 shows that N c is a function of both C/D and GSI for the given values of m i = 5, 10, 20, and 30. In general, N c increases nonlinearly as the cover depth C/D increases for all GSI values selected in the figures. The rate of increase is different for each GSI value. The curve becomes flattered (gradient decrease) as GSI decreases.

Conclusions
This short technical note has successfully studied the stability of spherical cavity in axisymmetric Hoek-Brown rock mass using the rigorous upper and lower bound finite element limit analysis. The solution was formulated to find the limit normalized surface pressures σ s /σ ci that is a function of four dimensionless parameters; namely the coverdepth ratio C/D, the Geological Strength Index GSI, the Hoek-Brown m i parameter, and the normalized uniaxial compressive strength σ ci /γD. A new stability criterion for predicting the stability of cavities in rock masses is developed by using a least square method of the computed solutions. The main findings of the present study are summarized as follows.
• The limit normalized surface pressures σ s /σ ci increases as the cover depth ratio C/D increases. The greater the values of GSI and m i , the larger the σ s /σ ci. In addition, the effect of σ ci /γD on σ s /σ ci is insignificant for all considered depth ratios in this study. • The failure mechanism of a cavity resembles a chimney type of vertical slippage when C/D is small. The lateral size of the failure mechanism extends when C/D increases. • The new cavity stability factors N c and N γ for the stability of cavities in rock masses are proposed in this paper, where N c is a function of C/D only while N γ is a function of C/D, GSI, and m i .  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.