# Unified Analytical Solutions of Circular Tunnel Excavated in an Elastic-Brittle-Plastic Rock Mass considering Blast-Induced Damage and Dead Weight Loading

Authors: Yue Cao; Liang Chen; Jinhai Xu; Chong Li; Wei Zhang
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

Stress and deformation around circular tunnel are crucial for optimizing the support system and evaluating the tunnel stability. The damage zone induced by blasting or mechanical excavation can dramatically influence the support design and methods because the self-weight of broken rock mass at the roof of the tunnel can exert a high pressure on the support system, leading to the support system instability due to the overload. This paper presents a new closed-form solution for analyzing the stress and deformation of deep circular tunnel excavated in elastic-brittle rock mass with the consideration of the rock gravity and damage zone by using the unified strength criterion. A new modified equilibrium equation in the fracture zone is used to determine the stress and the radius of fracture zone. The correctness of the solution is also verified by comparison with the numerical simulation results. The results illustrate that the rock gravity, damage zone radius, and intermediate principal stress have an extremely important influence on the ground response. The tunnel surface convergence and damage zone radius with the consideration of the gravity are obviously larger than those without consideration of the gravity. The rock gravity effect under the high intermediate principal stress gradually weakens, illustrating that the intermediate principal stress is beneficial to tunnel stability. Large deformation instability of the tunnel is dependent on the extension of damage zone. The larger the radius of damage zone, the larger both fracture range and tunnel surface deformation. The proposed solution in this study is novel and can be used to assess the ground convergence for different scenarios and to optimize the support system during the early design stage of the tunnel.


---

## Body

## 1. Introduction

During excavation of rock, the extent of damage zone varies with the tunnel shape, excavation method, and rock mass properties. The distribution of damage greatly impacts the ground response and support design. This damaged zone can be considered a rock failure zone in which the rock mass properties are seriously weakened due to the disturbance caused by the excavation. The extent of the damage zone near the tunnel surface in blasted tunnels is far more significant owing to the impact wave and stress redistribution [1, 2]. The damage can range from a few centimetres in tunnels with tunnel boring machine (TBM) to several decimetres and several meters in drilling and blasting [3, 4]. It is therefore necessary to adopt the blast-induced damage zone (BIDZ) in account for the response analysis of surrounding rock. It can be assumed that the damage to the material forms a cylindrical zone around the tunnel [2].Many researchers have mainly concentrated on determining the range of excavation- or blast-induced damage zones beyond the perimeter of tunnels excavated in strain-softened, elastic-brittle-plastic, and elastoplastic rock masses [2, 5–8]. Generally, three sets of criteria are used to analyze the stresses and deformation behaviour of tunnels: Hoek-Brown (H-B), Mohr-Coulomb (M-C), and generalized Hoek-Brown (GHB) criteria. A small number of true triaxial results from laboratory testing and numerical simulation have demonstrated that intermediate principal stress has a more significant effect on rock strength and can influence the rock failure behaviour after tunnel excavation [9–11]. However, the above three sets of criteria ignore this effect, leading to several limitations in engineering applications. Even though several researchers deem the Drucker-Prager criterion as the failure condition of rock strength when studying the response behaviour of a rock mass, they put the same amount of emphasis on it in studying the impact of both the minimum and maximum principal stress [12, 13]. Yu and He [14] have advanced a unified strength criterion (USC) that will be utilized to depict the reaction conduct of rock mass in this study. In recent years, some researchers have investigated the response behaviour of circular tunnel considering different rock materials and flow rules based on USC [15–19]. However, those researchers ignored the effect of rock gravity in the damage zone on stresses and deformation of the tunnel.The gravitational forces inside the damage zone have a significant influence on the response characteristic of the rock mass and existing support system. Furthermore, the relation between support pressure and tunnel convergence at the roof is unusual compared with that at the side wall because the gravity forces do not exist at the side wall. Several researchers have analyzed the impact of rock gravity on the ground reaction of circular tunnels excavated in M-C, H-B, or GHB rock mass [1, 2, 20–24]. They thought that the dead weight of the material in the broken zone can be represented by body forces acting radially towards the centre of the circular tunnel. This approximation was acceptable, since stability conditions for the roof are of primary concerns. Carranza-Torres and Fairhurst [21] improved the equilibrium equation of the failure zone by equating the dead weight of rock mass in the failure zone to body forces acting downward towards the centre of the opening at the tunnel roof and presented an analytical solution for rock material obeying Mohr-Coulomb criterion. The results showed that the gravity played an important role in roof instability, so it should be viewed in the equilibrium equation. Herein, the effect of rock material’s dead weight was considered in a similar manner to Carranza-Torres and Fairhurst [21] but provided an analytical solution of a deep tunnel that considers a plastic zone and BIDZ with the USC. However, the effect of the damage zone was ignored. Hedera and Weems [2] divided the surrounding rock of the tunnel into elastic, plastic undamaged, and plastic damaged zones and advanced a closed solution of opening which considers the effect of rock gravity with the generalized Hoek-Brown criterion. However, the previous study did not consider the effect of gravitational load in rock masses that follow a uniform strength criterion.Both intermediate principal stress and gravitational force exert considerable influence on the equilibrium and should be evaluated using the governing equations for accurately predicting the stresses and deformation around the tunnel. Combining the impact of the damage zone, rock gravity, and intermediate principal stress is crucial in the estimation of roof convergence and optimization of the support system but has not been investigated yet. In this paper, the modified equilibrium equation is used in the fracture zone in conjunction with consideration of the gravity effect and accounts for the intermediate principal stress and the damage zone. This analytical solution is beneficial for a more direct understanding of the response behaviour of the damaged zone.


## 2. Unified Strength Criterion (USC)

The USC takes into account the effect of all the stress components acting on an orthogonal octahedron based on a twin-shear element model [18] and can be expressed in terms of the following principal shear stress parameters [14, 15]:(1a)τ13+bτ12+βσ13+bσ12=Cτ12+βσ12≥τ23+βσ23,(1b)τ13+bτ23+βσ13+bσ23=Cτ12+βσ12≤τ23+βσ23,where τ12, τ13, and τ23 are the principal shear stresses acting on the orthogonal octahedron element, defined as τij=σi−σj/2, and subscripts “i” and “j” denote 1, 2, and 3. σ12, σ13, and σ23 represent the normal stresses acting on the orthogonal octahedron element corresponding to the principal shear stress, defined as σij=σi+σj/2. σ1, σ2, and σ3 are the maximum, intermediate, and minimum principal stresses, respectively. The coefficients of β and C are material parameters that can be determined experimentally. In practice, the USC can be also expressed by the following Mohr-Coulomb strength parameters [25]:(2a)σ1−1−sinφ1+b1+sinφbσ2+σ3=2ccosφ1+sinφσ2≤σ1+σ32+σ1−σ32sinφ,(2b)bσ2+σ11+b−1−sinφ1+sinφσ3=2ccosφ1+sinφσ2≥σ1+σ32+σ1−σ32sinφ.In equations (2a) and (2b), tensile stress is positive, and effective stress should be used under undrained conditions. As in previous reports [18, 26, 27], the parameter b is shown in equation (3), where σc, σt, and τ0 denote, respectively, uniaxial compressive strength, uniaxial tensile strength, and shear strength:(3)b=σc+σtτ0−σtσcσt−τ0σc.As shown in Figure1, when the value of the parameter b varies between 0 and 1, a series of yield criteria will be obtained, which consider the effect of intermediate principal stress to describe the failure behaviour of different rock materials subjected to a complex stress state. From the analysis above, the USC can be applied to materials with different properties and, as such, has a wider application range. In solving the problem of tunnel convergence, the compressive pressure is ordinarily taken to be positive, and the USC in equations (2a) and (2b) can be changed as follows [19, 25]:(4a)1−sinφ1+sinφσ1−bσ2+σ31+b=2ccosφ1+sinφσ2≤σ1+σ32−σ1−σ32sinφ,(4b)bσ2+σ11−sinφ1+b1+sinφ−σ3=2ccosφ1+sinφσ2≥σ1+σ32−σ1−σ32sinφ.Figure 1
Limiting loci of USC in deviatoric plane.


## 3. Tunnel Condition and Basic Equation

### 3.1. Definition of the Problem

In Figure2, the depicted circular tunnel of radius R0 is in an elastic state following excavation. After blasting, a cylindrical BIDZ of radius Rd will be developed around the tunnel with lower rock mass parameters. As a result of installing a lining, a uniform internal pressure is applied to the tunnel surface. The elastic-brittle-plastic model is used in this paper (see Figure 3), which means that the strength of the rock mass decreases suddenly after peak strength is attained, and the residual values for cohesion, internal friction angle, Poisson’s ratio, and elastic modulus (cr, φr, νr, and Er, respectively) are used to analyze the stresses and rock mass deformation in the plastic zone. It is worth noting that the analytical solution based on the elastic-perfectly-plastic model can be considered a special case in which the residual rock parameter values are equal to the initial values of c, φ, ν, and E). According to the range of the plastic zone and BIDZ, two different cases can be considered, as shown in Figure 2: Case 1 is that in which the BIDZ radius is less than that of the plastic zone, and Case 2 represents the 100% overlap between the BIDZ and plastic zone.Figure 2
Deep tunnel subjected to a hydrostatic stress field with gravity effect and cylindrical BIDZ: (a)Rp>Rd; (b) Rp=Rd.
(a)(b)Figure 3
Simplified mechanical behaviour of rock material in this study for both plastic zone and damaged zone. (a) Stress-strain relationship; (b) flow rule.
(a)(b)The axial stressσz is the intermediate principal stress, which can be expressed as σz=σ2=nσr+σθ/2, and 2ν≤n≤1 [18]. ν is the initial Poisson ratio of the rock mass. To simplify calculation, n=1 can be taken in the damage and plastic zones, and the relationship σr+σθ/2≥σr+σθ/2−σθ−σrsinφ/2 can be deduced. Thus, equation (4b) should be embraced and afterwards can be changed as follows:(5a)σθ=ξσr+Yfor peak strength,(5b)σθ=ξrσr+Yrfor plastic zone,where(6)ξ=2+b+2+3bsinφ2+b1−sinφ,Y=41+bccosφ2+b1−sinφ,ξr=2+b+2+3bsinφr2+b1−sinφr,Yr=41+bcrcosφr2+b1−sinφr.The failure parameters of damaged rock mass are denoted by subscript “D” (namely, cRd, φRd, νRd, and ERd).


### 3.2. Plastic Potential Functions

In deformation analysis of a circular tunnel, an appropriate plastic potential function needs to be determined according to the incremental theory of plasticity. Because the plastic strain increments depend on the plastic potential function, when the plastic potential function is consistent with the yield criterion, the function is called the associated flow rule; otherwise, it is called the nonassociated flow rule. It is believed, in general, that plastic deformation satisfies the nonassociated flow rule [6, 18, 19]. Thus, it is assumed that the plastic potential function can be expressed as(7)Fσθ,σr=σθ−βσr,where β=2+b+2+3bsinψ/2+b1−sinψ is the dilation parameter associated with the coefficient of intermediate principal stress b and the dilation angle ψ. It is generally assumed that the plastic strain increments are proportional to the stress gradient of the plastic potential function F [2]:(8a)dεθp=∂F∂σθdλ,(8b)dεrp=∂F∂σrdλ,where dεθp and dεrp are the circumferential and radial plastic strain increments, respectively; dλ is a nonnegative proportional constant. Submitting equations (7), (8a), and (8b), the contact between the plastic strain and dilation parameter is obtained:(9)εrp+βεθp=0.Most previous studies regardedβ=1+sinψ/1−sinψ as a dilation parameter to investigate the displacement behaviour of rock mass in the fracture zone while ignoring the influence of intermediate principal stress on the rock dilation effect, which is potentially unreasonable.


### 3.3. Equilibrium Equation considering Gravity’s Effect

The rock weight of the fracture zone (i.e., plastic and damaged zones) at the tunnel roof can apply higher additional pressure on the support system and significantly affect tunnel stability. Hence, the rock weight of the fracture zone should be considered in the equilibrium equation. Detournay [28] advanced a model describing the body forces acting towards the centre of an opening, the methods of which have been widely applied to elastic-plastic analysis of a circular tunnel at a roof location [2, 21, 22, 24, 29]. The same method is used in the present work. Assuming that the gravity load is located directly in the centre of the tunnel, the equilibrium equation is as follows [24]:(10)∂σr∂r+σr−σθr+γ=0,where γ is the rock mass unit weight in the fracture zone.


## 3.1. Definition of the Problem

In Figure2, the depicted circular tunnel of radius R0 is in an elastic state following excavation. After blasting, a cylindrical BIDZ of radius Rd will be developed around the tunnel with lower rock mass parameters. As a result of installing a lining, a uniform internal pressure is applied to the tunnel surface. The elastic-brittle-plastic model is used in this paper (see Figure 3), which means that the strength of the rock mass decreases suddenly after peak strength is attained, and the residual values for cohesion, internal friction angle, Poisson’s ratio, and elastic modulus (cr, φr, νr, and Er, respectively) are used to analyze the stresses and rock mass deformation in the plastic zone. It is worth noting that the analytical solution based on the elastic-perfectly-plastic model can be considered a special case in which the residual rock parameter values are equal to the initial values of c, φ, ν, and E). According to the range of the plastic zone and BIDZ, two different cases can be considered, as shown in Figure 2: Case 1 is that in which the BIDZ radius is less than that of the plastic zone, and Case 2 represents the 100% overlap between the BIDZ and plastic zone.Figure 2
Deep tunnel subjected to a hydrostatic stress field with gravity effect and cylindrical BIDZ: (a)Rp>Rd; (b) Rp=Rd.
(a)(b)Figure 3
Simplified mechanical behaviour of rock material in this study for both plastic zone and damaged zone. (a) Stress-strain relationship; (b) flow rule.
(a)(b)The axial stressσz is the intermediate principal stress, which can be expressed as σz=σ2=nσr+σθ/2, and 2ν≤n≤1 [18]. ν is the initial Poisson ratio of the rock mass. To simplify calculation, n=1 can be taken in the damage and plastic zones, and the relationship σr+σθ/2≥σr+σθ/2−σθ−σrsinφ/2 can be deduced. Thus, equation (4b) should be embraced and afterwards can be changed as follows:(5a)σθ=ξσr+Yfor peak strength,(5b)σθ=ξrσr+Yrfor plastic zone,where(6)ξ=2+b+2+3bsinφ2+b1−sinφ,Y=41+bccosφ2+b1−sinφ,ξr=2+b+2+3bsinφr2+b1−sinφr,Yr=41+bcrcosφr2+b1−sinφr.The failure parameters of damaged rock mass are denoted by subscript “D” (namely, cRd, φRd, νRd, and ERd).


## 3.2. Plastic Potential Functions

In deformation analysis of a circular tunnel, an appropriate plastic potential function needs to be determined according to the incremental theory of plasticity. Because the plastic strain increments depend on the plastic potential function, when the plastic potential function is consistent with the yield criterion, the function is called the associated flow rule; otherwise, it is called the nonassociated flow rule. It is believed, in general, that plastic deformation satisfies the nonassociated flow rule [6, 18, 19]. Thus, it is assumed that the plastic potential function can be expressed as(7)Fσθ,σr=σθ−βσr,where β=2+b+2+3bsinψ/2+b1−sinψ is the dilation parameter associated with the coefficient of intermediate principal stress b and the dilation angle ψ. It is generally assumed that the plastic strain increments are proportional to the stress gradient of the plastic potential function F [2]:(8a)dεθp=∂F∂σθdλ,(8b)dεrp=∂F∂σrdλ,where dεθp and dεrp are the circumferential and radial plastic strain increments, respectively; dλ is a nonnegative proportional constant. Submitting equations (7), (8a), and (8b), the contact between the plastic strain and dilation parameter is obtained:(9)εrp+βεθp=0.Most previous studies regardedβ=1+sinψ/1−sinψ as a dilation parameter to investigate the displacement behaviour of rock mass in the fracture zone while ignoring the influence of intermediate principal stress on the rock dilation effect, which is potentially unreasonable.


## 3.3. Equilibrium Equation considering Gravity’s Effect

The rock weight of the fracture zone (i.e., plastic and damaged zones) at the tunnel roof can apply higher additional pressure on the support system and significantly affect tunnel stability. Hence, the rock weight of the fracture zone should be considered in the equilibrium equation. Detournay [28] advanced a model describing the body forces acting towards the centre of an opening, the methods of which have been widely applied to elastic-plastic analysis of a circular tunnel at a roof location [2, 21, 22, 24, 29]. The same method is used in the present work. Assuming that the gravity load is located directly in the centre of the tunnel, the equilibrium equation is as follows [24]:(10)∂σr∂r+σr−σθr+γ=0,where γ is the rock mass unit weight in the fracture zone.


## 4. Analysis of Stresses and Deformations

### 4.1. Analysis Solution of Elastic Zone

The rock mass remains elastic outside the fracture zone, as shown in Figure2, and the solutions are similar for Cases 1 and 2 mentioned above. Under axisymmetric plane strain conditions, the elastic solution given by Lame can be used to calculate the radial displacement and the stresses of the elastic zone as follows [17, 24].For Case 1,(11a)σrr=σ0−σ0−σRpRp2r2,(11b)σθr=σ0+σ0−σRpRp2r2,(11c)urr=1+vEσ0−σRPRP2r.In Case 2, the stresses and deformation of the tunnel in the elastic zone can be obtained only by substitutingRp and σRp by Rd and σRd in equations (11a)–(11c). σRp and σRd are the radial boundary stresses at the interface between the elastic and fracture zones, respectively. The radial boundary stress should satisfy the initial yield condition (equation (5a)) at the elastic-fracture zone interface. Hence, the radial boundary stress can be obtained by substituting equations (11a) and (11b) into equation (5a):(12)σRp=σRd=2σ0−Y1+ξ.


### 4.2. Analysis Solutions of Fracture Zone

#### 4.2.1. Stresses and Radius Analysis

(1) For Case 1. As shown in Figure 2(b), the fracture zone mainly includes the plastic zone and BIDZ. Substituting equation (5b) into equation (10) and applying the boundary condition σr=σRp at r=Rp, the radial and circumferential stresses of the plastic zone can be expressed as follows:(13a)σrr=ArRpr1−ξr+Brr−Rp2−ξrrξr−1+Cr,(13b)σθr=ξrArRpr1−ξr+ξrBrr−Rp2−ξrrξr−1+Cr,where Ar=σRp−Cr, Br=γ/2−ξr, and Cr=Yr/1−ξr.By considering the blast-induced damage effect in the same way, substituting equation (5b) into equation (10), and using the boundary condition σr=σi at r=R0, the radial and circumferential stresses for the BIDZ can be obtained:(14a)σrr=ARdR0r1−ξRd+BRdr−R02−ξRdrξRd−1+CRd,(14b)σθr=ξRdARdR0r1−ξRd+ξRdBRdr−R02−ξRdrξRd−1+CRd,where ARd=σi−CRd, BRd=γ/2−ξRd, and CRd=YRd/1−ξRd.The closed solution is acquired by substitutingr=Rd into equations (13a) and (14a) and then equating the right-hand sides of equations (13a) and (14a):(15)ARdR0Rd1−ξRd+BRdRd−R02−ξRdRdξRd−1+CRd=ArRpRd1−ξr+BrRd−Rp2−ξrRdξr−1+Cr.From the above equation, the radius of the BIDZ (Rd) is known, and the approximate solutions of Rp can then be obtained using numerical iterations as in the bisection or Newton-Raphson methods. The bisection method is used herein to solve equation (15). Once the Rp values are determined from equation (15), substitute the value into equations (11a)–(11c) to get the stress and displacement in the elastic zone.(2) For Case 2. The fracture zone only consists of a BIDZ with a lower strength parameter of rocks. Just as equations (14a) and (14b), the radial and circumferential stress of the BIDZ in Case 2 is consistent with that in Case 1. By considering the continuity of radial stress at the elastic-BIDZ interface, the radius of the BIDZ can be determined. Then, substituting r=Rd into equation (14a), followed by equating the right-hand sides of equations (12) and (14a), gives(16)ARdR0Rd1−ξRd+BRdRd−R02−ξRdRdξRd−1+CRd=2σ0−Y1+ξ.The approximate solution forRd can also be calculated by using the method of bisection.


#### 4.2.2. Displacement Analysis

The total circumferential strainεθr and radial strain εrr at radius r for the axial-symmetric condition can be expressed by the radial displacement urr [1]:(17)εrr=∂ur∂r,εθr=εθe+εθp.Assuming that the total strain in the fracture zone consists of a plastic strain component and an elastic strain component, the total strain in the fracture zone can then be expressed as follows:(18a)εrr=εre+εrp,(18b)εθr=εθe+εθp.By adopting Hooke’s law, the elastic strain componentsεre and εθe in the fracture zone can be expressed as follows [1, 18]:(19a)εre=1+νrEr1−νrσrr−σ0−νrσθr−σ0,(19b)εθe=1+νrEr1−νrσθr−σ0−νrσrr−σ0.Substituting equations (9) and (17) into equations (19a) and (19b) gives(20)∂ur∂r+βur=εre+βεθe=fr.(1) For Case 1. Upon solving equation (20), the general solution for the radial displacement in the plastic zone can be deduced using the displacement continuity condition at the elastic-plastic interface; that is,(21)ur=1rβ∫Rprrβεre+βεθedr+uRpRprβ=1rβ∫Rprrβfrdr+uRpRprβ,where uRp=1+ν/Eσ0−σRpRp.Substituting equations (13a), (13b), (19a), and (19b) into equation fr, we get the following equation:(22)fr=1+νrErδa1Rpr1−ξr+δa2r−Rp2−ξrrξr−1+δa3,where δa1, δa2, and δa3 are constants. δa1=δ1Ar+δ2ξrAr, δa2=δ1Br+δ2ξrBr, δa3=δ1+δ2Cr+2νr−11+βσ0, δ1=1−νr−βνr, δ2=β−νr−βνr, and uRp=1+ν/Eσ0−σRpRp.Upon substituting equation22) into equation (21),(23)ur=1rβ1+νrErδa1ξr+βRd1−ξr−δa2ξr+βRp2−ξrrξr+β−Rpξr+β+δa22+βr2+β−Rp2+β+δa31+βr1+β−Rp1+β+uRpRprβ.As shown in equation (23), the radial displacement in the BIDZ can also be deduced:(24)ur=1rβ1+νRdERdδb1ξRd+βRd1−ξRd−δb2ξRd+βRd2−ξRdrξRd+β−RdξRd+β+δb22+βr2+β−Rd2+β+δb31+βr1+β−Rd1+β+uRdRdrβ,where δa1, δa2, and δa3 are constants. The radial displacement uRd at the plastic-BIDZ interface is obtained by substituting r=Rd into equation (23). δb1=δ3ARd+δ4ξRdARd, δb2=δ3BRd+δ4ξRdBRd, δb3=δ3+δ4CRd+2νRd−11+βσ0, δ3=1−νRd−βνRd, and δ4=β−νRd−βνRd.(2) For Case 2. The postpeak rock only consists of the BIDZ. As in equation (24), the radial displacement for the BIDZ can be calculated by substituting equations (14a), (14b), (19a), (19b), and (22) into equation (21) and considering the radial displacement continuity condition at the elastic-BIDZ interface.(25)ur=1rβ1+νRdERdδb1ξRd+βRd1−ξr−δb2ξRd+βRd2−ξrrξRd+β−RdξRd+β+δb22+βr2+β−Rd2+β+δb31+βr1+β−Rd1+β+uRdRdrβ,where uRd=1+ν/Eσ0−σRdRd.The tunnel roof convergenceu0 can be determined by substituting r=R0 into equations (24) and (25) for various cases, respectively.


### 4.3. Theoretical Value and Comparative Analysis

The solution for a circular opening consists of a series of simple formulas, which obviate the need to solve a complex differential equation. The solution takes into consideration rock gravity and the BIDZ, which can produce a relatively accurate result. Hence, the new solution has a wider engineering application.Zhang et al. [18] devised a closed-form solution and analyzed the stress and displacement of a circular opening based on the USC. However, their solution did not account for the effect of rock gravity, which is inconsistent with engineering practice. Chen et al. [19] also used the USC and presented a closed-form solution of a deep circular tunnel based on a new four-stage constitutive model. However, their solution only considered the cohesion degradation while ignoring the attenuation of other rock parameters. The rock gravity and BIDZ effects were also not considered in Chen’s solution. However, these solutions ignore the intermediate principal stress. The aggregate effect of rock gravity, BIDZ, and intermediate principal stress is, however, rationally considered in the new closed-form solution presented herein. The closed-form solution for Case 2, in which γ=0, cRd=cr, φRd=φr, νRd=νr, and ERd=Er, is identical to Zhang’s solution (n=0) [18]. Park’s solution, based on the M-C failure criterion, can be derived using γ=0, cRd=cr, φRd=φr, νRd=ν, ERd=E, and b=0 [6]. Zhang and Park’s solutions are special cases in this paper. Hence, this new solution includes a series of traditional results and is widely used in rock mass.


### 4.4. Further Analysis on Modified USC

#### 4.4.1. Numerical Verification

Cheng [30] embedded the ideal elastoplastic constitutive model into ABAQUS finite-element software based on the USC and then simulated the stress distribution around a circular tunnel under conditions of various b values. Table 1 lists the geometry and rock parameters in this simulation. The associated flow rule, that is, ψ=φ, was adopted to determine the dilatancy coefficient β. However, the numerical simulation does not consider the effect of intermediate principal stress on the dilatancy coefficient. Therefore, β=1+sinψ/1−sinψ is used in this analytical solution. Since the damaged zone is ignored in Pan’s numerical simulation, Case 2 is used to facilitate a comparison with Pan’s simulation results. The analysis results on Figure 4 are highly consistent with Pan’s simulation results, which further confirm the accuracy of the new solution.Table 1
The geometry and rock parameters by Pan et al. [30].
Geometry and rock parametersValueIn situ stress,σ0 (MPa)30Tunnel radius,R0 (m)10Rock gravity,γ (kN/m3)0Internal pressure,σi (MPa)0Dilatancy angle,ψ(°)30Initial elastic modulus,E (GPa)2Initial cohesive,c (MPa)6Initial friction angle,φ (°)30Initial Poisson’s ratio (ν)0.3Damage elastic modulus,ERd (GPa)2Damage cohesive,cRd (MPa)6Damage friction angle,φRd (°)30Damage Poisson’s ratio (νRd)0.3Figure 4
Comparison of stress distribution under (a)b = 0.0, (b) b = 0.5, and (c) b = 1.0.
(a)(b)(c)


#### 4.4.2. Parametric Study

A 3 m radius transportation roadway of Quandian coal mine in China is buried at about 650 m, with 15 MPa average in situ stress. The surrounding rock is mainly sandy mudstone. The special rock parameters of the roadway are shown in Table2. It should be noted that the 0.2 times of the initial strength parameters were used as the damage parameters, reflecting a serious damage phenomenon induced by blasting excavation.Table 2
The geometry and rock parameters used in this study.
Geometry and rock parametersValueIn situ stress,σ0 (MPa)15Tunnel radius,R0 (m)3Rock gravity,γ (kN/m3)27.5Dilatancy angle,ψ (°)0Initial elastic modulus,E (GPa)2.2Initial cohesive,c (MPa)5Initial friction angle,φ (°)30Initial Poisson’s ratio (ν)0.25Residual elastic modulus,Er (GPa)1.6Residual cohesive,cr (MPa)3.6Residual friction angle,φr (°)24Residual Poisson’s ratio (νr)0.25Damage elastic modulus,ERd (GPa)0.44Damage cohesive,cRd (MPa)1.0Damage friction angle,φRd (°)6Damage Poisson’s ratio (νRd)0.25


#### 4.4.3. The Effect of Intermediate Principal Stress and Gravity on Ground Response Curves

Case 2 is used here to investigate the effect of intermediate principal stress and gravity on the ground response curve (GRC). Figure5 shows the GRC of the tunnel ignoring (Figure 5(a)) and considering (Figure 5(b)) the gravity effect under various intermediate principal stresses. The main research results are as follows.Figure 5
GRC at different gravity condition and intermediate principal stress. (a) Ignoring the gravity effect. (b) Considering the gravity effect.
(a)(b)As shown in Table3, with increasing b value, the tunnel surface convergence gradually decreases. When the b value increases from 0 to 0.5 and 1, the tunnel surface convergence decreases by 71.97% and 85.57% with the consideration of the gravity as well as by 68.93% and 83.48% with ignoring the gravity, respectively. The above analysis reveals that high intermediate principal stress can effectively improve the bearing capacity of surrounding rock and restrain tunnel convergence.Table 3
Roof surface displacement of the tunnel under various gravity conditions.
ParametersConditionsb0.00.250.50.751.0u0 (m)Gravity included1.0880.5070.3050.2100.157Gravity ignored0.9140.4580.2840.1990.151Considering the gravity condition, the roof surface displacement of the tunnel increases by 19.04% forb = 0, 7.39% for b = 0.5, and 3.79% for b = 1.0 compared with ignoring the gravity, respectively. This means that the gravitational condition can accelerate tunnel convergence, which is not conducive to long-term tunnel stability. Thus, the tunnel’s support strength design should fully consider the rock gravity at the roof. Moreover, with increasing intermediate principal stress, the rock gravity effect at the roof gradually weakens, again indicating that intermediate principal stress is beneficial to tunnel stability.With increasingb value, the critical pressure at the elastic-damage zone interface gradually decreases, demonstrating that the higher the intermediate principal stress, the smaller the support resistance needed to maintain tunnel stability.


#### 4.4.4. The Effect of Gravity and Intermediate Principal Stress on Damage Zone Radius

Figure6 shows the trend of damage zone radius against the internal pressure at various intermediate principal stress and gravity conditions. The damage zone radii are listed in Table 4, from which the following can be seen.Figure 6
Evolution law of damage zone at various gravity conditions and intermediate principal stress.Table 4
Damage zone radius of the tunnel under various gravity conditions.
ParametersConditionsb0.00.250.50.751.0Rd/R0Gravity included3.7262.6142.0961.8031.616Gravity ignored3.4272.4962.0341.7651.591For the damage zone, as the internal pressure decreases gradually, the dimensionless valueRd/R0 nonlinearly increases at various b values, while the increase rates vary with the b value. The slope of the “σi vs. Rd/R0” curve for b = 0 is obviously less than that b = 1 at any gravity condition, which indicates that the extension process of the damage zone is closely dependent on the intermediate principal stress.Forσi=0MPa, with increasing the b value, the radius of the damage zone gradually decreases. When the b value increases from 0 to 0.5 and 1, the dimensionless value Rd/R0 decreases by 43.75% and 56.63% with the consideration of the gravity as well as by 40.65% and 53.57% with ignoring the gravity effect, indicating that the higher the intermediate principal stress, the smaller the damage zone radius. In practice, the damage zone is the main support area for the deep tunnel due to the lower rock strength. One of the most important functions of the bolt is to reinforce the damaged block and improve its residual strength. The support parameter design of bolts is dependent on the range of damage zones. At high intermediate principal stress, the range of damage zone is obviously smaller than that at lower intermediate principal stress. It means that we can optimize bolt support parameters by the methods of reducing its length and increasing its pretension for the tunnel under the high intermediate principal stress state to save its support cost.The roof rock gravity has an essential influence on the radius of damage zone. Considering the gravity effect, the dimensionless value ofRd/R0 increases by 8.72% for b = 0, by 3.05% for b = 0.5, and by 1.57% for b = 1.0 compared with ignoring the gravity, respectively. It means that the range of damage zone under the effect of gravity is larger than that without the influence of gravity. Hence, the roof support design of the tunnel should consider the effect of roof rock gravity for maintaining the long-term stability of the tunnel.


#### 4.4.5. The Effect of Damage Zone Range on the Stress and Displacement

The stress and deformation of the tunnel are dependent on the damage zone range. In this work, Case 1 is used to explore the effect of damage zone radius on the stress and deformation of the tunnel while ignoring the gravity effect (γ=0). Figure 7 plots the trend of hoop stress and radial stress against the damage zone radius. With increasing the damage zone radius, the location of the maximum hoop stress dramatically shifted away from the surface of the tunnel, while its value has only lightly changed, indicating that the range of damage zones has an essential influence on the stress state of the circular tunnel. Therefore, the larger the radius of the damage zone, the larger the looseness range of the tunnel and the higher the instability risk of the tunnel.Figure 7
Hoop stress and radial stress of tunnel under various damage ranges (γ=0 and b=0). (a) Hoop stress. (b) Radial stress.
(a)(b)Figure8 shows the relationship between radial displacement and radius. According to the figure, the extent of the damage zone has a profound influence on the large deformation instability of the tunnel. The original radial displacement of the tunnel is always at a lower level of the plastic zone-damage zone interface regardless of the existence of any damage zone range. Once the rock mass at the roof enters the damage zone, the radial displacement of the tunnel sharply increases in a nonlinear manner, especially under Rd=9m. Furthermore, with increasing damage zone radius, the surface displacement of the tunnel observably increases, indicating that a large damage zone range is not conducive to tunnel stability. This is because the larger the damage zone range, the larger the surface displacement of the tunnel and then the larger the deformation pressure acting on the support structure, which leads to the collapse of the support structure due to the overload. Hence, the support methods of the tunnel should fully consider the effect of BIDZ range on the ground response.Figure 8
Displacement change of tunnel under various damage ranges (γ=0 and b=0).


## 4.1. Analysis Solution of Elastic Zone

The rock mass remains elastic outside the fracture zone, as shown in Figure2, and the solutions are similar for Cases 1 and 2 mentioned above. Under axisymmetric plane strain conditions, the elastic solution given by Lame can be used to calculate the radial displacement and the stresses of the elastic zone as follows [17, 24].For Case 1,(11a)σrr=σ0−σ0−σRpRp2r2,(11b)σθr=σ0+σ0−σRpRp2r2,(11c)urr=1+vEσ0−σRPRP2r.In Case 2, the stresses and deformation of the tunnel in the elastic zone can be obtained only by substitutingRp and σRp by Rd and σRd in equations (11a)–(11c). σRp and σRd are the radial boundary stresses at the interface between the elastic and fracture zones, respectively. The radial boundary stress should satisfy the initial yield condition (equation (5a)) at the elastic-fracture zone interface. Hence, the radial boundary stress can be obtained by substituting equations (11a) and (11b) into equation (5a):(12)σRp=σRd=2σ0−Y1+ξ.


## 4.2. Analysis Solutions of Fracture Zone

### 4.2.1. Stresses and Radius Analysis

(1) For Case 1. As shown in Figure 2(b), the fracture zone mainly includes the plastic zone and BIDZ. Substituting equation (5b) into equation (10) and applying the boundary condition σr=σRp at r=Rp, the radial and circumferential stresses of the plastic zone can be expressed as follows:(13a)σrr=ArRpr1−ξr+Brr−Rp2−ξrrξr−1+Cr,(13b)σθr=ξrArRpr1−ξr+ξrBrr−Rp2−ξrrξr−1+Cr,where Ar=σRp−Cr, Br=γ/2−ξr, and Cr=Yr/1−ξr.By considering the blast-induced damage effect in the same way, substituting equation (5b) into equation (10), and using the boundary condition σr=σi at r=R0, the radial and circumferential stresses for the BIDZ can be obtained:(14a)σrr=ARdR0r1−ξRd+BRdr−R02−ξRdrξRd−1+CRd,(14b)σθr=ξRdARdR0r1−ξRd+ξRdBRdr−R02−ξRdrξRd−1+CRd,where ARd=σi−CRd, BRd=γ/2−ξRd, and CRd=YRd/1−ξRd.The closed solution is acquired by substitutingr=Rd into equations (13a) and (14a) and then equating the right-hand sides of equations (13a) and (14a):(15)ARdR0Rd1−ξRd+BRdRd−R02−ξRdRdξRd−1+CRd=ArRpRd1−ξr+BrRd−Rp2−ξrRdξr−1+Cr.From the above equation, the radius of the BIDZ (Rd) is known, and the approximate solutions of Rp can then be obtained using numerical iterations as in the bisection or Newton-Raphson methods. The bisection method is used herein to solve equation (15). Once the Rp values are determined from equation (15), substitute the value into equations (11a)–(11c) to get the stress and displacement in the elastic zone.(2) For Case 2. The fracture zone only consists of a BIDZ with a lower strength parameter of rocks. Just as equations (14a) and (14b), the radial and circumferential stress of the BIDZ in Case 2 is consistent with that in Case 1. By considering the continuity of radial stress at the elastic-BIDZ interface, the radius of the BIDZ can be determined. Then, substituting r=Rd into equation (14a), followed by equating the right-hand sides of equations (12) and (14a), gives(16)ARdR0Rd1−ξRd+BRdRd−R02−ξRdRdξRd−1+CRd=2σ0−Y1+ξ.The approximate solution forRd can also be calculated by using the method of bisection.


### 4.2.2. Displacement Analysis

The total circumferential strainεθr and radial strain εrr at radius r for the axial-symmetric condition can be expressed by the radial displacement urr [1]:(17)εrr=∂ur∂r,εθr=εθe+εθp.Assuming that the total strain in the fracture zone consists of a plastic strain component and an elastic strain component, the total strain in the fracture zone can then be expressed as follows:(18a)εrr=εre+εrp,(18b)εθr=εθe+εθp.By adopting Hooke’s law, the elastic strain componentsεre and εθe in the fracture zone can be expressed as follows [1, 18]:(19a)εre=1+νrEr1−νrσrr−σ0−νrσθr−σ0,(19b)εθe=1+νrEr1−νrσθr−σ0−νrσrr−σ0.Substituting equations (9) and (17) into equations (19a) and (19b) gives(20)∂ur∂r+βur=εre+βεθe=fr.(1) For Case 1. Upon solving equation (20), the general solution for the radial displacement in the plastic zone can be deduced using the displacement continuity condition at the elastic-plastic interface; that is,(21)ur=1rβ∫Rprrβεre+βεθedr+uRpRprβ=1rβ∫Rprrβfrdr+uRpRprβ,where uRp=1+ν/Eσ0−σRpRp.Substituting equations (13a), (13b), (19a), and (19b) into equation fr, we get the following equation:(22)fr=1+νrErδa1Rpr1−ξr+δa2r−Rp2−ξrrξr−1+δa3,where δa1, δa2, and δa3 are constants. δa1=δ1Ar+δ2ξrAr, δa2=δ1Br+δ2ξrBr, δa3=δ1+δ2Cr+2νr−11+βσ0, δ1=1−νr−βνr, δ2=β−νr−βνr, and uRp=1+ν/Eσ0−σRpRp.Upon substituting equation22) into equation (21),(23)ur=1rβ1+νrErδa1ξr+βRd1−ξr−δa2ξr+βRp2−ξrrξr+β−Rpξr+β+δa22+βr2+β−Rp2+β+δa31+βr1+β−Rp1+β+uRpRprβ.As shown in equation (23), the radial displacement in the BIDZ can also be deduced:(24)ur=1rβ1+νRdERdδb1ξRd+βRd1−ξRd−δb2ξRd+βRd2−ξRdrξRd+β−RdξRd+β+δb22+βr2+β−Rd2+β+δb31+βr1+β−Rd1+β+uRdRdrβ,where δa1, δa2, and δa3 are constants. The radial displacement uRd at the plastic-BIDZ interface is obtained by substituting r=Rd into equation (23). δb1=δ3ARd+δ4ξRdARd, δb2=δ3BRd+δ4ξRdBRd, δb3=δ3+δ4CRd+2νRd−11+βσ0, δ3=1−νRd−βνRd, and δ4=β−νRd−βνRd.(2) For Case 2. The postpeak rock only consists of the BIDZ. As in equation (24), the radial displacement for the BIDZ can be calculated by substituting equations (14a), (14b), (19a), (19b), and (22) into equation (21) and considering the radial displacement continuity condition at the elastic-BIDZ interface.(25)ur=1rβ1+νRdERdδb1ξRd+βRd1−ξr−δb2ξRd+βRd2−ξrrξRd+β−RdξRd+β+δb22+βr2+β−Rd2+β+δb31+βr1+β−Rd1+β+uRdRdrβ,where uRd=1+ν/Eσ0−σRdRd.The tunnel roof convergenceu0 can be determined by substituting r=R0 into equations (24) and (25) for various cases, respectively.


## 4.2.1. Stresses and Radius Analysis

(1) For Case 1. As shown in Figure 2(b), the fracture zone mainly includes the plastic zone and BIDZ. Substituting equation (5b) into equation (10) and applying the boundary condition σr=σRp at r=Rp, the radial and circumferential stresses of the plastic zone can be expressed as follows:(13a)σrr=ArRpr1−ξr+Brr−Rp2−ξrrξr−1+Cr,(13b)σθr=ξrArRpr1−ξr+ξrBrr−Rp2−ξrrξr−1+Cr,where Ar=σRp−Cr, Br=γ/2−ξr, and Cr=Yr/1−ξr.By considering the blast-induced damage effect in the same way, substituting equation (5b) into equation (10), and using the boundary condition σr=σi at r=R0, the radial and circumferential stresses for the BIDZ can be obtained:(14a)σrr=ARdR0r1−ξRd+BRdr−R02−ξRdrξRd−1+CRd,(14b)σθr=ξRdARdR0r1−ξRd+ξRdBRdr−R02−ξRdrξRd−1+CRd,where ARd=σi−CRd, BRd=γ/2−ξRd, and CRd=YRd/1−ξRd.The closed solution is acquired by substitutingr=Rd into equations (13a) and (14a) and then equating the right-hand sides of equations (13a) and (14a):(15)ARdR0Rd1−ξRd+BRdRd−R02−ξRdRdξRd−1+CRd=ArRpRd1−ξr+BrRd−Rp2−ξrRdξr−1+Cr.From the above equation, the radius of the BIDZ (Rd) is known, and the approximate solutions of Rp can then be obtained using numerical iterations as in the bisection or Newton-Raphson methods. The bisection method is used herein to solve equation (15). Once the Rp values are determined from equation (15), substitute the value into equations (11a)–(11c) to get the stress and displacement in the elastic zone.(2) For Case 2. The fracture zone only consists of a BIDZ with a lower strength parameter of rocks. Just as equations (14a) and (14b), the radial and circumferential stress of the BIDZ in Case 2 is consistent with that in Case 1. By considering the continuity of radial stress at the elastic-BIDZ interface, the radius of the BIDZ can be determined. Then, substituting r=Rd into equation (14a), followed by equating the right-hand sides of equations (12) and (14a), gives(16)ARdR0Rd1−ξRd+BRdRd−R02−ξRdRdξRd−1+CRd=2σ0−Y1+ξ.The approximate solution forRd can also be calculated by using the method of bisection.


## 4.2.2. Displacement Analysis

The total circumferential strainεθr and radial strain εrr at radius r for the axial-symmetric condition can be expressed by the radial displacement urr [1]:(17)εrr=∂ur∂r,εθr=εθe+εθp.Assuming that the total strain in the fracture zone consists of a plastic strain component and an elastic strain component, the total strain in the fracture zone can then be expressed as follows:(18a)εrr=εre+εrp,(18b)εθr=εθe+εθp.By adopting Hooke’s law, the elastic strain componentsεre and εθe in the fracture zone can be expressed as follows [1, 18]:(19a)εre=1+νrEr1−νrσrr−σ0−νrσθr−σ0,(19b)εθe=1+νrEr1−νrσθr−σ0−νrσrr−σ0.Substituting equations (9) and (17) into equations (19a) and (19b) gives(20)∂ur∂r+βur=εre+βεθe=fr.(1) For Case 1. Upon solving equation (20), the general solution for the radial displacement in the plastic zone can be deduced using the displacement continuity condition at the elastic-plastic interface; that is,(21)ur=1rβ∫Rprrβεre+βεθedr+uRpRprβ=1rβ∫Rprrβfrdr+uRpRprβ,where uRp=1+ν/Eσ0−σRpRp.Substituting equations (13a), (13b), (19a), and (19b) into equation fr, we get the following equation:(22)fr=1+νrErδa1Rpr1−ξr+δa2r−Rp2−ξrrξr−1+δa3,where δa1, δa2, and δa3 are constants. δa1=δ1Ar+δ2ξrAr, δa2=δ1Br+δ2ξrBr, δa3=δ1+δ2Cr+2νr−11+βσ0, δ1=1−νr−βνr, δ2=β−νr−βνr, and uRp=1+ν/Eσ0−σRpRp.Upon substituting equation22) into equation (21),(23)ur=1rβ1+νrErδa1ξr+βRd1−ξr−δa2ξr+βRp2−ξrrξr+β−Rpξr+β+δa22+βr2+β−Rp2+β+δa31+βr1+β−Rp1+β+uRpRprβ.As shown in equation (23), the radial displacement in the BIDZ can also be deduced:(24)ur=1rβ1+νRdERdδb1ξRd+βRd1−ξRd−δb2ξRd+βRd2−ξRdrξRd+β−RdξRd+β+δb22+βr2+β−Rd2+β+δb31+βr1+β−Rd1+β+uRdRdrβ,where δa1, δa2, and δa3 are constants. The radial displacement uRd at the plastic-BIDZ interface is obtained by substituting r=Rd into equation (23). δb1=δ3ARd+δ4ξRdARd, δb2=δ3BRd+δ4ξRdBRd, δb3=δ3+δ4CRd+2νRd−11+βσ0, δ3=1−νRd−βνRd, and δ4=β−νRd−βνRd.(2) For Case 2. The postpeak rock only consists of the BIDZ. As in equation (24), the radial displacement for the BIDZ can be calculated by substituting equations (14a), (14b), (19a), (19b), and (22) into equation (21) and considering the radial displacement continuity condition at the elastic-BIDZ interface.(25)ur=1rβ1+νRdERdδb1ξRd+βRd1−ξr−δb2ξRd+βRd2−ξrrξRd+β−RdξRd+β+δb22+βr2+β−Rd2+β+δb31+βr1+β−Rd1+β+uRdRdrβ,where uRd=1+ν/Eσ0−σRdRd.The tunnel roof convergenceu0 can be determined by substituting r=R0 into equations (24) and (25) for various cases, respectively.


## 4.3. Theoretical Value and Comparative Analysis

The solution for a circular opening consists of a series of simple formulas, which obviate the need to solve a complex differential equation. The solution takes into consideration rock gravity and the BIDZ, which can produce a relatively accurate result. Hence, the new solution has a wider engineering application.Zhang et al. [18] devised a closed-form solution and analyzed the stress and displacement of a circular opening based on the USC. However, their solution did not account for the effect of rock gravity, which is inconsistent with engineering practice. Chen et al. [19] also used the USC and presented a closed-form solution of a deep circular tunnel based on a new four-stage constitutive model. However, their solution only considered the cohesion degradation while ignoring the attenuation of other rock parameters. The rock gravity and BIDZ effects were also not considered in Chen’s solution. However, these solutions ignore the intermediate principal stress. The aggregate effect of rock gravity, BIDZ, and intermediate principal stress is, however, rationally considered in the new closed-form solution presented herein. The closed-form solution for Case 2, in which γ=0, cRd=cr, φRd=φr, νRd=νr, and ERd=Er, is identical to Zhang’s solution (n=0) [18]. Park’s solution, based on the M-C failure criterion, can be derived using γ=0, cRd=cr, φRd=φr, νRd=ν, ERd=E, and b=0 [6]. Zhang and Park’s solutions are special cases in this paper. Hence, this new solution includes a series of traditional results and is widely used in rock mass.


## 4.4. Further Analysis on Modified USC

### 4.4.1. Numerical Verification

Cheng [30] embedded the ideal elastoplastic constitutive model into ABAQUS finite-element software based on the USC and then simulated the stress distribution around a circular tunnel under conditions of various b values. Table 1 lists the geometry and rock parameters in this simulation. The associated flow rule, that is, ψ=φ, was adopted to determine the dilatancy coefficient β. However, the numerical simulation does not consider the effect of intermediate principal stress on the dilatancy coefficient. Therefore, β=1+sinψ/1−sinψ is used in this analytical solution. Since the damaged zone is ignored in Pan’s numerical simulation, Case 2 is used to facilitate a comparison with Pan’s simulation results. The analysis results on Figure 4 are highly consistent with Pan’s simulation results, which further confirm the accuracy of the new solution.Table 1
The geometry and rock parameters by Pan et al. [30].
Geometry and rock parametersValueIn situ stress,σ0 (MPa)30Tunnel radius,R0 (m)10Rock gravity,γ (kN/m3)0Internal pressure,σi (MPa)0Dilatancy angle,ψ(°)30Initial elastic modulus,E (GPa)2Initial cohesive,c (MPa)6Initial friction angle,φ (°)30Initial Poisson’s ratio (ν)0.3Damage elastic modulus,ERd (GPa)2Damage cohesive,cRd (MPa)6Damage friction angle,φRd (°)30Damage Poisson’s ratio (νRd)0.3Figure 4
Comparison of stress distribution under (a)b = 0.0, (b) b = 0.5, and (c) b = 1.0.
(a)(b)(c)


### 4.4.2. Parametric Study

A 3 m radius transportation roadway of Quandian coal mine in China is buried at about 650 m, with 15 MPa average in situ stress. The surrounding rock is mainly sandy mudstone. The special rock parameters of the roadway are shown in Table2. It should be noted that the 0.2 times of the initial strength parameters were used as the damage parameters, reflecting a serious damage phenomenon induced by blasting excavation.Table 2
The geometry and rock parameters used in this study.
Geometry and rock parametersValueIn situ stress,σ0 (MPa)15Tunnel radius,R0 (m)3Rock gravity,γ (kN/m3)27.5Dilatancy angle,ψ (°)0Initial elastic modulus,E (GPa)2.2Initial cohesive,c (MPa)5Initial friction angle,φ (°)30Initial Poisson’s ratio (ν)0.25Residual elastic modulus,Er (GPa)1.6Residual cohesive,cr (MPa)3.6Residual friction angle,φr (°)24Residual Poisson’s ratio (νr)0.25Damage elastic modulus,ERd (GPa)0.44Damage cohesive,cRd (MPa)1.0Damage friction angle,φRd (°)6Damage Poisson’s ratio (νRd)0.25


### 4.4.3. The Effect of Intermediate Principal Stress and Gravity on Ground Response Curves

Case 2 is used here to investigate the effect of intermediate principal stress and gravity on the ground response curve (GRC). Figure5 shows the GRC of the tunnel ignoring (Figure 5(a)) and considering (Figure 5(b)) the gravity effect under various intermediate principal stresses. The main research results are as follows.Figure 5
GRC at different gravity condition and intermediate principal stress. (a) Ignoring the gravity effect. (b) Considering the gravity effect.
(a)(b)As shown in Table3, with increasing b value, the tunnel surface convergence gradually decreases. When the b value increases from 0 to 0.5 and 1, the tunnel surface convergence decreases by 71.97% and 85.57% with the consideration of the gravity as well as by 68.93% and 83.48% with ignoring the gravity, respectively. The above analysis reveals that high intermediate principal stress can effectively improve the bearing capacity of surrounding rock and restrain tunnel convergence.Table 3
Roof surface displacement of the tunnel under various gravity conditions.
ParametersConditionsb0.00.250.50.751.0u0 (m)Gravity included1.0880.5070.3050.2100.157Gravity ignored0.9140.4580.2840.1990.151Considering the gravity condition, the roof surface displacement of the tunnel increases by 19.04% forb = 0, 7.39% for b = 0.5, and 3.79% for b = 1.0 compared with ignoring the gravity, respectively. This means that the gravitational condition can accelerate tunnel convergence, which is not conducive to long-term tunnel stability. Thus, the tunnel’s support strength design should fully consider the rock gravity at the roof. Moreover, with increasing intermediate principal stress, the rock gravity effect at the roof gradually weakens, again indicating that intermediate principal stress is beneficial to tunnel stability.With increasingb value, the critical pressure at the elastic-damage zone interface gradually decreases, demonstrating that the higher the intermediate principal stress, the smaller the support resistance needed to maintain tunnel stability.


### 4.4.4. The Effect of Gravity and Intermediate Principal Stress on Damage Zone Radius

Figure6 shows the trend of damage zone radius against the internal pressure at various intermediate principal stress and gravity conditions. The damage zone radii are listed in Table 4, from which the following can be seen.Figure 6
Evolution law of damage zone at various gravity conditions and intermediate principal stress.Table 4
Damage zone radius of the tunnel under various gravity conditions.
ParametersConditionsb0.00.250.50.751.0Rd/R0Gravity included3.7262.6142.0961.8031.616Gravity ignored3.4272.4962.0341.7651.591For the damage zone, as the internal pressure decreases gradually, the dimensionless valueRd/R0 nonlinearly increases at various b values, while the increase rates vary with the b value. The slope of the “σi vs. Rd/R0” curve for b = 0 is obviously less than that b = 1 at any gravity condition, which indicates that the extension process of the damage zone is closely dependent on the intermediate principal stress.Forσi=0MPa, with increasing the b value, the radius of the damage zone gradually decreases. When the b value increases from 0 to 0.5 and 1, the dimensionless value Rd/R0 decreases by 43.75% and 56.63% with the consideration of the gravity as well as by 40.65% and 53.57% with ignoring the gravity effect, indicating that the higher the intermediate principal stress, the smaller the damage zone radius. In practice, the damage zone is the main support area for the deep tunnel due to the lower rock strength. One of the most important functions of the bolt is to reinforce the damaged block and improve its residual strength. The support parameter design of bolts is dependent on the range of damage zones. At high intermediate principal stress, the range of damage zone is obviously smaller than that at lower intermediate principal stress. It means that we can optimize bolt support parameters by the methods of reducing its length and increasing its pretension for the tunnel under the high intermediate principal stress state to save its support cost.The roof rock gravity has an essential influence on the radius of damage zone. Considering the gravity effect, the dimensionless value ofRd/R0 increases by 8.72% for b = 0, by 3.05% for b = 0.5, and by 1.57% for b = 1.0 compared with ignoring the gravity, respectively. It means that the range of damage zone under the effect of gravity is larger than that without the influence of gravity. Hence, the roof support design of the tunnel should consider the effect of roof rock gravity for maintaining the long-term stability of the tunnel.


### 4.4.5. The Effect of Damage Zone Range on the Stress and Displacement

The stress and deformation of the tunnel are dependent on the damage zone range. In this work, Case 1 is used to explore the effect of damage zone radius on the stress and deformation of the tunnel while ignoring the gravity effect (γ=0). Figure 7 plots the trend of hoop stress and radial stress against the damage zone radius. With increasing the damage zone radius, the location of the maximum hoop stress dramatically shifted away from the surface of the tunnel, while its value has only lightly changed, indicating that the range of damage zones has an essential influence on the stress state of the circular tunnel. Therefore, the larger the radius of the damage zone, the larger the looseness range of the tunnel and the higher the instability risk of the tunnel.Figure 7
Hoop stress and radial stress of tunnel under various damage ranges (γ=0 and b=0). (a) Hoop stress. (b) Radial stress.
(a)(b)Figure8 shows the relationship between radial displacement and radius. According to the figure, the extent of the damage zone has a profound influence on the large deformation instability of the tunnel. The original radial displacement of the tunnel is always at a lower level of the plastic zone-damage zone interface regardless of the existence of any damage zone range. Once the rock mass at the roof enters the damage zone, the radial displacement of the tunnel sharply increases in a nonlinear manner, especially under Rd=9m. Furthermore, with increasing damage zone radius, the surface displacement of the tunnel observably increases, indicating that a large damage zone range is not conducive to tunnel stability. This is because the larger the damage zone range, the larger the surface displacement of the tunnel and then the larger the deformation pressure acting on the support structure, which leads to the collapse of the support structure due to the overload. Hence, the support methods of the tunnel should fully consider the effect of BIDZ range on the ground response.Figure 8
Displacement change of tunnel under various damage ranges (γ=0 and b=0).


## 4.4.1. Numerical Verification

Cheng [30] embedded the ideal elastoplastic constitutive model into ABAQUS finite-element software based on the USC and then simulated the stress distribution around a circular tunnel under conditions of various b values. Table 1 lists the geometry and rock parameters in this simulation. The associated flow rule, that is, ψ=φ, was adopted to determine the dilatancy coefficient β. However, the numerical simulation does not consider the effect of intermediate principal stress on the dilatancy coefficient. Therefore, β=1+sinψ/1−sinψ is used in this analytical solution. Since the damaged zone is ignored in Pan’s numerical simulation, Case 2 is used to facilitate a comparison with Pan’s simulation results. The analysis results on Figure 4 are highly consistent with Pan’s simulation results, which further confirm the accuracy of the new solution.Table 1
The geometry and rock parameters by Pan et al. [30].
Geometry and rock parametersValueIn situ stress,σ0 (MPa)30Tunnel radius,R0 (m)10Rock gravity,γ (kN/m3)0Internal pressure,σi (MPa)0Dilatancy angle,ψ(°)30Initial elastic modulus,E (GPa)2Initial cohesive,c (MPa)6Initial friction angle,φ (°)30Initial Poisson’s ratio (ν)0.3Damage elastic modulus,ERd (GPa)2Damage cohesive,cRd (MPa)6Damage friction angle,φRd (°)30Damage Poisson’s ratio (νRd)0.3Figure 4
Comparison of stress distribution under (a)b = 0.0, (b) b = 0.5, and (c) b = 1.0.
(a)(b)(c)


## 4.4.2. Parametric Study

A 3 m radius transportation roadway of Quandian coal mine in China is buried at about 650 m, with 15 MPa average in situ stress. The surrounding rock is mainly sandy mudstone. The special rock parameters of the roadway are shown in Table2. It should be noted that the 0.2 times of the initial strength parameters were used as the damage parameters, reflecting a serious damage phenomenon induced by blasting excavation.Table 2
The geometry and rock parameters used in this study.
Geometry and rock parametersValueIn situ stress,σ0 (MPa)15Tunnel radius,R0 (m)3Rock gravity,γ (kN/m3)27.5Dilatancy angle,ψ (°)0Initial elastic modulus,E (GPa)2.2Initial cohesive,c (MPa)5Initial friction angle,φ (°)30Initial Poisson’s ratio (ν)0.25Residual elastic modulus,Er (GPa)1.6Residual cohesive,cr (MPa)3.6Residual friction angle,φr (°)24Residual Poisson’s ratio (νr)0.25Damage elastic modulus,ERd (GPa)0.44Damage cohesive,cRd (MPa)1.0Damage friction angle,φRd (°)6Damage Poisson’s ratio (νRd)0.25


## 4.4.3. The Effect of Intermediate Principal Stress and Gravity on Ground Response Curves

Case 2 is used here to investigate the effect of intermediate principal stress and gravity on the ground response curve (GRC). Figure5 shows the GRC of the tunnel ignoring (Figure 5(a)) and considering (Figure 5(b)) the gravity effect under various intermediate principal stresses. The main research results are as follows.Figure 5
GRC at different gravity condition and intermediate principal stress. (a) Ignoring the gravity effect. (b) Considering the gravity effect.
(a)(b)As shown in Table3, with increasing b value, the tunnel surface convergence gradually decreases. When the b value increases from 0 to 0.5 and 1, the tunnel surface convergence decreases by 71.97% and 85.57% with the consideration of the gravity as well as by 68.93% and 83.48% with ignoring the gravity, respectively. The above analysis reveals that high intermediate principal stress can effectively improve the bearing capacity of surrounding rock and restrain tunnel convergence.Table 3
Roof surface displacement of the tunnel under various gravity conditions.
ParametersConditionsb0.00.250.50.751.0u0 (m)Gravity included1.0880.5070.3050.2100.157Gravity ignored0.9140.4580.2840.1990.151Considering the gravity condition, the roof surface displacement of the tunnel increases by 19.04% forb = 0, 7.39% for b = 0.5, and 3.79% for b = 1.0 compared with ignoring the gravity, respectively. This means that the gravitational condition can accelerate tunnel convergence, which is not conducive to long-term tunnel stability. Thus, the tunnel’s support strength design should fully consider the rock gravity at the roof. Moreover, with increasing intermediate principal stress, the rock gravity effect at the roof gradually weakens, again indicating that intermediate principal stress is beneficial to tunnel stability.With increasingb value, the critical pressure at the elastic-damage zone interface gradually decreases, demonstrating that the higher the intermediate principal stress, the smaller the support resistance needed to maintain tunnel stability.


## 4.4.4. The Effect of Gravity and Intermediate Principal Stress on Damage Zone Radius

Figure6 shows the trend of damage zone radius against the internal pressure at various intermediate principal stress and gravity conditions. The damage zone radii are listed in Table 4, from which the following can be seen.Figure 6
Evolution law of damage zone at various gravity conditions and intermediate principal stress.Table 4
Damage zone radius of the tunnel under various gravity conditions.
ParametersConditionsb0.00.250.50.751.0Rd/R0Gravity included3.7262.6142.0961.8031.616Gravity ignored3.4272.4962.0341.7651.591For the damage zone, as the internal pressure decreases gradually, the dimensionless valueRd/R0 nonlinearly increases at various b values, while the increase rates vary with the b value. The slope of the “σi vs. Rd/R0” curve for b = 0 is obviously less than that b = 1 at any gravity condition, which indicates that the extension process of the damage zone is closely dependent on the intermediate principal stress.Forσi=0MPa, with increasing the b value, the radius of the damage zone gradually decreases. When the b value increases from 0 to 0.5 and 1, the dimensionless value Rd/R0 decreases by 43.75% and 56.63% with the consideration of the gravity as well as by 40.65% and 53.57% with ignoring the gravity effect, indicating that the higher the intermediate principal stress, the smaller the damage zone radius. In practice, the damage zone is the main support area for the deep tunnel due to the lower rock strength. One of the most important functions of the bolt is to reinforce the damaged block and improve its residual strength. The support parameter design of bolts is dependent on the range of damage zones. At high intermediate principal stress, the range of damage zone is obviously smaller than that at lower intermediate principal stress. It means that we can optimize bolt support parameters by the methods of reducing its length and increasing its pretension for the tunnel under the high intermediate principal stress state to save its support cost.The roof rock gravity has an essential influence on the radius of damage zone. Considering the gravity effect, the dimensionless value ofRd/R0 increases by 8.72% for b = 0, by 3.05% for b = 0.5, and by 1.57% for b = 1.0 compared with ignoring the gravity, respectively. It means that the range of damage zone under the effect of gravity is larger than that without the influence of gravity. Hence, the roof support design of the tunnel should consider the effect of roof rock gravity for maintaining the long-term stability of the tunnel.


## 4.4.5. The Effect of Damage Zone Range on the Stress and Displacement

The stress and deformation of the tunnel are dependent on the damage zone range. In this work, Case 1 is used to explore the effect of damage zone radius on the stress and deformation of the tunnel while ignoring the gravity effect (γ=0). Figure 7 plots the trend of hoop stress and radial stress against the damage zone radius. With increasing the damage zone radius, the location of the maximum hoop stress dramatically shifted away from the surface of the tunnel, while its value has only lightly changed, indicating that the range of damage zones has an essential influence on the stress state of the circular tunnel. Therefore, the larger the radius of the damage zone, the larger the looseness range of the tunnel and the higher the instability risk of the tunnel.Figure 7
Hoop stress and radial stress of tunnel under various damage ranges (γ=0 and b=0). (a) Hoop stress. (b) Radial stress.
(a)(b)Figure8 shows the relationship between radial displacement and radius. According to the figure, the extent of the damage zone has a profound influence on the large deformation instability of the tunnel. The original radial displacement of the tunnel is always at a lower level of the plastic zone-damage zone interface regardless of the existence of any damage zone range. Once the rock mass at the roof enters the damage zone, the radial displacement of the tunnel sharply increases in a nonlinear manner, especially under Rd=9m. Furthermore, with increasing damage zone radius, the surface displacement of the tunnel observably increases, indicating that a large damage zone range is not conducive to tunnel stability. This is because the larger the damage zone range, the larger the surface displacement of the tunnel and then the larger the deformation pressure acting on the support structure, which leads to the collapse of the support structure due to the overload. Hence, the support methods of the tunnel should fully consider the effect of BIDZ range on the ground response.Figure 8
Displacement change of tunnel under various damage ranges (γ=0 and b=0).


## 5. Conclusions

A new elastic-plastic solution of deep circular tunnel excavated in elastic-brittle plastic rock mass is presented by using unified strength criterion. The new solution considers the effect of intermediate principal stress, rock gravity, and damage zone on ground response. Associated and nonassociated flow rules are also considered in the solution. The correctness of the solution is also verified by comparison with the numerical simulation results. This study demonstrates that the intermediate principal stress, rock gravity, and damage zone can dramatically influence the stress and deformation of surrounding rock. Some interesting meaningful conclusions are summarized as follows.The intermediate principal stress has a significant influence on the tunnel convergence, damage zone radius, and critical pressure at the interface of elastic-damage zone. With the increasing intermediate principal stress, the tunnel convergence, the radius of damage zone, and critical pressure decrease gradually under various gravity conditions, indicating that the intermediate principal stress is beneficial to the tunnel stability.The surface displacement and damage zone radius of the tunnel with the consideration of the gravity are obviously larger than those with ignoring the gravity, indicating that the rock gravity can accelerate roof instability of the tunnel. Furthermore, the rock gravity effect under the high intermediate principal stress gradually weakens, illustrating again that the intermediate principal stress is crucial for the tunnel stability.Large deformation instability of the tunnel is dependent on the extension of damage zone. Originally, both fracture range and tunnel convergence are in a lower level before damage zone. Once the rock mass enters the damage zone, both fracture range and tunnel convergence sharply and nonlinearly increase. The larger the radius of damage zone, the larger both fracture range and tunnel surface deformation. Hence, the support method design of the tunnel should fully consider the effect of the range of BIDZ on the ground response.The intermediate principal stress, rock gravity, and damage zone can significantly influence the ground response. Hence, it is suggested to perform parametric studies during the early design stage of the tunnel to evaluate the ground convergence for different scenarios.



---

*Source: 5537963-2021-04-24.xml*