Three dimensional numerical rock damage analysis under blasting
load
Ozgur Yilmaz ⇑,Tugrul Unlu
Department of Mining Engineering,Bülent Ecevit University,Zonguldak,Turkey
a r t i c l e i n f o Article history:
Received 5March 2012
Received in revised form 30January 2013Accepted 1July 2013
Available online 27July 2013Keywords:
Blasting damage Numerical analysis Strain rate effect
Anisotropic in situ stresses
a b s t r a c t
In this study,the behaviour of rock mass subjected to blasting load is investigated using three dimen-sional finite difference numerical modelling.In the analyses,Mohr-Coulomb failure criterion has been used for the characterisation of the rock mass strength.Stresses acting on the borehole boundary have been simulated by an exponential function which reaches its maximum within a short time and then falls to zero value in a considerable period.The strain rate effect on the mechanical properties of rock material has also been taken into account in the analyses.Different explosive and site conditions have been stud-ied to investigate the effects of loading rate and anisotropic high in situ stresses on blasting performance and blast induced damage zones.Results have shown that the most efficient explosive in
rock blasting will be the one with low frequency content but with a sufficiently high borehole wall pressure.In addi-tion,it has been verified that the directions and the magnitudes of major principle stresses affect the development of the crack zone around the borehole.Finally,it has been seen that proposed equation for the dynamic compressive strength of rock material fits very well to general suggestions.
Ó2013Elsevier Ltd.All rights reserved.
1.Introduction
In mining and construction operations,the use of explosives is probably the most widely used means of crushing rock,as well as the most cost-effective rock excavation method.At the under-ground openings which are excavated by the drilling and blasting method,predicting the blast induced damage zone is important for long term stability.However,rock crushing process is a very fast process under blasting,so it is extremely difficult to study rock blasting failure mechanism.Therefore,a clear description of a com-plete and accurate rock blasting damage theory is highly difficult phenomenon.
In the analysis of rock behaviour in underground excavations subjected to dynamic loads such as blasting,it is vital to know the dynamic mechanical properties of rock material.In addition,blast induced
damage zones around a blasthole are affected by not only the high anisotropic in situ stresses but also the distance between the hole and the free surface.On the other hand,investi-gation of rock blasting and its effects by scaled or full size experi-ment is highly expensive and inconvenient.In conjunction with the development of high-speed computer algorithms,numerical simulation has become one of the major means to study the dynamics of rock blasting.
The objective of this study is to recommend a simple blast damage analysis model by using three-dimensional numerical technique.A dynamic geo-technical software FLAC 3D (v.2.1)(Itasca,2003)has been used in order to conduct a numerical simulation of the detonation of a finite cylindrical charge.
2.Numerical modelling of blast induced damage analysis Two basic forms of energy are released when an explosive re-acts.The first type of energy is called as shock energy and the sec-ond one is called as gas energy (Kutter and Fairhurst,1971).This released energy generates a stress wave loading travels outward from the explosion cavity followed by a longer duration gas pres-surization loading.Recent studies have revealed that the stress wave is responsible for the initiation of the crushing zone and the surrounding radial fractures,while the gas pressure further ex-tends the fractures.Numerical results obtained by some research-ers indicated that modelling only the stress wave propagation could give reasonable prediction of rock mass response to blasting load (Donze et a
l.,1997;Hao et al.,2002;Cho and Kaneko,2004;Jong et al.,2005;Lee et al.,2005;Ma and An,2008;Shin et al.,2011).Thus,blasting damage to rock mass due to the pressure of the expanding gases is not considered in this investigation.
Blasting is an impulsive load problem requiring three-dimen-sional consideration.Shock waves propagate spherically and a complex interaction between detonation and explosion gas pres-sures induces strains in surrounding rocks.Therefore,main con-cern should be given to the modelling of borehole wall pressure.Calculation of the borehole wall pressure and the other numerical modelling details are given below.
0886-7798/$-see front matter Ó2013Elsevier Ltd.All rights reserved./10.1016/j.tust.2013.07.007
⇑Corresponding author.Tel.:+9037225740101525;fax:+903722574023.
E-mail addresses:oyilmaz@ (O.Yilmaz),unlutugrul@hotmail (T.Unlu).
2.1.Calculation of the borehole wall pressure and its time history
An important assumption in this study is that the damage zone radius is dependent upon the borehole
wall pressure.Borehole wall pressure describes the expansion work of the explosive during the rock breakage process and it is the most important information in the evaluation of the explosive performance and the prediction of blasting results.Despite the importance of this parameter,direct measurements of the borehole pressure have not adequately been carried out due to the lacking of feasible methods,instead,various empirical formulas or detonation theories are used to estimate it (Esen et al.,2003).The detonation pressure needs to be calculated in order to estimate the pressure acting on the borehole wall since detonation is the beginning phase of the explosive process.Deter-mination of the detonation pressure is very complex.Fickett and Davis(1979)and Henrych(1979)proposed the following simpli-fied equation,which is derived using the equation of state of ideal gases,to evaluate of maximum borehole wall pressure in an explo-sive wave propagating from a cylindrical charge:
P d¼q e V2d=ð1þcÞð1Þwhere P d is the detonation pressure(Pa),V d is detonation velocity (m/s),q e is density of the unreacted explosive(kg/m3)and c is spe-cific heat ratio.The assumption c=3give the well-known expres-sion of the detonation pressure for many explosives(Persson et al.,1993):
P d¼q e V2d=4ð2ÞThe detonation pressure generates a gas pressure in the reacted portion of the explosive column.The gas pressure,often called explosion pressure,is the pressure that is exerted on t
he borehole walls by the expanding gases after the chemical reaction has been completed.Explosion pressure is approximately one-half the deto-nation pressure(Konya and Walter,1991).The explosion pressure, thereby,can be estimated through following equation:
P e¼P d=2¼q e V2d=8ð3ÞThe explosion pressure expresses the pressure of a fully coupled charge completelyfilling the blasthole.The explosive pressure, thus,also called as full coupled borehole pressure.In thefield, however,the decoupling technique is used,which refers to leaving an empty space between an explosive column and the blasthole wall.For considering the decoupling effect,the following calibra-tion equation is used(Nie and Olsson,2000):
P b¼P e r2k
c
ð4Þ
where P b is the borehole wall pressure(Pa),r c is the ratio of explo-sive diameter to blasthole diameter(coupling ratio=d c/d h)and k is explosive’s adiabatic expansion constant.Assuming k=1.5for aver-age adiabatic exponent,Eq.(4)gets the simplest form as follows (Persson et al.,1993),
P b¼P eðd c=d hÞ3ð5ÞConsequently,the borehole wall pressure acting on the bound-ary of the borehole can be estimated as follows:
P b¼ðq e V2d=8Þðd c=d hÞ3ð6ÞThis pressure is oriented radially outward from the wall of the explosive charge.If the explosive charge was in intimate contact with the borehole wall(fully coupled conditions),the explosion pressure(Eq.(3))would be the borehole wall pressure.In this study,for simplicity,the explosive charge diameter is assumed to be equal to the borehole wall diameter.
Since the actual dynamic pressure acting on the borehole wall varies with the time,time history of the dynamic pressure has to be considered to accurately simulate the dynamic analysis for rock blasting.
There are several ways to obtain the pressure–time profile of the explosion pressure.Numerically,these borehole pressure pro-files can be approximated by one of three procedures(Sarahan and Mitri,2008):(1)Equation-of-State(EoS)of the expanding gases products,(2)direct input of dynamic pressure as a function of time or,(3)pressure-decay functions.
The EoS describes material behaviour in the high-rate intense pressure environment and the equation relates different material quantities as a single function regardless of prior history of defor-mation.The equation-of-state of the expanding gases products uses the ideal detonation theory.The accuracy of us
ing EoS is therefore questionable since the detonation process for many explosives is non-ideal.There are also some problems in using the EoS of expanding gases such as choosing the accurate EoS and determining the required parameters of preferred EoS.
The use of Gaussian function and triangular load function has also been attempted to approximate measured dynamic-pulse load.These procedures,however,are not close to the physical char-acteristics of the dynamic load and hence carry no physical mean-ing.The Gaussian function is mainly introduced to avoid numerical errors associated with the application of very high pressure,which is in the order of GPa,in a very short duration in the order of micro-seconds(Sarahan and Mitri,2008).
Despite the fact that the physical processes associated with the detonation of a column of explosives have been studied for various decades,most of the blasting applications use empirical attenua-tion laws to represent the stress wavefield,or apply discrete spherical source models to simulate continuous cylindrical detona-tion(Hustrulid,1999).Several stress wave patterns in the form of decay or window functions have also been employed to replicate the measured stress wave records in the close vicinity of an explo-sive source.
In this study,the pressure function which was originally intro-duced by Starfield and Pugliese(1968)and
modified by Jong et al. (2005)is employed as a stress wave function to replicate the time history of the dynamic pressure.According to this function that has been generally applied,stress wave is expressed as follows:
P t¼4P bðeÀb t=
ffiffi
2
p
ÀeÀ
ffiffi
2
p
b tÞð7Þwhere P t is the time history of the borehole pressure,P b is the bore-hole pressure(Pa),b is da
mping factor(1/s)and t is time(s).Explo-sively generated pressures are generally specified by rise rime,t r and a corresponding peak magnitude.Damping factor is determined according to peak pressure rising time.According to Eq.(7),the peak pressure occurs at the time t r¼À
ffiffiffi
2
p
lnð1=2Þ=b.So we can determine the damping factor with respect to peak pressure rise time as follows:
b¼À
ffiffiffi
2
p
lnð1=2Þ=t rð8ÞIn this study,explosives were chosen according to their detona-tion velocity(V d),in orde
r to generate ideal and non-ideal detona-tion.For this purpose,an emulsion for ideal detonation case and ANFO for non-ideal detonation have been chosen.Vanbrabant et al.(2002),conducted detonation velocity measurements for 10m long charges of ANFO and BLASTEX(emulsion)and provided values of4052m/s for ANFO and5582m/s for the emulsion.Se-lected explosive properties and calculated borehole wall pressures are given in Table1.
Many publications cite explosive peak pressure attainment time in the order of several milliseconds but there are also in borehole measurements suggesting that the peak pressure attainment time is between20and150l s,depending on explosive type and con-finement(Sarahan et al.,2006).Generally,harder rocks generate
O.Yilmaz,T.Unlu/Tunnelling and Underground Space Technology38(2013)266–278267
high peak magnitude with short rise time while the contrary is true for softer formations(Simha,1996).
In this study,the rise time of explosion pressure to its peak is varied between10and1000l 10,50,150and1000l s)to investigate the dynamic damage zone processes.The peak pressure P d is kept at constant value of1.6and4.9GPa according to explo-sive types.Fig.1shows the pressure–time histories used for the numerical simulations.As can be seen in Fig.1,observed periods (T)of these pres
sure waves are different in all of the considered cases.According to observed periods one can obtain the corre-sponding angular velocities(x)as:
x¼2p
observedperiod radians
second
ð9Þ
Based on this,frequency of vibration(cycles/s)can be written in terms of the angular velocity as:
f¼x=2pð10ÞThe terms like angular velocity,frequency and period are im-ported from the point of the shape of applied pressure waves.It is critical to know that how much time the explosive pressure acts to the borehole wall.Moreover,the ratio b/x is directly affects the shape of the wave(especially fall time).In this work,it is kept at constant value of1.4.Properties of the considered pressure wave loading rates and the other parameters)are listed in Table2.2.2.Effects of strain rate
The strain rate effect on the strength of various engineering materials such as concrete,mortar rock,cer
amics,silicon carbide, polymericfibre,and composite materials,has become an impor-tant factor in both the material model and the design of structures that may experience high strain rate in a range of applications when impact or blast loading is involved(Li and Meng,2003). Extensive experimental results have shown that there is an appar-ent increase of the dynamic strength when the engineering mate-rial is subjected to high strain rate(Lankford,1981;Malvern and Ross,1985;Bischoff and Perry,1991;Zhao,2000;Grote et al., 2001;Li et al.,2005;Fukui et al.,2004;Zhou and Hao,2008).Strain rates reported to be of relevance in rock mechanics range from 10À14to108sÀ1(Olsson,1991)and indicate the large variations in strains to which rocks can be subjected depending on the mode of fracturing or excavation.
The dynamic mechanical properties of rock are basic informa-tion in assessing the stability of rock structures under dynamic loads.Mechanical properties of rocks,including compressive strength,tensile strength,shear strength and toughness,are af-fected by the rate of loading(Chong et al.,1980;Blanton,1981; Zhao et al.,1999a;Cho et al.,2003;Fukui et al.,2004;Li et al., 2011).To analyse and design underground structures in rock,it is essential to obtain dynamic characteristics and to understand the rock dynamic properties under different strain rates.A consid-erable number of studies have been conducted in recent years to study the strain rate effect on rock strength.Studies have indic
ated that rock material strengths generally increase with increasing loading rate(Zhao et al.,1998,1999a,1999b;Zhao,2000;Li et al.,2005;Cai et al.,2007;Kubota et al.,2008).
Researches on the strength of rocks at strain rates higher than 100per second have shown that the dynamic compressive strength can be expressed as(Birkimer,1971;Grady and Kipp, 1979);冶金专业
1=3
Table1
Explosives properties(Vanbrabant et al.(2002))and corresponding borehole
pressures.
ANFO explosive Emulsion explosive
V d,m/s q e,kg/m3P b,GPa V d,m/s q e,kg/m3P b,GPa
4052780 1.655821250 4.9
Fig.1.Pressure–time histories used for the numerical simulation.
268O.Yilmaz,T.Unlu/Tunnelling and Underground Space Technology38(2013)266–278
r
cd
¼65 e1=3ð13ÞIt can be said that transition point from low to high strain rate sensitivity of strength of rock takes place between10and100sÀ1 strain rates.At the stage of low strain rates,dynamic compressive strength of rocks may increase up to20%.This probably takes place at a point of30sÀ1strain rate.After this stage,dynamic compres-sive strength may be described by the power of1/3which is also suggested by other researchers(Birkimer,1971;Grady and Kipp, 1979;Li,1994;Li et al.,2005).In Eq.(11),coefficient K is estimated as40%of the unconfined compressive strength of rock.As a result, the dynamic increasing factor of rock material may be calculated as:
r
cd
¼0:4r cð eÞ1=3for e>30sÀ1ð14Þand the variation of dynamic rock strength according to strain rates is given in Fig.3.It can be seen from the Fig.3,Eq.(14)gives com-patible results compared to suggestions
of the other researchers.
2.3.Dynamic mechanical properties of rock mass
In this study,two different rock types have been chosen in terms of their strength.This rock types are granite for very high strength rock and sandstone for the relatively lower strength rock. Experimental studies conducted on different rock materials have shown that dynamic deformation modulus is generally greater than static deformation modulus(Bieniawski,1984;Eissa and Kazi, 1988;Price et al.,1994).Eissa and Kazi(1988)have proposed the following equation for the relationship between static and dy-namic deformation modulus of rock material.
E dyn¼ðE s=100:02Þ1=0:77=q ið15Þwhere E s(GPa),static deformation modulus,E dyn(GPa),dynamic deformation modulus and q i(g/cm3)density of the rock material. Table3gives the summary of the typical static mechanical intact rock properties used for the simulation as well as references of these values and the calculated dynamic mechanical properties of these rock materials.As the strain rates of rock material near the charge hole may range from103to105sÀ1,the strain rate value of104sÀ1has been selected as an average value in Eq.(14).
Table4gives the dynamic mechanical properties of rock masses for two different rock types used in nu
merical analyses.In this study,Geological Strength Index(GSI)values have been chosen as the same value of85for the two different type of rock material to eliminate the differences on the rock mass structure and composition status of rock mass and surface condi-tions of discontinuities).The GSI value of85represents intact rock specimens or massive in situ rock with few widely spaced discon-tinuities in very good or good surface conditions.We selected this value to consider the rock mass in good conditions.This selection does not specify any real conditions.Since,FLAC3D(v.2.1)is used Mohr-Coulomb failure criteria like most of geo-technical code,rock mass properties such as internal friction angle and cohesion,in this case,must be calculated.In Table4,the dynamic rock mass defor-mation modulus(E rm)dyn,the dynamic rock mass tensile strength (r tm)dyn,the dynamic cohesion of rock mass(c m)dyn and the dy-namic internal friction angle of rock mass(u m)dyn were estimated from the dynamic unconfined compressive strength of the intact rock material(r ci)dyn,the dynamic intact deformation modulus
(E i)dyn and the GSI by using Eq.(16)(Hoek and Diederichs,2006),
(17)(Hoek and Brown,1988),(18)and(19)(Hoek et al.,2002).It is important to note that the failure envelope range(r3max)has a direct effect on the calculated Mohr-Coulomb parameters.In this study,r3max is selected as equal to r ci/4.This is based on the empirical observation that the stress ran
ge associated with brittle failure,occurs when r3max is less than about one-quarter of r ci (Hoek et al.,2002).This equations has originally been given for sta-tic conditions by the authors,but in this study these equations have been adapted for the dynamic conditions by replacing the sta-tic values with the dynamic ones.
E rm¼E i0:02þ
1ÀD=2
1þeðð60þ15DÀGSIÞ=11Þ
ð16Þ
Table2
Fig.2.Pressure–time history applied to the borehole wall in detonation velocity.
O.Yilmaz,T.Unlu/Tunnelling and Underground Space Technology38(2013)266–278269
r tm ¼a r ci
¼m b Àffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffim 2b þ4s
q ð17Þ
c ¼
r ci ½ð1þ2a Þs þð1þa Þm b r 3n ðs þm b r 3n Þa À1
ð1þa Þð2þa Þffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1þ6am b ðs þm b r 3n Þa À1=ðð1þa Þð2þa ÞÞ
q ð18Þ
u ¼sin À1
6am b ðs þm b r 3n Þ
a À1
2ð1þa Þð2þa Þþ6am b ðs þm b r 3n Þa À1
"
#
ð19Þ
m b ¼m i exp
GSI À100
28À14D
;s ¼exp GSI À100
9À ;
a ¼12þ
16
e ÀGSI 15e À203
ð20a Þr 3n ¼r 3max =r ci
ð20b Þ
where D is the factor to allow for the effects of blast damage and stress relaxation (06D 61)which chosen as 0in this study.2.4.Boundary conditions and mechanical damping
In dynamic problems of numerical models,boundary conditions cause the reflection of outward propagating waves back into the model and do not allow the necessary energy radiation.The use of a larger model can minimise the problem,since material damp-ing will absorb most of the energy reflected from outer boundaries.However,this solution leads to a large computational burden.To avoid this,viscous (or absorbing)boundaries can be used as an alternative.The viscous boundary developed by Lysmer and Kuhle-meyer (1969)was used in the study.
In time-domain programs,Rayleigh damping is commonly used to provide damping that is approximately frequency-independent over a restricted range of frequencies.Although Rayleigh damping embodies two viscous elements (in which the absorbed energy is dependent on frequency),the frequency-dependent effects are ar-ranged to cancel out at the frequencies of interest.In rock and soil,the natural damping is independent of frequency (Gemant and Jackson,1937).
Alternatively,the local damping embodied in FLAC 3D ’s static solution scheme may be used dynamically,but with a damping coefficient appropriate to wave propagation.The use of local damping is simpler than Rayleigh damping,because there is no need to specify a frequency.Thus local damping in dynamic problems is useful as an approximate way to include hysteretic damping (Itasca,2003).
It is known that damping has much less importance in control-ling the maximum response of a structure to impulsive loads than for periodic and harmonic loads (Clough and Penzien,1975).The maximum response to an impulsive load will be reached in very short time,before the damping forces can absorb much energy from the structure.However,if a linear failure criterion
like
Variation of dynamic increasing factor with strain rate for various materials.(See above-mentioned references for further Table 3
Selected static intact rock properties and calculated dynamic intact rock properties.Intact rock material properties
Granite Bukit
Timah (Singapore)Zhao (1996)Sandstone Zonguldak (Turkey)Gerçek and
Müftüog ˘lu (1993)Unconfined compressive strength (r ci )MPa
18698Deformation modulus (E i )GPa 84
31Poisson ratio (t )0.250.23Density (q )kg/m 3
26102638Brazilian tensile strength (r tB )MPa
À11.1À7.9Hoek-Brown strength arameter,m i
15.311.1Dynamic unconfined compressive strength (r ci )dyn (MPa)(Eq.(15))
1600
840
Dynamic deformation modulus (E i )dyn (GPa)(Eq.(16))
11431
Table 4
Dynamic mechanical properties of rock masses with GSI value of 85(r 3max =r ci /4).Dynamic mechanical properties of rock mass Granite Bukit Timah (Singapore)Sandstone
Zonguldak (Turkey)Rock mass deformation modulus (E rm )dyn GPa 10629Rock mass tensile strength (r tm )dyn (MPa)
À33.8À24.5Cohesion (c m )dyn (MPa)
17390Internal friction angle (u m )dyn
(°)
44
41
发布评论