In this paper, a multiscale modelling and simulation of destress blasting of rock mass is presented. The proposed and novel approach combines two separate 3D solutions: the first was obtained for the small-scale problem, face(s) blasting process, and the second for the global scale problem, seismic wave propagation within very large volumes of surrounding rock mass. Both the approaches were based on explicit dynamic modelling methodology using the central difference method. In the local case, the arbitrary Lagrangian–Eulerian (ALE) procedure with the Jones–Wilkins–Lee (JWL) equation defining an explosive material was used. For this purpose, a selected volume of a rock mass comprising a blasted mining face was modelled in detail. From the numerical simulation, pressure distribution over the modelled rock was obtained, which was the basis for an initial condition for the global 3D FE model. The peak particle velocity (ppv) distribution obtained from finite element analysis was compared with experimental outcomes. A reasonable agreement between both results was achieved; therefore, the adopted multiscale modelling approach confirmed its effectiveness and that it can be successfully implemented in further engineering practice.

1. Introduction

The flat-type copper ore body in Polish deep underground mines is basically excavated with drill-and-blast technology which seems to be relatively well suited to the hardness of the local rocks and to local mining/geological conditions. Besides the direct production potential, blasting works are also recognized as one of the most important active methods of seismic events prevention. This is achieved through strain-release massive simultaneous blasting engaging ten to forty production faces near potentially unstable main roof strata and/or pillars. The strain-release effect is in proportion to the rock mass ability for strain energy accumulation. A significant number of recorded seismic events may be directly explained by the blasting works’ inducing effects. However, at present time, these operations are conducted intuitively based on previous experience as well as on a trial-and-error approach rather than upon an intentional and scientifically justified approach. In result, they are not as effective as they could be.

In the increasingly difficult geological and mining conditions in which extraction of copper ore deposits in Poland is conducted, ensuring an effective and safe mining becomes a key task and a significant challenge for mine operators. Recently, new geomechanical hazards have been identified under these conditions, particularly those related to induced seismicity, which are mostly attributable to high stress values, lower deformability, and higher strength of rocks surrounding a deep deposit. Therefore, more attention should be paid to the dual role played by multiple-faced blasting operations since it is a technology which provides sustainable high production and it is applied for rock burst hazards prevention.

The overall goal of the research is to develop and implement firing patterns for the multifaced blasting production process which may generate the effect of wave amplification due to the interference of postblasting seismic waves [13]. This should increase the capability of inducing stress relief, manifested as a seismic event in a rock mass being mined. As a result, the development of an appropriate group winning blasting procedure may help to develop an effective method which would assist, better than methods used so far, in the stability control in underground workings as well as mitigate risks associated with the dynamic effects of the rock mass pressure. It is possible to achieve a precise delay time between the detonation of explosives in the individual mining faces with electronic detonators. It may be expected that the blasting method described above will provide an amplification of induced blasting wave in a specific location in the main roof strata in which the high likelihood of instability occurrence had been identified. Such amplification is assumed to be able to cause a local instability which in turn may allow for stress relief in the abovementioned high-stress concentration areas. To estimate blast wave propagation in the rock mass, a numerical model of the copper mine was developed and detonation in the mining face is simulated.

Numerical modelling of rock masses subjected to blast loading has long history. Back in 1986, Taylor and his coworkers [3] developed a continuum damage model Taylor–Chen–Kuszmaul (TCK) to characterize the dynamic fracture of rock in tension. In 1998, Maxwell and Young [4] modelled the damage zone using the continuum method. Fairlie [5] described a numerical simulation of street-channelled blast using AUTODYN-2D & 3D explicit codes designed for nonlinear dynamic problems, often referred as hydrocodes. He pointed out that ALE numerical processors, ranging from Lagrangian (grid moves with material) to Eulerian (grid fixed in space), are available for practical applications. AUTODYN-2D has been also utilized by Ma et al. [6] in modelling of wave propagation induced by underground explosion in continuous rock mass. The authors of references [7, 8] presented excellent reviews on the computer simulation of rock blasting conducted till 1990s. Articles [9, 10] and especially review [11] list the number of researches published from the 1990s until today. Review of the literature reveals that most of the research focuses on constitutive modelling or investigation of fracture patterns using simple specimens, or models very often simplified to the 2D domain. Global modelling of an entire bench or mining face is rare. References [12, 13] are examples of research dealing with large-scale models. In both cases, the authors present new approaches to conduct numerical simulations of rock blast. In the first case, the author proposes a connection of damage in rock with ppv. In the second case, the authors propose a new numerical method called the hybrid stress blasting model (HSBM). Reference [11] presents the hybrid discrete-continuous method (finite element method-discrete element method) (FEM-DEM) allowing researchers to model the effect of borehole blasts in a mining face. In [14], the authors simulate fragmentation of the large bench using mixed Lagrange–Euler approach. Due to inherent multiscale nature of fragmentation process, modelling of this phenomenon is extremely difficult. The same observation applies to the presented study, where the shock-induced waves within a single mining face are modelled with taking into account their influence on large rock mass behaviour.

Within the area of numerical methods of mechanics, problems that arise from multiscale character of investigated phenomena are usually addressed by the so-called “global-local” approach. Traditionally, it is done by local refinement of the model in critical areas [15], either manually or automatically (coupled approach) or by creating separate models for global and local analyses (submodelling) [16]. Regardless of adopted methodology, crucial for successful simulation is proper modelling of initial as well as boundary conditions allowing to exchange information between the global and local domain.

To override multiscale-induced problems, the authors present a novel numerical submodelling approach which shall provide a better understanding of the blasting process for practitioners. This approach includes modelling from the face firing to the induced seismic wave propagation in remote areas of surrounding rock mass. While the overall idea of a hybrid, explicit-explicit or explicit-implicit, approach is not new [17, 18], implementation of the submodelling technique and not the coupled approach, together with the scale of the global model, makes the presented research unique. As such, the multiscale technique was never used to model effects of blast wave propagation in the rock mass. This novel technique combines two separate 3D solutions. The small-scale model of the mining face representing the stress-relieving hole (SRH) is used to determine the history of blast pressure, while the global model is used to investigate wave propagation. The global model had to be large enough to allow for investigation of wave influence on the rock burst hazards. Thus, the global model was developed in geological scale, where mining faces are represented by a few finite elements.

The pressure distribution over the SRH, obtained from the small-scale model simulation, was used as an initial condition for the global 3D finite element (FE) model. Subsequently, the peak particle velocity distribution obtained from finite element analysis (FEA) using the global model was compared with experimental outcomes. A reasonable agreement between the results was observed; therefore, the adopted multiscale modelling method confirmed its effectiveness and that it can be successfully implemented in similar problems.

2. Description of Rock Mass Blasting Modelling and Simulation

Two different numerical models were used to conduct FEA of rock mass located around the given room-and-pillar panel in one of the KGHM mines:Case no. 1: small-scale (local) problem—time-dependent blasting process developing in the SRH fired at the chosen mining face of rock wall.Case no. 2: global scale problem—seismic wave propagation within very large volume of surrounding rock mass, following the SRH blasting.

Generally, the proposed coupling procedure is based on the following consecutive phases:(1)The history of the detonation pressure within each side of the SRH was determined. Since the utilizing time step size should cover an extremely short period of time, the model cannot be too large and therefore also cannot represent the rock mass located in more remote locations. However, after detonation, the ground velocity can rapidly drop down and it may stabilize on the level of 5000.0–6000.0 m/s at a distance of few meters from the fired SRH. Thus, one may use the results of the small-scale problem solutions as the initial conditions for the global-scale problem.(2)Taking into account that the main roof strata in Polish copper mines are mostly composed of competent sediment rocks (dolomite/anhydrite/sandstone), it is justified to provide them with equations of state (EOS) adequate for the elastic continuous body. Afterwards, using any standard 3D FEM dynamic code, one may model very large volumes of rock mass surrounding the considered blasted face using small time step size (e.g., Δt ≈ 10−5 s) and the appropriate time-dependent functions of the initial load acting as pressure applied to the correspondent elements (segments) representing sides of the fired mining face.

2.1. Simulation Procedure

Numerical studies of local and global blasting process were carried out using the LS-Dyna commercial code [19]. For the analysis, the authors adopted the FEM explicit central difference time integration, in which the velocity and displacement are calculated as follows:where ; , , and are nodal acceleration, velocity, and displacement vectors, respectively.

2.2. Local Submodelling Definition
2.2.1. Blast Simulation Modelling

The Arbitrary Lagrangian–Eulerian formulation (ALE) algorithm is adopted based on the previous paper [20] and the authors’ experience in modelling and simulation of blast phenomena. Moreover, the ALE method was confirmed to be an effective and reliable method for investigating blast pressure wave propagation and its interaction with various structures [10, 21, 22].

The coupling between fluid-like and solid domains was one of the most challenging tasks during the modelling process. Unfortunately, parameters responsible for the coupling are not universal and have to be estimated for a specific problem [23, 24]. In the presented study, the coupling was simulated using a penalty method, which can be characterized as a fictional spring connecting nodes of two interacting domains (Figure 1). Stiffness of the spring depends on penetration distance between the Lagrangian and Eulerian domains.

Generation of blast wave was modelled using the JWL equation of state. The equation describes behaviour of the detonation products has the following form [19]:where  = ρ0, ρ0 is the initial density of the high explosive; ρ is the actual density of the high explosive; E0 is the detonation energy per unit volume and initial value for E of the high explosive; and A, B, R1, R2, and ω are the empirical constants determined for a specific type of an explosive material based on the experiments [25] using Gurney energy, detonation pressure, and explosion heat.

All required parameters of an emulsion high explosive (HE) used in mining industry were defined based on the results of the so-called cylindrical test (Table 1). Detailed description of the test procedure can be found in [25].

2.2.2. Rock Material Properties

In the local approach, the rock mass was considered as isotropic, with the properties of one type of rock, e.g., dolomite. From the available material models in LS-Dyna, the Johnson-Holmquist II (JH-2) was used for predicting behaviour of the dolomite [26, 27]. The JH-2 model assumes that the strength of material, both intact and fractured, is dependent on pressure, strain rate, and damage.

The normalized intact strength of the JH-2 is given bywhereas the normalized fracture strength is described using the following formula:where A, B, C, M, and N are the material constants,  = P/PHEL is the normalized pressure (P: actual pressure, PHEL: pressure at HEL),  = T/PHEL is the normalized maximum tensile hydrostatic pressure (T–maximum tensile hydrostatic pressure), and is the dimensionless strain rate (: actual equivalent strain rate, : reference equivalent strain rate).

The accumulated damage resulted from the fracture of material is given bywhere is the plastic strain during a cycle of integration, is the plastic strain to fracture under constant pressure (D1, D2: damage constants).

In Table 2, the parameters for the JH-2 for the dolomite are presented. The process of their determination was not a straightforward task, and it required an iterative approach with the coupling of numerical simulations and experimental tests. However, in the presented study, the procedure of their determination is not described.

2.2.3. Model Definition

For the numerical model development in the local scale, a selected mining face of rock mass was modelled in detail. The cubical-like geometry of the rock mass had the following dimensions: width W = 6.0 m, length Lrock = 4.0 m, and height H = 3.5 m. The SRH consisted of four typical production blast holes used for rock mass destressing with the diameter D = 63.5 mm and length LHE = 3.0 m. In each hole, the HE with the mass of 3.5 kg was used. In the FEA, the quarter of the model was considered which drastically reduced computation time. The model consisted of 2,671,306 Lagrangian and Eulerian elements which represented a few very large elements in the structural model (Figure 2) For the coupling between the Eulerian and Lagrangian domains, the second-order advection method was used [19, 23]. A finer mesh was adopted within the area of direct HE-rock interaction for preventing undesirable leakage effect. On the outer walls of the air domain, the pressure flow was allowed by using nonreflecting boundary condition.

Based on the developed model, pressure histories were measured at the mining face sides, representing the finite elements (segments) of the mining face in the global model. Therefore, the pressure histories from the local approach were the initial loading condition in the global modelling.

2.3. Global Modelling Definition
2.3.1. Rock Material Properties

The proposed approach has been studied based on one of the sections in the Lubin mine where copper ore is exploited using the room-and-pillar technology. The copper ore deposit in the analysed area is located at the depth of 742.0 m below the ground and is horizontally flat. It is covered immediately by very thick and stiff main roof strata consisting of 20.0 m layer of dolomite overlaid by 137.0 m thick strong anhydrite plate and 233.0 m thick competent sandstone stratum (Figure 3). The averaged geological data over the considered area and the estimated rock mass parameters are given in Table 3, where σcm, σtm, and σsm are the compression, tensile, and shear strengths, respectively, in rock mass in which values were scaled down according to Hoek’s approach [28] (GSI: Geological Strength Index). It was assured that the overburden rock plates reflect the real lithology in the area and the technological and remnant pillars work effectively within a postcritical phase. The pillars materials reveal a linear-elastic characteristic, whereas copper ore rock mass behaviour is represented by elastic-plastic with a strain softening kind of mechanical model. The entire numerical model’s general boundary conditions were described by displacement-based relationships. More information on the applied solution method can be found in the paper [29, 30].

Except strength and deformation data, the rock mass must be characterized by dynamic parameters between which the structural damping coefficient and maximum frequency value within the ground response range are the most important.

As Preece [31] indicates, the damping coefficient Cd value depends on the type of rock and the degree of its structural fracturing. According to Preece, the minimum damping factor for intact rock which may be employed in blasting calculation is 0.17 while the maximum value Cd should be less than 0.5. For more jointed or even loose rocks, the damping factor may reach a value of 0.7 ÷ 1.0. On the other hand, another numerical modelling [32] revealed that the value of the overall structural damping coefficient assumed as Cd = 0.05 was too low to suppress the rock vibrations of the higher frequency (50 Hz) at the distance of 40–70 m. Based on these estimates and taking into account that roof strata in Polish copper mines are generally composed of competent measure rocks, the further global size numerical simulation utilizes a moderate value of Cd = 0.50.

2.3.2. Model Definition

As a global model for the problem, the multiplate overburden model was accepted with the following simplifying assumptions:(i)Overburden strata consist of several homogeneous rock plates reflecting the real lithology in the area(ii)Technological and remnant pillars work effectively within postcritical phase (elastic-plastic with strain softening behaviour)(iii)The value of carried loads depends on pillar size and actual extraction ratio

The FE model was developed based on the geometry of the selected mining section in the Lubin mine (Figure 4). The global modelling was compared with the actual measurements from the seismic station located at the distances of 380.0 m and 590.0 m in a straight line from the left-side faces and right-side faces, respectively. The global FE model covered the rock mass volume of 2000.0 × 2000.0 × 1725.0 m and was discretized using 13,150,548 brick elements. The main aim of the global-approach modelling was to verify the paraseismic wave propagation within the rock mass due to group production faces’ blasting. During the actual tests, the blasting was performed simultaneously on the 14 left-side faces (see Figure 4) and the paraseismic wave history was measured. This test was also repeated for the 11 right-side faces, with the identical measurements. The discussed tests were reproduced in the global FEA with the applied pressure time histories obtained from the submodelling presented earlier. A detailed description of the local-global coupling is included in the Results and Discussion section. The whole model was fixed at the bottom, gravity load was applied, and the residual stress state was considered as the initial displacement tensor. For the outer faces of the model, symmetry conditions were implemented (displacements in axis perpendicular to the faces were constrained). The discussed global FE model is presented in Figure 5.

3. Results and Discussion

3.1. Pressure Loading Characterization

From the local modelling simulation, the damage of the rock in the mining face was obtained (Figure 6). The red colour represents a fully fractured material: the largest damage can be observed in the area of the four blast holes. Moreover, tensile cracks can be also noticed. In Figure 7, pressure distribution for the selected time moments is presented. It should be pointed out that the range of their values was limited to 100.0 MPa for a better presentation of the results. To transfer from the local to the global domain, several locations corresponding to the discretization of the global model were chosen in the local domain (see Figure 8). Pressure time histories acquired in these locations (see Figure 9) were used as loads and were applied on the faces of the elements representing mining blast faces in the global model.

3.2. Results of the Global Simulation of Rock Mass Destress Blasting

The main objective of the paper was to compare the results obtained from numerical simulations of rock mass dynamic response of group-faces blasting with the actual underground measurements from the seismic station. Example of obtained results in the form of excited vertical velocities distribution is presented in Figure 10.

It is worth noticing that the precision of the nonelectrical detonators used is not high enough since the monitored multiface blasting effects are visibly divided into five isolated period high rock vibrations, which lasted for more than four seconds in total. This kind of blasting cannot be considered as an instant and simultaneous process as one theoretically could assume. Taking this into account, instead of similarity in the velocity distribution, the authors have focused on similarity among ppv. In Figure 11, vertical velocities vs. time characteristics obtained from numerical simulations are presented and compared with the experimental outcomes. The characteristic obtained from blasting faces is reproduced quite well: the oscillations begin nearly at the same time as in the experimental test, and their duration is also similar. It should be mentioned that the global model uses several simplifications, including elastic properties adopted almost in all layers of the modelled section of the Lubin mine.

4. Conclusion and Future Steps

The paper confirms the general ability of computer modelling for rock motion assessment due to dynamic phenomena such as multiple blasting of production faces. As shown above, a reasonable coupling of two different solutions based on a small and global interpretation of the multiface production blasting may produce very encouraging results which may be instantly applied in the prevention of the occurrence of seismic events in underground mines. However, the numerical model parameters (degree of damping, dominant frequency of the motion process, etc.) must be chosen wisely since some incorrect values may give unreasonable output data. The nonelectrical detonators, which somehow generate unpredictable delay times, are also very hard to introduce in the model. Due to the above shortcomings, in the authors’ opinion, a quantitative analysis of the results would be unjustified. On the other hand, the FEA results match well with the measured field data as far as the ppv is considered (Table 4).

Since ppv is widely used in the mining industry as one of the parameters to predict damage caused by blasting [3335], an excellent prediction of this parameter shows that the presented approach can be very useful for blast engineers. Results obtained above may suggest that similarity between computationally obtained and measured data in the time-distribution domain will always be poor unless higher precision of detonators is provided. This may suggest the need for a wider introduction of electronic detonators which can provide a more precise detonation regime for group face firing, which in turn is necessary to make it more effective as a seismic-events inducement approach. Outcomes of the FEA allowed to look into the phenomenon in a broader scope and indicated directions for further field trials where precisely described and selected time sequence could be implemented and checked on-site.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This paper has been prepared through the Horizon 2020 EU funded project on “Sustainable Intelligent Mining Systems (SIMS)” (Grant Agreement No. 730302). The research was also supported by the Interdisciplinary Centre for Mathematical and Computational Modeling (ICM), University of Warsaw (Grant Agreement No. GB73-19).