Abstract

Conventional smoothed particle hydrodynamics (SPH) methods suffer from disadvantages, such as difficult initial particle configuration, uneven distribution of generated particles, and low computational efficiency when applied to numerical simulation of shaped charge blasting. In this research, to overcome these problems, a modified SPH method that generates the particle configuration through self-adaptive optimization is developed by the combined application of MATLAB and LS-DYNA. The results presented in this paper demonstrate that the modified configuration method solves the problem of uneven distribution of particles in complex geometry domains by providing a more uniform smoothed particle distribution than the conventional SPH method. Furthermore, the results from the application of these two methods to the bidirectional-shaped charge blasting problem reveal that the defects in the particle configuration in the conventional SPH method lead to the development of main cracks in both the shaped and the unshaped directions. However, with the self-adaptive optimization method, the main cracks develop only in the shaped direction. In addition, the equivalent stress difference between the shaped and unshaped directions, 0.7 ms after detonation, is 120 MPa with the modified method. This is 85 MPa more than that with the conventional method.

1. Introduction

A comprehensive study on the principle of cumulative blasting and its effects began during World War II. Before the 1980s, the researchers in this field focused on the process of jet penetration mechanism and modeling of cumulative blasting [1, 2]. Since the 1980s, owing to the development of experimental methods and computer technologies, there have been abundant research studies on the subject of shaped charges and their energy release process [3, 4]. Studies on the cumulative blasting technology that enhances the blasting effect through changing the propagation medium of detonation wave and installing accessories on explosives, such as blasting by water-sealed explosive charge packages and notched blasting, have emerged subsequently [5, 6].

Due to the advances made in engineering software development, many researchers, nowadays, use commercial software to simulate the problems related to shaped charges. Kim et al. [7] and Donze et al. [8] used the discrete element method to simulate the behavior of rocks during blasting and studied the propagation behavior of stress waves in rock mass under different loading conditions. Yang et al. [9] and Fan and Li [10] used DYNA3D to study the presplitting with a decoupling charge of a binding energy tube. Their results showed that the shaped tube had a significant effect on the direction of energy propagation and protected the side rock wall. Ning et al. [11] used the discontinuous deformation analysis (DDA) method to simulate the motion behavior of a rock during the explosion process. Mussa et al. [12] used the arbitrary Lagrangian–Eulerian (ALE) technology in LS-DYNA to simulate the surface explosion of underground tunnels. Their simulation results proved the reliability of the peak particle velocity method and the single degree of freedom method for evaluating the damage criteria and degree of damage of box-shaped underground tunnels.

Due to the advantages offered by the SPH method over the conventional numerical methods in dealing with free surface and impact problems, more and more researchers are using the SPH method to analyze dynamic problems such as blasting [1316]. However, the accuracy of the SPH method depends on whether the particle distribution in the support domain is uniform. Although there are many ways to ensure the uniformity of particle distribution in large or simple geometric domains, for example [17, 18], the particles tend to distribute uniformly during the SPH simulation process if appropriate kernel function and smoothing length are adopted, but the traditional SPH method, typically, faces some difficulty in ensuring a uniform distribution of particles for models with complex geometries or those with a very small geometric domain. The existing repulsive force method [19] and the virtual particle method [20] primarily solve the problem of the boundary particles missing which happens due to these particles being “truncated” by the boundary in the support domain. The repulsive force method can easily form the initial disturbance; however, the repulsion force is an artificial force, and hence, the interaction between the real solid particles and the internal particles cannot be described. The virtual particle method, on the other hand, has difficulties in dealing with complex boundaries. The mass of the boundary particle in the unified semianalytical wall (USAW) boundary conditions method proposed by Ferrand et al. [21] depends on the angle between the boundary and the particle and, hence, cannot accurately simulate complex geometric domains. Moreover, if the particle distribution is not uniform, the higher-order continuity condition cannot be accurately satisfied, even if the support domain of the inner particles has not truncated the constants of the discrete form and the linear continuity conditions.

In summary, the uniformity of particle distribution has an important influence on the crack propagation of energy blasting. And different methods can be used to correct the continuity conditions of the discretization equations, but the positions of the particles are changing continuously during a blasting process [10, 22]. Considering the cost and efficiency, it is difficult to construct a smoothing function for the particles in each time step of the simulation. In this research, we have developed a MATLAB program to adaptively optimize the uniform generation of SPH particles. And when MATLAB completed the modeling process, LS-DYNA continued to complete the calculation and postprocessing. The computational efficiency is greatly improved, and the particle distribution is more uniform than that achieved with the conventional SPH method. In essence, the new configuration method solves the problem of SPH particle formation in complex geometric domains and the boundary forces calculation in the numerical simulation of shaped charge blasting.

2. Automatic Adaptive Optimization for SPH Particle Generation

2.1. Automatic Adaptive Optimization Model

The automatic adaptive optimization for particle generation is also based on the mesh generation methods, and the SPH particles are arranged at the positions of the conventional mesh nodes. The main difference between the SPH method and a conventional mesh generation method is that a very uniform mesh is generated in the former method. It is worth noting that the entire particle optimization process had been completed before the numerical model calculation was started.

For two-dimensional problems, the meshing is usually based on the Delaunay triangulation method. In general, any point set in the X-Y plane can be divided into a mesh of Delaunay triangles, assuming that each side of these triangular meshes is a member bar and each member bar corresponds to two endpoints. Assuming that the endpoint of every member bar can be moved, the movement function is , where is the expected length of the member bar and l is the current length of the member bar.

Now, let us suppose that the movement function is applied to each member bar in the form of a force. Each point beyond the calculated boundary point will be returned to the boundary by the external force applied perpendicular to the boundary. This force is just capable of preventing the point from escaping the boundary. The endpoints of all the member bars ultimately pass through the static balance equation of the bar system. The final expected result is that the length of each member bar reaches the same length, thereby achieving a more uniform meshing within the region and providing a carrier for the SPH particle points.

That is to say, firstly, the mesh was generated in the calculation domain by the triangulation method. And then, the movement function was used to adjust the length l between the adjacent nodes of the grid to l0. But, during the adjusting process, when the node jumps out of the calculation domain, a force that was perpendicular to the boundary of the calculation domain was applied to move the node to the boundary of the calculation domain. When the length between nodes is l0, a uniform mesh is obtained. Using the nodes of the uniform mesh as the carrier to arrange the particles, a uniformly distributed particle arrangement was obtained. The schematic diagram of the optimization process is shown in Figure 1.

For the mechanical equilibrium process that needs to be solved, the x and y directional components of the mesh nodes can be written as vector , where

Then, the force applied to each node contains both the x and y components as shown below:where contains the internal force term, which is the force exerted on each node inside the calculation area; is the external force applied to the node when the node is outside the calculation area, to pull it back into the calculation area. The first term of F contains the forces in the x and y directions.

It should be pointed out that F is not a real force, but rather a virtual force, whose value only depends on the geometric properties of all the bars in the calculation domain. In this method, the initial bar structure is also obtained by Delaunay triangular meshing and the mesh nodes are all in the domain. So initially, is equal to 0. Moreover, is not a continuous function of position , because with each Delaunay triangulation the node moves to a new position.

To make , the resultant force on each triangle node should be 0; in other words, each internal angle in the triangle reaches 60°, thereby achieving an equilateral triangle.

However, solving the function is very difficult, because not only is the node position moving but also the system is subjected to an external reaction force when the node moves out of the domain. To facilitate a solution for this function, artificial time-derivative terms are introduced. For , the partial differential equation of the system is assumed to be as follows:

If equation (3) has a steady-state solution, then there must be a such that . Equation (3) can be iteratively solved using the Euler method. Using time discretization, where , the approximate solution can be achieved by solving the following equation iteratively:

At the start of the iteration, , as well as the end position of each bar, is known. The external force will act when the node runs out of the computational domain and will pull back into the computational domain to the location .

There are many forms of the force function ; usually, the linear function is most commonly chosen:

The above bilinear force function form can yield satisfactory results. In general, the value of k is taken as l. It can be seen that when , where is a constant for all cells with uniform meshes.

2.2. SPH Particle Generation for Cumulative Blasting

As a first step, a set of randomly distributed points or equispaced mesh nodes based on a mesh generation method is generated in the solution domain. These nodes are stored in a matrix of size N rows and 2 columns. All the points falling outside the calculation domain are removed so that every point in the vector is in the calculation area. Then, as the solution of the system of equations progresses iteratively, the position vector is continuously updated. At the beginning of each iteration, Delaunay triangulation is performed on all the points in the calculation domain to determine the specific value of the force function. Then, the three sides of all the Delaunay triangles are composed into a matrix of three columns, with each row representing the force function value of the one vertex of the triangle. If a point escapes from the computational domain after the iteration, an external force is applied to return the point to the point closest to it in the computational domain. Note that the points can move on the boundary. The iterations are stopped when the target is reached. The adaptive optimization process after the initial particle allocation of the rock hole is shown in Figure 2.

In the directional blasting application, the hole spacing is generally about 500 mm. The diameter of the blast hole is usually 48 mm. The thickness of the shaped tube is 8 mm, and the diameter of the explosive package is only 24 mm. In the implementation of the SPH method for this problem, there should be an adequate number of SPH particles in the regions of the explosive package and the shaped tube, and the spacing between the particles should be uniform. The particle configuration diagram can be obtained as required, by using the adaptive optimization particle generation method.

The size of the hole of cumulative energy cannot be increased without limitation due to the limited thickness of the shaped tube. Therefore, as shown in Figure 3, the particle configuration result obtained by the modified method proposed in this paper is effective, and it has an adequate number and uniformity of particles in the computational domain.

3. Mechanical Model of the Bidirectional Cumulative Tensile Blasting

Bidirectional cumulative tensile blasting technology is a new directional fracturing control technique proposed by He et al. [23]. This technique takes advantage of the fact that the tensile strength of rock is much lower than its compressive strength. Here, the cumulative energy device is combined with a common mine explosive. A concentrated tensile force, larger than the rock’s tensile strength, is generated in two set directions so that the tensile failure of the rock mass occurs in these two set directions. Owing to its convenience of operation and the simple process, this technology has been successfully used in projects such as roadway retaining by roof cutting and chamber blasting.

The principle of bidirectional cumulative tensile blasting is shown in Figure 4. In this technology, the bidirectional cumulative device has an instantaneous inhibitory effect on the detonation product after the explosive detonation, and the two cumulative holes which are linearly distributed preferentially provide an instantaneous pressure relief space [22]. A high-energy jet is formed that is concentrated on the corresponding wall of each hole to create a radial initial crack. The high temperature, the pressure, and the gas velocity generated by the detonation are released from the cumulative hole, further driving the radial initial crack propagation. As the stress wave propagates in the radial direction, it generates a reverse tensile stress concentration in the vertical setting direction and also accelerates crack propagation along shaped energy direction. Thus, the rock mass is stretched and cracked along the energy accumulation direction. At the same time, due to the inhibition and buffering effect of the tube wall and the preferential unloading of the detonation products from the energy accumulation hole, the effect of stress in the direction of noncumulative energy decreases sharply. This, in turn, restricts the damage caused by the detonation products to the hole wall and inhibits the development of nondirectional cracks. When several such bidirectional cumulative blasting holes are simultaneously detonated, with the propagation of stress waves, a superimposed stress field is generated between the blast holes and the tensile stress between the blast holes increases. If the spacing of the blast holes is appropriate, the cracks in the shaped energy direction of the adjacent-blast holes can be cut through to form a smooth precracked surface for realizing the directional and accurate control of the blasting.

According to fracture mechanics theory [24], when the stress intensity factor at the crack tip is greater than the fracture toughness (KIC) of the rock, the crack initiation begins. Otherwise, the crack stops expanding. Considering the symmetrical crack of the hole wall of a bidirectional cumulative tensile blasting as an example (see Figure 5), when the shaped energy jet contacts the initial crack tip, the stress intensity factor at the crack tip is given bywhere is the initial cumulative crack length; is the explosive radius; is the pressure when the explosive particles fill the blast hole; and is the stress intensity factor correction coefficient, which is a function of the blast hole radius r and crack length a:

The function f is defined such that when , the value of changes greatly, and when , the value of changes little and approaches 1.

Furthermore, from the theory of fracture mechanics, the condition for crack initiation and crack propagation can be stated through equation (8) as follows:

When a crack initiation takes place, the “air wedge effect” of the subsequent explosive gas drives the crack to expand further, and the crack propagation causes the pressure of the explosive gas to drop. To ensure the continuous extension of the crack, the instantaneous pressure of the explosive gas must meet the condition stated by equation (9) as follows:

It can be seen from equations (3) and (4) that the tangential stress generated by the penetration of the shaped energy jet has the beneficial effect of reducing the pressure required for crack initiation and expansion in the direction of energy accumulation. Furthermore, the tangential stress is developed in such a way as to produce directional tensile cracks, which play an important role in effectively controlling the direction of crack propagation. In addition, the explosive gas preferentially enters the fissures in the direction of energy accumulation, which not only promotes the expansion of cracks in this direction but also controls the initiation and development of cracks in the noncumulative direction. Thus, the destruction of the rock mass in the noncumulative direction is effectively avoided.

4. Numerical Simulation of Cumulative Blasting

4.1. Material Model of Rock

The dynamic analysis software LS-DYNA was used to carry out a two-dimensional numerical simulation of the presplitting blasting of the cumulative tube. The rock selected in the numerical model was of granite. The Johnson–Holmquist-2 (JH-2) model was used to study the stress and crack propagation characteristics of the granite rock during the blasting process. The JH-2 strength model represents the equivalent stress of the material as a power function of the hydrostatic pressure and is related to the strain rate and the damage factor D. The dimensional-strength model is defined as

When the material is not damaged (i.e., D=0), the dimensional-equivalent stress can be expressed asAnd when the material is completely in failure (i.e., D = 1), the dimensional-equivalent stress is

In equations (10)–(12), is the strength of complete material and is the strength of completely broken material; A, B, C, M, and N are the material constants introduced in the material model; is the dimensional-hydrostatic pressure defined as , and is the pressure component of the material at the elastic limit of Hugoniot; and , where T is the maximum static tensile stress that the material can bear. is the ratio of real strain rate, and is the reference strain rate.

In the JH-2 model, plastic deformation depends on the stress level of the material. The equivalent plastic deformation resulting in failure is

In equation (13), and are material constants. With the accumulation of plastic deformation of the material, the total damage, D, accumulates and is expressed as

In equation (14), is the equivalent plastic deformation increment generated between two calculation steps.

Poisson’s ratio, elastic modulus, and yield strength of the granite needed for the simulations were obtained from uniaxial compression test data. The uniaxial tensile strength was obtained by the Brazilian splitting test. Velocities of the waves which is pressure wave and S which is shear wave were measured by wave velocity measuring equipment. The shear modulus was calculated from the results of the rock shear test. Table 1 shows all the static mechanical parameters of granite.

The material constants in the granite JH-2 model were calculated by the Hoek–Brown Failure Criterion, which can be stated as follows:

Equation (15) is the relationship between the granite samples obtained by Mahdi and Banadaki [25]. For normalized full strength and static pressure, and are the necessary parameters. can be obtained according to the method of Johnson and Holmquist [26, 27]:

Here, the Hugoniot elastic limit () for Kuru granite is 4.5 GPa.; represents the volumetric strain at , which is equal to 0.057 in this expression; represents the shear modulus; is assumed as 2.67 GPa; and is assumed as 2.74 GPa:

4.2. Explosive and Shaped Tube Material Model

In numerical simulations, the parameters of explosives were the same as those of emulsion explosives used in practice. The model used for the explosive material in LS-DYNA was MAT_HIGH_EXPLO_SIVE_BURN. The JWL EOS state equation with energy form is

In equation (18), is the pressure in the hole; , and are the explosion parameters; and Vr is the relative volume: V/V0. V is the specific volume of the detonation product and V0 is the initial specific volume of explosive before blasting: V0= 1/ρ0, and ρ0 is the density of explosive. E is the specific internal energy of detonation products. The specific parameters are shown in Table 2.

The cumulative material used in the shaped charge blasting simulations is a PVC pipe, and the constitutive equation is the plastic follow-up model in this paper. The parameters of this model are shown in Table 3. The plastic follow-up model is a hybrid model of isotropic-follow-up hardening or isotropic-follower hardening model, which is related to the strain rate, and can be considered for failure. Isotropic or follower hardening is selected by adjusting the hardening parameter between 0 (follow-up hardening only) and 1 (isotropic hardening only). For the strain rate, the Cowper–Symonds model is used, and the yield stress is related to the strain rate as follows:

In equation (19), is the yield stress of the shaped tube; is the strain rate of the shaped tube; and are the strain rate parameters; is the initial yield stress of the shaped tube; is the hardening parameter; is plastic hardening modulus, ; and finally, is the effective plastic strain.

5. Results and Discussion

Figures 6(a) and 6(b) show the results from the simulation of crack propagation using the conventional SPH and the automatic adaptive optimization SPH methods, respectively. We see from Figure 6(a), that when the conventional SPH method is used to simulate the blasting process, the crack propagation progresses not only in the direction of energy accumulation but in some other directions as well. This is due to the inability of the method to adaptively optimize the particle distribution during the particle motion. This result is quite different from the actual cumulative energy blasting. In Figure 6(b), which shows the result from the automatic adaptive optimization SPH method, the individual particles are moved through automatic adaptive optimization in line with the actual energy blasting process. In the initial stage after detonation, microcracks are formed near the shaped holes. With the continuous movement of the explosive particles, the cracks are driven to expand in the shaped energy direction, while the interfaces of the other direction are well preserved. The “air wedge action” generated by the particles entering the crack further accelerates the expansion of the horizontal crack and, eventually, presents a single horizontal crack along the shaped energy direction.

Figure 7(a) and 7(b) show a quantitative comparison of the simulation results for the crack propagation process by the conventional method and the automatic adaptive optimization method, respectively. When the conventional method is used, two cracks (Crack 1 and Crack 2) are generated in the shaped energy direction and other direction (Crack 3 and Crack 4) with a crack angle of 60° (see Figure 7(a)). The other-direction cracks are all in the upper part of the model, and the crack propagation speed in the other direction is slower than that in the shaped energy direction. The cracks are bifurcated in the shaped energy direction showing a plurality of Y-shaped cracks with different sizes (Crack 3 and Crack 4).

In Figure 7(b), when the automatic adaptive optimization method is used, two cracks are generated only in the direction of energy accumulation after detonation. As the explosive particles move continuously, the damage cracks also expand continuously, and finally, a smooth single horizontal crack appears along the direction of energy accumulation (Crack 5 and Crack 6), consistent with the actual blasting process. Furthermore, the crack propagation speed along the shaped energy direction is almost the same as that in the conventional method.

To prove that the simulation results of the optimized particle distribution model were more in line with the actual results than the unoptimized model, we conducted a site test in the gateway roof. Because the gateway roof is silty mudstone, compared to the crack propagation length, the silty mudstone is considered to be composed of uniformly distributed particles. After blasting, cracks occurred in both the orifice and the hole along the direction of energy accumulation. No cracks occurred in the orifice and the hole wall in the unshaped direction, which was in accordance with the numerical simulation results using the optimized model. The correctness of numerical simulation using the adaptive optimization model is verified, as shown in Figure 7(c).

At the time of the explosive detonation, the stress concentration appears in the rock mass in the shaped energy direction. Due to the action of the shaped tube, the rock mass in the shaped energy direction is subjected to tensile stress. In Figure 8(a), when the particle configuration is carried out by the conventional method, stress concentration occurs in both the direction of energy accumulation and the other direction, resulting in the occurrence of main blasting cracks in both these directions. The integrity of the rock mass in the other direction is destroyed, which is quite different from the actual blasting. In Figure 8(b), when the adaptive optimization method is used for particle configuration, two main cracks are generated only in the direction of energy accumulation after blasting, and no cracks are generated in the other direction. Thus, the numerical simulation results are more consistent with the actual blasting effect.

It can be seen from Figure 9(a) and 9(b) that when the particle configuration was carried out by the conventional SPH method, the maximum and the average equivalent stress in the direction of energy accumulation after the blasting were 320 MPa and 140 MPa, respectively. The same quantities, in the other direction, were 195 MPa and 105 MPa. Similarly, from Figure 9(c) when the adaptive optimization method for particle configuration was used, the maximum and the average equivalent stresses in the direction of energy accumulation were 380 MPa and 160 MPa, respectively, while in the other direction, the same quantities were 350 MPa and 40 MPa. Thus, in the modified method, the difference in the average equivalent stress between the shaped energy direction and other direction was 120 MPa, while this difference was only 35 MPa with the conventional method. Clearly, the numerical simulation results and the stress wave transfer process in the automatic adaptive optimization method were more in line with the actual energy blasting process.

As can be seen from the above analysis, the automatic adaptive optimization SPH method for the energy blasting numerical simulation provides better results than the conventional SPH method. The particle configuration, crack propagation process, and stress wave transfer process with the modified method were more in line with the actual energy blasting process.

6. Conclusions

The initial particle configuration in the SPH method has a great influence on the numerical simulation results and is a major drawback for applying the SPH method for the bidirectional-shaped charge blasting simulation. In this research, a new method for the particle configuration is presented which accomplishes automatic adaptive optimization of the SPH particles for the blasting simulation. This is achieved through the secondary development of a MATLAB program. The main conclusions are as follows:(1)The modified particle configuration method can provide a more uniform distribution of the initial particle field and enhance the computational efficiency of the SPH method. Thus, it helps overcome the shortcomings of the conventional method by solving a key issue associated with the generation of SPH particles.(2)Using the aforementioned method, a numerical simulation model for bidirectional cumulative tensile blasting was established. A comparison of results for this application using the conventional method and the modified configuration method was carried out to assess the modified method.(3)Due to the shortcomings of the conventional particle configuration method, the main cracks were generated in both the direction of energy accumulation and other directions in the simulation process. The automatic adaptive optimization method, on the other hand, primarily produced only the main crack in the shaped energy direction. Furthermore, the automatic adaptive optimization method yielded an equivalent stress difference of 120 MPa between the direction of energy accumulation and the other direction within 0.7 ms after detonation. This was 85 MPa larger than the corresponding difference when the conventional method was used. Thus, the results from the automatic adaptive optimization method are in line with the theory and engineering practice of cumulative blasting.

Data Availability

All data included in this study are available upon request from the corresponding author.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant nos. 51904188, 41602308, and 41702381) and the Open Fund of State Key Laboratory for GeoMechanics and Deep Underground Engineering (SKLGDUEK1821). In addition, thanks to all the people who contributed to the research and the manuscript. Finally, the authors would like to express their heartfelt thanks to the language editors from Elsevier for editing the manuscript language.