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

More deep tunneling projects will be constructed due to the increasing demand of underground energy and resource. +e 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. +e dynamic formation and development regularity of zonal disintegration in the deep tunnel are simulated by this subroutine.+e radial displacement, radial stress, and tangential stress show the oscillated variation of peaks and troughs alternately.+e 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. +e 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. +e method proposed in the current study can be utilized to provide a basis for the prediction and supporting design of fractured modes.


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. e 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]. e zonal disintegration phenomenon has been confirmed through many physical detection methods in the excavation of deep tunnels. e 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]. is phenomenon was observed by borehole television and an electrical resistivity survey in the Huainan coal mine in China ( Figure 2) [4]. e zonal disintegration phenomenon was reproduced in a model test. en 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].
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][12][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]. e numerical results indicated that the dynamic loads and high in situ stress are two dominant factors in the occurrence of zonal disintegration. e 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.
ere 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. e zonal disintegration calculation program was written based on the dynamic finite element equations and ABAQUS finite element software. e morphological characteristics of zonal disintegration under blasting excavation were obtained through numerical simulation.
is study reveals the morphological characteristics of zonal disintegration, which is used as a basis for the prediction and support control of fractured modes.

Introduction of High-Order Stress.
In the Toupin-Mindlin strain gradient theory [18][19][20], the total strain u i has two parts: the conventional Eulerian strain tensor ε ij and the high-order strain tensor η ijk , which are deduced as According to the second law of thermodynamics, the Helmholtz free energy function [18] can be expressed as 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 . e stress tensor σ ij and highorder stress tensor τ ijk can be derived from the free energy function corresponding to the strain tensor ε ij and strain gradient tensor η ijk , respectively: ijkl ; d is the damage tensor; E e ijkl is the elastic tensor; E p ijkl is the plastic tensor determined by plastic flow law and hardening law [21]; Λ ed ijklmn is the sixth-order elastic damage tensor considering the strain gradient, ; δ kn is 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].
H is 6 × 6 matrixH � H e − H p , and H e can be written as where E is the elastic modulus and υ is Poisson's ratio. H p is expressed as where F � F(σ, κ) � 0 is the yield function,κ is plastic work, dκ � σdε p , andQ is plastic potential function.
us, the elastoplastic damage constitutive equation based on the strain gradient can be expressed as

Damage Evolution Equation.
It is assumed that the damage of a rock mass can be represented by an isotropic damage variable d 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. e 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 Figure 4(a). e 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. e AB stage is the plastic damage stage. e 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] where ε is the generalized equivalent strain containing the strain gradient [24], ε � variation curve of damage variable d with equivalent strain ε is shown in Figure 4(b). us, equations (13) and (14) constitute an elastoplastic damage model based on the strain gradient.

Construction of C 1 Order 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. e interpolating shape functions of a conventional finite element generally have only first-order continuity, and the constructed elements have only C 0 order continuity. e second derivative after the interpolation of the displacement field is zero, and the influence of the strain gradient cannot be considered in such elements. erefore, it is necessary to consider finite elements with higher-order continuity, such as elements with C 1 order continuity. is 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 withC 1 order continuity is established in Figure 5. e geometric coordinates of the node i (i � 1; 2; 3; 4) are (x i , y i , z i ). Four-node displacements can be assumed as e area coordinates L i of any point P(x, y, z) within the element can be defined as where , and V � Vol(1234). en, the above formula can be expanded as follows: where where i, j, k, and l are in accordance with the rotation of 1 ⟶ 2 ⟶ 3 ⟶ 4 ⟶ 1 and meet the right-handed spiral rule. e displacement field u can be expressed as

Shock and Vibration 5
After introducing a shape function N, the displacement field in three dimensions is expressed as where Strain ε and strain gradients η take the following form: where L 1 and L 2 are differential operators that can be written as Figure 5: Tetrahedral elements withC 1 order continuity. where T . erefore, the geometric equation considering the strain gradient can be expressed as where B is the strain matrix. Finally, the element stiffness matrix of the C 1 order tetrahedral element is

Establishment of the Dynamic Finite Element Equation.
Similar to static finite element analysis, dynamic finite element analysis should establish the finite element equation. e space of the rock mass is discretized, and higher-order tetrahedral elements are constructed. e 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: where [M] is the element mass matrix; { } 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 is u � 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 T . erefore, the acceleration € u and velocity _ u at any point in the element are expressed, respectively, as € u � N € a and _ u � N _ a 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. erefore, the basic assumption based on the method can be expressed as [27] Equation (24) can be expanded as After subsisting equation (25) into equation (23) at time t1 + Δt, the expression can be arranged as 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: In the calculation, the elastic stiffness matrix [K e ] remains unchanged, and the plastic stiffness matrix [K p ] 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+Δt 0 solved by the elastic equivalent stiffness matrix [K e ] is taken as the initial displacement of the iterative calculation. e displacement value δ t1+Δt under the incremental load at each stage is solved by the continuous correction of the displacement. e basic iterative equation for the dynamic elastoplastic finite element method is K e u t1+Δt � ΔP p + K p u t1+Δt , where [K e ] is the elastic equivalent stiffness matrix.

Stress Decomposition.
Considering the complexity of deep underground engineering, a circular tunnel is used to carry out the mechanical analysis. e 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]. us, 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. erefore, the stress field of blasting excavation in Figure 6 is decomposed into two problems: (1) e elastoplastic stress field caused by the blasting load P b (t) (2) e elastoplastic stress field caused by the transient unloading of in situ stress P o (t)

Blasting Load.
e 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 (Figure 7). P b is the equivalent peak load of the blasting load on the excavation surface and t r and t d are the time of the load rising and the positive pressure, respectively.
Taking full section blasting in a tunnel as the research object, the influence of the blasting load on rock mass is taken into account. e 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 P b on the excavation surface is calculated as where ρ 0 is the density of explosives, D e is the detonation velocity of explosives, c is the isentropic exponent of explosives, which is generally taken as 3, d c is the diameter of the cartridge, d b is 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 velocity D e , the blasting load in the hole rises to the maximum, and the rising time t r can be written as e 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. e duration of the blasting load t d is calculated as where L 1 and L 2 are the respective lengths of the charging section and blockage section in the blast hole; C f is the average speed of expansion driven by crack explosion loading; C u1 is the unloading wave speed of gas explosion; and C u2 is the reflected wave velocity spread from the bottom to the top of the blast hole [32]. e time history relation P b (t) of the blasting load acting on the excavation surface can be expressed as follows:

Transient Unloading of In Situ
Stress. e 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. e in situ stress unloads at time t b , 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. t d − t b is the duration of transient unloading of in situ stress. e transient unloading curve of in situ stress coincides with the curve of the blasting load (Figure 8).
e time history relation of the in situ stress P 0 (t) on the excavation surface can be written as

Loading Method.
In the process of blasting excavation of an underground tunnel, the external loads P { } imposed on surrounding rock elements are mainly divided into blasting load F p , gravity loads F g , and excavation loads F e [34].
According to the stress generated by blasting on rock mass elements, the blasting load F p can be obtained as where F g is the vertical gravity loads caused by the removed elements and can be expressed as 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. F e , the transient unloading of in situ stress, can be expressed as where [B] is strain matrix and P 0 (t) is the in situ stress of the excavation unit, obtained from equation (33).

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. e maximum tensile strain criterion can be written as 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. If f ≥ 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 E C is given to the element where the tensile failure occurs; usually, E C � 0.05E. If f < 0, the tensile failure will not occur. e zonal disintegration energy damage failure criterion is used as the criterion of element failure.

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) Figure 9 shows the simplified stress-strain relation curve of rock under high stress. Elastic strain energy density (dW/dV) s , critical strain energy density (dW/dV) o , and P 0 ultimate strain energy density (dW/dV) u can be calculated as [36] dW dV e � ε e 0 σ ij dε ij , e damage failure criterion of zonal disintegration is judged by the change of the strain energy density.
(1) When (dW/dV) < (dW/dV) s � ε s 0 σ ij dε ij (ε s is 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) When (dW/dV) ≥ (dW/dV) s , the element enters into the plastic damage stage; microcracks in the rock mass appear and begin to expand.
(3) When (dW/dV) ≥ (dW/dV) D � ε o 0 σ ij dε ij (ε o is 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. is is called the plastic residual stage. (4) When (dW/dV) ≥ (dW/dV) u � ε u 0 σ ij dε 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.

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. Figure 10 shows the flowchart of the zonal disintegration procedure.

In Situ Monitoring Results.
e deep tunnel of the Dingji Mine in the Huainan area was selected as the research object. e depth of the tunnel is 910 m, and the monitoring section size is 5,000 × 3,880 mm. e 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. e zones between the fractured zones are relatively complete ( Figure 11) [37,38]. e average ranges of the fractured zones are listed in Table 1. e inner and outer diameters are the sides near and away from the tunnel in the distance, respectively.

Determination of Calculation Parameters.
e 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 is P 0 � 23.5 MPa; the rock density ρ � 2.45 g/cm 3 ; 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 . Figure 12 shows the layout of the hole in the tunnel excavation blasting. e 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 d b � 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 d c � 35 mm, the length of the charging section is L 1 � 1.4 m; the length of the blockage section is L 2    Relatively complete zones Fractured zones Figure 11: Failure distribution of model cavern after excavation.

Shock and Vibration 11
� 0.6 m; the density of the explosive is ρ 0 � 0.95-1.25 g/cm 3 ; the detonation velocity is D e � 3.0 × 10 3 m/s; the isentropic exponent is c � 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 P b � 52.75 MPa. It is concluded that the unloading wave speed of gas explosion C u1 is approximately equal to the reflected wave velocity C u2 , that is, C u1 � C u2 � (1-1.5) × 10 3 m/s. e average speed of expansion driven by crack explosion loading is C f � 0.25 C L , where C L is the longitudinal velocity of rock mass, C L � 4550 × 10 3 m/s. According to equations (30) and (31), it can be obtained that the rise time of the blasting load is t r � 470 μs, the duration of the blasting load is t d � 5700μs, and the unloading time of in situ stress is t b � 2130 μs. e 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. Figure 13 shows the three-dimensional computational mesh for the selected calculation range.

Establishment of the ABAQUS Numerical Model.
ere are 110,360 elements and 157,837 nodes.
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][40][41] is used as the boundary conditions of the computational model and the outward waves generated by the structure can be properly absorbed. Figure 14 shows the numerical nephogram obtained by numerical simulation.

Simulation Results and Comparative Analysis.
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). Figure 15 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 time t � 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. e 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. e radial displacement of the surrounding rock will also increase.
(2) At the time t � 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 time t � 2500 μs, the stress oscillation of the surrounding rock appears at a moderate and remote      , 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 time t � 6000 μs, the unloading of the dynamic load in the tunnel is completed. e radial stress changes from tension to compression at the near region. e 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. e radial displacement increases continuously in the original fractured zone and rises rapidly in the new fractured zone.
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 Figures 14 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. e numerical calculation of the fractured zones is compared with the in situ results in Figure 16. 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.
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.

Conclusion
e 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. e 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. e 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. e width and quantity of the fractured zone are in good agreement with the in situ observations. e 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.

Data Availability
e data used to support the findings of this study are available from the corresponding authors upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.