A Novel Continuous-Discontinuous Multi-Field Numerical Model for Rock Blasting

: During blasting, rock failure is caused by blasting wave and explosive gas pressure, as a multi-ﬁeld coupled process. Most numerical models focus on the effect of blasting wave where the gas pressure is commonly accounted for by empirical relations, ignoring the penetration and permeation of gas ﬂow in cracks. This can underestimate the failure region. In this work, a novel multi-ﬁeld model is developed in the framework of a continuous-discontinuous element method (CDEM), which is a coupled ﬁnite-discrete method with explicit integration strategy. The deformation and cracking of rock mass and the distribution of gas pressure are captured. The proposed method is veriﬁed by comparing the results to other results provided in published literature. Especially, by simulating the cases with blocked and unblocked blasting hole, we found that: (i) The fracture degree of the case with blocked blasting hole was 30% higher than that of the unblocked blasting hole. (ii) The radial main cracks in the fracture area are mainly caused by the explosive gas, and the tiny and dense cracks near the hole are induced by the explosion stress wave. (iii) The explosion crushing zone is mainly formed by the action of explosion stress wave, while the crack zone is formed by the combined action of the explosion stress wave and explosive gas. The proposed method provides a useful tool to properly simulate a rock blasting process.


Introduction
Considering rock blasting, the failure of rock mass is induced by blasting wave and successive gas pressure. When the former effect lasts for several microseconds, the latter effect can last for up to 1 millisecond. This whole process is short and accompanied by rapid energy transformation. The effect of blasting wave is widely studied, providing practical empirical formulas for applying equivalent time dependent dynamic loads. However, some empirical equations partly take into account the effect of explosive gas pressure. The gas flow can penetrate and permeate into cracks which promote the formation of main cracks. This process is similar to hydraulic fractures. In other words, ignoring this coupled hydro-mechanical process commonly leads to underestimation of damage regions.
Numerical simulations can assist the researchers to study the rock blasting efficiently and quantitatively. It is difficult for conventional finite element method (FEM) to capture initiations and propagations of cracks of strain-softening materials such as rocks. However, some novel numerical methods are proposed in the framework of continuum mechanics such as numerical manifold methods and extended finite element methods [1][2][3][4][5], cracking element methods (CEMs) [6][7][8][9][10], phase field methods [11][12][13], cracking particle methods [14][15][16], and peridynamics [17][18][19][20][21]. Blasting lasts for a very short period, which commonly results in local fragmentation and large displacements. The hybrid continuousdiscontinuous method with explicit integration can be advantageous considering the computing time and numerical stability. For example, the authors of Fakhimi and Lanari [22] simulated a rock blasting process based on a DEM-SPH coupling method, for which the SPH model is used to capture the gas flow and DEM is adopted to simulate rock fractures. The authors of Jayasinghe et al. [23] built an FEM-SPH coupling method to simulate rock blasting and found that the discontinuity and in situ stress of rock have significant influence on the blasting results. The authors of Wang et al. [24] used a hybrid continuous-discontinuous method to simulate the rock blasting induced by high-pressure gas and explosives, indicating that the former case provided fewer radial fractures.
The above research shows that the hybrid continuous-discontinuous method has unique advantages in the field of explosion simulation. In this work, in the framework continuous-discontinuous element method (CDEM) [25][26][27][28][29][30], we build a coupled hydromechanical multi-field model for simulating the rock blasting considering the coupled effects of blasting wave and gas pressure. The main features of our model include:

•
The domain is discretized into deformable blocks and interfaces which can capture the continuous-discontinuous process of rock blasting. • The explosive gas flow is simulated with Darcy's law and plate flow, where the opening of cracks directly increase the gas permeation coefficient. • A multi-field coupled iterative procedure is designed based on explicit integration framework for assuring efficiency and reliability.
By comparing our numerically-obtained results and experimental results provided by other literatures, our model is validated, which properly predicts the damage region of rock blasting considering blocked and unblocked blasting holes.
The paper is organized as follows: In Section 2, the framework of the continuous-discontinuous element method and multi-field coupling model are introduced, including control equations and numerical procedures. Subsequently, the equivalent loadings are described, including the dynamic loads induced by blasting wave and the quasi-static loads caused by the explosive gas pressure. In Section 3, a classical blasting test is used to verify the model. In addition, the single hole blasting of PMMA(Polymethyl Methacrylate) plate is simulated as the application, comparing to the experimental results considering blocked and unblocked blasting holes. Concluding remarks are given in Section 4.

Continuous-Discontinuous Element Method
The continuous-discontinuous element method (CDEM) is a numerical analysis framework with explicit integration built based on the Lagrange equation [31][32][33]. The elements adopted in this work include deformable block and interface elements, see Figure 1. In Figure 1, the solid black lines represent the real interface such as joints and other initial cracks. The dashed black lines represent the virtual interfaces, which will become activated when damage happens. Normal and tangential springs are arranged in the virtual interface. When these springs reach their strengths, they will break and the interfaces will open. The governing equation of CDEM method is as follows: • Block element The governing equation is established based on the Lagrangian energy system.
The energy function of element is The system external force includes damping and boundary external forces: Hence, Equation (1) can be rewritten as: Finally, the block element dynamic function is • Interface element An interface locates between two block elements. The normal and tangential contact forces of the two elements are the spring forces as The relative displacements between adjacent elements are The Mohr-Coulomb criterion is adopted; see Equation (9). When the normal stress of the spring is greater than the tensile strength, the normal spring will break and tensile failure will occur. Similarly, when the tangential spring reaches its strength, the tangential spring will break and shear failure will occur:

Permeation of Gas in Cracks
For obtaining the distribution of gas pressure, the finite volume method is used. When the springs connecting the adjacent elements break, Darcy's law and plate flow are considered to simulate the gas permeation in these cracks. The velocity of gas flow V is obtained with where k s is the gas permeation coefficient in crack. Considering the plate flow and cubic low, we have where ω is the normal crack opening, and µ g is the viscosity of gas. An iterative procedure is designed to trace the gas flow in the cracks. When the crack opens, the gas pressure of the nodes on the two crack surfaces is updated. The gas flow permeates further and results in increasing gas pressure in the deeper region. This iterative algorithm is shown in Figure 2.

Explosion Load
Yuan et al. [34] suggests a time-dependent explosion pressure curve with two peaks, see Figure 3a. The first peak corresponds to the blasting wave, and the second peak corresponds to the gas pressure. As mentioned before, comparing to the explosion stress wave which acts as dynamic loading, the explosive gas permeation process lasts longer, which is considered as a quasi-static loading. The first peak pressure P w can be determined base on C-J detonation theory, which is Furthermore, the action time t w According to [35], time of blast pressure ascent stage t wr = 0.1 t w and descent stage t w f = t w − t wr .
The second peak pressure of explosive gas P gmax was calculated based on Park et al. [36]: (14) where ζ is the uncoupling coefficient of blasting hole, ζ = r b / r e where r b is the radius of hole, and r e is the radius of charge. Based on [37], we assume a sine curve for the second peak. The complete time dependent explosion pressure curve used in this work is shown in Figure 3b.
A personal computer with an AMD Ryzen Threadripper 3960X 3.90 GHz processor and 32 GB memory was used to obtain the numerical results. A benchmark model with 10,000 triangular elements will take 2 h, and a complex case such as Section 3.2 of this paper (with 21,896 triangular elements and 32,972 interface elements) will take about 6 h.

Single Hole Blasting Test
To verify the proposed model, the single hole blasting test is considered. The model is a granite disk with radius 72 mm and a hole with a radius 3.225 mm, which is filled with PETN (Pentaerythritol Tetranitrate) [38], as a plane stress condition. The material parameters are taken based on [39]; see Table 1 and the parameter for explosive is listed in Table 2. In total, 39,480 triangular elements and 59,672 interface elements are created in the numerical model.  [38].

Material Properties Value
Charge density ρ w (kg/m 3 ) 1770 Specific internal energy of explosive U (J/m 3 ) 1.010 × 10 10 C-J pressure P CJ (GPa) 3.2 Detonation velocity V (m/s) 8300 The obtained cracking pattern is shown in Figure 4, compared to the MPM-CDEM results provided by [39] and experimental results provided in [38]. Figure 4 indicates that the explosion pulverized formed a pulverized zone around the blasting hole, and some radial cracks propagate to the edge of the model. Close to the edge of the model, some tangential cracks appear due to stress wave reflection, resulting in a nearly circular fractured region. Assuming the pulverized zone as the blasting funnel, the radius of the blasting funnel and number of main cracks are listed in Table 3, indicating the reliability of our model.

Blasting of the PMMA Plate
We used the proposed method to simulate the blasting test of the PMMA plate [40], see Figure 5. This case is also plane stress. We used Gmesh software to build a two-dimensional square model with a side length of 300 mm and set a blasting hole with a radius of 2 mm in the center of the model. We divided the model into 21,896 triangular elements and 32,972 interface elements, and the size of elements is 1-5 mm. To prevent the explosive stress wave reflection, we set a non-reflective boundary condition around the model. The material properties are listed in Table 4, which is determined based on [41]. The explosive is PbN 3 , the properties of which are listed in Table 5.   Table 5. PbN 3 main parameters.

Material Properties Value
Charge density ρ w (kg/m 3 ) 2560 Explosion heat Q (kJ/kg) 1524 C-J pressure P CJ (GPa) 3.59 Detonation velocity V (m/s) 2250 According to [40], two conditions are considered as (i) blocked blasting hole and the (ii) unblocked blasting hole, where the gas permeation in the second condition can be ignored in the experiment. The obtained results are shown in Figure 6. Figure 6a,b show the numerically-obtained results, and Figure 6c,d show the experimental results. It can be seen that, when the hole is not blocked, there is almost no radial main crack. Only dense cracks are generated near the hole, resulting in a crush zone caused by the blasting wave. When the hole is blocked, several radial main cracks are generated, which are caused by the gas permeation. The numerical results agree well with the experimental results, indicating the reliability of the proposed model. We would like to emphasize that we did not use special techniques to simulate explosion crushing area such as element dissolution/deletion because such techniques do not follow mass and energy conservation, which may bring other problems. Hence, in our results, a large number of micro-cracks appeared on the edge of the explosion crushing area in Figure 6. Considering the blocked condition, the evolution of gas pressure is shown in Figure 7. The damage region at t = 30 µs in Figure 7 is similar to the final damage region of the case with an unblocked blasting hole, indicating that the effect of blasting wave has finished. At t = 30.22 µs, the explosive gas begins to permeate into the crush zone with maximum gas pressure around 20 MPa. When t = 39.82 µs, the maximum gas pressure exceeds the strength of the PMMA plate. At t = 56.32 µs, the main cracks appear, and the gas begins to escape. The maximum gas pressure also drops. Finally, at t = 58.12 µs, the cracks become smooth, and the whole process ends.
The results of Mises stress are shown in Figure 8. In the beginning of the explosion, the blasting hole area was subjected to strong impacts where the stress peak reached 100 MPa. Before 30 µs, the explosive stress wave spread almost symmetrically. It can be found that the blasting wave almost reached the boundary of this region in 30 µs. With the propagating of the explosion wave, crushing and cracking zones were formed. When t > 30 µs, the explosive gas permeated into cracks, promoting the cracking processes. At t = 36.82 µs, the cracked zone was almost formed where the stress concentrations could be found at the tips of cracks. At t = 60 µs, the explosive gas began to dissipate. (The blue grid is the fracture grid which supports the explosive gas flow in the crack). When defining the ratio of the number of broken springs to the total number of springs as fracture degree, the history of fracture degree is shown in Figure 9. It can be seen that, in the first 30 µs, the fracture degree of the cases with blocked and unblocked blasting holes are almost the same. After this, the fracture degree of the case with a blocked hole continuously increases, which is 30% higher than the case with the unblocked blasting hole in the final stage. This conclusion is consistent with Ref. [40], and it also verifies the proposed model's rationality. However, due to the calculation efficiency, the proposed model is only suitable for the two-dimensional case, and the three-dimensional case will take too much time in calculating the explosive gas flow in the crack. This is also our future work.

Conclusions
A coupling method is proposed to effectively simulate explosive cracks in rock. The plate flow equation is used to describe the process of gas expansion and pressure propagation in explosion; CDEM is used to describe the process of blast stress wave shocking and crushing rock. By combining the dynamic and quasi-static loads, the progressive failure process of rock from continuous to discontinuous is reproduced. Compared to previous studies, the role of explosive gas in the explosion is considered, providing more reliable results of damage regions. The effectiveness and accuracy of the method in blasting are proven by numerical examples. The main conclusions of this paper are as follows: 1.
The whole process of rock blasting including the blasting stress wave, explosive gas permeation, and fracking processes are reproduced by the proposed model.

2.
For single-hole PMMA plate blasting, the established numerical model captures the crush and fracture zones. By comparing the two cases with blocked and unblocked blasting holes, it is found that the crushing zone in the explosion is mainly caused by the explosion stress wave. The tiny and dense cracks in edge of the fracture zone are caused by the reflection of the stress wave at the boundary. Meanwhile, the radial main cracks developed in the specimens are caused by the explosive gas. By comparing the fracture degrees, it can be found that the explosive gas can account for around 30% of the total damage. 3.
The explosion crushing zone is mainly formed by the action of the explosion stress wave, while the crack zone is formed by the combined action of the explosion stress wave and explosive gas.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
List of symbols