Numerical modelling of destress blasting – A state-of-the-art Numerical modelling of destress blasting – A state-of-the-art

Abstract As a proactive mine safety measure against the occurrence of rockburst, destress blasting has been applied to numerous mining conditions to precondition highly stressed rock mass to mitigate the risk of rockburst occurrence in deep mines as well as in deep underground constructions. However, the application of destress blasting mostly depends on engineering experience, while its mechanism and efficiency have not been well understood. Rapid advances in computer technology have made numerical simulation an economical and effective method to study the rock blasting effect. Enormous research efforts have been made to numerically investigate the blasting fracture mechanism, optimize blasting design, and assess the ef ﬁ ciency of destress blasting. This review focuses on the state-of-the-art progress in numerical modelling associated with destress blasting over the last two decades. Some commonly used modelling approaches for destressing blasting are compared and reviewed. Currently, two different ways of modelling based on static and dynamic modes are typically used to study the effect of blasting. In the static method, destress blasting is simulated by modifying the rock mass’s stiffness and strength properties to obtain the post-blast stress state in the destressed zone. The dynamic modelling technique focuses on the dynamic fracture process of coals and rock masses, during which the predetermination of the damage induced by blasting is not necessary. Moreover, the extent of damage zones around the blast hole can be precisely estimated in the dynamic modelling method by considering time-varying blast pressure and strain rate dependency on the strength of rock mass but at the cost of increased computation and complexity. Besides, different destress blasting modelling methods, generally classified into continuum-based, discrete-based, and coupled methods, are compared and reviewed. The fracture mechanism of blasting in the rock mass is revealed, and the destressing efficiency of the existing destress blasting design is assessed and compared with classical results. The factors that may affect the efficiency of destress blasting are summarized. Finally, the difficulties and challenges associated with the numerical modelling of destress blasting are highlighted briefly.


M ining in deep underground mines continues
to face the challenges of rockburst caused by high mining-induced stresses after excavation. Rockburst is a common disaster in deep underground mines, causing injury to mine operators and damage to underground work. Thus, protective measures are considered an effective way to deal with the adverse effects of excessive stress. Different attempts have been made to address the excessive stresses, and destress blasting is considered the proactive measure during underground mining [1e4].
The concept of destress blasting was reported in the early 1920s in the Nova Scotia province of Canada after a devastating disaster in Springhill Colliery. It has been viewed as an effective and important technique for managing strainburst hazards [5]. Destress blasting aims to migrate the stress concentration zone interior to the active mining face to reduce the potential of rockburst by generating blast-induced fracture networks (Fig. 1). The fracture generation in rock mass requires explosive detonation in a sequence of long holes ahead of the face or stressed pillars. Once the explosives are detonated in these boreholes carefully, damage and fractures would occur in the rock mass surrounding the blastholes. Thus, the highly fractured zone ahead of the working faces cannot carry the high magnitude of stress level anymore due to its low strength and stiffness, leaving a protective barrier between the working face and the stress concentration zone.
Up to now, destress blasting has been applied to numerous mining conditions to precondition highly stressed rock mass in deep mines. However, the design and implementation of destress blasting schemes extremely rely on engineering experience. Moreover, the field application of destress blasting in hard rock mines encountered many difficulties and challenges. The factors affecting destress blasting efficiency have not been grasped thoroughly. Suitable tools or methods are imperative for the optimization of blasting design in engineering site [6e8]. Blasting in the rock mass is a very fast and complex process, lasting for a few milliseconds, which involves a complex interaction between high-energy denotation products and rock mass. The dynamic response and fracture mechanism of rocks are difficult to fully understood based solely on theoretical and experimental studies. With the rapid development of computer technology and advanced numerical techniques, more detailed and reliable predictions of rock mass damage and fracturing under blast loading become available [9e11]. Numerical modelling provides a convenient, economical, and relatively accurate approach to studying destress blasting in the underground rock mass, especially for complicated cases where experiments and theoretical solutions are difficult to conduct [12e14].
Numerical simulation of destress blasting helps understand the rock dynamic response and destress mechanism under blasting loads, which are conducive to the field application of destress blasting. To accurately predict the blast-induced damage and assess the efficiency of destress blasting, numerous studies have been carried out based on different numerical tools and methods. This paper provides a comprehensive review of previously published literature on blasting modelling, involving the commonly used simulation methods, the characterization and representation of the denotation process, rock dynamical properties, and potential factors heavily affecting destress blasting performance and simulation modelling results. Moreover, based on the literature review, the fracture mechanism of rock mass subjected to blasting is first studied, and the destress efficiency of existing destress blasting schemes is further assessed.

Destress blasting simulation approach and method
In the past decades, the rapid advancement of computation capability has facilitated the use of computer technology to modelling complex blasting in rock masses [15]. Modelling of destress blasting requires considering several aspects to accurately reproduce the blasting process and rock response. Two different approaches, i.e. the static and dynamic approaches, have been used to perform destress blasting modelling and study the blasting effect. In terms of the traditional static method (i.e., equivalent modelling approach), destress blasting is simulated by modifying the mechanical properties and the post-blast stress state of the rock masses in the specified destressed zone. Note that the static modelling approach is not complicated and requires less time and computational capacity compared to the dynamic method. On the other hand, based on the developed dynamic modelling technique, the dynamic fracturing process of rock masses can be presented, in which the predetermination of the damage scope is not necessary. Moreover, the dynamic modelling method can give a more precise extent of fractured zones. These advantages make the dynamic modelling method more accurate than the static modelling approach. The comparison between the two approaches is summarized in Table 1.

Static approach
For the static approach, the blasting process is not taken into account, and the blast-induced damaged zone is assigned artificially when modelling. Therefore, the static approach is an equivalent blasting modelling approach. The numerical methods commonly used in the static approach are the finite element method (FEM) and finite difference method (FDM).
The scope and parameters of the weakened zone need to be obtained from actual blasting measurements or empirical formulas. Two parameters, namely rock fragmentation factor a and stress dissipation factor b (a and b ranging from 0 to 1), were proposed by Tang et al. [8] to weaken the rock mass in the blast-induced fracture zone. Thus, elastic modulus E of rocks is reduced to a E in the destressed zone, and stress is scaled by (1-b). Then, these scaled parameters were applied to simulate the effect of destress blasting. Moreover, suggested that a should be in the range 0.4e0.6, and the values of b should be 0.4 or more [16]. Using proper a and b, the equivalent blasting modelling approach was successfully applied to simulate the destress blasting in the mining panel [13,17]. It is worth noting that this approach assumed a uniformly distributed damage zone extending over the entire drift face, but whether these destress blasting practices can lead to full-face stress relief is questionable.

Dynamic approach
In the dynamic approach, the actual fracturing process of rock mass under blast loads can be obtained, and the assignment of the damage scope is not necessary. Numerous in-house and commercial software packages were utilized to model the blasting response of rocks. Existing numerical methods involved in dynamic modelling can be classified into continuous, discontinuous, and coupled methods [18].

Continuous method
The continuous method is suitable for problems whose system is a continuum with an infinite degree of freedom, whose behaviour is dominated by the governing differential equation and the continuity conditions at the interfaces between adjacent elements [19]. This type of method is mostly applied for non-fractured or slightly fractured rock masses. Continuum-based method for blasting modelling includes FEM and FDM. Several FEM codes, such as LS-DYNA [20,21], ABAQUS [22,23], AUTODYN [24,25], and DFPA [26,27], have been used to A fluid-solid coupling algorithm was adopted to analyseof the explosive detonation [25]. For example, Zhu et al. [28] used the AUTODYN code to study the fracture mechanism of rocks exposed to blasting and optimize blasting design. Ma and An [20] employed the LS-DYNA code to study the effect of loading rate, distance from a free face, in situ stress, and joint planes on the rock fracture patterns. With a combination of FEM and cellular automaton technique, a cellular automata software for engineering the rock mass fracturing process (CASRock) was developed to simulate the dynamic response of rock mass induced by blasting [29,30]. Fig. 2 shows the rock fracturing degree (RFD) under blasting loads modelled by CASRock, where RFD ! 0.8 means that the rock mass begins to damage, and RFD ¼ 2.0 means that the rock mass is completely unstable. The results showed that CASRock can well simulate the damage and failure process of the rock mass around the blasthole after blasting. Based on the FDM, the FLAC and FLAC3D packages developed by ITASCA consulting are also available for simulating blast-induced fracture zones by inputting various types of pressure functions [11,31]. For example, FLAC3D was adopted by Yilmaz and Unlu [31] to study the dimensions of the damaged zone in rock mass exposed to dynamic stress waves. Wang et al. [32] developed a finite difference code to explore the wave propagation in rock plates subjected to blast loads.

Discontinuous method
The discontinuous method includes the distinct element method (DEM), such as particle flow code (PFC), universal distinct element code (UDEC), discontinuous deformation analysis (DDA), and other discontinuous methods e.g. lattice model. The discrete element method (DEM) involves the definition of an assembly of indivisible elements for medium and proper interactions between elements for rock responses [33,34]. Finite displacements and rotations of discrete bodies, including complete detachment and automatically recognizing new contacts, are allowed in DEM [33]. Therefore, the properties of rock mass, such as joints and porosity, and dynamic responses characterization during blasting, such as the large deformation, fracturing, motion, and rotating or dispersing of elements, can be well modelled. For example, the particle flow code PFC3D was employed by Sarracino [35] and Potyondy et al. [36] to study the combined effect of stress waves and high-pressure gas flow. UDEC was adopted to simulate the dynamical response of rock masses subjected to stress waves [9,37]. DDA is an implicit discrete element method developed for modelling the static and dynamic behaviour of discrete block systems [38], which can solve the dynamic problem of fractured rock mass [39]. The lattice model is another type of discrete model. In the framework of the lattice model, point masses or particles interconnected through massless springs are adopted for rock materials subjected to dynamic loads [6,40].

Coupled method
Numerical modelling of the dynamic response of rock mass exposed to blasting loads involves several complicated physical processes, including the charge explosion, the crushed zone formation, the propagation of stress waves, and rock fracturing. Thus, it is quite difficult to include all using one method. Thus, the coupled method is a good choice to take advantage of two or more numerical methods. Wang and Konietzky [37], Chen and Zhao [41], and Deng et al. [42] simulated the charge explosive in the LS-DYNA or AUTODYN code to obtain the pressure-time history on the borehole wall. Then, the explosion history was imported to the UDEC to simulate the dynamic response and fracture patterns in rock masses. The coupling between the AUTODYN and PFC was also adopted by Ledoux [43] to model the fracture patterns of grout cylinders induced by the blast loads.
The hybrid method is another kind of coupled method by directly coupling the mesh of FEM or FDM with that of DEM, which becomes more and more popular in blasting modelling of the rock mass. Minchinton and Lynch [44] used the combined finite element-discrete element program MBM2D to model the rock fragmentation and heave during the blasting process. Mitelman and Elmo [45] employed a hybrid FEM-DEM code ELFEN to study blastinduced damage in a circular tunnel. The ELFEN [46], commercially developed and marketed by Rockfield Software Ltd, allowed realistic modelling of the transition of the rock mass from a continuum to discontinuum, with the explicit generation of stressinduced cracks. Fakhimi and Lanari [47] recommended a hybrid discrete element-smooth particle model (DEM-SPH) to simulate rock blasting. In their study, the rocks and the explosive gas were modelled by discrete circular disks and smooth particles, respectively, and their interaction was assumed to follow a perfectly plastic collision. The Hybrid Stress Blasting Model (HSBM) is a continuum and discontinuum hybrid approach, which combines the explicit finite difference code FLAC and a newly developed lattice scheme [48,49]. The code has been developed for several years through an international collaborative research project and seems to be promising in capturing many features of dynamic blasting, such as detonation, wave propagation, rock fragmentation, and muck pile formation [6,10].

Characterization of stress wave and explosive gas pressure
In the dynamic modelling approach, the characterization of blast loads is an important topic. The explosive blasting can be characterized by two distinct stress mechanisms:the shock wave pulse and the explosion gas pressure. The detonation of the explosive initiates a rapid chemical reaction and converts the solid mass of the explosive into gaseous products [50]. The denotation products impact the borehole wall and generate an intense, high amplitude compressive pulse. The dynamic stress field is generated when the shock (stress) wave propagates through the rock mass, which occurs in milliseconds. Then, the explosive gas would act on the fractured zone in rock mass induced by stress waves. The gas penetration is generally considered quasi-static, which occurs on the order of second [51]. The characterization of the pressure pulse and the explosive gas pressure is extremely essential in rock blasting modelling.

Pressure pulse on the borehole wall
Borehole wall pressure describes the expansion work of the explosive during explosive detonation, which is one of the most important indexes evaluating explosive performance. However, measuring the pressure applied to the blasthole wall is difficult due to the lack of feasible methods and tools. Instead, some empirical formulas or detonation theories are commonly utilized to estimate the pressure pulse.
The explosive detonation in rocks can be classified as ideal-and non-ideal detonation. For the ideal detonation, the explosive shock front propagates through the material at supersonic speed. If the detonation is non-ideal, the reaction zone is relatively long compared to that from the ideal detonation, where the reaction zone is assumed to be zero thickness [4,50]. Four commonly used approaches can be used to approximate the blasthole wall pressure time-history: (1) equation of state; (2) pressure decay function; (3) predefined wave function; (4) other methods.

Equation of state
Thermal-chemo-physical transformation occurs during the explosive denotation process [15]. Thus, the equation of state (EOS) has been a common approach for estimating the borehole wall pressure. The Jones-Wilkens-Lee (JWL) equation is the most commonly used due to its simple forms, experimental basis, and ease of calculation [7]. The JWL EOS is written as [52,53]: where P is the blasthole wall pressure, and A, B, R 1 , R 2 , and u are explosive constants. E 0 is the internal energy per initial volume (Pa), and V is the relative volume (m 3 ). The JWL EOS is a high-energy combustion model and is usually used for high explosives such as TNT (2,4,6-trinitrotoluene), RDX (hexahydro-1,3,5-trinitro-1,3,5-triazine), and HMX (octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine) [54], which has been widely applied to determine the pressure pulse [55,56]. Even though its a good performance in simulating the ideal detonation, the JWL is considered to be inaccurate for describing nonideal detonation of explosives, such as ANFO (full name: ammonium-nitrate-fuel-oil) [7,57]. Because the reaction zone at the detonation front of the non-ideal detonation is thick compared to that from the ideal denotation, Williamsburg equation of state is recommended for modelling non-ideal detonation [57,58]. However, few studies have used the equation of state to describe non-ideal detonation due to the complexity of non-ideal denotation.

Pressure decay function
The pressure decay function is another approach to generate the time-decaying pressure pulse.

Reference
Pressure decay function Remark Sharpe [59] P ¼ P P e Àat P P is the peak pressure applied on the borehole wall, a is constant, and t is the time Duvall [60] P ¼ P P ðe Àat À e Àbt Þ a and b are constants Jiang et al. [61] P ¼ P P t n e Àat n and a are constants, and the peak pressure occurs at the time t r ¼ n=a Cho and Kaneko [26] P ¼ P P xðe Àat À e Àbt Þ x ¼ 1=ðe Àatr À e Àbtr Þ, and the rise time to peak pressure t r ¼ ð1=ðb À aÞ)log(b=a) Trivino and Mohanty [62] PðtÞ  Table 2.

Predefined pressure function
Some specific functions with simple mathematical shapes were employed to approximate the pressure pulse from the explosive detonation. Donze et al. [34] employed a Gaussian time function to generate the pressure pulse. Zhu et al. [63] input triangle wave, sine wave, and Gaussian wave into the elastodynamic finite element model for back-analysis, and their results indicated that the Gaussian wave was more suitable for high-pressure air blasting modelling. Mitelman and Elmo [45] adopted a triangular-shaped pulse for ELFEN modelling. Based on a self-developed cellular automata software for engineering rockmass fracturing process (CASRock), Pan and Mei [30] used the triangular load function in their dynamic version to simulate the blast loads acting on borehole walls. Likewise, the triangular load function was also applied to simplify the blasting loads by Simons [64] and Yan et al. [65]. Resende et al. [66] employed a Ricker wavelet, a negative exponential function with no high-frequency corners, to carry out rock blasting modelling. However, these simple functions cannot provide a realistic representation of the pressure time-history from the detonation process [41]. Moreover, these simple functions carry no physical meaning. Even that, it provides an approach to simulate the blasting source, and their simulation results can qualitatively reveal the dynamic response of rock mass exposed to blast-induced stress wave.

Other approaches
Besides the above mentioned, there are some other ways to generate a pressure-time history. Xia et al. [67] used the typical velocity history recorded by the vibration sensors installed in the underground tunnel as input for blasting modelling. Saharan and Mitri [7] proposed the optimized pressure profiles to describe the borehole pressure. In their study, the peak borehole pressure and the pressure at different time scales were required to determine the optimized pressure-time profile. The peak borehole pressure (Pa) was estimated as follows [68], where r is the explosive density (kg m À3 ), VOD is the velocity of detonation m s À1 Â r c is the coupling ratio, i.e. the explosive diameter and borehole diameter ratio, and g is a constant, called the adiabatic exponent. In their study, the peak pressure of ideal and non-ideal detonation were calculated by properly importing the value of the VOD and r of emulsion explosive and ANFO explosive, respectively. Then, the rise time was also given for ideal and non-ideal detonation, respectively. Finally, the peak pressure and the rise time were applied to the proposed optimized pressure-time profiles. The authors stated that the optimized pressure-time profile gave a better approximation of the real pressure pulse than the Gaussian and triangular shape functions' profiles.

Explosive gas pressure
The explosive gas pressure plays an important role in rock blasting, but it is considered to be effective only when the rock mass has been preconditioned by the stress wave [51,69]. The explosive gas would squeeze into the blast-induced cracks or pre-existing fractures and expand its own space. The gas penetration process is coupled with the deformation and propagation of existing fractures or rock joints, resulting in more expansion and penetration of explosive gas into crack openings [70].
The gas flow in the blast-induced cracks and preexisting fractures is generally described with the gas flow model. Some studies utilized simplified gas models to calculate gas pressure. For example, Munjiza [71] assumed constant pressure distribution near the blast hole, neglecting the pressure change with the variable opening and length of different cracks. Thus, the gas pressure is only a function of the volume of the deformed hole and fractures [72]. considered explosion gas penetration during rock blasting by using a discontinuous deformation analysis method. The total volume for gas flow was tracked, and instant explosion gas pressure was derived using a simple polytropic gas pressure equation of state. Similarly, Yuan et al. [73] also utilized a simple polytropic equation of state to calculate the uniformly distributed gas pressure for the entire crack network. Sim et al. [74] considered two types of pressure distribution inside cracks, i.e. uniform and linear pressure distributions, to study the effect of gas pressure on crack propagation.
Several gas flow equations governed by conservation of mass and momentum have been widely employed to calculate explosive gas pressure. Cho et al. [27] and Goodarzi et al. [75] used the onedimension transient flow proposed by Nilson et al. [76] to describe the gas flow through fractures. Munjiza et al. [77] proposed a new approach to simulate the gas flow through a cracked medium, in which the gas flow through the cracks and fractures around the blasthole was modelled as an equivalent one-dimensional flow through several constant area ducts. From a different view, gas models based on porous media flow were also proposed. Preece and Taylor [78] employed an equivalent porous medium to simulate the blast-induced rock motion process. Mohammadi and Pooladi [79] proposed a new approach extending the original idea of the uniform gas flow in Munjiza et al. [77] to non-uniform isentropic gas flow. Later they also extended this idea further [70].

Overall consideration of stress wave and explosive gas pressure
Even though the respective contribution of denotative stress and explosive gas have been separately studied, numerical codes that can handle both the shock wave-induced dynamic rock fracture and the gas penetration are scarce [80]. Due to the different mechanisms and time scales involved, it is difficult to consider both the stress wave and explosive gas pressure during blasting modelling [70,81].
In recent years, some innovative studies have been performed to solve the coupled effect of stress wave and gas pressure in blasting modelling. The smooth particle hydrodynamic (SPH) method, pioneered by Russian scientists in the late 1950s, can handle large deformation as no real mesh exists in the SPH method [82]; thus, it has been used to perform blasting modelling over the years. Pramanik and Deb [83] studied the dynamic responses of rocks subjected to stress waves and explosive gases in the framework of SPH. In addition, SPH can also be combined with other methods to describe the blasting response of rock mass. For example, Fakhimi and Lanari [47] utilized a bonded particle system to mimic the dynamical reaction between the rocks and the denotation products. By modelling the borehole with and without a thin copper lining, Lanari and Fakhimi [84] used an improved hybrid SPH-BPM method to clarify the respective contribution of the shock wave and denotation gas on rock fracturing. Moreover, some other attempts have also been presented to simulate the whole blasting process, considering both the stress wave propagation and explosive gas penetration. An immersed-body method and a cohesive zone fracture model were also adopted by Yang et al. [81] to model the stress wave propagation, fracture extension, gasesolid interaction, and flying fragments in a practical complicated blasting engineering.

Direct method
In the direct method, rocks are idealized as a collection of structural units, such as springs, beams, etc., which are bonded together at their contact points. The breakage of individual structural units or bonds represents the formation of microcracks. The numerical codes for indirect approaches are various kinds of discrete element method (DEM) code.

Homogeneous and heterogeneous properties
Most numerical modelling studies generally treat rock masses as continuous and homogeneous mediums [24,25]. However, the rock mass contains many natural weaknesses, such as pores, grain boundaries, and pre-existing cracks, and inhomogeneity plays a significant role in the progressive failure process of rocks. Two approaches are commonly used to consider the inhomogeneity of rock materials in numerical modelling. One is to examine the comprehensive effects of discontinuities by using equivalent material properties of the rock mass, and the other is to examine the effects of a few discontinuities in the rock mass. Pisarenko and Lebedev [85] presumed that many fissures and defects exist in the rock mass, and their distribution is random. For example, the Weibull distribution has been widely used in numerical modelling to computationally produce heterogeneous rock mass [86e90]. When few joints with definite length and spacing exist in rock masses, proper modelling of those joints is critical [20,91].

Dynamic constitutive model
Rocks exposed to dynamic loads exhibit large strains, high strain rates, and high-pressure characteristics. A proper constitutive model should be used for modelling the rock fracturing under blasting loads. Currently, numerous dynamic constitutive models have been developed for dynamic responses of several materials [92e94]. In those models, commonly used constitutive models for rock dynamics modelling include the Johnson-Holmquist model series [95,96], the Holmquist-Johnson-Cook (HJC) model [97], and the Riedel-Hiermaier-Thoma (RHT) model [98].
The Johnson-Holmquist1(JH-1) was first proposed to account for large deformations [95]. Based on the JH-1 model, the Johnson-Holmquist 2 (JH-2) model was improved, and progressive damage with increasing deformation was considered [96]. In the JH-2 model, the dynamic behaviour of brittle materials such as rock and ceramics under the blast and impact loads are described through the strength model, polynomial EOS, and damage model. The natural logarithmic function of the equivalent strain rate is used to describe the strain rate strengthening effect. The JH-2 model has been widely adopted in rock blasting modelling [91,99,100].
The HJC model was originally proposed by Holmquist et al. [97] for impact computations and had relatively few parameters need to be determined. Thus, it has been considered a simple and effective model for impact, penetration, and blasting modelling. The HJC model is mainly valid under the compression-dominated stress states while not suitable for reproducing tension-dominated phenomena. Later, modifications of the HJC model were conducted by considering the compressive behaviour [101], strain-rate effects [101,102], and tensile behaviour [103].
RHT, introduced by Riedel et al. [98], is an advanced plasticity model capable of modelling the dynamic behaviour of concrete or rock-like materials when subjected to high strain rate loading conditions, such as impact and blasting [93,104]. The shear strength of the model is described through three limits, i.e. the initial elastic yield surface, the failure surface, and the residual friction surface. Thus, PHT is capable of considering pressure hardening, strain hardening, strain rate hardening, and strain softening.

Mechanism for destress blasting
Since the 1950s, many theories have been proposed to explain the fracture mechanism of rocks exposed to blast loads. It is believed that rock fracturing induced by the explosive explosion is the result of both actions of denotation wave and the explosion gases pressure. As illustrated by [51], the fracturing of rocks around the blastholes experiences three phases. Firstly, the shock wave causes a strong impact on rocks near the borehole and crushes the rocks. Then, stress waves attenuate from shock waves, propagate in the rocks in the radial direction and break the rocks with long radial cracks. Finally, explosive gases penetrate cracks and further extend radial fractures.

Role of the detonation wave
Numerous numerical studies focused on the dynamical response of rock mass subjected to blastinduced detonation wave [20,25,34,93]. For example, Xie et al. [93] adopted LS-DYNA to study the fracture process and mechanism of rocks exposed to blast loads (Fig. 3). Once an explosive charge is initiated, the detonation wave initially acts on the borehole wall working as a shock wave and then propagates in the radial direction of the borehole. The shock wave has a pressure much higher than the dynamic compressive strength of rocks, creating a crushed zone near the borehole (Fig. 3a). The shock wave is attenuated to the stress wave after a considerable amount of energy is consumed in the crushed zone. The stress wave poses the rocks to intense radial compression, which causes many radial cracks and forms a fractured zone outside the crushed zone (Fig. 3b). In the fractured zone, a large number of radial cracks are initiated from the boundary of the crushed zone, but only a few of them extend far due to the stress relaxation induced by the longer ones [105,106]. As a free face exists, the compressive stress wave would reflect from the free face as tensile and shear waves. If the reflected tensile wave is strong enough to exceed the dynamic tensile strength of rocks, the spalling phenomenon may occur back towards the interior of the rock, as shown in Fig. 3c. Some radial cracks may further extend when they interact with the reflected tensile wave (Fig. 3d). The simulated crack patterns match well with the experimental result in Jung et al. [107]. A similar fracture pattern under blast-induced stress wave has also been reproduced using the discrete element method and coupled method [40,73,108].

Role of explosion gases
To explain the contributions of the explosive gas in rock blasting, many researchers have adopted feasible workarounds. One common method in laboratory tests and numerical studies is comparing the blast-induced fracture patterns in a blasthole with and without a thin liner or tube [109,110]. Liners installed in the blasthole prevent gas penetration into the stress wave-induced cracks. Thus, the respective contribution of stress waves and explosive gases can be understood. Series of experiments have revealed that the additional fracturing would be generated in rock mass due to the explosive gas penetration [111e113]. The conclusion has also been validated in various numerical studies. For example, Lanari and Fakhimi [84] studied the respective contribution of stress waves and explosive gases by modelling rock blasting with and without using a thin copper lining. They found that most of the cracks in the rock masses are induced due to the stress wave, and the fracture increasing rate by the stress wave is much higher than that by the explosive gas penetration. Around 25% fewer cracks are observed in the rock when a liner is installed in the borehole. It also supports the earlier assertion that a large proportion of cracks is introduced by the stress wave. Based on the granular discrete element method, Yuan et al. [73] revealed that the explosive gas can enlarge the crushing area and increase the crack aperture. Some other numerical studies also indicated that the explosive gas can significantly extend existing fractures and cause fragmentation within a solid [74,75,79].

Optimization for destressing blasting design
Tang et al. [8] adopted the static approach to compare the destress efficiency of three different blasting patterns. The three blasting patterns, i.e. all-perimeter destressing, full-face-and-corner destressing, and wall-and-corner destressing, are examined and compared in a parametric to situate the full-face pattern commonly practiced in hard rock mines, which is based on a worldwide survey performed by Comeau et al. [114]. In the static approach, a uniformly distributed damage zone is assumed to extend over the entire drift face in the traditional static modelling approach. However, it is questionable that these destress blasting practices can lead to a full-face destress. Sainoki et al. [11] developed an alternative static-dynamic numerical modelling approach and compared the resulting destress blasting efficiency with that calculated by the traditional static approach. The magnitude of blast-induced damage zones from a dynamic model was used in the alternative numerical modelling approach. A typical destress blasting pattern was employed. The destressed zones used in the two approaches were presented. Then, the static methodology proposed by Tang and Mitri [115] was utilized to evaluate the effectiveness of destress blasting in two approaches. The numerical results in terms of major principle stress before and after destressing are presented. Compared to the result from the traditional modelling approach, the alternative modelling approach generates higher postdestressing stress. The traditional modelling approach suggests lower post-destress stresses, while the alternative modelling approach reduces stresses insignificantly. The study performed by Sainoki et al. [11] also indicated that common destress pattern may not lead to effective stress relief, and some improvement in the blasting scheme must be made to create uniform blast energy distribution ahead of the working face.
When carrying out blasting in rock masses, radial fractures extend in all directions in the absence of in situ stresses. However, the high in-situ stresses around the working face often limit blast-induced radial fracturing and confine preexisting natural fractures, which also inhibits the dilation of these fractures from gas penetration. Moreover, anisotropic stress tends to guide cracks to propagate preferentially along the maximum principal stress, which has been proven by several experiments and numerical studies (see section 6.9). Thus, the length and orientation of the blast-induced radial cracks are greatly affected by in situ stresses. For example, Saadatmand Hashemi [116] used LS-DYNA code and performed a numerical model on a destress blasting pattern containing two face holes with a spacing of 2.0 m and four corner holes at 30 to the direction of the face advancement with and without in situ stresses. The effective connection between radial cracks cannot be achieved due to the great distance between the corner-drilled holes and face holes caused, even without the in situ stresses. The condition would get worse when the in situ stresses are applied, and such damage packets might impose excessive stress concentration under high in-situ stresses and increase the probability of a rockburst [6,116]. As mentioned above, conventional destressing practices have limited efficiency in mitigating potential strain bursts around the working face.
Drover et al. [6] proved that the conventional destressing patterns provided limited fracturing and interaction between the damaged zones of these widely spaced blastholes, resulting in insufficient stress releases, especially in the presence of anisotropy stress fields. Saharan and Mitri [5] suggested that the shear failure mechanism is more beneficial for strain energy dissipation within rocks in destress blasting. Based on that, Drover et al. [6] proposed a conceptual destress blasting design that considered the likely length and directionality of radial fracture networks in the high-stress condition. The destress blasting pattern consists of a series of symmetrical rows of destressing charges, which are sub-parallel, yet almost oblique to the major principal stress (Fig. 4). Besides, the HSBM was used to optimize the design parameters, such as borehole diameter, spacing, and explosive type (Fig. 5). Note that the fracture interaction occurred efficiently between individual destressing charges, and series of parallel and continuous linear fracture planes were introduced throughout the face. It can be seen that the destressing concepts proposed by Drover et al. [6] provided a new solution for optimizing the destress blasting design under high-stress conditions.

Factors affecting destressing blasting
The fracturing process of rock material subjected to blast loads is complex, and the outcome of a blast Fig. 4. A novel destress blasting concept [6]. is easily affected by some parameters. The most important factors that may affect blast-induced rock fractures should be evaluated and understood. In destress blasting practice, the radial fracturing generated by a confined destressing charge is closely related to explosive properties, rock properties, natural joins, and in situ stress fields. Moreover, some other techniques initially developed for improving blasting efficiency or special purpose are also discussed here.

Rock materials
Rock-like materials with different mechanical properties and structural characteristics may have different blasting performances, depending on various factors such as physic-mechanical properties of rock, i.e., porosity, bulking, density, elasticity, plasticity, brittleness, rock strength, and other properties. Cho and Kaneko [26] examine the influence of rock inhomogeneity on fracture patterns by varying coefficients of uniformity. It has been suggested that extensions of the radial cracks around the borehole are related to rock inhomogeneity. On the other hand, rock compressive and shear strengths influence the hole expansion and the rock crushing, and the fracture toughness of rock mass determines radial crack extension or propagation [4,84]. Lanari and Fakhimi [84] studied the effect of rock strength on the blast-induced fracture patterns using a hybrid discrete elementsmooth particle model. It has been revealed that the blast-induced fractures in the stronger rock are more localized than those in the weaker rock. Moreover, the extensive damage and the resulting debris in the weaker rock would prevent the explosive gas from penetrating into cracks, which weakens the effect of the explosive gas. Wei et al. [55] numerically studied the effect of different RMR values on fracture patterns in the rock mass, and similar results were obtained. Mitelman and Elmo [45] employed the hybrid FEM-DEM code to simulate the fracturing of different rocks subjected to stress waves. It has been revealed that the larger extent of fracturing occurs in the weaker rock and dissipates more energy. Experimental results also showed that the properties of rock materials affect the transmitted pressure, blast-induced breakage and fracture patterns [110,117]. Moreover, rock mass permeability is considered the main control factor for the penetration of explosive gas [106].
The joint structure has a great effect on the blastinduced fracture extension. Wang and Konietzky [37] performed the blasting modelling in the joint rock mass, and two joint patterns, i.e. randomly polygonal joints and orthogonal joints, were considered. It has been found that the crack evolution exhibits strong anisotropy in the rock mass containing orthogonal joints compared to that containing randomly polygonal joints. The effect of different joint patterns on the damage of the rock mass subject to shock wave was also studied by Deng et al. [9]. Therefore, destressing blasting modelling and design should carefully consider the rock properties and joint conditions.

Charge properties
The mass of the explosive products is directly associated with the charge properties. The geometry of the explosive charge, such as charge shape, charge diameter, and charge length, significantly affect the wave interaction in the borehole. Srirajaraghavaraju [110] found that the transmitted pressure was affected by the charge shape in his experiments and transmitted pressure from a point source (concentrated charge) was significantly greater than that from the line source (detonating cord) as the same charge weight was utilized in two sources. However, the attenuation from the point source was greater due to the spherical propagation of the wavefront of the transmitted pulse compared to the conical spreading from the line source. At a given explosive charge, the increase of the charge length and charge diameter will increase the charge weight, resulting in increasing denotation products and explosive energy. Therefore, the peak pressure on the borehole wall increases, leading to an increased intensity of fracturing [109,110,118]. Both the numerical studies and experiments have confirmed the point [119e121]. The charge type influences explosive energy, explosive energy efficiency, and velocity of denotation (VOD), thus significantly affecting the blasting behaviour [122]. Experiments carried out by Khademian and Bagherpour [123] also indicated that the type of explosive has a great influence on rock fragments and inner microcracks.

Detonation initiation location
The detonation initiation location of the charge column affects rock fracture patterns. As the detonation wave spreads along the charge, the stress wave induced by detonation also propagates to the other side of the hole in the conical wavefront. Thus, the superposition of stress waves would cause a highenergy and high-stress zone in the hole far away from the detonation position [124,125]. Several detonation locations, such as bottom initiation, middle initiation, top initiation, and multi-point initiation, have been considered in blasting operations [4]. Several numerical modelling has been performed to study the effect of detonation initiation location on rock fracturing [124e126]. It has been found that the most serious damage occurs near the top for bottom initiation, while the most serious fractured zone occurs near the bottom for top initiation; the middle initiation takes advantage of the bottom and top initiation; the cumulative effects of stress wave in the top and bottom regions are all relatively good; the most fractured zone appeared at the centre area as the charge is initiated simultaneously at the top and bottom [28]. Numerical modelling in HSBM by Drover et al. [6] indicated limited vibrations at the face and greater rock mass damage occurred as the destressing charge is collar primed.

Decoupling and coupling
A fully coupled explosive is the case where the explosive charge fills the diameter of the blasthole. Conversely, the decoupled explosive means that gaps exist between the blasthole wall and the explosive charge. The coupling ratio is often defined as the charge diameter divided by the borehole diameter. The coupling ratio controls the blastinduced fracture patterns. Compared to the coupling under the same conditions (explosive type, charge weight, rock, etc.), the decoupling involves pressure pulse with lower peak pressure and long duration time. The explosive pressure transmitted to rock masses is reduced due to the air gap between the borehole wall and the explosive [4,119]. Thus, the size of the crushed and fractured zone decreases with the increase of the decoupling ratio [127e129]. The decoupling technique can be used to control blast-induced vibration and reduce the disturbed zone around holes. For example, Zhu et al. [28] adopted three types of coupling materials, i.e. water, sands, and air, to study the coupling medium on blast-induced damage. It has been found that air-coupling caused the least extended fractured zones, while water-coupling generated the greatest fractured zone because the water has the greatest ability to transmit stress waves. Therefore, the decoupling technique has the potential to enhance rock fragmentation. Experimental results also indicated that the fluid-filled borehole gives a better coupling with the borehole wall while preventing crushing of the wall as found in fully loaded holes [109,118,130].

Ideal and non-ideal detonation
Ideal or non-ideal may occur in an explosive detonation. For the ideal detonation, all the energy released within the reaction zone is used to drive the detonation shock front and drive it forward. However, the non-ideal detonation merely uses partial chemical energy to drive this detonation front, and some of the explosive material may react after the sonic plane [131,132]. An ideal detonation is of limited practical use due to the wide application of commercial explosives in mining engineering. Factors such as confining medium, coupling ratio, and explosive charge may affect the degree of non-ideal detonation [133]. Ideal detonation reaches its peak pressure in a short time and decays rapidly, while non-ideal detonation takes longer to reach its lower peak pressure and decay [134]. According to the characteristics, Saharan and Mitri [7] proposed optimized pressure-time profiles for ideal and nonideal detonation to investigate the rock dynamical response under blasting. It has been found that ideal detonation leads to more crushing around the borehole, which follows a large number of short radial cracks. In contrast, the non-ideal detonation results in a smaller crushing zone followed by a few long radial cracks. 290 6.6. Short delays Short delay blasting techniques are considered one of the potential methods for improving fracturing and fragmentation by overlapping waves and enhancing the tensile tail of the stress waves [104,135]. On the other hand, simultaneous initiation of massive blastholes in mines would release a large amount of energy that may potentially trigger a seismic event [6,136]. Thus, a delayed blast of those massive blastholes can overcome these problems. The nature of the delay blasting techniques is stress wave interaction emitted from different blastholes within the rock mass. If stress waves from two neighboring blastholes are effectively superimposed on each other, the final compressive wave is much longer than the single compressive wave, resulting in more fractured rocks and a rise of the final energy utilization [4]. The determination of the optimum delay time helps to optimize blast design in rock masses. Several numerical studies provided evidence that the interaction between the stress waves resulting from the detonation of neighbouring blastholes plays a great role in rock damage and fragmentation [104,137,138]. However, some studies also revealed that no distinct differences or high fragmentation improvements were caused by short delay blasting [139].

Air deck
Air deck is a novel technique and involves the replacement of some length of an explosive charge in a blasthole with air in a fully charged blasthole [140]. The air deck technique's motivation is the fact that the rock around a fully charged blasthole is broken too finely due to the absorption of much explosive energy. Thus, the air deck technique is developed to increase the explosion energy efficiency [140]. Air columns between adjacent individual charges provide space for detonation products expansion and shock wave collision, resulting in stress waves of lower amplitude and greater length propagating into rock masses [141]. As a result, the degree and uniformity of fragmentation, a decrease of average fragment size, and the amount of lump are increased, accompanied by the reduction of the specific consumption of explosives [140,142].
The mechanism involving air-deck in blasting practice has been further studied by numerical simulations. Liu and Katsabanis [143] performed numerical modelling on the air deck. It was found that a significant amount of the potential energy retained in the denotation was released and imparted to the rock mass due to the repetitive movement of the denotation products within the borehole with the air deck. As a result, a secondary loading wave was generated, increasing the amount of energy directed for a more even breakage of the entire rock material. The simulation results carried out by Lu and Hustrulid [144] also indicated that at least two stress waves are induced in the rock mass as air-deck exists.

Discontinuities
Rock often contains various scales of discontinuities, such as cracks, fissures, joints, cavities, bedding planes, faults, and other natural weaknesses. Discontinuities are frequently encountered in the stress wave travelling in rocks. When the stress wave encounters the joints, the transmitted portion may create some radial cracks in the rock block, while the reflected portion from the joint surface often causes some circumferential cracks [47,105]. Joint properties, such as width, spacing, orientation, stiffness, and infilling condition, has a great impact on the transmission and reflection characteristics of stress waves [9,105,145]. The joint sets potentially influence the destressing mechanism and should be considered in destressing pattern design. Drover et al. [6] suggested that collocating rows of destressing charges should be set along visible pre-existing discontinuities of suitable orientation to increase the gases penetrating and the interaction between the blast-induced fractures from adjacent destressing charges.

In situ stress field
The in situ stress around the blasthole has a great role in dynamic fracture extension and is essential for blasting design and application. Isotropic stress loading leads to an anisotropic fracture growth pattern in rock blasting. The existence of higher hydrostatic pressures significantly decreases the length of radial cracks and the size of the fractured zone, while the crushed zone near the blasthole has an increasing trend [34,128,146]. However, an anisotropic stress loading tends to guide cracks to extend preferentially along the direction of the maximum principal stress [93,94]. Besides, as shown in Fig. 6, the anisotropy of the fractured zone becomes much greater as the difference between two principal stresses increases [20,21,34]. A similar phenomenon is often observed in laboratory experiments [107]. The anisotropy of the failure pattern in the anisotropic in situ stress field has been reported in destress blasting applications [5].
Therefore, the effect of the in situ stress field on blasting should be considered in destress blasting design, which may affect the effectiveness of destress blasting.

Stemming
Proper stemming is beneficial for keeping more detonation energy and a longer detonation wave [147]. Consequently, the rock fracturing is more severe around borehole with stemming than that without stemming. Dally et al. [148] further confirmed the conclusion by performing two blasting experiments in two blastholes, one with stemming and another without stemming. It has been found that the maximum ratio of the crack radius at crack arrest to borehole radius was 7 for the model with the unstemmed charge, whereas with the stemmed charge this ratio was at least 37. On the other hand, the stemming column, as the least resistance path, also provides a path for the explosive energy ejection [5]. Therefore, the physical properties of stemming materials matter. To examine the effect of different stemming materials on the stress wave-induced damage in rock, a comparison of damage status for a cylindrical rock with different stemming materials, including intact rock, sand, and air, was made by Zhu et al. [28]. Simulation results showed that strong confinement introduced by rock stemming can efficiently block explosive products and produce more radial cracks than the other cases. Thus, proper application of stemming plugs may enhance explosive energy utilization to a good extent.

Boundary condition
The free boundary denotes a free surface that stress waves can partly reflect, while the transmit boundary allows full transmission of stress waves, in analogy to an infinitely large rock body. Therefore, the boundary condition should be carefully selected to obtain a realistic result in blasting modelling [31]. The effect of the boundary condition has been numerically studied. When a compressive stress wave propagates to a free surface, it will be reflected as tensile stress waves. As a consequence, rock fracturing close to the free surface will happen as long as the amplitude of the tensile wave reaches the dynamic tensile strength of the rock. These cracks developed well between the borehole and the free face are often named circumferential fractures or spalling, parallel to the free face [28,149]. In contrast, the transmit boundary does not cause any circumferential cracks. The effect of free boundary on spalling depends on several factors such as the geometrical form and amplitude of the tensile wave, the dynamic tensile strength of rocks, and the distance between the blasthole and the free surface [4,25,83]. The free surface can be viewed as one of the sources of tensile stresses in rock blasting, and the rocks are fractured efficiently when a nearby free surface exists.

Loading rate
The magnitude of the loading rate significantly affects the intensity of damage and failure of the rock masses [1]. Several numerical studies have focused on the effect of the loading rate on the fracturing during rock blasting [20,26,31]. Different waveforms with varying rising times are generated and applied on the borehole wall to study the effect of the loading rate. The rise time of explosion pressure to its peak varies strongly, depending on the blasthole diameterand confinement, rock strength, etc. [26,149]. It has been found that a high loading rate generally witnesses intensive crushing near the borehole and many short radial fractures, while a lower loading rate produces a relatively smaller crushed zone and fewer and longer radial fractures (Fig. 7). For a high strain rate, the crushed zone consumes too much Fig. 6. Effect of stress field on radial crack propagation in rock [20]. energy, while the rock in the relatively far-field is not sufficiently fractured, resulting in low energy efficiency [4]. Conversely, a lower loading rate always results in higher energy efficiency.

Conclusions
Static and dynamic approaches are two typically used approaches to modelling destress blasting of rock masses. Destress blasting in the static method (i.e., equivalent modelling approach) is achieved by modifying the mechanical properties and the postblast stress state of the rock masses in the specified destressed zone. On the other hand, the dynamic approaches can model the blasting process and give a more precise extent of fractured zones. For the dynamic approaches, the contribution of shock stress and denotative gas need to be considered, and their effect has been studied separately or jointly to reveal the blasting mechanism. Equation of state, pressure decay functions, predefined functions, etc., have been implemented within software packages to model the blasthole wall pressure. Lots of gas flow models have been developed to handle the interaction between the gas and fractures. Moreover, some innovative methods have been presented to solve the coupled effect of stress wave and gas pressure in rock blasting modelling. Modelling of rock masses can be classified into two categories: direct and indirect method. For the direct method, the microcracks can be directly modelled by the breakage of individual structural units or bonds. For the indirect method, a proper dynamic constitutive model should be selected to imitate the rock dynamic response in blasting practice. Commonly used dynamic constitutive models include JH2, JHC, RHT, and their modifications.
The fracture mechanism of blasting in rock mass has been widely studied, and similar fracture patterns are obtained. The shock wave causes a crushed zone near the boreholes and fractured zones with many radial cracks outside the crushed zone. The explosive gas would then penetrate the shock waveinduced fractures and further extend radial cracks. Shock waves are believed to have a major role in generating cracks compared with explosive gas penetration. According to numerous numerical validations, the destress efficiency of the commonly used blasting patterns is found to be overestimated, especially under high stress. Some novel destressing concept has been proposed, which provide a solution for optimizing the destress blasting design under high-stress conditions. The explosive in the rock masses is a complex process, and the blasting efficiency is dependent on numerous factors. Those factors which may affect the efficiency of destress blasting are briefly described. Understanding the mechanism of these factors should help improve the destress blasting design.

Ethical statement
We state that the research was conducted according to ethical standards.