Experimental testing and numerical simulations of blast-induced fracture of dolomite rock

In this paper, the Johnson-Holmquist II (JH-2) model with parameters for a dolomite rock was used for simulating rock fragmentation. The numerical simulations were followed by experimental tests. Blast holes were drilled in two different samples of the dolomite, and an emulsion high explosive was inserted. The first sample was used to measure acceleration histories, and the cracking pattern was analyzed to perform a detailed study of the blast-induced fracture to validate the proposed method of modelling and to analyze the capability of the JH-2 model for the dolomite. The second sample was used for further validation by scanning the fragments obtained after blasting. The geometries of the fragments were compared with numerical simulations to further validate the proposed method of modelling and the implemented material model. The outcomes are promising, and further study is planned for simulating and optimizing parallel cut-hole blasting.


Introduction
High explosives (HEs) are widely used in civil engineering and the mining industry for controlled fragmentation and rock removal [1,2]. However, dozens of tests are needed to analyze, fully understand and optimize the process, which is costly and timeconsuming. As a result, it is only possible to determine the burden. Available field data are insufficient, as several elements must be included in the process: number of drilled blast holes, their positions, delay intervals, amount of HE used or behavior and response of the rock itself [3]. Numerical simulations can be implemented to overcome these deficiencies and have proved robust and efficient for modelling almost every problem in the field of engineering [4][5][6][7][8][9][10][11].
Regardless of the method used to simulate the problem, the material behavior needs to be reproduced reliably in order to ensure that the results are correct and comparable to actual data. Therefore, a constitutive model that captures damage, cracking and other physical and mechanical properties as comprehensively as possible is needed [15, 34,35]. Several models are available for the response of brittle materials, especially rock. The most commonly used are the Johnson-Holmquist Concrete model (JHC model; sometimes referred to as the Holmquist-Johnson-Cook model) [36,37], the brittle damage model [22,38], the Federal Highway Administration (FHWA) soil model [13,39], the Riedel-Hiermaier-Thoma (RHT) model [3,40], the Continuous Surface Cap (CAP) model [38,41], and the Karagozian and Case Concrete (KCC) model [42,43]. Among constitutive models, the Johnson-Holmquist II (JH-2) model has been extensively used to reproduce the behavior of brittle and geomaterials under dynamic or impact loading conditions [15,18,[44][45][46][47][48][49][50][51].
In this study, the blast fragmentation of a dolomite rock was experimentally and numerically studied. The field tests were reproduced with numerical simulations adopting explicit LS-Dyna code. Two different samples of the dolomite were considered. In each sample, blast holes were drilled, filled with emulsion HE inserted in the bottom of the sample, and covered with stemming. The first sample was used to measure the acceleration histories to conduct a detailed study of the blast-induced fracture of the rock to validate the proposed method of modelling and to analyze the capability of the JH-2 model for the dolomite. Additionally, the damping properties of the rock material were evaluated. The second sample was used to thoroughly analyze fragmentation quantitatively and qualitatively using the scanned fragments. For the numerical analyses, the finite element method (FEM) and SPH were used to represent the dolomite rock and the HE with stemming, respectively.
2 Experimental procedure of the fragmentation test Equation Section (Next) The experimental tests presented in this paper are part of a wider project and only the results obtained from the two samples of dolomite are presented and discussed here. Sample no. 1 had a volume of approximately 0.364 m 3 and a mass of * 965.0 kg ( Fig. 1), whereas sample no. 2 had a volume of 0.472 m 3 and a mass of * 1240.0 kg (Fig. 2). In both samples, a blast hole with a diameter of 30.0 mm and depth of 430.0 mm was drilled near the center. In the bottom of the hole, an emulsion HE was inserted with a mass of 10.0 g. Emulsion explosives are the second most commonly used explosive type worldwide [52] and are commonly used in Polish underground mine operations [1,53]. Some of the physical properties of the HE are presented in Table 2. The blast hole was filled with stemming. In all cases, the blasting process was recorded using a Go-Pro camera. The first test was conducted to validate the proposed methodology of modelling as well as the JH-2 constitutive model for the dolomite. For this purpose, vertical acceleration histories were recorded using PCB piezoelectric gauges with a measurement range of up to 100,000 g. The gauges were attached to a steel ring, which was fixed to the rock sample using glue. Both sensors were placed in line from the blast hole: sensor A at a distance of 330 mm and sensor B at a distance of 530 mm from the hole. The sensors were connected to a SIRIUS HS measuring amplifier, and the cracking characteristics were analyzed and compared with numerical simulations. To further validate the modelling method, the second test was conducted. Fragments of rock sample no. 2 were scanned using a 3D laser scanner, and the field results were compared with the outcomes of finite element analysis (FEA).
3 Methodology for numerical simulation of the fragmentation test

Numerical integration procedure
Numerical studies of the blast-induced fracture and fragmentation of the dolomite were carried out using FEM and SPH with an explicit integration procedure. Massively parallel processing (MPP) LS-Dyna code, which uses explicit central difference time integration, was adopted. The dynamic equation solved during the calculations has the following form [54]: where M is the diagonal mass matrix; F ext n is the external and body forces; F int n is the stress divergence vector; C is the damping matrix; and x, _ x, € x are the nodal displacement, velocity and acceleration vectors.
Assuming that _ x % _ x nÀ1=2 , Eq. (3.1) is solved with numerical integration of acceleration € x n : Implementation of the central difference equations for velocity and displacement yields the following relations: The main disadvantage of this method is that it is conditionally stable and requires a time step to be limited according to the Courant-Friedrichs-Lewy (CFL) stability condition [54]: where Q is a function of the viscous coefficients C 0 and C 1 and is formulated as follows: where L E ¼ V E =A Emax is the characteristic length of the element; V E is the volume of the element; A Emax is the largest side of the element area; and c is the adiabatic speed of sound.

SPH method
The SPH method is a mesh-free particle method in which all computational properties, e.g., mass, position, and velocity, are associated with a particle. The description of the SPH continuum uses the principles of conservation of mass, momentum and internal energy as the governing equations in the Lagrangian description applied in FEM [55]: where r is stress; v is velocity; ij are indexes of components; E is internal energy; m is the mass of the particle; W is a kernel function; N is the number of particles within the smoothing length; and p, q are indexes denoting different particles. Equations (3.7) and (3.8) are solved by interpolation of a given function value \ f [ (i.e., density, velocity, energy, etc.) at a given point based on the known value of this function at the surrounding points (particles) using the following formula [55]: ð3:9Þ where\ f [ is the function interpolation; x is the vector defining the particle's position; h is the maximum distance between particles (smoothing length); and W is a kernel function.
In the LS-Dyna code, a cubic B-spline kernel with the following form is implemented [54]: 2ph 3 for the one-, two-and three-dimensional cases, respectively.
The smoothing length varies throughout the simulation. Therefore, the number of neighboring particles remains relatively constant and is obtained by recalculating the smoothing as follows: or by solving the continuity equation: where d 0 is the initial density, h 0 is the initial smoothing length, and p is the dimension of space.
4 Modelling and simulation of the fragmentation testEquation Section (Next)

Model definition
The geometries of the three rock fragments were achieved based on the 3D laser scanning process from which the point cloud was obtained. As a result, discrete models were developed consisting of mainly 1-point nodal tetrahedron elements; however, brick elements with one integration point were also used within the area of the blast hole. In addition, SPH particles were used to describe the HE and stemming. During parallel studies, this FEM-SPH coupled model proved to be sufficiently effective and gave identical results for the test case compared with ALE modelling, which has been widely implemented for analyzing blast wave interactions with different types of structures [27,56,57]. The main advantage of the meshless approach is the significantly shorter computational time. In the LS-Dyna code, two fundamental approaches are available in which particles from different SPH parts can interact with each other. The default strategy is to use particle approximation to compute the inter-part particle interaction. However, when interpolation functions are used, numerical instabilities can arise. To treat SPH parts separately and simulate the interactions between them in parallel, penalty-based node-to-node contact can be introduced. Such an approach was implemented in the present study. On the other hand, a typical node-tosurface contact based on the penalty formulation was used between the SPH parts (HE and stemming) and the dolomite rock [58].
For the discretization of each rock, an average 5.0 mm mesh size was assumed for brick elements and tetragonal elements. This mesh was selected based on parallel numerical simulations of the rock dolomite subjected to blast loading. Ultimately, 20,332,016 and 21,350,000 elements were used to represent sample no. 1 and no. 2 of the dolomite rock, respectively. The models with the initial boundary conditions are presented in Figs. 3 and 4. In the case of rock sample no. 1, the vertical acceleration histories were measured in nodes representing the sensors in the experimental tests.

Rock dolomite
In all cases, the dolomite was modelled using the JH-2 model, which is based on the relation between normalized values of equivalent stress and pressure. The model is described with three surfaces that represent intact, damaged and fractured states of the material (Fig. 5). The normalized intact strength of the material presented in Fig. 5 is described using the following formula [44,45]: where rÃ I ¼ r I =r HEL is the normalized intact strength (r I is the current equivalent stress, and r HEL is the equivalent stress at the Hugoniot elastic limit (HEL)); A, N are intact material constants; C is the strain rate coefficient; PÃ ¼ P=P HEL is the normalized hydrostatic pressure (P is the current hydrostatic pressure, and P HEL is the pressure at the HEL); TÃ ¼ T=P HEL is the normalized maximum tensile hydrostatic pressure; and _ eÃ ¼ _ e= _ e 0 is the dimensionless strain rate ( _ e is the current equivalent strain rate, and _ e 0 ¼ 1:0 s À1 is the reference strain rate).
The damaged material is represented with a dashed line in Fig. 5 and is determined by [44,45]: where D is a damage factor and takes a value between 0 and 1. The normalized fractured strength of the material is represented by rÃ F and is described using the following formula [44,45]: where B, M are the fractured material constants. In Fig. 5, rÃ FMax , which represents the maximum value of rÃ F , can be observed. This parameter makes it possible to control the upper limit of the fracture strength. The formulas above characterize the behavior of the material described with the JH-2 constitutive model. When the value of equivalent stress is larger than rÃ I , the material starts to deform plastically (Eq. 4.1). Then, damage to the material starts to accumulate, and its strength decreases until it reaches the damage surface rÃ D (Eq. 4.2). In this state, the material is partially damaged (0 B D B 1). The softening continues until the material is fully damaged, which means that the value of damage parameter D equals 1.0 and the material is characterized by the fractured material surface (Eq. 4.3). The damage is calculated as follows [44,45]: where De P is the increment of the equivalent plastic strain during a calculation cycle and e F P is the equivalent fracture plastic strain given by e F P ¼ D 1 P Ã þTÃ ð Þ D 2 (D 1 , D 2 -damage constants). The procedure for determining JH-2 constants is not a straightforward task and is not the main topic of the study; thus, only a brief discussion is included. Please see [59] for a detailed description of the procedure for determining JH-2 constants. In Table 1, the JH-2 parameters for the dolomite are presented. In the present paper, a different dolomite was studied than in the authors' previous paper [27], and therefore some of the constants were modified.

Stemming and high-explosive material models
The HE was described using the Jones-Wilkins-Lee (JWL) equation of state and the High Explosive Burn constitutive model (HEB). The equation describing the behavior of detonation products has the following form [58]: where V = q 0 /q; q 0 is the initial density of the HE; q is the actual density of the HE; E 0 is the detonation energy per unit volume and the initial value of E of the HE; and A HE , B HE , R 1 , R 2 , and x are empirical constants determined for a specific type of explosive material based on experiments [60] using the Gurney energy, detonation pressure and explosion heat. All required parameters of the HE used in the mining industry were defined based on the results of a so-called cylindrical test (Table 2). A detailed description of the test procedure can be found in [60]. The stemming was simulated using the FHWA soil model with the reported parameters of unsaturated sand (Table 3) [61] and supplemented with values from the literature [13,39]. The FHWA soil model is efficient for modeling the behavior of soil or/and sand with the consideration of strain softening, strain rate effects, moisture content and kinematic hardening [13,39]. Moreover, the FHWA material model is stable even if unconfined conditions occur.

Results and discussion
To investigate rock blasting and fragmentation, two different dolomite samples were considered. The first sample was used to measure acceleration histories for a detailed study of blast-induced fracture in order to validate the proposed modeling method and analyze the capability of the JH-2 model for the dolomite. Furthermore, damping properties were evaluated to confirm the proper dynamic response of the material. The second sample was used to analyze the fragmentation characteristics in detail using the scanned fragments to further validate the proposed method of modelling, the JH-2 material model and the correlated damping coefficients.

Rock sample no. 1
The main objective of the first stage was to compare the dynamic response of the rock material with the measurements from the field tests. Figure 6 presents the vertical acceleration histories for sensors A and B. The analysis of the results revealed that sensor A (closer to the blast hole), together with the steel ring, was detached from the rock mass immediately after the detonation of the HE (it was found next to the rock sample). Therefore, the gauge did not transfer any acceleration, and the results were rejected for further analysis. However, the signal from acceleration sensor B (farther from the blast hole) had a damped vibration characteristic. The first peaks, which occurred for 75 ls after detonation, were interference and were not taken into consideration for further study. Ultimately, the acceleration history from sensor B was used for comparison with the outcomes of the numerical simulations. However, before the validation analyses, the damping properties of the tested dolomite were calibrated.

Calibration of damping properties
For dynamic analysis of rock blasting, the damping mechanism and its magnitude are crucial for obtaining reliable outcomes with a correct material response [62]. The damping in the numerical simulations should represent the energy loss in the system during blast wave propagation [63]. Given the mentioned phenomena, the Rayleigh classical approach, which is commonly used in FE dynamic analyses, was introduced in the present numerical simulations of blastinduced fracture of the dolomite [64,65]. According to Rayleigh theory, the damping matrix is a linear ombination of the mass and stiffness matrixes: where C is the damping matrix, M is the diagonal mass matrix; K is the stiffness matrix; a is the damping coefficient for mass-weighted damping; and b is the damping coefficient for stiffness-weighted damping.
Determining the values of damping coefficients is quite difficult for rock-like and geomaterials due to the uncertain boundary conditions, especially initial cracks, heterogeneity and residual stress. Therefore, in the present paper, several coefficients of stiffnessweighted damping were tested for calibration purposes. Values of b of 0.05, 0.1, 0.15, 0.20 and 0.25 were taken into consideration, and the obtained vertical accelerations for each case were compared with the actual outcomes (Fig. 7). The results show that the stiffness-weighted damping coefficient significantly influences the acceleration oscillations and energy loss caused by damping.
For a deeper analysis of the effect of damping on the obtained results, the Fast Fourier Transform (FFT) was used to quantitatively compare all waveforms (Fig. 8). The same sampling frequency and window size were used for each characteristic, including the experimental one. To assess feasible values of the stiffness-weighted damping coefficient, the maximum FFT magnitude for each b was taken from Fig. 8 and compared with the actual outcomes (Fig. 9). The magnitudes were selected for frequencies less than 20 kHz. However, between 20 and 30 kHz, another dominant wave frequency occurred, which is clearly seen in the close-up view included in Fig. 8. This wave frequency was attributed to shock wave reflection from the ground on which the rocks were placed. In the numerical simulations, the ground was omitted. Moreover, the gauges, steel rings and glue were not modelled, which could also affect the differences in the obtained waveforms. Ultimately, a satisfactory correlation of the maximum amplitude was obtained for b equal to 0.20, and this value was used for the validation and further numerical simulations. The mass-weighted damping a influences lower frequencies, and its impact on the results was insignificant in the cases in the present study. Consequently, the a coefficient was set to 1.0 in all FEA.   Fig. 10, the vertical acceleration vs. time characteristics obtained from the field test and the FEM-SPH simulation with b equals to 0.20 are presented and compared. The characteristics obtained from acceleration sensor B are reproduced quite well: the oscillations begin nearly at the same time as in the experimental test, and their duration is also similar (energy dissipation is comparable). In the numerical simulations, several simplifications were considered: isotropic, homogenic properties of the dolomite were used, and the ground and the sensors with steel rings were modelled.
An example of the pressure distribution within the rock for selected time moments is presented in Fig. 11. For a better visualization of the outcomes, the fringe was limited to -25.0 MPa and ?25.0 MPa, and the model was cut. During the initial moments, the shock wave travelled into the rock, resulting in compressive damage in the area close to the blast hole. Next, the shock wave travelled back and forth and was reflected from the outer surfaces of the rock.
The reproduction of the actual test was qualitatively analyzed based on the comparison of the characteristics of the cracks and fragments in the experimental test with those from the numerical simulations. Figure 12 presents the results of the simulation and actual test at the selected stage of the blasting process. The cracking characteristic was satisfyingly reproduced in the numerical modelling. The sudden expansion of the HE products resulted in crack formation. Although the major damage was caused by radial cracks, transversal cracking was also observed. Ultimately, dolomite rock sample no. 1 was divided into five main pieces. However, some small pieces of rock were observed in the area directly under the blast hole.

Rock sample no. 2
In the next stage of the study, the proposed method of modelling, the implemented JH-2 constitutive model, and the correlated damping properties were validated based on the blasting fragmentation test conducted with the second rock sample. The crack pattern is shown in Fig. 13. The cracks are represented in red color, which corresponds to a damage index of 1.0; this value indicates that the material is fully damaged. The results are presented for 0.01 s of simulation time. The rock was fragmented into two main pieces, one of which was smaller than the other. In the field test, the formation of cracks resulted in a loss of material continuity. This phenomenon was not reflected in the numerical simulations because an erosion criterion was not used (numerical method for removing finite elements from the model). To compare the geometries of the fragments obtained in the numerical simulations and experimental tests, 3D laser scanning was used. The resulting point cloud was used to calculate the volumes of the real-world fragments and the fragments from the numerical simulations. Figure 14 presents the numerical simulations and field tests of rock sample no. 2 after blasting. To better analyze the results, the fully damaged finite elements were blanked. Moreover, fragment no. 2 was manually Based on 3D laser scanning, the volumes of the obtained fragments were calculated and compared with their discrete representatives ( Table 4). The volumes of fragment no. 1 obtained from the numerical simulations and actual test were nearly identical. However, higher errors were obtained for larger fragments. In general, the field results were satisfactorily reproduced. The discrepancies are attributable to the following factors: • the smaller fragments visible in Fig. 14 were not measured. • The fully damaged finite elements, which were blanked for volume calculations, were not taken into consideration. • The fragments within the area under the HE were treated as parts of fragments no. 1 and no. 2.

Conclusions
This paper presents experimental and numerical studies of blast fragmentation of two samples of a Based on the outcomes, the following conclusions can be drawn: • Knowledge of the damping mechanism is crucial in order to obtain reliable results from numerical simulations. In this study, the Rayleigh classical approach, which is commonly used in FE dynamic analyses, was considered. Due to the difficulty of determining damping properties for rock-like materials, several values of coefficients of stiffness-weighted damping were tested. b significantly influenced the acceleration oscillations and energy loss caused by damping. By contrast, the massweighted damping a only influenced lower frequencies and thus did not significantly impact the results in the present analysis. Consequently, a was set to 1.0 in all FEA. • The closest agreement between the acceleration history measured in FEA and the maximum FFT magnitude was obtained at b equal to 0.20. The numerical simulations satisfactorily reproduced the cracking characteristics and fragmentation of both rock samples. Furthermore, the fragments from rock sample no. 2 were scanned using a 3D  laser scanner to compare the volumes and geometries of the real-world and simulated fragments. • The implemented JH-2 constitutive model with parameters determined for dolomite proved to be efficient and accurate for simulating the blasting and subsequent fracture and fragmentation of the rock. Nevertheless, further studies are planned to validate the model in penetration and small-scale tests. The influence of the mesh will be also analyzed in future work. • The FEM-SPH coupled method of modelling was sufficiently effective, and the results of FEA were very similar to those of the actual tests. The main advantage of the adopted meshless approach for simulating the HE and stemming is the significantly shorter computational time compared with the ALE approach.
In future studies, the JH-2 model and the FEM-SPH method will be implemented in further numerical simulations of parallel cut-hole blasting, where a proper representation of fracture, cracking and fragmentation is crucial. Moreover, the investigations will be extended to optimize the cut-hole pattern and increase the effectiveness of the blasting method and rock fragmentation. These analyses will implement a methodology in which the fragment area is calculated based on 2D sections of the FE model and the Swebrec function is employed to fit the fragment area distribution [19,58].