Three-Dimensional Numerical Analysis of Blast-Induced Damage Characteristics of the Intact and Jointed Rockmass

This article reports numerical results investigating the damage evolution and spatial distribution characteristics of intact and jointed rockmass subjected to blast loading. The behaviors of rock material are described by the HolmquistJohnson-Cook (HJC) constitutive model incorporated in the finite element software LS-DYNA. Results indicate that the damage distribution shows a reverse S-shape attenuation with the increase of the distance from borehole, and a better goodness of fit with the Logistic function is observed. In the single-hole blasting of jointed rockmass, there are two types of regions around the intersection of borehole and joint in which the damage degree is quite different. The crushing damage develops in a Ψ-shape path along the joint. In the radial direction, the crushing damage and cracking damage of rock show different distribution forms with the increase of joint dip angle. As for the double-hole blasting, due to the superposition of the blast waves, the damage degree in the region between the two boreholes of intact rockmass is significantly large. For jointed rockmass, the joint has local enhancement or inhibition effect on the blast damage in the region between the two boreholes.


Introduction
Blasting technique has been widely employed in the areas of excavation and construction for many decades and it remains a popular method of rock breakage to the present day. The primary objective of drilling and blasting in practical activities is to achieve significant rock breakage without affecting the stability of retained rock-structure or its surrounding environment. However, due to the complexity and low controllability of blasting process, unplanned damage on the quality and integrity of the blasting projects is more or less induced, which has been drawing more and more attention [Zhang (2016)]. The damage area around the borehole is usually divided into three categories from near to far: crushed zone, cracked zone and elastic vibration zone [XaHykaeB (1980); Dai (2013)]. The crack problems were investigated by other researchers [Areias, Reinoso, Camanho et al. (2018); Areias, Msekh and Rabczuk (2017); Areias, Rabczuk and Camanho (2014)]. Also, many numerical methods for crack propagation simulation have been increasingly developed [Rabczuk and Belytschko (2007); ; ; Areias and Rabczuk (2017)]. In addition, the methodology for stochastic modelling of the fracture in polymer/particle nanocomposites [Hamdia, Silani, Zhuang et al. (2017)] may provide a new inspiration for rock cracking simulation, and a sensitivity analysis toolbox provided by Vu-Bac et al. [Vu-Bac, Lahmer, Zhuang et al. (2016)] can help determine the key parameters in blast damage analysis. It is all known that extensive research activities in the field of blast-induced damage have been carried out in recent years. Lu et al. [Lu, Leng, Chen et al. (2016)] explored the sizes of fractured zone and crushed zone and found that the size of the crushed zone ranged from 1.2 to 5.0 times the borehole radius considering the influence of the in-situ stress, and the size of the fractured zone was greatly affected by in-situ stress. Based on blasting tests, Esen et al. [Esen, Onederra and Bilgin (2003)] introduced a dimensionless parameter called the crushing zone index (CZI) and showed that the radius of crushing was a function of the CZI and the blasthole radius. Onederra et al. [Onederra, Furtney, Sellers et al. (2013)] developed the hybrid stress blasting model and stated that the model could adequately predict both the shape and scope of the damage zone considering the influence of point-of-initiation and free boundary. Using the stress intensity factor (SIF), Sim et al. [Sim, Cho and Song (2017)] evaluated the two-dimensional fragmentation zone of rock induced by gas pressure and concluded that the fragmentation zone in weathered rock was wider than in hard rock. In the research of crack evolution mechanism, Zhou et al. [Zhou, Zhuang and Rabczulk (2018) ;Zhou, Zhuang, Zhu et al. (2018); Zhou, Rabczulk and Zhuang (2018)] introduced a phase-field modeling of quasi-static and dynamic crack propagation in poroelastic media. As is well known, rockmass is usually broken up by joints into rock elements which are continuous and may be regarded as elastic bodies. The properties of rockmass are mainly determined by the properties of the rock elements and the joints. The propagation of blast wave in rockmass is significantly affected by the joint-filling material, joint thickness, dip angle, and the initial in-situ stress [Zhang, Peng, Fan et al. (2016); Qu and Liu (2015); Wei, Zhu, Bai et al. (2016)]. Fu et al. [Fu, Chen, Wen et al. (2018)] developed a singular boundary method (SBM) and applied it to analyse the wave propagation at low and moderate wave numbers in periodic structures. Zhao et al. [Zhao, Zhang and Ranjith (2017)] established a constitutive model considering dynamic compression and tensile failure to analyze the fracture propagation in coal using LS-DYNA. They reported that the transmission coefficient of stress wave at bedding plane decreased with the distance from borehole. Chen et al. [Chen and Zhao (1998)] analyzed the influence of joint on the propagation of stress wave based on the universal distinct element modelling, finding that the joint could attenuate the stress wave. Zeng et al. [Zeng, Wang, Sun et al. (2018)] proposed the geological mechanical models for the intact and jointed rockmass with the Froude Similarity principle, indicating that the peak velocity, and the transmission and reflection coefficients increased with the increase of joint dip angle. Ma et al. [Ma and An (2008)] applied the Johnson-Holmquist model to simulate both compressive crushed zone and tensile radial fractures. They claimed that loading rate had a significant effect on the fracture pattern and the existence of charge holder helped to control the initiation and propagation of fractures.
To date, most studies on the blast excavation of rockmass are based on two dimensional (2D) analysis. Although the 2D method can simplify the calculation process, it cannot well illustrate the spatial distribution of blast damage. In addition, it is lack of a complete study of the effect of joint skew intersection with borehole on the blast damage. Due to the widespread distribution of joints, the phenomenon of joint skewing with borehole is inevitable in practice, which must be paid enough attention to. In this study, 3D numerical simulations of single-hole and double-hole blasting are carried out for intact rockmass and jointed rockmass. The features of damage distribution of single-hole blasting in intact rockmass are investigated. Then, the propagation of blast wave and the evolution of rock blast damage are analyzed in detail considering the influences of joint with different dip angles. The characteristics of rock damage under double-hole blasting are also discussed. Finally, some useful conclusions and remarks are given.

Numerical simulation of blast damage
The finite element code LS-DYNA with an explicit integration scheme is used for the modeling of three-dimensional blast. It can deal with the nonlinear dynamic response of structures. Generally, the Lagrangian algorithm is applied for the solid materials, while the arbitrary Lagrangian-Eulerian (ALE) algorithm is ideally suited to simulate the problems of fluid and gas flow. Additional details on fluid-structure interaction of finite method between sub-grids, such as Lagrangian-Lagrangian and Lagrangian-Eulerian coupling can be found in Refs. Anghileri et al. [Anghileri, Castelletti and Tirelli (2005); Alia and Souli (2006); Ozdemir, Moatamedi, Fahjan et al. (2009);Rabczuk, Gracie, Song et al (2010)]. To make use of the advantages of both the algorithms, the multimaterial fluid-structure interaction method is adopted in this study. The Lagrangian algorithm is adopted for the rock and stemming, and the ALE algorithm is adopted for the air and explosive.

Material model 2.1.1 Rock
Rock is typically non-homogeneous in composition and may possess flaws such as cracks and voids. Strictly speaking, the linear-elastic model can only describe its initial deformation. Under blast loading, the flaws may cause elastoplastic damage in the rock medium. The elastoplastic damage is irreversible and dissipates the blast-wave energy. This study adopts the Holmquist-Johnson-Cook (HJC) model which has been specifically developed for rock-like materials, e.g., concrete, subjected to high pressures, large strains, and high strain rates [Holmquist, Johnson and Cook (1993); Ren, Wu, Fang et al. (2017)]. As seen from Fig. 1, this model has three components: strength criterion with strain rate effect, nonlinear hydrostatic compression, and damage evolution with plastic deformation. The equation of the yield surface (see Fig. 1(a)) is [Holmquist, Johnson and Cook (1993) where σ*=σ/fc is the normalized equivalent stress, σ is the actual equivalent stress and fc is the quasi-static uniaxial compressive strength; P*=P/fc is the normalized pressure, where P is the actual hydrostatic pressure; where ε  is the actual strain rate and 0 ε  =1.0 s -1 is the reference strain rate; D is the damage parameter; A, B, N and C are the constants determined by the material strength. The equation of state is presented in Fig. 1(b) [Holmquist, Johnson and Cook (1993)]. The pressure-volume response is separated into three parts. The first part is the linear elastic stage. The second part is the transition stage, which is characterized by the progressive compression of voids and plastic deformation. The third part is the compaction stage, which is usually expressed as where K1, K3, K3 are material constants. Modified volumetric strain is defined as In Fig. 1(b), Pcrush and μcrush are the crushing pressure and crushing volumetric strain, Plock and μlock are the locking pressure and locking volumetric strain, ρ0 and ρ are the initial and current densities, respectively.
The damage model of HJC describes accumulation damage from both equivalent plastic strain and plastic volumetric strain (see Fig. 1(c)). Its expression is given as [Holmquist, Johnson and Cook (1993) where ΔεP and ΔμP are the increments of the equivalent plastic strain and plastic volumetric strain, respectively; The term of f f P P ε µ + is the plastic strain at material fracture under a constant pressure; T*=T/fc is the normalized maximum tensile hydrostatic pressure where T represents the actual maximum tensile hydrostatic pressure; D1 and D2 are damage constants. It is noted that the damage parameter D varies in the range of [0, 1]. Based on impact tests and theoretical analysis, Fang et al. [Fang, Kong, Wu et al. (2014)] proposed a set of constitutive parameters of the HJC model for limestone, as listed in Tab. 1. These parameters are adopted in this study to simulate the blast damage in rockmass.

Joint
The influence of joint on blast damage is considered in this study. The filled material and joint thickness remain unchanged during the simulation. A soft rock is considered as the filling material which can be characterized by kinematic hardening model (MAT_PLASTIC_KINEMATIC) [Zhang, Peng, Fan et al. (2016)]. The joint thickness d is set as 1.0 cm. The main parameters of the joint are as follows: density ρ=1160 kg.m -3 , Young's modulus E=20.0 GPa, Poisson's ratio ν=0.3 and shear modulus G=2.5 GPa [Zhang, Peng, Fan et al. (2016)].

FE model
Owing to symmetry, a half model of single-hole blasting in intact rockmass is established, as shown in Fig. 2. Its dimensions are 140 cm (length)×70 cm (width)× 200 cm (depth). The borehole is placed in the middle of the rock with a radius rb of 2 cm and a length of 175 cm. The stemming length is 25 cm, and the stemming material is assumed to be the same as the rock. The coupling charge mode and the instantaneous ignition are adopted. The free boundary condition is set for the top surface of the model, and the non-reflecting boundary condition is set for the side and bottom surfaces. The free boundary represents a regular free surface in which stress wave will reflect back; the non-reflecting boundary means that the stress wave transmits the boundary and is not allowed to reflect back, which essentially is an analog of the semi-infinite rockmass. The normal displacement is constrained on the symmetry plane (i.e., XOY plane). In the fluid-structure interaction algorithm, the coupling medium is air, of which the radius is taken as 10 times the borehole radius. The total number of eight-node SOLID164 elements is about 600 000.  An oblique joint is added in the model for the simulation of jointed rockmass blasting, as shown in Fig. 3. The contact between joint and rock is simulated through the co-nodes of the elements. In the model, α is the angle between the joint and the X-axis, and O is the intersection point of the joint and the borehole. To present the distribution of blast damage in the radial direction of borehole, four cross-sections (S0, S1, S2 and S3) are selected (see Fig. 3(b)). On the symmetry plane of the model, O, O1, O2 and O3 are the intersection points of each cross-section with the joint, and the distances from these points to the borehole axis are 0 cm, 10 cm, 20 cm and 30 cm, respectively. Similar model setup is adopted for the simulation of double-hole blasting in intact and jointed rockmass, as shown in Fig. 4. Fig. 5 shows the blast damage contours of intact rockmass along the axial and radial direction of borehole. It can be seen that a crushed zone is detected near the borehole. In the crushed zone the rock completely fails (D=1). Outside the crushed zone, the cracked zone emerges with radial fractures. The radius of the crushed zone and cracked zone are denoted as R1 and R2, respectively. For the convenience of description, a circle with the radius of l is drawn from the center of the borehole. At the periphery of the circle, the rock damage is extracted from the points at an interval of 9° (with a total number of 40 points), as shown in Fig. 6(a). Then, the average of these values is taken as the damage value of rock at the distance of l from the borehole. Fig. 6(b) shows that the rock damage demonstrates a reverse S-shape attenuation as a function of the distance l from the borehole center. where A1, A2, B and l0 are constants. The best fitting result is obtained and shown in Fig. 6(b). The coefficient of determination R 2 is found to be 0.997, which indicates a promising fitting result. It can be concluded that the Logistic function is suited to characterize the distribution of rock damage. According to Eq. (9), when D=1, the radius of the crushed zone is 7 cm. The ratio of the crushed zone radius R1 to the charge radius r is 3.5. When D is less than 0.05, the damage is negligible. Thus, the ratio of the cracked zone radius R2 to the charge radius r is about 14.5. the damage gradually decreases with the increase of distance from the borehole. It can be seen that the ratio of the crushed zone radius to the charge radius is about 3.6, which agrees fairly well with that from this study. Based on a large amount of experimental data, Esen et al. [Esen, Onederra and Bilgin (2003)] proposed a theoretical formula as follows, which can be used to predict the radius of rock crushed zone:

Single-hole blasting in intact rockmass
where CZI represents the crushing zone index and is expressed as where Pb is the borehole pressure, K is rock stiffness, σc is unconfined compressive strength.
Under the condition of cylindrical coupling charge, Pb can be obtained from the following formulas [Dai (2013)]: where P0 is the detonation pressure of explosive (Pa); ρr, ρe, Cr, Ce and γ are, respectively, the density of rock (kg.m -3 ), the density of explosive (kg.m -3 ), the acoustic velocity of rock (m.s -1 ), the detonation velocity of explosive (m.s -1 ), and the adiabatic expansion coefficient of detonation products. The typical value for γ is 3. It is assumed that the material within the crushing zone is homogeneous and isotropic, so the rock stiffness K can be given by Esen et al. [Esen, Onederr and Bilgin (2003)]: where Ed is the dynamic Young's modulus and νd is the dynamic Poisson's ratio. If the Ed and νd are not available, it can be estimated from the following relationship [Dai (2013); Eissa and Kazi (1988) where Est=2(1+νst)G is the static Young's modulus, νst is the static Poisson's ratio, G is the shear modulus. For the rock considered in this study, Cr =3430 m . s -1 and νst =0.26 [Dai (2013)]. Other parameters are defined in Section 2.1. From Eqs. (10)~(16), it is found that the ratio of crushed zone radius to charge radius is about 4.2. Tab. 2 compares the radii of crushed zone and cracked zone. Obviously, the computed results are very close to those of other literature, which indicates that the numerical results in this study are reliable.  Fig. 8 shows the damage contour on the symmetry plane of the jointed rockmass when the joint dip angle α equals to 45°. Near the joint, the damage zones can be categorized as: (i) the "damage inhibition zone" (region I), and (ii) the "damage enhancement zone" (region II). Due to the dispersion of blast wave, the rock damage decreases and a damage inhibition zone emerges. On the other hand, the damage enhancement zone becomes visible because of the energy convergence of blast wave. In the damage enhancement zone, the crushing damage develops along the joint first and then gradually extends along the axis direction of the borehole, and finally connects with the crushed zone at the far end of the borehole. The intersection lines between the cross-sections (S0, S1, S2 and S3) and the model symmetry plane are taken as the monitoring lines (see Fig. 3(b)). Fig. 9 shows the relationship between the peak particle velocity (PPV) and the distance from the borehole center to the monitoring position. It can be seen that the PPVs on the cross-sections S1, S2 and S3 decrease sharply at the points O1, O2, and O3 compared with the PPV on the cross-section S0. This is due to energy dissipation as the blast wave propagates across the joint.

Figure 9
: PPV attenuation curves on the symmetry plane of jointed rockmass Fig. 10(a) illustrates the pressure contour of the cross-section S1 at t=0.1 ms. It is seen that the reflection and transmission of the pressure wave happen at the joint. Fig. 10(b) demonstrates the pressure and displacement-time history curves of the rock elements on the left and the right sides of point O1. Obviously, when t=1 ms, the pressure and displacement of the elements tend to be constant, which indicates that the calculation results have converged. The peak pressure and maximum displacement of the left element (No. 84703) are much larger than that of the right one (No. 491828) at point O1. This means that the energy of blast wave is greatly reduced after the wave passing through the joint. In addition, the pressure wave of right element at point O1 has two pressure peaks, which is mostly attributable to the repetitive reflections of blast wave in the joint.
(a) Pressure contour (t=0.1 ms) (b) Pressure and displacement-time history curves Figure 10: Pressure contour and dynamic response history curves on cross-section S1 Fig. 11 shows the final damage contours on the symmetry plane of jointed rockmass with different joint dip angles. The damage distribution of jointed rockmass is similar to that of intact rockmass (see Fig. 5(a)) when α is small (α=15°). As the dip angle rises, the damage distribution approximates to the simulation observations shown in Fig. 8.
Fig. 12 presents the blast damage distributions on the cross-sections S0, S1, S2 and S3 under different dip angles. It is seen that: (i) when the dip angle is small (for example, α=15°), the damage on cross-section S0 shows nearly uniform evolution in all directions.
With the increase of α, the range of cracked zone increases in the X direction, and the crushed zone develops prominently in the Z direction; (ii) on the cross-sections S1 and S2, the range and degree of rock damage on the right side of the joint decrease with the increase of dip angle. The damage increases in the region between the joint and the borehole due to the refection of blast wave at the joint. However, when the distance from the joint to borehole is larger (such as S3), this feature becomes indistinct; and (iii) when α rises to 60°, the shape of fracture is obvious, and the developments of main fractures in breadth and depth are significant. Fig. 13 shows the blast damage distribution on cross-section S4, which is parallel to the joint and close to the interface between rock and joint (see Fig. 3(b)). It is seen that the crushing damage on section S4 has a unique evolution path. Taking α=60° as an example, the crushing damage develops mainly towards the right side of the borehole and extends along a Ψ-shape path. This feature is more obvious under a larger dip α.
Figure 13: Damage distribution in the joint

Double-hole blasting in rockmass
The numerical results of double-hole blasting in intact and jointed rockmass are discussed in this section. The joint dip angle α of 45° is simulated, as shown in Fig. 4. Simultaneous ignitions of the explosives in the boreholes are assumed. It is seen from Fig. 14(a) that the degree of damage in the region between the two boreholes of intact rockmass is significant due to the superposition effect of the wave. The blast damage around a single borehole in the jointed rockmass still follows the way of damage evolution described above. However, the rock damage between the two boreholes shows some differences locally, as illustrated in region III and region IV of Fig. 14(b). For convenience, two lines L1-L1'and L2-L2' in Fig. 4 are selected. The damage distributions in the longitudinal direction, i.e. along line L1-L1' and line L2-L2', are shown in Fig. 15. It is noticed that the damage along line L2-L2' presents a "double-peak" feature, and the two peaks are almost all of the same magnitude. The damage value of line L2-L2' is greater than that of line L1-L1' in region IV (y=0.4~0.8 m and 1.2~1.6 m), but less than that of line L1-L1' in region III (y=0.8~1.2 m). These results indicate that the oblique joint intersecting with double boreholes has the dual effects of enhancing or inhibiting the blast damage compared to the double-hole case in the intact rockmass. distribution of blast damage were analyzed. Based on the numerical results, the following conclusions are drawn: (1) Under the single-hole blasting, the damage around the borehole in intact rockmass shows a reverse S-shape attenuation with the increase of the distance from the borehole center, which can be well described by the Logistic function.
(2) On the symmetric plane of the model, the crushed damage evolves first along the joint and then towards the direction of the borehole axis. Near the joint surface, the crushing damage evolves along a Ψ-shape path. The variation of joint dip angle appreciably affects the evolution path and space distribution of damage.
(3) For the double-hole blasting scenario, the rock damage in the region between the two boreholes increases significantly in the intact rockmass. However, in the jointed rockmass, the local enhancement and inhibition effects of joint on the rock damage may exist simultaneously in the region between the two boreholes.