1 Introduction

Stemming is an operation or a separation structure that separates explosive and atmosphere at the collar of a borehole or separates adjacent cartridges within a hole (Hustrulid et al. 2001; Liu et al. 2018b, 2019). It is one of the most crucial factors influencing the energy utilization of explosive detonation and thus controlling rock fragmentation in rock blasting of civil and mining engineering (Hosseini et al. 2023; Hustrulid et al. 2001; Langefors and Kihlström 1963; Zhang et al. 2022). However, up to now, although the stemming is employed, the energy utilization in rock blasting is low (Ouchterlony et al. 2004; Sanchidrián et al. 2007; Zhang 2016), and undesirable results of rock excavation and rock fragmentation such as under/overbreak, crushed rock particles and oversize fragmentation are frequently obtained (Cheng et al. 2021b; Hong et al. 2023; Li et al. 2023d, e; Rezaeineshat et al. 2020; Yi et al. 2017a, b). Furthermore, it is still unclear how much stemming affects rock cracking and rock fragmentation in blasting operations (Zhang 2016, 2019). Therefore, it is necessary to investigate rock fracturing and rock disintegration under blasting with different stemming conditions for increasing the utilization efficiency of explosion energy and improving the performance of rock fragmentation.

In practical blasting, sand or drill cutting is a common material poured into the borehole as the stemming. In most cases, sand is stemmed at the borehole collar to prevent air blasts and reduce gas escaping from the hole collar, leading to more explosion energy consumed in fragmenting rock mass (Brinkmann 1990; Zhang 2016). Moreover, sand is also employed to fill the gap between adjacent cartridges in a hole to form more shock wave interactions and prolong the duration of detonation-induced stress waves within the blasthole as well as reduce explosive consumption. This stemming structure is called a sand deck, and it is popularly employed in production blasting with deep-hole in mining, which aims to disintegrate rock mass efficiently with a low specific charge (Chen et al. 2021; Hustrulid et al. 2001; Zhang 2016). Sometimes, the air is used as the separation material filling the space between two cartridges in a hole or between explosive and stemming when the borehole collar is stemmed, which can be commonly found in perimeter blasts such as smooth blast and pre-split blast (Gao et al. 2023b; Mel’Nikov 1940; Zhang 2016). During the blasting, the air space induces efficient shock wave collisions within a blasthole, generating longer explosion stress waves with smaller pressure peak (Cheng et al. 2021b; Jhanwar et al. 2000). This air space structure is generally named an air deck (usually achieved using a thin paper tube), and similar to the sand deck, it was originally developed to reduce the crushed rock particles near the borehole, the explosive consumption and ground vibration, etc. (Zhang 2016). In actual rock blasting, the abovementioned stemming structures are all widely used to improve the performance of rock excavation. The stemming at the borehole collar is primarily applied in rock blasting to contain more explosion energy in borehole while the sand deck and air decks are mainly employed to adjust explosion stress and energy distribution within the hole.

It is generally accepted that a large part of detonation energy escapes through the borehole collar when no or too little stemming is used. With an unstemmed borehole, the explosion energy consumed in rock fracturing and rock fragmenting is much less than the energy generated by explosive burning. In contrast, when the stemming is properly employed, there will be no or little energy escapes so that much more explosion energy is dissipated in fragmenting rock (Hagan 1983; Zhang 2016; Zhang et al. 2020b, 2021). Laboratory experiment on the influence of stemming on blast-induced fracturing conducted by Dally et al. (1975) showed that the stemmed borehole leads to an increase in the energy carried in the explosion-induced pressure and produces longer and more radial fractures compared with an unstemmed hole. Site measurement of blast-formed gas pressure at blasthole collar in an underground gold mine implemented by Brinkmann (1990) indicated that without stemming, the high-speed explosion gases escape from hole collar to the atmosphere, which carries away more than 50% of detonation energy. In comparison, when sufficient stemming is used, the energy carried by explosion gases escaping from the hole collar is reduced by over 50%. Meanwhile, the duration of the explosion gas elapsing is extended by over 300%, which is very effective in creating and extending rock cracks. In recent years, Zhang et al. (2020b, 2021) carried out model tests on the influences of stemming on rock fragmentation generated by blasting, and the results showed that at a specific charge, the unstemmed blasting results in a waste of at least 25% of explosive energy in rock fragmentation when compared to stemmed blasting, and therefore blast with full stemming produces much smaller fragments than unstemmed and partially stemmed blasts. Most recently, experiments on blast-induced rock fracturing with iron ore samples considering different stemming structures were conducted by Ma et al. (2022) to study the effect of stemming on blast-formed internal fracturing within the sample using the combination of computed tomography scanning, digital image processing and 3D crack reconstruction. The experimental findings indicated that the initiation of a stemmed hole leads to the generation of more internal fractures compared to unstemmed hole detonation, resulting in a 33.6% increase in internal damage.

Since stemming can improve the efficiency of explosion energy dissipation in rock blasting, as indicated in the abovementioned studies, borehole is commonly stemmed in blasts of the open pit, quarries, raising, etc., for better rock fragmentation and reasonable excavation performance. However, stemming is generally not employed in upward-drilled blastholes of subsurface mining and it is also frequently absent in boreholes during tunnelling (Zhang 2016, 2019). The most important reason is that the findings in previous research are insufficient to convince blasting operators to use stemming for improving rock fragmentation performance. Till now, it is a lacking of quantitative investigation on clarifying how much stemming influences the blast-induced fragment size and fragment size distribution (FSD). Besides, few relevant studies were conducted about the dependence of rock fracturing and rock disintegration on the combination of stemming, air/sand deck and explosive, i.e., it is not clear yet how the variation of rock fragments is with distinct stemming structures. Therefore, it is necessary to quantitatively investigate the effects of the stemming condition on blast-produced rock fragmenting with different combinations of explosive, air/sand deck and stemming.

In the current study, the dependence of rock fragmentation on stemming is studied using combined finite element modelling in LS-DYNA (LSTC 2015) and image processing in ImageJ (Durda et al. 2015). First, the parameters of material models and equations of state (EOSs) employed in the developed numerical model are verified against the rock fracturing and FSD in three blasts with different combinations of stemming, air deck and explosive newly conducted by Zhang et al. (2021). Then, using the verified constitutive material models, a cubic single-hole model is constructed to model blast-produced pressure attenuation and rock cracking under different stemming conditions. The resulting crack pattern is obtained by post-processing the simulated numerical model in LS-PREPOST and is later image-processed to obtain the size data of blast-generated rock fragments. Thereafter, the corresponding FSDs determined by combining numerical modelling and image processing are fitted using a three-parameter generalized extreme value function, and the variations in blast-induced FSDs under different stemming conditions are quantitatively examined and discussed. This investigation provides some novel insights into the mechanisms of rock fragmenting induced by blasting with distinct stemming structures and supplies some guidance in improving rock fragmentation performance in blast engineering.

2 Numerical model and verification

In the present modelling, the finite element method with a widely used commercial hydrocode LS-DYNA is employed because of the advanced constitutive material models in LS-DYNA library. Using LS-DYNA, one can reasonably simulate the pressure during explosive burning and the rock damage behaviour which is crucial in simulating rock blasts (Gao et al. 2023a; Himanshu et al. 2022; Kucewicz et al. 2022; Li et al. 2023a, b, c; Pan et al. 2021; Pu et al. 2021; Zhao et al. 2019). Before the modelling, the developed finite element model is verified based on the rock cracking and FSD in blasting tests implemented by Zhang et al. (2021). The blast testing, material models and EOSs, and the numerical calibration results are introduced as follows:

2.1 Blast testing

In blast testing, cylinder granite specimens, 240 mm in diameter and 300 mm in height, were adopted. A central borehole with about 17 mm diameter was axially drilled from the top surface of the granite sample with a length of about 210 mm. All granite samples had no visible fracture or flaw. The granite had a density of 2650 kg/m3, Young’s modulus of 42.6 GPa, Poisson’s ratio of 0.23, and its longitudinal-wave velocity was 4200 m/s and compressive strength was 120 MPa. The explosive used in testing was Pentaerythritol tetranitrate (PETN), and it was placed at the bottom of borehole and initiated using an electrical detonator, as depicted in Fig. 1. The basic parameters for PETN are the density of 930 kg/m3, velocity of detonation (VOD) of 5300 m/s and explosion heat of 5925 kJ/kg. The charge length Le at borehole bottom was 15 mm for each sample except for sample S7 (20 mm). Sand was used as stemming material and stemmed at the collar of borehole whereas there was an air deck between stemming and explosive. During the blasting, the process of rock cracking was recorded using a high-speed camera. After running the test, the blast-induced FSDs were obtained by collecting, sieving/measuring and weighing the blasted fragments. Three blast tests of samples S6, S7 and S8 are chosen for calibrating the numerical model, and the geometry configuration of the selected samples is listed in Table 1.

Fig. 1
figure 1

Rock sample in blasting tests: a photograph and b diagrammatic sketch

Full size image
Table 1 Detailed geometry of samples S6, S7 and S8
Full size table

2.2 Material models and EOSs

2.2.1 Material model for rock

In rock blasting modelling, it is extremely important to select a constitutive material model for rock (Cheng et al. 2021a; Li et al. 2021; Liu et al. 2020). There are several constitutive models commonly used in modelling rock dynamic response under blasting such as the Concrete_Damage_Rel3 model (Liu et al. 2019, 2018b), HJC model (Luo et al. 2021; Wang et al. 2021), CSCM model (Tao et al. 2020) and RHT model (Hashemi and Katsabanis 2021; Jang et al. 2018; Katsabanis 2020; Leng et al. 2021; Li et al. 2020, 2022; Liu et al. 2018a, 2023). In this study, the advanced tensile-compressive damage model, RHT model, is chosen for modelling the rock behaviour under blasting because it takes strain rate for both compression and tension, strain hardening during material damage, and damage softening in compression, into account simultaneously. In the RHT model, the strength criteria are described using three stress limit surfaces including yield surface, failure surface and residual friction surface. Meanwhile, the pressure is described in Mie-Greisen form, embedding polynomial Hugoniot and pressure vs. porosity curve (Borrvall and Riedel 2011), which can precisely capture the real rock cracking behaviours under blasting. During the blast loading process, this model performs elastic before its pressure is lower than pore crush pressure, and once the material pressure is over the pore crushing pressure, pore collapse happens and results in a reduction in effective bulk modulus. As the stress state in rock continuously rises and exceeds its yield limit, the material deforms plastically. Subsequently, when this stress state reaches the failure surface, the material begins to accumulate damage gradually, accompanied by an increasing effective plastic strain. Finally, the stress state hits the last stress limit surface, i.e., residual surface, and the material is completely damaged.

When using this model, 37 parameters should be determined, and except for some parameters that have been given in this model, they can be obtained by laboratory testing, theoretical calculation and from reference as follows: the shear modulus of granite G can be determined by G = E/(2 × (1 + υ)), in which υ and E correspond to Poisson’s ratio and Young’s modulus, respectively, G = 17.32 GPa. Meanwhile. the rock tensile strength ft can be empirically calculated with ft = fc/10 (Li 2014b), in which fc denotes the uniaxial compressive strength of granite, ft = 12 MPa. Besides, the pore crushing pressure PEL is set as 2fc/3 according to Riedel et al. (2009), PEL = 80 MPa. In addition, the initial porosity for intact granite α0 = 1.016 is set referring to other works of Zhang et al. (2020a).

In the RHT model, the difference in strain rate effect between compression and tension is considered, and the dynamic strength increase in material is expressed using the strain strength factor Fr as (LSTC 2015)

$$F_{r} (\dot{\varepsilon }_{p} ) = \left\{ {\begin{array}{ll} {(\dot{\varepsilon }_{p} /\dot{\varepsilon }_{0}^{c} )^{{\beta _{c} }} } & {{\text{ }}(P \ge f_{c} /3)} \\ {\frac{{P + f_{t} /3}}{{f_{c} /3 + f_{t} /3}}(\dot{\varepsilon }_{p} /\dot{\varepsilon }_{0}^{c} )^{{\beta _{c} }} - \frac{{P - f_{c} /3}}{{f_{c} /3 + f_{t} /3}}(\dot{\varepsilon }_{p} /\dot{\varepsilon }_{0}^{c} )^{{\beta _{t} }} } & {( - f_{t} /3 < P < f_{t} /3{\text{ }})} \\ {(\dot{\varepsilon }_{p} /\dot{\varepsilon }_{0}^{t} )^{{\beta _{t} }} } & {(P \le - f_{t} /3)} \\ \end{array} } \right.$$
(1)

where \(\dot{\varepsilon }_{p}\) denotes strain rate. The reference strain rates for compression \(\dot{\varepsilon }_{0}^{c}\) and tension \(\dot{\varepsilon }_{0}^{t }\) are 3.0 × 10–8 s−1 and 3.0 × 10–9 s−1, respectively, which are given in the model. Meanwhile, the break compressive \(\dot{\varepsilon }^{c}\) and tensile \(\dot{\varepsilon }^{t}\) strain rates are also given, \(\dot{\varepsilon }^{c}\) = 3.0 × 1022 s−1 and \(\dot{\varepsilon }^{t}\) = 3.0 × 1022 s−1. P is the pressure, P = (σ1 + σ2 + σ3)/3. βc and βt refer to strain rate dependence exponents of compression and tension, respectively, and they are calculated by βc = 4/(20 + 3fc) = 0.0143 and βt = 1/(20 + fc) = 0.0105, respectively (Borrvall and Riedel 2011).

Meanwhile, the yield surface is described using fc, the normalized yield function σy.* and Willam–Warnke function R3 as (LSTC 2015)

$$\sigma_{y} (P^{*} ,\dot{\varepsilon }_{p} ,\varepsilon_{p}^{*} ) = f_{c} \sigma_{y}^{*} (P^{*} ,F_{r} (\dot{\varepsilon }_{p} ),\varepsilon_{p}^{*} )R_{3} (\theta ,P^{*} )$$
(2)

where P* = P/fc refers to the normalized pressure, \({\varepsilon }_{p}^{*}\) represents the effective plastic strain, θ denotes the lode angle. Moreover, the failure surface is expressed with relative tensile ft* and shear fs* strength, and failure surface parameters of A and N as (LSTC 2015)

$$\sigma _{f}^{*} (P^{*} ,F_{r} ) = \left\{ {\begin{array}{ll} {A(P^{*} - F_{r} /3 + (A/F_{r} )^{{ - 1/N}} )^{N} } & {(3P^{*} \ge F_{r} )} \\ {F_{r} f_{s}^{*} /Q_{1} + 3P^{*} (1 - f_{s}^{*} /Q_{1} )} & {(F_{r} > 3P^{*} \ge 0)} \\ {F_{r} f_{s}^{*} /Q_{1} - 3P^{*} (1/Q_{2} - f_{s}^{*} /Q_{1} f_{t}^{*} )} & {(0 > 3P^{*} \ge 3P_{t}^{*} )} \\ 0 & {(3P_{t}^{*} \ge 3P^{*} )} \\ \end{array} } \right.$$
(3)

where Q1 = R3(π/6,0) and Q2 = Q(P*), Pt* refers to the failure cut-off pressure, σf* denotes the normalized strength, σf* = σf / fc, in which σf = {[(σ1−σ2)2 + (σ2−σ3)2 + (σ1−σ3)2)]/2}0.5. Besides, fs* = τ/fc where \(\tau = 1.32f_{c} ((\sigma_{n} - f_{t} )/f_{c} )^{0.57}\) according to the work of Xie et al. (2017), fs* = 0.355, and ft* = ft/fc = 0.1. The A and N can be obtained by substituting the P* and σf* under various confining pressures into Eq. (3), and the P* and σf* can be estimated using Hoek–Brown criterion. With an intact granite, the Hoek–Brown criterion can be derived as σ1 = σ3 + fc (24σ3 / fc + 1)0.5, and thereby a mechanical parameter table of granite can be determined, as given in Table 2. Using P* and σf* under biaxial stress conditions (10 and 100 MPa), A = 2.76 and N = 0.72 are determined.

Table 2 Mechanical parameters of granite under confining pressure
Full size table

When the stress state in the RHT model hits its stress limit surface of failure, the material damage in the model is described with (Borrvall and Riedel 2011):

$$D = \sum \left( {\Delta \varepsilon_{p} /\varepsilon_{f} } \right)$$
(4)

where ∆εp and εf refer to accumulated effective plastic strain and failure strain, respectively, which can be used to identify rock cracking by critical value of D. This εf can be given as εf = D1 − (P* − (1 − D) Pt*)D2, in which D1 and D2 are damage constants, and D2 is given in the model as 1.0.

Besides, the pressure in p-α compaction EOS is written as (Borrvall and Riedel 2011)

$$P_{EOS} = \frac{1}{\alpha }((B_{0} + B_{1} \mu )\alpha_{0} \rho_{0} e + A_{1} \mu + A_{2} \mu^{2} + A_{3} \mu^{3} ) \, \mu > 0$$
(5)

where B0 and B1 are material constants and B0 = B1 = 2k − 1, in which k is also a material constant and can be obtained in reference (Meyers 1994). e represents the energy per unit mass, α is the porosity with an initial value of α0; ρ is the granite density with an initial value of ρ0, μ denotes the volumetric strain, and A1, A2, and A3 are the polynomial coefficients which can be obtained using the following equations (Xie et al. 2017).

$$A_{1} = \alpha_{0} \rho_{0} c^{2}$$
(6a)
$$A_{2} = \alpha_{0} \rho_{0} c^{2} (2k - 1)$$
(6b)
$$A_{3} = \alpha {}_{0}\rho_{0} c^{2} [(3k - 1)(k - 1)]$$
(6c)

where c denotes the longitudinal-wave velocity. So, A1 = 47.49 GPa, A2 = 57.94 GPa and A3 = 12.17 GPa can be determined. Added to that, when the material is under tension, A1 = T1 and T2 = 0 where T1 and T2 are also polynomial EOS parameters. Furthermore, it is reported by Li (2014a) that the lode angle dependence factor B and the tensile yield surface parameter gt* have almost no influence on the simulation results, and B = 0.0105 and gt* = 0.7 are set according to his recommendation. In addition, the remained parameters use the defaults, and the used parameters for granite are presented in Table 3.

Table 3 Parameters for granite
Full size table

2.2.2 Material model for PETN

To model PETN detonation, the combined High_Explosive_Burn and Jones-Wilkins-Lee EOS is employed, which can predict the process of explosive burning with high fidelity and is widely used in rock blasting modelling. Especially, this EOS can precisely express the pressure during the expansion of explosion products Pe as (Lee et al. 1968)

$$P_{e} = A\left( {1 - \frac{\omega }{{R_{1} V}}} \right)e^{{ - R_{1} V}} + B\left( {1 - \frac{\omega }{{R_{2} V}}} \right)e^{{ - R_{2} V}} + \frac{{\omega E_{e} }}{V}$$
(7)

where A, B, R1, R2 and ω are constants for PETN, which can be determined with regression equations according to the EOS data of PETN with various densities (Banadaki 2010). Ee denotes the detonation energy with an initial energy of Ee0, and Ee0 can be calculated by explosion heats, Ee0 = 5.51 × 109 J/m3. Additionally, the Chapman-Jouguet pressure (PCJ) of explosive is determined by the regression equation according to the pressure measurements of PETN (Green and Lee 2006). The parameters for PETN are given in Table 4.

Table 4 Parameters for PETN explosive
Full size table

2.2.3 Material model for air

The air deck in the blast testing is modelled by combing Mat_Null and Linear_Polynomial EOS. This pressure in air Pa under dynamic loading is expressed as (LSTC 2015)

$$P_{a} = C_{0}^{{}} + C_{1} \mu + C_{2} \mu_{2} + C_{3} \mu_{3} + (C_{4} + C_{5} + C_{6} \mu_{2} )E_{a}$$
(8)

where Ea is the internal energy. In this EOS, the air compression is expressed as μ = (ρa/ρa0)-1, in which ρa and ρa0 are the present and initial air density, respectively. C0 to C6 are material constants, while C4 and C5 are defined as equal to γ − 1, in which γ denotes the ratio of specific heats. The parameters for air are well calibrated in the previous studies, and thus they are directly determined with the references. The parameters for air are given in Table 5.

Table 5 Parameters for air
Full size table

2.2.4 Material model for sand

The Soil_and_Foam is widely applied in modelling the behaviour of sand, and the deviatoric perfectly plastic yield function ϕ for sand is written as (LSTC 2015)

$$\phi = J_{{2}} - \left( {a_{0} + a_{{1}} P_{s} + a_{{2}} P_{s}^{{2}} } \right)$$
(9)

where J2 is the second invariant, a0, a1 and a2 are the yield function constants. The sand has a density of 1950 kg/m3 in the blast testing. The remaining parameters for sand are well verified in the previous investigations. According to Wang et al. (2013), a0, a1 and a2 are 3.4 × 107, 6.387 × 103, 0.3, respectively, and the shear modulus, bulk modulus for unloading and pressure cutoff for tensile fracture are 63.8, 134 and −0. 005 MPa, respectively. Furthermore, the parameters for sand EOS were measured by Xiong et al. (2010) and listed in Table 6. The Vs and Vs0 are the present volume and initial volume, respectively.

Table 6 Parameters for EOS of sand (Xiong et al. 2010)
Full size table

2.3 Numerical verification

2.3.1 3D finite element model for blast tests

Finite element models that have the same dimensions and material as the samples S6, S7 and S8 are constructed according to Table 1, as shown in Fig. 2. The electrical detonator is not built because of its low explosive amount. In these models, common nodes are shared at the interfaces between the explosive, air, sand and rock. The Arbitrary-Lagrange-Euler (ALE) formulation is employed for Parts of explosive, air and sand to calculate the large deformation in the vicinity of borehole whereas the Lagrangian algorithm is applied for rock Part. Moreover, the multi-material ALE algorithm is invoked to model the material mixing among explosive, air and sand.

Fig. 2
figure 2

Numerical model for blast test modelling

Full size image

Before running the 3D models, a mesh sensitivity test is implemented using various mesh sizes for modelling the blasting test-S6, and the difference in simulated damage pattern with various mesh sizes is checked and compared. The mesh size sensitivity test results indicate that the simulated damage cracks are nearly the same with varying element sizes within 1–2 mm whereas the model with a finer element size produces more small fractures and takes much more calculation time, as presented in Fig. 3. At last, considering the modelling reliability and calculation time, the element size of 2 × 2 × 2/mm is used, and these 3D models are discretised into about 1,800,000 hexahedral elements.

Fig. 3
figure 3

Simulated damage patterns after detonation 0.4 ms with different mesh sizes (T—computational time using MPP solver version-810 with 28 cores)

Full size image

2.3.2 Determining damage threshold for rock cracking

In the present modelling, no element erosion or element deletion method is employed because the element erosion/deletion is an energetically non-conservative method that deletes rock element when it still has a residual strength and may thus result in transmission distortion of explosion stress waves. To model rock fracturing, an alternative method using post-processing in LS-PREPOST is adopted. The crack pattern is determined by blanking the damaged rock elements whose damage level D exceeds the critical value of generating rock cracks.

The key of this manner is to determine the damage threshold of crack formation by comparing the simulated rock cracking and fracture network in the physical test. Herein, the simulated crack patterns of S6 generated by blanking damaged elements with D = 0.2, 0.3, 0.4, 0.5 and 0.6 are shown in Fig. 4 (Specific operation in PRE-POST: Fcomp → Misc → History variable 4 → Range → Show → (0–0.4)), and the corresponding rock cracking and fragment morphology in physical testing are presented in Fig. 5a, b, respectively. It can be seen that the simulated main cracks develop more realistic when D ≤ 0.4, whereas many intermittent cracks are formed when D > 0.4. Meanwhile, too large cracks and too few fractures are created, when D < 0.4 and D > 0.4, respectively, compared to the physical testing. Thus, a threshold of D = 0.4 is determined for the formation of cracks. Following, the fracturing networks in simulation for S6, S7 and S8 are shown through blanking the damaged rock elements D exceeding 0.4 in LS-PREPOST.

Fig. 4
figure 4

Crack patterns of S6 with different damage levels

Full size image
Fig. 5
figure 5

Results of blast testing: a rock fracturing during blasting, b rock fragmentation after blasting

Full size image

2.3.3 Comparing results from blast testing and numerical modelling

It can be found in the results of physical testing (Fig. 5) that after blasting, the rock samples were disintegrated into a large number of fragments of various sizes. The high-speed photography results in Fig. 5a show that many relatively large fragments were produced on the upper part of the sample. Besides, the fragment morphology in Fig. 5b suggests that most of the relatively large fragments had an elongated shape, which implies that the formation of these large fragments was mainly controlled by the development of radial cracks. Furthermore, the amount of large fragments (individually placed) and fine particles (sieved and piled up) increases and decreases, respectively, in a sequence of S7, S8 and S6, indicating the blast-induced rock fragmentation performing worse in this order.

Figure 6 shows the evolution process of simulated rock fracturing, and it can be found that an ellipsoidal damage zone is instantly created after initiating the PETN. Then, the damage zone expands and the radial cracks radiate outward. During this stage, the spalling damage (tensile crack) is induced in the vicinity of the bottom surface and at the below part of cylindrical surface since compressive explosion stress waves reflect and transfer into tensile stress waves at these surfaces. Next, the blast-induced damage extends upwards, i.e., both the radial cracks and spalling cracks develop upward, and the radial cracks finally reach the top surface. As a consequence, the radial cracks and spalling cracks interconnect and disintegrate the rock sample, creating plenty of rock fragments of different sizes. The rock fragments in the upper part of sample are larger than those in the below part. It can be found by comparing the rock fracturing in simulation (Fig. 6) and blast testing (Fig. 5a) that a similar fracturing pattern is obtained when using numerical and experimental methods. After blast loading, the breakage extent in the simulation of S7 is larger than that of S6 and S8 because of more explosive. Meanwhile, the fracturing degree in S8 is slightly larger compared with that in S6 due to the shorter air deck and longer stemming, which agrees well with the physical testing.

Fig. 6
figure 6figure 6

Simulated damage crack evolution under blasting of S6, S7 and S8

Full size image

The final damage patterns with mesh and fracture networks in numerical modelling of samples S6, S7 and S8 are presented in Fig. 7. As can be found, the number and size of the relatively large fragments in the simulation are approximate their counterparts in blast testing (Fig. 5b). Then, the blast-induced FSD in simulation is obtained by image processing 5 uniformly distributed cut surfaces at x = 0 mm, y = 0 mm, z = 75, 150 and 225 mm of the model using ImageJ (Specific operation in PRE-POST: Splane → Input coordinates and vectors → Cut → Acen → Screenshot). During the image processing, the borders of blast-induced fragments are automatically identified using an edge-detection algorithm, and in this way, the simulated fragments are delineated. The fragment profile is subsequently fitted with an ellipse, and the fragment size data such as the long and short axes of the fitted ellipse can be determined. The fragment size data of these cut sections are merged and considered as the results of overall fragmenting in numerical modelling. A comparison of FSDs determined through combining finite element modelling and image processing and from blast test is presented in Fig. 8. It should be noted that no parametric study (no parameter adjusting) for constitutive models is involved in numerical calibration.

Fig. 7
figure 7

Simulated rock breakage: a simulated damage pattern, b cracking patterns using blanking damaged rock elements

Full size image
Fig. 8
figure 8

FSD comparison between blast testing and numerical modelling

Full size image

It can be observed in Fig. 8 that good agreements between simulation and blast testing are obtained such as the overall fragment size increasing in a sequence of S7, S8 and S6 and the similar size of large fragments (> 100 mm) in FSDs. Note that there is a small deviation in FSD between numerical and experimental results, especially when the fragment size is smaller than 30 mm. This deviation might be induced by the pre-existing microcracks or discontinuities in the granite which is activated and developed under instantaneous blast loading, resulting in more and finer particles in physical tests. Unfortunately, the pre-existing microcracks in the physical specimens are impossible to be built and modelled on account of the limitation of element size. However, although some differences exist between the blasts of the numerical model and physical sample, the good agreements of blast-induced fracturing and FSD indicate that the currently developed numerical model is reasonable to predict the rock fragmenting induced by blasting with various combinations of explosive, air deck and sand-stemming.

3 Numerical modelling of blasting with different stemming structures

3.1 Computational model

A cubic single-hole model, with a side length of 2 m and 100-mm-diameter borehole centrally drilled thoroughly in the model, is built to numerically investigate the rock fragmenting in blasting with distinct stemming structures. There are two charge modes considered with the same amount of explosive in those models, as shown in Fig. 9a. The first is the concentrated charge, i.e., the explosive is located in the model centre and sandwiched by stemming with an air deck (similar to the stemming structure in blast testing and defined as top stemming) or this explosive is directly sandwiched by stemming (in the following section, it is called as bottom stemming). When stemming locates above the explosive, the bottom stemming can be achieved by gravity, whereas it can be achieved by an air column (bamboo commonly used in practice) placed at the explosive bottom when stemming locates below the explosive. The second is the deck charge, i.e., two cartridges placed in the hole and separated by a central air or sand deck while the two ends of borehole are stemmed. When using the first charge mode, the effect of stemming ratio (defined as the stemming length over the uncharged length in a hole) on rock fragmentation is numerically analysed with the concentrated charge model as presented in Fig. 9b. Using the second charge mode, the influence of air/sand deck length on rock disintegration is numerically investigated with deck charge model as presented in Fig. 9c. Considering the rule-of-thumb with respect to specific charge in controlled blasting which involves various stemming structures, a 0.24-m-long cartridge is constructed when adopting concentrated charge. In the subsequent modelling, the explosive is ignited at its centre. The detailed parameters of the computational models are tabulated in Table 7.

Fig. 9
figure 9

Computational model: a sketch of different stemming structures in a hole, b concentrated charge model, c deck charge model

Full size image
Table 7 Detailed parameters of computational models
Full size table

To balance the computational accuracy and time, an element size of 20 × 20 × 20/mm is adopted after a mesh convergence test, and the model is discretised into 996 300 hexahedral elements. Meanwhile, common nodes are also used at the borders of the Parts of explosive, air, sand and rock, and the same algorithms as the blast test modelling are applied. The simulation concludes when there are no further extensions of existing cracks and no new cracks are formed, and the simulation duration is 2.0 ms. The modelling is performed using an MPP solver (810 version) on a system with 32 cores (28 cores are invoked) and 256 GB of RAM.

3.2 Modelling of rock blasting with various stemming ratios

Using the verified constitutive models, the blast-induced rock fragmenting with concentrated charge and various stemming ratios is first modelled. Herein, two stemming modes, i.e., top-stemming and bottom-stemming are involved, which are both widely adopted in controlled blasting. The explosion pressure in rock and the evolution processes of damage crack induced by blasting without stemming, 40% top stemming, 40% bottom stemming and full stemming are presented in Fig. 10. As can be observed, the blast-induced explosion pressure spherically radiates outwards after the detonation of concentrated charge, and an ellipsoid damage zone is simultaneously formed due to this high-intensity loading. The blast-induced compressive stress waves then arrive at free surfaces and shift into tensile stress waves, and thus spalling damage (crack) generates in the vicinity of surfaces. Thereafter, the outward propagating radial cracks gradually extend to free surfaces and interact with spalling cracks, and they eventually coalesce with each other, forming the final crack pattern and disintegrating rock into various sizes of fragments.

Fig. 10
figure 10

Simulated explosion pressure and crack evolution in blasting using different stemming ratios

Full size image

Besides, it can be observed that in the cases of bottom and full stemming, the high-pressure zone around borehole is much larger than those of top and no stemming. Due to this pressure difference, the damage extent of blasted rock with top and no stemmed boreholes is smaller. In addition, the rock damage extent increases with the increase of stemming ratio, indicating that stemming length plays a role in promoting blast-induced rock fracturing.

To reveal the mechanisms mechanisms responsible for variation in blast-produced rock damage with distinct stemming structures and lengths shown in Fig. 10, the time-history curves of explosion pressure at two monitoring points (in borehole and rock) are recorded during the simulation and presented in Fig. 11. It's worth noting that, in the case of top and no stemming, the monitored point is located in the air deck, whereas in the bottom and full stemming cases, it is located in the sand. As can be seen in Fig. 11a, the explosion pressure in sand (bottom or full stemming) is much higher than that in air (top or no stemming). Similar observations can be found in Fig. 11b that the pressure in rock with bottom and full stemming is larger than those with top and no stemming, indicating that the explosion-induced pressure is essentially dependent upon the stemming mode. This is because the difference in wave impedance between the explosive and air is extremely high whereas it is much lower between explosive and sand. This impedance difference allows higher explosion pressure transmitted into sand, leading to higher explosion pressure transmitting into rock in the blasting using bottom and full stemming.

Fig. 11
figure 11

Time-history curves of detonation pressure with different stemming ratios: a borehole pressure at z = 1.22, x = 0, y = 0, i.e., 0.1 m up to explosive, b pressure in rock at z = 1.0, x = 0.5, y = 0

Full size image

Besides, in Fig. 11, the same value of maximum pressure can be found at monitoring points in blasting with bottom and full stemming, and top and no stemming, respectively. This implies that the variation of stemming ratio does not influence the pressure peak induced by explosive detonation. With the same stemming mode, the pressure difference in blasthole can be found in pressure tail only. Accordingly, the variation in blast-induced pressure tail caused by changes in stemming ratio leads to a difference in blast-induced rock cracking. Full stemming leads to a higher pressure tail in the hole than 40% bottom stemming because of the best confinement. Meanwhile, a larger explosion pressure tail is induced by top stemming than no stemming owing to the reflection of explosion stress waves. With an air deck (top stemming), reflected detonation stress waves apparently pass through the monitoring point in hole twice, which is caused by the stress wave reflections at the interfaces of air and sand stemming near and away from the monitoring point. Consequently, the differences in pressure peak and pressure tail, arising from different stemming modes and various stemming lengths, respectively, induce the discrepancy in rock damage observed in Fig. 10.

The final damage cracks and fracturing patterns in blasting using various stemming ratios with top and bottom stemming are shown in Fig. 12a, b, respectively. In consistency with Figs. 10 and 11, The extent of rock fracturing induced by blasting with bottom stemming is greater than that with top stemming, and it increases with an increase in stemming ratio. This clearly indicates that changing from top stemming to bottom stemming and increasing the stemming length are significant factors in improving blast-induced rock fracturing. Using bottom stemming and longer stemming lengths increases the potential for the formation of smaller rock fragments.

Fig. 12
figure 12

Damage cracks and fracturing patterns in blasting using various stemming ratios

Full size image

3.3 Modelling of rock blasting with various deck lengths

Then, the rock fragmenting in blasting using air or sand deck is simulated with different deck lengths (the stemming length also changes accordingly). The simulated explosion pressure and damage pattern with air/sand deck lengths of 0.4 and 1.0 m are presented in Fig. 13. When using the deck charge, the blast-induced damage expands from different explosive positions and interacts, which is dominated by the colliding and interaction of detonation stress waves. Under this condition, the distance and medium between two cartridges have a great influence on the distribution of explosion stress, further dominating the performance of blast-induced rock fragmentation.

Fig. 13
figure 13

Simulated detonation pressure and damage crack evolution with various deck lengths

Full size image

It can be seen in Fig. 13 that the high-pressure zone around the borehole with a sand deck is remarkably larger than that with an air deck, and more rock fractures are generated in blasting using a sand deck. This is also attributed to the difference in wave impedance between sand and explosive smaller than that between air and explosive. More explosion energy travels through rock in blasting with sand deck charge than air deck one. In addition, the high-pressure zone induced by blasting becomes smaller as the deck length increases due to the rapid attenuation of explosion stress waves in air and sand. This ultimately leads to the formation of fewer cracks and less damage.

The time-history curves of detonation pressure within borehole and in rock with various deck lengths are plotted in Fig. 14. It can be observed that the maximum pressures in the borehole and rock when using a sand deck in blasting are significantly larger than those when using an air deck. This pressure peak decreases with an increase in deck length, as explained earlier, due to the impedance difference among explosive, air, and sand and the rapid decay of explosion pressure in sand and air. Besides, in Fig. 14a, The pressure in the borehole rises at different times due to the difference in velocity of detonation stress waves transmitting in sand and air, as well as varying deck lengths. A similar reason accounts for the discrepancy in the rising time of the explosion-induced pressure in rock (Fig. 14b). After the explosion stress waves collide at the monitoring points, the explosion pressure within the borehole oscillates for a period. This oscillation is more noticeable in blasting with an air deck because the oscillating pressure in the air deck hole is closer to its pressure peak. However, referring to Fig. 13, it should be pointed out that the rock damage in blasting using deck charge is mainly controlled by the pressure peak rather than this pressure oscillation in the hole.

Fig. 14
figure 14

Time-history curves of detonation pressure with various deck lengths: a pressure at the midpoint of the centre-to-centre line of two charges (z = 1.0, x = 0, y = 0), b pressure in rock at z = 1.0, x = 0.5, y = 0

Full size image

The simulated damage patterns and fracture networks with different air/sand deck lengths are shown in Fig. 15. It is evident that the extent of cracking induced by blasting with a sand deck is greater than that with an air deck, particularly with a relatively large deck length. Additionally, the fracturing extent decreases as the deck length increases. This implies that better rock fragmentation performance can be achieved by switching from an air deck to a sand deck and reducing the deck length. These results are to be expected since a reduction in deck length implies a longer stemming length, which consequently leads to the creation of more blast-induced rock fractures, as demonstrated in Sect. 3.2. Besides, the sand delivers more explosion energy into the system of explosive, stemming and rock than air does.

Fig. 15
figure 15

Damage cracks and fracturing patterns in blasting using various air/sand deck lengths

Full size image

4 Discussions

The current simulation has demonstrated that the rock fragmentation performance improves with longer stemming, shorter air/sand deck length, altering top stemming to bottom stemming, and changing air deck to sand deck. In this section, the quantitative discussions associated with the effect of the combination of explosive, air/sand deck and stemming on blast-induced rock fragmentation are given as follows.

After obtaining the blast-induced fracture networks, the calculated model is uniformly cut at 9 planes (x = −1, 0, 1, y = −1, 0, 1, z = 0.5, 1, 1.5 m), and the cracking patterns in cut surfaces are image-processed in ImageJ and merged to obtain the results of overall fragmenting. Following, the fragment size data of blast-indued rock fragmentation are characterized using a three-parameter generalized extreme value function (Hogan et al. 2012)

$$F(d;\mu ,\sigma ,\xi ) = \exp \left\{ { - \left[ {1 + \xi \left( {\frac{d - \mu }{\sigma }} \right)} \right]^{ - 1/\xi } } \right\}$$
(10)

where μ denotes the location parameter which is commonly considered as the average fragment size, σ refers to the scale parameter concerning the FSD range, and ξ is a parameter about fragment shape. This function was reported that it is suitable for characterizing the rock FSD data generated by dynamic load (Hou et al. 2015; Shen et al. 2017; Tao et al. 2020).

The FSDs and corresponding fitting curves under blasting with various stemming ratios are given in Fig. 16. As can be observed, the fit curve rises earlier with the increase of stemming ratio. In the pictures, the fitting curve rising earlier means an increasing amount of small fragments, which further leads to a smaller overall fragment size with a broader range of FSD. So, the fragment size decreases and the FSD span broadens as the stemming length increases. In the case of top stemming, the FSD changes gradually with the increase of stemming ratio, and thus the overall fragment size produced by rock blasting has a relatively large difference with various stemming lengths. In contrast, a dramatic change happens in FSD produced by blasting when using bottom stemming replacing no stemming. In other words, the use of bottom stemming can significantly improve the performance of rock fragmentation, even though only 20% stemming is applied. Moreover, when the length of bottom stemming continuously increases, the FSDs have a small difference. Therefore, in practical blasting, the bottom stemming is recommended since it is efficient in intensifying rock fragmentation and the blast performance is robust to stemming length, which can significantly decrease fragment size in rock blasting with the simple and highly fault-tolerant operation.

Fig. 16
figure 16

FSDs and corresponding fitting curves (R2 > 0.94) under blasting with various stemming ratios: a Top stemming, b bottom stemming

Full size image

Moreover, the FSDs and corresponding fitting curves under blasting with different air/sand deck lengths are plotted in Fig. 17. It can be found that the fitting curves rise later as the deck length increases, while they rise earlier in blasting with a sand deck. These observations are in agreement with the above findings that the blast-produced rock fragmenting with sand deck performs better and the rock fragmentation performance worsens with the increase of air/sand deck length.

Fig. 17
figure 17

FSDs and corresponding fitting curves (R2 > 0.95) under blasting with various deck lengths: a air deck, b sand deck

Full size image

The parameter variations in μ (regarded as the average fragment size) and σ (parameter of FSD range) with the change of stemming ratio and deck length are presented in Figs. 18 and 19, respectively. As can be observed, the parameter of average fragment size decreases and increases, respectively, with the increase in stemming ratio and deck length, while the FSD range parameter conversely varies with the change in length of stemming and air/sand deck. Meanwhile, the fragment size parameter and FSD range parameter in top-stemming blasting are larger and smaller, respectively, compared with the counterparts under bottom-stemming one, and they are higher and lower in blasting with air deck than those with sand deck, respectively, which again indicates that blasting with top stemming yields coarser fragments and narrower FSD span compared with that using bottom stemming whereas blasting employing sand deck creates finer fragments and broader FSD range than air deck. It is noteworthy that the average fragment size in blasting with 20% bottom stemming and full stemming is reduced by 41% and 49%, respectively, compared to blasting without stemming. The average fragment size in full stemming blasting is only about half of that in blasting without stemming. As a result of these changes in average fragment size, the span of FSD significantly increases when changing from no stemming to bottom stemming. In blasting with 20% bottom stemming and full stemming, the FSD range parameter increases by 26% and 31%, respectively, when compared to the case of no stemming.

Fig. 18
figure 18

Variations of fragment parameter μ with various a stemming ratios and b deck lengths

Full size image
Fig. 19
figure 19

Variations of fragment parameter σ with various a stemming ratios and b deck lengths

Full size image

Besides, it can be seen that the average fragment size under blasting using top stemming shows a linear reducing relationship with the increase of stemming ratio, and it performs a linear increasing correlation with the increase of air/sand deck length under blasting using deck charge. Similarly, the range of FSD with top stemming shows a linear increasing tendency with the stemming ratio, and meanwhile, it almost linearly decreases as the deck length increases. The opposite trend between fragment size and FSD span with the variation of stemming ratio and deck length can be explained as follows: when transmitting higher explosion energy into the rock by increasing stemming length, reducing air/sand deck length, replacing top stemming to bottom stemming and altering air deck to sand deck, finer rock fragments are formed, and thus the average fragment size decreases. Under this condition, the discrepancy in fragment size between the large fragments and fine particles increases, generating a broader span of FSD. Conversely, lower explosion energy transmitted into rock results in coarser rock fragmentation and thus leads to a smaller difference in size between big and small fragments which generates a narrower FSD range.

5 Conclusions

This study provides numerical insights into the influence of stemming on blast-produced rock fragmenting. In the current numerical modelling, the explosion pressure and rock fragmentation are numerically investigated with distinct combinations of explosive, air/sand deck and stemming by combining finite element modelling in LS-DYNA and image processing in ImageJ. According to the above findings and corresponding discussions, the subsequent summaries can be obtained:

After determining the parameters of granite, explosive, air and sand based on basic physical and mechanical parameters in blast testing, theoretical calculations and references, the good agreement between blast testing and simulation suggests that the currently developed numerical model is reasonable to model rock fracturing and rock fragmenting in blasting using different combinations of explosive, air deck and stemming. Following the simulation, the crack pattern can be obtained by blanking the rock elements whose damage level is higher than the critical value of crack formation, and the rock fragment size data can be determined using image processing.

The combination of explosive, air/sand deck and stemming plays a critical role in blast-induced rock fragmentation. The size of rock fragments during blasting decreases with an increase of stemming length and a decrease of air/sand deck lengths, changing from top stemming to bottom stemming and altering from air deck to sand deck. Conversely, the FSD range shows an opposite tendency. The average fragment sizes in blasting with 20% bottom stemming and full stemming are reduced by 41% and 49%, respectively, compared to blasting without stemming. The bottom stemming is recommended for practical blasting due to its high efficiency in improving rock fragmentation performance and high robustness related to stemming length. Furthermore, the average fragment size in blasting with top stemming and deck charge (both air and sand deck) presents a linear reduction correlation with an increase in stemming length and a linear increase relationship with an increase in deck length. In comparison, the FSD range in the corresponding case shows an opposite linear tendency with the change of stemming ratio and air/sand deck length.