# Simulating the Formation of Blasting-Excavation-Induced Zonal Integration in Deep Tunnels with an Elastoplastic Damage Model

Authors: Qiang Gao; Chuanxiao Liu; Jian Zhang; Guangtan Cheng
Journal: Shock and Vibration
 (2021)

---

## Abstract

More deep tunneling projects will be constructed due to the increasing demand of underground energy and resource. The zonal disintegration phenomena are frequently encountered with the surrounding rock of deep tunnels. To explain the mechanisms underlying the formation of zonal disintegration, an elastoplastic damage model and failure criterion are proposed in this study based on the strain gradient theory and the damage property of rock mass. A coupling calculation subroutine is thereafter developed by the ABAQUS code. The dynamic formation and development regularity of zonal disintegration in the deep tunnel are simulated by this subroutine. The radial displacement, radial stress, and tangential stress show the oscillated variation of peaks and troughs alternately. The coupling effect of the blasting load and the initial geostress transient unloading leads to the variation of alternation oscillation in the surrounding rock stress field, which is an important reason for the zonal disintegration of the surrounding rock. The morphological characteristics of fractured zones and nonfractured zones obtained from numerical simulations are in good agreement with the results from the in situ observations, which confirm the correctness and feasibility of the damage and numerical approach. The method proposed in the current study can be utilized to provide a basis for the prediction and supporting design of fractured modes.


---

## Body

## 1. Introduction

With rapid economy development and increasing demand for energy, many underground engineering applications, such as mining, traffic engineering, largescale water conservancy, and hydroelectric engineering, have entered a high-geostress and complex geological environment. The mechanical characteristics, deformation, and failure modes of deep rock masses differ from those of shallow rock masses. Zonal disintegration is a unique failure phenomenon in the deep surrounding rock, in which ruptured and nonruptured zones alternately appear in the surrounding rock of a deep tunnel [1].The zonal disintegration phenomenon has been confirmed through many physical detection methods in the excavation of deep tunnels. The periodic rupture phenomenon was discovered through spectroscopic observations in South African 2300 m gold mine [2]. A resistivity method was used to detect the zonal disintegration in the deep mine of Oktyabrskil (Figure 1) [3]. This phenomenon was observed by borehole television and an electrical resistivity survey in the Huainan coal mine in China (Figure 2) [4].Figure 1
Zonal disintegration of Oktyabrskil mine in Russia.Figure 2
Distribution of cracked zones in Huainan mine of China.The zonal disintegration phenomenon was reproduced in a model test. Then 2D and 3D equivalent material model tests were carried out and the zonal disintegration of surrounding rock was observed [5]. By using cement mortar as a similar material, the model test was performed and the alternate ruptured and nonruptured zones were generated in the surrounding rock of the model tunnel under high axial pressure [6]. A 3D triaxial loading geomechanical model was used to realistically reproduce the zonal disintegration phenomenon (Figure 3) in which the variation characteristics of strain and displacement were obtained using a variety of monitoring methods [7].Figure 3
Zone disintegration of the model test.For static numerical simulation, 3D realistic failure process analysis (RF-PA3D) was used to simulate the zonal disintegration phenomenon of the underground tunnel with the increase of axial load [8, 9]. Flac3D, a finite difference method software application, was implemented to conduct a numerical simulation of zonal disintegration in deep tunnels with a strain-softening joint model [10]. For dynamic numerical simulation, the unloading of tunnel excavation was regarded as a dynamic process. A corresponding dynamic model was set up and the Flac3D numerical code was used to simulate zonal disintegration in the surrounding deep tunnel rock [11–13]. Based on the continuous cap model suitable for dynamic properties of rock, LS-DYNA finite element software was used to simulate zonal disintegration, and it was proposed that static stress gradient and dynamic loading were the two prerequisites for the far-field fracture zone [14]. A constitutive model based on the general particle dynamics (GPD) method was proposed to investigate the zonal disintegration mechanism of isotropic rock masses around a deep circular tunnel subjected to dynamic unloading [15]. The numerical results indicated that the dynamic loads and high in situ stress are two dominant factors in the occurrence of zonal disintegration. The above dynamic load only refers to the in situ stress transient unloading and ignores the effect of blast loading on the rock mass surrounding the tunnel. There has been no lack of research on the coupling dynamic effect of the blast loading and the in situ stress transient release on the mechanical analysis of the zonal disintegration phenomenon.Zonal disintegration in a deep tunnel is completely different from that of a shallow cavity. It is difficult to explain zonal disintegration using elastic-plastic mechanics in a continuum mechanics frame. Zonal disintegration is a special and regular strain localization phenomenon [16]. Since the strain gradient plays a controlling role in the strain localization of rock [17], its influence should be considered in the study of zonal disintegration.An elastoplastic damage model and failure criterion were proposed based on the strain gradient theory. The zonal disintegration calculation program was written based on the dynamic finite element equations and ABAQUS finite element software. The morphological characteristics of zonal disintegration under blasting excavation were obtained through numerical simulation. This study reveals the morphological characteristics of zonal disintegration, which is used as a basis for the prediction and support control of fractured modes.


## 2. Elastoplastic Damage Model Based on Strain Gradient

### 2.1. Introduction of High-Order Stress

In the Toupin–Mindlin strain gradient theory [18–20], the total strain ui has two parts: the conventional Eulerian strain tensor εij and the high-order strain tensor ηijk, which are deduced as(1)εij=ui,j+uj,i2,ηijk=∂i∂juk=uk,ij.According to the second law of thermodynamics, the Helmholtz free energy function [18] can be expressed as(2)ψε,η=12εij:Eijklepd:εkl+12ηijk⋮Λijklmned⋮ηlmn,where σij is a second-order Cauchy stress tensor, which is the conjugate of the Eulerian strain tensor εij, and τijk is set as the third-order stress tensor which is conjugate to the strain gradient tensor ηijk. The stress tensor σij and high-order stress tensor τijk can be derived from the free energy function corresponding to the strain tensor εij and strain gradient tensor ηijk, respectively:(3)σij=∂ψ∂εij=Eijklepdεkl,τijk=∂ψ∂ηijk=Λijklmnedηlmn,where Eijklepd is the elastoplastic damage tensor, Eijklepd=Eijkle−Eijklpd; d is the damage tensor; Eijkle is the elastic tensor; Eijklp is the plastic tensor determined by plastic flow law and hardening law [21]; Λijklmned is the sixth-order elastic damage tensor considering the strain gradient, Λijklmned=1−dl2Eijlmeδkn; δknis the Kronecker symbol; and l is the internal length parameter of the material, which is closely related to the microcracks and microdefects inside the material, whose value is generally the average aggregate particle diameter of the corresponding material [22].


### 2.2. Establishment of Elastoplastic Damage Constitutive Equation

Cauchy stressσijand Eulerian strain εij are expressed as an incremental matrix:(4)dσ=1−ω+δijω3Hdε,where dσ=dσxxdσyydσzzdσxydσxzdσyzT and dε=dσxxdσyydσzzdσxydσxzdσyzT.H is 6×6 matrixH=He−Hp, and He can be written as(5)He=E1−υ1+υ1−2υEυ1+υ1−2υEυ1+υ1−2υEυ1+υ1−2υE1+υ1+υ1−2υEυ1+υ1−2υ000Eυ1+υ1−2υEυ1+υ1−2υE1−υ1+υ1−2υ000000E21+υ000000E21+υ000000E21+υ,where E is the elastic modulus and υ is Poisson’s ratio. Hpis expressed as (6)Hp=He∂Q/∂σ∂F/∂σHe∂F/∂σHe∂Q/∂σ−∂F/∂σσ∂Q/∂σ,where F=Fσ,κ=0 is the yield function,κ is plastic work, dκ=σdεp, andQ is plastic potential function.The high-order stress tensorτijk and strain gradient tensor ηijk are expressed as an incremental matrix:(7)dτ=1−ωΛdη,where(8)dτ=dτxxxdτyyxdτzzxdτxyxdτyzxdτxzxdτxxydτyyydτzzydτxyydτyzydτxzydτxxzdτyyzdτzzzdτxyzdτyzzdτxzzT,dη=dηxxxdηyyxdηzzxdηxyxdηyzxdηxzxdηxxydηyyydηzzydηxyydηyzydηxzydηxxzdηyyzdηzzzdηxyzdηyzzdηxzzT.Λ is 18×18 matrix,Λ=Εl2/81+υD1D2D3D4, and D1,D2,D3,D4 can be written as(9)D1=21660000006114000000641100000000010002310000700000000010000000400116400060062160002004611,D2=600000006400000002200000004001000100000402000020213000000000020000000060000000040,D3=321000000000002000000010000000004000000004000000006000000010000000200132312000000,D4=1000000002010012300000700010002011460000304116000060662100000100070000000001002000000010.Thus, the elastoplastic damage constitutive equation based on the strain gradient can be expressed as(10)dσdτ=Ddεdη=1−ω+δij3H001−ωΛdεdη.


### 2.3. Damage Evolution Equation

It is assumed that the damage of a rock mass can be represented by an isotropic damage variabled with the use of analytical methods from continuous damage mechanics. Since the damage variable d is an internal variable, equation (10) cannot form a complete damage constitutive model. The damage criterion and damage evolution law should also be determined.According to the mechanical properties of deep rock under high stress, the stress-strain relationship is simplified properly, as shown in Figure4(a). The OA stage is the linear elastic stage, in which no damage occurs to the rock mass. εs is the yield stress and σs is the yield strain. The AB stage is the plastic damage stage. The original cracks are pressed, and new cracks are produced after the rock mass enters the plastic stage. εu is the ultimate strain. Once the deformation reaches the ultimate strain, the rock mass is completely destroyed. To describe the strain-softening of deep rock under high stress, the damage evolution of rock is determined as [23](11)d=dε˜=0,ε<εs,εuε˜ε˜−εsεu−εs,εs≤ε<εu,1,ε≥εu,where ε˜ is the generalized equivalent strain containing the strain gradient [24], ε˜=2/3εijεij+l2ηijkηijk. The variation curve of damage variable d with equivalent strain ε˜ is shown in Figure 4(b).Figure 4
Nonlinear damage evolution law of rock under high stress. (a)σ−ε curve. (b)d−ε curve.
(a)(b)Thus, equations (13) and (14) constitute an elastoplastic damage model based on the strain gradient.


## 2.1. Introduction of High-Order Stress

In the Toupin–Mindlin strain gradient theory [18–20], the total strain ui has two parts: the conventional Eulerian strain tensor εij and the high-order strain tensor ηijk, which are deduced as(1)εij=ui,j+uj,i2,ηijk=∂i∂juk=uk,ij.According to the second law of thermodynamics, the Helmholtz free energy function [18] can be expressed as(2)ψε,η=12εij:Eijklepd:εkl+12ηijk⋮Λijklmned⋮ηlmn,where σij is a second-order Cauchy stress tensor, which is the conjugate of the Eulerian strain tensor εij, and τijk is set as the third-order stress tensor which is conjugate to the strain gradient tensor ηijk. The stress tensor σij and high-order stress tensor τijk can be derived from the free energy function corresponding to the strain tensor εij and strain gradient tensor ηijk, respectively:(3)σij=∂ψ∂εij=Eijklepdεkl,τijk=∂ψ∂ηijk=Λijklmnedηlmn,where Eijklepd is the elastoplastic damage tensor, Eijklepd=Eijkle−Eijklpd; d is the damage tensor; Eijkle is the elastic tensor; Eijklp is the plastic tensor determined by plastic flow law and hardening law [21]; Λijklmned is the sixth-order elastic damage tensor considering the strain gradient, Λijklmned=1−dl2Eijlmeδkn; δknis the Kronecker symbol; and l is the internal length parameter of the material, which is closely related to the microcracks and microdefects inside the material, whose value is generally the average aggregate particle diameter of the corresponding material [22].


## 2.2. Establishment of Elastoplastic Damage Constitutive Equation

Cauchy stressσijand Eulerian strain εij are expressed as an incremental matrix:(4)dσ=1−ω+δijω3Hdε,where dσ=dσxxdσyydσzzdσxydσxzdσyzT and dε=dσxxdσyydσzzdσxydσxzdσyzT.H is 6×6 matrixH=He−Hp, and He can be written as(5)He=E1−υ1+υ1−2υEυ1+υ1−2υEυ1+υ1−2υEυ1+υ1−2υE1+υ1+υ1−2υEυ1+υ1−2υ000Eυ1+υ1−2υEυ1+υ1−2υE1−υ1+υ1−2υ000000E21+υ000000E21+υ000000E21+υ,where E is the elastic modulus and υ is Poisson’s ratio. Hpis expressed as (6)Hp=He∂Q/∂σ∂F/∂σHe∂F/∂σHe∂Q/∂σ−∂F/∂σσ∂Q/∂σ,where F=Fσ,κ=0 is the yield function,κ is plastic work, dκ=σdεp, andQ is plastic potential function.The high-order stress tensorτijk and strain gradient tensor ηijk are expressed as an incremental matrix:(7)dτ=1−ωΛdη,where(8)dτ=dτxxxdτyyxdτzzxdτxyxdτyzxdτxzxdτxxydτyyydτzzydτxyydτyzydτxzydτxxzdτyyzdτzzzdτxyzdτyzzdτxzzT,dη=dηxxxdηyyxdηzzxdηxyxdηyzxdηxzxdηxxydηyyydηzzydηxyydηyzydηxzydηxxzdηyyzdηzzzdηxyzdηyzzdηxzzT.Λ is 18×18 matrix,Λ=Εl2/81+υD1D2D3D4, and D1,D2,D3,D4 can be written as(9)D1=21660000006114000000641100000000010002310000700000000010000000400116400060062160002004611,D2=600000006400000002200000004001000100000402000020213000000000020000000060000000040,D3=321000000000002000000010000000004000000004000000006000000010000000200132312000000,D4=1000000002010012300000700010002011460000304116000060662100000100070000000001002000000010.Thus, the elastoplastic damage constitutive equation based on the strain gradient can be expressed as(10)dσdτ=Ddεdη=1−ω+δij3H001−ωΛdεdη.


## 2.3. Damage Evolution Equation

It is assumed that the damage of a rock mass can be represented by an isotropic damage variabled with the use of analytical methods from continuous damage mechanics. Since the damage variable d is an internal variable, equation (10) cannot form a complete damage constitutive model. The damage criterion and damage evolution law should also be determined.According to the mechanical properties of deep rock under high stress, the stress-strain relationship is simplified properly, as shown in Figure4(a). The OA stage is the linear elastic stage, in which no damage occurs to the rock mass. εs is the yield stress and σs is the yield strain. The AB stage is the plastic damage stage. The original cracks are pressed, and new cracks are produced after the rock mass enters the plastic stage. εu is the ultimate strain. Once the deformation reaches the ultimate strain, the rock mass is completely destroyed. To describe the strain-softening of deep rock under high stress, the damage evolution of rock is determined as [23](11)d=dε˜=0,ε<εs,εuε˜ε˜−εsεu−εs,εs≤ε<εu,1,ε≥εu,where ε˜ is the generalized equivalent strain containing the strain gradient [24], ε˜=2/3εijεij+l2ηijkηijk. The variation curve of damage variable d with equivalent strain ε˜ is shown in Figure 4(b).Figure 4
Nonlinear damage evolution law of rock under high stress. (a)σ−ε curve. (b)d−ε curve.
(a)(b)Thus, equations (13) and (14) constitute an elastoplastic damage model based on the strain gradient.


## 3. Implementation of the Model of a Deep Tunnel under Dynamic Excavation

### 3.1. Construction ofC1Order Element considering Strain Gradient

To implement the elastoplastic damage model with the finite element method, the influence of the strain gradient should be considered first. The interpolating shape functions of a conventional finite element generally have only first-order continuity, and the constructed elements have onlyC0order continuity. The second derivative after the interpolation of the displacement field is zero, and the influence of the strain gradient cannot be considered in such elements. Therefore, it is necessary to consider finite elements with higher-order continuity, such as elements with C1order continuity. This requires that the normal displacement and first-order partial differential maintain continuity so that the strain gradient can enter into the finite element equation.Taking the most versatile tetrahedral element as an example, a three-dimensional general tetrahedral element withC1order continuity is established in Figure 5. The geometric coordinates of the node i (i = 1; 2; 3; 4) are xi,yi,zi. Four-node displacements can be assumed as(12)a=a1a2a3a4T,ai=uu,xu,yu,zvv,xv,yv,zww,xw,yw,ziT.Figure 5
Tetrahedral elements withC1order continuity.The area coordinatesLi of any point Px,y,z within the element can be defined as(13)Li=AiV,i=1,2,3,4,where A1=VolP234A2=VolP134A3=VolP214A4=VolP231, and V=Vol1234.Then, the above formula can be expanded as follows:(14)Li=−1i+1ai+bix+cty+diz6v,where(15)V=161x1y1z11x2y2z21x3y3z31x4y4z4,ai=xjyjzjxkykzkxlylzl,bi=−1yjzj1ykzk1ylzl,ci=1xjzj1xkzk1xlzl,d=−1xjyj1xkyk1xlyl,where i,j,k,andl are in accordance with the rotation of 1⟶2⟶3⟶4⟶1 and meet the right-handed spiral rule.The displacement fieldu can be expressed as(16)u=∑i=14αiLi3+α5L12L2+α6L12L3+α7L12L4+α8L22L3+α9L22L4+α10L22L1+α11L32L4+α12L32L1+α13L32L2+α14L42L1+α15L42L2+α16L42L3.The values of the 16 coefficientsαii=1,…16can be determined from the displacement ui1=1,…4of the four vertices and the gradient ui,x,ui,y,ui,zi=1,…4 of the three directions.After introducing a shape functionN, the displacement field in three dimensions is expressed as(17)u=uvw=N1000N1000N1|N2000N2000N2|N3000N3000N3|N4000N4000N4a1a2a3a4=Na,where(18)Ni=Ni0NixNiyNiz=Li3+3Li2AiAjLj−AkLk+AlLl6VAiLi2pklLj+pljLk+pjkLl6VAiLi2qklLj+qljLk+qjkLl6VAiLi2rklLj+rljLk+rjkLlT,A1=b2b3b4c2c3c4d2d3d4,A2=b1b3b4c1c3c4d1d3d4,A3=b1b2b4c1c2c4d1d2d4,A4=b1b2b3c1c2c3d1d2d3,pjk=ckdj−cjdk,qjk=bjdk−bkdj,rjk=bkcj−bjck,01×4=0000.Strainε and strain gradients η take the following form:(19)ε=B1a=L1Na,η=B2a=L2Na,where L1 and L2 are differential operators that can be written as(20)L1=∂∂x∂∂y∂∂z∂∂y∂∂x∂∂z∂∂z∂∂y∂∂x,L2=M000M000M,where M=∂2/∂x2∂2/∂y2∂2/∂z2∂2/∂x∂y2∂2/∂y∂z2∂2/∂x∂z2T,01×6=000000T.Therefore, the geometric equation considering the strain gradient can be expressed as(21)εη=Ba=B1B2a,where B is the strain matrix.Finally, the element stiffness matrix of theC1order tetrahedral element is(22)K=∭BTDBdxdydz.


### 3.2. Establishment of the Dynamic Finite Element Equation

Similar to static finite element analysis, dynamic finite element analysis should establish the finite element equation. The space of the rock mass is discretized, and higher-order tetrahedral elements are constructed. The shape function and stiffness matrix of the elements are deduced. According to the dynamic equilibrium equations of each node, the Newmark-β integral method [25, 26] is adopted to solve the following motion equation:(23)Mu¨+Cu˙+Ku=P,where M is the element mass matrix; M=∫VρNTNdV,ρ is the density of the rock mass; C is the element damping matrix, C=∫VcNTNdV, c is damping coefficient; K is the element stiffness matrix; u¨,u˙, and u are the acceleration, velocity, and displacement matrices at any point in the element; and P is the external loads.Displacement at any point in the element isu=Na, where N is a function only of position and is not related to time t, and a is related to time and can be expressed as a=a1ta2ta3ta4tT. Therefore, the acceleration u¨ and velocity u˙ at any point in the element are expressed, respectively, as u¨=Na¨ and u˙=Na˙ where · indicates a derivative with respect to time.In the Newmark-β method, when the control parameters satisfy α=1/2, β=1/4, the method changes to an average constant acceleration method with second-order accuracy, which satisfies the engineering requirements. Therefore, the basic assumption based on the method can be expressed as [27](24)u˙t1+Δt=u˙t1+Δt2u¨t1+u¨t1+Δt,ut1+Δt=ut1+u˙t1Δt+Δt24u¨t1+u¨t1+Δt.Equation (24) can be expanded as(25)u¨t1+Δt=4Δt2ut1+Δt−4Δt2ut1−4Δtu˙t1−u¨t1,u¨t1+Δt=u˙t1+Δt2u¨t1+u¨t1+Δt=2Δtut1+Δt−2Δtut1−u˙t1.After subsisting equation (25) into equation (23) at time t1+Δt, the expression can be arranged as(26)K¯ut1+Δt=P¯t1+Δt,where K¯ is the equivalent stiffness matrix, K¯=4/Δt2M+2/ΔtC+K;P¯t1+Δt is the equivalent load, P¯t1+Δt=Pt1+Δt+M4/Δt2ut1+4/Δtu¨t1+u¨t1+C2/Δtut1+u˙t1.In iterative computation, the process of each time step are as follows:(1)
Calculate the equivalent loadP¯t1+Δt according to the operation state and external load of the node.(2)
Calculate the equivalent stiffness matrixK¯ according to the node.(3)
Solve equation set (24). When the structure is in the elastic state, the equivalent stiffness matrix K¯ does not change with the displacement of the node, which is a linear equation set and can be solved directly by triangular decomposition. When the structure is in the plastic state, the equivalent stiffness matrix K¯ is related to the displacement of the node and must be solved through displacement iteration.(4)
Calculate the acceleration and velocity at timet1+Δt by equation (25).


### 3.3. Dynamic Incremental Variable Plastic Stiffness Iterative Method

When the structure enters the plastic stage described by the Mohr–Coulomb plasticity criterion, equation set (26) is nonlinear. The equivalent stiffness matrix K¯ is related to displacement and must be solved iteratively. The dynamic incremental variable plasticity stiffness iteration method is used to solve the equation set [28].The equivalent stiffness matrixK¯ of dynamic calculation is made up of mass matrix M, damping matrix C, and stiffness matrix K. The plastic stiffness matrix KP is formed from K when the element enters plasticity.To guarantee the convergence of iteration, according to the incremental plastic stiffness iteration method, the plastic load of each time step is applied in multiple stages. In each stage of plastic load application, the overall stiffness matrix of the structure can be decomposed into two parts: elastic and plastic:(27)K=∭vBTDe−DpBdv=Ke−Kp.In the calculation, the elastic stiffness matrixKe remains unchanged, and the plastic stiffness matrix Kp is changed according to the critical stress state imposed after each stage load. In the process of each stage of the plastic load iteration, according to the basic idea of the incremental load method, the plastic stiffness matrix is kept constant, and the displacement δt1+Δt0 solved by the elastic equivalent stiffness matrix K¯e is taken as the initial displacement of the iterative calculation. The displacement value δt1+Δt under the incremental load at each stage is solved by the continuous correction of the displacement. The basic iterative equation for the dynamic elastoplastic finite element method is(28)K¯eut1+Δt=ΔP¯p+Kput1+Δt,K¯e=4Δt2M+2ΔtC+Ke,where K¯e is the elastic equivalent stiffness matrix.


### 3.4. Calculation Method of External Load

#### 3.4.1. Stress Decomposition

Considering the complexity of deep underground engineering, a circular tunnel is used to carry out the mechanical analysis. The blasting excavation process of a deep tunnel is accompanied by the blasting load and in situ stress transient unloading, and blasting excavation is the precondition for transient unloading of in situ stress [29]. Thus, on the one hand, the stress wave produced by the blasting load and detonation gas will cause blasting damage of surrounding rock within a certain range; and on the other hand, transient unloading of stress in a deep tunnel will cause stress redistribution of surrounding rock. Therefore, the stress field of blasting excavation in Figure 6 is decomposed into two problems:(1)
The elastoplastic stress field caused by the blasting loadPbt(2)
The elastoplastic stress field caused by the transient unloading of in situ stressPotFigure 6
Forces on circular tunnel under dynamic unloading.


#### 3.4.2. Blasting Load

The blasting load is actually the shock wave pressure acting on the surrounding rock after the explosion. In this paper, the change process of the blasting load is described by a triangle function (Figure7). Pb is the equivalent peak load of the blasting load on the excavation surface and tr and tdare the time of the load rising and the positive pressure, respectively.Figure 7
Explosion load curve.Taking full section blasting in a tunnel as the research object, the influence of the blasting load on rock mass is taken into account.The detonation CJ model [30, 31] is used to calculate the peak load of the blasting load. For cylindrical radial uncoupled loading conditions, the equivalent peak load Pb on the excavation surface is calculated as(29)Pb=ρ0De221+γdcdb2γdbS,where ρ0 is the density of explosives, De is the detonation velocity of explosives, γ is the isentropic exponent of explosives, which is generally taken as 3, dcis the diameter of the cartridge, dbis the diameter of the blast hole, and Sis the distance of adjacent blast holes.When the explosion wave is propagating to the blasting hole with velocityDe, the blasting load in the hole rises to the maximum, and the rising time tr can be written as(30)tr=L1D.The blasting load is a transient load that changes with time during blasting. With explosive gas escaping, the blasting load decreases further. When the pressure of the blast hole drops to atmospheric pressure, the synchronous unloading of the excavation load is completed. The duration of the blasting loadtd is calculated as(31)td=L1D+1/4S2+L22Cf+L1+L2Cu1+L1+L2Cu2,where L1 and L2 are the respective lengths of the charging section and blockage section in the blast hole; Cf is the average speed of expansion driven by crack explosion loading; Cu1is the unloading wave speed of gas explosion; and Cu2 is the reflected wave velocity spread from the bottom to the top of the blast hole [32].The time history relationPbt of the blasting load acting on the excavation surface can be expressed as follows:(32)Pbt=0,t≤0ort≥td,Pbttr0<t≤tr,Pbtd−ttd−trtr<t<td.


#### 3.4.3. Transient Unloading of In Situ Stress

The unloading of in situ stress on the excavation surface [33] should satisfy the continuous stress condition, so the initial time and duration of transient unloading are determined by the blasting process. The in situ stress unloads at time tb, defined as the moment when the blasting load decreases to the in situ stress. When the blasting load is further reduced to zero, the synchronous unloading of in situ stress is completed. td−tb is the duration of transient unloading of in situ stress. The transient unloading curve of in situ stress coincides with the curve of the blasting load (Figure 8).Figure 8
Coupling effect curves of blasting load and transient release of in situ stress.The time history relation of the in situ stressP0t on the excavation surface can be written as(33)P0t=P0,t≤tb,−P0t−tdtd−tb,tb<t<td,0,t≥td.


#### 3.4.4. Loading Method

In the process of blasting excavation of an underground tunnel, the external loadsP imposed on surrounding rock elements are mainly divided into blasting load Fp, gravity loads Fg, and excavation loads Fe [34].According to the stress generated by blasting on rock mass elements, the blasting loadFp can be obtained as(34)Fp=∭vBTPbtdv,where Fg is the vertical gravity loads caused by the removed elements and can be expressed as(35)Fg=∭vNTρgNdv,where ρ is the density of the rock mass, g is the acceleration of gravity, and N is the shape function matrix which can be obtained by Hermit interpolation function. Fe, the transient unloading of in situ stress, can be expressed as(36)Fe=∭vBTP0tdv,where B is strain matrix and P0t is the in situ stress of the excavation unit, obtained from equation (33).


### 3.5. Zonal Disintegration Damage Failure Criterion

#### 3.5.1. Maximum Tensile Strain Criterion

Tensile failure of surrounding rock is generally the result of strain development in the direction perpendicular to the surface of the wall. The maximum tensile strain criterion can be written as(37)f=εtmax−εtu=0,where εtmax is the maximum tensile strain of a calculation element [35] and εtu is the ultimate tensile strain, which is obtained by a uniaxial tensile test or Brazilian disk test.Iff≥0, then tensile failure will occur and the element is destroyed. To maintain the integrity and continuity of the entire calculation during the numerical simulation, a small residual elastic modulus EC is given to the element where the tensile failure occurs; usually, EC=0.05E. If f < 0, the tensile failure will not occur. The zonal disintegration energy damage failure criterion is used as the criterion of element failure.


#### 3.5.2. Zonal Disintegration Energy Damage Failure Criterion

Considering the influence of the strain gradient term, the formula for calculating the strain energy density of an element is written as(38)dWdV=∫0εijσijdεij+18Gl24ηijjηikk+4ηiikηkjj+4ηiikηjjk+7ηijkηijk+2ηijkηkji.Figure9 shows the simplified stress-strain relation curve of rock under high stress. Elastic strain energy density dW/dVs, critical strain energy density dW/dVo, and ultimate strain energy density dW/dVu can be calculated as [36](39)dWdVe=∫0εeσijdεij,dWdVo=∫0εoσijdεij,dWdVu=∫0εuσijdεij.Figure 9
Stress-strain relation curve of rock under high stress.The damage failure criterion of zonal disintegration is judged by the change of the strain energy density.(1)
WhendW/dV<dW/dVs=∫0εsσijdεij (εsis the yield strain in the loading process), then the element is in the linear elastic stage, and damage within the element will not occur.(2)
WhendW/dV≥dW/dVs, the element enters into the plastic damage stage; microcracks in the rock mass appear and begin to expand.(3)
WhendW/dV≥dW/dVD=∫0εoσijdεij (εois the critical strain in the loading process), the microcracks in the rock mass develop into a macroscopic crack, and there is still bearing capacity in the rock body. This is called the plastic residual stage.(4)
WhendW/dV≥dW/dVu=∫0εuσijdεij (εu is the ultimate strain in the loading process), the element is damaged to failure. In the case of tensile failure, a residual elastic modulus value is given to the failed element to maintain the integrity and continuity of the entire calculation.


### 3.6. Numerical Simulation Procedure

Based on the elastoplastic damage model and the energy damage failure criterion, the code for the zonal disintegration calculation was programmed with the finite element software ABAQUS and the dynamic incremental variable plastic stiffness iterative method. Figure10 shows the flowchart of the zonal disintegration procedure.Figure 10
Flow chart of zonal disintegration subroutine.


## 3.1. Construction ofC1Order Element considering Strain Gradient

To implement the elastoplastic damage model with the finite element method, the influence of the strain gradient should be considered first. The interpolating shape functions of a conventional finite element generally have only first-order continuity, and the constructed elements have onlyC0order continuity. The second derivative after the interpolation of the displacement field is zero, and the influence of the strain gradient cannot be considered in such elements. Therefore, it is necessary to consider finite elements with higher-order continuity, such as elements with C1order continuity. This requires that the normal displacement and first-order partial differential maintain continuity so that the strain gradient can enter into the finite element equation.Taking the most versatile tetrahedral element as an example, a three-dimensional general tetrahedral element withC1order continuity is established in Figure 5. The geometric coordinates of the node i (i = 1; 2; 3; 4) are xi,yi,zi. Four-node displacements can be assumed as(12)a=a1a2a3a4T,ai=uu,xu,yu,zvv,xv,yv,zww,xw,yw,ziT.Figure 5
Tetrahedral elements withC1order continuity.The area coordinatesLi of any point Px,y,z within the element can be defined as(13)Li=AiV,i=1,2,3,4,where A1=VolP234A2=VolP134A3=VolP214A4=VolP231, and V=Vol1234.Then, the above formula can be expanded as follows:(14)Li=−1i+1ai+bix+cty+diz6v,where(15)V=161x1y1z11x2y2z21x3y3z31x4y4z4,ai=xjyjzjxkykzkxlylzl,bi=−1yjzj1ykzk1ylzl,ci=1xjzj1xkzk1xlzl,d=−1xjyj1xkyk1xlyl,where i,j,k,andl are in accordance with the rotation of 1⟶2⟶3⟶4⟶1 and meet the right-handed spiral rule.The displacement fieldu can be expressed as(16)u=∑i=14αiLi3+α5L12L2+α6L12L3+α7L12L4+α8L22L3+α9L22L4+α10L22L1+α11L32L4+α12L32L1+α13L32L2+α14L42L1+α15L42L2+α16L42L3.The values of the 16 coefficientsαii=1,…16can be determined from the displacement ui1=1,…4of the four vertices and the gradient ui,x,ui,y,ui,zi=1,…4 of the three directions.After introducing a shape functionN, the displacement field in three dimensions is expressed as(17)u=uvw=N1000N1000N1|N2000N2000N2|N3000N3000N3|N4000N4000N4a1a2a3a4=Na,where(18)Ni=Ni0NixNiyNiz=Li3+3Li2AiAjLj−AkLk+AlLl6VAiLi2pklLj+pljLk+pjkLl6VAiLi2qklLj+qljLk+qjkLl6VAiLi2rklLj+rljLk+rjkLlT,A1=b2b3b4c2c3c4d2d3d4,A2=b1b3b4c1c3c4d1d3d4,A3=b1b2b4c1c2c4d1d2d4,A4=b1b2b3c1c2c3d1d2d3,pjk=ckdj−cjdk,qjk=bjdk−bkdj,rjk=bkcj−bjck,01×4=0000.Strainε and strain gradients η take the following form:(19)ε=B1a=L1Na,η=B2a=L2Na,where L1 and L2 are differential operators that can be written as(20)L1=∂∂x∂∂y∂∂z∂∂y∂∂x∂∂z∂∂z∂∂y∂∂x,L2=M000M000M,where M=∂2/∂x2∂2/∂y2∂2/∂z2∂2/∂x∂y2∂2/∂y∂z2∂2/∂x∂z2T,01×6=000000T.Therefore, the geometric equation considering the strain gradient can be expressed as(21)εη=Ba=B1B2a,where B is the strain matrix.Finally, the element stiffness matrix of theC1order tetrahedral element is(22)K=∭BTDBdxdydz.


## 3.2. Establishment of the Dynamic Finite Element Equation

Similar to static finite element analysis, dynamic finite element analysis should establish the finite element equation. The space of the rock mass is discretized, and higher-order tetrahedral elements are constructed. The shape function and stiffness matrix of the elements are deduced. According to the dynamic equilibrium equations of each node, the Newmark-β integral method [25, 26] is adopted to solve the following motion equation:(23)Mu¨+Cu˙+Ku=P,where M is the element mass matrix; M=∫VρNTNdV,ρ is the density of the rock mass; C is the element damping matrix, C=∫VcNTNdV, c is damping coefficient; K is the element stiffness matrix; u¨,u˙, and u are the acceleration, velocity, and displacement matrices at any point in the element; and P is the external loads.Displacement at any point in the element isu=Na, where N is a function only of position and is not related to time t, and a is related to time and can be expressed as a=a1ta2ta3ta4tT. Therefore, the acceleration u¨ and velocity u˙ at any point in the element are expressed, respectively, as u¨=Na¨ and u˙=Na˙ where · indicates a derivative with respect to time.In the Newmark-β method, when the control parameters satisfy α=1/2, β=1/4, the method changes to an average constant acceleration method with second-order accuracy, which satisfies the engineering requirements. Therefore, the basic assumption based on the method can be expressed as [27](24)u˙t1+Δt=u˙t1+Δt2u¨t1+u¨t1+Δt,ut1+Δt=ut1+u˙t1Δt+Δt24u¨t1+u¨t1+Δt.Equation (24) can be expanded as(25)u¨t1+Δt=4Δt2ut1+Δt−4Δt2ut1−4Δtu˙t1−u¨t1,u¨t1+Δt=u˙t1+Δt2u¨t1+u¨t1+Δt=2Δtut1+Δt−2Δtut1−u˙t1.After subsisting equation (25) into equation (23) at time t1+Δt, the expression can be arranged as(26)K¯ut1+Δt=P¯t1+Δt,where K¯ is the equivalent stiffness matrix, K¯=4/Δt2M+2/ΔtC+K;P¯t1+Δt is the equivalent load, P¯t1+Δt=Pt1+Δt+M4/Δt2ut1+4/Δtu¨t1+u¨t1+C2/Δtut1+u˙t1.In iterative computation, the process of each time step are as follows:(1)
Calculate the equivalent loadP¯t1+Δt according to the operation state and external load of the node.(2)
Calculate the equivalent stiffness matrixK¯ according to the node.(3)
Solve equation set (24). When the structure is in the elastic state, the equivalent stiffness matrix K¯ does not change with the displacement of the node, which is a linear equation set and can be solved directly by triangular decomposition. When the structure is in the plastic state, the equivalent stiffness matrix K¯ is related to the displacement of the node and must be solved through displacement iteration.(4)
Calculate the acceleration and velocity at timet1+Δt by equation (25).


## 3.3. Dynamic Incremental Variable Plastic Stiffness Iterative Method

When the structure enters the plastic stage described by the Mohr–Coulomb plasticity criterion, equation set (26) is nonlinear. The equivalent stiffness matrix K¯ is related to displacement and must be solved iteratively. The dynamic incremental variable plasticity stiffness iteration method is used to solve the equation set [28].The equivalent stiffness matrixK¯ of dynamic calculation is made up of mass matrix M, damping matrix C, and stiffness matrix K. The plastic stiffness matrix KP is formed from K when the element enters plasticity.To guarantee the convergence of iteration, according to the incremental plastic stiffness iteration method, the plastic load of each time step is applied in multiple stages. In each stage of plastic load application, the overall stiffness matrix of the structure can be decomposed into two parts: elastic and plastic:(27)K=∭vBTDe−DpBdv=Ke−Kp.In the calculation, the elastic stiffness matrixKe remains unchanged, and the plastic stiffness matrix Kp is changed according to the critical stress state imposed after each stage load. In the process of each stage of the plastic load iteration, according to the basic idea of the incremental load method, the plastic stiffness matrix is kept constant, and the displacement δt1+Δt0 solved by the elastic equivalent stiffness matrix K¯e is taken as the initial displacement of the iterative calculation. The displacement value δt1+Δt under the incremental load at each stage is solved by the continuous correction of the displacement. The basic iterative equation for the dynamic elastoplastic finite element method is(28)K¯eut1+Δt=ΔP¯p+Kput1+Δt,K¯e=4Δt2M+2ΔtC+Ke,where K¯e is the elastic equivalent stiffness matrix.


## 3.4. Calculation Method of External Load

### 3.4.1. Stress Decomposition

Considering the complexity of deep underground engineering, a circular tunnel is used to carry out the mechanical analysis. The blasting excavation process of a deep tunnel is accompanied by the blasting load and in situ stress transient unloading, and blasting excavation is the precondition for transient unloading of in situ stress [29]. Thus, on the one hand, the stress wave produced by the blasting load and detonation gas will cause blasting damage of surrounding rock within a certain range; and on the other hand, transient unloading of stress in a deep tunnel will cause stress redistribution of surrounding rock. Therefore, the stress field of blasting excavation in Figure 6 is decomposed into two problems:(1)
The elastoplastic stress field caused by the blasting loadPbt(2)
The elastoplastic stress field caused by the transient unloading of in situ stressPotFigure 6
Forces on circular tunnel under dynamic unloading.


### 3.4.2. Blasting Load

The blasting load is actually the shock wave pressure acting on the surrounding rock after the explosion. In this paper, the change process of the blasting load is described by a triangle function (Figure7). Pb is the equivalent peak load of the blasting load on the excavation surface and tr and tdare the time of the load rising and the positive pressure, respectively.Figure 7
Explosion load curve.Taking full section blasting in a tunnel as the research object, the influence of the blasting load on rock mass is taken into account.The detonation CJ model [30, 31] is used to calculate the peak load of the blasting load. For cylindrical radial uncoupled loading conditions, the equivalent peak load Pb on the excavation surface is calculated as(29)Pb=ρ0De221+γdcdb2γdbS,where ρ0 is the density of explosives, De is the detonation velocity of explosives, γ is the isentropic exponent of explosives, which is generally taken as 3, dcis the diameter of the cartridge, dbis the diameter of the blast hole, and Sis the distance of adjacent blast holes.When the explosion wave is propagating to the blasting hole with velocityDe, the blasting load in the hole rises to the maximum, and the rising time tr can be written as(30)tr=L1D.The blasting load is a transient load that changes with time during blasting. With explosive gas escaping, the blasting load decreases further. When the pressure of the blast hole drops to atmospheric pressure, the synchronous unloading of the excavation load is completed. The duration of the blasting loadtd is calculated as(31)td=L1D+1/4S2+L22Cf+L1+L2Cu1+L1+L2Cu2,where L1 and L2 are the respective lengths of the charging section and blockage section in the blast hole; Cf is the average speed of expansion driven by crack explosion loading; Cu1is the unloading wave speed of gas explosion; and Cu2 is the reflected wave velocity spread from the bottom to the top of the blast hole [32].The time history relationPbt of the blasting load acting on the excavation surface can be expressed as follows:(32)Pbt=0,t≤0ort≥td,Pbttr0<t≤tr,Pbtd−ttd−trtr<t<td.


### 3.4.3. Transient Unloading of In Situ Stress

The unloading of in situ stress on the excavation surface [33] should satisfy the continuous stress condition, so the initial time and duration of transient unloading are determined by the blasting process. The in situ stress unloads at time tb, defined as the moment when the blasting load decreases to the in situ stress. When the blasting load is further reduced to zero, the synchronous unloading of in situ stress is completed. td−tb is the duration of transient unloading of in situ stress. The transient unloading curve of in situ stress coincides with the curve of the blasting load (Figure 8).Figure 8
Coupling effect curves of blasting load and transient release of in situ stress.The time history relation of the in situ stressP0t on the excavation surface can be written as(33)P0t=P0,t≤tb,−P0t−tdtd−tb,tb<t<td,0,t≥td.


### 3.4.4. Loading Method

In the process of blasting excavation of an underground tunnel, the external loadsP imposed on surrounding rock elements are mainly divided into blasting load Fp, gravity loads Fg, and excavation loads Fe [34].According to the stress generated by blasting on rock mass elements, the blasting loadFp can be obtained as(34)Fp=∭vBTPbtdv,where Fg is the vertical gravity loads caused by the removed elements and can be expressed as(35)Fg=∭vNTρgNdv,where ρ is the density of the rock mass, g is the acceleration of gravity, and N is the shape function matrix which can be obtained by Hermit interpolation function. Fe, the transient unloading of in situ stress, can be expressed as(36)Fe=∭vBTP0tdv,where B is strain matrix and P0t is the in situ stress of the excavation unit, obtained from equation (33).


## 3.4.1. Stress Decomposition

Considering the complexity of deep underground engineering, a circular tunnel is used to carry out the mechanical analysis. The blasting excavation process of a deep tunnel is accompanied by the blasting load and in situ stress transient unloading, and blasting excavation is the precondition for transient unloading of in situ stress [29]. Thus, on the one hand, the stress wave produced by the blasting load and detonation gas will cause blasting damage of surrounding rock within a certain range; and on the other hand, transient unloading of stress in a deep tunnel will cause stress redistribution of surrounding rock. Therefore, the stress field of blasting excavation in Figure 6 is decomposed into two problems:(1)
The elastoplastic stress field caused by the blasting loadPbt(2)
The elastoplastic stress field caused by the transient unloading of in situ stressPotFigure 6
Forces on circular tunnel under dynamic unloading.


## 3.4.2. Blasting Load

The blasting load is actually the shock wave pressure acting on the surrounding rock after the explosion. In this paper, the change process of the blasting load is described by a triangle function (Figure7). Pb is the equivalent peak load of the blasting load on the excavation surface and tr and tdare the time of the load rising and the positive pressure, respectively.Figure 7
Explosion load curve.Taking full section blasting in a tunnel as the research object, the influence of the blasting load on rock mass is taken into account.The detonation CJ model [30, 31] is used to calculate the peak load of the blasting load. For cylindrical radial uncoupled loading conditions, the equivalent peak load Pb on the excavation surface is calculated as(29)Pb=ρ0De221+γdcdb2γdbS,where ρ0 is the density of explosives, De is the detonation velocity of explosives, γ is the isentropic exponent of explosives, which is generally taken as 3, dcis the diameter of the cartridge, dbis the diameter of the blast hole, and Sis the distance of adjacent blast holes.When the explosion wave is propagating to the blasting hole with velocityDe, the blasting load in the hole rises to the maximum, and the rising time tr can be written as(30)tr=L1D.The blasting load is a transient load that changes with time during blasting. With explosive gas escaping, the blasting load decreases further. When the pressure of the blast hole drops to atmospheric pressure, the synchronous unloading of the excavation load is completed. The duration of the blasting loadtd is calculated as(31)td=L1D+1/4S2+L22Cf+L1+L2Cu1+L1+L2Cu2,where L1 and L2 are the respective lengths of the charging section and blockage section in the blast hole; Cf is the average speed of expansion driven by crack explosion loading; Cu1is the unloading wave speed of gas explosion; and Cu2 is the reflected wave velocity spread from the bottom to the top of the blast hole [32].The time history relationPbt of the blasting load acting on the excavation surface can be expressed as follows:(32)Pbt=0,t≤0ort≥td,Pbttr0<t≤tr,Pbtd−ttd−trtr<t<td.


## 3.4.3. Transient Unloading of In Situ Stress

The unloading of in situ stress on the excavation surface [33] should satisfy the continuous stress condition, so the initial time and duration of transient unloading are determined by the blasting process. The in situ stress unloads at time tb, defined as the moment when the blasting load decreases to the in situ stress. When the blasting load is further reduced to zero, the synchronous unloading of in situ stress is completed. td−tb is the duration of transient unloading of in situ stress. The transient unloading curve of in situ stress coincides with the curve of the blasting load (Figure 8).Figure 8
Coupling effect curves of blasting load and transient release of in situ stress.The time history relation of the in situ stressP0t on the excavation surface can be written as(33)P0t=P0,t≤tb,−P0t−tdtd−tb,tb<t<td,0,t≥td.


## 3.4.4. Loading Method

In the process of blasting excavation of an underground tunnel, the external loadsP imposed on surrounding rock elements are mainly divided into blasting load Fp, gravity loads Fg, and excavation loads Fe [34].According to the stress generated by blasting on rock mass elements, the blasting loadFp can be obtained as(34)Fp=∭vBTPbtdv,where Fg is the vertical gravity loads caused by the removed elements and can be expressed as(35)Fg=∭vNTρgNdv,where ρ is the density of the rock mass, g is the acceleration of gravity, and N is the shape function matrix which can be obtained by Hermit interpolation function. Fe, the transient unloading of in situ stress, can be expressed as(36)Fe=∭vBTP0tdv,where B is strain matrix and P0t is the in situ stress of the excavation unit, obtained from equation (33).


## 3.5. Zonal Disintegration Damage Failure Criterion

### 3.5.1. Maximum Tensile Strain Criterion

Tensile failure of surrounding rock is generally the result of strain development in the direction perpendicular to the surface of the wall. The maximum tensile strain criterion can be written as(37)f=εtmax−εtu=0,where εtmax is the maximum tensile strain of a calculation element [35] and εtu is the ultimate tensile strain, which is obtained by a uniaxial tensile test or Brazilian disk test.Iff≥0, then tensile failure will occur and the element is destroyed. To maintain the integrity and continuity of the entire calculation during the numerical simulation, a small residual elastic modulus EC is given to the element where the tensile failure occurs; usually, EC=0.05E. If f < 0, the tensile failure will not occur. The zonal disintegration energy damage failure criterion is used as the criterion of element failure.


### 3.5.2. Zonal Disintegration Energy Damage Failure Criterion

Considering the influence of the strain gradient term, the formula for calculating the strain energy density of an element is written as(38)dWdV=∫0εijσijdεij+18Gl24ηijjηikk+4ηiikηkjj+4ηiikηjjk+7ηijkηijk+2ηijkηkji.Figure9 shows the simplified stress-strain relation curve of rock under high stress. Elastic strain energy density dW/dVs, critical strain energy density dW/dVo, and ultimate strain energy density dW/dVu can be calculated as [36](39)dWdVe=∫0εeσijdεij,dWdVo=∫0εoσijdεij,dWdVu=∫0εuσijdεij.Figure 9
Stress-strain relation curve of rock under high stress.The damage failure criterion of zonal disintegration is judged by the change of the strain energy density.(1)
WhendW/dV<dW/dVs=∫0εsσijdεij (εsis the yield strain in the loading process), then the element is in the linear elastic stage, and damage within the element will not occur.(2)
WhendW/dV≥dW/dVs, the element enters into the plastic damage stage; microcracks in the rock mass appear and begin to expand.(3)
WhendW/dV≥dW/dVD=∫0εoσijdεij (εois the critical strain in the loading process), the microcracks in the rock mass develop into a macroscopic crack, and there is still bearing capacity in the rock body. This is called the plastic residual stage.(4)
WhendW/dV≥dW/dVu=∫0εuσijdεij (εu is the ultimate strain in the loading process), the element is damaged to failure. In the case of tensile failure, a residual elastic modulus value is given to the failed element to maintain the integrity and continuity of the entire calculation.


## 3.5.1. Maximum Tensile Strain Criterion

Tensile failure of surrounding rock is generally the result of strain development in the direction perpendicular to the surface of the wall. The maximum tensile strain criterion can be written as(37)f=εtmax−εtu=0,where εtmax is the maximum tensile strain of a calculation element [35] and εtu is the ultimate tensile strain, which is obtained by a uniaxial tensile test or Brazilian disk test.Iff≥0, then tensile failure will occur and the element is destroyed. To maintain the integrity and continuity of the entire calculation during the numerical simulation, a small residual elastic modulus EC is given to the element where the tensile failure occurs; usually, EC=0.05E. If f < 0, the tensile failure will not occur. The zonal disintegration energy damage failure criterion is used as the criterion of element failure.


## 3.5.2. Zonal Disintegration Energy Damage Failure Criterion

Considering the influence of the strain gradient term, the formula for calculating the strain energy density of an element is written as(38)dWdV=∫0εijσijdεij+18Gl24ηijjηikk+4ηiikηkjj+4ηiikηjjk+7ηijkηijk+2ηijkηkji.Figure9 shows the simplified stress-strain relation curve of rock under high stress. Elastic strain energy density dW/dVs, critical strain energy density dW/dVo, and ultimate strain energy density dW/dVu can be calculated as [36](39)dWdVe=∫0εeσijdεij,dWdVo=∫0εoσijdεij,dWdVu=∫0εuσijdεij.Figure 9
Stress-strain relation curve of rock under high stress.The damage failure criterion of zonal disintegration is judged by the change of the strain energy density.(1)
WhendW/dV<dW/dVs=∫0εsσijdεij (εsis the yield strain in the loading process), then the element is in the linear elastic stage, and damage within the element will not occur.(2)
WhendW/dV≥dW/dVs, the element enters into the plastic damage stage; microcracks in the rock mass appear and begin to expand.(3)
WhendW/dV≥dW/dVD=∫0εoσijdεij (εois the critical strain in the loading process), the microcracks in the rock mass develop into a macroscopic crack, and there is still bearing capacity in the rock body. This is called the plastic residual stage.(4)
WhendW/dV≥dW/dVu=∫0εuσijdεij (εu is the ultimate strain in the loading process), the element is damaged to failure. In the case of tensile failure, a residual elastic modulus value is given to the failed element to maintain the integrity and continuity of the entire calculation.


## 3.6. Numerical Simulation Procedure

Based on the elastoplastic damage model and the energy damage failure criterion, the code for the zonal disintegration calculation was programmed with the finite element software ABAQUS and the dynamic incremental variable plastic stiffness iterative method. Figure10 shows the flowchart of the zonal disintegration procedure.Figure 10
Flow chart of zonal disintegration subroutine.


## 4. Engineering Application

### 4.1. In Situ Monitoring Results

The deep tunnel of the Dingji Mine in the Huainan area was selected as the research object. The depth of the tunnel is 910 m, and the monitoring section size is 5,000 × 3,880 mm. The damaged scopes of the surrounding rock in the boreholes are depicted on the map, and the fractured areas within 0.5 m in different holes are connected to form fractured zones with different depths. The zones between the fractured zones are relatively complete (Figure11) [37, 38].Figure 11
Failure distribution of model cavern after excavation.The average ranges of the fractured zones are listed in Table1. The inner and outer diameters are the sides near and away from the tunnel in the distance, respectively.Table 1
In situ test results of fractured zones.
Fractured zones numberField measured resultsInner diameter (m)Outer diameter (m)12.504.9925.756.5537.748.7849.8710.40


### 4.2. Determination of Calculation Parameters

The calculation parameters of the tunnel’s surrounding rock are as follows: the tunnel shape is a semicircular arch; the section size is 5 m × 3.88 m; the simulation range of the calculation area is 30 m × 30 m × 30 m; the initial geostress isP0 = 23.5 MPa; the rock density ρ  = 2.45 g/cm3; the compressive strength σc  = 88.55 MPa; the elastic modulus E = 77.82 GPa; Poisson's ratio υ  = 0.286; the cohesion force c = 9 MPa; the internal friction angle φ  = 43°; the internal material length l = 0.01 m; the yield strain εs  = 0.725 × 10−3; and the ultimate strain εu = 4.65 × 10−3.Figure12 shows the layout of the hole in the tunnel excavation blasting. The blasting design of the underground tunnel in the actual project is simplified, and the full section blasting is used. Only the effect of the blasting load on surrounding rock at the excavation radius (periphery holes and bottom holes) is considered. For tunneling blasting, the diameter of the blast hole is db  = 42 mm and the distance of adjacent blast holes is S = 0.3 m; the explosive is grade-three water-gel explosive for permissible coal mine; the diameter of the cartridge is dc  = 35 mm, the length of the charging section is L1  = 1.4 m; the length of the blockage section is L2  = 0.6 m; the density of the explosive is ρ0  = 0.95–1.25 g/cm3; the detonation velocity is De  = 3.0 × 103 m/s; the isentropic exponent is γ  = 3; and the equivalent blasting peak load on the surface of the blasting excavation can be calculated by equation (33), in which the value is Pb  = 52.75 MPa.Figure 12
The schematic diagram of blast hole arrangement.It is concluded that the unloading wave speed of gas explosionCu1 is approximately equal to the reflected wave velocity Cu2, that is, Cu1  =  Cu2 = (1–1.5) × 103 m/s. The average speed of expansion driven by crack explosion loading is Cf  = 0.25 CL, where CLis the longitudinal velocity of rock mass, CL  = 4550 × 103 m/s. According to equations (30) and (31), it can be obtained that the rise time of the blasting load is tr  = 470 μs, the duration of the blasting load is td  = 5700μs, and the unloading time of in situ stress is tb  = 2130 μs.The calculation step is determined by the demanded precision, so the initial space step is 0.01 m, and the initial time step is 10μs.


### 4.3. Establishment of the ABAQUS Numerical Model

Figure13 shows the three-dimensional computational mesh for the selected calculation range. There are 110,360 elements and 157,837 nodes.Figure 13
Three-dimensional computational mesh of rock mass.In the dynamic calculation, the outwardly propagating stress wave will be reflected back to the model interior at the constraint boundary of the model, causing the energy not to be dissipated outwards. As the time increases, the model elements will be completely destroyed and the result will be completely distorted. To solve this problem, a free-field boundary [39–41] is used as the boundary conditions of the computational model and the outward waves generated by the structure can be properly absorbed.


### 4.4. Simulation Results and Comparative Analysis

Figure14 shows the numerical nephogram obtained by numerical simulation.Figure 14
Fractured zone nephogram of surrounding rock. (a)t=250μs. (b)t=500μs. (c)t=2500μs. (d)t=6000μs. (e) Ultimate time.
(a)(b)(c)(d)(e)At the time of 250μs, the radial displacement of the tunnel reaches 4.91 cm, and a damage zone is generated within the range of about 2.05 m around the tunnel as in Figure 14(a).At the time of 500μs, the damage zone around the tunnel increases significantly, and the radial displacement increases to 6.85 cm as in Figure 14(b).At the time of 2500μs, there is a sudden increase in radial displacement at 5.5 m, where the second fractured zones occur as in Figure 14(c).At the time of 6000μs, the ultimate time, a third fractured zone occurs about 7.8 m away from the tunnel periphery as in Figures 14(d) and 14(e).Figure15 shows the stresses and displacement of the surrounding rock in a deep tunnel caused by blasting excavation. A wave phenomenon can be seen in the radial displacement, tangential stress, and radial stress of the surrounding rock. From the curves in Figure 15, the following trends can be observed.(1)
At timet=250μs, the blasting load is at the rising stage; the radial stress and tangential stress in the surrounding rock near the blasting area are relatively high. With the increase of the radius, the stress decreases rapidly, and the attenuation speed decreases with the increase of the radius. The stress in surrounding rock at moderate and blasting areas tends to be gentle. When the tangential tensile stress in the surrounding rock is greater than the tensile strength, the radial tensile cracks will occur and lead to the fracture phenomenon. The radial displacement of the surrounding rock will also increase.(2)
At the timet=500μs, the blasting load is around the peak load, and the stress and displacement values of the surrounding rock increase to a certain extent.(3)
At the timet=2500μs, the stress oscillation of the surrounding rock appears at a moderate and remote distance under the interaction of blasting load and in situ stress unloading. As shown in Figures 8(a) and 8(b), the radial and tangential stress are both in the trough state at a radius of 5.65 m, and the variation of radial displacement shows a rapidly increasing trend from a radius of 5 m and reaches the peak at the radius of 5.8 m. It is known from the previous analysis that the surrounding rock in this region is damaged and the fractured zone is produced.(4)
At the timet=6000μs, the unloading of the dynamic load in the tunnel is completed. The radial stress changes from tension to compression at the near region. The tangential arc crack will be generated when the stress exceeds the tensile strength of the rock mass. In the moderate and remote area of the tunnel, the stresses in the surrounding rock present a more obvious changing mode. The radial displacement increases continuously in the original fractured zone and rises rapidly in the new fractured zone.Figure 15
The displacement and stress changes of tunnel surrounding rock in different times. (a) Radial displacement. (b) Radial stress. (c) Tangential stress.
(a)(b)(c)It can be realized from the values at the ultimate time that the displacement and stresses present an oscillating mode, which is reflected in the alternate appearance of the wave peak and trough. By comparative analysis of Figures14 and 15, the surrounding rock stresses at the relatively complete zones are in the peak ranges while the displacements are in the trough ranges; on the contrary, the surrounding rock stresses in the fractured zones are in the trough ranges while the displacements are in the peak ranges.The numerical calculation of the fractured zones is compared with the in situ results in Figure16. It is worth mentioning that the average radius is the mean value of the inner and outer diameter of each fractured zone, and the average width is the mean value of the difference between the internal and external diameters of each fractured zone.Figure 16
Width and quantity comparison in surrounding rock of tunnel.We can conclude that the range, width, and position of the fractured zones obtained by numerical simulation are basically consistent with the field measurement results, indicating that the numerical calculation method is feasible and reliable.


## 4.1. In Situ Monitoring Results

The deep tunnel of the Dingji Mine in the Huainan area was selected as the research object. The depth of the tunnel is 910 m, and the monitoring section size is 5,000 × 3,880 mm. The damaged scopes of the surrounding rock in the boreholes are depicted on the map, and the fractured areas within 0.5 m in different holes are connected to form fractured zones with different depths. The zones between the fractured zones are relatively complete (Figure11) [37, 38].Figure 11
Failure distribution of model cavern after excavation.The average ranges of the fractured zones are listed in Table1. The inner and outer diameters are the sides near and away from the tunnel in the distance, respectively.Table 1
In situ test results of fractured zones.
Fractured zones numberField measured resultsInner diameter (m)Outer diameter (m)12.504.9925.756.5537.748.7849.8710.40


## 4.2. Determination of Calculation Parameters

The calculation parameters of the tunnel’s surrounding rock are as follows: the tunnel shape is a semicircular arch; the section size is 5 m × 3.88 m; the simulation range of the calculation area is 30 m × 30 m × 30 m; the initial geostress isP0 = 23.5 MPa; the rock density ρ  = 2.45 g/cm3; the compressive strength σc  = 88.55 MPa; the elastic modulus E = 77.82 GPa; Poisson's ratio υ  = 0.286; the cohesion force c = 9 MPa; the internal friction angle φ  = 43°; the internal material length l = 0.01 m; the yield strain εs  = 0.725 × 10−3; and the ultimate strain εu = 4.65 × 10−3.Figure12 shows the layout of the hole in the tunnel excavation blasting. The blasting design of the underground tunnel in the actual project is simplified, and the full section blasting is used. Only the effect of the blasting load on surrounding rock at the excavation radius (periphery holes and bottom holes) is considered. For tunneling blasting, the diameter of the blast hole is db  = 42 mm and the distance of adjacent blast holes is S = 0.3 m; the explosive is grade-three water-gel explosive for permissible coal mine; the diameter of the cartridge is dc  = 35 mm, the length of the charging section is L1  = 1.4 m; the length of the blockage section is L2  = 0.6 m; the density of the explosive is ρ0  = 0.95–1.25 g/cm3; the detonation velocity is De  = 3.0 × 103 m/s; the isentropic exponent is γ  = 3; and the equivalent blasting peak load on the surface of the blasting excavation can be calculated by equation (33), in which the value is Pb  = 52.75 MPa.Figure 12
The schematic diagram of blast hole arrangement.It is concluded that the unloading wave speed of gas explosionCu1 is approximately equal to the reflected wave velocity Cu2, that is, Cu1  =  Cu2 = (1–1.5) × 103 m/s. The average speed of expansion driven by crack explosion loading is Cf  = 0.25 CL, where CLis the longitudinal velocity of rock mass, CL  = 4550 × 103 m/s. According to equations (30) and (31), it can be obtained that the rise time of the blasting load is tr  = 470 μs, the duration of the blasting load is td  = 5700μs, and the unloading time of in situ stress is tb  = 2130 μs.The calculation step is determined by the demanded precision, so the initial space step is 0.01 m, and the initial time step is 10μs.


## 4.3. Establishment of the ABAQUS Numerical Model

Figure13 shows the three-dimensional computational mesh for the selected calculation range. There are 110,360 elements and 157,837 nodes.Figure 13
Three-dimensional computational mesh of rock mass.In the dynamic calculation, the outwardly propagating stress wave will be reflected back to the model interior at the constraint boundary of the model, causing the energy not to be dissipated outwards. As the time increases, the model elements will be completely destroyed and the result will be completely distorted. To solve this problem, a free-field boundary [39–41] is used as the boundary conditions of the computational model and the outward waves generated by the structure can be properly absorbed.


## 4.4. Simulation Results and Comparative Analysis

Figure14 shows the numerical nephogram obtained by numerical simulation.Figure 14
Fractured zone nephogram of surrounding rock. (a)t=250μs. (b)t=500μs. (c)t=2500μs. (d)t=6000μs. (e) Ultimate time.
(a)(b)(c)(d)(e)At the time of 250μs, the radial displacement of the tunnel reaches 4.91 cm, and a damage zone is generated within the range of about 2.05 m around the tunnel as in Figure 14(a).At the time of 500μs, the damage zone around the tunnel increases significantly, and the radial displacement increases to 6.85 cm as in Figure 14(b).At the time of 2500μs, there is a sudden increase in radial displacement at 5.5 m, where the second fractured zones occur as in Figure 14(c).At the time of 6000μs, the ultimate time, a third fractured zone occurs about 7.8 m away from the tunnel periphery as in Figures 14(d) and 14(e).Figure15 shows the stresses and displacement of the surrounding rock in a deep tunnel caused by blasting excavation. A wave phenomenon can be seen in the radial displacement, tangential stress, and radial stress of the surrounding rock. From the curves in Figure 15, the following trends can be observed.(1)
At timet=250μs, the blasting load is at the rising stage; the radial stress and tangential stress in the surrounding rock near the blasting area are relatively high. With the increase of the radius, the stress decreases rapidly, and the attenuation speed decreases with the increase of the radius. The stress in surrounding rock at moderate and blasting areas tends to be gentle. When the tangential tensile stress in the surrounding rock is greater than the tensile strength, the radial tensile cracks will occur and lead to the fracture phenomenon. The radial displacement of the surrounding rock will also increase.(2)
At the timet=500μs, the blasting load is around the peak load, and the stress and displacement values of the surrounding rock increase to a certain extent.(3)
At the timet=2500μs, the stress oscillation of the surrounding rock appears at a moderate and remote distance under the interaction of blasting load and in situ stress unloading. As shown in Figures 8(a) and 8(b), the radial and tangential stress are both in the trough state at a radius of 5.65 m, and the variation of radial displacement shows a rapidly increasing trend from a radius of 5 m and reaches the peak at the radius of 5.8 m. It is known from the previous analysis that the surrounding rock in this region is damaged and the fractured zone is produced.(4)
At the timet=6000μs, the unloading of the dynamic load in the tunnel is completed. The radial stress changes from tension to compression at the near region. The tangential arc crack will be generated when the stress exceeds the tensile strength of the rock mass. In the moderate and remote area of the tunnel, the stresses in the surrounding rock present a more obvious changing mode. The radial displacement increases continuously in the original fractured zone and rises rapidly in the new fractured zone.Figure 15
The displacement and stress changes of tunnel surrounding rock in different times. (a) Radial displacement. (b) Radial stress. (c) Tangential stress.
(a)(b)(c)It can be realized from the values at the ultimate time that the displacement and stresses present an oscillating mode, which is reflected in the alternate appearance of the wave peak and trough. By comparative analysis of Figures14 and 15, the surrounding rock stresses at the relatively complete zones are in the peak ranges while the displacements are in the trough ranges; on the contrary, the surrounding rock stresses in the fractured zones are in the trough ranges while the displacements are in the peak ranges.The numerical calculation of the fractured zones is compared with the in situ results in Figure16. It is worth mentioning that the average radius is the mean value of the inner and outer diameter of each fractured zone, and the average width is the mean value of the difference between the internal and external diameters of each fractured zone.Figure 16
Width and quantity comparison in surrounding rock of tunnel.We can conclude that the range, width, and position of the fractured zones obtained by numerical simulation are basically consistent with the field measurement results, indicating that the numerical calculation method is feasible and reliable.


## 5. Conclusion

The effect of strain gradient should be considered in the study of the zonal disintegration of deep tunnels.Based on the strain gradient theory and the damage characteristics of rock mass, an elastoplastic damage model was developed. The zonal disintegration calculation program was compiled with the dynamic finite element equation and the dynamic incremental variable stiffness iterative method, which is embedded in the finite element software ABAQUS. According to the failure criterion, the number and depth of fractured zones are obtained under the action of blasting excavation. The failure process of the surrounding rock in the deep tunnel can be observed dynamically, and the formation mechanism of the zonal disintegration phenomenon in the deep tunnel can be clearly identified.In the numerical calculation, the radial displacement, radial stress, and tangential stress of the surrounding rock appear in the changing law of the wave peak and trough alternation. The width and quantity of the fractured zone are in good agreement with the in situ observations. The applicability of the numerical simulation method is analyzed.When the tunnel is disturbed by blasting excavation, the coupling effect of the blasting load and the initial geostress transient unloading leads to the variation of alternation oscillation in the surrounding rock stress field, which is an important reason for the zonal disintegration of the surrounding rock.



---

*Source: 9991251-2021-05-12.xml*