Blast damage zone strength reduction method for deep cavern excavation and its application

Abstract

The drill and blast (D&B) method is widely used to excavate underground spaces, but explosions generally cause damage to the rock. Still, no blast simulation method can provide computational accuracy and efficiency. In this paper, a blast equivalent simulation method called the blast damage zone strength reduction (BDZSR) method is proposed. This method first calculates the range of the blast-induced damage zone (BDZ) by formulae, then reduces the strength and deformation parameters of the rock within the BDZ ahead of excavation, and finally calculates the excavation damage zone (EDZ) for the D&B method by numerical simulation. This method combines stress wave attenuation, rock damage criteria and stress path variation to derive the BDZ depth calculation formulae. The formulae consider the initial geo-stress, and the reliability is verified by numerical simulations. The calculation of BDZ depth with these formulae allows the corresponding numerical simulation to avoid the time-consuming dynamic calculation process, thus greatly enhancing the calculation efficiency. The method was applied to the excavation in Jinping Class II hydropower station to verify its feasibility. The results show that the BDZSR method can be applied to blast simulation of underground caverns and provide a new way to study blast-induced damage.

Graphical abstract

Keywords

attenuation of stress wave / deep-buried underground building / drill and blast / in situ stress

1 Introduction

Underground buildings are important to different types of production and transportation and are mainly used for infrastructures such as underground parking, metros, railway tunnels, mining, hydropower station diversion tunnels and powerhouses. The most important characteristic of underground buildings is that the rock is its natural building material. The rock provides support (bearing capacity) for underground buildings, so the quality of the rock plays a decisive role in the stability of underground construction. The drill and blast (D&B) method is an economical and efficient method of underground excavation that is widely used in practical engineering. However, underground blasting usually leads to cracks in the rock and thus deteriorates the physical and mechanical properties of the rock and even creates a larger excavation damage zone (EDZ) in the rock mass.

Rock in underground engineering is subjected to geo-stress loading, especially in the case of deeper underground engineering, where the magnitude of in situ stress is generally higher. It is difficult to consider the initial stress of the rock in laboratory blasting experiments, so blast-induced damage is generally studied by numerical simulation methods. Li et al. [1] simulated the bursting process of a single blasthole with the aid of the Jones-Wilkins-Lee equation of state for explosion in LS-DYNA software and investigated wave propagation and blast-induced damage in a rock mass. The method provides a good simulation of the near-zone response of the explosion and has widely been used to study blast damage mechanisms [2–5]. However, this method applicable to only single-hole blasting scenarios is difficult to apply to numerical calculations for large and complex numerical simulations. During excavation in an underground space, the EDZ is formed primarily due to in situ stress redistribution and blast loading, especially for deeply buried underground engineering with a complex distribution of in situ stress [6,7]. To predict the damage zone of blasting excavation, Hu et al. [8] modeled the blasthole in a numerical simulation of the slope and applied the time-history curve of an explosion load on the blasthole. However, this method creates a large number of meshes for calculation, which makes it difficult to extend its application to numerical simulations of large and complex underground engineering projects. To overcome this limitation, many scholars [9–11] have applied the equivalent load of an explosion to a wider boundary, such as the excavation face of a cavern, avoiding the need to model each blasthole individually. The method can simulate the medium- and far-field responses of explosions very well and can be applied to numerical simulations of blasting in large-scale underground buildings. Although the problem of huge number of meshes can be solved by the above method, time-consuming dynamic calculations remain necessary, and this wastes considerable computational resources.

Additionally, some scholars have tried to calculate the depth of a blast-induced damage zone (BDZ) directly by deriving a formula. Leng et al. [12] derived a formula for the radius of the crush zone of rock blasted under column charge detonation conditions based on elastoplastic mechanics theory and found that the range of the crush zone is usually 1.2 to 5.0 times the radius of the blasthole. Wei et al. [13] derived the relationship between the radii of different BDZs and a blasthole based on the rock strain rate. Dai [14] combined the von Mises criterion and stress wave attenuation theory to derive the depth of the rock compressive damage zone (CDZ) and fracture zone due to blasting. However, these researchers ignored the effects of in situ stress and lateral stress coefficients [15].

In addition to the above studies, some authors have focused their research on the propagation and attenuation of blast-induced stress waves in rock masses [16]. Chen et al. [17] studied the blasting dynamic response of underground caverns by modeling the explosion wave propagation in rock masses with the universal distinct element code (UDEC). Hao et al. [18] investigated the attenuation of stress waves in jointed rock and its relationship with the direction of the joint. Wang and Li [19] studied the propagation of blast stress waves under uncoupled charging by outdoor blasting experiments and obtained an exponential relationship between the blast stress wave amplitude and distance of the wave front from the blasting source in the case of four different uncoupled coefficients. Deng et al. [20] proposed the UDEC-AUTODYN method to simulate the blast process, and the generated shock waves, to study the blast wave effects on underground structures and stress wave attenuation at different distances from the blasthole. Zhu et al. [21] studied the propagation of blast-induced waves in jointed rock and soil cover by the finite element method-discrete element method (FEM-DEM) and found that the ground soil cover can significantly attenuate blast-induced waves. Xie et al. [22] proposed the conceptualization of engineering-disturbed rock dynamics and investigated wave propagation, attenuation and superposition in rock. Zhu et al. [23] simulated the blasting process with UDEC and found that the burial depth of underground openings could significantly influence the propagation of stress waves and dynamic responses of underground openings under dynamic extension. The relationship between the attenuation of the stress wave and the distance from the blasthole has been validated, and in turn, the depth of a BDZ can be calculated by the attenuation of the stress wave.

Considering the above limitations of numerical blasting simulation and calculation formulae, a new blast simulation method called blast damage zone strength reduction (BDZSR) is proposed, which does not require dynamic calculation in the numerical simulation and takes into account the strength deterioration of the damaged rock. In this method, the rock damage process caused by blast loading and that caused by in situ stress redistribution are viewed as two separate processes [24]. In other words, the EDZ consists of the BDZ and the stress redistribution-induced damage zone (SRDZ). First, a new research approach combining attenuation of stress wave propagation theory, rock failure criteria and stress path variation is proposed to derive formulae for BDZ depth calculation, so that the designers can quickly determine the approximate depth of rock damage caused by different blasting schemes on site. Secondly, damage variables are introduced to simulate the reduction in the physical and mechanical properties of the rock in the BDZ. Finally, the distribution of the EDZ due to the D&B method for underground engineering is calculated by numerical simulation software, and the method is applied to calculate the extent of the damage zone for the D&B excavation of the auxiliary cavern in Jinping II hydropower station. A comparison of these results with data measured in the field and determined with the equivalent load method shows that the BDZSR method has good performance in terms of calculation efficiency and accuracy, providing an important reference for the design of blasting and support schemes.

2 Simulation method of blast load and rock damage criteria

2.1 Blast load history and peak pressure

During the blasting process, a blast-induced stress wave acts on the surrounding rock across the diameter of the blasthole: the surrounding rock is compacted, and the stress wave decays rapidly after reaching the peak value [25]. Therefore, the blasting load in geotechnical engineering numerical simulation is generally simplified into a triangular load form [26,27], as shown in Fig.1. P₀ is the peak of the explosion load, t_r is the rising time of the explosion load, and t_d is the duration of the explosion load.

Under the conditions of coupled charging, the rock near the blasthole is subjected to impact load due to the explosion of the column charge. According to acoustic principles, the explosion load p₀ is

(1)p_0=\frac{1}{1+\gamma }{\rho }_0{V}_{\mathrm{D}}^2,

where ρ₀ is the explosive density of the column charge, kg/m³; V_D is the velocity of detonation, m/s; and γ is the ratio of the specific heats of the detonation gases, usually taken as γ = 3.

In the case of uncoupled charging, the explosion load p is

(2)p=\frac{1}{2}{p}_0{\left(\frac{d_{\mathrm{b}}}{d_{\mathrm{c}}}\right)}^{-2\gamma }{l}_{\mathrm{e}}n,

where d_b and d_c are the radii of the blasthole and the charge, respectively, mm; l_e is the coefficient of axial charges; and n is the pressure increase factor, usually 8–11.

2.2 Rock damage criteria

The damage variables defined by existing damage criteria are mostly functions of crack density, which is mainly determined through probability distributions, making it difficult to apply to numerical simulations using models with complex in situ stress conditions. This paper defines a damage model for blasting processes based on the Mohr‒Coulomb criterion, combined with the first strength criterion and the maximum compressive stress criterion, including both tensile and compressive intercepts, which can be expressed in the (σ₁, σ₃) plane, as shown in Fig.2. The Mohr‒Coulomb criterion is selected because it is widely recognized and the parameters are easy to calibrate. This model is coded with FLAC^3D’s built-in FISH language. Therefore, the plastic damage state can be divided into three domains, i.e., domain 1 (compressive damage), domain 2 (shear damage) and domain 3 (tensile damage). FLAC^3D is robust enough to handle complex operations of a large numerical model. It has a wide range of commonly used built-in constitutive models and can carry out dynamic calculations. Blast simulation research based on FLAC^3D can provide the basis for subsequent application to larger numerical simulation models.

The damage criterion from point A to point B in Fig.2 is based on the maximum compressive stress criterion f_c = 0.

(3)f_{\mathrm{c}}={\sigma }_1-{\sigma }_{\mathrm{c}},

where σ₁ is the maximum principal stress and σ_c is the compressive strength of the rock.

The damage criterion from point B to point C is based on the Mohr‒Coulomb criterion f_s = 0.

(4)f_{\mathrm{s}}={\sigma }_1-{\sigma }_3{N}_{\text{φ}}+2c\sqrt{N_{\text{φ}}},

where N_{\text{φ}} is:

(5)N_{\text{φ}}=(1+\sin \phi )/(1-\sin \phi ),

where σ₃ is the minimum principal stress, φ is the internal friction angle, and c is the cohesion.

From point C to point D the damage criterion is based on the tensile strength criterion f_t = 0.

(6)f_{\mathrm{t}}={\sigma }_3-{\sigma }_{\mathrm{t}},

where σ_t is the tensile strength of the rock.

The potential function can be represented by three functions, namely, g^c, g^s, and g^t, where g^c is the compressive plastic flow. The function g^c corresponds to the associated law and is written as:

(7)g^{\mathrm{c}}={\sigma }_1.

The function g^s corresponds to a non-associated law and has the form:

(8)g^{\mathrm{s}}={\sigma }_1-{\sigma }_3{N}_{\text{ψ}},

where N_ψ is:

(9)N_{\text{ψ}}=(1+\sin \psi )/(1-\sin \psi ),

where ψ is the dilation angle.

The function g^{\mathrm{t}} corresponds to an associated law and is written as follows:

(10)g^{\mathrm{t}}={\sigma }_{\mathrm{t}}.

h₁(σ₁, σ₃) = 0 is defined as the value of σ₁ on the Mohr‒Coulomb envelope in the (σ₁, σ₃) plane when σ₁ = σ_c. h₂(σ₁, σ₃) = 0 is defined as the diagonal line between f_s = 0 and f_t = 0 in the plane of (σ₁, σ₃). The plastic damage state is determined by the positive and negative domains of the discriminant function, as shown in Fig.2, with the following expressions:

(11)h_1=\left({\sigma }_{\mathrm{c}}+2c\sqrt{N_{\text{φ}}}\right)/N_{\text{φ}},

(12)h_2={\sigma }_3-{\sigma }_{\mathrm{t}}+a^{\mathrm{p}}\left({\sigma }_1-{\sigma }^{\mathrm{p}}\right),

where a^p and σ^p are constants that are defined as:

(13)a^{\mathrm{p}}=\sqrt{1+N_{\text{φ}}^2}+N_{\text{φ}},{\sigma }^{\mathrm{p}}={\sigma }_{\mathrm{t}}{N}_{\text{φ}}-2c\sqrt{N_{\text{φ}}}.

When a point breaks the above yield criterion, the (σ₁, σ₃) plane in Fig.2 indicates that the elastic point falls within domains 1, 2 or 3. If the point falls within domain 1, i.e., when f_c < 0 and h₁ < 0, then the point is in a compressive damage state, for which Eq. (7) applies. If the point falls within domain 2, i.e., when h₁ > 0, when f_s < 0 and when h₂ < 0, then the point is in a shear damage state, and Eq. (8) applies. If the point falls within domain 3, i.e., when h₂ > 0 and f_t < 0, then the point is in a tensile damage state, and Eq. (10) applies.

3 Blast damage zone strength reduction method and application

3.1 Distribution of the blast-induced damage zone and numerical calculation model

3.1.1 Distribution of the blast-induced damage zone

After the explosion of the charge in the rock, the strength of the blast-induced stress wave exceeds the compressive strength of the rock by several times, and the rock near the blasthole is thus crushed. As the propagation distance increases, the blast-induced stress wave gradually decreases in strength to the level of the shear and tensile strength of the rock, causing many radial cracks in the rock outside the crush zone, under the influence of shear and tensile stress [28,29]. Existing blasting theory generally divides zones into a crush zone, a fracture zone and an elastic vibration zone [30]. Combined with the blast strength criterion presented in Subsection 2.2, in this study the blast damage distribution is further divided into a CDZ, a shear damage zone (SDZ), a tensile damage zone (TDZ) and an elastic zone, as shown in Tab.1 and Fig.3.

3.1.2 Numerical calculation model

To investigate the effects of blasting on the surrounding rock, a numerical model was introduced to aid the study. The size of the calculation model is 10 m × 10 m, with a blasthole diameter of 40 mm and a charge diameter of 32 mm. The blasthole is located in the center of the model, as shown in Fig.4.

The mechanical parameters used for the numerical calculation are shown in Tab.2. The uncoupling factor of the blasthole is 1.25, the explosive is emulsified explosive, the density is 1000 kg/m³, the velocity of detonation is 3600 m/s and the coefficient of axial charges is 1. Then, according to Eqs. (1) and (2), the average detonation pressure is 1620 MPa, and the blasthole wall load is 424.7 MPa. The time-history function for the blast load on the blasthole wall is chosen as a triangular function, as shown in Fig.1, where the rise time t_r = 0.5 ms, the positive pressure time t_d = 5.0 ms and the peak pressure is 424.7 MPa.

3.2 Calculation of the blast-induced damage zone considering in situ stress

3.2.1 Effect of in situ stress on the depth of the blast-induced damage zone

Without considering the effect of in situ stress, the starting point of the stress path at any elastic point in Fig.2 is (σ₁, σ₃) = (0, 0). Considering the effect of in situ stress, assuming in situ stress is σ₀ (for underground engineering, σ₀ is negative), the starting point of the stress path at any elastic point becomes (σ₁, σ₃) = (σ₀, σ₀), as shown in Fig.5, and the starting point of any elastic point moves from point O to point E. The elastic point moves in a direction away from the damage criterion line ABCD, compared to the starting point of the stress path, without considering in situ stress. This indicates that in blasting engineering, in situ stress enhances the yield limit of the stress state to some extent and that stress waves generated by explosions need to overcome in situ stress to propagate. This phenomenon has also been confirmed by many scholars [31–33], with different research methods.

3.2.2 Depth of compressive damage zone

As seen from Subsection 2.2, compressive damage occurs when any elastic point in the numerical simulation falls into domain 1, so the most important foundation of calculating the BDZ is to investigate the relationship between the maximum principal stress and the minimum principal stress at any elastic point under blast loading. The blasting problem in rock can be reduced to an axisymmetric plane problem in polar coordinates, and the general expression in polar coordinates for axisymmetric stress is shown below:

(14)\{\begin{array}{rl}& {\sigma }_{\mathrm{r}}=\frac{A}{r^2}+B(1+2\ln r)+2C,\\ & {\sigma }_{\text{θ}}=-\frac{A}{r^2}+B(3+2\ln r)+2C,\\ & {\tau }_{\mathrm{r}\text{θ}}={\tau }_{\text{θ}\mathrm{r}}=0.\end{array}

Because the ring is connected at multiple locations, it can be determined from the displacement singular value condition that B = 0. At the inner wall of the blasthole, i.e., when r = r_b, the value of the radial stress σ_p caused by explosion is p. In the process of increasing radial stress, any elastic point in the CDZ first reaches the Mohr‒Coulomb envelope for plastic damage, and then compressive plastic damage occurs as the maximum principal stress gradually increases until it reaches compressive strength; this can be expressed in Fig.5 as the point moving in the direction of envelope CB and finally reaching point B, where compressive damage occurs. Therefore, at r = R₁, Eq. (4) is only just applicable. Adding the above conditions into Eq. (14), the undetermined coefficient can be solved.

(15)\{\begin{array}{rl}& A=\frac{\left(N_{\text{φ}}-1\right)r_{\mathrm{b}}^2{R}_1^2}{R_1^2\left(N_{\text{φ}}-1\right)+\left(N_{\text{φ}}+1\right)r_{\mathrm{b}}^2}\left(p-\frac{2c\sqrt{N_{\text{φ}}}}{N_{\text{φ}}-1}\right),\\ & 2C=\frac{p\left(N_{\text{φ}}+1\right)r_{\mathrm{b}}^2+2c\sqrt{N_{\text{φ}}}{R}_1^2}{R_1^2\left(N_{\text{φ}}-1\right)+\left(N_{\text{φ}}+1\right)r_{\mathrm{b}}^2}.\end{array}

At the border of the CDZ, i.e., R = R₁, the radial stress σ_pr caused by the explosion is equal to the dynamic compressive strength [σ_c], i.e., σ_pr = [σ_c]. Combining the above condition with Eqs. (14) and (15), the CDZ depth R₁ can be expressed as

(16)R_1=\sqrt{\frac{2p{N}_{\text{φ}}-\left[{\sigma }_{\mathrm{c}}\right]\left(N_{\text{φ}}+1\right)-2c\sqrt{N_{\text{φ}}}}{\left[{\sigma }_{\mathrm{c}}\right]\left(N_{\text{φ}}-1\right)-2c\sqrt{N_{\text{φ}}}}}{r}_{\mathrm{b}},

where the dynamic compressive strength [σ_c] can be calculated by the loading strain rate \dot{\epsilon }. According to available papers, the dynamic compressive strength [σ_c] is proportional to the loading strain rate \dot{\epsilon }^1/3 [13]:

(17)\left[{\sigma }_{\mathrm{c}}\right]={\sigma }_{\mathrm{c}}{\dot{\epsilon }}^{1/3},

where \dot{\epsilon } is the loading strain rate, s⁻¹. The loading strain rate is usually taken as 10⁰–10⁴ s⁻¹.

However, Eq. (16) does not take into account the effect of in situ stress. In the CDZ, the rock first undergoes shear damage, followed by compressive damage as the maximum principal stress gradually increases. This indicates that any elastic point in the case of in situ stress needs to pass farther through the stress path in the (σ₁, σ₃) plane for the maximum principal stress to reach dynamic compressive strength. Therefore, Eq. (16) can be improved by combining numerical simulation experiments. Numerical simulations of the depth of the damage zone were carried out using different working conditions with a vertical in situ stress of 20 MPa and lateral stress coefficients of 0.7–1.0, and the results are shown in Fig.6, where the numerical calculation model parameters are shown in Subsection 3.1.

As shown in Fig.6, at the same maximum initial stress, different lateral stress coefficients have almost no effect on the depth of the CDZ. The result of the numerical calculation is 0.1072 m, which is less than the 0.1313 m calculated by Eq. (16). Because the stress wave generated by the explosion needs to overcome in situ stress to propagate through the rock, in situ stress increases the compressive strength of the rock. Based on this idea and combined with the results of numerical calculation, Eq. (16) can be improved. When considering the effect of in situ stress, the depth of the CDZ is

(18)R_1=\sqrt{\frac{2p{N}_{\text{φ}}-\left[{\sigma }_{\mathrm{c}}\ast \right]\left(N_{\text{φ}}+1\right)-2c\sqrt{N_{\text{φ}}}}{\left[{\sigma }_{\mathrm{c}}\ast \right]\left(N_{\text{φ}}-1\right)-2c\sqrt{N_{\text{φ}}}}}{r}_{\mathrm{b}},

where \left[{\sigma }_{\mathrm{c}}\ast \right] is the equivalent dynamic compressive strength under the effect of in situ stress, defined as:

(19)\left[{\sigma }_{\mathrm{c}}\ast \right]=\left({\sigma }_{\mathrm{c}}+{\sigma }_{01}\right){\dot{\epsilon }}^{1/3},

where σ₀₁ is the initial maximum principal stress.

The depth of the CDZ calculated according to Eq. (18) is R₁ = 0.1098 m, which is closer to the actual depth calculated by numerical simulation. To verify the feasibility of the improved formula, four work conditions with initial maximum principal stresses of 10, 15, 25, and 30 MPa were selected for blasting simulation, and the comparison of the formula calculation results with the numerical simulation is shown in Tab.3.

It can be seen that Eq. (18) can be effectively applied to calculate the CDZ depth when considering the in situ stress and lateral stress coefficient.

3.2.3 Depth of shear damage zone

As shown in Fig.6, at the edge of the SDZ, r = R₂, the point stress state reaches the Mohr‒Coulomb criterion, and the maximum and minimum principal stresses at this point satisfy the relationship in Eq. (20).

(20){\sigma }_1-N_{\text{φ}}{\sigma }_3+2c\sqrt{N_{\text{φ}}}=0.

It is possible to determine the radial stress at the time of plastic damage by studying the variation in the stress path at any elastic point outside the CDZ under blast loading for different in situ stresses and lateral stress coefficients. From Fig.6, unlike the extent of the CDZ, the extent of the SDZ is strongly influenced by the lateral stress coefficient λ and cannot be reduced to an axisymmetric plane problem in polar coordinates. The λ (λ = σ₀₁/σ₀₃) is defined as the ratio of the initial maximum principal stress σ₀₁ to the initial minimum principal stress σ₀₃ in the cross-section of the blasthole. This shows that the pattern of the stress path is different under different lateral stress conditions. To verify the above hypothesis, numerical simulations were carried out for a single blasthole explosion considering a lateral pressure coefficient of 0.7–1.0 and an initial maximum principal stress of 20 MPa. The stress path was recorded at locations 0.15 m from the blasthole in the vertical direction, and the results are shown in Fig.7. The line AB in the figure is the compressive strength criterion, and BC is the shear strength criterion, corresponding to Eqs. (3) and (4), respectively.

As seen from Fig.7, for the same initial maximum principal stress but with different lateral stress coefficient λ, the change in stress paths at the same location is entirely different, and there is a large difference in the rate of change k close to the strength criterion line. The closer λ is to 1.0, the smaller the absolute value of the gradient of stress; in addition, the closer it is to 1.0, the more slowly σ₃ changes with σ₁. Based on the results of the numerical simulation, the gradient of stress k (k = Δσ₃/Δσ₁) for different lateral pressure coefficients prior to plastic damage can be defined as

(21)

The relationship between Δσ₃ and Δσ₁ for any elastic point outside the CDZ moving in the (σ₁, σ₃) plane toward the strength criterion line can be derived from Eq. (21), i.e., Δσ₃ = kΔσ₁. Then, any elastic point in the (σ₁, σ₃) plane taking into account in situ stress and lateral stress coefficients can be expressed as

(22)k\left({\sigma }_1-{\sigma }_{01}\right)={\sigma }_3-{\sigma }_{03},

where σ₀₁ is the initial maximum principal stress and σ₀₃ is the initial minimum principal stress.

Outside the CDZ, in situ stress can be seen as a reverse force preventing the propagation of the blast stress wave, so that in Eq. (22) σ₁ can be expressed as σ₁ = σ_pr ‒ σ₀₁. Equation (22) can be expressed as

(23){\sigma }_3=k\left({\sigma }_{\mathrm{p}\mathrm{r}}-2{\sigma }_{01}\right)+{\sigma }_{03}.

When in situ stress is not considered, the maximum principal stress at the boundary between the CDZ and the SDZ is the dynamic compressive strength [σ_c]. When considering in situ stress, it can be assumed that the in situ stress inhibits the propagation of the stress wave onward to the next damage zone, so that when considering propagation from the CDZ boundary to the next damage zone, the maximum principal stress value at the CDZ boundary is σ_R1:

(24){\sigma }_{\mathrm{R}1}=\frac{\left({\sigma }_c-{\sigma }_{01}\right){\dot{\epsilon }}^{1/3}}{b}.

At the boundary of the SDZ, the stress state at this location is exactly changed from the elastic state to the shear damage state. Thus, taking Eqs. (23) and (24) into the Mohr‒Coulomb strength criterion,

(25)\left({\sigma }_{\mathrm{p}\mathrm{r}}-{\sigma }_{01}\right)-k\left({\sigma }_{\mathrm{p}\mathrm{r}}-2{\sigma }_{01}\right)N_{\text{φ}}-{\sigma }_{03}{N}_{\text{φ}}+2c\sqrt{N_{\text{φ}}}=0,

where σ_pr is the radial stress caused by blasting.

The transmitted shock wave front in the rock propagates outward continuously, and the radial stress σ_pr caused by explosion attenuation is continuous. σ_pr at any point in the rock can be expressed as

(26){\sigma }_{\mathrm{p}\mathrm{r}}={{\sigma }_{\mathrm{r}}|}_{R=r}{\overline{r}}^{-\alpha },

where σ_r |_{R = r} is the radial stress at R = r; \overline{r} is the distance ratio, \overline{r}=r^{\prime }/r; r^{\prime } is the distance from any point to the center of the blasthole; and α is the load propagation attenuation index, which can be determined from the impedance of the rock. The equation for the attenuation of the stress wave in column charging expressed in terms of wave impedance is [34]:

(27)\alpha =-2.86\times {10}^{-7}\rho C_{\mathrm{p}}+1.89,

where ρC_p is the impedance of the rock, g/cm³·cm/s; ρ is the density of the rock; and C_p is the velocity of the rock in the longitudinal direction.

Inserting Eqs. (24) and (26) into Eq. (25) gives the depth of SDZ R₂ as

(28)R_2={\left(\frac{{\sigma }_{\mathrm{R}1}\left(1-k{N}_{\text{φ}}\right)}{\left(1-2k{N}_{\text{φ}}\right){\sigma }_{01}+{\sigma }_{03}{N}_{\text{φ}}-2c\sqrt{N_{\text{φ}}}}\right)}^{1/\alpha }{R}_1.

To verify the feasibility of Eq. (28), the numerical simulation was set at 20 working conditions with an initial maximum principal stress ranging from 10 to 30 MPa and a lateral stress coefficient ranging from 0.7 to 1.0. The numerical calculation results are shown in Fig.9‒Fig.13. The comparison of the calculation results of the two methods is listed in Tab.4, and the comparison of the relative errors is presented in Fig.8. The results indicate the following.

1) The relative error between the depth of the plastic damage zone calculated by the formulae and the numerical simulation is small, ranging from 1.4% to 15%, with an average relative error of 6.1%. This shows that the depth of the SDZ can be accurately calculated by Eq. (28) when the in situ stress and lateral pressure coefficients are taken into account.

2) The depth of the CDZ is hardly influenced by the variation in the lateral pressure coefficient, and its depth is mainly influenced by the initial maximum principal stress. The higher the initial maximum principal stress is, the smaller the extent of the CDZ.

3) Different lateral pressure coefficients affect the depth of the damage zone under the same initial maximum principal stress. The BDZ will develop in the direction of the maximum principal stress and perpendicularly to the axial direction of the blasthole, which is consistent with the conclusions of other scholars [35–39]. The depth of the damage zone perpendicular to the direction of the maximum principal stress is related to the lateral stress coefficient. Based on the numerical simulation results, the damage depth R_2^{\text{'}} perpendicular to the direction of the initial maximum principal stress can be expressed as:

(29)R_2^{\prime }=\frac{\lambda }{2-\lambda }{R}_2.

4) The TDZ was not observed in the results posted in Fig.9–Fig.13. This indicates that in deeply buried underground engineering where initial stress is considered, single-hole blasting of intact rock may not create a TDZ, i.e., the stress state of the zone in Fig.2 has difficulty falling into domain 3. If the initial stress is not considered, a TDZ is likely to be induced [28].

3.3 Parameter reducing method for damaged rock

With the formulae in Subsection 3.2, it is possible to calculate the depth of damage to the surrounding rock caused by column charge blasting for different blast loadings, for different maximum initial principal stresses, and for lateral stress coefficients, directly, without modeling the blasthole. After calculations by Eqs. (28) and (29), all zones within the BDZ can be identified as zones that are damaged by blasting. The blasting process lasts only a few milliseconds, which is almost negligible compared to the long-term deformation of the surrounding rock after excavation, so it can be considered that the damage of this area has already occurred before the static equilibrium calculation of the excavation in the numerical simulation.

After the damage of the rock occurs, significant deterioration of the rock mechanical parameters also occurs. In the blast damage model, the equivalent physical parameters of the surrounding rock in the damaged state are related to those in the undamaged state as follows:

(30)K^{\prime }=K(1-D),

where K^{\prime } and K are the equivalent bulk moduli of the damaged surrounding rock and the undamaged surrounding rock, respectively, and D is the damage variable.

At construction sites, acoustic testing is used to obtain the P-wave velocity V_pm of the rock mass, which can reveal the quality and intactness of the surrounding rock mass. The intact rock coefficient K_v can be expressed in terms of acoustic velocity according to the recommendation of the Standard for Engineering Classification of Rock Masses (GB/T 50218-2014) in China.

(31)K_{\mathrm{v}}={\left(\frac{V_{\mathrm{p}\mathrm{m}}}{V_{\mathrm{p}\mathrm{r}}}\right)}^2,

where V_pr is the intact rock wave velocity.

The damage variable D defined by the intact rock coefficient K_v can be expressed by the following [40]:

(32)D=1-K_{\mathrm{v}}.

When K_v is greater than 0.75, the rock mass can be considered to be intact rock. Therefore, in this work, the rock was considered to be damaged when D ≥ 0.25 and not damaged when D < 0.25 [41].

In some engineering, it is difficult to carry out acoustic testing on all rock masses, so an estimation method can be used to calculate the speed of acoustic wave propagation in the rock. The speed of sound wave propagation in rock is closely related to the mechanical properties of the surrounding rock, and the velocity of longitudinal wave propagation in rock can be expressed as [8]:

(33)V_{\mathrm{p}}=\sqrt{\frac{3K(1-v)}{\rho (1+v)}},

where V_p is the velocity of the longitudinal wave, ν is Poisson’s ratio, and ρ is the density of the rock.

According to Lemaitre’s hypothesis of strain equivalence, the strain induced by stress in damaged rock is equivalent to the strain induced by an equivalent stress in intact rock, which is expressed as follows [42–44]:

(34)\epsilon =\frac{{\sigma }^{\text{*}}}{E}=\frac{\sigma }{E^{\text{*}}}=\frac{\sigma }{E(1-D)},

where ε is the material strain, {\sigma }^{\text{*}} is the equivalent stress, σ is the apparent stress, E is the Young’s modulus of the undamaged rock and E^{\text{*}} is the Young’s modulus of the damaged rock. According to Eq. (34), the relationship between the equivalent stress and the apparent stress can be expressed as

(35){\sigma }_{ij}^{\ast }=\frac{{\sigma }_{ij}}{1-D}.

Taking Eq. (35) into Eq. (4) gives a new Mohr‒Coulomb damage criterion for the damage state, as follows:

(36)f_{\mathrm{s}}=\frac{{\sigma }_1}{(1-D)}-\frac{{\sigma }_3}{(1-D)}{N}_{\text{φ}}+2c\sqrt{N_{\text{φ}}}.

Equation (36) can be applied to rock damaged by blasting to simulate the effect of blasting on the mechanical properties of the surrounding rock.

3.4 Process of the blast damage zone strength reduction method

The process of the BDZSR method is shown in Fig.14. After inversion of in situ stress, the extent of BDZ is calculated by formulae in Subsection 3.2, then the mechanical properties of rock in BDZ are reduced by Eq. (36). Finally, the extent of the SRDZ is calculated by numerical simulation software. The detailed process is listed as follows.

Step 1: Inversion of in situ stress in the engineering area. Information on the magnitude and direction of the maximum principal stress in the zone near the blasthole is read.

Step 2: The depth of damage to the surrounding rock for each blasthole explosion is calculated according to Eqs. (28) and (29). The surrounding rock within the damage zone is the BDZ, which is defined as a separate group of zones called damage in FLAC^3D.

Step 3: The damage zones defined in Step 2 are set to the plastic damage state, and Eq. (36) is applied to the damaged zone to simulate the deterioration of the mechanical properties of the surrounding rock.

Step 4: Excavation and static balance calculations are carried out according to the design of the project and the form of the SRDZ. Then, Step 2 is repeated after each level of excavation calculations according to the designed excavation sequence, until the final level of excavation is completed.

3.5 Application and validation of the blast damage zone strength reduction method

3.5.1 Application of the blast damage zone strength reduction method

To verify the reliability of the BDZSR method, the auxiliary cavern of the deeply buried underground cavern group in the Jinping II Hydropower Station was selected as the research object. The auxiliary cavern was constructed by the D&B method, with a diameter of 5.50–6.25 m. The size of the auxiliary cavern used for the numerical simulation is 5.5 m × 5.7 m. The blastholes of the auxiliary cavern were 45 mm in diameter, with a 32 mm diameter charge and 500 mm spacing. The emulsion explosive used in the package had a density of 1000 kg/m³ and a burst speed of 3600 m/s. using Eq. (2), a blasthole wall load of 209.5 MPa was calculated. The numerical simulation model is shown in Fig.15, and the location of the blasthole is shown in Fig.16.

According to the data provided by Luo et al. [11], the in situ stress in the auxiliary cavern is 50.8 MPa in the direction of the cavern axis, 43.9 MPa in the horizontal direction and 38.5 MPa in the vertical direction. The rock mechanical parameters used in the numerical simulation model of the auxiliary cavern are shown in Tab.5. The rock mechanical parameters were justified in Ref. [11]. The location of the acoustic test hole for the auxiliary cavern is shown in Fig.17.

The depth of the BDZ of the surrounding rock around the blasthole under the influence of explosion was calculated by Step 2 in Subsection 3.4, as shown in Fig.18. According to Step 3, the BDZ is set to a plastic damage state, and the mechanical parameters of the rock in this area are reduced by Eq. (36).

The wave velocity reduction rate of 25% was taken from the field acoustic test provided by Luo et al. [11]. The auxiliary cavern was then excavated to obtain the damage zone caused by blasting, the result of which is shown in Fig.19(a). To compare the effect of blasting on the damage zone of the surrounding rock, the result of the cavern excavation without considering blasting is shown in Fig.19(b). The black points are the measured damage zone depth for each measurement hole and the red curve is obtained by connecting each measured point by interpolated curves. A comparison of the damage zone of cavern excavation considering blasting effects and ignoring blasting effects and the results of acoustic tests from the field are shown in Tab.6.

The extent of the damage zone is small when blasting is not considered, with an average depth of 1.16 m and an average relative error of 14.7% from the measured depth of the damage zone. The average depth of damage for the condition simulating the blast using the BDZSR method is 1.30 m, which is close to the measured values in the field, with a relative error of 8.9%. The depth of damage for measurement hole #6 in the top arch of the auxiliary cavern was 1.33 m, which is closer to the measured depth than the 1.13 m calculated without the consideration of blasting. The depth of damage calculated by the BDZSR method for measurement hole #1 below the auxiliary hole was 1.83 m, which is closer to the measured depth than the 1.70 m without consideration of blasting. In summary, the depth of the damage zone calculated by the BDZSR method is in good agreement with the measured values. This means that the BDZSR method can well simulate the damage to the surrounding rock by blasting.

3.5.2 Comparison of the blast damage zone strength reduction method with the equivalent load method

The equivalent load method is a widely used method in the simulation of blasting effects in geotechnical engineering. This method does not require modeling of the blastholes and simulates the damage to the surrounding rock by applying a discounted equivalent normal load between the blastholes, as shown in Fig.20. In the figure, P_e is the applied equivalent pressure, r_b is the radius of the borehole and a is the borehole spacing, which can be expressed as

(37)P_{\mathrm{e}}=\frac{2r_{\mathrm{b}}}{a}P,

where P is the transmitted wave pressure of the burst of gases in the rock, calculated by Eq. (2).

In this study, the equivalent load P_e applied on the excavation face is 18.9 MPa. The time curve for the blasting dynamics calculations is shown in Fig.4, and the results are shown in Fig.21.

A comparison of the results of the equivalent load method and the BDZSR method is shown in Tab.7, and a comparison of the calculation efficiency of the two methods when using the same computer is shown in Fig.22. The results indicate the following.

1) The simulation results of the BDZSR method at measurement holes #1, #2, #3, #5, #8, and #9 were all better than those of the equivalent load method in terms of the measured depth. The average depth of the BDZ calculated by the BDZSR method was 1.17 m, with a relative error of 8.9%, while the average depth of the damage zone calculated by the equivalent load method was 1.43 m, with a relative error of 19.2%.

2) Under the same computer configuration, the calculation time for the equivalent load method was 5127 s, while the total calculation time for the BDZSR method was 2388 s, a reduction in the calculation time of 53.4%. This indicates that the BDZSR has a faster calculation speed while ensuring the accuracy of BDZ depth prediction and can be applied to more complex large underground plant cavern groups.

4 Conclusions

A new blasting simulation method called BDZSR is proposed to solve current problems in underground engineering. These problems are that the calculation formula of blast-induced damage depth ignores the influence of in situ stress, and that most blasting simulation methods require long calculation times in the application of the D&B method. The BDZSR method is verified by practical engineering measurement data. The following conclusions can be drawn.

1) A new research method combining attenuation of stress wave propagation theory, rock failure criteria and stress path variation is proposed to investigate the BDZ. The formulae for calculating the depth of the BDZ are derived, which take into account the in situ stress and lateral pressure coefficient. These formulae allow a builder to determine the depth of the BDZ directly at the construction site and can be easily applied to secondary development in numerical simulation.

2) When considering in situ stress, the lateral pressure coefficient significantly influences the development of the SDZ but not the CDZ. The development of the CDZ is mainly controlled by the magnitude of the initial stress.

3) In contrast to most existing methods, the BDZSR method does not need to perform dynamic calculations in excavation simulations, thus greatly enhancing the calculation efficiency. The method is applied to the excavation simulation of the Jinping II auxiliary cavern, and it is found that the method improved the calculation accuracy by 10.3% and saved 53.4% of the calculation time compared with the traditional equivalent pressure load method.

4) The approach in this paper is based on rock damage criteria, ignoring the effect of intermediate principal stress σ₂ on rock damage. A true tri-axial stress damage criterion should be proposed in subsequent studies to predict the EDZ more accurately for the D&B method.

References

[1] Li H B, Xia X, Li J C, Zhao J, Liu B, Liu Y Q. Rock damage control in bedrock blasting excavation for a nuclear power plant. International Journal of Rock Mechanics and Mining Sciences, 2011, 48(2): 210–218

[2] Zhu W C, Wei J, Zhao J, Niu L L. 2D numerical simulation on excavation damaged zone induced by dynamic stress redistribution. Tunnelling and Underground Space Technology, 2014, 43: 315–326

[3] Xie L X, Lu W B, Zhang Q B, Jiang Q H, Wang G H, Zhao J. Damage evolution mechanisms of rock in deep tunnels induced by cut blasting. Tunnelling and Underground Space Technology, 2016, 58: 257–270

[4] Chen M, Ye Z W, Lu W B, Wei D, Yan P. An improved method for calculating the peak explosion pressure on the borehole wall in decoupling charge blasting. International Journal of Impact Engineering, 2020, 146: 103695

[5] Hu Y G, Yang Z W, Huang S L, Lu W B, Zhao G. A new safety control method of blasting excavation in high rock slope with joints. Rock Mechanics and Rock Engineering, 2020, 53(7): 3015–3029

[6] Diederichs M S, Kaiser P K, Eberhardt E. Damage initiation and propagation in hard rock during tunnelling and the influence of near-face stress rotation. International Journal of Rock Mechanics and Mining Sciences, 2004, 41(5): 785–812

[7] Golshani A, Oda M, Okui Y, Takemura T, Munkhtogoo E. Numerical simulation of the excavation damaged zone around an opening in brittle rock. International Journal of Rock Mechanics and Mining Sciences, 2007, 44(6): 835–845

[8] HuY GLuW BJinX HChenMYanP. Numerical simulation for excavation blasting dynamic damage of rock high slope. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(11): 2204−2213 (in Chinese)

[9] Lu W B, Yang J H, Chen M, Zhou C B. An equivalent method for blasting vibration simulation. Simulation Modelling Practice and Theory, 2011, 19(9): 2050–2062

[10] Lu W, Yang J, Yan P, Chen M, Zhou C, Luo Y, Jin L. Dynamic response of rock mass induced by the transient release of in-situ stress. International Journal of Rock Mechanics and Mining Sciences, 2012, 53: 129–141

[11] LuoSYanPLuW BChenMWangG H. Research on the simulation of blasting damage and its mechanism of deep tunnel excavation. Chinese Journal of Rock Mechanics and Engineering, 2021, 40(Sup 1): 2760−2772 (in Chinese)

[12] LengZ DLuW BChenMYanPHuY G. Improved calculation model for the size of crushed zone around blasthole. Explosion and Shock Waves, 2015, 35(1): 101−107 (in Chinese)

[13] WeiDChenMYeZ WLuW BLiT. study on blasting failure zone based on rate-dependent dynamic characteristics of rock mass. Advanced Engineering Sciences, 2021, 53(1): 67−74 (in Chinese)

[14] DaiJ. Dynamic Behaviors and Blasting Theory of Rock. Beijing: Metallurgical Industry Press, 2002 (in Chinese)

[15] Luo S, Yan P, Lu W B, Chen M, Wang G H, Lu A, Liu X. Effects of in-situ stress on blasting damage during deep tunnel excavation. Arabian Journal for Science and Engineering, 2021, 46(11): 11447–11458

[16] Zhao J, Zhou Y X, Hefny A M, Cai J G, Chen S G, Li H B, Liu J F, Jain M, Foo S T, Seah C C. Rock dynamics research related to cavern development for Ammunition storage. Tunnelling and Underground Space Technology, 1999, 14(4): 513–526

[17] Chen S G, Cai J G, Zhao J, Zhou Y X. Discrete element modelling of an underground explosion in a jointed rock mass. Geotechnical and Geological Engineering, 2000, 18(2): 59–78

[18] Hao H, Wu Y K, Ma G W, Zhou Y X. Characteristics of surface ground motions induced by blasts in jointed rock mass. Soil Dynamics and Earthquake Engineering, 2001, 21(2): 85–98

[19] WangWLiX C. Experimental study of propagation law of explosive stress wave under condition of decouple charge. Rock and Soil Mechanics, 2010, 31(6): 1723−1728 (in Chinese)

[20] Deng X F, Chen S G, Zhu J B, Zhou Y X, Zhao Z Y, Zhao J. UDEC–AUTODYN Hybrid modeling of a large-scale underground explosion test. Rock Mechanics and Rock Engineering, 2015, 48(2): 737–747

[21] Zhu J B, Li Y S, Wu S Y, Zhang R, Ren L. Decoupled explosion in an underground opening and dynamic responses of surrounding rock masses and structures and induced ground motions: A FEM-DEM numerical study. Tunnelling and Underground Space Technology, 2018, 82: 442–454

[22] Xie H P, Zhu J B, Zhou T, Zhang K, Zhou C T. Conceptualization and preliminary study of engineering disturbed rock dynamics. Geomechanics and Geophysics for Geo-Energy and Geo-Resources, 2020, 6: 1–34

[23] Zhu J B, Li Y S, Peng Q, Deng X F, Gao M Z, Zhang J G. Stress wave propagation across jointed rock mass under dynamic extension and its effect on dynamic response and supporting of underground opening. Tunnelling and Underground Space Technology, 2021, 108: 103648

[24] Yan P, Lu W B, Chen M, Hu Y G, Zhou C B, Wu X X. Contributions of in-situ stress transient redistribution to blasting excavation damage zone of deep tunnels. Rock Mechanics and Rock Engineering, 2015, 48(2): 715–726

[25] GuoDXiaoW HGuoDLuY. Numerical simulation of surface vibration propagation in tunnel blasting. Mathematical Problems in Engineering, 2022: 3748802

[26] YangJ HLuW BChenM. Determination of the blasting load variation in borehole. In: Proceedings of the 2nd National Academic Conference on Engineering Safety and Protection. Beijing: Chinese Society for Rock Mechanics & Engineering, 2010, 773–777

[27] Zhang Y B, Yang F, Huang L Y. Study on the dynamic responses of tunnel linings with small net distance under bias pressures caused by blasting. Geotechnical and Geological Engineering, 2022, 40(2): 605–615

[28] Li X D, Liu K W, Yang J C, Song R T. Numerical study on blast-induced fragmentation in deep rock mass. International Journal of Impact Engineering, 2022, 170: 104367

[29] Lu W B, Leng Z D, Chen M, Yan P, Hu Y. A modified model to calculate the size of the crushed zone around a blast-hole. The Journal of The Southern African Institute of Mining and Metallurgy, 2016, 116(5): 413–422

[30] Esen S, Onederra I, Bilgin H A. Modelling the size of the crushed zone around a blasthole. International Journal of Rock Mechanics and Mining Sciences, 2003, 40(4): 485–495

[31] Aydan Ö. In situ stress inference from damage around blasted holes. Geosystem Engineering, 2013, 16(1): 83–91

[32] Yang J H, Yao C, Jiang Q H, Lu W B, Jiang S H. 2D numerical analysis of rock damage induced by dynamic in-situ stress redistribution and blast loading in underground blasting excavation. Tunnelling and Underground Space Technology, 2017, 70: 221–232

[33] Yilmaz O, Unlu T. Three dimensional numerical rock damage analysis under blasting load. Tunnelling and Underground Space Technology, 2013, 38: 266–278

[34] ZhangJ H. Study on the attenuation law of explosion stress wave in rock of cylinder charge. Thesis for the Master’s Degree. Taiyuan: North University of China, 2005 (in Chinese)

[35] BaiY. Blasting damage model and numerical test of rock under effect of in situ stress. Dissertation for the Doctoral Degree. Shenyang: Northeastern University, 2014 (in Chinese)

[36] BaiYZhuW CWeiC HWeiJ. Numerical simulation on two-hole blasting under different in-situ stress conditions. Rock and Soil Mechanics, 2013, 34(Sup 1): 466−471 (in Chinese)

[37] ChenMHuY GLuW BYanPZhouC B. Numerical simulation of blasting excavation induced damage to deep tunnel. Rock and Soil Mechanics, 2011, 32(5): 1531−1537 (in Chinese)

[38] MuZ MPanF. Numerical study on the damage of the coal under blasting loads coupled with geostatic stress. Chinese Journal of High Pressure Physics, 2013, 27(3): 403−410 (in Chinese)

[39] YangJ H. Coupling effect of blasting and transient release of in situ stress during deep rock mass excavation. Dissertation for the Doctoral Degree. Wuhan: Wuhan University, 2014 (in Chinese)

[40] Shen Y J, Yan R X, Yang G S, Xu G L, Wang S Y. Comparisons of evaluation factors and application effects of the new [BQ]GSI system with international rock mass classification systems. Geotechnical and Geological Engineering, 2017, 35(6): 2523–2548

[41] Zhang R, Xie H P, Ren L, Deng J H, Gao M Z, Feng G, Zhang Z T, Li X P, Tan Q. Excavation-induced structural deterioration of rock masses at different depths. Archives of Civil and Mechanical Engineering, 2022, 22(2): 81

[42] CaoW GZhaoM HLiuC X. A study on damage statistical strength theory for rock. Chinese Journal of Geotechnical Engineering, 2004, 26(6): 820−823 (in Chinese)

[43] CaoW GFangZ LTangX J. A study of statistical constitutive model for soft and damage rocks. Chinese Journal of Rock Mechanics and Engineering, 1998, 166(6): 628−633 (in Chinese)

[44] ZhangDLiuE LLiuX YZhangGYinXSongB T. A damage constitutive model for frozen sandy soils based on modified Mohr–Coulomb yield criterion. Chinese Journal of Rock Mechanics and Engineering, 2018, 37(4): 978−986 (in Chinese)

RIGHTS & PERMISSIONS

Higher Education Press