Modelling of blast-induced damage in tunnels using a hybrid fi nite-discrete numerical approach

Amichai Mitelman,Davide Elmo

NBK Institute of Mining Engineering,University of British Columbia,Vancouver,Canada

Modelling of blast-induced damage in tunnels using a hybrid fi nite-discrete numerical approach

Amichai Mitelman,Davide Elmo*

NBK Institute of Mining Engineering,University of British Columbia,Vancouver,Canada

A R T I C L E I N F O

Article history:

Received 21 August 2014

Received in revised form

23 September 2014

Accepted 25 September 2014

Available online 29 October 2014

Blast damage

This paper presents the application of a hybrid f i nite-discrete element method to study blast-induced damage in circular tunnels.An extensive database of f i eld tests of underground explosions above tunnels is used for calibrating and validating the proposed numerical method;the numerical results are shown to be in good agreement with published data for large-scale physical experiments.The method is then used to investigate the in f l uence of rock strength properties on tunnel durability to withstand blast loads.The presented analysis considers blast damage in tunnels excavated through relatively weak (sandstone)and strong(granite)rock materials.It was found that higher rock strength will increase the tunnel resistance to the load on one hand,but decrease attenuation on the other hand.Thus,under certain conditions,results for weak and strong rock masses are similar.

?2014 Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.Production and hosting by Elsevier B.V.All rights reserved.

1.Introduction

Underground excavations in rock(e.g.civil road tunnels,mine drifts)can be subjected to blast loads associated with excavation methods,or may be even subjected to blast loads caused by external sources.Blast damage includes both cracking of the rock mass material and induced damage to structural reinforcing elements(concrete liner).Even in the case where the blast may not necessarily reduce the load capacity of the engineered excavation, there is the potential for fragments of rock material and/or structural elements to be ejected with large velocities,thus imposing a signi f i cant hazard to either humans or equipment.

Due to the complexity of the mechanisms involved,blast design in construction projects and mining largely relies on simpli f i ed empirical approaches.Most commonly,peak particle velocity(PPV) attenuation is estimated based on f i eld tests,and compared to PPV based damage thresholds(Dowding,1996).Different authors used numerical simulations to study the response of underground structures to blasting(e.g.Jiang and Zhou,2012;Deng et al.,2014), and rock mass damage is de f i ned based on the observation of the plastic zones created by the blast and/or by measuring PPV results (Wei and Zhao,2008).

In this paper,a fracture mechanics based f i nite-discrete element approach(FEM-DEM)is adopted,using the proprietary code ELFEN (Rock f i eld,2007).In the models,blast-induced cracking and spalling of the rock material are simulated using a Rankine rotating crack failure criterion.The hybrid FEM-DEM approach allows for an immediate and explicit simulation of the damage caused to the tunnel walls.Blast load generated by the explosive detonation is initially estimated using the ANSYS Autodyn software(ANSYS, 2013)and subsequently inserted into ELFEN.

For calibrating and validating the proposed numerical method, observations and results from extensive f i eld tests conducted by Engineering Research Associates are used,hereafter referred to as the ERA tests(ERA,1953).In these tests,single delay charges in the range of 145-145,000 kg of TNT were detonated above unlined tunnels with diameters of 2-10 m in sandstone and granite.All charges were buried and fully coupled to the ground.Four damage zones from total collapse to light damage were empirically de f i ned as a function of the scaled distance of the charge to the tunnel.

The in f l uence of rock mass strength on tunnel durability to withstand dynamic loads is debated amongst authors.For instance, Rozen et al.(1988)proposed an empirical correction factor to the empirical guidelines established in the ERA tests based on the Rock Quality Designation(RQD)(Deere and Miller,1966)of the rock mass.They found that for lower RQDs the PPV induced by blasting increases.In contrast,Wu et al.(1998)emphasized the role of discontinuities in the rock mass on wave attenuation,implying that a weaker and heavily jointed rock mass is favourable in terms of tunnel dynamic resistance.Once the proposed method of numerical simulation is found to be consistent with the ERA tests,the method is then extended to attempt to determine the overall impact of rock mass strength on tunnel dynamic strength.

2.The f i nite-discrete element method

As discussed in Hamdi et al.(2014),two main approaches are used for the numerical modelling of rock mass behavior,based on the concept that the deformation of a rock mass subjected to applied external loads can be considered to be either continuous or discontinuous.The main differences between the continuum and discontinuum analysis techniques lie in the conceptualization and modelling of the fractured rock mass and the subsequent deformation that can take place in it.A continuum model re f l ects mainly material deformation of the system,while a discontinuum model re f l ects the movement component of the system.The continuum approach may circumvent some of the dif f i culties associated with the discrete method,in terms of complexity of the model and impracticality of modelling every fracture in a deterministic way. However,an intrinsic limitation of the equivalent continuum approach is that the stress acting on a speci f i c fracture is usually not the same as that deduced from the overall stress,because it depends on the stiffness of the fracture itself and on the stiffness of the fracture’s surrounding matrix(Cai and Horii,1993).

Hybrid fi nite-discrete element(FEM-DEM)codes combine the aspects of both fi nite elements and discrete elements,and also allow for the incorporation of fracture-mechanics principles to allow for the realistic simulation of brittle fracture-driven processes and a full consideration of the failure kinematics(Pine et al.,2006; Mahabadi et al.,2012;Hamdi et al.,2014).In FEM-DEM model,the fi nite element-based analysis of continua is merged with discrete element-based transient dynamics,contact detection,and contact interaction solutions(Munjiza,2004).FEM-DEM based numerical analysis of fracturing processes in rock considers that such problems are often highly dynamic,with rapidly changing domain con f i gurations,thus requiring suf f i cient resolution and allowing for multiphysics phenomena.Such problems are typically simulated employing time-integration schemes of an explicit nature(Owen et al.,2004).Application of dynamic explicit time-integration schemes to multifracturing solids,particularly to those involving high nonlinearity and complex contact conditions,has increased notably in recent years(e.g.Owen et al.,2004;Jaini and Feng,2011).

There are advantages in employing a hybrid FEM-DEM approach to model blast-induced damage,including:

(1)A better description of the physical processes involved,accounting for diverse geometrical shapes and effective handling of large numbers of contact entities with speci f i c interaction laws.

(2)The implementation of speci f i c fracture criteria and propagation mechanisms allows the simulation of the progressive fracture process within both the f i nite and discrete elements.

Among the different hybrid FEM-DEM codes currently available, the code ELFEN(Rock f i eld,2007)incorporates a coupled,elastoplastic,fracture-mechanics constitutive criterion that allows realistic modelling of the transition from a continuum to a discontinuum,with the explicit generation of stress-induced cracks.

As an FEM/DEM code,ELFEN has the capability of modelling preexisting discontinuities.In the current paper,the rock mass is modeled as an equivalent continuum.The effect of joints on wave propagation has been investigated by different authors(Cai and Zhao,2000;Chen et al.,2000).Work is being carried out to test the proposed approach with the addition of discontinuities preinserted in the model.

Within the ELFEN code,the constitutive behavior used to simulate multi-fracturing of brittle materials is achieved by employing a fracture energy approach controlled by designated constitutive fracture criteria.In this paper,the rotating crack model is used to simulate crack formation under tensile conditions within the initially continuum-meshed geometry.

The Rankine rotating crack failure criterion is based on the concept of Mode I fracturing studied in fracture mechanics.Once the maximum principal stress reaches the tensile strength limit, tensile softening is initiated and the elastic modulus is degraded in the direction of the major principal stress invariant.Finally,the mesh topology is updated and when new surfaces and/or bodies are formed they interact with each other according to the discrete contact properties assigned(Rock f i eld,2007).The yield surface and softening curve for the Rankine rotating crack failure criterion are shown in Fig.1.

3.Assessment and characterization of blast load

To the authors’knowledge,there is no available software program that can model all stages of blast-induced damage(i.e. explosive detonation,wave propagation,fracturing and spalling). Therefore,it was decided to simulate the load generated by the explosive detonation using ANSYS Autodyn(ANSYS,2013),similar to the work carried out by Chen and Zhao(1998).Autodyn is a fi nite-difference software,specially designed to solve a wide variety of non-linear problems including the resultant stresses emanated from explosive materials using the empirical Jones-Wilkins-Lee(JWL)equations.

Fig.1.Yield surface and softening curve for the Rankine rotating crack in ELFEN(from ELFEN user’s manual(Rock f i eld,2007)),whereσ1andσ2are the tensile strengths,ftis the elastic limit,σnnand εnnare the principal stress and strain,and E and Edare the elastic and residual Young’s moduli.

Three material models are used for this stage:TNT,air,and the rock material.The TNT material properties are based on the inhouse Autodyn material library and presented in Table 1.The input parameters for the air material are also obtained from the Autodyn material library.The equation of state(EOS)for air is the ideal gas equation.The internal energy corresponding to the atmospheric pressure is assigned to the air material as an initial condition.

The mean values of the rock material properties from different sites of the ERA tests are listed in Table 2.The sandstone properties from the ERA tests refer to intact rock material.It is interesting to note that no rating of the rock mass with respect to the presence of discontinuities is available from the ERA tests since the tests were conducted before the development of the rock mass classi f i cation systems most commonly used today,i.e.RMR(Bieniawski,1989),Q-index(Barton et al.,1974),GSI(Hoek et al.,1995).Considering that rock mass conditions for the ERA tunnels were such to allow for unsupported tunnels to be excavated and that the presence of joints results in a reduction in rock mass strength,the massive rock mass conditions were assumed for the test conditions(e.g.GSI and/or RMR greater than 80).The GSI index has the advantage of being related to the Hoek-Brown failure criterion for rock masses,which is widely accepted in geotechnical and rock engineering applications.In the current analysis,the Hoek-Brown parameters and the fi tting Mohr-Coulomb parameters are initially estimated (assuming GSI of 80 for the rock mass)and then converted to the generalized Drucker-Prager strength model that is available in Autodyn.The Hoek-Brown and Mohr-Coulomb parameters for the sandstone material used for the Autodyn analysis are listed in Table 3.

The EOS for a material de f i nes the relationship between the hydrostatic pressure,density,and internal energy.This relationship is affected by the strain rate of the load applied,implying that the bulk modulus increases with increasing strain rate(Chen et al., 2000;Zhao,2000).As there is no information regarding the dynamic properties of the rock from the ERA tests,the bulk and shear static moduli from Table 3 are inserted into the Autodyn model and gradually increased in order that PPV results are consistent with the ERA tests.An increase by a factor of 25%is found to be appropriate for the EOS of the sandstone material.

The model geometry is shown in Fig.2.In the axisymmetric model,the TNT material is located in proximity of the model origin (bottom-left corner in Fig.1).History points Nos.1,5 and 9 are set at a distance of 2.4 m from the TNT location to provide the necessary stress input for the ELFEN models.The integrated Autodyn-ELFEN modelling approach is designed to take advantage of the strengths of Autodyn for blast modelling.The large deformations that the supersonic region undergoes due to the extreme pressures generated by blast are accommodated by the Autodyn coupled Euler-Lagrange solver.In the supersonic region,the shock wave conforms to the Rankine-Hugoniot relations which are calculated by Autodyn.Therefore the stress input for the ELFEN model is de f i ned using the history points at some distance away from the blast location rather than using the stress history of the TNT-rock material boundary.The procedure is schematically illustrated in Fig.3.The remaining history points are used to con f i rm that the simulated PPV attenuation results are in agreement with the f i eld tests.

Table 1TNT material properties.

Table 2Sandstone averaged intact properties from the ERA tests.

Table 3Mohr-Coulomb f i t to Hoek-Brown failure criterion for the sandstone material.

Kuhlmeyer and Lysmer(1973)found that in order to obtain realistic results in dynamic simulations,the mesh size must be no larger than one-eighth of the minimal wavelength.This criterion is used for selecting the mesh size for both the Autodyn and ELFEN models.The hypothesis is also validated by performing a mesh sensitivity analysis and results converge at a mesh size of 12 cm.

Extensive tests of ground shock effects in soils carried out by Drake and Little(1983)found that the pressure-time histories consist of a compressive pulse with a short rise time followed by a negligible negative tensile pulse.Dominant frequency has been found to be an additional important parameter for assessment of blasting impact(Dowding,1996).In the current analysis,a triangular shaped pulse with the frequency and peak compressive stress determined by the Autodyn model is selected for the subsequent ELFEN modelling.The load function used in the ELFEN model is shown in Fig.3,based on the Autodyn results.The averaged peak stress of the compressive pulse and pulse duration from the Autodyn results are 320 MPa and 0.8 ms,respectively.

4.Modelling of damage zones based on ERA test results

4.1.ELFEN modelling methodology and set-up

Fig.2.Model geometry and history points layout for the Autodyn model.

As introduced in Section 2,ELFEN is a hybrid FEM-DEM code for 2D and 3D modelling and is extensively applied to modelling ofrock mechanics problems(e.g.Cai and Kaiser,2004;Stead et al., 2004;Elmo et al.,2008).Note that in the literature,generally the coupled Rankine-Mohr-Coulomb failure criterion is used in ELFEN to simulate quasi-static rock behavior.However,the Rankine rotating crack failure criterion is suf f i cient as spalling failure caused by tensile re f l ection from the tunnel roof is governed by the tensile strength alone.

The parameters required for this criterion are the tensile strength of the rock and the fracture energy.The tensile strength of the sandstone material used is taken from the intact strength measured in the ERA tests.The fracture energy GFis the integral of the stress over the strain in the softening portion of the stressstrain curve,and is related to the critical stress intensity factor KIC(or fracture toughness)and elastic modulus E(Rock f i eld,2007):

Fig.3.Procedure for blast load function estimation from Autodyn to ELFEN models.

The fracture toughness KICcan be estimated using an empirical relationship with the rock tensile strength(Zhang,2002):

Therefore it is possible to estimate the fracture energy using Eqs. (1)and(2)solely based on the knowledge of the rock elastic modulus and tensile strength.The Rankine rotating crack input parameters used in the ELFEN models for the sandstone material are listed in Table 4.

Four damage zones are identi f i ed in the ERA test results based on the scaled distance from the exploding charge(the scaled distance is de f i ned as the distance from the exploding charge divided by the cube root of the charge weight).The four zones are categorized as follows:

Table 4The input parameters for Rankine rotating crack failure criterion.

(1)Zone 1:total collapse of the tunnel;

(2)Zone 2:heavy damage;

(3)Zone 3:moderate damage;

(4)Zone 4:light damage.

Each zone is characterized by the damage area,de f i ned as the difference in cross-section before and after the blast,and the maximum velocity of the broken rock fragments(spalls)ejected into the tunnel.Tunnels in distances greater than Zone 4 exhibited no damage at all.

Waves that propagate through a medium attenuate due to geometrical spreading and material damping.PPV attenuation is commonly simpli f i ed to the form of:

where R is the radial distance from the blast,W is the weight of the explosive,b in the case of a spherical charge refers to the cube root, and K and n are site-speci f i c constants.There is no analytical method of prediction of these constants based on the materialproperties;therefore for the numerical simulations arti f i cial global damping must be applied in order to obtain the attenuation rates measured in f i eld.As shown in Fig.4,two damping regions are de f i ned in the ELFENmodels.Region 1 is the area around the loaded arc.According to the modelling results,when cracking is simulated, the wave attenuates rapidly.Therefore a lower damping percentage is required within Region 1 to yield a wave attenuation that agrees with the physical measurements from the ERA tests.

The following procedure is implemented in the ELFEN models:

(1)An initial ELFEN model is created that corresponds to a speci f i c distance of the damage zones in the ERA tests.

(2)The model above is calibrated so that the damping percentage for Regions 1 and 2 in the model yields PPV results that are consistent with the PPV results measured in the ERA tests.

(3)The other models(Zones 2 and 3,and model of no damage)are then assigned the same damping percentage as the calibrated model.

(4)Results of damage area and velocity of the spalled rock in all models are then compared to those of the ERA tests.

The distance of the center of the charge to the blast for the different models is listed in Table 5.The model set-up for the Zone 4 model used for calibration is shown in Fig.4.Note that Zone 1 damage,de f i ned as the distance in which complete collapse of the tunnel occurs,is not considered in the current analysis,since the complex phenomena that occur in the supersonic zone may not be properly captured by the current modelling procedure.An additional model representing a zone with no damage is included as well in the analysis.

A tunnel subjected to a spherically propagating blast wave is essentially a 3D problem.However,it is assumed that 2D plane strain modelling is acceptable.The 2D simpli f i cation is probably somewhat conservative as it implies that the maximum blast load acts in f i nitely.Works done by others(Chen et al.,2000;Deng et al., 2014)have shown that 2D models can be successful in predicting blasting outcomes tested in f i eld.

The tunnel diameter is 10 m.The model boundaries are set to a distance suf f i ciently large to avoid wave re f l ections.Deng et al. (2014)found the initial stresses to have little in f l uence in depths of up to 200 m.The tunnels in the ERA tests were located at shallower depths;therefore the effects of initially induced stresses associated with the tunnel excavation are not included in the models.

Fig.4.Model geometry for the Zone 4 model.

Table 5Vertical distance of blast charge to tunnel roof for models of the ERA damage zones.

4.2.Simulated results of damage zones

A qualitative description of explosive interaction with rock was given by Brady and Brown(2004),who stated that subsequently to blasting,the region surrounding the blasthole would experience expansion and dense fracturing,while farther away from the blast, radial cracks would appear.Finally,if the wave is re f l ected from a free face,additional fractures in the form of spalling may appear. Fig.5 shows the fracture pattern generated in the calibrated ELFEN model for Zone 3,which matches the qualitative description given by Brady and Brown(2004).

Simulated results of the damaged area above the tunnel roof and the velocity of the spalled rock fragments fromthe three models are presented in Table 6 along with the corresponding results measured in the ERA tests.

For the Zone 4 model,fracturing occurs above the tunnel boundary and does not cause the rock to detach and f l y into the tunnel space.This is consistent with the ERA observations,as in Zone 4,pieces of rocks were found to fall into the tunnel with no initial velocity,and the detachment of these pieces is attributed topre-existing joints.No fracturing occurred in the model at a distance greater than Zone 4.

In the ERA tests it was observed that,at closer distances to the blast,rocks are ejected radially,and,as the distance increases,the spalling tends to occur only above the tunnel roof.In addition,it was found in the ERA tests that spalling begins with small particle size with high velocity ejections from the tunnel perimeter.The subsequent deeper ejections involve comparatively larger particle size with low ejection velocity.Both these phenomena are captured in the results of the ELFEN models and can be observed in Figs.6 and 7.

Note that the models are simulated for a duration of 10 ms which is suf f i cient for capturing the complete kinematics of the failed blocks.In order to display spalling in a vivid manner such as the image displayed in Fig.6,relatively long computing times are required so that the detached pieces gain distance from their initial location.

Fig.5.(a)Rock fracturing due to an explosive charge,modi f i ed from Brady and Brown (2004).(b)Results of Zone 3 model.

Table 6Comparison of results for damage zones from the ERA tests and ELFEN models.

Fig.6.Progressed spalling in Zone 3 model.

5.Effects of varying rock materials and properties on tunnel durability

5.1.Modelling set-up

When considering the in f l uence of rock strength properties on the problem of tunnel subjected to blast loads,it is important to consider the interrelated effect of strength properties on both tensile failure(cracking)and wave attenuation.According to Zhou (2011),for massive rocks the attenuation will be less rapid than in weak and fractured rocks.Typical values of the attenuation coef f i cient n from Eq.(3)are 1.5 for hard and strong rocks and greater than 2 for weak and soft rocks.

In an elastic analysis,the wave will attenuate only due to geometrical spreading,independent of the material properties. When the Rankine fracturing failure criterion is used,different rock types(with different tensile strengths and fracture energy)will respond differently to the load and the resultant attenuation will vary due to energy dissipation through fracturing.

Wave attenuation is considered to be site-speci f i c,and even for tests undertaken at the same site,a large scatter in results is often encountered(Dowding,1996).Due to these limitations,it is assumed that true wave attenuation rates cannot be predicted by means of numerical simulation alone.An attempt is herein made todetermine whether the simulated results in ELFEN are consistent with typical attenuation values encountered in f i eld and published in the literature.Subsequently,a preliminary comparison between rock types with different strength properties is made.

Models with material properties typical of limestone,basalt and granite are considered in the analysis,as listed in Table 7.These types of rocks were chosen as their strength spans from weak to strong,and the objective of the modelling is to capture the overall in f l uence of rock strength on the durability of the tunnel subjected to blast loading.

All models are set-up in an identical con f i guration to the models calibrated to the ERA tests discussed in the previous section,and with a distance of charge to tunnel of 14 m corresponding to“Zone 3”type damage.

Fig.7.Results for models of the damage zones:(a)Zone 2;(b)Zone 3;(c)Zone 4;and(d)Zone with no damage.The distance of the tunnel roof to the center of the charge is denoted as r.

5.2.Simulated results for different rock types

Observing the simulated results of fracturing presented in Fig.8, it can be seen that there is only minor decrease of the fractured area above the tunnel as the rock material is stronger.On the other hand it can be noted that for the stronger rock types,less fracturing occurs in the area of the blast.

These results are in agreement with the ERA tests for sandstone and granite,where the extentof damage measured in the tunnels in granite compared to those in sandstone was found to be only slightly more favourable for the granite.The discrepancy between the difference in tensile strength of the granite and sandstone to the similar amount of damage can therefore be attributed to the difference in attenuation:for the stronger rocks attenuation is less rapid,and the arriving pulse is larger.Fig.9 shows the best f i t curves for the peak stress as a function of the scaled distance from the charge for recorded from the different models.

The two contrasting effects of the rock strength,i.e.wave attenuation and rock resistance to fracturing even out each other.It is apparent that tunnel durability subjected to blast loads cannot be estimated in a straightforward manner and both effects should be weighed carefully.

PPV is recorded for varying distances for each of the different rock types and the attenuation constants n and K are found by curve fi tting.The results are presented in Table 8 and show that,as the rock material is weaker,the attenuation is more rapid.The numerical explanation for this is that the larger extent of fracturing that occurs around the area of the blast for the weaker rock dissipates more energy.

Results show that the models using the Rankine fracturing failure criterion in ELFEN yield attenuation rates for that follows the trend encountered in fi eld measurements.As is stated,this is not to say that numerical modelling is capable of predicting actualattenuation rates,as these are highly site-speci f i c.However,it is believed that the ELFEN models capture the general trend of the different rock types and can be used for gaining an initial estimation and understanding of the tunnel response to blasting in different rock types.In-situ measurements can later be used to modify the damping percentage,so the attenuation rate is better fi tted to the actual fi eld behavior.

Table 7Properties used in the models investigating the effects of varying rock types.

Fig.8.Fracturing results for a 10 m diameter tunnel in different rock types:(a)limestone;(b)sandstone;(c)basalt;and(d)granite.

6.Conclusions

An integrated approach is used in this paper to simulate the response of a tunnel subjected to blasting.The code Autodyn is used to estimate the blast load,and the hybrid FEM-DEM code ELFEN is used to simulate the fracturing of the rock resulting from the blast wave.Observations and results from f i eld tests are used to calibrate the wave attenuation and to compare the results of the damaged area and rock spall velocity.Altering the distance of the charge to the tunnel,the simulated results are found to beconsistent with the test measurements.Overall,the modelling results show that the proposed method is a reliable tool for analyzing spalling damage to tunnels induced by blasting.

Table 8PPV attenuation constants for the rock type models.

The modelling shows that a tensile failure criterion with explicit fracturing simulation is capable of capturing the rock response to blasting in terms of the fracture pattern.Varying rock material properties shows that fracturing simulation yields attenuation rates that are consistent with measurements encountered in f i eld. Furthermore,the modelling results show that due to the contrasting effects of the rock strength(i.e.tensile strength and wave attenuation)the tunnel response to a given blast load in weak and strong rocks is similar.

This method can be further extended to investigate the in f l uence of other varying conditions such as tunnel shape,in-situ stresses,tunnel support,and multiple blast loads.The proposed numerical approach has the capability of modelling pre-existing discontinuities.Although this is not considered in the current paper,the effect of joints on wave propagation is the subject of ongoing research.

Con f l ict of interest

The authors wish to con f i rm that there are no known con f l icts of interest associated with this publication and there has been no signi f i cant f i nancial support for this work that could have in f l uenced its outcome.

ANSYS.AUTODYN user manual 2014,revision 3.0.San Ramon,California:Century Dynamics;2013.

Barton NR,Lien R,Lunde J.Engineering classi f i cation of rock masses for the design of tunnel support.Rock Mechanics 1974;6(4):189-239.

Bieniawski ZT.Engineering rock mass classi f i cations.New York:Wiley;1989.

Brady BHG,Brown ET.Rock mechanics for underground mining.3rd ed.Dordrecht, Netherlands:Springer;2004.

Cai M,Kaiser PK.Numerical simulation of the Brazilian test and the tensile strength of anisotropic rocks and rocks with preexisting cracks.International Journal of Rock Mechanics and Mining Sciences 2004;41(Supp.1):478-83.

Cai JG,Zhao J.Effects of multiple parallel fractures on apparent attenuation of stress waves in rock masses.International Journal of Rock Mechanics and Mining Sciences 2000;37(4):661-82.

Cai M,Horii H.A constitutive model and FEM analysis of jointed rock masses.International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 1993;30(4):351-9.

Chen SG,Zhao J.A study of UDEC modelling for blast wave propagation in jointed rock masses.International Journal of Rock Mechanics and Mining Sciences 1998;35(1):93-9.

Chen SG,Cai JG,Zhao J,Zhou YX.Discrete element modelling of an underground explosion in a jointed rock mass.Geotechnical and Geological Engineering 2000;18(2):59-78.

Deere DU,Miller RP.Engineering classi f i cation and index properties for intact rock. Tech Rept.No.AFWL-TR-65-116.New Mexico:Air Force Weapons Lab,Kirkland Air Force Base;1966.

Deng XF,Zhu JB,Chen SG,Zhao ZY,Zhou YX,Zhao J.Numerical study on tunnel damage subject to blast-induced shock wave in jointed rock masses.Tunnelling and Underground Space Technology 2014;43:88-100.

Dowding CH.Construction vibrations.Englewood Cliffs:Prentice-Hall;1996.

Drake JL,Little CD.Ground shock from penetrating conventional weapons.In: Proceedings of the symposium of the interaction of non-nuclear munitions with structures.Colorado,USA:US Air Force Academy;1983.

Elmo D,Vyazmensky A,Stead D,Rance JR.Numerical analysis of pit wall deformations induced by block-caving mining:a combined FEM/DEM-DFN synthetic rock mass approach.In:Proceedings of the 5th international conference and exhibition on mass mining,June 2008.Lulea,Sweden:Lulea University of Technology Press;2008.p.1073-82.

Engineering Research Associates(ERA).Underground explosion test program f i nal report,vol.2.Sacramento District,USA:Sandstone Rock;1953.

Hamdi P,Stead D,Elmo D.Damage characterization during laboratory strength testing:a 3D-f i nite-discrete element approach.Computers and Geotechnics 2014;60:33-46.

Hoek ET,Kaiser PK,Bawden WF.Support of underground excavations in hard rock. Rotterdam:A.A.Balkema;1995.

Jaini ZM,Feng YT.Computational modelling of damage behaviour of reinforced concrete slabs subjected to blast loading.In:Proceedings of the 8th international conference on Structural dynamics,EURODYN 2011,Leuven,Belgium;4-6 July,2011.

Jiang N,Zhou C.Blasting vibration safety criterion for a tunnel liner structure. Tunneling and Underground Space Technology 2012;32:52-7.

Kuhlmeyer RL,Lysmer J.Finite element method accuracy for wave propagation problems.Journal of the Soil Mechanics and Foundations Division 1973;99(5): 421-7.

Mahabadi OK,Lisjak A,Munjiza A,Grasselli G.Y-Geo:new combined f i nite-discrete element numerical code for geomechanical applications.International Journal of Geomechanics 2012;12(6):676-88.

Munjiza A.The combined f i nite-discrete element method.Chichester,England,UK: John Wiley&Sons;2004.

Owen DRJ,Feng YT,de Souza Neto EA,Cottrell MG,Wang F,Andrade Pires FM.The modelling of multi-fracturing solids and particulate media.International Journal for Numerical Methods in Engineering 2004;60(1):317-39.

Pine RJ,Coggan JS,Flynn ZN,Elmo D.The development of a new numerical modelling approach for naturally fractured rock masses.Rock Mechanics and Rock Engineering 2006;39(5):395-419.

Rock f i eld.ELFEN user’s manual.Swansea,UK:Rock f i eld Software Ltd.;2007.

Rozen A,Johnson WJ,Rizzo PC.Explosions in soils:the effects of soil properties on shock attenuation.In:Proceedings of the 23rd department of defense explosives safety seminar,Atlanta,GA,USA;1988.p.1831-41.

Stead D,Coggan JS,Eberhardt E.Realistic simulation of rock slope failure mechanisms:the need to incorporate principles of fracture mechanics.International Journal of Rock Mechanics and Mining Sciences 2004;41(3):466.

Wei X,Zhao Z.Response characteristics of underground rock cavern subjected to blast load.In:Proceedings of the World Tunnel Congress 2008-underground facilities for better environment and safety,India;2008.p.665-73.

Wu YK,Hao H,Zhou YX,Chong K.Propagation characteristics of blast-induced shock waves in a jointed rock mass.Soil Dynamics and Earthquake Engineering 1998;17(6):407-12.

Zhang ZX.An empirical relation between mode I fracture toughness and the tensile strength of rock.International Journal of Rock Mechanics and Mining Sciences 2002;39(3):401-6.

Zhao J.Applicability of Mohr-Coulomb and Hoek-Brown strength criteria to the dynamic strength of brittle rock.International Journal of Rock Mechanics and Mining Sciences 2000;37(7):1115-21.

Zhou Y.Explosion loading and tunnel response.Chapter 17 in Advances in rock dynamics and applications.London,UK:CRC Press;2011.

Amichai Mitelman holds a bachelor degree in civil engineering,and he has 8 years of experience working as a civil engineer.Currently,Mr.Mitelman is completing his Master of Applied Science degree(MASc)at the NBK Institute of Mining Engineering of the University of British Columbia.

Davide Elmo holds a B.Eng(Hons)in Engineering Geology from Portsmouth University(UK)and a Ph.D.in Geomechanics from Camborne School of Mines(University of Exeter,UK).He is an assistant professor of rock mechanics at the NBK Institute of Mining Engineering of the University of British Columbia.He is a Fellow of the Geological Society of London,Fellow of the Institute of Materials, Minerals and Mining,Graduate Member of the Institute of Civil Engineers(UK)and Member of the Canadian Institute of Mining.He has over 13 years of experience in rock engineering and engineering geological studies,numerical modelling studies,mining and civil projects.He has carried out work related to numerical modelling of block cave mining,mechanical behavior of hard rock pillars,slope stability analysis,applications of synthetic rock mass modelling and discrete fracture network modelling.

*Corresponding author.Tel.:+1 604 822 9304.

E-mail address:delmo@mining.ubc.ca(D.Elmo).

Peer review under responsibility of Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.

1674-7755?2014 Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.Production and hosting by Elsevier B.V.All rights reserved.

http://dx.doi.org/10.1016/j.jrmge.2014.09.002

Tunnels

Numerical analysis

Finite-discrete method