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

Zhiliang Wang ,Youpeng Huang and Feng Xiong

Abstract: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 Holmquist- Johnson-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.

Keywords:Rockmass,joint,blast damage,spatial distribution,numerical analysis.

1 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);Rabczuk,Zi,Bordas et al.(2010);Rabczuk,Bordas and Zi (2010);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 newinspiration 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 e 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 skewintersection 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.

2 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.

2.1 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σ*=σ/fcis the normalized equivalent stress,σis the actual equivalent stress andfcis the quasi-static uniaxial compressive strength;P*=P/fcis the normalized pressure,wherePis the actual hydrostatic pressure;=/is the dimensionless strain rate,whereis the actual strain rate and=1.0 s-1is the reference strain rate;D is the damage parameter;A,B,N and C are the constants determined by the material strength.

Figure1:HJC constitutive model

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 aswhereK1,K3,K3are material constants.Modified volumetric strain is defined as

In Fig.1(b),Pcrushandμcrushare the crushing pressure and crushing volumetric strain,Plockandμlockare the locking pressure and locking volumetric strain,ρ0andρa(bǔ)re 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 ΔεPand ΔμPare the increments of the e equivalent plastic strain and plastic volumetric strain,respectively;The term ofis 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 e HJC model for limestone,as listed in Tab.1.These parameters are adopted in this study to simulate the blast damage in rockmass.

Table1:Parameters of HJC model for limestone

2.1.2 Air

The air is usually described by the ideal gas law (i.e.,Gamma law).This polytropic equation of state is given as Alia et al.[Alia and Souli (2006)]

For ideal gas,C0=C1=C2=C3=C6=0,C4=C5=γ-1.Then,Eq.(6) can be simplified as

where the specific volumeV=ρ/ρ0-1,ρ0andρa(bǔ)re the initial and current densities of air;Eis the specific internal energy;γis the coefficient of adiabatic expansion.For the diatomic molecules comprising air,γ=1.4.At timet=0,the initial densityρ0=1.293×10-3g.cm-3and the initial specific internal energyE0=0.25 MPa [Alia and Souli (2006)].

2.1.3 Explosive

The “Jones-Wilkins-Lee” (JWL) equation of state is here used to model the pressure generated by chemical energy in an explosive as below [Alia and Souli (2006)]:

whereA,B,R1,R2andωare JWL parameters.The computational parameters for explosive are as follows [Zhao,Zhang and Ranjith (2017)]:the density of explosive isρe=1000 kg.m-3,the detonation velocity isCe=4000 m.s-1,A=541 GPa,B=9.4 GPa,R1=4.5,R2=1.1,ω=0.35,E0=8 GPa,and the Chapman-Jouget pressurePCJ= 6.1 GPa.

2.1.4 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 thicknessdis set as 1.0 cm.The main parameters of the e joint are as follows:densityρ=1160 kg.m-3,Young’s modulusE=20.0 GPa,Poisson’s ratioν=0.3 and shear modulusG=2.5 GPa[Zhang,Peng,Fan et al.(2016)].

2.2 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 radiusrbof 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.

Figure2:Numerical model of single-hole blasting in intact rockmass

Figure3:Numerical model of single-hole blasting in jointed rockmass

Figure4:Numerical models of double-hole blasting in intact and jointed rockmass

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 e 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,S2and S3) are selected (see Fig.3(b)).On the symmetry plane of the model,O,O1,O2and O3are 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.

3 Results and discussion

3.1 Single-hole blasting in intact rockmass

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 e crushed zone and cracked zone are denoted asR1andR2,respectively.

For the convenience of description,a circle with the radius oflis 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 oflfrom the borehole.Fig.6(b) shows that the rock damage demonstrates a reverse S-shape attenuation as a function of the distancelfrom the borehole center.

Figure5:Damage contours of intact rockmass after single-hole blasting

Figure6:Distribution law of rock blast damage with the increase of distance from borehole

According to this characteristic,the Logistic function is selected and applied to fit it.The fitting formula is given as follow:

whereA1,A2,Bandl0are constants.

The best fitting result is obtained and shown in Fig.6(b).The coefficient of determinationR2is 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),whenD=1,the radius of the e crushed zone is 7 cm.The ratio of the crushed zone radiusR1to the charge radiusris 3.5.WhenDis less than 0.05,the damage is negligible.Thus,the ratio of the cracked zone radiusR2to the charge radiusris about 14.5.

Figure7:Blast damage zones of PMMA sheet in single-hole blasting [Zhang (2016)]

Zhang [Zhang (2016)] investigated the crushing and cracking phenomenon of rock based on the blast test of PMMA sheet.The experimental observation is shown in Fig.7 where 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:

whereCZIrepresents the crushing zone index and is expressed as

wherePbis the borehole pressure,Kis rock stiffness,σcis unconfined compressive strength.

Under the condition of cylindrical coupling charge,Pbcan be obtained from the following formulas [Dai (2013)]:

whereP0is the detonation pressure of explosive (Pa);ρr,ρe,Cr,Ceandγ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 stiffnessKcan be given by Esen et al.[Esen,Onederr and Bilgin (2003)]:

whereEdis the dynamic Young’s modulus andνdis the dynamic Poisson’s ratio.If theEdandνdare not available,it can be estimated from the following relationship [Dai (2013);Eissa and Kazi (1988)]:

whereEst=2(1+νst)Gis the static Young’s modulus,νstis the static Poisson’s ratio,Gis the shear modulus.For the rock considered in this study,Cr=3430 m.s-1andν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.

Table2:The radii of crushed zone and cracked zone from different researches

3.2 Single-hole blasting in jointed rockmass

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.

Figure8:Damage contour of jointed rockmass after single-hole blasting

The intersection lines between the cross-sections (S0,S1,S2and 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,S2and S3decrease sharply at the points O1,O2,and O3compared with the PPV on the cross-section S0.This is due to energy dissipation as the blast wave propagates across the joint.

Figure9:PPV attenuation curves on the symmetry plane of jointed rockmass

Fig.10(a) illustrates the pressure contour of the cross-section S1att=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 e rock elements on the left and the right sides of point O1.Obviously,whent=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 O1has two pressure peaks,which is mostly attributable to the repetitive reflections of blast wave in the joint.

Figure10: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.

Figure11:Damage contours of jointed rockmass with different joint dip angles in single-hole blasting

Fig.12 presents the blast damage distributions on the cross-sections S0,S1,S2and S3under different dip angles.It is seen that:(i) when the dip angle is small (for example,α=15°),the damage on cross-section S0shows 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 S1and 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 S4has 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α.

Figure13:Damage distribution in the joint

3.3 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 e 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.

Figure14:Damage contours of intact and jointed rockmass in double-hole blasting

Figure15:Comparison of damage distribution at the position of monitoring lines

4 Concluding remarks

Three-dimensional simulations of rock damage during single-hole and double-hole blasting were performed.The effects of oblique joint on the evolution and spatial 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.

Acknowledgments:This study was financially supported by the National Natural Science Foundation of China (51579062,51379147),which is gratefully appreciated.