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.

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

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

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 [1113]. 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 [1820], the total strain has two parts: the conventional Eulerian strain tensor and the high-order strain tensor , which are deduced as

According to the second law of thermodynamics, the Helmholtz free energy function [18] can be expressed aswhere is a second-order Cauchy stress tensor, which is the conjugate of the Eulerian strain tensor , and is set as the third-order stress tensor which is conjugate to the strain gradient tensor . The stress tensor and high-order stress tensor can be derived from the free energy function corresponding to the strain tensor and strain gradient tensor , respectively:where is the elastoplastic damage tensor, ; d is the damage tensor; is the elastic tensor; is the plastic tensor determined by plastic flow law and hardening law [21]; is the sixth-order elastic damage tensor considering the strain gradient, ; 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].

2.2. Establishment of Elastoplastic Damage Constitutive Equation

Cauchy stress and Eulerian strain are expressed as an incremental matrix:where and .

is matrix, and can be written aswhere is the elastic modulus and is Poisson’s ratio. is expressed as where is the yield function, is plastic work, , and is plastic potential function.

The high-order stress tensor and strain gradient tensor are expressed as an incremental matrix:where

is matrix,, and can be written as

Thus, the elastoplastic damage constitutive equation based on the strain gradient can be expressed as

2.3. 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. 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 Figure 4(a). The OA stage is the linear elastic stage, in which no damage occurs to the rock mass. is the yield stress and 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. 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], . The variation curve of damage variable d with equivalent strain is shown in Figure 4(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 of 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. The interpolating shape functions of a conventional finite element generally have only first-order continuity, and the constructed elements have only order 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 order 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 withorder continuity is established in Figure 5. The geometric coordinates of the node i (i = 1; 2; 3; 4) are . Four-node displacements can be assumed as

The area coordinates of any point within the element can be defined aswhere , and .

Then, the above formula can be expanded as follows:wherewhere are in accordance with the rotation of and meet the right-handed spiral rule.

The displacement field can be expressed as

The values of the 16 coefficients can be determined from the displacement of the four vertices and the gradient of the three directions.

After introducing a shape function N, the displacement field in three dimensions is expressed aswhere

Strain and strain gradients take the following form:where L1 and L2 are differential operators that can be written aswhere

Therefore, the geometric equation considering the strain gradient can be expressed aswhere B is the strain matrix.

Finally, the element stiffness matrix of the order tetrahedral element is

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:where is the element mass matrix; , is the density of the rock mass; is the element damping matrix, , c is damping coefficient; is the element stiffness matrix; ,, and are the acceleration, velocity, and displacement matrices at any point in the element; and is the external loads.

Displacement at any point in the element is , where is a function only of position and is not related to time t, and is related to time and can be expressed as . Therefore, the acceleration and velocity at any point in the element are expressed, respectively, as and where · indicates a derivative with respect to time.

In the Newmark- method, when the control parameters satisfy , , 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]

Equation (24) can be expanded as

After subsisting equation (25) into equation (23) at time , the expression can be arranged aswhere is the equivalent stiffness matrix, ; is the equivalent load, .

In iterative computation, the process of each time step are as follows:(1)Calculate the equivalent load according to the operation state and external load of the node.(2)Calculate the equivalent stiffness matrix according to the node.(3)Solve equation set (24). When the structure is in the elastic state, the equivalent stiffness matrix 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 is related to the displacement of the node and must be solved through displacement iteration.(4)Calculate the acceleration and velocity at time 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 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 matrix of dynamic calculation is made up of mass matrix , damping matrix , and stiffness matrix . The plastic stiffness matrix is formed from 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:

In the calculation, the elastic stiffness matrix remains unchanged, and the plastic stiffness matrix 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 solved by the elastic equivalent stiffness matrix is taken as the initial displacement of the iterative calculation. The displacement value 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 iswhere 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 load (2)The elastoplastic stress field caused by the transient unloading of in situ stress

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 (Figure 7). is the equivalent peak load of the blasting load on the excavation surface and and 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.

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 on the excavation surface is calculated aswhere is the density of explosives, is the detonation velocity of explosives, is the isentropic exponent of explosives, which is generally taken as 3, is the diameter of the cartridge, is the diameter of the blast hole, and is the distance of adjacent blast holes.

When the explosion wave is propagating to the blasting hole with velocity , the blasting load in the hole rises to the maximum, and the rising time can be written as

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 load is calculated aswhere and are the respective lengths of the charging section and blockage section in the blast hole; is the average speed of expansion driven by crack explosion loading; is the unloading wave speed of gas explosion; and is the reflected wave velocity spread from the bottom to the top of the blast hole [32].

The time history relation of the blasting load acting on the excavation surface can be expressed as follows:

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 , 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. 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).

The time history relation of the in situ stress on the excavation surface can be written as

3.4.4. Loading Method

In the process of blasting excavation of an underground tunnel, the external loads imposed on surrounding rock elements are mainly divided into blasting load , gravity loads , and excavation loads [34].

According to the stress generated by blasting on rock mass elements, the blasting load can be obtained aswhere is the vertical gravity loads caused by the removed elements and can be expressed aswhere is the density of the rock mass, is the acceleration of gravity, and is the shape function matrix which can be obtained by Hermit interpolation function. , the transient unloading of in situ stress, can be expressed aswhere is strain matrix and 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 aswhere is the maximum tensile strain of a calculation element [35] and is the ultimate tensile strain, which is obtained by a uniaxial tensile test or Brazilian disk test.

If , 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 is given to the element where the tensile failure occurs; usually, . 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

Figure 9 shows the simplified stress-strain relation curve of rock under high stress. Elastic strain energy density , critical strain energy density , and ultimate strain energy density can be calculated as [36]

The damage failure criterion of zonal disintegration is judged by the change of the strain energy density.(1)When (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 , the element enters into the plastic damage stage; microcracks in the rock mass appear and begin to expand.(3)When (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. This is called the plastic residual stage.(4)When ( 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. Figure 10 shows the flowchart of the zonal disintegration procedure.

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 (Figure 11) [37, 38].

The average ranges of the fractured zones are listed in Table 1. The inner and outer diameters are the sides near and away from the tunnel in the distance, respectively.

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 is P0 = 23.5 MPa; the rock density  = 2.45 g/cm3; the compressive strength  = 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  = 0.725 10−3; and the ultimate strain  = 4.65 10−3.

Figure 12 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  = 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  = 35 mm, the length of the charging section is  = 1.4 m; the length of the blockage section is  = 0.6 m; the density of the explosive is  = 0.95–1.25 g/cm3; the detonation velocity is  = 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  = 52.75 MPa.

It is concluded that the unloading wave speed of gas explosion is approximately equal to the reflected wave velocity , that is,  =  = (1–1.5) × 103 m/s. The average speed of expansion driven by crack explosion loading is  = 0.25 , where is the longitudinal velocity of rock mass,  = 4550 × 103 m/s. According to equations (30) and (31), it can be obtained that the rise time of the blasting load is  = 470 , the duration of the blasting load is  = 5700, and the unloading time of in situ stress is  = 2130 .

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 .

4.3. Establishment of the ABAQUS Numerical Model

Figure 13 shows the three-dimensional computational mesh for the selected calculation range. There 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 [3941] 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

Figure 14 shows the numerical nephogram obtained by numerical simulation.

At the time of 250 , 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 , 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 , 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 , 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 , 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 time , 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 , 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 time , 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.

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.

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

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.

Data Availability

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

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This work was financially supported by the Key Development Program for Research of Shandong Province, under Grant/Award no. 2018GNC110023, and National Natural Science Foundation of China, under Grant/Award no. 51574156. The authors are deeply grateful for the support.