A Study on the Impact of Different Delay Times on Rock Mass Throwing and Movement Characteristics Based on the FEM–SPH Method

Abstract

Burst morphology is a crucial indicator for evaluating the effectiveness of blasting, as it directly reflects the actual state of the blasting results. The results of rock displacement following blasting partially reflect the effectiveness of throw blasting, while the rock ejection process serves as the macroscopic manifestation of the blasting method. To accurately assess the impact of different delay times on burst formation, this study addressed the issues of rock movement and ejection in underground blasting. Using three-dimensional modeling, we constructed a FEM–SPH model and utilized LS-DYNA numerical simulation software to investigate the movement patterns of rock in precise delayed blasting scenarios underground. This study explored the spatiotemporal evolution characteristics of rock movement post-blasting. Digital electronic detonators were used to set precise inter-row delay times of 25 ms, 50 ms, and 75 ms. The results revealed that the ejection distance of blasted rock in underground mining increased with longer inter-row delay times, while the slope angle of the blasted muck pile decreased as the delay time increased. Furthermore, at a micro level, the study found that a 75 ms delay created new free surfaces, providing effective compensation space for subsequent blasts, thereby improving blasting outcomes. Analysis of the 25 ms and 50 ms delay periods indicated a clamping effect on rock movement. Field comparisons of blasting results were conducted to validate the influence of precise delay times on the movement patterns and spatiotemporal evolution characteristics of blasted rock.

1. Introduction

In the study of rock movement during blasting, throwing techniques are primarily applied in high-bench blasting in open-pit mining. This technique involves directly throwing part of the rock to a more ideal area or goaf through blasting. The earliest application of high-bench blasting and throwing techniques dates back to the 1960s, when Chironis et al. [1] discovered that approximately 40% of the blasted rock was directly thrown into the goaf, sparking interest and further research by scientists worldwide. Subsequent research and experimental results from various countries have demonstrated that this technique can significantly enhance the economic benefits of mining operations. However, the process of rock blasting and throwing can also lead to issues such as equipment damage, scattering of blasted rock piles, and increased costs. Asri et al. [2] studied the application of cast blasting in Moroccan phosphate mines. Zhang et al. [3] determined the key factors influencing the effectiveness of casting blasting based on gray relational analysis theory, subsequently optimizing the blasting parameters for open-pit mining. Relevant studies indicate that the impact of shock waves on the borehole wall generates radial and circumferential cracks, and the explosive gases further propagate these cracks, leading to rock fragmentation. The residual energy of the explosive gases propels the fragmented rock in the direction of least resistance. Li et al. [4], through studying the impact of explosive charge ratio on high-bench casting blasting, found that the casting distance increases with the average explosive charge ratio. Zhang et al. [5] employed a locally weighted linear regression method to determine the formation of the muck pile resulting from casting blasting. Huang et al. [6] established calculation methods for the projection velocity and rock-breaking projection energy under two blasting methods: spherical charges and cylindrical charges. By these methods, predicted velocity characteristics of fractured rocks can be obtained.

Most researchers have focused on rock throwing in open-pit mines, aiming to improve the efficiency of shoveling, loading, and transportation by controlling the rock throwing distance to achieve the desired blast pile morphology. However, few have examined the impact of rock movement on fragmentation and blasting effects in underground settings, where such movement is rarely studied. In underground blasting, rock movement is constrained by tunnel layouts, making it difficult to disperse blast piles, which is often overlooked in the practical production stages. Nonetheless, in the confined spaces of underground operations, blast pile morphology can significantly affect the efficiency of mechanical shoveling, loading, and transportation. Therefore, this paper aims to investigate the influence of rock movement patterns on fragmentation and throwing in underground blasting.

Given that the blasting process occurs instantaneously, existing technologies struggle to capture the complete process. International scholars were the first to conduct blasting research using numerical simulations. Early studies in this area were pioneered by these researchers. For continuous medium analysis methods, as early as the 1970s, R. L. Taylor and J. R. Rice used the finite element method (FEM) to analyze crack propagation and stress wave transmission issues. This method has since been widely applied to numerical simulations of rock blasting. For instance, Ma et al. [7] simulated the rock fragmentation process under different inter-hole delay conditions using the FEM. Zhou et al. [8] applied the FEM to study the explosive effects of coupling interactions in different charge structures. Additionally, the Finite Difference Method (FDM) has been used in blasting engineering, primarily to assess the impact of blasting activities on other targets [9]. However, continuous medium analysis methods are limited in simulating rock mass movement and the subsequent fragmentation process. To simulate the entire rock fragmentation process, discontinuous medium analysis methods are needed. In this context, Zhao et al. [10] employed the Discontinuous Deformation Analysis (DDA) program to simulate tunnel excavation blasting. Han et al. [11] used a coupled FEM and discrete element method (DEM) model for bench blasting in deep tunnels. Pan et al. [12] combining finite and discrete element methods, investigated the effects of mineral particle composition and size on the macroscopic blasting pattern of rocks by studying the micromechanical properties of rock mineral particles. Jayasinghe et al. [13] developed a coupled model of Smoothed Particle Hydrodynamics (SPH) and the FEM to study damage and crack patterns in rock blasting. Moreover, the Universal Distinct Element Code (UDEC) program [14] and the Particle Flow Code (PFC) [15] are also capable of simulating the rock blasting process.

Smoothed Particle Hydrodynamics (SPH) is a numerical method that integrates mesh-free techniques, Lagrangian methods, and particle flow methods. This approach uses particle-like elements to carry mechanical quantities, tracking their motion and accurately obtaining corresponding mechanical information through specific mathematical relationships. Due to algorithmic differences, the SPH method requires strict computational conditions and extended time compared with traditional finite element methods. It also has shortcomings in addressing boundary conditions. Therefore, to comprehensively and accurately simulate the blasting impact process, an SPH–FEM method is proposed, which employs the SPH method in the near field and the finite element method in the far field. In the numerical simulation of rock blasting processes, the SPH–FEM coupling method combines the advantages of both approaches. SPH is used to simulate rock fragmentation and large deformations, while the conventional FEM ensures stability and efficiency in the far field of the blast. At the coupling interface, SPH particles and finite elements transfer mechanical information of the same material through solid contact, ensuring displacement coordination between the two via a point-to-surface bonding approach. Additionally, this method can fully utilize finite element boundary constraint methods to accurately analyze specific issues such as rock blast damage, thereby improving computational efficiency. Given the wide application of the LS-DYNA (R14) program in blasting engineering and its integration of the SPH algorithm, it effectively simulates the throwing and collision of rock fragments post-blast. This study employs LS-DYNA software and the FEM–SPH coupling simulation method.

The underground ore body of the Karatongk Mine in Xinjiang contains disseminated copper–nickel ore, primarily comprising mafic complex rocks that mostly exhibit a massive structure. Aside from the brecciated rock group, other ore groups display good stability, with the main structure of the ore rock being dense and hard, and the Protodyakonov coefficient (f) ranging between 8 and 10. The lithology is complex, consisting of rocks with massive structures such as norite, gabbro-norite, diabase, and norite. In addition to the mafic complex rocks, there are strata dominated by tuffaceous limestone and intermediate to acidic dikes. The geological structure is well-developed, with frequent occurrences of fault crush zones within the rock mass, and it is situated in an area active with neotectonic movements. The blasting height of the ore rock is 4.5 m, with a width of 4.6 m and an advance of 2.3 m. A total of six rows of blast holes are arranged. The distance from the first row of hole collars to the roof is 1.15 m. The row spacing between the first to second to third rows of the blast holes is 0.6 m, and the row spacing between the third to fourth to fifth to sixth rows of the blast holes is 0.65 m. The spacing between blast holes in each row is 1.1 m. The angle between the first row of blast holes and the driving face is 64°, while the angle between the second row of blast holes and the driving face is 78°. The drilling layout, initiation network, and parameters are shown in Figure 1a,b and

Table 1. Blasting parameters.

Number	Index	Value
1	Number of holes (pcs)	30
2	Hole depth (m)	First row 2.4 m, remaining rows 2.3 m
3	Number of electronic detonators (rounds)	30
4	Charge rolls per hole (pcs)	First row 5, remaining rows 4
5	Length of charge per hole (m)	1500
6	Total charge per cycle (kg)	37.5
7	Feet per cycle (m)	2.3
8	Cubic breakage per cycle (m3)	47.61
9	Explosives consumption per cycle (kg/t)	0.252

, respectively.

Based on the actual engineering situation, an FEM–SPH model was established using LS-DYNA software, as shown in Figure 2. The model dimensions are 3.0 m in length, 9.6 m in width, and 10.0 m in height. It consists of 266,076 SPH (Smoothed Particle Hydrodynamics) particles and 204,595 FEM (finite element method) elements. The thickness of the outer layer FEM mesh is 1.5 m, and the ore rock is modeled using SPH particles.

2. Materials and Methods

2.1. Mechanics of the RHT Model

LS-DYNA is a dynamic finite element method program that offers various models for simulating rock or rock-like materials. Among the commonly used models is the Holmquist–Johnson–Cook (HJC) model, one of two proposed by Timothy Johnson and Gordon Holmquist [16]. The second is the JH-2 constitutive model [17,18,19]. Another well-known material model in the literature is the Riedel–Hiermaier–Thoma (RHT) model [20]. The RHT model is widely used to assess the damage in concrete or rock materials by considering three stress limit surfaces to account for strength reduction and strain rate effects, as detailed by Borrvall and Riedel [21]. It describes the continuous accumulation of strain from elastic to plastic, where the accumulation of plastic strain causes damage, leading to model softening. Michał Kucewicz et al. [22] found that the JHC model could not initiate tensile cracking and only fully damaged materials with high compressive residual strength. Kong et al. [23] analyzed the advantages and disadvantages of the JHC and RHT models. The JHC model uses two stress invariants to represent the current damage surface but cannot capture the shear dilatancy behavior during the transition from low to high pressure. While the RHT model better describes material failure characteristics, it has more complex constitutive parameters than the HJC model, requiring significant effort to determine these parameters. To more accurately investigate the rock fragmentation process, this study employs the RHT model combined with a coupled FEM–SPH (finite element method–Smoothed Particle Hydrodynamics) approach to study rock movement during blasting. This methodology allows for a more precise simulation of the complex behaviors and interactions during the blasting process, providing deeper insights into the mechanisms of rock movement.

2.2. Calibration of RHT Model Parameters

The RHT model consists of 38 parameters, with 4 defaulted by the software and 34 requiring calibration.

2.2.1. Material Mechanical Parameters

Mechanical parameters were determined through mechanical testing, combined with data from the 1985 “Xinjiang Karatongk Copper–Nickel Mine Rock Mechanics Test Research” by the China Nonferrous Engineering Design and Research Institute, as shown in

Table 2. Basic physical parameters of the mineral rock model.

Item Name	ρ (kg/m3)	P (MPa)	Pcut (MPa)	E (GPa)	μ	v0 (m/s)	α0 %
Mineral rock	3121	116.3	6.06	54.4	0.25	4591	1.1

.

2.2.2. $p-\alpha$ Compaction Equation of State

In the $p-\alpha$ equation of state, the shear and pressure components are coupled. Pressure is expressed using the Mie–Gruneisen form, incorporating a polynomial Hugoniot curve and an A-B compaction relationship [24].

$$P_R={\displaystyle \frac{1}{{\alpha }_0}}\left(\left(B_0+B_1\mu \right){\alpha }_0{\rho }_{0}e+A_1\mu +A_2{\mu }^2+A_3{\mu }^3\right)$$

(1)

where P_R is the EOS pressure in the RHT model; α₀ represents the initial porosity parameter of the ore rock; A₁, A₂, and A₃ are the polynomial coefficients; B₀ and B₁ are material constants; μ is volumetric strain; and ${\rho }_0$ and $e$ are the initial density and the internal energy per unit mass of the rock, respectively. When T₂ = 0, B₀, B₁, A₁, A₂, and A₃ are derived from the Rankine–Hugoniot and Mie–Gruneisen equation of state [25].

$$B_0=B_1=2s-1,$$

(2)

$$A_3={\rho }_0{v}_0^2,$$

(3)

$$A_3={\rho }_0{v}_0^2(2s-1),$$

(4)

$$A_3={\rho }_0{v}_0^2(3s^2-4s+1),$$

(5)

where s is an empirical constant specific to the rock. Meyers, in his monograph [26], provides values for the material parameter s. For the magmatic intrusion mine, the value of s = 1.14.

2.2.3. Strain Rate Parameter Calibration

In the RHT model, the hydrostatic pressure is represented through Hugoniot polynomial coefficients and the p − α compression relationship. Under quasi-static conditions, the failure strength of the rock can be obtained by triaxial compression tests. The dynamic failure strength is obtained from the static breaking strength, and the failure surface stress is described as the yield surface through the compressive strength, regularized yield function Willam–Warnke function [27]:

$${\sigma }_y\left(P_0^{\mathrm{*}},{\dot{\epsilon }}_p,{\epsilon }_p^{\mathrm{*}}\right)=f_c{\sigma }_y^{\mathrm{*}}\left(P_0^{\mathrm{*}},F_r\left({\dot{\epsilon }}_p\right),{\epsilon }_p^{\mathrm{*}}\right)R_3\left({\theta }_l,P_0^{\mathrm{*}}\right),$$

(6)

where $f_c$ is the unconfined uniaxial compressive strength; ${\sigma }_y^{*}$ is the normalized yield function; R₃ is the Lode angle function, which accounts for the reduced strength of shear and tensile meridians [13]; ${\theta }_l$ is the Lode angle; Fr is the dynamic strain rate increase factor; $p_0^{*}$ is the normalized pressure (the parameters with a superscript * in the paper are the result of the normalization process); ${\stackrel{.}{ \epsilon }}_p$ is the strain rate; ${\mathsf{\epsilon }}_p^{\mathrm{*}}$ is the effective plastic strain; and $p_0^{*}=p_0/p_c$, $p_0$ is the hydrostatic pressure.

$$F_r\left({\dot{\epsilon }}_p\right)\left\{\begin{array}{c}{({\dot{\epsilon }}_p/{\dot{\mathsf{\epsilon }}}_0^{\mathrm{c}})}^{{\beta }_c} P\ge f_c/3\\ {\displaystyle \frac{P+f_t/3}{f_c/3+f_t/3}}{({\dot{\epsilon }}_p/{\dot{\mathsf{\epsilon }}}_0^{\mathrm{t}})}^{{\beta }_c}-{\displaystyle \frac{P-f_t/3}{f_c/3+f_t/3}}{({\dot{\epsilon }}_p/{\dot{\mathsf{\epsilon }}}_0^{\mathrm{c}})}^{{\beta }_t} -f_t/3<P<f_c/3 \\ {({\dot{\epsilon }}_p/{\dot{\mathsf{\epsilon }}}_0^{\mathrm{t}})}^{{\beta }_t} P\le -f_t/3\end{array},\right.$$

(7)

$${\beta }_c={\displaystyle \frac{4}{20+3f_c}},$$

(8)

$${\beta }_t={\displaystyle \frac{2}{20+f_c}},$$

(9)

where ${\dot{\mathsf{\epsilon }}}_0^{\mathrm{c}}$ is the reference strain rate under compression, ${\dot{\mathsf{\epsilon }}}_0^{\mathrm{t}}$ is the reference strain rate under tension, f_t is the tensile strength, ${\mathsf{\beta }}_{\mathrm{c}}$ and ${\mathsf{\beta }}_{\mathrm{t}}$ are the material constants for compression and tension, ${\dot{\mathsf{\epsilon }}}_0^{\mathrm{c}}=3.0\times {10}^{-5} {\mathrm{s}}^{-1},{\dot{\mathsf{\epsilon }}}_0^{\mathrm{t}}=3.0\times {10}^{-6} {\mathrm{s}}^{-1},{\dot{\mathsf{\epsilon }}}^{\mathrm{c}}=3.0\times {10}^{25} {\mathrm{s}}^{-1},{\dot{\mathsf{\epsilon }}}^{\mathrm{t}}=3.0\times {10}^{25} {\mathrm{s}}^{-1}$, respectively [28,29].

2.2.4. Material Damage Model Parameters

The model damage curve is defined as:

$${\sigma }_f^{*}(P_0^{*},F_r)=\left\{\begin{array}{c}A{(P_0^{*}-F_r/+{(A/F_r)}^{-1/N})}^N\\ F_r{F}_s^{*}/Q_1+3P_0^{*}(1-F_s^{*}/Q_1)\\ F_r{F}_s^{*}/Q_1-3P_0^{*}(1/Q_2-F_s^{*}/(F_t^{*}{Q}_1\left)\right)\\ 0\end{array}\right.\begin{array}{cc}& 3P_0^{*}\ge F_r\\ & F_r>3P_0^{*}\ge 0\\ & 0>3P_0^{*}\ge 3P_t^{*}\\ & 3P_t^{*}>3P_0^{*}\end{array},$$

(10)

When the material is under quasi-static loading conditions, $F_r$ takes the value of 1. Parameters A and N can be obtained from the expression for the damage surface at $3P_0^{*}\ge F_r$. Thus, the failure surface expression when $3P_0^{*}\ge F_r$ is as follows:

$${\sigma }_f^{*}\left(P_0^{*}\right)=A{(P_0^{*}-1/3+{\left(A\right)}^{-1/N})}^{N}3P_0^{*}\ge F_r,$$

(11)

where ${\mathsf{\sigma }}_{\mathrm{f}}^{\mathrm{*}}$ is the normalized strength, ${\mathsf{\sigma }}_{\mathrm{f}}^{\mathrm{*}}={\mathsf{\sigma }}_{\mathrm{f}}/{\mathsf{\sigma }}_{\mathrm{c}}$, and A and N are the failure surface parameters. The results of normalized p and Y are obtained from the experimental data of triaxial rock mechanics in

Table 3. Triaxial rock mechanics experiments and normalization parameters.

σ2 = σ3/MPa	σ1/MPa	σc/MPa	P*	σ*
0	116.3	116.3	0.333	1
5	170.1	116.3	0.516	1.419
6	178.9	116.3	0.547	1.487
7	187.8	116.3	0.578	1.554
10	213.5	116.3	0.669	1.749

* denotes the normalized (regularized) dimensionless parameter to eliminate the effect of units in the ontological relationship.

, with P as x-axis and Y as y-axis, and the curve of mode 3 is obtained by linear fitting to obtain A = 2.4, N = 0.82, as shown in Figure 3.

The p and y in Table 3 are obtained from the following Equations (12) and (13) [30]:

$$P_0^{*}=({\mathsf{\sigma }}_1+{\mathsf{\sigma }}_3)/3f_{\mathrm{c}},$$

(12)

$$P_0^{*}=({\mathsf{\sigma }}_1-{\mathsf{\sigma }}_3)/f_{\mathrm{c}},$$

(13)

In the RHT model, the damage variable is defined as the accumulation with plastic strain ${\epsilon }_p$:

$$D=\sum {\displaystyle \frac{d{\epsilon }_p}{{\epsilon }_p^f}}$$

(14)

where ${\epsilon }_p^f$ is the plastic strain at failure.

When the stress state reaches the material’s ultimate strength on the damage surface, damage accumulates gradually during inelastic deformation or plastic strain. In this model, the plastic strain at the point of damage is defined as [31]:

$${\epsilon }_p^f=\left\{\begin{array}{c}{D}_1{[P_0^{*}-(1-D\left)P_t^{*}\right]}^{D_2}\\ {\epsilon }_p^m\end{array}\right.\begin{array}{cc}& P_0^{*}\ge (1-D)P_t^{*}+{\left({\epsilon }_p^m/D_1\right)}^{1/D_2}\\ & P_0^{*}<(1-D)P_t^{*}+{\left({\epsilon }_p^m/D_1\right)}^{1/D_2}\end{array},$$

(15)

where ${\epsilon }_p^m$ is the minimum damaged residual strain and D₁ and D₂ are the damage constants for the RHT model. The variables are based on the findings in the literature [32], D₁ = 0.04, D₂ = 1.0.

In the RHT model, the maximum reduction in strength at the stretching meridian is described by the factor of the following equation:

$$\mathrm{Q}={\mathrm{Q}}_0+{\mathrm{B}}_{\mathrm{Q}}{\mathrm{P}}_0^{\mathrm{*}},$$

(16)

where Q₀ and B_Q are introduced as lode angle-dependence factors to describe the maximum reduction in strength on the tensile meridian. Thus, the tensile strength can be obtained by multiplying Equations (1) and (2). Q₀ and B_Q values can be obtained from the results of biaxial compression or triaxial tension. In this study, Q₀ = 0.68 and B_Q = 0.0105 were determined using the method of Li [33].

In summary, the mineral rock RHT parameters are calibrated as shown in

Table 4. RHT model parameters.

Parameters	Values	Parameters	Values
Mass density RO (kg/m3)	3100	Elastic shear modulus SHEAR (GPa)	21.76
Compressive strength FC (MPa)	116.3	Erosion plastic strain EPSF	2.0
Compressive strain rate dependence exponent BETAC	0.011	Relative tensile strength FT*	0.052
Compressive yield surface parameter GC*	0.53	Relative shear strength FS*	0.187
Crush pressure PEL (MPa)	77.53	Residual surface parameter AF	1.62
Compaction pressure PCO (GPa)	6.0	Residual surface parameter AN	0.62
Damage parameter D1	0.04	Reference compressive strain rate EOC	3.0 × 10−5
Damage parameter D2	1.0	Reference tensile strain rate EOT	3.0 × 10−6
Tensile strain rate dependence exponent BETAT	0.014	Break compressive strain rate EC	3.0 × 10−25
Tensile yield surface parameter GT*	0.7	Break tensile strain rate ET	3.0 × 10−25
Hugoniot polynomial coefficient A1 (GPa)	65.78	Failure surface parameter A	2.40
Hugoniot polynomial coefficient A2 (GPa)	84.19	Failure surface parameter N	0.823
Hugoniot polynomial coefficient A3 (GPa)	22.28	Lode angle dependence factor Q0	0.68
Parameter for polynomial EOS B0	1.28	Lode angle dependence factor B	0.0105
Parameter for polynomial EOS B1	1.28	Shear modulus reduction factor XI	0.5
Parameter for polynomial EOS T1 (GPa)	65.78	Minimum damaged residual strain EPM	0.015
Parameter for polynomial EOS T2	0.0	Gruneisen gamma GAMMA	0.0
Porosity exponent NP	3.0	Initial porosity ALPHA	1.0
Volumetric plastic strain fraction in tension PTF	0.001

.

2.3. The Explosive Model

In the simulation, *MAT_HIGH_EXPLOSIVE_BURN was applied to model the explosives. The (JWL) EOS was used to characterize the relationship between the change in pressure and the relative volume of rock during blasting. The JWL-EOS equation of state is expressed as [34]:

$$P=A(1-{\displaystyle \frac{\omega }{R_{1}V}})e^{-R_{1}V}+B(1-{\displaystyle \frac{\omega }{R_{2}V}})e^{-R_{2}V}+{\displaystyle \frac{\omega E}{V}}$$

(17)

where P is the burst pressure, E is the internal energy per unit volume, V is the relative volume, and A, B, R₁, R₂ and ω are material constants.

Table 5. Explosive material parameters.

Density/(kg·m−3)	Velocity of Donation/(m·s−1)	Pcj/GPa	A/GPa	B/GPa	R1	R1	ω	E0 (GPa)
1320	6690	16	586	21.6	5.81	1.77	0.282	7.38

shows the parameters used in the simulations, which were taken from the literature [35].

2.4. The Blasting Rock Mass Throwing Model

The underground blasting space was relatively confined, and the blasting operations in mine tunnels imposed multiple constraints on rock mass throwing. Research on the movement of blasted rock in underground environments was relatively scarce, with most studies focusing on rock mass movement in open-pit bench blasting. Many scholars had used the principle of spherical charge to approximate cylindrical charges as multiple spherical charges for the study of blasting issues. In this paper, based on the principle of spherical charge blasting, a model for underground blasting and rock throwing was constructed (Figure 4). In a two-dimensional cross-section of the blast hole, the radial lines emanating outward from the initiation point of the spherical charge divided the underground rock mass into multiple triangular finite elements. Assuming that the momentum imparted by the explosive energy was the same on different finite elements, according to the principle of spherical charge blasting, the movement velocity and direction of an individual element could be determined by the vector summation of the initial throw velocities imparted by multiple spherical charges. The direction vector of the initial throw velocity acting on a specific element, generated by the explosive energy of several spherical charges, determined its movement direction. This method allowed for the determination of the blasting rock throwing velocity direction for any rock element within the underground rock mass.

As shown in Figure 4, the underground blasting charge segment (BC) in this study is equivalently represented by N spherical charges, spaced equally along the axial direction. The mining rock mass in the stope is then divided into multiple triangular finite elements. The i-th rock mass element on the x-axis and the j-th rock mass element on the y-axis are denoted as ${\mathrm{M}}_{\mathrm{i}\mathrm{j}}$, with an initial throw velocity of ${\mathrm{V}}_{\mathrm{i}\mathrm{j}}$. The initial throw velocity generated by the n-th spherical charge at the rock mass element ${\mathrm{M}}_{\mathrm{i}\mathrm{j}\mathrm{n}}$ is denoted as ${\mathrm{V}}_{\mathrm{i}\mathrm{j}\mathrm{n}}$. At point ${\mathrm{M}}_{\mathrm{i}\mathrm{j}}$, the angle between the line connecting the rock mass element and the n-th spherical charge with the x-axis is ${\mathsf{\theta }}_{\mathrm{i}\mathrm{j}\mathrm{n}}$. Thus, the initial throw velocities in the x-axis and y-axis directions are denoted as ${\mathrm{V}}_{\mathrm{M}\mathrm{i}\mathrm{j}\mathrm{x}}$ and ${\mathrm{V}}_{\mathrm{M}\mathrm{i}\mathrm{j}\mathrm{y}}$, respectively.

After equivalently representing the spherical charges, the initial throw velocities of the rock mass element ${\mathrm{M}}_{\mathrm{i}\mathrm{j}}$ in the x-axis and y-axis directions can be expressed as shown in Equations (18) and (19):

$${\mathrm{v}}_{{\mathrm{M}}_{\mathrm{i}\mathrm{j}}\mathrm{X}}={\sum }_{\mathrm{n}=1}^{\mathrm{N}}{\mathrm{v}}_{\mathrm{i}\mathrm{j}\mathrm{n}}\cos {\mathsf{\theta }}_{\mathrm{i}\mathrm{j}\mathrm{n}},$$

(18)

$${\mathrm{v}}_{{\mathrm{R}}_{\mathrm{i}\mathrm{j}}\mathrm{y}}={\sum }_{\mathrm{n}=1}^{\mathrm{N}}{\mathrm{v}}_{\mathrm{i}\mathrm{j}\mathrm{n}}\mathrm{s}\mathrm{i}\mathrm{n}{\mathsf{\theta }}_{\mathrm{i}\mathrm{j}\mathrm{n}} ,$$

(19)

The vector resultant velocity in the x-axis and y-axis directions, which represented the initial throw velocity direction, is given by Equation (20). The throw angle, representing its throw inclination, is expressed through the inverse trigonometric function in Equation (21).

$${\mathrm{v}}_{{\mathrm{R}}_{\mathrm{i}\mathrm{j}}}=\sqrt{{\mathrm{v}}_{{\mathrm{R}}_{\mathrm{i}\mathrm{j}}\mathrm{y}}^2+{\mathrm{v}}_{{\mathrm{R}}_{\mathrm{i}\mathrm{j}}\mathrm{x}}^2},$$

(20)

$${\mathsf{\theta }}_{\mathrm{i}\mathrm{j}}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left({\displaystyle \frac{{\mathrm{v}}_{{\mathrm{R}}_{\mathrm{i}\mathrm{j}}\mathrm{y}}}{{\mathrm{v}}_{{\mathrm{R}}_{\mathrm{i}\mathrm{j}}\mathrm{x}}}}\right),$$

(21)

3. Results

In blasting engineering, the collision characteristics of rock movement following multi-row blast-hole detonation in underground mines are rarely studied. In underground mining blasting projects, the size of rock fragments are influenced not only by the energy of the initial blast but also by subsequent collisions that occur as the rock moved. These collisions could cause the rock mass to fragment again (once or multiple times). A numerical calculation and analysis of the model were conducted using the FEM–SPH method. To study the movement characteristics of the rock under precise delay time conditions, this research investigated the collision patterns of rock movement under precise delay-controlled detonation networks with delay times of 25 ms, 50 ms, and 75 ms, respectively.

3.1. Rock Displacement Patterns Under Different Delay Time Conditions

The Ls-PrePost (R14) software was utilized to perform box cutting on the numerical model, obtaining the post-blast half-region model of the rock mass. This approach was employed to study the internal rock mass movement and collision patterns under the influence of blasting.

Under the 25 ms delay time condition, as depicted in Figure 5, it was observed that before the detonation of the second row of blast holes, the farthest movement distance of the rock mass was 0.41 m (25 ms (a1)). After the detonation of the second row of blast holes, the rock mass exhibited an upward and outward movement trend, with the overall outward movement being 1.19 m, which was only slightly larger than the 1.09 m movement observed after the first row of blast holes with a 50 ms delay, and significantly less than the 2.307 m movement under a 75 ms delay. This indicated that the short inter-row delay resulted in a weaker outward throwing effect on the rock mass. Additionally, it inhibited the movement of the internal rock mass, leading to a more pronounced restrictive effect (25 ms (b1)). Following the detonation of the third row of blast holes, the rock mass in the areas of both the first and second rows of blast holes did not effectively throw outward, creating a movement trend that restricted the outward movement of the rock mass in the area of the third row of blast holes (25 ms (c1)). After the detonation of the fifth row of blast holes, since the rock in the areas of the previous multiple rows of blast holes did not effectively throw outward, the internal rocks collided with each other and the restrictive effect became more evident, and the rock at the bottom of the upper holes was not effectively thrown out (Figure 5, 25 ms (e1)). Finally, after the detonation of the bottom holes, the throwing effect of the rock in the compensation space was worse than under the 50 ms and 75 ms delay conditions, forming a concave groove at the bottom, indicating a root phenomenon (Figure 5, 25 ms (f1)).

Under the 50 ms delay time condition, as illustrated in Figure 5, it was evident that prior to the detonation of the second row of blast holes, the farthest movement distance of the rock mass was 1.09 m, exhibiting a stronger effect than the 25 ms delay but weaker than the 75 ms delay (50 ms (a2)). Upon the detonation of the second row of blast holes, the rock mass demonstrated a movement trend inclined upward at a 30° angle. However, compared with the 75 ms delay, the rock mass in the area of the first row of blast holes did not completely throw out, resulting in a squeezing and restrictive effect with the rock mass moving outward after the detonation of the second row of blast holes (50 ms (b2)). Following the detonation of the third row of blast holes, the movement distance of the rock mass in the area of the second row of blast holes was significantly less than that observed under the 75 ms delay, and the upper rock mass did not fully throw out, leading to a squeezing and restrictive effect with the rock mass in the area of the second row of blast holes (50 ms (c2)). After the detonation of the fourth row of blast holes, the restrictive phenomenon was pronounced but less pronounced than the 25 ms delay and more pronounced than the 75 ms delay (50 ms (d2)). Upon the detonation of the fifth row of blast holes, the rock mass at the bottom of the uppermost holes was not thrown out (50 ms (e2)). Following the detonation of the bottom holes, the throwing effect of the rock mass within the compensation space was considerably stronger than the 25 ms delay but weaker than the 75 ms delay, and a clear root phenomenon was discernible (50 ms (f2)).

Under the 75 ms delay blasting condition, as illustrated in Figure 5, it was evident that before the detonation of the second row of blast holes, the farthest throw distance of the rock mass was 2.307 m, which was significantly greater than the throw distances observed under the 25 ms and 50 ms delays, providing a larger compensation space for the detonation of subsequent rows (75 ms (a3)). After the detonation of the second row of blast holes, the rock mass at the excavation face had been thrown out, and the compensation space in the area of the second row of blast holes was larger compared with the 25 ms and 50 ms delays (75 ms (b3)). Following the detonation of the third row of blast holes, the rock mass from the first two rows had the farthest throw distance among the three delay conditions, and it did not affect the outward movement of the rock mass in the area of the third row of blast holes (75 ms (c3)). After the detonation of the fourth row of blast holes, the outward movement trajectory of the rock mass within the upper compensation space was clear, and it did not impose a restrictive effect on the movement of the lower rock mass, with the movement distance of the rock mass near the roof being larger than that under 25 ms and 50 ms delays (75 ms (d3)). Upon the detonation of the fifth row of blast holes, the rock mass in the upper area had been thrown out, providing a compensation space for the movement of the lower rock mass, and no restrictive effect occurred (75 ms (e3)). After the detonation of the bottom holes, a concave groove was also formed, indicating a root phenomenon, but the amount of rock mass in the upper area was significantly less than that under 25 ms and 50 ms delays, with more rock mass having been thrown beyond the free face (75 ms (f3)).

Based on the study of the internal movement patterns of rock mass, a demonstration analysis of the overall movement trend of blasted rock mass was conducted, as illustrated in Figure 6a–e. These figures depicted the outward throwing movement cloud maps of the rock mass after the detonation of the second to sixth rows of blast holes, respectively, reflecting the spatiotemporal evolution pattern of rock mass movement at the excavation face post-blasting. From these cloud maps, it was observed that the rock mass at the excavation face first fell to the bottom of the blast heap, exhibiting a top-down collapse of the rock mass at the excavation face and a bottom-up stacking effect of the blast heap. In the simulation, tensile fracture phenomena appeared in the FEM elements at the upper roof. Subsequent field blasting tests confirmed that this manifested as scaling on both sides of the tunnel at the roof, consistent with the tensile fracture observed in the FEM elements at the roof in the numerical simulation, as shown in Figure 6f–h.

3.2. Stress and Damage Characteristics of Ore Rocks Under Different Delay Time Conditions

To further analyze the movement patterns of the rock mass after blasting, six monitoring lines were set at elevations of 4.0 m, 3.35 m, 2.7 m, 2.05 m, 1.4 m, and 0.75 m in the FEM–SPH model. Points were sampled every 0.46 m, from the bottom of the borehole to the excavation face, labeled as A, B, C, D, E, and F, to monitor the movement patterns of the rock mass.

From Figure 7, it was evident that under a 25 ms delay condition, at the 4.0 m elevation, except for the rock mass 0.46 m from the excavation face which moved upward, the rest of the rock mass moved downward, with the smallest movement distance at 0.46 m from the bottom of the borehole. At the 3.35 m elevation, the rock mass within 0.46 m of the excavation face moved upward, and the angle of upward movement at the excavation face was greater than that in the interior regions, with the smallest movement distance at the bottom of the borehole. At the 2.7 m elevation, the rock mass at the bottom of the borehole did not move outward, but the rock mass 0.46 m from the bottom of the borehole moved further than before, showing a complete movement trend. The movement of the rock mass at lower elevations maintained the same trend without significant changes.

Under a 50 ms delay condition, the rock mass at 4.0 m moved downward as a whole, with the rock mass 0.46 m from the bottom of the borehole first moving downward and then upward, and the movement distance was larger than at 25 ms. The movement pattern at 3.35 m was the same as at 25 ms, but the angle of upward movement within 0.46 m of the excavation face was greater than at 25 ms. At the 2.7 m elevation, both the rock mass at the bottom of the borehole and the excavation face moved upward, with more intense upward movement in the region of the borehole bottom. At the 2.05 m elevation, only the rock mass at the excavation face moved upward, while the rest moved downward. The movement patterns at lower elevations were similar to those at 25 ms.

Under a 75 ms delay condition, the movement of the rock mass in the upper regions was more complex, showing collision effects and trajectory shifts during outward movement at the 4.0 m and 3.35 m elevations. At the 2.7 m elevation, the outward movement distance of the rock mass was greater than at 25 ms and 50 ms, but the upward movement capability was less than at 25 ms and 50 ms. At the 2.05 m elevation, the rock mass inclined upward more noticeably. At the 1.4 m elevation, the rock mass movement was more inclined towards the horizontal direction, and at the lowest 0.75 m elevation, the upward movement trend was stronger than at 25 ms and 50 ms. However, in general, the upper rock mass moved further after blasting and fragmentation, and as the elevation decreased, the movement distance of the rock mass also decreased, with the outward movement distance increasing with the increase in delay time between rows. To thoroughly analyze the movement patterns of rock mass post-blasting, SPH particles were used as monitoring points at each row of boreholes in the vertical direction at the excavation face within the numerical analysis model. The velocity time–history curves under 25 ms, 50 ms, and 75 ms delay conditions were obtained, as shown in Figure 8.

Under the 25 ms delay condition, the initiation velocities of the second and fourth rows of boreholes increased most significantly, exhibiting a secondary velocity increase due to the apparent effects of subsequent blasting. For the first, third, fifth, and sixth rows, the velocity increase due to subsequent row explosions was not significant, displaying a slow and steady increase. Post-initiation, a velocity decrease was observed across all boreholes, with a more pronounced decrease in the second and fourth rows, while the other rows showed minor decreases. This phenomenon might have been attributed to the short delay time, where the detonation of subsequent boreholes continuously applied explosive forces to the preceding rock masses, causing collisions and sustained velocity increases. Additionally, the second row displayed minimal velocity changes post-explosion, indicating difficulty in outward rock mass ejection.

Under the 50 ms delay condition, each subsequent borehole initiation promoted a velocity increase in the previously initiated rock masses. The velocities of rock masses in the first, second, and third rows continued to rise due to subsequent borehole explosions. After the initiation of the fourth, fifth, and sixth rows, the velocities showed significant decreases, but the velocities in the fourth and fifth rows exhibited abrupt increases after the subsequent borehole detonations. Each borehole’s velocity showed a more substantial drop post-initiation than at 25 ms, indicating a reduction in rock mass density and collision frequency in a confined space. Overall, rock mass velocities continued to accumulate, suggesting ongoing collisions leading to continuous velocity accumulation, showing a clear and consistent pattern. At this delay, the second row also exhibited minimal velocity changes post-explosion, but the movement speed was higher than at 25 ms, indicating an enhanced ability for rock mass ejection compared with 25 ms.

Under the 75 ms delay condition, the velocity time–history curves showed a different pattern compared with the 25 ms and 50 ms delays. The subsequent borehole initiations did not significantly enhance the velocities of the previously exploded rock masses. Post-initiation velocities of the first and second rows decreased continuously, indicating sparse rock mass density in the area, with most rock masses already ejected to the free face, preventing sustained velocity accumulation from subsequent explosions. The third and fourth rows exhibited slight velocity recoveries after initial decreases, suggesting some impact from subsequent explosions, but overall the rock mass density was low, and collision impacts were weak, with most rock masses effectively ejected. The fifth and sixth rows showed immediate velocity decreases post-initiation, without significant velocity fluctuations. Therefore, the velocity trend under the 75 ms delay indicated that the number of fragmented rock masses in the confined blasting space was sparse, with most already ejected to the free face. The rock masses in the second row, which initially moved slowly, showed continuous velocity increases post-explosion, indicating an enhanced ability for outward movement and ejection.

Huang Yonghui’s casting model experiments revealed that approximately 20 ms after explosive detonation, the acceleration of the fragmented material was completed [36]. This study found that in underground mine blasting, the rock completed its initial acceleration within approximately 15 ms following the detonation of the explosives. Cui Xinnan et al. [37] observed that the ejection velocity of fragmented rock in bench blasting exhibited a pattern of acceleration, constant velocity, secondary acceleration, and deceleration. This study also identified similar patterns; however, these patterns were not evident with a short inter-row delay time of 25 ms but became more pronounced with a longer inter-row delay time of 75 ms.

4. Comparative Verification of Rock Mass Movement and Ejection Blasting Experiments Under Different Delay Times

In the field blasting comparison and verification experiments for rock mass movement, the confined space within underground tunnels and the constraints of objective conditions led to a relatively concentrated pile of blasted material. Nevertheless, the ejected material predominantly accumulated within 10 m of the excavation face under different delay time conditions as follows.

Under the 25 ms delay condition, the blasting ejection distance was more concentrated near the excavation face, as illustrated in Figure 9 (25 ms (a)). The ejection distance was shorter, with the peak accumulation of the blasted pile located approximately 4.2 m from the excavation face, and the top width of the pile was about 1.8 m, with a curved distribution downward.

Under the 50 ms delay condition, the blasting ejection distance was farther than at 25 ms, but the concentration degree of the blasted pile slightly decreased. The pile’s accumulation pattern was similar to that at 25 ms, with the peak accumulation distance approximately 5.1 m from the excavation face and the top width about 1.2 m, displaying a curved distribution downward, as shown in Figure 9 (50 ms (c)).

Under the 75 ms delay condition, the distribution characteristics of the blasted pile differed from those at 25 ms and 50 ms. The primary feature was a significant decrease in concentration, and no steep “peak” wave-like shape appeared near the excavation face. Instead, it exhibited a gentle “slope” shape, with distribution characteristics shown in Figure 9 (75 m s(e)).

The overall ejected and accumulated pile characteristics under the three different delay time conditions are illustrated in Figure 9 (25 ms (b), 50 ms (d), and 75 ms (f)). The ejection distances of the piles were measured during the field experiments. Under the 25 ms delay condition, the pile was distributed within a 9.6 m range; under the 50 ms delay condition within a 13 m range; and under the 75 ms delay condition within a 17.2 m range. Hence, within a certain delay time range, the ejection distance of the blasted pile increased as the inter-row delay time increased, showing a transition from a concentrated accumulation near the excavation face to a more dispersed distribution towards the farther end. Additionally, a shorter inter-row delay resulted in a steeper pile, while a longer inter-row delay produced a gentler slope. The slope angle of the pile decreased with increasing inter-row delay time within a certain range. Li, X.L et al. [38] focused on the control of large fragments and concluded that the inter-row delay time for high bench cast blast should be between 63.5 ms and 102 ms, and should not exceed 100 ms. The findings of this study, indicating that a 75 ms inter-row delay time yielded better blasting results, were consistent with this conclusion.

From the field blasting experiments, it was evident that the rock mass ejection patterns under different delay times were consistent with the trends observed in the numerical simulation analysis.

5. Discussion

(1)

The geological conditions underground in the study area of this paper are exceedingly complex, with varying degrees of joint and fissure development among different ore bodies and disparities in lithological weaknesses. The numerical simulations in this study assume a homogeneous, monolithic ore rock body, providing conclusions of general regularity. However, these conclusions may exhibit distinct characteristics under varying geological conditions. Future research could focus on ore rock bodies where joint and fissure development results in suboptimal blasting effects.

(2)

Constrained by computational demands, particularly the high processing power required by the FEM–SPH method, this study omits the post-fracture collision scenarios between ore rock and tunnels. Even under simplified boundary conditions, the computation extended beyond 200 h. Accounting for such collisions would render solutions infeasible on standard computational equipment, as the real scenario involves complex interactions where ore rock rebounds and accumulates into a muck pile, challenging accurate quantification.

(3)

The dynamics of ore rock movement are intricately linked to the size of fragments resulting from blasting. Future research could develop fragment size prediction models grounded in big data analytics theory, utilizing more extensive field-measured data to furnish more advanced predictive techniques for blasting operations.

6. Conclusions

This paper, based on practical engineering, focuses on the movement and ejection of rock mass during underground mining blasts. A FEM–SPH numerical model was constructed through three-dimensional modeling, and LS-DYNA numerical simulation software was employed to investigate the movement patterns of rock mass under precise delay blasting in underground environments. The study aimed to explore the spatiotemporal evolution of rock mass movement post-blasting and verify the influence mechanism of precise delay times on rock mass movement patterns through comparative field blasting experiments. The dynamic response characteristics of rock mass movement after underground access stope blasting were revealed through precise delay initiation using digital electronic detonators under 25 ms, 50 ms, and 75 ms delay conditions.

(1)

Numerical simulations of multi-row borehole blasting using the SPH method revealed the movement characteristics of rock mass. The results indicated that within a certain delay time range, the ejection distance of rock mass increased with increasing inter-row delay time, exhibiting a pattern where the accumulation was concentrated near the excavation face and dispersed towards the farther end.

(2)

Research on rock mass movement patterns under different delay times showed that at 25 ms and 50 ms delays, due to the short inter-row delay, the rock mass already affected by the explosion could not be effectively ejected to the free face during the delay period, resulting in significant restraint and inhibiting the creation of new free faces, which could lead to poor blasting outcomes. At 75 ms delay, the blasting effect could effectively eject the rock mass during the inter-row delay period, creating new free faces and providing effective compensation space for subsequent blasting, which had a beneficial impact on the blasting outcome.

(3)

Studies on the spatiotemporal evolution of rock mass movement under different delay times revealed that the pile morphology was significantly influenced by the delay time. Within a certain delay time range, the slope angle of the pile decreased with increasing inter-row delay time. Blasting piles with shorter inter-row delays were steeper, while those with longer inter-row delays were gentler.