#### 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.

#### 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]:where , , and are the principal shear stresses acting on the orthogonal octahedron element, defined as , and subscripts “” and “” denote 1, 2, and 3. , , and represent the normal stresses acting on the orthogonal octahedron element corresponding to the principal shear stress, defined as . , , and are the maximum, intermediate, and minimum principal stresses, respectively. The coefficients of and are material parameters that can be determined experimentally. In practice, the USC can be also expressed by the following Mohr-Coulomb strength parameters [25]:

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 is shown in equation (3), where , , and denote, respectively, uniaxial compressive strength, uniaxial tensile strength, and shear strength:

As shown in Figure 1, when the value of the parameter 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]:

#### 3. Tunnel Condition and Basic Equation

##### 3.1. Definition of the Problem

In Figure 2, the depicted circular tunnel of radius is in an elastic state following excavation. After blasting, a cylindrical BIDZ of radius 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 (, , , and , 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 , , , and ). 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.

**(a)**

**(b)**

**(a)**

**(b)**

The axial stress is the intermediate principal stress, which can be expressed as , and [18]. is the initial Poisson ratio of the rock mass. To simplify calculation, can be taken in the damage and plastic zones, and the relationship can be deduced. Thus, equation (4b) should be embraced and afterwards can be changed as follows:where

The failure parameters of damaged rock mass are denoted by subscript “*D*” (namely, , , , and ).

##### 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 aswhere is the dilation parameter associated with the coefficient of intermediate principal stress and the dilation angle . It is generally assumed that the plastic strain increments are proportional to the stress gradient of the plastic potential function [2]:where and are the circumferential and radial plastic strain increments, respectively; is a nonnegative proportional constant. Submitting equations (7), (8a), and (8b), the contact between the plastic strain and dilation parameter is obtained:

Most previous studies regarded 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]: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 Figure 2, 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,

In Case 2, the stresses and deformation of the tunnel in the elastic zone can be obtained only by substituting and by and in equations (11a)–(11c). and 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):

##### 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 at , the radial and circumferential stresses of the plastic zone can be expressed as follows:where , , and .

By considering the blast-induced damage effect in the same way, substituting equation (5b) into equation (10), and using the boundary condition at , the radial and circumferential stresses for the BIDZ can be obtained:where , , and .

The closed solution is acquired by substituting into equations (13a) and (14a) and then equating the right-hand sides of equations (13a) and (14a):

From the above equation, the radius of the BIDZ () is known, and the approximate solutions of 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 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 into equation (14a), followed by equating the right-hand sides of equations (12) and (14a), gives

The approximate solution for can also be calculated by using the method of bisection.

###### 4.2.2. Displacement Analysis

The total circumferential strain and radial strain at radius for the axial-symmetric condition can be expressed by the radial displacement [1]:

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:

By adopting Hooke’s law, the elastic strain components and in the fracture zone can be expressed as follows [1, 18]:

Substituting equations (9) and (17) into equations (19a) and (19b) gives

*(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,where .

Substituting equations (13a), (13b), (19a), and (19b) into equation , we get the following equation:where , , and are constants. , , , , , and .

Upon substituting equation 22) into equation (21),

As shown in equation (23), the radial displacement in the BIDZ can also be deduced:where , , and are constants. The radial displacement at the plastic-BIDZ interface is obtained by substituting into equation (23). , , , , and .

*(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.where .

The tunnel roof convergence can be determined by substituting 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 , , , , and , is identical to Zhang’s solution () [18]. Park’s solution, based on the M-C failure criterion, can be derived using , , , , , and [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, 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.

**(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 Table 2. 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.

###### 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). Figure 5 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.

**(a)**

**(b)**

As shown in Table 3, 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.

Considering the gravity condition, the roof surface displacement of the tunnel increases by 19.04% for *b* = 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 increasing *b* 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

Figure 6 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.

For the damage zone, as the internal pressure decreases gradually, the dimensionless value nonlinearly increases at various *b* values, while the increase rates vary with the *b* value. The slope of the “ vs. ” 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 , 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 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 of 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 (). 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.

**(a)**

**(b)**

Figure 8 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 . 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.

#### 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.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (nos. 51874277 and 51874278) and the Natural Science Foundation of Jiangsu Province of China (no. BK20181357).