The free surface influence on blast vibration

International Journal of Rock Mechanics & Mining Sciences 77 (2015) 182–191

Contents lists available at ScienceDirect

International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms

The free surface influence on blast vibration D.P. Blair n TNL Consultants, 20 Cinnabar Place, Carine, WA 6020, Australia

art ic l e i nf o

a b s t r a c t

Article history: Received 25 September 2014 Received in revised form 23 March 2015 Accepted 19 April 2015 Available online 16 May 2015

When a charged blasthole initiates it typically radiates P waves and shear-vertical (SV) waves. Furthermore, dependent upon the velocity of detonation (VoD), P-Mach and/or SV-Mach waves will form, which will contain a significant proportion of the blast vibrational energy. It is important to understand the influence of the free surface on such wave radiation when considering near field vibration monitoring for material damage/fragmentation and far field surface vibration monitoring for statutory compliance. The reflection characteristics of plane P- and SV-waves incident on a free surface are given by the well-known Zoeppritz equations. However, it is much more difficult to derive the reflection characteristics for the non-planar waves radiated from a blasthole. This work describes an analytical model in which the non-planar incident radiation from a blasthole is given by the Heelan solution and its reflection from the free surface is approximated using the Zoeppritz equations. Predictions of the free surface influence using this analytical model are in good agreement with the predictions from a dynamic finite element model. The analytical model also shows that as the VoD increases, the reflected vibrational energy is channelled further away from the near field. Thus the model predicts that the explosive VoD can be altered to control damage in specific near-field regions. In this regard, the present analytical method produces a complete near field vibration damage model in which the VoD and the free surface play crucial roles. At high VoD, the radiation approaches grazing incidence at the surface in the far field. This phenomenon is relevant to monitoring far field vibration. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Blast vibration Free surface Analytical/numerical models Near field vibration damage model

1. Introduction In many mine-blasting applications, vibration monitoring usually plays an important experimental role in any attempt to determine the rock response to explosive loads. There are two broad distance zones of interest for monitoring blast-induced vibration—the near field where rock damage can result from excessive vibration, and the far field where the levels of vibration are often subject to statutory compliance. Although the vibration monitor is typically located on the surface for compliance monitoring, it is often located beneath the surface for near field monitoring. Thus it is important to investigate the influence of the free surface for monitors placed at surface and sub-surface locations. Any surface, whether free, fixed or partly constrained, will disturb the wave field of incident vibration. For example, body waves at non-vertical incidence on a free boundary will generate reflected and mode-converted body waves as well as surface waves. Although the exact analytical solutions for plane wave reflection and mode-conversion at a free boundary are given by the

n

Tel.: þ 61 458132139. E-mail address: [email protected]

http://dx.doi.org/10.1016/j.ijrmms.2015.04.006 1365-1609/& 2015 Elsevier Ltd. All rights reserved.

well-known Zoeppritz equations [1], it is not possible to give an exact analytical solution for the surface waves. In the most simplistic form of blast wave reflection in the near field, Persson et al. [2] assume that the influence of a nearby surface is to double all motion close to that surface. Unfortunately, this assumption is only valid for vertically incident waves, and is in serious error for waves at oblique incidence and for all reflected waves at depth beneath the surface. Jiang et al. [3] also considered blast wave reflection in the near field, and gave a realistic treatment for the reflection of planar and non-planar waves at all angles of incidence. However their treatment is limited solely to surface motion under the influence of an incident Pwave only. In this regard, analytical models of wave radiation from a blasthole show that a significant amount of energy is transported in the shear-vertical (SV) mode [4,5], and so it is required to examine the reflection of S-waves from the boundary. These analytical models also predict that the radiation from an element of charge (used to construct the full column of explosive) is non-planar and thus the resulting radiation from the full column is also non-planar. With regard to the reflection of non-planar waves, it is possible to express a spherical wave as a superposition of plane waves using the Weyl integral, then apply the Zoeppritz equations to each plane wave component to ultimately determine the reflected spherical wave [1]. However, in the present situation there are two main problems with this approach.

D.P. Blair / International Journal of Rock Mechanics & Mining Sciences 77 (2015) 182–191

First such a calculation would require a significant computing overhead, and secondly (and more importantly), the radiation from a small element of charge is non-planar and also non-spherical, and thus it is not clear how to derive the exact reflection characteristics of these particular non-planar P- and SV-waves. It is for such reasons that the plane wave approximation is used. Use of this approximation is further justified by the fact that Jiang et al. [3] have shown for a finite length of charge there is not a significant difference in the reflection characteristics between plane and curved P-wavefronts, and this difference also decreases significantly, even for moderate source depths. Rossmanith and Uenishi [6] used fracture mechanics and analytical models of wave propagation to investigate blast wave reflection in scaled laboratory samples of cubes and cylinders. The models and experiments demonstrated the significant influence of wave reflection and mode conversion due to the surrounding boundaries, which is entirely consistent with the approaches given in [1]. However, unlike the more realistic solutions presented in [3,5], the solution used in [6] assumed an infinite length of charge. The main aim of the present work is to develop an approximate analytical model in which non-planar incident radiation of P- and SV-waves is given by the Heelan solution [4] and the reflection of these waves at a free surface is calculated using the plane wave Zoeppritz equations. Surface waves will be neglected due to their difficult analytical nature. The first part of this work describes the ray path treatment for the surface reflection of incident waves produced by a single element of charge given by the Heelan solution [7,4]. The second part of this work describes the incident and reflected radiation produced by a full column of explosive modelled as a series of Heelan elements, appropriately delayed according to the explosive velocity of detonation. The final part of this work is to compare the predictions of the approximate analytical model with the numerical results obtained from a dynamic finite element model (DFEM) of the radiation from a full column of explosive.

2. Wave reflection for a single element of charge 2.1. Ray path treatment for wave reflection The vibration due to an element of explosive charge at S(rS, zS) and detected at a monitoring point M(rM, zM) beneath the surface is due to direct waves (P, S), reflected waves (PP, SS) and mode converted waves (PS, SP). The ray paths for these six wave fields are illustrated in Fig. 1, in which the source lies at a depth b-zS and the monitor lies at a depth b-zM. It appears that a solution involving all these waves to form a total vibration wave at such a monitor has not been previously reported in the literature. This total waveform obviously depends on the angles shown in Fig. 1, which must be evaluated. The following treatment makes use of the ray parameter, p, as defined by Aki and Richards [1] in which i is the angle between a P-wave and the vertical and j is the angle between an S-wave and the vertical. All angles are measured positive anti-clockwise. If

Fig. 1. Ray paths for direct, reflected and mode-converted waves; Poisson solid.

183

rM 4 rS (as in Fig. 1) then according to the sign convention, angles

γ1, i1, ϕ and j2 are inherently positive whereas the remainder are inherently negative. However, if the source lies to the right of the monitor (rM orS) then all angles change sign. If α and β are the Pand S-wave velocities, respectively, then a generalised ray parameter, p, is defined as p¼

j sin ij

α

¼

j sin jj

ð1Þ

β

Thus the angles for the mode-converted waves are defined by








sin i1
¼ αp ;
sin j
¼ βp ;
sin i2
¼ αp ;
sin j
¼ β p 1 1 1 2 2 2 ð2Þ From Fig. 1 the angle quantities, as follows: tan ϕ ¼

ϕ can be obtained directly from known

rM  rS 2b  zS  zM

ð3Þ

The following function f(ϕ) can then be defined in terms of known quantities: f ðϕÞ ¼ ð1 þ tan 2 ϕÞfðb  zM Þ2 þ ðb  zS Þ2 g þ 2ðb  zM Þðb  zS Þ tan 2 ϕ ð4Þ Using the geometry of Fig. 1, it can be shown that the solutions for p1 and p2 are given by the following four equations: f ðϕÞ  gðp1 Þ ¼ 0

ð5Þ

f ðϕÞ  hðp2 Þ ¼ 0

ð6Þ

gðp1 Þ ¼

ðb  zM Þ2 ðb zS Þ2 2αβ p21 ðb zM Þðb zS Þ þ þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1  α2 p21 1  β 2 p2 ð1  α2 p2 Þð1  β p2 Þ

ð7Þ

ðb  zS Þ2 2αβp22 ðb  zM Þðb  zS Þ þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1  α2 p22 ð1  α2 p2 Þð1  β p2 Þ

ð8Þ

1

hðp2 Þ ¼

ðb  zM Þ2 1  β p22 2

þ

1

2

1

2

Brent’s method was used to solve Eqs. (5) and (6) for p1 and p2, which must also lie in the following ranges:







sin ϕ
o αp o
sin γ
;
sin ϕ
o αp o
sin γ j ð9Þ 1 2 1 2 Using Eq. (2), the relevant angles shown in Fig. 1 can be pffiffiffi calculated; in this particular case α =β ¼ 3 was assumed (i.e. a Poisson solid). 2.2. The Heelan element of charge As illustrated in Fig. 1, the direct components (P, S) travelling from the source to the monitor are obviously independent of the surface and can be calculated using any analytical solution in which the P- and Swaves can be separately identified (i.e. decoupled), such as the Heelan solution, for example. Although Blair [4] has shown that this solution overestimates the true wave motion at the high frequencies in the near field of a blast, it can be used to provide a relative assessment of the wave motion, except for regions very close to the blasthole axis. The exact solution given by Blair [5], which does not suffer from the Heelan limitations, could also be used. However, the P- and S-waves are not always decoupled in this solution, especially in regions close to the blasthole, and this natural phenomenon makes such a solution difficult to implement. Furthermore, this exact solution is significantly more numerically intensive than the Heelan solution. It is for these reasons that the Heelan solution is used in the present analysis. Assuming an infinite medium, Fig. 2 shows the geometry applicable to the Heelan solution [7] for a single element of charge located at S(rS, zS), for which the P- and S-wave particle velocity components at the monitor are given by the orthogonal components vR(R, ψ) and vψ(R, ψ), respectively; arrows indicate the

184

D.P. Blair / International Journal of Rock Mechanics & Mining Sciences 77 (2015) 182–191

SS) and mode-converted (PS, SP) wave fields can be determined by the use of three fictitious sources placed above the horizontal surface as illustrated in Fig. 3, for which SAA ¼SA, SBB ¼SB and SCC ¼SC. The r and z components of these fictitious sources are given by their appropriate Heelan solutions. The total solution at the monitor station, M, requires an additional three general Heelan solutions given by vr(RB þBM, ψB) and vz(RB þBM, ψB); vr(RA þAM, ψA) and vz(RA þ AM, ψA); vr(RC þBM, ψC) and vz(RC þ CM, ψC). These Heelan solutions can be obtained from the known quantities  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rM  rS RM ¼ ðr M  r S Þ2 þ ðzM  zS Þ2 ; ψ M ¼ π  atan 2b  zM  zS ð11aÞ

Fig. 2. Geometry for the Heelan solution in an infinite medium.

positive direction for all relevant variables. These Heelan solutions are decoupled, and their (r, z) components for P- and S-waves radiated from a blasthole of radius a can be derived for the present geometry from the modified form given in [4]:  2   β α f P ðζ Þ ¼ 1  2 cos 2 ζ ; f S ðζ Þ ¼ ð10aÞ sin 2ζ

α

β

RB þ BM ¼

ð2b  zS  zM Þ ; cos ϕ

RA þ AM ¼

b zS b  zM þ ; cos j1 cos i1

ψB ¼ ϕ

ð11bÞ

ψ A ¼ i1

ð11cÞ

alternatively, RA þ AM ¼

ðb  zS Þð1  β =αÞ ðr M  r S Þðβ=αÞ  cos j1 sin j1

ð11dÞ

b zS b  zM þ ; cos i2 cos j2

ð11eÞ

vrP ¼ gðR; t  t P Þf P ðψ Þ sin ψ

ð10bÞ

vzP ¼  gðR; t  t P Þf P ðψ Þ cos ψ

ð10cÞ

RC þ CM ¼

vrS ¼ gðR; t  t S Þf S ðψ Þ cos ψ

ð10dÞ

vzS ¼ gðR; t  t S Þf S ðψ Þ sin ψ

ð10eÞ

In the case of near-vertical incidence (i1, j1-0), Eq. (11c) is required, and in the case when i1-π/2 (i.e. j1-arcsin[β/α]), Eq. (11d) is required; outside these limiting cases, either equation may be used. The travel times from source to monitor for the direct and reflected waves are

gðR; tÞ ¼ khðtÞ

a  R

ð10f Þ

The functions fP(ζ) and fS(ζ) are the generalised Heelan radiation patterns for P-wave and S-wave radiation, respectively, at any angle ζ. The function h(t) is a suitably normalised waveform shape as a function of time, t, and is dependent on the pressure-time history of the explosive (see Appendix A); tP and tS are the arrival times of the P- and S-wave, respectively, k is a constant dependent on the material shear modulus as well as the peak (Von Neumann) pressure produced by the explosive. 2.3. Wave reflection from a Heelan element of charge The method of images is the traditional analytical approach used to account for the surface reflection due to any buried source. In this method, the particle motions due to the real source are superimposed on the motions due to its image source symmetrically located above the surface. The calculation of a normal stress field,  σzz, to be applied to the boundary, is also required in order to balance the stress field, σzz, that exists in the absence of any free surface. This image method, which does not involve the direct calculation of wave reflection coefficients, was used by Lee [8] to account for the surface reflection of a dilatational point source. However, lengthy and time-consuming calculations are required to calculate the required stress field, even for such a simple source. It would clearly be prohibitive to attempt this approach for the more complicated Heelan source. Thus a solution is sought using a modified method of images that requires three fictitious sources and calculation of wave reflection coefficients, but does not require the lengthy calculations for a balancing stress field. This modified method is now discussed. Under the modified method of images, the amplitude of the horizontal (r) and vertical (z) components for the Heelan radiation in the presence of a free surface can be given by considering the components of the six wave fields (P, S, PP, PS, SP, SS) reaching the monitor. The exact directions and amplitudes of the reflected (PP,

tP ¼

RM

α

;

tS ¼

RM

β

;

t PP ¼

ψ C ¼ j2

2b zS  zM ; α cos ϕ

t SS ¼

2b  zS zM β cos ϕ

ð12aÞ

The travel times from source to monitor for mode converted waves are t PS ¼

b zS

þ

b  zM

α cos i2 β cos j2

;

t SP ¼

b  zS

þ

b zM

β cos j1 α cos i1

ð12bÞ

alternatively, t SP ¼

ðb  zS Þf1  ðβ=αÞ2 g ðr M  r S Þðβ =αÞ2  β cos j1 β sin j1

ð12cÞ

Fig. 3. Modified method of images used to evaluate the total wave field at M.

D.P. Blair / International Journal of Rock Mechanics & Mining Sciences 77 (2015) 182–191

185

If the wave reflection coefficients for the PP, SS, PS and SP waves are ΩPP, ΩSS, ΩPS and ΩSP, respectively, then the (r, z) components of reflected waves arriving at the monitor are given by vrPP ¼ gðRB þ BM ; t  t PP ÞΩPP f P ðϕÞ sin ϕ

ð13aÞ

vzPP ¼  gðRB þ BM ; t  t PP ÞΩPP f P ðϕÞ cos ϕ

ð13bÞ

vrSS ¼ gðRB þ BM ; t  t SS ÞΩSS f S ðϕÞ cos ϕ

ð13cÞ

vzSS ¼ gðRB þ BM ; t  t SS ÞΩSS f S ðϕÞ sin ϕ

ð13dÞ

The mode-converted waves arriving at the monitor are given by vrPS ¼ gðRC þ CM ; t  t PS ÞΩPS f P ðj2 Þ cos j2

ð14aÞ

vzPS ¼ gðRC þ CM ; t  t PS ÞΩPS f P ðj2 Þ sin j2

ð14bÞ

vrSP ¼ gðRC þ CM ; t  t SP ÞΩSP f S ði1 Þ sin i1

ð14cÞ

vzSP ¼ gðRC þ CM ; t  t SP ÞΩSP f S ði1 Þ cos i1

ð14dÞ

Fig. 4. Critical behaviour of SP rays arriving at 20 near-surface stations.

The reflection coefficients are given in Appendix A, and it should be noted that all these coefficients are real with the exception of ΩSS, which may be complex. This completes the approximate solution for a single element of charge. A valuable insight to the approximations involved in the present model can be obtained by considering a monitor at various depths, where S-waves at critical incidence play a role. A useful check on some of the equations can also be obtained by considering P-waves incident on a monitor placed precisely at the surface. 2.4. Monitors at various depths For the ray model illustrated in Fig. 1, only the six rays, P, S, PP, SS, PS and SP, are considered for each monitoring point. However, ray theory demands that other distinct waves will also be generated at each surface point. These waves can be described by lower case nomenclature. For example, the incident S-wave at A also generates a reflected S-wave at A that may be termed ss. Likewise there are ps and sp waves generated at B due to the incident P- and S-waves, respectively, and a reflected P-wave, pp, is also generated at the point C. The basic assumption of the present model is that these 4 wave types, pp, ss, ps and sp, are directed away from the point M and so do not contribute to the total vibration at this monitoring station. In other words, this assumption implies that M must be sufficiently below the surface. If the monitor is located at a great depth then all the angles j1, i1, ϕ, i2 and j2, are small and there is no S ray at critical incidence at the surface. As the monitor approaches the surface along a fixed vertical line all these angles increase and ϕ is the first angle to go critical. This occurs when ϕ ¼arcsin(β/α) and then the incident Swave at the surface point B produces the critically refracted Pwave (sp) that travels along the surface. Fig. 1 was constructed assuming β/α ¼1/√3, yielding a critical angle of 35.3 degrees. The angle ϕ as shown is 53.1 degrees, which is much larger than critical, yet despite this, all the 6 rays remain distinct and the ray sp does not influence the vibration at M. However, when ϕ does exceed the critical angle, as in Fig.1, the reflection coefficient, ΩSS, becomes complex (see Appendix A) but nevertheless can still be modelled for ϕ approaching π/2 (i.e. the source approaching the surface). As the monitor approaches the surface, ϕ approaches arctan[(rM rS)/(b  zS)] and j1 approaches the critical angle. In order to illustrate this critical incidence, Fig. 4 shows a line of 20 monitors all at a very small fixed depth, d. The mode-converted rays, SP and PS, arriving at each station from the source, S, are also shown. The critical behaviour of the SP rays, whose surface intersection point is r ¼rC, defines a critical incident angle j1 ¼arcsin(β/α), for which i1-π/2. Thus, as illustrated by Fig. 4,

the basic assumption of the present model is that the angle j1 can never exceed arcsin(β/α). 2.5. P-waves incident on a surface monitor If the monitor is located beneath the surface, there is a direct Pwave ray path to the monitor as well as 2 distinct P-wave ray paths incident at the surface points B and C (Fig. 1). However, if the monitor is on the surface then B and C merge to a single point, in other words there is no differentiation between the pp and PP waves and no differentiation between the ps and PS waves. The ray paths, SB, SC and SM, from the source to the monitor are then identical. Under this condition,

   β i2 ¼  ϕ; j2 ¼ arcsin sin ϕ ð15aÞ

α

ϕ ¼ arctan RB þ BM ¼



rM  rS b  zS



¼ π ψ

b  zS ¼ RC þ CM ¼ R cos ϕ

t PP ¼ t PS ¼ t P

ð15bÞ

ð15cÞ ð15dÞ

Eq. (15) can then be used to calculate the r and z components of the P-wave amplification factors, χPr and χPz, for surface monitors. These amplification factors are defined for a given ξ-direction (r or z) as the ratio of the total surface motion to the incident motion. Thus,

χ ξP ¼

ðvξP þ vξPP þvξPS Þ v ξP

ð16Þ

Using Eqs. (10) and (13)–(16) it can be shown that

χ rP ¼ 1 þ ΩPP þ ΩPS

f P ðj2 Þ cos j2 f P ðϕÞ sin ϕ

ð17aÞ

χ zP ¼ 1  ΩPP þ ΩPS

f P ðj2 Þ sin j2 f P ðϕÞ cos ϕ

ð17bÞ

If plane wave (rather than Heelan) radiation is assumed, i.e. fP(ζ) ¼1 for all ζ, then Eq. (17) reduce to those given by Jiang et al. [3] for the amplitude of reflected plane P-waves. The vector amplification factor, χVP, for P-waves is given by

χ 2VP ¼ χ 2rP sin 2 ϕ þ χ 2zP cos 2 ϕ

ð18Þ

Fig. 5 shows plots of this vector amplification factor as a function of the surface incident angle, ϕ. If Poisson’s ratio, ν, is 0.5, then fP(j2)/fP(ϕ)¼ 1, and so the χVP plot for ν ¼ 0.5 is identical to that for the plane wave results shown in [3].

186

D.P. Blair / International Journal of Rock Mechanics & Mining Sciences 77 (2015) 182–191

If ϕ o ϕC :

If ϕ Z ϕC :

t SP -t S

ð19fÞ

2rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3  2 β 1  6 α tan ϕ7 7 t SP -ðb zS Þ6 þ 4 α 5 β

The S-wave amplification factors, χSr and arbitrarily close to the surface are given by

χ ξS ¼

ðvξS þ vξSS þvξSP Þ v ξS

ð19gÞ

χSz, for a monitor ð20Þ

The S-wave amplifications are somewhat more difficult to evaluate than the P-wave amplifications because of the possibility of reflection beyond the critical angle. However, if ϕ (and hence j1) is less than ϕC, then Eqs. (10), (13), (14) and (20) yield

Fig. 5. The P-wave vector amplification versus incident angle.

χ rS ¼ 1 þ ΩSS þ ΩSP

f S ði1 Þ sin i1 f S ðϕÞ cos ϕ

ð21aÞ

χ zS ¼ 1  ΩSS þ ΩSP

f S ði1 Þ cos i1 f S ðϕÞ sin ϕ

ð21bÞ

The vector amplification factor,

χ

Fig. 6. Magnitude of S-wave vector amplification versus incident angle.

2 VS

¼ χ cos ϕ þ χ 2 rS

2

2 zS

χVS, for S-waves is given by

sin ϕ 2

ð22Þ

However, when ϕ exceeds ϕC, the reflection coefficient ΩSS becomes complex, producing a change in shape of the elemental waveform, h(t), of Eq. (10f); this aspect is detailed in Appendix A. The real (ΩRSS) and imaginary (ΩISS) components of ΩSS are defined by Eq. (A13). The amplification factors χrS, χzS and χVS in Eq. (22) are then all complex, and the magnitude of the S-wave vector amplification is given by the magnitude of the complex square root


qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi





χ
¼
χ 2 cos 2 ϕ þ χ 2 sin 2 ϕ
ð23Þ rS zS
VS


Fig. 6 shows the magnitude of this vector amplification factor as a function of the surface incident angle, ϕ. The sharp features in most plots occur when the S-wave incident angle equals the critical angle, which, in terms of Poisson’s ratio, ν, is defined by sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi β 1  2ν sin ϕC ¼ ¼ ð24Þ 2ð1  νÞ α

Fig. 7. Ray paths from a 5 m column of explosive arriving at a monitor, M.

2.6. S-waves incident on a near-surface monitor Here it is necessary to distinguish between S-waves below and above critical incidence. The critical angle, ϕC, is defined by sinϕC ¼ (β/α), and if the monitor depth, d, approaches zero, then the following conditions hold

   α If ϕ o ϕC : i1 -arcsin sin ϕ ; j1 -  ϕ ð19aÞ

β

π

If ϕ Z ϕC :

i1 - ; 2

If ϕ o ϕC :

RA þ AM-

3. An extended column of explosives

j1 -  ϕC

t SS ¼ t S

ð19bÞ

b  zS b  zS ¼R cos j1 cos ϕ 2

If ϕ Z ϕC :

If ν-0.5, then (β/α)-0 and all incident angles are critical according to Eq. (24). If ν ¼0.0, then (β/α)¼ 0.707 and ϕ is critical for all angles above 45 degrees. Figs. 5 and 6 show that the surface motion is doubled for vertically incident waves, and this is the basis of the simplistic form of blast wave reflection applied by Persson et al. [2] for near-field vibration. However, Figs. 5 and 6 also show that for blast waves at non-vertical incidence (which would certainly represent the vast majority of near-field incident waves) the surface motion is much more complex than that assumed in [2].

ð19cÞ 3

6 1  β =α 7 RA þ AM-ðb zS Þ4qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ tan ϕ5 2 1  ðβ =αÞ

ð19dÞ

ð19eÞ

The ray path solution is now applied to each of the N Heelan elements assumed to constitute a column of explosives of length L; the radiation from each element is represented by a single ray. The element length is fixed at 0.1a, thus for an 89 mm diameter blasthole there are N ¼1123 elements along a 10 m charge column. The waveforms previously calculated for a single, nth, element are appropriately delayed according to the velocity of detonation, VD, and then superimposed to form a total vibration due to the entire column. Fig. 7 illustrates the ray paths for twenty-one elements within a 5 m column of explosive, showing the twenty-one beams

D.P. Blair / International Journal of Rock Mechanics & Mining Sciences 77 (2015) 182–191

187

given by V ξW ðtÞ ¼

N X n¼1

vnξW ¼

N X

g RnW ; t  t W ðn 1ÞL=V D hnξW ðθn Þ

n¼1

ð27Þ Hence the vector sum wave at M due to all wave types (P, S, PP, SS, PS, SP) and all elements of the charge is given by the following sum over all W wave types vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 12 0 12 u0 u X X u ð28Þ V rW ðtÞA þ @ V zW ðtÞA VðtÞ ¼ t@ W

Fig. 8. Snapshot of wave fronts; head wave 15 m from base of charge.

Fig. 9. Snapshot of wave fronts; head wave 23 m from base of charge.

due to each wave type arriving at the monitor M. Of course, interference will occur between beams of differing wave types for regions where they overlap in space and time; such interference is neglected in the present model except where the beams meet at M. In fact, waveform calculations are only made at the point M. Eqs. (10a–10f), (13a–13d) and (14a–14d) can be used to describe the vibration from the nth Heelan element of charge, which can be cast in its most general form as vnξW ¼ gðRnW ; t  t nW ÞhnξW ðθn Þ

ð25Þ

where ξ denotes the component (r or z), W denotes the wave type (P, S, PP etc), RnW and tnW denote the distance and travel time, respectively, for the nth ray from its element to monitor for wave type W, and hnξW(θn) is a function dependent upon the n th ray and an associated angle θn. For example, assuming an SP wave, this function is given by analogy to Eq. (14c) as hnξW ðθn Þ ¼ ΩnSP f S ðin1 Þ sin in1

ð26Þ

where ΩnSP is the SP-wave reflection coefficient for the ray from the nth element, and in1 is the ray’s refraction angle. It should also be appreciated that the general form of the function hnξW(θn) still applies for the direct (P, S) waves. However in this case, unlike Eq. (26), a reflection coefficient is not involved and θn is only an incident angle. Due to the finite velocity of detonation, the time delay of the nth element is given by (n  1)L/VD. Thus, according to Eq. (25) the superposition of all n elemental charges produces a ξ-component (r or z) of vibration at M due to wave type W that is

W

If there are k monitoring stations specified over a two dimensional grid, then the vector waveform of Eq. (28) can be evaluated for all times, t, and all stations, k, and stored in a two dimensional array V(k,t). This array can then be used to determine the peak value of the vector wave over all times for each monitor station. The peak value is then stored in the one dimensional array VMAX(k), which can be used to contour the vector peak particle velocity (VPPV) over the grid region if required. Furthermore, the array V(k,t) can also be used to determine the value of the vector wave itself at a specified time, tS, for all monitors. This wave value is stored in the one dimensional array V(k,tS), which can be used to display a snapshot of the waveforms over the entire grid region at the particular time tS. In the present analytical model, these snapshots were obtained from vector waveforms stored at each point in a 301  301 grid covering the 15  15 m region of interest. The standard model under consideration in this study is a 10 m column of explosive having a diameter of 89 mm and a velocity of detonation, VD, of 5500 m/s. As noted previously, the column consists of 1123 elements and the base of the charge is 15 m below the free surface. The rock properties assumed are α ¼4000 m/s and ν ¼0.25, hence (α/β)2 ¼ 3 and the Mach number (VD/α) is 1.375. Because VD 4 α there will be a P-Mach (PM) wave as well as an S-Mach (SM) wave. Fig. 8 shows a snapshot of the wavefronts from the explosive column, taken at the time instant, tS, when the head wave has just reached the surface, i.e. has travelled 15 m from the initiation point (base of charge). The head wave is defined as the fastest moving wave of the system, and since VD 4 α its travel time is given by the time taken for the velocity of detonation to travel the column length plus the time taken for the P-wave to travel the extra 5 m, thus tS ¼10/VD þ5/ α ¼3.07 ms. At this time both conical Mach waves have decoupled from the blasthole (i.e. have degenerated), yet still contain most of the vibrational energy. Spherically spreading P- and S-waves radiate off the base (loading phase, L) and top (unloading phase, U) of the explosive column. Thus, as shown in Fig. 8, there are 6 individual wavefronts radiating from the explosive column, and each of these wavefronts has the possibility of generating a reflected and mode converted wave at the surface. In other words, an explosive column near a free surface can generate up to 18 particular wavefronts. Fig. 9 shows a snapshot of the wavefronts when the head wave has travelled 23 m from the base of the charge and is obviously outside the 15  15 m grid. The degenerate P-Mach wave can still be seen, along with its reflection (PMP) and mode converted (PMS) waves as well as the degenerate S-Mach wave along with its reflection (SMS) and mode converted (SMP) waves. The elemental shape for the SMS wave is phase-dependent, as given by Eq. (A15) and illustrated in Fig. A1, and thus is different to all other elemental shapes. Multiple snapshots are now used to show the directivity of wavefronts in time. In this regard, Fig. 10a shows the results for 130 snapshots in a large-scale model (60  15 m); the separation in head wave travel distance between each snapshot is 1 m, and the Mach number remains at 1.375. Fig. 10b shows

188

D.P. Blair / International Journal of Rock Mechanics & Mining Sciences 77 (2015) 182–191

of Blair [9], who did not model reflection. The present work shows that this energy at near-vertical incidence is efficiently reflected directly downwards as a dominant PP wave. It is obvious from Fig. 10a and b that the VoD (i.e. Mach number) controls the directivity of blast waves reflected from the surface, and as the VoD increases, the reflected vibrational energy is channelled further away from the near field. Thus the model can be used to predict near-field vibration damage in which the explosive VoD can be altered to control damage in specific near-field regions. At high VoD, the radiation approaches grazing incidence at the surface in the far field. This phenomenon is relevant to monitoring far field vibration.

4. A numerical model of the free surface influence

Fig. 10. (a) Progression of wavefronts; Mach number 1.375. (b) Progression of wavefronts; Mach number 1.00.

Fig. 11. Axisymmetric DFEM with dashpots attached to an outer 45  45 m grid.

Fig.12. DFEM Snapshot of wave fronts; head wave 23 m from base of charge.

the results for a slightly smaller model (45  15 m) having a Mach number of 1.0 (i.e. there is no P-Mach wave). In this case, the S-Mach wave is incident at the surface in the nearer field and also reflects from the surface at a steeper angle. Furthermore, for this Mach number, there is significant energy radiated from the blasthole directly above the source, which is consistent with the results

An axisymmetric dynamic finite element model (DFEM) was set up to compare its numerical predictions with the analytical results for the previous 15  15 m model. It is well known that the size and non-uniformity of the finite elements within a given dynamic numerical model introduce spurious effects that are not present in the continuum from which the discrete model is derived [10,11]. In particular, any DFEM prediction is dependent on the relative size of the finite elements compared to the shortest wavelength propagated throughout the grid. In this regard, small, uniform element sizes are the optimum arrangement for high frequency propagation. In conducting a study of such propagation throughout an axisymmetric DFEM grid, Blair [12] found that between six and twelve elements per wavelength were required in order to avoid any spurious effects. In this regard it has been shown [5] that applied pressure loads from realistic explosives have onset rise times less than 10 μs. This is consistent with the induced vibration wave shape shown in Fig. A1, which has a rise time of 8 μs. An estimation of the shortest wavelength of interest is given by the lowest velocity of propagation (VS) multiplied by the rise time of the pressure load. This gives a minimum wavelength of 23 mm. Using this information, three tests were run for a DFEM grid using square elements of side 20 mm, (barely 1 element per wavelength), 10 mm and 5 mm (approximately 5 elements per wavelength). It was too time-consuming to use smaller elements. As expected, the models with 10 mm and 20 mm elements produced a significant level of spurious noise, and thus the square element size was fixed at 5 mm. Fig. 11 shows the DFEM central 15  15 m grid having nine million square elements. In order to reduce unwanted reflections from the finite boundaries, the model size was extended to a 45  45 m grid using additional graded elements to the outer boundaries, with energy absorbing dashpots attached to these boundaries. The complete model consisted of 11.56 million elements, and, like the analytical solution, the contouring of wavefronts was done over 301  301 grid points in the 15  15 m zone. Fig. 12 shows the DFEM predicted snapshot of the wavefronts when the head wave has travelled 23 m from the base of the charge. A comparison of Figs. 12 and 9 shows that the analytical and DFEM results are in good agreement with regard to all wavefronts interacting with the free horizontal surface. The significant variation of the SMS component across its wavefront is due to the influence of S-wave incidence beyond the critical angle, and this variation is clearly seen in the analytical and numerical solutions. However, it should be noted that the DFEM wavefronts (Fig. 12) are somewhat broader than the analytical wavefronts (Fig. 9). This is consistent with the results of Blair [4] who showed that realistic wave shapes (as assessed by either DFEM or full-field solutions) are broader than the Heelan wave shapes if the mean dimensionless frequency of the applied pressure load is sufficiently high, as in the present case. Ironically, this limitation of the

D.P. Blair / International Journal of Rock Mechanics & Mining Sciences 77 (2015) 182–191

Fig. 13. (a) Blasthole radiation reflected by a presplit; Mach number 1.286. (b) Blasthole radiation reflected by a presplit; Mach number 1.833.

Heelan wave shapes improves the resolution of the wavefronts, which, in turn, provides an improved insight to the behaviour of blast waves interacting with a free surface. It should also be noted that the Heelan solution cannot predict the interaction of the blast waves with the blasthole itself, whereas the DFEM naturally includes this realistic phenomenon. This is clearly seen in Fig. 12 as disturbances in the wave field for regions directly above the explosive column.

5. Discussion As noted previously, surface waves have been neglected in the present analytical model. However, the free-surface DFEM would naturally include such waves. The fact that these waves are not apparent in the DFEM results of Fig. 12 is evidence that they are not significant within the 15  15 m zone considered. This zone would constitute the near field, and thus the present analysis suggests that neglecting surface waves in the analytical model is reasonable when assessing blast vibration in the near field of an explosive column. Thus the model can be used to predict nearfield vibration damage in which the explosive VoD can be varied to produce damage in specific regions throughout the near field as indicated by Fig. 10a and b. The present analytical model is one of axisymmetric radiation from a blasthole interacting with a single free planar surface. However, although only vertical blastholes have been considered

189

so far, there is no restriction on the orientation of the surface and/ or the blasthole. The solution only requires an axis of symmetry along the blasthole, and all the blast vibration components can be rotated through any desired angle to mimic wave radiation from a blasthole of any orientation and reflecting off a planar surface whose normal has any orientation. The one restriction here is that the solution is only valid for the plane containing the blasthole axis and the normal to the surface (i.e. the plane of the paper). It is not possible to construct this particular axisymmetric model using a numerical dynamic code such as DFEM. These codes require a single axis of symmetry for the entire model, which includes any structures, boundaries and free surfaces. For example, if the blasthole is precisely vertical, then a horizontal free surface can be specified in any axisymmetric DFEM as illustrated in Fig. 11b. However, if the blasthole is not vertical, then it is impossible to have a horizontal planar surface; to model this arrangement would require a DFEM that is fully three-dimensional. This is a severe disadvantage compared to the present analytical approach, because such a numerical model would require at least 15 m length in the third dimension to avoid wave reflections from those extra boundaries. Furthermore, this 15 m span, too, would require at least 3000 elements of 5 mm size to avoid spurious contamination of the wave propagation (as previously noted). Thus, at least 27 billion three-dimensional elements would be needed to achieve the required spatial resolution of uncontaminated wavefronts. It is for this reason that current three-dimensional numerical models of near-field blasting could be questioned, simply because they either use unrealistically low frequency propagation and/or use unrealistically large elements. Another distinct advantage of the analytical model is its flexibility for isolating individual wave types. For example, the detailed behaviour of the SS-wave is easily achieved by suppressing all other wave types in the analytical solution; this was invaluable in determining the nature of the sharp cut-off associated with the SMS wavefront as shown in Figs. 9 and 12. Suppressing individual wave types is not possible using DFEM, which solves the entire problem with all waves naturally included. As examples of the capability of the present analytical model, Fig. 13a and b show the radiation from an angled buffer blasthole, of diameter 114 mm, close to a presplit, which is relevant to wall control blasting. The presplit is considered to act as a free surface (i.e. an open joint) and thus it will not transmit any radiation; there is no horizontal surface in this model. In both cases a Poisson solid (ν ¼0.25) is assumed; in Fig. 13a α ¼3500 m/s, VD ¼ 4500 m/s hence the Mach number (VD/α) is 1.286 and in Fig. 13b α ¼3000 ms/s, VD ¼5500 m/s hence the Mach number is 1.833. All the relevant wave types are identified. The relatively low Mach number describing the radiation in Fig. 13a produces a relatively large angle of incidence for the PM wave at the presplit, which results in an insignificant PMP wave (not labelled). Thus for this particular incident angle almost all the P-Mach energy is converted to S-waves, in other words the reflection coefficient, ΩPP, approaches zero. The higher Mach number describing the radiation in Fig. 13b produces a relatively small angle of incidence, which, in turn, results in a reasonable amount of energy in the PMP wave.

6. Conclusions An analytical model has been developed for blast wave radiation interacting with a free surface. There are two main limitations of the model. First, the model assumes plane wave reflection coefficients for the surface rather than the actual situation of curved wave reflection coefficients. Second, the model cannot account for the complex surface waves. In this regard, any reasonable DFEM solution would naturally include the complex influence of curved wavefront reflection from the surface as well as the

190

D.P. Blair / International Journal of Rock Mechanics & Mining Sciences 77 (2015) 182–191

formation and propagation of surface waves. Despite these limitations of the analytical model, the evidence via a comparison with an appropriate DFEM solution shows that such limitations are not significant, at least within 15 m (the near field) of the blast. Thus the present model can be considered as a reasonable engineering solution that yields insights on the distribution of blast damage. In this regard, the results of Fig. 10a and b imply that the explosive velocity of detonation has a significant influence on the distribution of blast damage around a blasthole close to a free surface. In the far field, however, the predictions from the analytical model could be problematic simply because this is the region where the surface waves are expected to dominate. Nevertheless, the evidence of Fig. 10a and b still suggests that increasing the VoD will direct more energy towards the surface. It is reasonable to assume that such a process will also increase the surface wave contribution to the far field. Although the DFEM solution has distinct advantages over the analytical model, it suffers from the problems of mesh dependency and reflections from the model boundary. Clearly the analytical model has no such limitations. Furthermore, any numerical approach, such as DFEM, faces severe practical limitations when attempting to solve fully three-dimensional problems, even a problem as simple as an angled blasthole beneath a horizontal surface. The fact that this specific problem can be easily achieved with the analytical model is a clear indication of its power and flexibility. Blair [9] has also raised the limitations with current three-dimensional numerical models of blast wave propagation and, furthermore, suggested that problems might still arise with the resolution of two-dimensional numerical models. The resolution of the analytical model is much easier and simpler to control than that of any mesh-dependent numerical model because it does not suffer from the problem of spurious oscillations.

Appendix A. The reflection coefficients

Ω3 ΩPS ¼ ; Ω1 þ Ω2

Ω  Ω2 Ω4 ΩSS ¼ 1 ;Ω ¼ Ω1 þ Ω2 SP Ω1 þ Ω2

ðA1Þ

Ω1 ðpÞ ¼

1

β2

 2p2

Ω3 ðp; iÞ ¼ 4p

1

β2

Ω2 ðp; i; jÞ ¼ 4p2

; !

 2p2

cos i ;

β

cos i cos j

α

Ω4 ðp; jÞ ¼ 4p

β

!

1

β2

 2p2

cos j

However, for the present analysis, these reflection coefficients must be given in terms of the incident angles j1, i2 and ϕ, illustrated in Fig. 1. Thus the various wave types are characterised by sin ϕ ;

i ¼  ϕ;

α

sin ϕ ; SS wave : p ¼ PS wave : p ¼ SP wave : p ¼

j sin i2 j

α

j sin j1 j

β

j ¼ arcsin½ðβ = αÞ sin ϕ

i ¼ arcsin½ðα=β Þ sin ϕ;

β

ΩSS ¼

ΩSP ¼

ðA5Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

cos 2 2ϕ  2ðβ=αÞ sin 2ϕ sin

ϕ 1  ðα=βÞ2 sin 2 ϕ

cos 2 2ϕ þ 2ðβ=αÞ sin 2ϕ sin

ϕ 1  ðα=βÞ2 sin 2 ϕ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ðA6Þ

2j sin 2j1 cos 2j1 j qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  ðα=βÞ2 sin 2 j1

ðA7Þ

ðα=βÞ cos 2 2j1 þ 2 sin 2j1 sin j1

It should be noted in Eqs. (A4) to (A7), that for any angle ζ r π/2, terms such as sin2 ζsin ζ ( ¼2sin2 ζcos ζ) Z0. Also, for angles approaching vertical incidence, sin ζ approaches ζ for all angles of incidence and reflection, then neglecting all terms in sin2 ζ, it is simple to show that

ΩPP -  1:0;





ΩPS -4
i2
β=α;

ΩSS -1:0;





ΩSP -4
j1
β=α

ðA8Þ

Eqs. (A8) and (18a) to (21b) can then be used to determine all surface amplification factors for P- and S-waves at vertical incidence.

χ rP -4β=α; χ rS -2;

χ zP -2;

‘ χ V P -2

ðA9Þ

χ zS -4β=α;

‘ χ V S -2

ðA10Þ

Grazing incidence occurs when the source and monitor are at the same depth, and this depth approaches zero. Under these conditions:













 














ϕ
o
i1
;
ϕ
-π =2;
i2
¼
i1
-π =2;
j1
¼
j2
-arcsin β ;















j ¼ ϕ

i ¼ i2 ;

;



i ¼ i1 ¼ arcsin½ðα=βÞ
sin j1
;

Then it can be shown that

ΩPP -  1:0;

ΩPS -0;

ΩSS -1:0;

ΩSP -

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 1  ðβ=αÞ2 ðα=β Þ2  2

ðA12Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

cos 2 2ϕ  2iðβ=αÞ sin 2ϕ sin

ϕ ðα=βÞ2 sin 2 ϕ  1

cos 2 2ϕ þ 2iðβ=αÞ sin 2ϕ sin

ϕ ðα=βÞ2 sin 2 ϕ  1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ΩRSS þ iΩISS

The magnitude and phase of ΩSS are given by

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  


ΩSS
¼ Ω2 þ Ω2 ¼ 1:0; pha ¼ arctan ΩISS SS RSS ISS

Ω

ðA13Þ

ðA14Þ

RSS

The non-zero phase beyond the critical angle implies a change in shape of the elemental waveform, h(t) of Eq. (10f). If HI is the Hilbert transform operator, then the elemental waveform for the SS wave beyond critical incidence is given by [1] hSS ðtÞ ¼ cos ðphaSS ÞhðtÞ þ sin ðphaSS ÞH I ½hðtÞ

ðA15Þ

The elemental function, h(t), has the form [4]



j ¼ j2 ¼ arcsin½ðβ =αÞ
sin i2


;

ðA11Þ

ΩSS ¼

α ðA2Þ

PP wave : p ¼

ðA4Þ   2 2ðβ=αÞj sin 2i2 j 1  2ðβ=αÞ2 sin i2 ΩPS ¼  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 2ðβ=αÞ2 sin 2 i2 þ 2ðβ=αÞ3 sin 2i2 sin i2 1 ðβ=αÞ2 sin 2 i2

Another case of interest also arises when ϕ exceeds ϕC, i.e. when ϕ 4arcsin(β/α); then Eq. (A6) may be rewritten as the complex quantity

where !2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2  1  2ðβ=αÞ2 sin 2 ϕ þ 2ðβ =αÞ3 sin 2ϕ sin ϕ 1 ðβ =αÞ2 sin 2 ϕ ΩPP ¼  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1  2ðβ=αÞ2 sin 2 ϕ þ 2ðβ=αÞ3 sin 2ϕ sin ϕ 1  ðβ=αÞ2 sin 2 ϕ

α

Aki and Richards [1] give the reflection coefficients for incident planar waves as a function of the general ray parameter, p, as defined by Eq. (1):  Ω1 þ Ω2 ΩPP ¼ ; Ω1 þ Ω2

yields

hðtÞ ¼ ½t η  2ηt η  1 þ ηðη  1Þt η  2 e  γ t

j ¼ j1 ðA3Þ

The reflection coefficients for the ray paths as shown in Fig. 1 can be obtained by substituting Eqs. (A2) and (A3) into (A1). This

ðA16Þ

where γ is a decay parameter describing the explosive pressure load on the borehole wall and η is a constant; in the present model η ¼ 6 and γ ¼200,000. Fig. A1 shows a plot of the elemental waveform, hSS(t), as a function of the time t, for increments of π/8 of the phase function,

D.P. Blair / International Journal of Rock Mechanics & Mining Sciences 77 (2015) 182–191

Fig. A1. Wave shapes for hSS(t) at π/8 increments in phaSS.

phaSS. Thus the waveform for phaSS ¼ 0 represents the original shape, h(t); the shape for phaSS ¼ π/2 represents its Hilbert transform. In the present work, this Hilbert transform is evaluated upfront, using the Fast Fourier Transform, and stored in an array for continual accessing. References [1] Aki K, Richards P. Quantitative seismology. second ed. San Francisco: University Science Books; 2002. [2] Persson PA, Holmberg R, Lee J. Rock blasting and explosives engineering. Boca Raton, Fla: CRC; 1993.

191

[3] Jiang J, Baird GR, Blair DP. Polarization and amplitude attributes of reflected plane and spherical waves. Geophys J Int 1998;132:577–83. [4] Blair DP. A comparison of Heelan and exact solutions for seismic radiation from a short cylindrical charge. Geophysics 2007;72:E33–41. [5] Blair DP. Seismic radiation from an explosive column. Geophysics 2010;75: E55–E65. [6] Rossmanith HP, Uenishi K. On size and boundary effects in scaled model blasts spatial problems. Int J Frag Blast 2005;9:139–74. [7] Heelan PA. Radiation from a cylindrical source of finite length. Geophysics 1953;18:685–96. [8] Lee I. Transient ground motion in an elastic homogeneous halfspace to blasting loading. Soil Dynam Earthq Eng 1996;15:151–9. [9] Blair DP. Blast vibration dependence on charge length, velocity of detonation and layered media. Int J Rock Mech Min Sci 2014;65:29–39. [10] Celep Z, Bazant ZP. Spurious reflection of elastic waves due to gradually changing finite element size. Int J Numer Method 1983;19:631–46. [11] Holmes N, Belytschko T. Post processing of finite element transient response calculations by digital filters. Comp Struct 1976;6:211–6. [12] Blair DP. Acoustic pulse transmission in half-spaces and finite-length cylindrical rods. J Acoust Soc Am 1985;50:1676–83.