Abstract

Shaped energy blasting has been widely used in the field of geotechnical engineering because of its good orientation and high energy utilization. However, the bifurcation of cracks in the direction of energy accumulation seriously affects the precracking effect in the direction of energy accumulation. In order to study the influence of the shaped energy angle on the crack propagation and bifurcation in the direction of energy accumulation, this paper used theoretical analysis and numerical simulation to study the influence of the energy angle on the crack propagation law in the energy-concentration direction. It was found that the energy release rate in the direction of energy accumulation after blasting was the main determinant of crack propagation and bifurcation in the direction of energy accumulation, and it decreased with the increase of the shaped energy angle. When the energy release rate in the direction of energy absorption exceeded a certain critical value, the stress intensity factor K at the crack tip would be affected by the impact load more than the bifurcation toughness KB, resulting in bifurcation of the crack in the direction of the energy. The SPH method was used to simulate and analyze the energy blasting of four different shaped energy angles. The results show that as the shaped energy angle increases when the shaped energy angle is greater than or equal to 35°, the cracks in the direction of energy accumulation after blasting are bifurcated, two cracks at the crack tip. When the shaped energy angle is less than 24°, only one horizontal crack is generated in the direction of shaped energy, which is in good agreement with the theoretical analysis. The research in this paper will provide a certain research basis for the design of the blasting device and the optimization of the blasting effect.

1. Introduction

With the widespread application of coal-free column cutting and roadway technology and the improvement of underground engineering, blasting excavation requirements, the requirements for the connectivity of cracks, and the linearity of cracks in the direction of shaped energy in shaped energy blasting are also increasing. The research and application of shaped energy blasting technology originated from the First World War and began to be mainly used for weapons and equipment. With the development of blasting technology, shaped energy blasting is gradually applied to geotechnical engineering, blasting molding, and other fields. In recent years, scholars at home and abroad have proposed the technology of blasting, using energy-concentrating hoods, energy-collecting granules, and strip-shaped and point-shaped PVC energy-concentrating tubes to collect energy blasting and focusing on grooving parameters. The influence parameters and mechanism of energy device parameters on the shaped charge jet, crack propagation, and stress wave propagation have been deeply studied, and rich results have been obtained.

In the shaped energy blasting, the design parameters of the shaped energy device are one of the important factors affecting the blasting effect of the shaped energy and directly determine the effect of the directional precracking. Guo et al. [14] studied the principle of crack generation after the explosion and the energy accumulation of various parameters of shapes energy blasting. In the shaped energy blasting, the detonation mode, different charging positions, different charging coupling coefficients, blasting control hole spacing, and so on, the parameters have different effects on the generation of blasting cracks. Wu et al.’s study [5] investigated the crack propagation in materials similar to rocks and mutual influences thereon at five different notch angles (45°, 60°, 90°, 120°, and 135°) by using acoustic emission (AE) and static strain test system based on the soundless cracking technology. The notch angle was negatively correlated with the fracture degree as evinced by specimens with prefabricated notches showing the best overall cracking effect and the fragmentation of sample when the notch angle was 45°. Li et al. [6] studied the stress state near the kerf angle in the grooving blasting. It is considered that the crack growth zone formed by the slanting direction stress zone and the crack suppression zone formed by the nongrooving direction compressive zone are the same. The cracks in the grooving direction are cracked, are expanded, and finally form a fracture surface. Xiao et al. [7] and Zhang et al. [8, 9] believe that the peak of the tangential strain of the crack along the grooving direction increases gradually during the increase of the groove angle from 30° to 60° and slightly decreases between 60° and 90°. The reason is that as the kerf angle increases, the grooving becomes wider or the grooving volume increases so that more explosive gas enters the crack, which provides more energy for crack development and strains along the grooving direction. The attenuation is small. When the angle of the groove is small, there is less explosive gas entering the crack, so the tangential strain in the direction of the groove is attenuated quickly. Yang et al. [10] studied the effects of three different materials of slitting bags on the explosion impact and explosive gas in grooving. In 2003, Academician He et al. proposed a bidirectional cumulative tensile blasting technology with a point-shaped PVC energy collecting tube as a collecting device for the complex manufacturing, construction and construction of existing shaped energy devices, high cost, difficult installation, and adaptability to geotechnical engineering [11]. This technology fully considers the rock’s high compressive strength and low compressive strength. It produces tensile stress concentration in the set direction and direct fracture of the rock mass. It is mainly used for precise directional precracking of rock. The bidirectional cumulative tensile blasting combines the energy collecting device with the ordinary explosives smartly and only needs to load the ordinary explosive into the shaped energy device on the spot and then insert the collecting device into the blast hole, the process is simple, the construction is convenient, the construction speed is fast, and the construction cost is low.

In fact, the shaped energy angle of the collecting holes on the shaped energy device has a great influence on the blasting effect. However, the influence of the concentrating holes on the blasting effect is not very clear, which brings some troubles to the further optimization of the concentrating device. It also causes some problems with the widespread application of bidirectional cumulative tensile blasting technology. In shaped energy blasting, the type of explosive, rock properties, and charge structure all affect the crack propagation. In order to study the crack propagation and bifurcation law of shaped energy blasting individually, this paper used LS-DYNA software and applied the smooth particle flow method to simulate the concentrating blasting with the shaped energy angles of 24°, 35°, 45°, and 60°, respectively, and the evolution law of the crack propagation and bifurcation of the energy in the direction of shaped energy at different shaped energy angles was analyzed, providing a certain research basis for the optimization and wide application of bidirectional cumulative tensile blasting technology.

2. Bidirectional Cumulative Tensile Blasting Rupture Mechanism

2.1. Principle of Concentrated Energy Blasting

The dynamic process of bidirectional cumulative tensile blasting is mainly in the following stages: in the initial stage, after the explosive is detonated, high temperature and high-pressure gas are rapidly generated, and the shock wave formed by the explosion directly impacts the rock mass through the shaped energy hole so that the surrounding of the shaped energy hole the rock is fractured. In the second stage, the shaped energy holes on the shaped energy tube cause the compressive stress formed by the contact of the shock wave with the wall of the hole to be rapidly converted into the tensile stress along the shaped energy direction. Because the tensile strength of the rock is much smaller than the decompression strength, it affects the crack propagation in the direction of the shaped energy. In the third stage, the interaction of the air wedge action and the stress wave, with the blasting progress, the effect of the stress wave is gradually attenuated, but the explosive gas enters the crack due to the air wedge action, increasing the tensile stress of the rock mass at the front end of the fracture, and the crack is further expanded. The use of shaped energy tubes and their parameters plays a very important role in the propagation of explosive energy [1215]. Therefore, the bidirectional cumulative tension blasting is quite different from the ordinary blasting process.

The shaped energy angle is the preset angle of the shaped energy hole of the shaped energy device, as shown in Figure 1. The size of the shaped energy angle has a great influence on the formation of the initial crack. The increase of the shaped energy angle will directly lead to the increase of the area of the energy accumulation hole and thus will bear the load released by more explosives. As shown in Figure 2, the stress intensity factor change increases with the increase of the load, which will eventually increase the energy release rate.

2.2. Influence of Energy Angle on Energy Blasting

In the bidirectional cumulative tensile blasting, the shaped energy angle has an important influence on the crack propagation of rock directional blasting. The key to the bidirectional cumulative tensile blasting technology is the shaped energy tube with point-shaped energy holes. Therefore, the optimization of the shaped energy angle parameters in the energy blasting has important scientific significance for improving the effect of the blasting crack propagation (Figure 3).

The bidirectional cumulative blasting is mainly the action of stress wave and explosive gas. In the third stage of blasting, the stress wave is gradually weakened, and the explosive gas collides with the quasistatic gas pressure in the direction of the collecting hole. As a result, the explosion cracks in the direction of the shaped energy hole and the peripheral position continue to expand forward. As the crack develops, the quasistatic gas pressure in the hole will become smaller and smaller. When the crack tip strength factor is greater than the rock tensile strength, the crack will stop expanding [16]. The crack propagation of explosive gas is mainly reflected in its quasistatic effect. It is assumed that there is no quasistatic pressure of the explosive gas in the nonpolymerization direction, and a quarter model is used for analysis. As shown in Figure 4, the explosive gas is set. The quasistatic pressure on the rock mass is the load Q.

The crack generated by the two-way shaped tensile blasting is a type I crack, and the stress intensity factor at the crack tip iswhere a is the crack length, KI is the stress intensity factor at the crack tip, d1 is the inner diameter of the shaped hole, Q is the load at which the explosive gas acts on the crack tip, and r is the radius of the blast hole.

In the third stage of the shaped energy blasting, the explosion crack of the concentrating hole is expanded by the quasistatic pressure of the gas in the horizontal direction and is subjected to a circular area having an area approximately d1/2 as a radius, and it is obtained as follows:where P is the static pressure of the explosive gas at the abscissa x.

The inner radius of the collecting hole changes with the angle of the collecting angle, sowhere z is the thickness of the collecting tube and ω is the angle of convergence.

Simultaneous (1)–(3) can be obtained as follows:

It can be obtained from the above formula that, in the interval where ω belongs to 0° to 180°, the stress intensity factor KI increases with the increase of ω, and the stress intensity factor increases with the increase of the shaped energy angle. Therefore, it is possible to change the magnitude of the shaped energy angle to affect the stress intensity factor at the crack tip.

When the energy intensity blasting reaches a high level, the crack will bifurcate during the expansion process, and the defects in front of the main crack develop into microcracks. The interaction of these microcracks determines the main crack propagation path and bifurcation behavior [17].

The energy release rate is always expressed in the form of the total energy or strain energy differently, and the energy of the crack propagation is provided by the strain energy released into the structure. In the plane stress problem, when the crack is a type I crack:where G is the energy release rate, E is the elastic modulus, and KI is the stress intensity factor.

It can be seen that the energy release rate is related to KI. When KI is greater than the fracture toughness KIC of the rock, the crack begins to crack. When K = KM = 1.4 KIC, the crack will likely appear as microscopic bifurcation (minor branch cracks from the main crack). When K reaches a certain critical value KB, the crack will be microscopically bifurcated (grown from the main crack and two or more cracks are growing at almost the same speed) [18]. For high-speed extended cracks, when the energy release rate is much larger than the crack propagation resistance, the excess energy will generate new cracks, causing the crack to bifurcate. Such a phenomenon will seriously affect the effect of directional blasting of rock [19] (Figure 5).

2.3. Crack Bifurcation Angle Calculation

Based on the J-integral theory, the maximum energy release rate of the pure I-type crack bifurcation is calculated, and the crack tip is divided into four quadrants to calculate the maximum energy release component of each quadrant. Formula [18] that extends toward the maximum energy release rate is

When the energy release rate reaches the maximum, that is, when the Gb is maximum, the crack will start to propagate along the angle of α, then it can be obtained from formula (6) that

Let

One has

Substituting (8) to (9) into (7) giveswhere Gb is the energy release rate, υ is Poisson’s ratio, α is the bifurcation angle, and E is the elastic modulus.

To get the maximum value of Gb, so . We know that . Therefore, when , the energy release rate is the highest. When the crack bifurcates, its bifurcation angle is 32.48°.

3. Numerical Simulation of Concentrated Energy Blasting

3.1. Establishment of the Numerical Model

The numerical simulation software selects the full-featured nonlinear finite element numerical simulation software LS-DYNA. The algorithm is mainly the Lagrange algorithm and contains a thermal analysis function. The main analysis object is nonlinear dynamic analysis, which is very suitable for the explosive explosion, vibration and shock, and other issues. Considering that the initial ground stress is much smaller than the compressive stress formed by the detonation gas, the initial ground stress is neglected during the numerical simulation. In the numerical simulation process, the smooth particle flow method is used to carry out a two-dimensional numerical simulation of the presplitting blasting of the energy-concentrating tube. In the process of simulating material large deformation in the numerical method based on the grid, it is easy to terminate the calculation due to grid distortion, and it is difficult to simulate the problem with a strong impact discontinuity. Explosion, penetration, and high-speed collision involve large deformation, the strong impact, and the interaction process of the material interface; the smooth particle flow method makes up for these problems based on the grid computing method [20, 21]. Explosives use mineral emulsion explosives. The model of explosives uses MAT_HIGH_EXPLO_SIVE_BURN. The pressure-volume relationship of explosives during detonation is expressed by the JWL equation of state aswhere P is the pressure of detonation products, A, B, R1, R2, and omega are the material constants determined by experiments, V is the relative volume of detonation products, and E0 is the initial internal energy density of detonation products. Parameters of explosives and the JWL equation of state are shown in Table 1.

The polyenergy material is selected from the PVC pipe, and the specific parameters of the constitutive equation using the plastic follow-up model are shown in Table 2. The plastic follow-up model is a hybrid model of isotropic, creep hardening or isotropic and follow-up hardening, which is related to strain rate and can be considered for failure. Isotropic or follower hardening is selected by adjusting the hardening parameter β between 0 (synchronized hardening only) and 1 (isotropic hardening only). The strain rate is considered in the Cowper-Symonds model, and the yield stress is expressed by the factor associated with the strain rate as follows:where σ0 is the initial yield stress. C and p are copper-Symonds strain rate parameters. is the effective plastic strain. Ep is plastic hardening modulus, , Etan is the tangent modulus. Specific parameters of shaped charge materials are shown in Table 2.

The rock material in the model uses the Johnson-Holmquist II (JH2) model to study the crack propagation of rock mass blasting process under blast loading. Fine rock sandstone is selected as the reference to the rock parameters. The specific parameters are selected as shown in Table 3.

3.2. Numerical Simulations of Single Hole Blasting with Different Shaped Energy Angles

In this paper, LS-DYNA software was used to simulate the influence of different shaped energy angles on the blasting, and the length and angle of the crack were analyzed. The JH2 model used in the rock materials in the numerical model is a nonhardened material model. Its yield surface and failure surface coincide. Generally, the equivalent stress is used to express the elastic limit of the failure criterion. The yield surface expression in the model includes the normalized pressure and the maximum normalized tensile hydrostatic pressure. Therefore, in numerical simulation, the damage is the result of the combined effect of the two. In addition, the damage itself is a microscopic description of the crack, and the macrocrack is also caused by the action of a large number of microcracks. Therefore, the damage in the numerical simulation results would be used to describe macroscopic cracks.

In order to verify the influence of the shaped energy angle on the effect of the shaped energy blasting, a single blast hole with different angles of shaped energy angle was numerically simulated with reference to ordinary blasting. As shown in Figure 6(a), the blasting simulation results in the case of no shaped energy device: in the initial stage, a shear stress failure zone was formed around the hole, the crack was uniformly cracked in six directions at an angle of 60°, and Figure 6(b) can be seen (the development process of the equivalent stress in the blasting process). After forming the shear stress failure zone at 0.1 ms as shown in Figure 6(a), six cracks at equal angles were formed, and stress concentration occurs at the crack tip. As the blasting progresses, the crack tip was continuously unstable and the crack was continuously caused, expanding until the crack runs through.

As shown in Figure 7(a), the blasting simulation results with a shaped energy angle of 24° were obtained: in the initial stage, microcracks were formed near the shaped energy holes, and as the explosive particles move, the cracks along the collecting energy were driven. The direction was cracked while presenting a single horizontal crack (crack I; crack II). It could be seen from Figure 7(b) that, during the development of the equivalent stress in the blasting process, stress concentration occurred near the shaped energy hole immediately after the explosive was detonated, and the nearby rock was subjected to the tensile stress. As the blasting progresses, the crack tip was continuously destabilized, causing the crack to expand continuously until the crack penetrates, presenting a single horizontal crack pattern.

As shown in Figure 8(a), the blasting simulation results with a shaped energy angle of 35° resulted in the formation of microcracks near the shaped energy holes in the initial stage, but with the continuous movement of the blasting particles, four and horizontal cracks appear. Symmetric cracks (crack I, crack II, crack III, and crack IV) at an angle of 30° are shown in Figure 8(a) and begin to dominate the cracking direction of the crack. Figure 8(b) shows the evolution of the equivalent stress during the blasting process. As shown in Figure 8(b), immediately after the detonation of the explosive, stress concentration occurs near the shaped energy hole, but four stress concentration regions symmetric with the direction of shaped energy appear immediately, and the nearby rock is subjected to tensile stress. It can be seen from the figure that the stress in the direction of the shaped energy hole is negative. As the blasting progresses, the crack tip is continuously destabilized, causing the crack to expand continuously until the crack penetrates, and the equivalent stress tends to zero, forming a symmetric crack.

Figures 9(a) and 10(a) show the results of blasting damage development at the shaped energy angles of 45° and 60°. In the initial stage, microcracks are formed near the shaped energy holes, but with the continuous movement of the explosive particles, four symmetric cracks (crack I, crack II, crack III, and crack IV) appearing at different angles to the horizontal crack appear, and with time, the symmetric crack begins to dominate the cracking direction of the crack. Figures 9(b) and 10(b) show the evolution of the equivalent stress during the blasting process. Immediately after the detonation of the explosive, stress concentration occurred near the shaped energy hole, but four stress concentration areas symmetrically with the direction of the energy accumulation appeared immediately, and the nearby rock was subjected to tensile stress, which can be seen from the figure. The hole direction stress is negative. As the blasting progresses, the crack tip is continuously destabilized, causing the crack to expand continuously until the crack penetrates, and the equivalent stress tends to zero, forming a symmetric shear crack.

3.3. Crack Propagation Process of Different Shaped Energy Angles

The relationship between the crack propagation angle and the crack length during the crack propagation process of the 24° shaped energy angle blasting can be obtained from Figure 11(a). It can be seen from the figure that two cracks are 180° apart. In Figure 11(b), in the middle, the peak lengths of the two cracks are shown. After 0.1 ms, the crack tip cracking directions are constant at 0° and 180°, respectively. In Figure 11(c), the length of crack propagation can be seen as a function of time. The initial expansion speed of the crack is the same, and the slope is consistent, but at the end of the latter, the crack propagation speed is greater than the speed of the right crack. The shaped energy device can be obtained with a satisfactory blasting effect.

The relationship between the crack propagation angle and the crack length during the crack propagation process of the 35° shaped energy angle blasting can be obtained from Figure 12(a). From the figure, we can see the four crack length peaks, corresponding to four cracks. In Figure 12(b), the peak lengths of the four cracks are shown. After 0.1 ms, the crack tip cracking directions are constant at 30°, 170°, 190°, and 350°, respectively. And 35°, in Figure 12(c), it can be seen that the length of crack propagation with time and the initial expansion speed of the four cracks are basically consistent.

From Figures 13(a) and 14(a), the relationship between crack propagation angle and crack length during the crack propagation process of 45° and 60° energy angle blasting is obtained. The crack length also increases with time. Four packs of crack length are obtained from the figure, indicating that four cracks appear at four angles, and four crack length peaks are shown in Figures 13(b) and 14(b). The crack tip cracking direction is 0.1 ms later. The constant lengths of crack propagation over time can be seen in Figures 13(c) and 14(c) at constant temperatures of 30°, 160°, 200°, and 340°, respectively, and the initial expansion speeds of the four cracks are substantially the same.

Numerical simulation results showed that a single horizontal crack (crack I; crack II) appeared in the single-hole blasting with a shaped energy angle of 24°, which expanded in the direction of energy accumulation, and the two cracks were 180°. At the single hole blasting with shaped energy angles of 35°, 45°, and 60°, four cracks (crack I, crack II, crack III, and crack IV) at an angle of about 30 degrees from the horizontal appeared. The symmetry crack began to dominate the cracking direction of the crack. In the theoretical analysis, the increase of the shaped energy angle is the increase of the stress intensity factor and the energy release rate at the crack tip. When the energy release rate is far greater than the crack propagation resistance, excess energy will generate new cracks, causing cracks to bifurcate. In the theoretical analysis, the J-integral theory is used to calculate the 32.48o bifurcation angle when the crack is bifurcated.

It could be known from Table 4 and Figure 15 that when the shaped energy angle was 24°, there was no bifurcation after the blasting. When the shaped energy angle was 35°, after the blasting, the main crack was one-side bifurcation. From the main crack, two bifurcation cracks which were 38° and 40° from the main crack are forked. When the shaped energy angle was 45°, after the blasting, the main crack was both sides bifurcation. From the main crack, four bifurcation cracks with 31°, 31°, 31°, and 30° were branched from the main crack. When the shaped energy angle was 60°, after the blasting, the main crack was both sides bifurcation. From the main crack, four bifurcation cracks with 26°, 33°, 29°, and 27° are branched from the main crack. The above angle values were measured from the shaped energy blast damage development map and were close to the theoretical value of 32.48°. Therefore, the numerical simulation agrees with the theoretical analysis conclusion.

4. Field Experiment

4.1. Project Overview

The 1105 working face of Hecaogou No. 2 Coal Mine adopts the inclined longwall coal mining method. The working face is mainly a near-level, gentle monoclinic stratum, and the coal seam inclination angle is 1° to 3°, with an average of 2°. The length of the 1105 working face is 120 m, and the trend length is 1140 m, which is shown in Figure 16. The main coal seam thickness of this working face is 0.83∼0.85 m, the average thickness is 0.84 m, and the mining height is 1.2 m. The average coal seam depth is 63 m. The test roadway selected the 1105 working face to return to the wind, and the roadway width × height = 3.5 m × 2.4 m. Roadway support design is shown in Figure 17.

4.2. Scheme Design

Combined with the results of numerical simulations, when the shaped energy angle is greater than or equal to 35 degrees, the crack would produce a bifurcation phenomenon. In order to ensure the linearity of the crack, the minimum shaped energy angle was designed to be 0°. The field experiment used the method of interval blasting. The distance between the gun holes was 500 mm, the depth of the gun holes was 3.5 m, and the angle between the holes and the plumb line was 20°. During the experiment, one shaped tube was used for each blast hole, and three rolls of mining emulsion explosives were installed in the middle of the shaped tube. The charging structure diagram is shown in Figure 18. The diameter of the hole was 48 mm, the diameter of the shaped tube was 36 mm, the diameter of the drug roll was 32 mm, and the uncoupling coefficient of the charging structure was 1.5.

4.3. Experimental Results

Since the cracks formed between the blast holes cannot be visually observed after the blasting, the CXK6 mine intrinsically safe drilling imager was used to peek at the site, as shown in Figure 19. It can be seen that after the explosives exploded, obvious continuous cracks were generated from the interval holes to the bottom of the holes, and the cracks extended along the line of the gun holes, indicating that continuous cut surfaces can be formed between the holes. It is concluded that the linearity of crack propagation can be guaranteed when the energy accumulation angle is 0°.

5. Conclusion

This paper used LS-DYNA software to adopt smooth particle flow method to simulate the bidirectional cumulative tensile blasting when the shaped energy angles are 24°, 35°, 45°, and 60°, respectively, to explore the blasting crack growth, bifurcation and equivalent stress development law in the direction of energy accumulation during blasting which were studied in detail, and the effect of energy accumulation devices with different energy angles on the directional blasting effect of rocks.(1)According to the theoretical analysis, during the blasting, the increase of the energy-gathering angle of the energy-gathering device will increase the explosive gas bursting into the fracture to a certain extent, thereby releasing more energy when the fracture expands. The stress intensity factor of the crack tip will increase with the increase of the energy concentration angle, and the increasing trend will become more and more significant. It will also release more energy in the energy concentration direction. For high-speed extended cracks, when the energy release rate is much larger than the crack propagation resistance, the excess energy will generate new cracks, causing the crack to bifurcate. The crack and main crack will be 32.48° according to the maximum energy release theory.(2)Using the LS-DYNA software and the smooth particle flow method, a numerical model with an energy accumulation angle of 24°/35°/45°/60° was established. The upper crack is bifurcated. When the energy accumulation angle is less than 35°, that is, when the energy accumulation angle is 24°, the crack in the energy accumulation direction is ideal, and no bifurcation occurs. It is known from the field experiments that the linearity of crack propagation can be guaranteed when the energy accumulation angle is 0°. The numerical results, field experiments, and theoretical analysis results agree well.(3)According to the equivalent stress evolution law of the four models in the numerical simulation process, it can be known that as the shaped energy angle increases, the stress near the energy accumulation angle continues to increase, thereby releasing more energy in the energy accumulation direction. When the shaped energy angle is equal to or greater than 35°, the results show that as the shaped energy angle increases when the shaped energy angle is greater than or equal to 35°, the cracks in the direction of energy accumulation after blasting are bifurcated, two cracks at the crack tip, and with time, the symmetrical cracks begin to dominate the crack’s cracking direction.

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 (Grants nos. 51904188 and 41602308), the Open Fund of State Key Laboratory for GeoMechanics and Deep Underground Engineering (Grant no. SKLGDUEK1821), and Research and Development Project of Guizhou University of Engineering Science (Grant no. G2018016). In addition, thanks are due to all the people who contributed to the research and the manuscript.