Analytical Solution Implemented by an Incremental Procedure for the Ground Reaction Based on Hoek–Brown Failure Criterion in the Tunnel Convergence-Confinement Method

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]. 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 [3,4]. 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 [5,6,7,8,9].

The CCM focuses on achieving three main objectives (shown in Figure 1): (1) 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 [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 the convergence and enhancing tunnel stability [16,17,18,19,20,21,22,23]. (3) 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 confinement loss can be obtained by the relationship between convergence measurement and unsupported distance to the tunnel face. It could be a positive scalar, an incremental value, or a function for considering the effect of advancing excavation of the tunnel face [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 ensure the stability and safety of underground excavations [29,30].

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 characteristic of the GRC is investigated and implemented by using the non-linear failure criterion of generalized Hoek-Brown [31,32,33]. 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 [34,35,36,37]. 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 [38,39,40]. 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 objective of this paper is, therefore, to present a comprehensive discussion on the non-linear failure criterion of generalized Hoek-Brown on a ground reaction in the convergence-confinement method for the tunnel support design, providing the full set of equations and the details of their derivation, verification, comparison, together with illustrative diagrams of the distribution of stresses/displacements at the intrados and on the cross-sections of the tunnel.

2. Problem description

2.1. Relationship between confinement loss and tunnel advancing excavation effect

As a result of the advancing excavation of the tunnel face, a stress field change occurs in the surrounding rock, which in turn produces the results of stress redistribution and convergence of the excavation face. In other words, the effect of advancing excavation of the tunnel face is the consideration of a change of stresses and displacements caused by tunnel sequential excavation at certain measurement points of an observation section. Applying the tunnel convergence measurement, establishing the functions concerning the effect of advancing excavation of the tunnel face, and providing the procedure to estimate the confinement loss are mainly proposed in this study. In addition, the function of the confinement loss curve (CLC) can be shown in Figure 2, and at a certain distance (z) to the working face, the definition of confinement loss (λ_z) can be given as [5,6]

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

(1)

Concerning the assumption of the function concerning the effect of advancing excavation of the tunnel face, f(z/R), which can also be called the advancing effect function, can be proposed 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

The convergence-confinement method is based on the analysis of the stress and strain state that developed in the rock mass around a tunnel. The stresses at the periphery of the tunnel are altered due to the excavation of the tunnel. Therefore, the variation of stresses can be mainly interpreted by the stress gradient which is a difference between the far-field stress and the near-field stress around a tunnel. By using the hypothesis of increments in the numerical analysis, the increment of stresses could take into account the confinement loss λ as a fraction of the stress gradient. This is an important principle concept used by the convergence-confinement method. The variation of stresses can be described by the following representations [5,6]:

(1) Stresses in the far-field (the initial in-situ stresses, r –› ∞): Before opening the tunnel, the stresses in the rock mass are in equilibrium and are the same as the lithostatic stresses and constant in the surrounding medium of the tunnel. Hence, the initial anisotropic stress in the rock mass can be shown as shown in Figure 3, and indicated by the following equation,

$${\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 in the polar coordinates (r, θ ), respectively.

(2) Stresses in the near-field (the final stresses or the boundary stresses at the intrados of the tunnel, R ≤ r < ∞ ): After opening the tunnel, the stresses at any point around the tunnel in the rock mass are altered and can be obtained by the Kirsch solutions as follows,

$${\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 }_r^f$ and ${\sigma }_{\theta }^f$ are the final radial stress and the final tangential stress around the tunnel proximity, respectively.

Assuming that the stress gradient due to the advancing excavation of the tunnel face is a difference between the far-field stress and the near-field stress, a difference between the radial stress and the tangential stress can be expressed by the following equations,

$${\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)

Using the hypothesis of increments in the numerical analysis, the increment of stresses can be expressed by the multiplication of the stress gradient with the confinement loss λ which is a fraction of the stress gradient, and can be obtained as the followings,

$$∆{\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)

Finally, the variation of stresses due to the advancing excavation of the tunnel (0 ≤ λ ≤ 1) at a certain point of tunnel proximity (r, θ ) can be represented by the followings,

$${\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 the Kirsch 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 to be aware that these equations must be satisfied by the boundary conditions proposed by the convergence-confinement method, and also can be verified by the following interpretations: (1) When λ = 0, it means that the tunnel is not excavated yet, and stresses in the rock mass are the same as the initial in-situ stresses ${\sigma }_r^i$ and ${\sigma }_{\theta }^i$ (or the far-field stress, e.g. r –› ∞), and (2) When λ = 1, it means that the tunnel is already excavated, and stresses in the rock mass are the same as Kirsch’s stresses ${\sigma }_r^f$ and ${\sigma }_{\theta }^f$ (or the near-field stress, e.g. r = R).

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

In the more general case of rock mass with elastic-perfectly-plastic behavior, the confinement loss in the elastic limit situation (λ_e) appears when the elastic limit of the rock mass is reached, for a certain radial stress (σ_r) and with the decreasing of r, and when the stress state reaches the limit conditions defined by the strength criterion. Such a value of r is called the plastic radius (R_p) which is not only a function of the peak strength parameters of the rock mass but also a function of the confinement loss in the elastic limit situation (λ_e). When the confinement loss continuously increases, the plastic radius increases, but the radial stress decreases in the plastic region; therefore this radial stress is obviously a function of the plastic radius. For this reason, the radial stress in the plastic region can be obtained by the following representations.

The generalized Hoek-Brown failure criterion [33] 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, 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$$

(16)

where σ_ci is the uniaxial compression strength (UCS) of the intact rock, and the empirical constants m_b and s are determined based on laboratory testing, field observations, and experience. 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)$$

(17)

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

(18)

where m_b is a reduced value of the material constant. m_i, and 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. In this study, it’s supposed that a = 0.5, 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}$$

(19)

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)}$$

(20)

Moreover, with the 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}$$

(21)

where the stability number N can be defined and given as

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

(22)

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

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

(23)

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

(24)

Substituting equations (23) and (24) into the Hoek-Brown strength criterion equation (21), can be shown the relationship as the following,

$$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}$$

(25)

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

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

(26)

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

3.2. Derivation of the plastic radius

To study the plastic radius in the plastic region, the differential equation of equilibrium for the axisymmetric problem can be expressed as

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

(27)

By substituting the Hoek-Brown failure criterion, equation (21) into this equation (27), then one

$${\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}$$

(28)

can obtain the relationship as follows,

By integrating the above equation, and therefore the plastic radius around a tunnel can be obtained as

$$\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\}$$

(29)

where the R_p is the radius of the elastic-plastic interface or so-called plastic radius, and is also a function of the peak strength parameters of the rock mass (σ_ci , m_b , s ), and the initial vertical stress (σ_v). In particular, it is also dependent on the confinement loss (λ) which is an important key 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 (21), into equation (27), it can be expressed as

$${\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}$$

(30)

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

$$\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)$$

(31)

In addition, the tangential stress can be also 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)$$

(32)

It should be noted that the above-expressed stresses can be represented by the functions as σ =f(r, θ, λ) or σ =f(r, θ, Z). In addition, the condition of the elastic-plastic interface needs to be satisfied, and it is verified by setting r = R_p, then the above equations become the equations (22) and (23).

3.4. Derivation of displacement in the plastic region

In the elastic region, the radial and tangential strains, ${\epsilon }_r^e$ and ${\epsilon }_{\theta }^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\}$$

(33)

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

For the small strain problem, the relationship between strain and displacement at any point in the rock mass can be expressed by the equations of compatibility as

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

(34)

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

(35)

Additionally, to determine the displacement field in the plastic region, a plastic flow rule is needed. By assuming that the elastic strains are relatively small in comparison to the plastic strains and that a non-associated flow rule is valid as shown in Figure 5, the plastic parts of the radial and tangential strains may be related to the plane strain condition.

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

(36)

where the coefficient of passive pressure K_p can be expressed as

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

(37)

where $\psi$ is the dilation angle of the intact rock. The above equation with 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$$

(38)

where the subscripts e and p represent the elastic and plastic parts, respectively. Using Eq. (34), (35), and (38) 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}$$

(39)

By substituting the equation of constitutive law, equation (33), 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\}$$

(40)

This differential equation may be solved by engineering mathematics with the homogeneous solution and the particular solution, and by using the following boundary condition for the radial displacement (u_r) in the plastic region.

(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$$

(41)

By integrating the above equation with the boundary condition (r = R_p and $u_r=u_{R_p}$), it can be represented as

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

(42)

Finally, the radial displacement with the homogeneous solution can be expressed as

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

(43)

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)$$

(44)

To consider in the 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}$$

(45)

(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 (31) and (32), it 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]$$

(46)

By rearranging the above equation, therefore it can be represented as

$$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]$$

(47)

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

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

(48)

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

(49)

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

(50)

By combining the homogeneous solution equation (43) and the particular solution equation (47), 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}$$

(51)

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]$$

(52)

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}$$

(53)

$$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]$$

(54)

$$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]$$

(55)

Finally, the closed-form analytical solution for the radial displacement in the plastic region can be obtained. For validation of this equation, one can examine the radial displacement at the elastic-plastic interface (r = R_p), then the same result is obtained as shown by the equation (44).

4. Implementation of incremental procedure for the analytical solution

Calculation steps in the incremental procedure method

The incremental procedure used by this study is also called an explicit procedure which can realize and transform the analytical solutions into an executable computation that can be directly calculated by using a simple calculation spreadsheet (e.g., Microsoft Excel) as shown in Figure 6. The particular procedure proposed by this paper is called the explicit analysis method (EAM), which can deal with using the confinement loss as an incremental step to simulate the effect of advancing excavation of the tunnel, calculating the stresses/displacements in each step, and drawing the ground reaction curve, the stress path at the intrados of the tunnel, and the distribution of stresses/displacements on the cross-sections of the tunnel as shown in Figure 7.

The steps of calculation for the EAM are described by the following representations.

(1)

Input data: This relates to the data of the initial in-situ stress, geometry of the tunnel, material properties, and unsupported distance.

(2)

To estimate the confinement loss λ_z at a certain distance z from the tunnel working face, one can use a given value of λ_z as a chosen effect of the unsupported distance of tunnel excavation. Therefore, λ_z can be determined from the Eq. (1).

(3)

Dividing the confinement loss λ_z by n segments, the incremental step λ can be expressed as

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

(56)

(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}$$

(57)

(5)

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

(6)

According Eq. (26), estimates the confinement loss in the elastic limit situation ( e).

(7)

If ${\lambda }_{i+1}<{\lambda }_e$ , it means that the stress state is in the elastic region and that the radial and tangential stresses/displacements can be calculated with Eq. (11), (12), and (13)

(8)

If ${\lambda }_{i+1}\ge {\lambda }_e$, it means that the stress state is in the plastic region and that the plastic radius Rp can be calculated with Eq. (29). Once one obtained this value Rp, the procedure automatically substitutes into Eq. (31), Eq. (32), and Eq. (52) for the radial and tangential stresses, and the radial displacement, respectively.

(9)

Recording the calculated data, relates to 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)

When i < n-1, repeat steps (4) through (10).

(11)

When i = n-1, the process is not repeated, and the data from each step is recorded.

(12)

Drawing the distribution of stresses/displacements at the intrados and on the cross-section 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) [41] and Rocksupport (2004) [42]. 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. 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 by using the typical cases of hard rock (Case I) and soft rock (Case II~IV) [41,42]. Table 1 show 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

The behavior of rock mass at the intrados of the tunnel is altered by the advancing excavation of the tunnel face and may be interpreted by the ground reaction curve and the stress path at any point of the intrados of the tunnel. In this paper, the dimensionless radial displacement (2Gu_r/Rσ_v), and radial and tangential stresses normalized by the vertical stress (σ_v), in the elastic region and plastic region are plotted in Fig 8-15.

As shown in Figure 8, the results obtained by the EAM using the data of Case I [41] 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 8(a)), and the dimensionless radial and tangential stresses are equal to 0 and 2, respectively (Figure 8(b)).

In addition, the characteristic behavior of the soft rocks (Case II~IV) is described in Figure 9, Figure 10 and Figure 11. It indicates that the behavior of the soft rocks presents the 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 the characteristic non-linearity by the curve BG in Figures 9(b)–11(b).

According to the analysis results between EAM and the listed articles (shown in Table 2), 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

Around a circular tunnel excavation in an isotropic stress state, one is interested in the stress/displacement distribution over the tunnel cross-section. 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 the horizontal line on the right side in Figure 12(a); (2) When 0$\le \lambda \le {\lambda }_e$ ($R_p\le r\le \infty$), it describes that the stresses are in the elastic region, and the radial and tangential stresses increase with the increase of the confinement loss, and both stresses are symmetrically separated along the horizontal axis (r/R); (3) When $\lambda ={\lambda }_e$ ($r=R_p$), the plastic radius appears and 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 Figures 13a–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 radial stress and tangential stress being decreased steeply; and (5) Until $\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.

Figure 12. Results obtained by the EAM using the data of Case I [41], (a) stress distribution, and (b) radial displacement distribution on the cross-section of the tunnel (R_p/R = 1.0).

Figure 12. Results obtained by the EAM using the data of Case I [41], (a) stress distribution, and (b) radial displacement distribution on the cross-section of the tunnel (R_p/R = 1.0).

Figure 13. Results obtained by the EAM using the data of Case II [41], (a) stress distribution, and (b) radial displacement distribution on the cross-section of the tunnel (R_p/R = 1.65).

Figure 13. Results obtained by the EAM using the data of Case II [41], (a) stress distribution, and (b) radial displacement distribution on the cross-section of the tunnel (R_p/R = 1.65).

Figure 14. Results obtained by the EAM using the data of Case III [42], (a) stress distribution, and (b) radial displacement distribution on the cross-section of the tunnel (R_p/R = 2. 59).

Figure 14. Results obtained by the EAM using the data of Case III [42], (a) stress distribution, and (b) radial displacement distribution on the cross-section of the tunnel (R_p/R = 2. 59).

Figure 15. Results obtained by the EAM using the data of Case IV [42], (a) stress distribution, and (b) radial displacement distribution on the cross-section of the tunnel (R_p/R =5.38).

Figure 15. Results obtained by the EAM using the data of Case IV [42], (a) stress distribution, and (b) radial displacement distribution on the cross-section of the tunnel (R_p/R =5.38).

As shown in Table 2, according to the analysis results between EAM and the listed articles, 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 a large error value.

6. Conclusions

Through the rigorous theoretical analysis (CCM) of nonlinear mechanical behavior (Hoek-Brown model) caused by tunnel advancing excavation, the derivation of mechanical partial differential equations, the establishment and practice of incremental calculation program (EAM), and the comparison of published results, the following conclusions can be drawn:

(1)

A theoretically consistent closed-form analytical solution was derived for the elastic-perfectly-plastic analysis of a circular tunnel in a rock mass governed by the Hoek-Brown non-linear failure criterion and subject to isotropic in-situ stress.

(2)

The confinement loss in the elastic limit situation is a function of the peak strength parameters of the rock mass, and the initial vertical stress.

(3)

The plastic radius is also a function of the peak strength parameters of the rock mass, the confinement loss in the elastic limit situation, and is also dependent on the confinement loss. As the confinement loss increases, the plastic radius increases and consequently leads to both the radial stress and tangential stress being decreased, even as the radial displacement increases simultaneously in the plastic region.

(4)

The incremental procedure of the explicit analysis method (EAM) is particularly proposed to deal with using the confinement loss as an incremental step to simulate the effect of advancing excavation of the tunnel face, calculating the stresses and displacements in each step, and drawing the ground reaction curve, the stress path at the intrados of the tunnel, and the distribution of stresses/displacements on the cross-sections of the tunnel.

(5)

The agreement between the published results and the proposed closed-form solutions with the explicit procedure was found to be excellent in elastic-perfectly-plastic media.

(6)

The variation of stresses due to the advancing excavation of the tunnel face can be mainly interpreted by the stress gradient which is a difference between far-field stress and near-field stress around a tunnel and can be obtained employing the hypothesis of increments in the numerical analysis. The increment of stresses can be considered with the confinement loss as a fraction of the stress gradient.

(7)

The proposed incremental procedure method can deal with the effect of the non-linear failure criterion of the rock mass, not only it can be a useful tool for the analysis of a circular tunnel in an isotropic stress state, but it perhaps may be applied to the simulation of the behavior of the tunnel under the condition of an anisotropic stress in the next stage study.