An Improved Incremental Procedure for the Ground Reaction Based on Hoek-Brown Failure Criterion in the Tunnel Convergence-Confinement Method

Abstract

The purpose of this study is to introduce the Hoek–Brown nonlinear failure criterion into the convergence–confinement method (CCM), which is based on the linear failure criterion, and implant it into the ground reaction curve (GRC) in this theory. Based on consistent and rigorous theoretical analysis and according to the mechanical relationships, the principle of equilibrium, the material constitutive law, and the strain/displacement compatibility equation are used to deduce the closed-form analytical solution of the stress/displacement of the rock mass in the plastic region and to complete the simulation and interpretation of the nonlinear behavior of surrounding rock due to the advancing excavation of a circular tunnel. In this paper, confinement loss is regarded as the concept of increment in numerical analysis, and the calculation process and analysis steps of the explicit analysis method (EAM) are proposed. This method can not only realize the calculation of analytical solutions but can also be performed using simple spreadsheets. In addition, to verify the feasibility of this calculation method, this study uses published data to verify the validity of the analytical solution implemented by the incremental procedure and to study the influence of nonlinear failure criteria on the GRC, especially the mechanical behavior at the intrados of the tunnel and the distributions of tress/displacement on the cross-section of the tunnel. A comparison between the published results and the proposed solution in this study reveals consistent and favorable trends.

1. Introduction

The convergence–confinement method (CCM) is a systematic approach used in the design and construction of tunnels in rock masses. It recognizes that rock masses tend to deform and converge around an excavated tunnel due to the redistribution of stress. The method aims to control and manage this convergence by providing adequate confinement or support to ensure the stability and safety of the tunnel [1,2,3]. The method provides a framework for understanding and managing the convergence behavior of the rock mass, while the support design translates this understanding into practical measures to enhance tunnel stability [4,5]. The fundamental concept of the CCM is based on the interaction between the primary and secondary stress fields induced by tunnel excavation. When a tunnel is excavated, the primary stress in the rock mass is redistributed, causing the rock to deform and converge toward the opening. This convergence can lead to displacements, deformations, and potentially unstable conditions if not properly controlled [6,7,8,9].

The CCM focuses on achieving three main objectives, as shown in Figure 1: (1) The ground reaction curve (GRC) aims to limit the amount of convergence and control the rate at which it occurs. By controlling convergence, the risk of tunnel deformation and failure is reduced. The level of acceptable convergence depends on factors such as the tunnel purpose, ground conditions, support system, and safety requirements. Much of the research literature is devoted to this approach to developing mathematical and empirical expressions for ground reaction curves under different behavioral assumptions [10,11,12,13,14,15]. (2) Support confining curve (SCC) refers to the application of external support to the tunnel walls to counteract the inward movement and maintain stability. This is achieved through the installation of various support systems such as rock bolts, shotcrete, steel ribs, or other reinforcement elements. The support elements help confine the rock mass and distribute the stresses, reducing convergence and enhancing tunnel stability. Many studies on SCC have dealt with the limitations of nonlinearity of structural support, time dependence, progressive hardening, and transient conditions. Under different assumptions, the current research results have achieved a solution for the support design between SCC and GRC in the final equilibrium state (as shown at point E in Figure 1) [16,17,18,19,20,21,22,23].

(3) The confinement loss curve (CLC) emphasizes the importance of an observational approach during tunnel construction. This means continuously monitoring the behavior of the rock mass and adjusting the support measures accordingly. The value of the confinement loss can be obtained through the relationship between the unsupported distance from the tunnel face and the convergence measurement data of the tunnel. This value can be an incremental value, a positive scalar, or a function that takes into account the effect of advancing excavation on the tunnel face. This CLC principle shows that during the continuous advancement of the tunnel face, the confinement loss increases, resulting in a decrease in the radial stress around the tunnel and an increase in the radial displacement. Recently, much research literature has mentioned the effect of confinement losses around tunnels and constitutive models using numerical procedures to improve the accuracy of CCM analysis [24,25,26,27,28]. By implementing the convergence–confinement method, engineers can effectively manage the convergence behavior of rock masses in tunneling projects. It provides a systematic and proactive approach to ensuring the stability and safety of underground excavations [29,30]. In recent years, many scholars have proposed related numerical simulation research, such as the hybridized FSM and DDM stress analysis of the discontinuous surrounding rock of the tunnel, the mixed indirect boundary element method, and the numerical simulation of rock failure analysis in the crushing zone [31,32,33].

The ground reaction curve (GRC) is a concept used in the convergence–confinement method to understand the response of the rock mass and support system to tunnel excavation. It represents the relationship between the convergence (displacement) of the tunnel wall and the applied support pressure. In addition, the failure criterion provides insight into the strength and behavior of the rock mass, which influences the shape and characteristics of the ground reaction curve. In this study, the characteristics of the GRC are investigated and implemented by using the non-linear failure criterion of generalized Hoek–Brown [34,35,36]. The criterion allows engineers to assess the stability of rock slopes or tunnels under different stress conditions, taking into account the strength and quality of the rock mass. This criterion helps engineers evaluate the stability of the rock mass and determine the support requirements during tunneling [37,38,39,40]. By using the generalized Hoek-Brown failure criterion, engineers can estimate the rock mass strength and predict the potential failure mechanisms that may occur during tunneling. This information is crucial in selecting appropriate support systems and determining the required support pressure to control convergence and maintain tunnel stability [41,42,43,44]. Finally, the generalized Hoek–Brown failure criterion provides a basis for understanding the strength and failure behavior of the rock mass, while the ground reaction curve provides a means to evaluate the actual response of the rock mass and support system during tunneling. The failure criterion informs the design of the support system, and the ground reaction curve validates its effectiveness. Together, these concepts contribute to the overall assessment and management of convergence and confinement in the convergence-confinement method for tunneling projects.

The novelty and improvement of the GRC in this paper is to adopt the Hoek–Brown nonlinear failure criterion and insert it into the CCM theoretical analysis to replace the traditional Mohr–Coulomb linear failure criterion. Its advantage is that it can more realistically reflect the nonlinear mechanical behavior of the rock mass in the failure stage, especially in response to the plastic behavior of the surrounding rock. The purpose of this study is to use the theoretical analysis and basic principles of the convergence–confinement method to comprehensively discuss the mechanical behavior of tunnel excavation-induced ground reaction by considering the generalized Hoek–Brown nonlinear failure criterion. It not only provides a complete set of equations and details of their derivation, verification, and comparison but also illustrates the influence and availability of this nonlinear failure criterion based on the analysis results of the stress/displacement distribution on the cross-section and at the intrados of the tunnel.

2. Related Problem Description

2.1. Relationship between Tunnel Advancing Excavation Effect and Confinement Loss

Due to the advancing excavation of the tunnel face, the surrounding rock converges and the stress redistributes. That is to say, the impact of the unsupported distance on the advancing excavation of the working face is to consider the changes in stress and displacement caused by a certain observation section in the front after the continuous excavation of the tunnel. Based on the relevant basic research [6,7], using the tunnel convergence measurement data, a statistical regression analysis is used to establish a function related to the effect of tunnel face advancing excavation and to estimate the confinement loss. The function of this confinement loss curve (CLC) is shown in Figure 2. At a certain distance (z) from the working face, the definition of confinement loss (λ_z) can be given as

$${\lambda }_z={\lambda }_0+\left(1-{\lambda }_0\right)f\left(\frac{z}{R}\right)$$

(1)

For the assumption of the effect function of advancing excavation in the tunnel face, f(z/R) can also be called the advancing effect function, which can be expressed by the following:

$$f\left(\frac{z}{R}\right)=1-{\left(\frac{m}{m+\frac{z}{R}}\right)}^2$$

(2)

where R is the tunnel excavation radius, λ₀ is a confinement loss at the working face (z = 0, point A in Figure 1), and can be given as:

$${\lambda }_0=1-{\left(\frac{m}{m+\frac{d}{R}}\right)}^2$$

(3)

2.2. Stress Variation around a Circular Tunnel

In the analysis of the convergence–confinement method, the stress changes in the surrounding rock are generated due to the advancing excavation of the tunnel. The change in this stress can be described by the concept of stress gradient, which is expressed as the difference between the far-field stress and the near-field stress around the tunnel. The assumption of increment in numerical analysis is to treat the confinement loss as part of an increment or as a fraction of the stress gradient. This is an important principle and concept used in the theory of the convergence-confinement method. The change in this stress can be described by the following expression [6,7]:

(1) The stress state of the far-field (r –› ∞): Before the excavation of the tunnel, the stress of the surrounding rock is the initial in situ stress or the initial isotropic stress (as shown in Figure 3). Its stress can be expressed as follows:

$${\sigma }_r^i={\sigma }_{\theta }^i={\sigma }_v$$

(4)

where ${\sigma }_r^i$, ${\sigma }_{\theta }^i$, and ${\sigma }_v$ are the initial radial stress, the initial tangential stress, and the vertical stress.

(2) The stress state of the near-field (R ≤ r < ∞): After the tunnel is excavated, the stress change of the surrounding rock can be obtained through Kirsch’s solution, as shown in the following:

$${\sigma }_r^f=\left(1-\frac{R^2}{r^2}\right){\sigma }_v$$

(5)

$${\sigma }_{\theta }^f=\left(1+\frac{R^2}{r^2}\right){\sigma }_v$$

(6)

where ${\sigma }_{\theta }^f$ and ${\sigma }_r^f$ are the final tangential stress and the final radial stress around the tunnel proximity, respectively.

Assuming that the stress gradient produced by the advancing excavation of the tunnel face is the difference between the far-field stress and the near-field stress, the changes in the tangential stress and radial stress can be expressed as follows:

$${\overline{\sigma }}_r={\sigma }_r^f-{\sigma }_r^i=\left(1-\frac{R^2}{r^2}\right){\sigma }_v-{\sigma }_v$$

(7)

$${\overline{\sigma }}_{\theta }={\sigma }_{\theta }^f-{\sigma }_{\theta }^i=\left(1+\frac{R^2}{r^2}\right){\sigma }_v-{\sigma }_v$$

(8)

Regarding the assumption of increment in numerical analysis, the stress increment can be expressed by multiplying the stress gradient by the confinement loss, and its relationship can be expressed as follows:

$$∆{\sigma }_r=\lambda {\overline{\sigma }}_r=-\lambda \frac{R^2}{r^2}{\sigma }_v$$

(9)

$$∆{\sigma }_{\theta }=\lambda {\overline{\sigma }}_{\theta }=+\lambda \frac{R^2}{r^2}{\sigma }_v$$

(10)

Therefore, at a certain point around the tunnel, the stress change due to the advancing excavation of the tunnel can be obtained, as shown in the following equations:

$${\sigma }_r={\sigma }_r^i+∆{\sigma }_r=\left(1-\lambda \frac{R^2}{r^2}\right){\sigma }_v$$

(11)

$${\sigma }_{\theta }={\sigma }_{\theta }^i+∆{\sigma }_{\theta }=\left(1+\lambda \frac{R^2}{r^2}\right){\sigma }_v$$

(12)

In addition, from the results obtained by Kirsch’s solution, the radial displacement (u_r) in the elastic region can be given as the following:

$$\frac{2G}{{\sigma }_v}\frac{u_r}{R}=\mathsf{\lambda }\left(\frac{R}{r}\right)$$

(13)

where G is the shear modulus of the rock mass.

It is worth noting here that these equations must satisfy the boundary conditions and can be verified by the following instructions: (1) if λ = 0, it means that the tunnel has not been excavated, and the stress state of its surrounding rock is the same as the initial in-situ stress in the far-field; (2) if λ = 1, it indicates that the unsupported tunnel has been completely excavated, and the stress of its surrounding rock is the same as Kirsch’s stress in the near-field.

3. Derivation of Stress/Displacement Equations in the Plastic Regions

3.1. Derivation of the Confinement Loss in Elastic Limit with the Non-Linear Failure Criteria

Assuming that the mechanical behavior of the rock mass is elastic-perfectly-plastic, when the stress of the surrounding rock of the tunnel meets the limit conditions defined by the failure criterion, a disturbance zone, or the so-called plastic radius, gradually appears around the tunnel. It is not only a function of the peak strength parameter of the surrounding rock but also a function of the confinement loss in the elastic limit situation. When the value of confinement loss increases, the plastic radius increases, but the radial stress in the plastic zone decreases. Therefore, the radial stress is a function of the plastic radius.

The generalized Hoek–Brown failure criterion [34] is a widely used non-linear criterion in geotechnical engineering to assess the rock mass strength and predict its failure behavior. This failure criterion considers the influence of several parameters to capture the complex behavior of rock masses. These parameters include the uniaxial compressive strength (UCS), the geological strength index (GSI), and other geotechnical parameters. The UCS is the strength of the intact rock material, and the GSI represents the quality and condition of the rock mass. This criterion is used in various geotechnical applications, including tunneling, slope stability analysis, and rock engineering design. It provides a more realistic representation of the strength and failure behavior of rock masses compared to linear failure criteria such as the Mohr–Coulomb criterion shown in Figure 4. The mathematical expression for the generalized Hoek–Brown failure criterion is as follows,

$$f\left({\sigma }_3,{\sigma }_1\right)={\sigma }_1-{\sigma }_3-{\sigma }_{ci}{\left(m_b\frac{{\sigma }_3}{{\sigma }_{ci}}+s\right)}^a$$

(14)

where the empirical constants m_b and s are determined based on laboratory testing, field observations, and experience, and σ_ci is the uniaxial compression strength (UCS) of the intact rock. These constants vary depending on the rock type and the condition of the rock mass. The constants can be obtained from published charts or tables specific to the rock types and rock mass classifications. This equation describes the relationship between the stress state and the strength of the rock mass, and the parameters can be expressed as:

$$m_b=m_{i}exp\left(\frac{GSI-100}{28-14D}\right)$$

(15)

$$\mathrm{s}=exp\left(\frac{GSI-100}{9-3D}\right)$$

(16)

where m_b is a reduced value of the material constant. m_i, s, and a are constants for the rock mass. D is a factor that depends upon the degree of disturbance to which the rock mass has been subjected by blast damage and stress relaxation. It varies from 0 for undisturbed, in situ rock masses to 1 for very disturbed rock masses. Since the parameter a is variable and must be obtained from the experimental data by statistical regression analysis, in order to facilitate the calculation and derivation of the closed-form analytical solution of the differential equation, it is recommended to use a = 0.5 in this study. Then the representation of the equation can be given as follows:

$$f\left({\sigma }_3,{\sigma }_1\right)={{\sigma }_1-\sigma }_3-\sqrt{m_b{\sigma }_3{\sigma }_{ci}+s{\sigma }_{ci}^2}$$

(17)

In addition, the above equation can be normalized by the vertical stress (σ_v) for dimensionless, therefore:

$$f\left(\frac{{\sigma }_3}{{\sigma }_v},\frac{{\sigma }_1}{{\sigma }_v}\right)=\frac{{\sigma }_1}{{\sigma }_v}-\frac{{\sigma }_3}{{\sigma }_v}-\sqrt{m_b\left(\frac{{\sigma }_3}{{\sigma }_v}\right)\left(\frac{{\sigma }_{ci}}{{\sigma }_v}\right)+s\left(\frac{{\sigma }_{ci}^2}{{\sigma }_v^2}\right)}$$

(18)

Moreover, with consideration of the effect of stability, this equation can be redefined by the following:

$$f\left(\frac{{\sigma }_3}{{\sigma }_v},\frac{{\sigma }_1}{{\sigma }_v}\right)=\frac{{\sigma }_1}{{\sigma }_v}-\frac{{\sigma }_3}{{\sigma }_v}-\sqrt{2N{m}_b\left(\frac{{\sigma }_3}{{\sigma }_v}\right)+4s{N}^2}$$

(19)

where the stability number N can be defined and given as:

$$N=\frac{{\sigma }_{ci}}{{2\sigma }_v}$$

(20)

In the previous description of the stresses at the intrados of the tunnel ($r=R$), one can obtain the radial and tangential stresses as follows:

$${{\sigma }_3=\sigma }_r=\left(1-{\lambda }_e\right){\sigma }_v$$

(21)

$${{\sigma }_1=\sigma }_{\theta }=\left(1+{\lambda }_e\right){\sigma }_v$$

(22)

Substituting Equations (23) and (24) into the Hoek–Brown strength criterion Equation (21), the relationship can be shown as follows:

$$f\left(\frac{{\sigma }_3}{{\sigma }_v},\frac{{\sigma }_1}{{\sigma }_v}\right)=2{\lambda }_e-\sqrt{2N{m}_b\left(1-{\lambda }_e\right)+4s{N}^2}$$

(23)

Finally, the confinement loss in the elastic limit situation (λ_e) can be obtained as follows:

$${\lambda }_e=\frac{N}{4}\left(\sqrt{m_b^2+8\frac{m_b}{N}+16s}-m_b\right)$$

(24)

It is worth noting that the confinement loss in the elastic limit situation (λ_e) is a function of the initial vertical stress (σ_v) and the peak strength parameters of the rock mass (σ_ci, m_b, s).

3.2. Derivation of the Plastic Radius

To study the calculation of the plastic radius in the plastic zone, this paper adopts the assumption of axisymmetry, and its differential equation in the equilibrium state can be expressed as:

$$\frac{d{\sigma }_r}{dr}+\frac{{\sigma }_r{-\sigma }_{\theta }}{r}=0$$

(25)

By substituting the Hoek–Brown failure criterion, Equation (19) into this Equation (25), one can obtain the relationship as follows:

$${\int }_{{\sigma }_r^{-}}^{{\sigma }_r^{+}}\frac{d{\sigma }_r}{\sqrt{2N{m}_b{\sigma }_r+4s{N}^2}}={\int }_R^{R_p}\frac{dr}{r}$$

(26)

By integrating the above equation, the plastic radius around a tunnel can be obtained as follows:

$$\left(\frac{R_p}{R}\right)=exp\left\{\frac{1}{N{m}_b}\left[\left(2{\lambda }_e\right)-\sqrt{2N{m}_b\left(1-\lambda \right)+4s{N}^2}\right]\right\}$$

(27)

where R_p is the radius of the elastic-plastic interface, or so-called plastic radius, and is also a function of the initial vertical stress (σ_v) and the peak strength parameters of the rock mass (σ_ci, m_b, s). In particular, it is also dependent on the confinement loss (λ) which is an important factor in the incremental procedure to simulate the effect of advancing excavation of the tunnel face in the convergence–confinement method. In addition, the condition of the elastic-plastic interface needs to be satisfied, and it is verified by setting λ = λ_e, therefore R_p/R = 1. At the same time, the radial and tangential stresses are also satisfied by the Hoek-Brown non-linear failure criterion.

3.3. Derivation of Stress in the Plastic Region

To study the radial stress in the plastic region by substituting the Hoek-Brown failure criterion, Equation (19), into Equation (25), it can be expressed as follows:

$${\int }_{{\sigma }_r^{+}}^{{\sigma }_r^{-}}\frac{d{\sigma }_r}{\sqrt{2N{m}_b{\sigma }_r+4s{N}^2}}={\int }_{R_p}^r\frac{dr}{r}$$

(28)

Thus, the radial stress in the plastic region can be obtained as follows:

$$\frac{{\sigma }_r}{{\sigma }_v}=\left(1-{\lambda }_e\right)+\left(2{\lambda }_e\right)ln\left(\frac{r}{R_p}\right)+\frac{N{m}_b}{2}{ln}^2\left(\frac{r}{R_p}\right)$$

(29)

In addition, the tangential stress can also be given as:

$$\frac{{\sigma }_{\theta }}{{\sigma }_v}=\left(1+{\lambda }_e\right)+\left(2{\lambda }_e+N{m}_b\right)ln\left(\frac{r}{R_p}\right)+\frac{N{m}_b}{2}{ln}^2\left(\frac{r}{R_p}\right)$$

(30)

Therefore, the above stress can also be expressed as a function of σ = f(r, θ, λ) or σ = f(r, θ, Z). In addition, the condition of the elastic-plastic interface needs to be satisfied, and if it is verified by setting r = R_p, then the above equations become Equations (20) and (21).

3.4. Derivation of Displacement in the Plastic Region

In the elastic region, the tangential and radial strains, ${\epsilon }_{\theta }^e$ and ${\epsilon }_r^e$, can be expressed by the equations of constitutive law as:

$$\left\{\begin{array}{c}{\epsilon }_{\theta }^e\\ {\epsilon }_r^e\end{array}\right\}=\frac{1}{2G}\left[\begin{array}{cc}1-\nu & \nu \\ \nu & 1-\nu \end{array}\right]\left\{\begin{array}{c}{\sigma }_{\theta }\\ {\sigma }_r\end{array}\right\}$$

(31)

where ν is the Poisson’s ratio and G is the shear modulus of the rock mass.

Regarding the assumption of the small strain problem, the relationship between displacement and strain at any point in the surrounding rock of the tunnel can be expressed by a compatibility equation as follows:

$${\epsilon }_r=\frac{\partial u_r}{\partial r}$$

(32)

$${\epsilon }_{\theta }=\frac{\partial u_{\theta }}{r\partial \theta }+\frac{u_r}{r}\cong \frac{u_r}{r}$$

(33)

In this study, the non-associated plastic flow rule and the normality rule are used to determine the flow amount and direction of the displacement in the plastic region, as shown in Figure 5. Via the assumption that the plane strain and the elastic strain are relatively smaller than the plastic strain, the plastic parts of the radial and tangential strains can therefore be expressed as follows:

$${\epsilon }_r^p+{K_{\psi }\epsilon }_{\theta }^p=0$$

(34)

where the coefficient of passive pressure, K_p, can be expressed as:

$$K_{\psi }={tan}^2\left(45+\frac{\psi }{2}\right)$$

(35)

where $\psi$ is the dilation angle of the intact rock. The above equation in another form can be expressed as the following,

$$\left({\epsilon }_r-{\epsilon }_r^e\right)+K_{\psi }\left({\epsilon }_{\theta }-{\epsilon }_{\theta }^e\right)=0$$

(36)

where the subscripts p and e represent the plastic and elastic parts, respectively. Using Equations (32), (33), and (36) leads to the following differential equation:

$$f\left(r\right)={\epsilon }_r^e+{K_{\psi }\epsilon }_{\theta }^e=\frac{\partial u_r}{\partial r}+K_{\psi }\frac{u_r}{r}$$

(37)

By substituting the equation of constitutive law, Equation (31), into the above equation, then:

$$f\left(r\right)=\frac{\partial u_r}{\partial r}+K_{\psi }\frac{u_r}{r}=\frac{1}{2G}\left\{\left[\left(1-\upsilon \right)K_{\psi }-\nu \right]\left({\sigma }_{\theta }\right)+\left[\left(1-\upsilon \right)-\nu K_{\psi }\right]\left({\sigma }_r\right)\right\}$$

(38)

To obtain the analytical solution of the radial displacement (u_r) in the plastic region, the homogeneous solution and the particular solution can be obtained by using the engineering mathematics method in consideration of the boundary conditions. This differential equation can be expressed as:

(1)

Homogeneous solution ($f\left(r\right)=0$):

$$f\left(r\right)=\frac{\partial u_r}{\partial r}+K_{\psi }\frac{u_r}{r}=0$$

(39)

Integrating the above equation with the boundary conditions (r = R_p and $u_r=u_{R_p}$), it can be expressed as:

$${\int }_{u_r}^{u_{R_p}}\frac{\partial u_r}{\partial r}=-K_{\psi }{\int }_r^{R_p}\frac{\partial r}{r}$$

(40)

Therefore, the homogeneous solution of the radial displacement can be expressed as:

$$\frac{u_r}{r}=\frac{1}{r^{K_{\psi }+1}}{u}_{R_p}{\left(R_p\right)}^{K_{\psi }}$$

(41)

where $u_{R_p}$ is the radial displacement at the elastic-plastic interface (r = R_p) and can be given as:

$$u_{R_p}={\lambda }_e\left(\frac{{\sigma }_v}{2G}\right)\left(R_p\right)$$

(42)

To consider it in its normalized form, it can be obtained as:

$$\frac{2G}{{\sigma }_v}\frac{u_r}{R}={\lambda }_e\left(\frac{r}{R}\right){\left(\frac{R_p}{r}\right)}^{K_{\psi }+1}$$

(43)

(2)

Particular solution ($f\left(r\right)=\frac{\partial u_r}{\partial r}+K_{\psi }\frac{u_r}{r}$):

According to the stresses obtained in the plastic region, Equations (29) and (30), can be written in a new form as:

$$f\left(r\right)=\frac{\partial u_r}{\partial r}+K_{\psi }\frac{u_r}{r}=\frac{1}{2G}\left[\left(1-v-\nu K_{\psi }\right)\left({\sigma }_r\right)+\left(K_{\psi }-v{K}_{\psi }-\nu \right)\left({\sigma }_{\theta }\right)\right]$$

(44)

By rearranging the above equation, it can be represented as follows:

$$f\left(r\right)=\frac{1}{2G}\left[D_1+D_{2}ln\left(\frac{r}{R_p}\right)+D_3{ln}^2\left(\frac{r}{R_p}\right)\right]$$

(45)

where the coefficients D₁, D₂, and D₃ can be given as:

$$D_1=\left(K_{\psi }-1\right)\left({\lambda }_e\right)$$

(46)

$$D_2=\left(1+K_{\psi }\right)\left(1-2v\right)\left(2{\lambda }_e\right)+\left(K_{\psi }-v{K}_{\psi }-\nu \right)N{m}_b$$

(47)

$$D_3=\left(1+K_{\psi }\right)\left(1-2v\right)\left(\frac{N{m}_b}{2}\right)$$

(48)

By combining the homogeneous solution Equation (41) and the particular solution Equation (45), the radial displacement in the plastic region can be expressed as:

$$\frac{u}{r}=\frac{1}{r^{K_{\psi }+1}}{\int }_{R_p}^r{r}^{K_{\psi }}f\left(r\right)dr+\frac{u_r}{r}$$

(49)

In addition, the radial displacement in the plastic region can be obtained as:

$$\frac{u}{r}=\frac{1}{2G}\frac{1}{r^{K_{\psi }+1}}\left[D_1{f}_1\left(r\right)+D_2{f}_2\left(r\right)+D_3{f}_3\left(r\right)+2G{u}_{R_p}{\left(R_p\right)}^{K_{\psi }}-D_1{f}_1\left(R_p\right)-D_2{f}_2\left(R_p\right)-D_3{f}_3\left(R_p\right)\right]$$

(50)

where $f_1\left(r\right)$, $f_2\left(r\right)$ and $f_3\left(r\right)$ can be given as:

$$f_1\left(r\right)=\int r^{K_{\psi }}dr=\frac{r^{K_{\psi }+1}}{K_{\psi }+1}$$

(51)

$$f_2\left(r\right)=\int r^{K_{\psi }}ln\left(\frac{r}{R_p}\right)dr=\frac{r^{K_{\psi }+1}}{K_{\psi }+1}\left[ln\left(\frac{r}{R_p}\right)-\frac{1}{K_{\psi }+1}\right]$$

(52)

$$f_3\left(r\right)=\int r^{K_{\psi }}{ln}^2\left(\frac{r}{R_p}\right)dr=\frac{r^{K_{\psi }+1}}{K_{\psi }+1}\left[{ln}^2\left(\frac{r}{R_p}\right)-\frac{2}{K_{\psi }+1}ln\left(\frac{r}{R_p}\right)+\frac{2}{{\left(K_{\psi }+1\right)}^2}\right]$$

(53)

From the above series of mathematical operations, the closed-form analytical solution of the radial displacement in the plastic zone is finally obtained. For validation of this equation, one can examine the radial displacement at the elastic-plastic interface (r = R_p), and the same result is obtained as shown by Equation (42).

4. Implementation of Incremental Procedure for the Analytical Solution

Calculation Steps and Analysis Process of the Incremental Method

The calculation steps and analysis process of the incremental method used in this study can be called an explicit calculation program. It can realize the calculation of the analytical solution obtained above and convert it into an executable calculation program, such as using a simple calculation spreadsheet to perform direct calculation, and the results can be shown in Figure 6. Therefore, the procedure proposed in this paper can be called an explicit analysis method (EAM), which can consider the confinement loss as an incremental step, thereby simulating the effect of tunnel advancing excavation, calculating the stress/displacement in each step, and drawing the ground reaction curve, the stress path of the intrados of the tunnel, and the stress/displacement distribution on the cross-section of the tunnel. The results can be seen in Figure 7.

The calculation steps and analysis process of the explicit analysis method (EAM) can be interpreted one by one through the following analysis steps:

(1)

Data input stage: This stage is the input action of calculation data, including data such as tunnel excavation radius, unsupported distance, overburden depth, physical and mechanical parameters of surrounding rock, etc.

(2)

Calculation stage of confinement loss: To estimate the confinement loss λ_z at a certain distance z from the tunnel working face, it can be calculated by Equation (1).

(3)

Calculation stage of the incremental step: Dividing the confinement loss by n segments yields the incremental step, which can be expressed as follows:

$$∆\lambda =\frac{1}{n}\left(1-{\lambda }_z\right)$$

(54)

(4)

Calculating each step value of λ as

$$\begin{array}{cc}{\lambda }_{i+1}={\lambda }_z& i=0\\ {\lambda }_{i+1}={\lambda }_i+i\mathsf{\Delta }\lambda & i=1~\left(n-1\right)\end{array}$$

(55)

(5)

Attaining the final value ${\lambda }_n={\lambda }_{i+1}$

(6)

According to Equation (24), it estimates the confinement loss in the elastic limit situation (λ_e).

(7)

Calculation stage of stress/displacement in the elastic zone: If ${\lambda }_{i+1}<{\lambda }_e$, it indicates that the stress state is in the elastic zone, and the radial stress, tangential stress, and radial displacement can be calculated according to Equations (11), (12), and (13), respectively.

(8)

Calculation stage of stress/displacement in the plastic zone: If ${\lambda }_{i+1}\ge {\lambda }_e$, it indicates that the stress state is in the plastic zone, and the plastic radius can be calculated according to Equation (27). Once this value is obtained, the program will automatically substitute Equations (29), (30), and (50) to calculate the radial stress, tangential stress, and radial displacement, respectively.

(9)

Data recording stage: In this stage, the stress/displacement are recorded respectively, which can be expressed as the representation of the distribution of stresses/displacements ($\frac{r}{R},\frac{{\sigma }_r}{{\sigma }_v}$), ($\frac{r}{R},\frac{{\sigma }_{\theta }}{{\sigma }_v}$), and ($\frac{r}{R},\frac{{2Gu}_r}{R{\sigma }_v}$) on the cross-section of the tunnel and ($\frac{{2Gu}_r}{R{\sigma }_v},\frac{{\sigma }_r}{{\sigma }_v}$) at the intrados of the tunnel.

(10)

Judgment stage 1 of the calculation ends: When i < n − 1, repeat steps (4) through (10).

(11)

Judgment stage 2 of the calculation ends: When i = n − 1, stop performing calculations and record all the obtained data.

(12)

Drawing production stage: Drawing the distribution diagram of stress/displacement on the cross-section and at the intrados of the tunnel.

5. Comparison of Results Obtained between Published Data and This Study

Case studies comparing this study with the EAM calculations include published research such as Sharan (2003) [44] and Rocksupport (2004) [16]. First, to analyze the ground reaction caused by tunnel-advancing excavation under unsupported conditions, the input data required for numerical calculation are shown in

Table 1. Input data of the computation of EAM by using the published data [16,44].

Reference	Sharan (2003) [44]		Rocksupport (2004) [16]
Parameter	Case I	Case II	Case III	Case IV
E (MPa)	40,000	5500	3530	2100
ν	0.20	0.25	0.3	0.3
mi	7.5	7.5	10	12
GIS	100	80	22	17
D	0	0	0	0
σci (MPa)	300	30	5	4
Kψ	0	0	0	0
σv (MPa)	108	30	1.62	2.02
R (m)	4.0	5.0	6.0	5.0

. To investigate the effect of parameters used by the different case studies on the stress/displacement, the dimensionless radial and tangential stresses and the dimensionless radial displacement in the elastic region and the plastic region are compared using the typical cases of hard rock (Case I) and soft rock (Cases II–IV) [16,44]. Table 1 shows the properties of the rocks, and it can be observed that the dilation angle (K_ψ) and the degree of disturbance (D) of the rocks are assumed to be zero in this study. In addition, the values of GSI and UCS of the hard rock used in the calculation are much larger than those of the soft rock.

5.1. Stress/Displacement at the Intrados of the Tunnel

Due to the advancing excavation of the tunnel face, the mechanical behavior of the surrounding rock is altered. This behavior can be explained by the ground reaction curve and stress path at any point on the interior of the tunnel. For the presentation of stress/displacement in the elastic and plastic regions, the radial displacement is expressed in a dimensionless (2Gu_r/Rσ_v) manner, while the radial and tangential stresses are expressed in a normalized vertical stress (σ_v) manner. The obtained results are shown in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15.

As shown in Figure 8, the results obtained by the EAM using the data from Case I [44] indicate that the behavior of the hard rock presents linear elasticity in the ground reaction curve and the stress path. The calculated result is in agreement with the theoretical value that the dimensionless radial displacement is equal to 1.0 (Figure 8a), and the dimensionless radial and tangential stresses are equal to 0 and 2, respectively (Figure 8b).

In addition, the characteristic behavior of the soft rocks (Cases II–IV) is described in Figure 9, Figure 10 and Figure 11. It indicates that the behavior of the soft rocks presents non-linear elastic-perfectly-plasticity in the ground reaction curve. The interface between the elastic region (line AB) and the plastic region (curve BC) is represented by point B. For the stress path, the Hoek-Brown failure criterion is obviously shown to have characteristic non-linearity by the curve BG in Figure 9b, Figure 10b and Figure 11b.

According to the analysis results between EAM and the listed articles (shown in

Table 2. Comparison of results between EAM and published studies.

Published Studies	Radial Displacement, uR (mm)	Plastic Zone Radius, Rp (m)	EAM Radial Displacement, uR (mm) (Error * %)	EAM Plastic Zone Radius, Rp (m) (Error * %)
Case I	12.52	N/A	12.96 (3.51%)	N/A
Case II	56.05	8.12	57.30 (2.23%)	8.25 (1.67%)
Case III	12.0	13.77	12.06 (0.5%)	15.55 (12.93%)
Case IV	65.5	26.30	66.80 (1.98%)	26.9 (2.28%)

* Error (%) = |(published data − EAM)/published data| × 100%.

), the percentage error of the radial displacement is from 0.5% to 3.51% for the hard and soft rocks. The results obtained by the calculation of EAM in this study show a consistent trend with the data in the references.

5.2. Distribution of Stress/Displacement on the Cross-Section of the Tunnel

Under the isotropic stress condition, the distribution of stress/displacement on the cross-section of the tunnel caused by the excavation of the circular tunnel is investigated in this study. According to the stress distribution results of the tunnel, the following five stress states can be used to describe and verify the stress changes caused by the advancing excavation of the tunnel face: (1) When λ = 0 ($r\to \infty$), it indicates that the stresses are in the initial stress state as shown by the horizontal line on the right side in Figure 12a; (2) When 0 $\le \lambda \le {\lambda }_e$ ($R_p\le r\le \infty$), it indicates that the stress state is in the elastic region, and at this time the tangential stress and radial stress increase with the increase of the confinement loss, and the two stresses are symmetrically separated up and down along the horizontal axis (r/R); (3) When $\lambda ={\lambda }_e$ ($r=R_p$), the plastic radius appears, and two stresses are at the elastic-plastic interface. The radial stress begins to change the curvature, and the tangential stress attains the maximum value as shown in Figure 13a, Figure 14a and Figure 15a; (4) When ${\lambda }_e\le \lambda \le 1$ ($R\le r\le R_p$), it indicates that the stresses are in the plastic region, and this leads to both the tangential stress and radial stress being decreased steeply; and (5) When $\lambda =1$ ($r=R$), the radial stress becomes zero and the tangential stress is equal to the coefficient $\sqrt{4s{N}^2}$ proposed by the Hoek–Brown failure criterion.

According to results obtained by the explicit analysis method (EAM) and the analysis results provided by the listed literature (as shown in Table 2), the percentage error of the plastic zone radius is from 1.67% to 2.28% for the soft rocks. However, the plastic radius obtained in Case III shows an error value of up to 12%, which may occur during the recursive process of the plastic radius solution.

6. Conclusions

Through the nonlinear mechanical behavior (Hoek–Brown model) caused by the advancing excavation of a circular tunnel, this study conducts rigorous theoretical analysis (CCM), derivation of mechanical differential equations, establishment and practice of incremental calculation programs (EAM), as well as a comparison with published literature results. The following conclusions can be drawn: (1) The closed-form analytical solution for the behavior of ground based on the Hoek–Brown nonlinear failure criterion due to the advancing excavation of a circular tunnel in an isotropic stress state is investigated. (2) The method of judging whether the stress state is in the elastic region or the plastic region depends on the calculation of the confinement loss in the elastic limit situation, which is affected by the peak strength parameters of the rock mass and the initial vertical stress. (3) It is found that when the confinement loss increases, the plastic radius also increases, resulting in a gradual decrease in the radial stress and tangential stress while the radial displacement increases. The plastic radius is also affected by factors such as confinement losses at the elastic limit and rock mass peak strength parameters. (4) In particular, this study proposes an incremental procedure of the explicit analysis method (EAM), which takes the confinement loss as an incremental step to simulate the effect of advancing excavation on the tunnel face, calculates the stress and displacement at each step, and plots the ground reaction curves, stress paths at the intrados of the tunnel, and the distributions of stress/displacement on the cross-section of the tunnel. (5) The comparison results between the analytical solution of the explicit program and the published literature show a consistent trend in elastoplastic media. (6) The stress change caused by the advancing excavation of the tunnel face can be explained by the stress gradient, which is the difference between the far-field stress and the near-field stress around the tunnel. The increase in stress can be accounted for as part of the stress gradient using confinement losses. (7) The proposed incremental procedure can deal with the influence of the nonlinear failure criterion of rock mass and can not only be a useful tool for the analysis of circular tunnels in an isotropic stress state but may also be applied to the simulation of the ground behavior of the tunnel under anisotropic stress conditions in the next stage of research. However, this method also has its limitations; namely, the CCM analysis can only be used for the excavation of circular deep tunnels. As for the analysis of non-circular shallow tunnels, it cannot meet the assumptions of this theory.