Numerical Investigation on Crack Formation and Penetration Mechanism between Adjacent Blastholes

A deep understanding of the crack formation mechanism between blastholes had great significance for improving the energy utilization rate of explosives. From the perspective of static mechanics of stress wave theory, this paper theoretically derived the stress distribution on the blasthole connecting. It was proven that the stress at the midpoint of the blasthole connecting was not the maximum. The analysis results contradicted the original theoretical results. Moreover, the finite element software LS-DYNA was used to numerically simulate the crack formation between adjacent blastholes in infinite rock media. The fluid-solid coupling method was used to simulate the effect of stress wave and detonation products on rocks. The simulation results were consistent with the actual situation and they showed that cracks were formed in the blastholes wall firstly. Stress wave superposition was not the main reason for crack penetration. Stress wave reflection stretching at the crack tip was vital for crack penetration.


Introduction
It had great significance for the design of blasting parameters in mining and roadway driving to have a deep understanding of the variation law of explosive stress and strain field between adjacent blastholes and explain the mechanism of crack penetration between blastholes correctly [1,2]. For blasting engineering problems, the biggest concern was how to qualitatively describe and quantitatively calculate the deformation and failure of rock and soil media [3,4]. Previous studies on rock and soil blasting were mostly based on theoretical analysis and experimental research. Wang used laser dynamic caustics to analyze the dynamic fracture effect of detonation crack of different charge structures and to examine the directional fracture-controlled blasting mechanism of slotted cartridge [5]. Using caustics method, Qiu found that vertical P 1 P waves tended to cause crack arrest, while horizontal P 2 P waves assisted in crack opening [6]. In the same experiment method, an optically geometrical superposition of light deflections from the running crack and blasting waves was proposed. e stress field near the running crack tip in crack-wave interaction was obtained [7]. Using slit charge blasting method, Yang concluded that crack propagation was suppressed and the length of crack was short when the static stress was perpendicular to the slit direction [8]. Complete theoretical analysis was only applicable to some simple and ideal problems. Wang obtained the cumulative damage model of rocks mass through theoretical derivation. Experimental study indicated that blast activated and then extended the initial cracks in rock mass, leading to accumulation of rock mass damage. e rock mass damage accumulation can be conveniently quantified using the proposed damage variable. Meanwhile, the instantaneous nature of explosion process, the complexity of acted medium, and the limitations of test technology made researchers unable to obtain sufficient spatial and temporal data [9]. Up to now, the determination of blasting parameters in engineering calculation mainly relied on experimental formulas, mostly obtained by sorting out experimental results on the basis of energy principle and geometric similarity law [10,11].
Numerical simulation technology provided a new way for the study of explosion problem in rock and soil. is technology was used to predict physical behavior of the deformable objects with very high level of accuracy, such as the finite element method (FEM) [12][13][14][15][16]. Compared with theoretical analysis and experimental study, this method could ignore the limitation of geometric size of the research object and visually showed the phenomenon which was difficult to observe in the experiment. Numerical simulation had become an indispensable means to study the problems in explosion field. Seelig investigated the crack trajectory under the influence of stress waves by adding new boundary elements of constant length at the running crack tip [17]. Yi carried out numerical simulation on the stress wave under the condition of differential blasting of adjacent blastholes [18]. Yang used particle flow code method to obtain the stress wave gradient law, and the effects of blasting rock breaking, crack extension, and explosion effect outside the blasting crater area were determined with different initiation modes [19]. Yuan studied the effect of shock wave and detonation gas on fracture formation by using numerical simulation method. e results showed that the shock wave played a major role in crushing the rock near the wall of the blasthole and manufacturing a few primary radius cracks further away from the blasthole, whereas the detonation gas further extended the fractures to increase the crushed area [20]. Cho coupled the dynamic finite element and finite difference method to analyze the dynamic fracture process and gas flow through the fractures [21]. Lanari [22] used a hybrid two-dimensional finite element-discrete elementsmoothed particle program to model rock blasting. By controlling whether the gas penetrated into cracks, they clarified that the shock wave was the cause of most cracks, whereas the gas mostly separated the fragments. Mohammadi [23] combined the finite element method and finite difference method to investigate the solid behavior during blasting. e results indicated that gas flowed into opening cracks and enlarged the fracture networks. Yi studied the stress wave propagation between the double blastholes by numerical simulation. It was believed that the stress wave superposition did not completely determine the formation of cracks between blastholes [18].
Based on the above narrative, this paper used LS-DYNA software and adopted fluid-solid coupling method to study the response process of rock under the action of doubleblastholes explosion load in infinite medium.

Traditional Stress Wave eory.
e traditional stress wave theory believed that the rock elements were subjected to radial compressive stress and tangential tensile stress under the action of stress wave generated by the explosive in blasthole. Figure 1(b) is the same as (a) in the process of stress wave's propagation. Figures 1(a) and 1(b) show the stress state of the rock element in the case of single-hole charging. For Figure 1(c), after the double blastholes were simultaneously detonated, the stress wave propagated along the blastholes connecting and converged at the midpoint of the connecting, where the tensile and compressive stresses reached the maximum. When the tensile stress reached the rock allowable value, the rock produced cracks. is theory held that when two adjacent blastholes were detonated at the same time, the crack was initially produced in the middle of the two blastholes' connecting.
By analyzing the traditional stress wave theory, it can be found that the core of this theory was a simple stress superposition. After explosion of an explosive, rock in the middle of two blastholes was more stressed than blasthole wall due to the simultaneous action of the two blastholes. Although the dynamic propagation process of stress wave was mentioned in traditional stress wave theory, the variation of stress wave in this process was not considered in practice. erefore, the theory was essentially formed from the perspective of statics. e traditional stress wave theory was verified in the view of stress superposition.

Statics eory Calculation
Research. According to the knowledge of elastic mechanics, when uniform load q was applied to the blasthole wall in the infinite medium, the stress distribution in the surrounding medium can be expressed by the following formulas: In the formulas, a is the radius of the blasthole, r is the distance from the center of the blasthole, σ r is the radial stress, and σ θ is the tangential stress. If there were two blastholes and the hole spacing was b, it can be concluded that the stress distribution after superimposition on the double-blastholes connecting was When x � 0, that is, the selected position was the midpoint of the blastholes connecting, at this time, Usually b � (16 ∼ 24)a; it can be seen that the stress at the midpoint of the double-blastholes connecting was still much smaller than the stress at the blasthole wall. Figure 2 is a plot of the stress distribution on the blasthole connecting at b � 20a.
It can be seen from Figure 2 that when the doubleblastholes detonated at the same time, the superposition of intermediate stresses among double blastholes was very weak, and the intermediate stress was still much smaller than the wall of blasthole. e above derivation results were contrary to the view that the stress in the midpoint of the blasthole connecting was largest and the crack occurred firstly in this position in the traditional stress wave theory. is derivation was based on the static thinking of the traditional stress wave, so it can prove that the traditional stress wave theory itself had contradictions.

Physical Model.
Since the cylindrical charge was used in the actual blasting operation and the charge diameter was far less than its length, the model can be simplified as a problem of finite element at plane without considering the special blasting effect at both ends of the charge.
All parts of the model adopted solid units. e size of the whole calculation field was 300 cm × 200 cm × 0.5 cm, the diameter of the blasthole was 5 cm, and the diameter of the explosive was 3 cm. e centers of the two blastholes were located at 125 cm and 175 cm in the x direction, respectively, and the explosives in the two blastholes were detonated simultaneously. It was assumed that the blasthole had sufficient depth, so that the material at the cross section away from the orifice was hardly displaced and deformed in the depth direction of the hole, so only a single-layer grid needed to be divided during the model building process, and all nodes in the model were constrained parallel to the freedom of the blastholes. No reflective boundary conditions were applied around the rock. e model is shown in Figure 3. e peripheral marks of the rock represented displacement constraints. In the center of the figure are two blastholes 50 cm apart.

Algorithm and Material Model.
is problem belonged to contact explosion. For contact detonation, although the Lagrange algorithm could be used to define the contact between the explosive and the rock to consider the interaction, the explosive element was prone to serious distortion during the explosion. Fluid-solid coupling mechanics was a branch of mechanics generated by the intersection of fluid mechanics and solid mechanics. It was a science to study the interaction of deformed solids under the action of flow field and the influence of solid shape on flow field. In this paper, the fluid-solid coupling method was used to calculate the effect of explosive on rock. e model contained three materials: explosives, air, and rock. e air and explosive adopted Euler algorithm, while the rock adopted Lagrange algorithm. Because of the high strain rate of air and explosives, the mesh was fixed in Euler algorithm, which can reduce the influence of mesh distortion on the calculation process. In the Lagrangian algorithm, the meshes moved together with the rock units, which can display the mechanical behavior of the rock more accurately. During modeling, air should overlap with rock near the blastholes. In this paper, the air domain diameter of each blasthole was taken as 3 times the diameter of the blasthole, and a small width of air domain was established to connect the two blastholes. Figures 4 and 5 are, respectively, the distribution diagram of air domain and the enlarged grid division diagram. e mesh of rock and air elements in a small area around the blastholes was obtained by sweeping method. e mesh of other parts and explosive elements were acquired by mapping method. e element type was SOLID164.
Explosives, air, and rock were coupled by * CONSTRAINED_LAGRANGE_IN_ SOLID [24]. In LS-DYNA software, there were two main methods to simulate the generation and propagation of cracks. e first was to generate cracks in the structure through the failure of elements, and the other was to generate cracks through the defined node constraint failure. In this section, the first method was used to study the failure form of rock between the blastholes by defining the keyword * MAT_ADD_EROSION. Element destruction criteria Shock and Vibration were added in * MAT_ADD_EROSION. When the stress and strain of the element in the finite element model exceeded the set value, the failed element was removed and cracks were formed in the model. e HJC model usually adopted the equivalent strain failure f S failure mode, but this mode cannot display the failure caused by stretching and shearing. Better simulation results were achieved by adding the maximum compressive failure strain and maximum shear strain failure criteria in the * MAT_ADD_EROSION keyword. When using this algorithm, in order to improve the calculation accuracy, the mesh of the overlapping part of the air and rock was encrypted, as shown in Figure 5.
Explosives were described by using the HIGH_-EXPLOSIVE_BURN material model and the JWL equation of state. JWL equation was expressed as where E � ρ 0 e is the internal energy of the initial volume of the unit. A, B, R 1 , R 2 , ω are constants determined by experiments. P is the pressure generated by the explosive explosion. V is the specific volume. e material parameters of explosives are shown in Table 1. ρ is the material density. D and P CJ are detonation velocity and CJ explosion pressure of explosive, respectively. Air was described by using the MAT_NULL material model and the EOS_ LINEAR_ POLYNOMIAL equation of state. e equation is expressed as where C 0 ∼ C 6 are coefficients of equation of state. E 0 is the initial internal energy in unit volume. V 0 is the initial relative volume. Material parameters of air are shown in Table 2.
e rock adopted the HJC model, which was mainly applied to the simulation of concrete and rock under high  strain rate and large deformation. is model was suitable for Lagrange and Euler algorithms and had been widely used in numerical simulation. e strength of HJC model was described by normalized equivalent stress: where σ ′ � (σ/f c ′ ) is the ratio of actual equivalent stress to static yield strength (where σ ′ is equivalent stress, σ is actual equivalent stress, and f c ′ is the quasi-static uniaxial compressive strength). P ′ � (P/f c ′ ) is dimensionless pressure (P refers to actual pressure). ε ′ � (ε/ε 0 ) is dimensionless strain rate (where ε is the actual strain rate and ε 0 � 1.0 s − 1 is the reference strain rate). D is the degree of damage. A, B, N, and C are the strength parameters of the material. Damage factor D (0 ≤ D ≤ 1) was obtained by the accumulation of equivalent plastic strain and plastic volume strain.
where Δε P is equivalent plastic strain increment. Δμ P is equivalent volume strain increment. ε f p and μ f p are equivalent plastic strain and plastic volume strain of crushing under atmospheric pressure, respectively. f(P) � ε f P + μ f P � D 1 (P * + T * ) D 2 is the plastic strain of material fracture under normal pressure P. P * and T * are the normalized pressure and the maximum normalized tensile hydrostatic pressure that the material can bear. D 1 and D 2 are the damage constant. In LS-DYNA, the HJC model was defined as * MAT_JOHNSON_HOLMQUIST_CONCRETE. rough extensive literature review, the most widely used rock parameters were selected. e parameters of rock materials used in this paper are shown in Table 3. Figure 6 shows the distribution of cracks around the blasthole at different times during the crack propagation process.

Process of Cracks Propagation.
As can be seen from Figure 6, the crack propagation of double blastholes in infinite rock medium can be divided into three stages: (1) formation of crushing zone at the blasthole wall; (2) the formation of fracture zone around the blasthole; (3) the penetration and expansion of cracks between blastholes.
e simulation results agreed with the observed results in reality.

Analysis of Mechanical Parameters of Blasthole Connecting.
is paper studied the stress involved in the theory of strength of materials mechanical. e four rock elements A, B, C, and D were selected at equal intervals from the blasthole wall to the midpoint of the connection line of the blasthole to observe various stress changes. Figure 7 is the location diagram of the taken element. Figures 8 and 9 are the time distribution curves of the maximum principal stress σ 1 and the minimum principal stress σ 3 of the three elements. e simulation in this section used the cm-g-μs unit system, the time unit was μs, and the stress unit was 10 11 Pa. According to the stipulations in material mechanics stress theory, tensile stress was positive and compressive stress was negative. e peak value of each curve in Figure 8 is the maximum tensile stress value of the corresponding element. e peak value of each curve in Figure 9 is the maximum compressive stress of the corresponding element. It could be seen from Figures 8 and 9 that the maximum tensile stress and maximum compressive stress at four elements A, B, C, and D satisfy A > B > D > C. Analysis of the absolute value of the peaks of the curves in Figures 8 and 9 showed that the maximum compressive stress value of each element was about 2-3 times the maximum tensile stress.
Observing the measured curves of each element in Figure 8, it could be found that when the stress wave propagated to each selected element on the blasthole   Table 2: Material parameters of air.
Shock and Vibration 5 connecting, the maximum principal stress on the element was of negative value firstly and then turned to positive value after a few microseconds. erefore, it could be known that after the explosive exploded, the rock on the connection line of the blasthole was in a three-way compression state during the propagation and action of the stress wave, and then the maximum principal stress quickly converted into tensile stress within a few microseconds. In Figure 9, during the stress wave propagation process, the minimum principal stress of each element first showed a negative peak and then rose to fluctuate around zero value.
Earlier, the paper analyzed the maximum principal stress σ 1 and the minimum principal stress σ 3 of the rock on the blasthole connecting, according to the maximum shear stress formula: It was not difficult to find that the maximum shear stress on the four elements A, B, C, and D also satisfied A > B > D > C. Figure 10 shows the maximum shear stress on the four measurement points of A, B, C, and D.
rough the above analysis, no matter what kind of stress, there is no maximum stress at the midpoint of the blastholes connection. So the crack would not appear first at the midpoint of the line. e simulation results showed that the peak values of various stress parameters about element C were the smallest among the four measurement elements.
is was because the stress wave propagated in the rock and decayed rapidly. After reaching the midpoint, the stress wave and the stress wave generated by the adjacent blasthole were superimposed. e superposition of stress waves at the midpoint caused the stress in D element to be slightly greater than that in C element. Figure 11 is a graph showing the distribution of the absolute value of the stress peak on the blasthole connecting.

Stress Analysis of Penetrating Cracks.
We selected three elements, A, B, and C, near the midpoint of the throughcrack according to Figure 12 for research. Figures 13 and 14 are the equivalent stress and maximum shear stress curves, respectively. In Figure 13, the three measuring points were closer to the midpoint of the blasthole connecting, so when the stress wave initially swept through the measuring points, there was only a small fluctuation in the equivalent stress value. At 65 μs, the two stress waves met and superimposed at the midpoint of the blasthole connecting. e equivalent stress value at the measuring point rose sharply to the maximum peak, but the element had not failed. After 65 μs, the two stress waves continued to move towards each other. After 100 μs, the equivalent stresses of the three measuring points under the action of the stress wave generated by the adjacent blastholes had shown a large upward trend and had successively failed.
It was found from Figures 13 and 14 that the failure time of the three elements A, B, and C near the midpoint of the crack was not around 65 μs. e equivalent stress and maximum shear stress of these three elements when they failed were less than the peak stress value when the stress wave met at the midpoint of the blasthole connecting. is showed that the reason why the cracks penetrate each other was not that the stress waves met and superimposed at the midpoint of the blasthole connecting. Shearing force was not the main factor for cracks penetration.
In the past, it was generally believed that the mechanism of cracks formation between blastholes was related to the tangential tensile stress. erefore, the maximum principal stress of the three elements A, B, and C was compared with the tangential (Y-direction) tensile stress. Figures 15 and 16 are the maximum principal stress curve and the Y-direction stress curve of elements A, B, and C, respectively.
It can be seen from Figures 15 and 16 that when the stress wave is transmitted to the selected element, the Y-direction stress and the maximum principal stress of the element are both compressive stress and then convert into tensile stress. In the first 100 μs time period, the stress curves in Figures 15  and 16 were almost completely consistent, indicating that the Y-direction stress acting on these three elements during this time period was the maximum principal stress. After 100 μs, the maximum principal stress and the Y-direction stress curve both maintained an upward trend, but the maximum principal stress increased significantly more than  is showed that, after 100 μs, the element has tensile stress in other directions besides the tensile stress in the Y-direction. e combination of the two stresses caused the maximum principal stress on the element to be greater than the Y-direction stress.
It can be seen from Figure 17 that after the stress wave was transmitted to each element, the X-direction of each element was in a compressed state for a long time. From 95 to 100 μs, the stresses in the X-direction of elements A, B, and C were successively converted into tensile stresses.
is was because the stress wave generated by a blasthole propagated to the tip of crack generated by adjacent blasthole to produce a radial reflection stretching effect. Combined with the above analysis, the generation of this radial tensile stress caused the maximum principal stress at the crack tip to exceed the Y-direction tensile stress.
Comparing the stress values in Figures 15 and 16, the maximum principal stress increased by an average of 33.3% compared to the Y-direction stress value when the elements failed. e tensile stresses in the Y-direction (tangential direction) when the elements A, B, and C failed in respectively. e X-direction (radial direction) tensile stresses were 2.1 MPa, 1.9 MPa, and 1.8 MPa, respectively. Although the latter was smaller than the former, the two were in the same order of magnitude. e maximum principal stress at the tip of the crack was increased by 33.3% due to reflective stretching. erefore, according to the first strength theory, radial reflective stretching could not be ignored during the penetration of the crack between the blastholes.

Conclusion
Based on the existing stress wave theory, theoretical derivation showed that the midpoint of the blasthole connecting was not the largest stress, and the result contradicted its own assumption. e fluid-solid coupling method could be used to simulate the crack propagation of rock mass under blast loading, which could be consistent with the actual results.
is method was suitable for the discussion and research of blasting theory. Based on the double blastholes model of the dimensional parameters given in this paper, the study found that, at the midpoint of the blasthole connecting, there was no maximum value of stress regardless of the type. e article analyzed the penetration mechanism of cracks between blastholes. When the two stress waves met at the midpoint of the blasthole connecting, the cracks between the blastholes did not penetrate each other. e stress wave continued to propagate to the tip of the crack generated by the adjacent blasthole and produced a reflective stretching effect. It made the maximum principal stress value of the crack tip increase by 33.3% compared with the tangential tensile stress. erefore, the reflection of stress waves at the tip of the crack generated by the adjacent blasthole also plays an important role in the process of cracks penetration.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.