# Study on the Influences of Empty Hole and In Situ Stress on Blasting-Induced Rock Damage

Authors: Ningkang Meng; Jian-biao Bai; Wenlong Shen; Yong Chen; Xiangyu Wang; Wenda Wu; Bowen Wu
Journal: Geofluids
 (2021)

---

## Abstract

Blasting has been used extensively in deep underground excavation. In this paper, a numerical simulation LS-DYNA was carried out to investigate the influence of in situ stress and empty hole on rock fracture development. Based on the simulation process and existing laboratory data, the micromechanical properties of the Riedel-Hiermaier-Thoma (RHT) model were calibrated. By then, this model was used to study the effects of in situ stress and empty hole on the degree of damage and fracture development in the rock. Simulation results revealed that the difference between horizontal and vertical principal stresses can influence the degree of damage as well as the fracture initiation and development, which in turn leads to the fracture development towards a high-stress concentration area. Numerical simulation also showed that the total area of fracture development gradually decreases with the increasing in situ stress magnitudes. Under the same stress conditions, empty hole can facilitate the fracture development at 90°, while preventing the subsequent fracture development at 60° and 120°. This would result in further fracture propagation at 90°, such that the control of fracture direction can be achievable.


---

## Body

## 1. Introduction

As the mineral at shallow depth is rapidly extracted, it is inevitable to conduct mining activities at deeper location, including coal mines [1]. Comparing with shallow depth of cover, in situ stress magnitudes are substantially higher in deeper locations. Therefore, to effectively and efficiently extract coal underground, it is critical to consider the in situ stress magnitudes while implementing blasting and developing innovative techniques.D&B are not only used in coal mining but also in many fields, such as mining other mineral resources, slope treatment of opencast mines, tunnel excavation, and building demolition. Drill and blast is a traditional and effective method for coal extraction, which has been used widely in deep roadway development and rock fracturing. Prior to blasting, rock is subjected to in situ stress condition. For surface or shallow underground mining, the buried depth of coal and rock is small, and the thickness of overlying strata is relatively small. Therefore, for shallow coal and rock excavation projects, the in situ stress is low and can be neglected during blasting. However, in situ stress magnitudes in deep locations are more significant, which can be close or higher than the uniaxial compressive strength. Due to the high-stress conditions, the mechanical properties of rock at deeper locations are different than those at subsurface [2]. Lu et al. [3] found that the blasting-induced fracture development was heavily influenced by high in situ stress. They argued that when the magnitudes of in situ stress surpass 10–12 MPa, the blasting design for subsurface is unreasonable.Many researchers have investigated the fracture evolution and mechanism of rock under blasting. Based on the finite element method (FEM), Bendezu et al. [4] proposed an analysis technique to simulate the blasting-induced fracture development in rock. The study used three approaches to simulate rock fracturing process. Ma and An [5] used the 2D finite element method to investigate the effects of joint, loading rate, and in situ stress on blasting-induced rock fracture. Results revealed that the induced fracture developed towards the maximum principal stress direction. Jayasinghe et al. [6] investigated the effects of discontinuity persistence and high in situ stresses on the evolution of blasting-induced damage and revealed that discontinuity persistence and spatial distribution of rock bridges have a significant influence on the evolution of blasting-induced damage. Shang [7] studied the rupture of veined granite in polyaxial compression via a discrete element method model. Yang et al. [8] studied the effect of blasting sequencing on the damage zone of surrounding rock. Based on the simulation result through FLAC3D, the influences of loading rate on blasting performance and area were investigated. Zhao et al. [9] used discontinuous deformation analysis software to analyze the rock response and fracture process, which was informative on the influence of dynamic loading due to blasting on rock fracture and fracture development. Based on the randomness of rock fracture development from conventional blasting, some researchers suggested that the use of empty hole can reduce the degree of damage, while inducing direction-controlled fracture initiation and development [10–13]. The existence of empty hole can change the stress distribution and damage extend in the rock, which has been successfully implemented in preconditioning and roadway excavation. Cho et al. [14] conducted laboratory experiments on polymethyl methacrylate (PMMA) to study the empty-hole effect on fracture development. The study described the fracture development pattern was due to the stress condition. Based on the utilization of a high-speed camera and digital image, He and Yang [15] studied the fracture evolution between adjacent blasting holes. Chen et al. [16] used LS-DYNA to simulate the effect of superposition of stress and reflected waves on the effective stress distribution in the rock. Based on the same software and the investigation of the effect of empty hole, Meng et al. [17] monitored the change of tensile and compressive failure elements in the surrounding rock using Fortran language. They found that empty hole can significantly increase the stress concentration in the surrounding rock and the numbers of tensile fractures were considerably higher in the rock near the empty hole. In addition, Li et al. [18] used PMMA to conduct experiments to study the fracture development with empty hole. However, there are limited studies on the influence of in situ stress (static loading) on the mechanism of blasting-induced damage, especially the combined effect of in situ stress and empty hole on the fracture development.In the real situation, natural rock mass exists many preexisting discontinuities and homogeneity. To facilitate to study the effect of in situ stress on rock mass crack growth, the rock mass is simplified as an isotropic homogeneous continuum in this paper. A static-dynamic coupling numerical model was constructed using LS-DYNA; the calibration process was carried out against existing laboratory data. Based on the simulation results, the influences of in situ stress on rock damage and fracture development were discussed. Thereby, the combined effect of in situ stress and empty hole on the damage of rock and fracture development was analyzed.


## 2. Simulation Software Calculation Scheme

### 2.1. Lagrange, Euler, and ALE

LS-DYNA is a nonlinear software which simulates complex nonlinear dynamic problems. It has great advantages in analyzing geometric nonlinearity, material, and contact nonlinearity [19]. There are main calculation schemes: Lagrange (LF), Euler (EF), and Arbitrary Lagrangian-Eulerian (ALE) calculation schemes. EF and ALE calculation schemes are generally applied to solve nonlinear dynamic problem with large deformation [20]. Considering the effect of in situ stress, it is important to set the in situ stress in the model before performing the dynamic analysis. The implicit-explicit calculation method in LS-DYNA can solve the static-dynamic coupling problem.


### 2.2. Implicit-Explicit Analysis Method

The implicit-explicit calculation method in LS-DYNA can solve the problem of static loading prior to dynamic analysis. In numerical simulation, the explicit calculation method is more suitable for dynamic problem [21]. However, ANSYS/LS-DYNA is advantageous for implicit calculation. Hence, the implicit-explicit calculation method in LS-DYNA was chosen to simulate the in situ stress (static)–blasting (dynamic) load coupling effect. The procedures are listed as follows:
(a)
Conducting implicit analysis in ANSYS and set in situ stress in the model in advance. In implicit analysis, the chosen element type was entity 185. By then, the boundary conditions were set up on the prestress model. Subsequently, the simulation is processed(b)
Converting implicit analysis in ANSYS to explicit analysis in LS-DYNA. To generate in situ stress in LS-DYNA, entity 185 was converted as entity 164 with eight nodes in LS-DYNA. During this process, a dynamic relaxation file was generated and the displacement in the file was used to generate the in situ stress in LS-DYNA(c)
Conducting explicit analysis. Once the conversion was completed, the displacement and constitutive models in the file were used to generate in situ stress. Finally, the blasting load was simulated and analyzed


## 2.1. Lagrange, Euler, and ALE

LS-DYNA is a nonlinear software which simulates complex nonlinear dynamic problems. It has great advantages in analyzing geometric nonlinearity, material, and contact nonlinearity [19]. There are main calculation schemes: Lagrange (LF), Euler (EF), and Arbitrary Lagrangian-Eulerian (ALE) calculation schemes. EF and ALE calculation schemes are generally applied to solve nonlinear dynamic problem with large deformation [20]. Considering the effect of in situ stress, it is important to set the in situ stress in the model before performing the dynamic analysis. The implicit-explicit calculation method in LS-DYNA can solve the static-dynamic coupling problem.


## 2.2. Implicit-Explicit Analysis Method

The implicit-explicit calculation method in LS-DYNA can solve the problem of static loading prior to dynamic analysis. In numerical simulation, the explicit calculation method is more suitable for dynamic problem [21]. However, ANSYS/LS-DYNA is advantageous for implicit calculation. Hence, the implicit-explicit calculation method in LS-DYNA was chosen to simulate the in situ stress (static)–blasting (dynamic) load coupling effect. The procedures are listed as follows:
(a)
Conducting implicit analysis in ANSYS and set in situ stress in the model in advance. In implicit analysis, the chosen element type was entity 185. By then, the boundary conditions were set up on the prestress model. Subsequently, the simulation is processed(b)
Converting implicit analysis in ANSYS to explicit analysis in LS-DYNA. To generate in situ stress in LS-DYNA, entity 185 was converted as entity 164 with eight nodes in LS-DYNA. During this process, a dynamic relaxation file was generated and the displacement in the file was used to generate the in situ stress in LS-DYNA(c)
Conducting explicit analysis. Once the conversion was completed, the displacement and constitutive models in the file were used to generate in situ stress. Finally, the blasting load was simulated and analyzed


## 3. Model Calibration

Based on the modified HJC model, Riedel et al. [22] proposed a tensile and compressive damage model, i.e., Riedel-Hiermaier-Thoma (RHT) model. Comparing with the existing models, this model considers the influences of confining pressure, strain rate, strain hardening, and damage softening on rock during blasting and dynamic loading. To describe the effect of compaction and hardening, Mie-Gruneisen’s EOS and p-α models were embedded into LS-DYNA [23].RHT is an appropriate constitutive model for dynamic mechanical response of rock from blasting in deep mining locations [24, 25]. Although this model has been widely applied in the damage evolution simulation of concrete, it is rarely used for rock. In addition, this model has not been used extensively in LS-DYNA. To embed the RHT model to LS-DYNA and simulate damage evolution, it is critical to determine the mechanical properties based on the existing laboratory testing [26]. Based on the comparison with blasting results, the calibration and examination of the model were carried out.


### 3.1. RHT Model

In the RHT model, pressure is represented by Mie-Greisen, which is a polynomial function of the Hugoniot curve and has relationship withp-α compaction [23]. Figure 1 is p-α EOS illustration. In the compaction model, when the applied pressure is lower than the pore crushing pressure, the model will be always elastic. Once the pore crushing pressure is surpassed, the stiffness of the material will be reduced due to pore collapse and subsequently leads to a decrease in bulk modulus. As the relationship between pressure and volume is nonlinear, plastic volumetric deformation would occur when the pressure is high than the pore crushing pressure. To simulate pore collapse in the model, the internal parameter α was used to represent porosity. This parameter is the fraction between matrix material and porous material α decreases with the increasing pressure, and the loading path cannot be reversed. When the pressure reaches pore collapse pressure and compaction pressure, a equals 0 and 1, respectively.Figure 1
p-α EOS illustration in the RHT model.There are three stress limit surfaces in RHT to explain the strength and strain rate, which include yield surface, failure surface, and residual surface [23]. Figure 2 depicts a typical loading path, where the arrow represents the elastic state prior to the yield surface. The material exhibits linear hardening behaviour under plastic deformation outside this surface. When the material enters the yielding stage, damage starts to happen until the value reaches 1. The plastic strain and hardening properties of the material are applied to form an effective yield surface, which is created from an interpolation between the initial yield surface and failure surface. When the residual strength is reached, the material is deemed as fully damaged, whereas the residual strength is determined by residual surface. Borrvall and Riedel [23] provided a detailed discussion on the RHT model.Figure 2
Stress limit surfaces and loading path illustration.


### 3.2. Parameter Determination and Calibration

#### 3.2.1. Yield Surface Equation

Yield surface,σf, describes the maximum stress that the material can withstand without any damage, as expressed in Equation (1) [23]:
(1)σf=fcσf∗p∗,Frε̇p,p∗R3θ,p∗,in which θ is the Lode angle, εṗ is the plastic strain rate, and p∗ is the normalized pressure under uniaxial compressive strength fc. R3 is introduced to represent the change in the shear and tensile failure plane. Fr is the growth rate of dynamic strain rate, as defined in Equations (2)–(5) [23]:
(2)Frε̇p,p∗=Frc,3p∗≥Frc,Frc−3p∗−FrcFrc+Frtft∗Frt−Frc,Frc>3p∗≥−Frtft∗,Frt,−Frtft∗>3p∗,(3)Frc,tε̇p=ε̇pε0ċ,tβc,t,ε̇p≤εpċ,t,γc,tε̇p3,ε̇p>εpċ,t,(4)logγc=6βc−0.492,(5)logγt=7βt−0.492.In the above equations,c and t are compression and tension, βc and βt are compressive and tensile strain indexes, and ft∗ is the normalized tensile strength.


#### 3.2.2. Failure Surface Equation

By submittingFe and Fc into Equation (1), σel can be calculated:
(6)σel=fcσf∗p∗Fe,Frε̇p,p∗R3θ,p∗Fep∗Fcp∗.When the yield surface is reached, the damage from plastic deformation starts to accumulate:(7)D=∑Δεpεpf,εpf=D1p∗−1−Dpt∗D2,where D is the damage, ranging from 0 (undamaged) to 1 (fully damaged), ∆εp is the accumulated plastic strain, εpf is the equivalent plastic strain at failure, pt∗ is the normalized failure cutoff pressure, and D1 and D2 are constants.


#### 3.2.3. Residual Surface Equation

Due to the accumulation of damage, yield surface starts to soften and eventually becomes residual face:(8)σr∗p∗=Afp∗nf,p∗>0,0,p∗≤0,where Af and nf are constants.


#### 3.2.4. Strain Rate

This paper uses regression to represent these relationships, as discussed by Zhang and Zhao [27]. As displayed in Figure 3, the relationships between tensile and compression strengths and strain rate were obtained. Equations (9) and (10) [27] represent the growth rate of strain rates under compression and tension:
(9)Frcε̇p=0.509εp0.26,ε̇p>10s−1,ε̇pε̇0c0.026,ε̇p≤10s−1,where Frcεṗ is the rate of strain rate under compression, εṗ is the strain rate, ε0ċ is the standard strain rate under compression, and ε0ċ=3.0×10−5s−1.
(10)Frtε̇p=2.4εp0̇.3,ε̇p>10−1s−1,ε̇pε̇0t0.007,ε̇p≤10−1s−1,in which Frtεṗ is the rate of strain rate under tension, εṗ is the strain rate, ε0t is the standard strain rate under tension, and ε0ṫ=3.0×10−6s−1.Figure 3
Curve fitting on growth rates of strain rate: (a) dynamic compression (data from Zhang and Zhao [27]); (b) dynamic tension (data from Zhang and Zhao [27]).
(a)(b)Rock strength is influenced by strain rate. Because of this, the strain rate used in the rock tests is different from the reference strain rate, and the rock strength must be converted to a value that corresponds to the reference strain rate. The strain rate in the compression tests is controlled at4.63×10−5s−1, and the equivalent strain rate is 2.59×10−5s−1. To better estimate rock strength under high in situ stresses, the empirical equation of Hoek and Brown is used to obtain the strength of rock under various confining stresses.
(11)σ1=σ3+167.124σ3167.1+10.5,where σ1 and σ3 are the maximum and minimum effective stresses at failure.The mechanical parameters of the rock under various confining pressures can be obtained from Equations (9), (10), and (11).


#### 3.2.5. Jones-Wilkins-Lee (JWL) EOS

The Jones-Wilkins-Lee model is a high-energy combustion model, which has been implemented widely in numerical simulation because of its successful prediction on the blasting-induced pressure range. This study used JWL EOS to investigate the blasting process, which the equation can be expressed as(12)Pb=A1−ωR1Ve−VR1+B1−ωR1Ve−VR2+ωEiV,where Pb is the pressure applied to the borehole wall by a unit explosion product, V is the relative volume of the explosion product, Ei is the initial internal energy density of the explosion product, and A, B, R1, R2, and ω are the material constants determined by blasting experiments. In the paper, the PETN explosive is selected. However, due to the limitation of the measurement instruments, the JWL EOS parameters of the explosive are referred to Banadaki [28], as listed in Table 1.Table 1
Properties of explosive and JWL EOS.
ρ0 (g·cm-3)D (m·s-1)A (GPa)B (GPa)R1R2ωEi (GPa)1320669058621.65.811.770.2827.38


### 3.3. Simulation and Comparison with Laboratory Results

To verify the RHT model and utilize it for blasting analysis, a calibration process was conducted against the laboratory results from Banadaki [28]. Banadaki [28] carried out experiments on two types of cylindrical granite samples. In this paper, the results of Barre granite samples were used for comparison purposes. The diameter and height of the sample are 144 mm and 150 mm, with a predrilled 6.45 mm diameter blasting hole in the middle, as displayed in Figure 4(a). A 1.2 mm diameter copper tube was installed in the blasting hole to prevent gas leaking into fractures.Figure 4
Illustration of laboratory and numerical setup: (a) laboratory setup; (b) numerical setup.
(a)(b)A synthetic rock specimen was constructed in LS-DYNA with the same dimensions. In this study, an algorithm was implemented in LS-DYNA to model the detonation of explosive, and the Lagrangian meshes were used to model the rock and copper tube as shown in Figure4(b). Based on the sensitivity analysis, the model parameters were adjusted. The RHT model parameters of Barre are presented in Table 1. Tables 2–4 show the explosive, gas, and copper tube parameters, respectively.Table 2
Parameters for the RHT model for rock.
ParameterValueParameterValueDensity2660 kg/m3Reference compressive strain rate E0C3.0e−5Shear modulus17 GPaReference compressive strain rate E0T3.0e−6Erosion plastic strain2.0Break compressive strain rate EC3.0e25Parameter for polynomial EOS B01.22Break tensile strain rate ET3.0e25Parameter for polynomial EOS B11.22Compressive strain rate dependence exponent0.025Parameter for polynomial EOS T143.87 GPaTensile strain rate dependence exponent0.045Parameter for polynomial EOS T20Pressure influence on plastic flow in tension0.001Hugoniot polynomial coefficient A143.87 GPaCompressive yield surface parameter GC∗0.85Hugoniot polynomial coefficient A249.40 GPaTensile yield surface parameter GT∗0.4Hugoniot polynomial coefficient A311.62 GPaShear modulus reduction factor0.25Failure surface parameter A2.50Damage parameter D10.025Failure surface parameter N0.85Damage parameter D21.0Compressive strength150 MPaMinimum damaged residual strain0.01Relative shear strength0.07Residual surface parameter AF2.5Relative tensile strength0.05Residual surface parameter AN0.85Lode angle dependence factor Q00.72Gruneisen gamma0Lode angle dependence factor B0.01Crush pressure133 MPaCompaction pressure6 GPaPorosity exponent3Initial porosity1.0Table 3
Parameters for air.
ρ0 (g·cm-3)C4C51.29e−030.040.04Table 4
Parameters for copper.
Density (kg/m3)Young’s modulus (GPa)Poisson’s ratioYield stress (MPa)Tangent modulus (MPa)βC (s-1)P89301170.3540010001.34e65.28Figure5 shows the experimental results of blasting-induced fracture development and numerical results from this study. Contours were used to represent the damage from 0 to 1, where blue indicates 0 (not damaged) and red indicates 1 (fully damaged). The other colours represent different degrees of damage in the rock.Figure 5
Comparison between laboratory and numerical results on fracture development due to blasting: (a) laboratory results; (b) numerical results.
(a)(b)Figure5 revealed the good alignment between experimental and numerical results. After the explosion, a fractured zone can be observed around the blasting hole and the radial fractures developed outwards from the hole. Thereby, there were some fractures that occurred in tangential direction around the rock boundary. To further examine the model setup, stress conditions at different locations from the blasting hole were monitored, which can be seen in Figure 6. In this paper, the damage threshold is set to 0.85, and the rationality of the threshold setting is determined by numerical simulation and experimental results. The expansion of rock mass fissures is by removing the unit body with a damage level greater than 0.85 to form a blank area (cracks). Based on the comparison with experimental results, it can be seen that there is an overall agreement between simulation and experiments, suggesting that the model was calibrated properly for rock damage.Figure 6
Experimental and numerical results of stress conditions at different locations from the blasting hole.


## 3.1. RHT Model

In the RHT model, pressure is represented by Mie-Greisen, which is a polynomial function of the Hugoniot curve and has relationship withp-α compaction [23]. Figure 1 is p-α EOS illustration. In the compaction model, when the applied pressure is lower than the pore crushing pressure, the model will be always elastic. Once the pore crushing pressure is surpassed, the stiffness of the material will be reduced due to pore collapse and subsequently leads to a decrease in bulk modulus. As the relationship between pressure and volume is nonlinear, plastic volumetric deformation would occur when the pressure is high than the pore crushing pressure. To simulate pore collapse in the model, the internal parameter α was used to represent porosity. This parameter is the fraction between matrix material and porous material α decreases with the increasing pressure, and the loading path cannot be reversed. When the pressure reaches pore collapse pressure and compaction pressure, a equals 0 and 1, respectively.Figure 1
p-α EOS illustration in the RHT model.There are three stress limit surfaces in RHT to explain the strength and strain rate, which include yield surface, failure surface, and residual surface [23]. Figure 2 depicts a typical loading path, where the arrow represents the elastic state prior to the yield surface. The material exhibits linear hardening behaviour under plastic deformation outside this surface. When the material enters the yielding stage, damage starts to happen until the value reaches 1. The plastic strain and hardening properties of the material are applied to form an effective yield surface, which is created from an interpolation between the initial yield surface and failure surface. When the residual strength is reached, the material is deemed as fully damaged, whereas the residual strength is determined by residual surface. Borrvall and Riedel [23] provided a detailed discussion on the RHT model.Figure 2
Stress limit surfaces and loading path illustration.


## 3.2. Parameter Determination and Calibration

### 3.2.1. Yield Surface Equation

Yield surface,σf, describes the maximum stress that the material can withstand without any damage, as expressed in Equation (1) [23]:
(1)σf=fcσf∗p∗,Frε̇p,p∗R3θ,p∗,in which θ is the Lode angle, εṗ is the plastic strain rate, and p∗ is the normalized pressure under uniaxial compressive strength fc. R3 is introduced to represent the change in the shear and tensile failure plane. Fr is the growth rate of dynamic strain rate, as defined in Equations (2)–(5) [23]:
(2)Frε̇p,p∗=Frc,3p∗≥Frc,Frc−3p∗−FrcFrc+Frtft∗Frt−Frc,Frc>3p∗≥−Frtft∗,Frt,−Frtft∗>3p∗,(3)Frc,tε̇p=ε̇pε0ċ,tβc,t,ε̇p≤εpċ,t,γc,tε̇p3,ε̇p>εpċ,t,(4)logγc=6βc−0.492,(5)logγt=7βt−0.492.In the above equations,c and t are compression and tension, βc and βt are compressive and tensile strain indexes, and ft∗ is the normalized tensile strength.


### 3.2.2. Failure Surface Equation

By submittingFe and Fc into Equation (1), σel can be calculated:
(6)σel=fcσf∗p∗Fe,Frε̇p,p∗R3θ,p∗Fep∗Fcp∗.When the yield surface is reached, the damage from plastic deformation starts to accumulate:(7)D=∑Δεpεpf,εpf=D1p∗−1−Dpt∗D2,where D is the damage, ranging from 0 (undamaged) to 1 (fully damaged), ∆εp is the accumulated plastic strain, εpf is the equivalent plastic strain at failure, pt∗ is the normalized failure cutoff pressure, and D1 and D2 are constants.


### 3.2.3. Residual Surface Equation

Due to the accumulation of damage, yield surface starts to soften and eventually becomes residual face:(8)σr∗p∗=Afp∗nf,p∗>0,0,p∗≤0,where Af and nf are constants.


### 3.2.4. Strain Rate

This paper uses regression to represent these relationships, as discussed by Zhang and Zhao [27]. As displayed in Figure 3, the relationships between tensile and compression strengths and strain rate were obtained. Equations (9) and (10) [27] represent the growth rate of strain rates under compression and tension:
(9)Frcε̇p=0.509εp0.26,ε̇p>10s−1,ε̇pε̇0c0.026,ε̇p≤10s−1,where Frcεṗ is the rate of strain rate under compression, εṗ is the strain rate, ε0ċ is the standard strain rate under compression, and ε0ċ=3.0×10−5s−1.
(10)Frtε̇p=2.4εp0̇.3,ε̇p>10−1s−1,ε̇pε̇0t0.007,ε̇p≤10−1s−1,in which Frtεṗ is the rate of strain rate under tension, εṗ is the strain rate, ε0t is the standard strain rate under tension, and ε0ṫ=3.0×10−6s−1.Figure 3
Curve fitting on growth rates of strain rate: (a) dynamic compression (data from Zhang and Zhao [27]); (b) dynamic tension (data from Zhang and Zhao [27]).
(a)(b)Rock strength is influenced by strain rate. Because of this, the strain rate used in the rock tests is different from the reference strain rate, and the rock strength must be converted to a value that corresponds to the reference strain rate. The strain rate in the compression tests is controlled at4.63×10−5s−1, and the equivalent strain rate is 2.59×10−5s−1. To better estimate rock strength under high in situ stresses, the empirical equation of Hoek and Brown is used to obtain the strength of rock under various confining stresses.
(11)σ1=σ3+167.124σ3167.1+10.5,where σ1 and σ3 are the maximum and minimum effective stresses at failure.The mechanical parameters of the rock under various confining pressures can be obtained from Equations (9), (10), and (11).


### 3.2.5. Jones-Wilkins-Lee (JWL) EOS

The Jones-Wilkins-Lee model is a high-energy combustion model, which has been implemented widely in numerical simulation because of its successful prediction on the blasting-induced pressure range. This study used JWL EOS to investigate the blasting process, which the equation can be expressed as(12)Pb=A1−ωR1Ve−VR1+B1−ωR1Ve−VR2+ωEiV,where Pb is the pressure applied to the borehole wall by a unit explosion product, V is the relative volume of the explosion product, Ei is the initial internal energy density of the explosion product, and A, B, R1, R2, and ω are the material constants determined by blasting experiments. In the paper, the PETN explosive is selected. However, due to the limitation of the measurement instruments, the JWL EOS parameters of the explosive are referred to Banadaki [28], as listed in Table 1.Table 1
Properties of explosive and JWL EOS.
ρ0 (g·cm-3)D (m·s-1)A (GPa)B (GPa)R1R2ωEi (GPa)1320669058621.65.811.770.2827.38


## 3.2.1. Yield Surface Equation

Yield surface,σf, describes the maximum stress that the material can withstand without any damage, as expressed in Equation (1) [23]:
(1)σf=fcσf∗p∗,Frε̇p,p∗R3θ,p∗,in which θ is the Lode angle, εṗ is the plastic strain rate, and p∗ is the normalized pressure under uniaxial compressive strength fc. R3 is introduced to represent the change in the shear and tensile failure plane. Fr is the growth rate of dynamic strain rate, as defined in Equations (2)–(5) [23]:
(2)Frε̇p,p∗=Frc,3p∗≥Frc,Frc−3p∗−FrcFrc+Frtft∗Frt−Frc,Frc>3p∗≥−Frtft∗,Frt,−Frtft∗>3p∗,(3)Frc,tε̇p=ε̇pε0ċ,tβc,t,ε̇p≤εpċ,t,γc,tε̇p3,ε̇p>εpċ,t,(4)logγc=6βc−0.492,(5)logγt=7βt−0.492.In the above equations,c and t are compression and tension, βc and βt are compressive and tensile strain indexes, and ft∗ is the normalized tensile strength.


## 3.2.2. Failure Surface Equation

By submittingFe and Fc into Equation (1), σel can be calculated:
(6)σel=fcσf∗p∗Fe,Frε̇p,p∗R3θ,p∗Fep∗Fcp∗.When the yield surface is reached, the damage from plastic deformation starts to accumulate:(7)D=∑Δεpεpf,εpf=D1p∗−1−Dpt∗D2,where D is the damage, ranging from 0 (undamaged) to 1 (fully damaged), ∆εp is the accumulated plastic strain, εpf is the equivalent plastic strain at failure, pt∗ is the normalized failure cutoff pressure, and D1 and D2 are constants.


## 3.2.3. Residual Surface Equation

Due to the accumulation of damage, yield surface starts to soften and eventually becomes residual face:(8)σr∗p∗=Afp∗nf,p∗>0,0,p∗≤0,where Af and nf are constants.


## 3.2.4. Strain Rate

This paper uses regression to represent these relationships, as discussed by Zhang and Zhao [27]. As displayed in Figure 3, the relationships between tensile and compression strengths and strain rate were obtained. Equations (9) and (10) [27] represent the growth rate of strain rates under compression and tension:
(9)Frcε̇p=0.509εp0.26,ε̇p>10s−1,ε̇pε̇0c0.026,ε̇p≤10s−1,where Frcεṗ is the rate of strain rate under compression, εṗ is the strain rate, ε0ċ is the standard strain rate under compression, and ε0ċ=3.0×10−5s−1.
(10)Frtε̇p=2.4εp0̇.3,ε̇p>10−1s−1,ε̇pε̇0t0.007,ε̇p≤10−1s−1,in which Frtεṗ is the rate of strain rate under tension, εṗ is the strain rate, ε0t is the standard strain rate under tension, and ε0ṫ=3.0×10−6s−1.Figure 3
Curve fitting on growth rates of strain rate: (a) dynamic compression (data from Zhang and Zhao [27]); (b) dynamic tension (data from Zhang and Zhao [27]).
(a)(b)Rock strength is influenced by strain rate. Because of this, the strain rate used in the rock tests is different from the reference strain rate, and the rock strength must be converted to a value that corresponds to the reference strain rate. The strain rate in the compression tests is controlled at4.63×10−5s−1, and the equivalent strain rate is 2.59×10−5s−1. To better estimate rock strength under high in situ stresses, the empirical equation of Hoek and Brown is used to obtain the strength of rock under various confining stresses.
(11)σ1=σ3+167.124σ3167.1+10.5,where σ1 and σ3 are the maximum and minimum effective stresses at failure.The mechanical parameters of the rock under various confining pressures can be obtained from Equations (9), (10), and (11).


## 3.2.5. Jones-Wilkins-Lee (JWL) EOS

The Jones-Wilkins-Lee model is a high-energy combustion model, which has been implemented widely in numerical simulation because of its successful prediction on the blasting-induced pressure range. This study used JWL EOS to investigate the blasting process, which the equation can be expressed as(12)Pb=A1−ωR1Ve−VR1+B1−ωR1Ve−VR2+ωEiV,where Pb is the pressure applied to the borehole wall by a unit explosion product, V is the relative volume of the explosion product, Ei is the initial internal energy density of the explosion product, and A, B, R1, R2, and ω are the material constants determined by blasting experiments. In the paper, the PETN explosive is selected. However, due to the limitation of the measurement instruments, the JWL EOS parameters of the explosive are referred to Banadaki [28], as listed in Table 1.Table 1
Properties of explosive and JWL EOS.
ρ0 (g·cm-3)D (m·s-1)A (GPa)B (GPa)R1R2ωEi (GPa)1320669058621.65.811.770.2827.38


## 3.3. Simulation and Comparison with Laboratory Results

To verify the RHT model and utilize it for blasting analysis, a calibration process was conducted against the laboratory results from Banadaki [28]. Banadaki [28] carried out experiments on two types of cylindrical granite samples. In this paper, the results of Barre granite samples were used for comparison purposes. The diameter and height of the sample are 144 mm and 150 mm, with a predrilled 6.45 mm diameter blasting hole in the middle, as displayed in Figure 4(a). A 1.2 mm diameter copper tube was installed in the blasting hole to prevent gas leaking into fractures.Figure 4
Illustration of laboratory and numerical setup: (a) laboratory setup; (b) numerical setup.
(a)(b)A synthetic rock specimen was constructed in LS-DYNA with the same dimensions. In this study, an algorithm was implemented in LS-DYNA to model the detonation of explosive, and the Lagrangian meshes were used to model the rock and copper tube as shown in Figure4(b). Based on the sensitivity analysis, the model parameters were adjusted. The RHT model parameters of Barre are presented in Table 1. Tables 2–4 show the explosive, gas, and copper tube parameters, respectively.Table 2
Parameters for the RHT model for rock.
ParameterValueParameterValueDensity2660 kg/m3Reference compressive strain rate E0C3.0e−5Shear modulus17 GPaReference compressive strain rate E0T3.0e−6Erosion plastic strain2.0Break compressive strain rate EC3.0e25Parameter for polynomial EOS B01.22Break tensile strain rate ET3.0e25Parameter for polynomial EOS B11.22Compressive strain rate dependence exponent0.025Parameter for polynomial EOS T143.87 GPaTensile strain rate dependence exponent0.045Parameter for polynomial EOS T20Pressure influence on plastic flow in tension0.001Hugoniot polynomial coefficient A143.87 GPaCompressive yield surface parameter GC∗0.85Hugoniot polynomial coefficient A249.40 GPaTensile yield surface parameter GT∗0.4Hugoniot polynomial coefficient A311.62 GPaShear modulus reduction factor0.25Failure surface parameter A2.50Damage parameter D10.025Failure surface parameter N0.85Damage parameter D21.0Compressive strength150 MPaMinimum damaged residual strain0.01Relative shear strength0.07Residual surface parameter AF2.5Relative tensile strength0.05Residual surface parameter AN0.85Lode angle dependence factor Q00.72Gruneisen gamma0Lode angle dependence factor B0.01Crush pressure133 MPaCompaction pressure6 GPaPorosity exponent3Initial porosity1.0Table 3
Parameters for air.
ρ0 (g·cm-3)C4C51.29e−030.040.04Table 4
Parameters for copper.
Density (kg/m3)Young’s modulus (GPa)Poisson’s ratioYield stress (MPa)Tangent modulus (MPa)βC (s-1)P89301170.3540010001.34e65.28Figure5 shows the experimental results of blasting-induced fracture development and numerical results from this study. Contours were used to represent the damage from 0 to 1, where blue indicates 0 (not damaged) and red indicates 1 (fully damaged). The other colours represent different degrees of damage in the rock.Figure 5
Comparison between laboratory and numerical results on fracture development due to blasting: (a) laboratory results; (b) numerical results.
(a)(b)Figure5 revealed the good alignment between experimental and numerical results. After the explosion, a fractured zone can be observed around the blasting hole and the radial fractures developed outwards from the hole. Thereby, there were some fractures that occurred in tangential direction around the rock boundary. To further examine the model setup, stress conditions at different locations from the blasting hole were monitored, which can be seen in Figure 6. In this paper, the damage threshold is set to 0.85, and the rationality of the threshold setting is determined by numerical simulation and experimental results. The expansion of rock mass fissures is by removing the unit body with a damage level greater than 0.85 to form a blank area (cracks). Based on the comparison with experimental results, it can be seen that there is an overall agreement between simulation and experiments, suggesting that the model was calibrated properly for rock damage.Figure 6
Experimental and numerical results of stress conditions at different locations from the blasting hole.


## 4. Simulation Results and Discussion

### 4.1. Model Construction

In the real situation, natural rock mass exists many preexisting discontinuities and homogeneity. To facilitate to study the effect of in situ stress on rock mass crack growth, the rock mass is simplified as an isotropic homogeneous continuum in this simulation. The magnitudes of in situ stress generally increase with the depth of cover, and it is difficult to excavate tunnel or extract minerals under high in situ stress conditions [29–31]. In this paper, a plane model which considers the static-dynamic coupling effect was implemented, which can represent the blasting-induced rock damage and fracture development in deep locations [32]. A numerical model with dimensions of 600cminlength×600cm (in width) is established by the numerical modelling code LS-DYNA (Figure 7). The borehole diameter is 5.0 cm, the roll diameter is 4.0 cm, and the empty borehole diameter is 5.0 cm. To eliminate the effects of explosion and rock boundary on simulation, the model boundary is set to a nonreflective boundary. Model grids are based on a Cartesian coordinate system using “mapping” uniform division, and the average size of explosives is 1.8 mm. The number of explosive units is 697. The minimum and maximum sizes of air and rock are 2.4 mm and 6.7 mm, respectively. The spacing ratio of adjacent elements is 1.15. The number of air and rock units is 2153276 and 196723, respectively. To accurately reflect the situation of explosion and rock mass damage, the model adopts the ALE algorithm, explosives and rocks adopt the Lagrange algorithm, and air adopts the Euler algorithm.Figure 7
Blasting model without empty hole.There were four different stress conditions used in this study. By changing the magnitude of vertical stress, the effect of in situ stress on rock damage and fracture development was investigated. The horizontal stress was set as 0 and 30 MPa, whereas the vertical stress was set at 0, 30, and 70 MPa. Each case is shown as follows:Case 1.
σh=σv=0 MPa. The blasting model does not have any in situ stress.Case 2.
σh=30 MPa, σv=0 MPa. To investigate low-stress conditions on rock damage.Case 3.
σh=30 MPa, σv=30 MPa. To compare with Cases 1 and 2 and further explore the influence of in situ stress.Case 4.
σh=30 MPa, σv=70 MPa. To study the complex in situ stress conditions.Figures7 and 8 show the simulation without and with empty holes. The simulation procedures include the following: (a) applying σv and σh to the boundary; (b) using control_dynamic_relaxation in LS-DYNA to apply initial in situ stress on the rock boundary for stress initialisation; (c) after the initialisation, the blasting simulation was carried out to study the rock damage and fracture development.Figure 8
Blasting model with empty hole.


### 4.2. Effect of In Situ Stress on the Degree of Damage of Rock

Figure9 shows the blasting-induced rock damage in four scenarios. Based on the figure, it can be seen that the rock damage was more severe without in situ stress, where the fracture development was more obvious. The fracture development general decreases with the increasing in situ stress; this is due to the increasing confinement that prevents further fracture development. Based on the comparison between Figures 9(a) and 9(b), rock damage and fracture development decreased at θ=0° and 180° and increased along 90°and 270°. Based on the comparison between Figures 9(b) and 9(d), when the σv was greater than σh, rock damage and fracture development increased at θ=0° and 180° and decreased along 90° and 270°. By comparing Figures 9(b) and 9(c), as σv increases, the fracture development gradually diminishes between 60° and 120°, while it increased at 0° and 180°. When there is difference between σv and σh, the anisotropic stress field led to the fracture developing towards a high-stress concentration zone, whereas the fracture development in a low-stress concentration zone was reduced.Figure 9
Fracture development under different stress conditions. (a)σh=σv=0 MPa; (b) σh=0 MPa, σv=30 MPa; (c) σh=30 MPa, σv=30 MPa; (d) σh=70 MPa, σv=30 MPa.
(a)(b)(c)(d)Table5 summarises the total fracture areas (blank areas) under four scenarios. The blank areas can be indicative to rock damage extent. From Table 5, Case 1 has the highest blank area whereas Case 4 has the lowest blank area. This suggests that high in situ stress prevented the fracture development.Table 5
Area of blanked elements (unit: m2).
Case 1Case 2Case 3Case 42.041.871.591.06


### 4.3. Effect of Empty Hole on the Degree of Damage of Rock

Figure10 shows the simulation results at σh=σv=0 MPa. When the shock wave is attenuated into stress wave, the tensile stress is higher than the tensile strength, such that tensile cracks occurred in the rock. When the stress wave transmits to the empty borehole free surface, the compressive wave reflects at the free surface and forms the tensile wave. The compressive wave superimposes on the tensile wave nearby the empty borehole, leading to a notable increase of the tensile stress component. If the superimposed tensile stress exceeds the tensile strength of the rock mass, the superimposed tensile wave induces tensile damage and generates the fractures propagating toward the rock mass. Finally, the fractures develop at 90 degrees.Figure 10
Fracture development atσh=σv=0 MPa.
(a)
Without empty hole(b)
With empty holeBy comparing the crack propagation process of the secondary crack and the primary crack in Figures10(a) and 10(b), it can be found that under the influence of empty hole, the directions of the secondary crack gradually changed towards the side of empty hole during the propagation. This indicates that empty hole can change the propagation direction of fracture and achieve controlled fracture development.


### 4.4. The Combined Effect of In Situ Stress and Empty Hole on the Degree of Damage of Rock

To study the combined effect of in situ stress and empty hole on the evolution of rock blasting damage, Figure11 shows the distribution of rock damage at 500 μs under four conditions. Figure 11 shows the fracture development was influenced by in situ stresses and subsequently resulted in different fracture patterns. When σv<30 MPa, the rock damage in horizontal direction was more severe than that in vertical direction, vice versa for σv>30 MPa. In addition, if the differential stress between the two principal stresses was high, the fracture lengths in vertical and horizontal directions were considerably different. Based on the comparison between Figures 9 and 11, the empty hole resulted in fracture prevention at 60° and 120° under the same stress condition. Thereby, the secondary crack propagation direction gradually deviated to 90°, while the primary crack developed significantly in 90° direction. This suggests that empty hole can facilitate fracture development at 90° during blasting process. It is also possible to control the fracture direction as the secondary fracture development direction was deviated towards the empty hole (at 90°).Figure 11
Fracture development under various in situ stress and empty holes: (a)σh=σv=0 MPa; (b) σh=30 MPa, σv=0 MPa; (c) σh=30 MPa, σv=30 MPa; (d) σh=30 MPa, σv=70 MPa.
(a)(b)(c)(d)


## 4.1. Model Construction

In the real situation, natural rock mass exists many preexisting discontinuities and homogeneity. To facilitate to study the effect of in situ stress on rock mass crack growth, the rock mass is simplified as an isotropic homogeneous continuum in this simulation. The magnitudes of in situ stress generally increase with the depth of cover, and it is difficult to excavate tunnel or extract minerals under high in situ stress conditions [29–31]. In this paper, a plane model which considers the static-dynamic coupling effect was implemented, which can represent the blasting-induced rock damage and fracture development in deep locations [32]. A numerical model with dimensions of 600cminlength×600cm (in width) is established by the numerical modelling code LS-DYNA (Figure 7). The borehole diameter is 5.0 cm, the roll diameter is 4.0 cm, and the empty borehole diameter is 5.0 cm. To eliminate the effects of explosion and rock boundary on simulation, the model boundary is set to a nonreflective boundary. Model grids are based on a Cartesian coordinate system using “mapping” uniform division, and the average size of explosives is 1.8 mm. The number of explosive units is 697. The minimum and maximum sizes of air and rock are 2.4 mm and 6.7 mm, respectively. The spacing ratio of adjacent elements is 1.15. The number of air and rock units is 2153276 and 196723, respectively. To accurately reflect the situation of explosion and rock mass damage, the model adopts the ALE algorithm, explosives and rocks adopt the Lagrange algorithm, and air adopts the Euler algorithm.Figure 7
Blasting model without empty hole.There were four different stress conditions used in this study. By changing the magnitude of vertical stress, the effect of in situ stress on rock damage and fracture development was investigated. The horizontal stress was set as 0 and 30 MPa, whereas the vertical stress was set at 0, 30, and 70 MPa. Each case is shown as follows:Case 1.
σh=σv=0 MPa. The blasting model does not have any in situ stress.Case 2.
σh=30 MPa, σv=0 MPa. To investigate low-stress conditions on rock damage.Case 3.
σh=30 MPa, σv=30 MPa. To compare with Cases 1 and 2 and further explore the influence of in situ stress.Case 4.
σh=30 MPa, σv=70 MPa. To study the complex in situ stress conditions.Figures7 and 8 show the simulation without and with empty holes. The simulation procedures include the following: (a) applying σv and σh to the boundary; (b) using control_dynamic_relaxation in LS-DYNA to apply initial in situ stress on the rock boundary for stress initialisation; (c) after the initialisation, the blasting simulation was carried out to study the rock damage and fracture development.Figure 8
Blasting model with empty hole.


## 4.2. Effect of In Situ Stress on the Degree of Damage of Rock

Figure9 shows the blasting-induced rock damage in four scenarios. Based on the figure, it can be seen that the rock damage was more severe without in situ stress, where the fracture development was more obvious. The fracture development general decreases with the increasing in situ stress; this is due to the increasing confinement that prevents further fracture development. Based on the comparison between Figures 9(a) and 9(b), rock damage and fracture development decreased at θ=0° and 180° and increased along 90°and 270°. Based on the comparison between Figures 9(b) and 9(d), when the σv was greater than σh, rock damage and fracture development increased at θ=0° and 180° and decreased along 90° and 270°. By comparing Figures 9(b) and 9(c), as σv increases, the fracture development gradually diminishes between 60° and 120°, while it increased at 0° and 180°. When there is difference between σv and σh, the anisotropic stress field led to the fracture developing towards a high-stress concentration zone, whereas the fracture development in a low-stress concentration zone was reduced.Figure 9
Fracture development under different stress conditions. (a)σh=σv=0 MPa; (b) σh=0 MPa, σv=30 MPa; (c) σh=30 MPa, σv=30 MPa; (d) σh=70 MPa, σv=30 MPa.
(a)(b)(c)(d)Table5 summarises the total fracture areas (blank areas) under four scenarios. The blank areas can be indicative to rock damage extent. From Table 5, Case 1 has the highest blank area whereas Case 4 has the lowest blank area. This suggests that high in situ stress prevented the fracture development.Table 5
Area of blanked elements (unit: m2).
Case 1Case 2Case 3Case 42.041.871.591.06


## 4.3. Effect of Empty Hole on the Degree of Damage of Rock

Figure10 shows the simulation results at σh=σv=0 MPa. When the shock wave is attenuated into stress wave, the tensile stress is higher than the tensile strength, such that tensile cracks occurred in the rock. When the stress wave transmits to the empty borehole free surface, the compressive wave reflects at the free surface and forms the tensile wave. The compressive wave superimposes on the tensile wave nearby the empty borehole, leading to a notable increase of the tensile stress component. If the superimposed tensile stress exceeds the tensile strength of the rock mass, the superimposed tensile wave induces tensile damage and generates the fractures propagating toward the rock mass. Finally, the fractures develop at 90 degrees.Figure 10
Fracture development atσh=σv=0 MPa.
(a)
Without empty hole(b)
With empty holeBy comparing the crack propagation process of the secondary crack and the primary crack in Figures10(a) and 10(b), it can be found that under the influence of empty hole, the directions of the secondary crack gradually changed towards the side of empty hole during the propagation. This indicates that empty hole can change the propagation direction of fracture and achieve controlled fracture development.


## 4.4. The Combined Effect of In Situ Stress and Empty Hole on the Degree of Damage of Rock

To study the combined effect of in situ stress and empty hole on the evolution of rock blasting damage, Figure11 shows the distribution of rock damage at 500 μs under four conditions. Figure 11 shows the fracture development was influenced by in situ stresses and subsequently resulted in different fracture patterns. When σv<30 MPa, the rock damage in horizontal direction was more severe than that in vertical direction, vice versa for σv>30 MPa. In addition, if the differential stress between the two principal stresses was high, the fracture lengths in vertical and horizontal directions were considerably different. Based on the comparison between Figures 9 and 11, the empty hole resulted in fracture prevention at 60° and 120° under the same stress condition. Thereby, the secondary crack propagation direction gradually deviated to 90°, while the primary crack developed significantly in 90° direction. This suggests that empty hole can facilitate fracture development at 90° during blasting process. It is also possible to control the fracture direction as the secondary fracture development direction was deviated towards the empty hole (at 90°).Figure 11
Fracture development under various in situ stress and empty holes: (a)σh=σv=0 MPa; (b) σh=30 MPa, σv=0 MPa; (c) σh=30 MPa, σv=30 MPa; (d) σh=30 MPa, σv=70 MPa.
(a)(b)(c)(d)


## 5. Conclusion

In situ stress is high and complex in deep locations. This paper used numerical simulation technique to compare the existing laboratory data with simulation results as well as verified the numerical model. By then, the influences of in situ stress and empty hole on rock damage and fracture development were investigated. The following conclusions were drawn based on the analysis:(1)
The magnitude and state of the in situ stress in the rock can influence the damage evolution and the crack initiation and development substantially, which in turn affect the rock breakage. As the in situ stress continued to increase, the degree of fracture development gradually decreased. Rock mass failure is caused by a combination of shock waves, stress waves, stress wave superposition, and reflected stress waves. The length of damage cracks gradually decreases with increasing hydrostatic pressure. With the increasing maximum principal stress, the main damage cracks of the rock mass gradually deviate to the direction of the maximum principal stress during production blasting, and the length of the maximum principal crack gradually decreases(2)
Asσv increased, the empty hole resulted in the prevention of fracture propagation at 60° and 120°, and fracture developed significantly along 0° and 180°. Given the differential stress between the two principal stresses was high, the fracture lengths in vertical and horizontal directions were considerably different. The rock is affected by the anisotropic stress field, which caused the cracks to propagate mainly in the direction of high-stress concentration, and the crack propagation of rock in the direction of low in situ stress was suppressed(3)
Under the same stress condition, the empty hole resulted in the prevention of fracture development at 60° and 120°. Thereby, the secondary crack propagation direction gradually deviated to 90°, while the primary crack developed significantly in 90° direction. This suggests that empty hole can facilitate fracture development at 90° and prevent secondary fracture development along 60° and 120°. Therefore, it is also possible to control the fracture direction as the secondary fracture development direction was deviated towards the empty holeWith the increase of vertical in situ stress, cracks mainly propagate along the vertical direction, which has a negative effect on roof cracking. Therefore, the spacing of the holes is reduced accordingly. At the same time, empty holes are arranged in the middle of the hole to suppress the crack propagation in the vertical direction.



---

*Source: 6643042-2021-05-08.xml*