Abstract

Based on the particle flow code (PFC2D) within the discrete element method (DEM), the rock mass model was established according to the site rock conditions and the rock mass blasting was simulated by the explosive particle expansion method. The influence of various parameters (the peak pressure action coefficient of the borehole wall, explosive particle buried depth, and charge mode) in the explosive particle expansion method on blasting effect was investigated. Furthermore, the relationship between the various parameters and the geometry size of the blasting crater was obtained. By comparing the size of blasting crater in the field blasting test and numerical simulation example, the reliability of rock mass blasting simulated by the explosive particle expansion method using PFC is verified. The result shows that this paper provides a reliable new numerical simulation method for rock mass blasting and can be used to guide field blasting.

1. Introduction

With the development of numerical simulation and computer technology, the numerical simulation technology has been gradually applied in the field of rock blasting. In recent years, blasting researchers have studied the rock damage mechanism under blast loading by computer simulation, predicted the blasting funnel forming effect, and selected reasonable blasting parameters. The widely used numerical methods can be roughly divided into two categories: finite element method (FEM) based on continuum mechanics and discrete element method (DEM) based on Newton's classical mechanics [1].

Jing [2] reviewed the technology, progress, problems, and possible future development of numerical modelling of rock mechanics and introduced in detail the technical status and progress related to the main methods. As a tool of numerical modelling, discrete element software (PFC) has its own advantages and can simulate large deformation rock problems [3]. Park et al. [4] used the discrete element software (PFC3D) to simulate rock joints and used the code to carry out a series of direct shear tests. It was proved that the contact bond model is feasible to reproduce rock joints, and the effects of geometric characteristics and microscopic properties of joints on their shear behavior are studied. Based on the agglomerated particle model, Rong et al. [5] established four rock samples with different particle shapes and studied the influence of particle shape on rock mechanics.

Deng et al. [68] studied the dynamic response of underground and ground structures under the action of underground explosion by using the verified finite element code AUTODYN and DEM-based UDEC program. The joints in rock mass have great influence on the propagation of blast stress wave. Hao et al. [9, 10] studied the influence of rock joints with different dip angles on stress wave propagation. Fakhimi et al. [11] used a smooth particle hydrodynamics method (SPH) to simulate the mechanical behavior of rock using a bonded particle model. It assumed the gas particles can interact with rock particles. The mixed model has a good simulation of the process of rock crack growth under gas-rock interaction during rock blasting. Singh [12] studied the influence of explosion design parameters on rock fragmentation. Kunihisa et al. [13] analyzed the internal stress changes in materials caused by a free surface blasting and simulated the cracks and crater formation process. José Sanchidrián et al. [14] recorded the initial velocity of the seismic wave and blasting rock surface using a seismograph and high-speed camera and determined the seismic wave energy, kinetic energy, and fracture energy during the blasting, respectively. Meysam Lak et al. [15] deduced the general Green’s function solution of elastic wave propagation caused by rock blasting, used Navier’s equation as the control equation, and obtained a general two-dimensional elastic Green function in displacement. The strain and stress fields related to the displacement Green function are also obtained by using elasticity theory. Then, the general two-dimensional problem including explosion in a blast hole is solved by using the analytical solution obtained. Zhu et al. [16] used AUTODYN 2D code to carry out numerical blasting of a circular rock model with a decomposition explosion center. The influence of various factors on rock dynamic fracture is studied. An et al. [17] used the hybrid finite-discrete element method (FEM-DEM) to simulate rock fracture and debris accumulation. The process of rock fragmentation during blasting is reproduced. Liu et al. [18] analyzed the rock breaking mechanism and effect of confined blasting, through blasting impact dynamics, fluid mechanics, fracture mechanics, and blasting test. The crack propagation model of the surrounding rock under the condition of confined blasting is established. He et al. [19] used the field test and numerical simulation to study the dynamic failure process between adjacent boreholes. A high-speed camera (HS) and digital image correlation (DIC) were used in the field test, and the strain field and crack growth process were obtained. Hu et al. [20] introduced the smooth particle fluid dynamics (SPH) algorithm and the improved blasting damage model into the explicit finite element program LS-DYNA with subroutine interface and developed the coupled SPH-DAM-FEM numerical simulation program. The modified damage model is used to couple the SPH technique to simulate the near zone in the blasting process, and the finite element method is used to solve the far-field response of the rock. Daniel Johansson [21] studied the effect of the delay time on the rock breaking effect.

Rock belongs to discrete materials, and rock fracture caused by blasting is a problem of large deformation of materials. The PFC method in existing simulation methods is not limited by material deformation, which can effectively simulate discontinuous phenomena such as medium cracking and separation and can well simulate the blasting effect of rock mass under explosive load, so it has good theoretical feasibility. The geometric elements of the blasting crater were studied through numerical simulation, and the conclusions obtained can better guide the blasting construction, which has great engineering significance.

2. Parameter Calibration of the Particle Flow Discrete Element Rock Mass Model

2.1. Basic Mechanical Parameters of Limestone

Limestone samples were taken from the test site of Langshan limestone mine. The basic mechanical parameters of limestone in the site were determined by the rock uniaxial compressive strength experiment. The bedrock of the mining area is bare, the limestone is compact and hard, the fractures are not developed, and the integrity is good. Samples with good integrity and homogeneity were selected during the test, and cylindrical specimens with a height-to-diameter ratio of 1 : 2 (50 mm × 100 mm) were prepared. The obtained mechanical parameter data are shown in Table 1.

2.2. Calibration of Rock Model Mesoscopic Parameters

The PFC method is used to study various rock problems, and the biggest difficulty encountered is the selection of mesoscopic parameters. The particle flow code model of uniaxial and biaxial compression is used to get the macroscopic parameters of rock materials, which is the main method to calibrate the mesoscopic parameters of the model [22].

First, a numerical experiment container was constructed by using PFC. The size of the container is the same as the sample size. When numerical simulation is carried out, the loading wall should have sufficient strength to ensure that it will not be destroyed by particles. The effective modulus of the loading wall was set to 100 Gpa.

Finally, according to uniaxial compressive strength, uniaxial tensile strength, and elastic modulus of the rock obtained in the laboratory, all the calibrated mesoscopic parameters of limestone are shown in Table 2. The numerical model of rock mass (8 m long and 6 m wide) is produced by using the mesoscopic parameter data as the follow-up numerical test model.

3. Numerical Simulation of Rock Blasting Based on the Particle Flow Discrete Element Method

3.1. Boundary Condition

This paper assumes that the rock mass model is the infinite medium model and ignores the reflection of stress waves at the infinite boundary of the model. Therefore, it is necessary to set the stress wave dispersion boundary in the PFC model to absorb the incident wave energy generated inside, so as to simulate the infinite medium. One method is to build a large-scale model to dissipate the energy generated by explosion through internal friction, local damping, and bonding contact of the model. However, this method requires a large number of particles, and it is difficult to solve the model. The second method is to fix the velocity and displacement of the particles at the boundary, so that the reflection of the stress wave at the boundary makes the calculation result of the model inconsistent with the actual situation. The third method is to apply the boundary force to the boundary particles to satisfy the absorbing stress wave energy.

In this paper, the viscous boundary proposed by Kouroussis and Verlinden [23] and the dispersion effect of stress wave propagation at rock mass boundary proposed by Shi [24] are considered. The boundary conditions of rock mass are established, as shown in Figure 1.

The relationship between boundary force and particle motion velocity iswhere r is the particle radius, ρ is the rock medium density, C is the wave velocity, and is the particle velocity.where ξ and η are the correction coefficients of P-wave and S-wave dispersion effect, respectively; Cp and Cs are the P-wave velocity and S-wave velocity of stress wave propagating in the model medium, respectively; and are the normal and tangential velocity of particles at the boundary, respectively.

3.2. Expansion Loading of Explosive Particle

When the explosive particle expands, it will produce superposition amount with surrounding rock particle system. According to the particle contact principle in PFC, the explosive particle will generate radial pressure on the surrounding rock particles when the explosive particle expands to the borehole wall. The peak pressure action coefficient of the hole wall α during explosive particle expansion was defined as the ratio of the hole wall peak pressure pm to the normal stiffness of explosive particle kn, i.e., α=pm/kn [25]. When the contact stiffness and explosion pressure are known, the peak value of particle radius variation iswhere dm is the expansion amount of explosive particle, r0 is the explosive particle initial radius, pm is the peak pressure of the hole wall, kn is the normal stiffness of explosive particle, and α is the peak pressure action coefficient of the blast hole wall.

The explosive loading propagates outward in the form of spherical wave, which is equivalent to pulse stress wave. In this paper, it simplified as a semisinusoidal wave with the same rising and falling time, and its expression iswhere f is the frequency of semisinusoidal wave, t is duration of blasting action, and p(t) is expansion pressure. The duration of blasting action is taken as 10 ms. Similarly, the semisinusoidal wave has a period of 10 ms, namely, 1/f = 10 ms.

The limestone numerical model adopted in the following factor analysis research is shown in Figure 1. Also, the mesoscopic parameters of the model are taken from the data shown in Table 2. According to equations (3) and (4), the explosive particle expansion method in PFC was used to simulate rock mass blasting.

3.3. Influence of Peak Pressure Action Coefficient on Blasting Effect

The buried depth of explosive particle is 0.5 m, the diameter of explosive particle is 32 mm, and the diameter of borehole is 50 mm. The peak pressure action coefficient α changed from 0.1 to 5.0. The radius and depth of the blasting crater were recorded after each blasting operation. The relationship between the peak pressure action coefficient α and the depth of the blasting crater is shown in Figure 2. Through linear regression analysis, the relationship between the depth of blasting crater h and the peak pressure action coefficient α is h = 0.6143α + 0.2775 and r2 = 0.839.

In addition, the relationship between the peak pressure action coefficient α and the radius of blasting crater R is shown in Figure 3. Furthermore, through regression analysis, the relationship between the above two parameters is R = 0.3107α + 0.4517 and r2 = 0.959.

Another important parameter often used in engineering blasting is blasting action index n, which is the ratio of blasting crater radius R to minimum burden W, which is expressed as n = R/W. The relationship between the peak pressure action coefficient α and blasting action index n is shown in Figure 4. Through linear regression analysis, the relationship between blasting action index n and the peak pressure action coefficient α is n = 0.1793α + 0.928 and r2 = 0.895.

With the increase in the peak pressure action coefficient, the blasting action index n also increases. The form of blasting action changes from weakened throwing blasting to normal throwing blasting and finally to enhanced throwing blasting. Figure 5 shows the numerical simulation of the blasting crater under three throwing blasting actions.

3.4. Influence of Charging Mode on Blasting Effect

The multiparticle expansion loading method was used to simulate the explosion of cylindrical explosive in rock mass, and the damage and blasting vibration of surrounding rock by two charge methods were compared in this part. The cylindrical charge uses five explosive particles. The buried depths of explosive particles are 2 m, 1.85 m, 1.7 m, 1.55 m, and 1.4 m, respectively. The diameter of explosion particle is 32 mm. The borehole diameter is 50 mm. The peak pressure action coefficient α of single explosive particle is 0.55. Relatively, the depth of explosive particle in the spherical charge is 2 m, the diameter of this explosive particle is 32 mm, and the diameter of the blast hole is 50 mm.

To simulate the difference of blasting effect between spherical charge and cylindrical charge in numerical simulation under the condition of the same charge quantity, the peak pressure action coefficient of spherical charge is 5 times that of single explosive particle in cylindrical charge. Therefore, the peak pressure action coefficient of spherical charge is taken as 2.75. The cylindrical charge structure in the rock mass model is shown in Figure 6.

The stress wave propagation process of cylindrical charge obtained by numerical simulation is shown in Figure 7, showing an obvious columnar wave.

Cylindrical explosive forms a crater depth of 1.99 m and radius of 1.25 m, while spherical explosive forms a crater depth of 2.4 m and radius of 2.08 m. Figures 8 and 9 show the blasting crater and rock crack of spherical and cylindrical explosives. It can be seen that, under the same charge quantity, the spherical charge forms more cracks and damage to the surrounding rock mass.

The particle position is shown in Figure 10. Figure 11 shows the curve of particle velocity in x direction with time at a distance of 3.9 m from the explosive particle. Under the same charge quantity, the peak velocity of the cylindrical charge is lower than the spherical charge velocity.

It is reasonable to adopt the method of simulating cylindrical charge by the multiexplosive particle expansion loading, which accords with the theory and practical engineering.

3.5. The Influence of Explosive Particle Buried Depth on Blasting Effect

In this part, the influence of buried depth on blasting effect is studied by adjusting the buried depth of explosive particle. The diameter of explosive particle is 32 mm, the borehole diameter is 50 mm, and the peak pressure action coefficient α is 0.75. The buried depth of explosive particle is changed. The radius and depth of the blasting crater and the number of discrete blocks generated after each blasting operation are recorded. The specific values are shown in Table 3. In this paper, the rock crushing index δ is defined aswhere N is the number of discrete blocks produced by blasting, R is the radius of blasting crater, h is the depth of the blasting crater, and A is the sectional area of the blasting crater.

Figure 12 shows the variation curve of the number of generated blocks N with the buried depth h of explosive particle. When the buried depth is 0.7 m, the number of generated blocks reaches the largest.

Figure 13 shows the variation curves of the sectional area of the blasting crater A and the crushing index δ with the buried depth h of explosive particle. The figure is divided into three parts: (1) 1.0 m–0.7 m is the impact crushing zone. When other conditions are unchanged in this area, the sectional area of the blasting crater increases with the decrease in the buried depth of explosive particle. When the depth reaches 0.7 m, the area size of the blasting crater and the crushing index reach the maximum of the whole curve. The optimum depth of the explosive particle buried depth is 0.7 m. (2) When the depth of blasting particle is gradually reduced, that is, 0.7 m–0.3 m, the rock is in the fracture zone. When the sectional area of the blasting crater is reduced, the fracture index is first decreased and then increased. The crushing index reaches the maximum of the whole curve at the buried depth of 0.3 m, and the area reduces to the minimum value of the fracture zone. The energy is consumed in rock crushing and throwing. At this time, the buried depth of the explosive particle is the transition depth. (3) When the explosive particle is buried shallowly, that is, 0–0.3 m, the particle is in the air explosion zone, and the little energy is used to break the rock. Below 1.0 m, it is an elastic deformation zone. The blasting action occurs only inside the rock mass and does not cause damage to the surface rock.

According to Hausser’s formula, the charge quantity Q of standard throwing blasting is calculated:where q is the specific charge and W is the minimum burden. Figure 12 shows the explosive optimal burial depth is 0.7 m. It can be considered that the blasting carter when the explosive explodes at a depth of 0.7 m is the standard throwing blasting crater. Therefore, it can be approximately considered that the minimum burden is 0.7 m and W = 0.7 m.

In terms of the relationship between critical depth and charge quantity in Livingston blasting funnel theory,where We is the explosive critical depth and Eb is the rock deformation energy coefficient. According to site rock conditions, the specific charge q = 1.3 kg/m3 and the deformation energy coefficient Eb is 1.3. Therefore, through equations (6) and (7), it can be appropriate that when the minimum burden is 0.7 m, the charge quantity Q is 0.4459 kg and the critical depth is 0.993 m. The critical depth obtained by theoretical calculation is close to that obtained in Figure 13.

When the buried depth of explosive particle changes, the blasting effect also changes. With the decrease in buried depth, the deformation and failure modes of rock blasting are divided into four parts, which is consistent with theoretical research. The numerical simulation of particle flow code method can be used to determine the depth of explosive according to the characteristics of rock, which has important engineering significance for improving the blasting effect.

4. Example

According to the peak pressure action coefficient α, the linear regression equation of the blasting crater depth is

By taking the minimum burden equal to the buried depth of the explosive particle, W = h, we obtain

The charge quantity Q is 60 g and 80 g. Because the charge quantity is small, it can be equivalent to spherical charge. According to the on-site mine blasting design scheme, the specific charge q is 0.35 kg/m3. According to equation (9), when the charge quantity is 60 g and 80 g, the peak pressure action coefficient α is 0.45 and 0.54, respectively.

The blasting crater and rock cracks after numerical simulation are shown in Figures 14 and 15. The specific values of the blasting effect are shown in Table 4.

5. Field Blasting Experiment

5.1. Blasting Effect Test Results of Rock Mass

The field blasting experiment of single hole was carried out on the site. The blasting effect of field test is shown in Figure 16. After the blasting action, the depth and diameter of the blasting crater were measured using a tape measure. The specific data are shown in Table 5.

5.2. Comparison of Field Experiment and Numerical Simulation

The specific calculation methods involved in the comparative analysis in this part are as follows:(1)The calculation method of blasting crater measured average value in field test is defined as follows:where is the measured average value, is the measured value, and G is the number of experimental groups.(2)In order to compare the results of field test and numerical simulation, the calculation method of error rate is defined as follows:where Er is the error rate between the measured value and numerical value and S and F are the values obtained by numerical simulation and field test, respectively.

According to equations (10) and (11), the values of different charge amount experiments can be calculated, and the analysis results are shown in Table 6.

According to the data in Table 6, it can be concluded that the error rate between the field test and numerical simulation results is controlled within ±10%. Therefore, it shows that the PFC method of explosive particle expansion loading has good reliability and can be used to guide filed blasting construction.

6. Conclusions

In this work, a numerical model of limestone rock mass was established by PFC. The explosive particle expansion loading method was used to simulate rock blasting. Then, the blasting effects of peak pressure action coefficient, charging mode, and explosive buried depth were investigated. Finally, the numerical simulation results were compared with the size of blasting crater obtained by field test. The main conclusions are as follows:(1)In the numerical simulation, with the decrease in buried depth, rock deformation and failure caused by blasting can be divided into four parts (elastic deformation zone, impact fracture zone, fracture zone, and air explosion zone), which is consistent with the theoretical research. Moreover, the explosive particle loading method is used to simulate cylindrical charge, which confirms to the theory and practical engineering.(2)Combining engineering blasting theory with the PFC method in DEM theory, explosive particle expansive loading method is adopted to simulate rock mass blasting. The correctness and reliability of the numerical simulation are verified by the error within ±10% between the field blasting test results and the numerical simulation results. Furthermore, compared with the FEM, the DEM using PFC is more in line with the actual rock mass. This study provides a new reliable numerical simulation method for rock mass blasting. According to rock characteristics and numerical simulation, the reasonable explosive buried depth and charge quantity can be determined, which is of great engineering significance for improving blasting.(3)The numerical simulation of blasting crater size in this paper is in good agreement with the field results. In the numerical simulation, the blasting effect under different charging quantity is explored indirectly by adjusting the action coefficient of the peak pressure. However, there is no theoretical relationship between the action coefficient of peak pressure and charge quantity. So, there are still some limitations in practical engineering. In the future research work, it is still indispensable to improve the explosive particle expansion loading algorithm in PFC procedure, so that the algorithm can directly establish the corresponding relationship with the actual engineering charge. In addition, the number of field tests should be increased to make this numerical simulation more convenient to be applied in engineering practice and to guide blasting design scheme better.

Data Availability

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

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (no. 51874189).