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.

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

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 with p-α 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, equals 0 and 1, respectively.

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.

3.2. Parameter Determination and Calibration
3.2.1. Yield Surface Equation

Yield surface, , describes the maximum stress that the material can withstand without any damage, as expressed in Equation (1) [23]: in which is the Lode angle, is the plastic strain rate, and is the normalized pressure under uniaxial compressive strength . is introduced to represent the change in the shear and tensile failure plane. is the growth rate of dynamic strain rate, as defined in Equations (2)–(5) [23]:

In the above equations, and are compression and tension, and are compressive and tensile strain indexes, and is the normalized tensile strength.

3.2.2. Failure Surface Equation

By submitting and into Equation (1), can be calculated:

When the yield surface is reached, the damage from plastic deformation starts to accumulate: where is the damage, ranging from 0 (undamaged) to 1 (fully damaged), is the accumulated plastic strain, is the equivalent plastic strain at failure, is the normalized failure cutoff pressure, and and are constants.

3.2.3. Residual Surface Equation

Due to the accumulation of damage, yield surface starts to soften and eventually becomes residual face: where and 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: where is the rate of strain rate under compression, is the strain rate, is the standard strain rate under compression, and . in which is the rate of strain rate under tension, is the strain rate, is the standard strain rate under tension, and .

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 at , and the equivalent strain rate is . 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. where and 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 where is the pressure applied to the borehole wall by a unit explosion product, is the relative volume of the explosion product, is the initial internal energy density of the explosion product, and , , , , 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.

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.

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 Figure 4(b). Based on the sensitivity analysis, the model parameters were adjusted. The RHT model parameters of Barre are presented in Table 1. Tables 24 show the explosive, gas, and copper tube parameters, respectively.

Figure 5 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 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.

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 [2931]. 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 (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.

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.  MPa. The blasting model does not have any in situ stress.

Case 2.  MPa,  MPa. To investigate low-stress conditions on rock damage.

Case 3.  MPa,  MPa. To compare with Cases 1 and 2 and further explore the influence of in situ stress.

Case 4.  MPa,  MPa. To study the complex in situ stress conditions.

Figures 7 and 8 show the simulation without and with empty holes. The simulation procedures include the following: (a) applying and 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.

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

Figure 9 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 and 180° and increased along 90°and 270°. Based on the comparison between Figures 9(b) and 9(d), when the was greater than , rock damage and fracture development increased at and 180° and decreased along 90° and 270°. By comparing Figures 9(b) and 9(c), as increases, the fracture development gradually diminishes between 60° and 120°, while it increased at 0° and 180°. When there is difference between and , 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.

Table 5 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.

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

Figure 10 shows the simulation results at  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.

By comparing the crack propagation process of the secondary crack and the primary crack in Figures 10(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, Figure 11 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  MPa, the rock damage in horizontal direction was more severe than that in vertical direction, vice versa for  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°).

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 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 hole

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

Data Availability

Get the paper data from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Authors’ Contributions

M.N. conceived the research and wrote the paper, J.B. and W.X. directed the paper structure. Y.C. and W.B. analyzed the data and contributed numerical simulation. W.W. and S.W. participated in collecting the data and language editing. All authors have read and approved the final manuscript.


This research has been supported by the National Natural Science Foundation of China (52074239 and 51927807), the Fundamental Research Funds for the Central Universities (2018ZDPY02 and 2019XKQYMS63), and the independent research project of State Key Laboratory of Coal Resources and Safe Mining, CUMT (SKLCRSM19X0015).