A coupled method to study blast wave propagation in fractured rock masses and estimate unknown properties

A coupled method to study blast wave propagation in fractured rock masses and estimate unknown properties

Computers and Geotechnics 49 (2013) 134–142 Contents lists available at SciVerse ScienceDirect Computers and Geotechnics journal homepage: www.elsev...

1MB Sizes 1 Downloads 20 Views

Computers and Geotechnics 49 (2013) 134–142

Contents lists available at SciVerse ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

A coupled method to study blast wave propagation in fractured rock masses and estimate unknown properties Nima Babanouri a, Hamid Mansouri a,⇑, Saeed Karimi Nasab a, Mojtaba Bahaadini b a b

Department of Mining Engineering, Shahid Bahonar University of Kerman, Iran School of Mining Engineering, The University of New South Wales, Sydney, Australia

a r t i c l e

i n f o

Article history: Received 1 September 2012 Received in revised form 10 November 2012 Accepted 26 November 2012 Available online 3 January 2013 Keywords: Blast wave propagation Fractured rock mass Linear superposition method Three-dimensional discrete element method Simulated annealing Back analysis

a b s t r a c t The prediction of blast-induced ground vibrations across fractured rock masses is of great concern to rock engineers in evaluating the stability of rock slopes in open pit mines. Hence, many researchers have studied the wave propagation through rock slopes and tried to simulate its complexities using different methods. At the present work, the ‘linear superposition method’ was used at the Gol-e-Gohar iron ore mine to simulate the production blast seismograms based upon measurements of single-hole shot vibrations, carried out at a distance of 39 m from the blast. The production blast seismograms were then used as the input to predict particle velocity time histories in the mine wall using the three-dimensional discrete element method. By back analysis of measured blast waveforms at a distance of approximately 300 m away, mechanical properties of in-filled faults of the study area were determined. In combination with the numerical model, a ‘simulated annealing’ search algorithm was employed to find the optimum values of unknown parameters. The final results of time histories of particle velocity show a good agreement with the measured time histories. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction In general, surface rock structures are more sensitive than underground ones to dynamic loads such as earthquake and blast waves. Wave propagation through these structures and their dynamic response are very complicated especially when the rock masses are regarded as discontinuous media containing various types of fractures (in the form of micro-cracks to faults). When a blast-induced shock wave propagates in an intact medium, it is attenuated simultaneously by material damping and geometrical spreading. However, when a shock wave propagates in a jointed rock mass, wave attenuation is controlled by the characteristics of joints surfaces, the rock wave impedance of both walls of each discontinuity, and the angles between discontinuities and the direction of propagation. Hence, many researchers have studied effects of blast-induced wave propagation through rock slopes and tried to simulate the complexities using analytical, empirical, artificial intelligence, and numerical methods. Analytical solutions [1,2] encounter serious limitations when the problem geometry is complicated and/ ⇑ Corresponding author. Address: Department of Mining Engineering, Shahid Bahonar University of Kerman, Kerman 76196-37147, Iran. Tel.: +98 341 2112764; fax: +98 341 2121003. E-mail addresses: [email protected] (N. Babanouri), [email protected] ac.ir, [email protected] (H. Mansouri), [email protected] (S.K. Nasab). 0266-352X/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compgeo.2012.11.008

or a greater number of influencing parameters are to be incorporated. Furthermore, based on the results of blasts monitoring, many empirical relationships have been presented which predict the peak particle velocity (PPV) at any points considering mostly two parameters namely the maximum charge per delay and distance from the blast to the measuring point [3–9]. Geological and geotechnical parameters, blast geometry, and also characteristics of explosive charge have not yet been incorporated in this type of relationships. In reality, the data used to determine the PPV– scaled distance relations show large scatter most of the time [10,11]. Because the number of influencing parameters is too high, artificial neural networks and a small number of artificial intelligence methods were employed to predict rock blasting vibrations. Many researchers used neural network and support vector machine to estimate the PPV and air blast [12–16]. Empirical and artificial intelligence methods allow only an estimation of PPV and provide no information about the complete seismic waveform. The linear superposition method, which takes into account all the influencing parameters, models complete seismic waveforms produced by a blast. In the linear superposition method, blast vibrations are simulated at the desired location around the blast based on the vibrations of a single hole shot measured at that location. This method was used and validated by Hinzen [11]. The linear superposition method was then used to control the blasts carried out in the Val d’Azergues cement factory in France [17], and validated to simulate blasting vibrations at the Sarcheshmeh copper

N. Babanouri et al. / Computers and Geotechnics 49 (2013) 134–142

mine in Iran [18]. Toraño et al. employed the linear superposition method in analysis of ground vibrations produced by blasting at a Spanish quarry [19]. Numerical methods are able to model seismic waveforms of production blast at any distances and directions from the blast provided that the mechanical properties of media are known. Toraño et al. developed a finite element model including randomness to study the blasting vibrations [20]. In comparison with continuum-based methods such as the finite element method, the discrete element method [21] is considered more appropriate for analysis of the wave propagation through fractured rock masses. Eberhardt and Stead [22] and Stead et al. [23] used discrete element method to dynamic analysis of a natural rock slope in western Canada; the dynamic input was applied along the bottom boundary of the model as a sine wave of a specified amplitude, frequency, and duration. Bhasin and Kaynia [24] and Kveldsvik et al. [25] applied discrete element method to analysis of, respectively, 700 m and 800 m high Norwegian rock slopes and estimated the potential rock volumes associated with potential catastrophic failures. Liu et al. [26] modeled the dynamic response of a jointed rock slope to a blast by UDEC (universal distinct element code). Their results indicated that the simulated velocity history of the slope toe agreed well with in situ measurements. Fan et al. studied the effect of incident boundary on wave propagation in jointed rock masses using discrete element method [27]. For the problem of 2-D compressive wave propagation in a rock mass with parallel joints in the radial direction normal to the joints, Lei et al. performed parametric studies on the transmission ratio and the maximum rebound ratio in UDEC [28]. Zhao et al. carried out a calibration work of UDEC modeling on P-wave propagation across multiple nonlinearly deformable fractures [29]. Chen and Zhao [30,31], and Zhao and Cai [32] used UDEC to model the blast wave propagation in fractured rock masses. It was demonstrated that the rock masses containing multiple fractures have greater and faster attenuation of the blast wave than the continuous media. The above-mentioned numerical studies are all two-dimensional, while there are a smaller number of three-dimensional analyses of wave propagation within jointed rock masses. Kulatilake et al. [33] studied the effect of intermitted joints on the deformation of jointed rock mass


through 3DEC (three-dimensional distinct element code, developed by Cundall [34–36]). Lichun et al. [37] investigated the threshold control of safety blasting vibration velocity by coupling 3DEC and a search algorithm. In the regions near to the blast at large open-pit mines, where vibration intensity is too high, it is not possible to measure total blast vibrations by conventional equipment. Therefore, it is difficult to provide real dynamic inputs for numerical models. Coupling two methods, linear superposition and numerical, is likely to overcome this problem. In this article, linear superposition method and three-dimensional discrete element method were first coupled to simulate blast wave propagation through the rock masses of northern wall of the Gol-e-Gohar iron ore mine. To any unknown parameter, instead of a specified value, a range was assigned within which optimum value of each parameter was back calculated by iterative comparison of simulated blast waveforms with the real seismograms. A simulated annealing (SA) search algorithm was responsible for supervising the change of unknown parameters values in the numerical model at each run. In this research, the emphasis has been placed on estimating properties of in-filled faults of the study area and the numerical model’s additional damping. 2. Study area The Gol-e-Gohar Iron ore mine is located southeast of Iran (55 km southwest of Sirjan and 325 km northeast of Shiraz), between 55°150 4000 E and 55°220 3300 E longitudes and 29°030 1000 N and 29°070 0400 N latitudes, at an altitude of 1740 m above sea level. This area has about 1135 million tons of geological iron ore distributed within six anomalies. Among these anomalies, the mine No. 1 with 270 million tons of iron ore is under mining. Fig. 1 is a photograph of the northern sector of mine No. 1. As shown on the figure, the upper part consists of soil which is highly compacted deposit of clayey sandy gravel with cobbles, the middle part is quartz–schist, and the lower part includes magnetite as iron ore. At the present time, the northern wall has not reached its final status. Several local failures have occurred at the Gol-e-Gohar mine No. 1 which were mostly structural failure and were controlled by

Fig. 1. Photo of the study area at the Gol-e-Gohar iron ore mine.


N. Babanouri et al. / Computers and Geotechnics 49 (2013) 134–142

Fig. 2. Contours of pole concentrations along with major planes plot.

the faults. Faults of the area are mostly clay-filled with high thickness. Site investigations showed that thickness of clay filling varies within 10–50 cm. At present, comprehensive slope stability studies are in progress to redesign the Gol-e-Gohar mine walls. The aim of these studies is to prevent large scale failures which threaten the stability of the Gol-e-Gohar mine. The determination of properties of these in-filled faults is one of the most challenging issues in the slope stability studies. Three nearly mutually orthogonal major joint sets could be identified throughout quartz–schist and magnetite rock masses at the northern wall of the pit as shown in Fig. 2. During blasting operations, a staggered pattern is used and the explosive charged in dry blast holes is ANFO (ammonium nitratefuel oil). The diameter and depth of blast holes are 0.25 m and 17 m, respectively. Drilling cuttings are used as stemming materials. The number of rows and the number of holes per each row are respectively 2–6 and 10–20, and the inter-row delay is 25–80 ms. Blast holes are vertical and, the drilling operation is performed carefully with a low deviation.

3.1. Geometry and constitutive models Fig. 3 shows the 3DEC model of northern sector of the Gole-Gohar iron ore mine with a width of 300 m. As can be seen, there are three types of materials: soil, quartz–schist, and magnetite. The blast was located at magnetite region, and the monitoring spot was located in soil at a distance of approximately 300 m away from the blast. Soil medium has no joint and is only cut by the faults. At quartz–schist and magnetite regions, in addition to the faults, there are three nearly mutually orthogonal joint sets which were incorporated into the model. Fig. 4 shows a two-dimensional cross section of the 3DEC model (the third joint set, Jt, is parallel with this section and therefore does not appear here). The joint sets Jn and Js are specified and

3. Modeling procedure The discrete element method (DEM) [21], which was presented by Cundall in 1979, was used to simulate the response of discontinuous media to dynamic loadings. The jointed rock masses are considered as the composition of discontinuities and rock blocks which can have rotations and large displacements along discontinuities. Hence, the characterization of nonlinear and large deformations in jointed rock masses is simulated by DEM more realistically than other existing numerical methods. Because of the three-dimensional geometry of the problem, 3DEC was selected, which is a three-dimensional distinct element code which utilizes a Lagrangian calculation scheme to model large movements and deformations of a blocky system [34–36]. In this section, details of the dynamic discrete element analyses combined with the linear superposition method are described.

Fig. 3. Geometry of the 3DEC model of study area showing the three different regions, monitoring spot, and local coordinate system.


N. Babanouri et al. / Computers and Geotechnics 49 (2013) 134–142

Table 3 Dependence of normal and shear stiffness coefficients of faults on clay-filling thickness [38]. Clay-filling thickness (mm)

Kn (GPa/m)

Ks (GPa/m)


50–100 10–20 <1

0.1–0.5 0.5–2.0 >0.5

0.01–0.1 0.1–0.6 >1.0

6 3.6 5

Table 4 Laboratory properties of the clay filling. Kn (GPa/m)

Ks (GPa/m)

U (°)

c (kPa)





Fig. 4. A cross section of the 3DEC model showing faults and joint sets (faults are shown as bold lines).

faults are shown as bold lines in this figure. As can be deduced from Figs. 3 and 4, the study area is kinematically infeasible to slip. However, as exploitation operations proceed, rock blocks can be kinematically active. In order to consider the effect of intact rock deformation on wave propagation, blocks were assigned to be deformable. The Mohr–Coulomb elasto-plastic constitutive model was selected for blocks; values of required properties for different blocks are given in Table 1. To the faults and joints, the Coulomb slip model was assigned. Properties of joint sets along with their spacing are listed in Table 2. All mechanical properties given in Tables 1 and 2 were achieved by a large number of experiments at rock mechanics laboratory of the Gol-e-Gohar iron ore mine. In contrast, no test was performed on the faults. These faults, which are of particular importance in slope stability, have clay filling with a thickness ranging within 10–50 cm. Infanti and Kanji [38] reported the values of normal and shear stiffness of joints in relation to clay in-filling thickness at shallow depths (Table 3). The faults of the study area have differing thicknesses of clay filling, which are mainly greater than those listed in Table 3. On the other hand, since the thickness of the filling material was much greater than the amplitude of surface roughness, shear behavior of the fault is completely controlled by the in-filled material. Therefore, it is expected that the faults behave plastically and has very low stiffness coefficients. Hence, laboratory samples of infilled clay were prepared and then underwent mechanical tests. The results are presented in Table 4. In this study, based on values available in Tables 3 and 4, following ranges of search were considered for mechanical properties of faults (Table 5). 3.2. Mesh generation The blocks are meshed into triangular elements. The size of elements has a great influence on the computational results. A finer

Table 5 Search ranges of faults properties. Kn (GPa/m)

Ks (GPa/m)

U (°)

c (kPa)





mesh size results in higher accuracy. Kuhlemeyer and Lysmer [39] showed that for the accurate representation of wave transmission through a model, the spatial element size, Dl, must be smaller than approximately one-tenth to one-eighth of the wavelength, k, associated with the highest frequency component of the input wave. Alternatively, in order to accurately simulate the transmission of shock wave in rock joints and eliminate the wave distortion, the maximum input frequency of the wave should be less than a maximum frequency of fmax which is determined by the following equation [36,39]:

fmax ¼

c c ¼ ; k 10Dl


where c refers to the smaller value between Cp and Cs which are velocity of ‘P’ and ‘S’ waves, respectively, and their values can be obtained by the following formulae:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi K þ 4G 3 Cp ¼ ;


sffiffiffiffi G ; Cs ¼




where K, G and q are the bulk modulus, shear modulus and density, respectively. In this study, a uniform mesh of tetrahedral elements with the average edge length of 5 m was used. The elements size was selected fine enough to make sure that the waves with frequencies up to 20 Hz could be transmitted without distortion.

Table 1 Properties of blocks used in the 3DEC model. Material type

Density (kg/m3)

K (MPa)

G (MPa)

Elastic modulus (MPa)

Poisson’s ratio

U (°)

c (kPa)

Soil Quartz–schist Magnetite

2000 2660 4500

125 15,600 32,100

57.7 7200 14,800

150 18,720 38,485

0.3 0.3 0.3

40 36.2 36.9

15 11,000 20,500

Table 2 Characteristics of rock joints used in the 3DEC model. Rock type

Joint set

Kn (GPa/m)

Ks (GPa/m)

U (°)

c (kPa)

Spacing (m)

Quartz–schist Magnetite

Jn, Js, Jt Jn, Js, Jt

1.3 1.5

0.635 0.734

27 30

376 384

28, 22, 22 14, 11, 20


N. Babanouri et al. / Computers and Geotechnics 49 (2013) 134–142

3.3. Boundary conditions In static analyses, fixed or elastic boundaries can be realistically placed at some distance from the region of interest. In dynamic problems, however, such 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 minimize the problem, since material damping will absorb most of the energy of the waves reflected from distant boundaries. However, this solution leads to a large computational burden [36]. One alternative solution is to use absorbing boundaries for which several formulations have been proposed; the viscous boundary developed by Lysmer and Kuhlemeyer [39] was used in our simulations. It is based on the use of independent dashpots in the normal and shear directions at the model boundaries. The dashpots provide viscous normal and shear tractions given by [36]:

t n ¼ qC p v n ;


t s ¼ qC s v s ;


where vn and vs are the normal and shear components of the particle velocity at the boundary, q is the mass density, and Cp and Cs are the P- and S-wave velocities [36]. The bottom and side boundaries of the model were set to viscous boundaries to reduce the effects of reflected waves by dynamic load and, fixed in direction normal to each boundary. One restriction when applying particle velocity input to model boundaries is that this boundary condition cannot be applied to the same boundary along with a viscous boundary condition. To overcome this, a stress boundary condition can be used instead. A velocity wave may be converted to a stress wave using following formulae:

rn ¼ 2qC p v n ;


rs ¼ 2qC s v s ;


where rn and rs are applied normal and shear stresses. It should be noted that the factor of two in these equations is adapted to account for the fact that the applied stress must be doubled to overcome the effect of the viscous boundary [36]. In this study, the particle velocity histories in the proximity of the explosive source, calculated by the linear superposition method (as explained in the following section), after conversion to stress time histories were applied to the left side of the 3DEC model as the input dynamic loads. 3.4. Application of linear superposition method to determine dynamic input The determination of dynamic loads is the basic information and a difficult problem to study the dynamic response of rock slopes under explosion. The commonly used method is treating the load as a triangular pulse by giving the peak value, duration and rise time, using existing empirical formulae. However, the results acquired by this method are very different from the actual in situ conditions [26]. On the other hand, measurement of vibrations in the proximity of production blasts, where vibration level is very high, cannot be taken by conventional devices. In this case, one solution is to measure vibrations at a short distance from a single-hole shot and then model production blast vibrations using the linear superposition method. The linear superposition method is based on this principle: the measured seismogram at a given point is the result of linear superposition of the single-hole seismograms in time domain. This principle was initially developed in order to predict vibration waveforms radiated from a single column of explosive and to predict the total vibrations of a full-scale blast [11]. Due to linearity of

the problem and the superposition principle in which distributed sources can be described as the sum of multiple point sources, there is no more difficulty in modeling even a very complicated source. In addition to spectral amplitudes, all phase effects from the superposition are included in the synthetic seismogram. The part of the solution that put into relation the force distribution at the source with the displacements at the receiver is referred to as the elasto-dynamic Green’s function. Derivation of the Green’s function is the key step of the synthetic seismogram calculation. Indeed, this function must take into account all of the elastic properties of the material and the appropriate boundary conditions. Green’s function, G(x, t), gives the displacement at point x, that results from the unit force function applied at point x0 [40].

ui ðx; tÞ ¼¼ Gij ðx; t; x0 ; t0 Þfj ðx0 ; t 0 Þ;


Eq. (8) gives the displacement u from a realistic source, with the force vector or source time function, f, of row blast, synthesized from the displacement produced by the simplest possible source which is a unidirectional unit impulse localized precisely in space and time [40]. The ground velocity can be achieved by a differentiation of Eq. (8). If the displacement field of the row blast is a linear superposition of the displacements produced by the individual holes of the blast pattern, the source function f can be separated into two parts (by ignoring mathematical dependency of x, x0, t and t0) as follows:

f ¼ fs  fR ; fR ¼ ai dðt  ti Þ i ¼ 1; . . . ; N;

ð9Þ ð10Þ

where N is the number of charges and ti is the firing time of ith charge. In the above formulae, the shape of the measured parameter (displacement or velocity) is assumed to be identical for all individual holes. The term of fs is the source time function of a single hole, fR is an impulse series, and d is the delta function. The superposition of the individual signals is mathematically expressed as a convolution. Convolution with a delta function leaves the original function unchanged. The delta function may act to produce a time shift in the original time series. The amplitudes ai of the impulse in Eq. (10) are scaling factor for seismic effects of individual holes. The arriving times at the observation point of the body waves from detonation of the individual holes are expressed by the impulse position in the series. The displacement time history from a single blast at a specific location can be written as:

us ¼ fs  G;


The displacement time history can be measured by a field test. Combining Eqs. (2) and (3), the displacement of a complete row blast can be obtained:

u ¼ fs  fR  G;


Convolution is commutative and associative, then:

u ¼ fs  G  fR ;


and by using Eq. (11):

u ¼ us  fR ;


The impulse series, fR, can be calculated and the convolution in Eq. (14) combines field measurement and computer simulation. By Eq. (14) and based upon measuring of a single-hole motions, ground motions of a production blast can be calculated without calling the Green’s function. The procedure starts by drilling a single hole and loading it by a charge similar to the holes of the actual blast pattern. The displacement or velocity time history is then measured at locations for which the ground vibrations are to be

N. Babanouri et al. / Computers and Geotechnics 49 (2013) 134–142

Fig. 5. Simulated blasting waveforms by linear superposition method at a distance of 39 m from the blast.

predicated. The next step is calculation of the impulse series for each geophone position. This series is convolved with us to simulate the seismogram of the blast at the specific locations [41]. At the Gol-e-Gohar iron ore mine, available devices can record the particle velocity up to 254 mm/s. Therefore, a single-hole blast was carried out with the same parameters of production blasts and the time histories of particle velocities were measured at a distance of 39 m from the blast (at a nearer distance from the blast, vibration level may exceed 254 mm/s). Then, the linear superposition method was used to simulate the production blast time histories of particle velocities. For the blast of concern, there were four15hole rows, the inter-row delay was 40 ms, and the holes in each row were shot simultaneously. The explosive charge per each hole was 450 kg and the explosive charge per delay was then 6.75 tonnes. Simulated waveforms of the production blast are shown for each wave component in Fig. 5. The simulated time histories of particle velocity were then filtered at 20 Hz, and after conversion to the stress time histories, were applied to the left side of the model where is located at the same distance of 39 m from the actual location of the blast. Since the original waveforms had very low amplitudes at frequencies higher than 20 Hz, the filtration made the time histories of particle velocities slightly smooth such that the filtered and non-filtered waveforms were nearly the same.

3.5. Additional damping Natural dynamic systems contain some degree of damping of the vibration energy within the system; otherwise, the system would oscillate indefinitely when is subjected to driving forces. Damping is due, in part, to energy loss as a result of internal friction in the intact material and slipping along interfaces, if these are present. For a dynamic analysis, the damping in the numerical simulation should be reproduced in magnitude and form of the energy losses in the natural system subjected to a dynamic loading. In soil and rock, natural damping is mainly hysteretic, i.e., independent of frequency. It is difficult to reproduce this type of damping numerically because of at least two problems. First, many simple hysteretic functions do not damp all components equally when several waveforms are superimposed. Second, hysteretic functions lead to path-dependence which makes results difficult to interpret. However, if a constitutive model is found that contains an adequate representation of the hysteresis which occurs in a real material, then no additional damping would be necessary. The existing models are not considered to model dynamic hysteresis well enough to omit additional damping completely [36,42].


In time-domain programs, Rayleigh damping is commonly used to provide damping that is approximately frequency-independent over a restricted range of frequencies. The idea is to try to get the right damping for the important frequencies in problem. The required input parameters to specify Rayleigh damping in 3DEC are the central frequency, fc, which accounts for important frequency range, and the critical damping ratio, nc, which varies within 2–5% for geological materials and within 2–10% for structural systems [20,34,40,43]. In analyses which use one of the plasticity constitutive models (e.g., Mohr–Coulomb), a considerable amount of energy dissipation can occur during plastic flow. Thus, for many dynamic analyses involving large strain, only a minimal percentage of damping (e.g., 0.5%) may be required. Furthermore, dissipation will increase with amplitude for stress/strain cycles that involve plastic flow [36,42]. In this research, the central frequency was determined from smoothed Fourier amplitude spectra of recorded waveforms. However, assigning the correct value of additional damping is one of the most challenging steps in the dynamic analysis with numerical codes. The value of 4 Hz was determined for fc, and a search range of 0.5–5% was considered for nc. 4. Blast back analysis Each blast could be regarded as a very large scale dynamic loading test that valuable data on unknown characteristics and uncertainties of rock masses could be acquired by blast monitoring and back analysis of the recorded responses of rock masses. At the present study, by iterative comparison of the blast waveforms simulated by 3DEC and the measured seismograms under the supervision of a SA search algorithm, optimum values of the unknown parameters of the in-filled faults and the numerical model were obtained. In this study, all the simulations were carried out on a PC with Intel Core i7 @3.4 GHz CPU and 8 GB RAM. Depending on chosen values for the unknown properties of the numerical models, the run time of each model varied from 1 to 3 h. 4.1. SA search algorithm SA is a probabilistic search technique, which is based on the annealing/cooling process of metals [44]. Physical annealing is a process in which a solid is heated until all particles randomly arrange themselves in the liquid state, followed by a slow cooling, spending a relatively long time to reach the freezing point. At higher temperature, the liquidized solid is allowed to move freely and reach thermal equilibrium. As the temperature is decreased slowly rather than quickly, nature is able to find the minimum energy level (freezing point) [45]. The traditional optimization methods are similar to the quickly cooled process for a liquid metal, greedily searching for the minimum in the downhill direction from the initial starting points. The problem is treated as searching for values of the variables so that the cost function (C) has its minimum value. In a simulated annealing process, if DC, defined as C(vk+1)  C(vk) in two consecutive random trials, is less than zero, vk will be replaced with vk+1. However in order to escape local minima, the algorithm must be allowed uphill moves. In this way, if DC > 0, vk+1 may still be accepted based on an auxiliary judgment. A random number, rk, is generated in (0, 1) and compared to the value exp[DC/T], where T is a specified parameter simulating temperature. If rk < exp[DF/T], vk will be replaced with vk+1 [45]. The parameter T slowly controls this searching process and plays the role of the temperature similar to that in physical annealing. The process starts from an initial temperature, T0, which is set sufficiently high to allow uphill moves away from local minimal


N. Babanouri et al. / Computers and Geotechnics 49 (2013) 134–142

[44]. However, too high a temperature will result in deteriorating random steps that lower the efficiency. The starting temperature T0 is then defined by:

T 0 ¼ F max  F min ;

was searched within the range of 0.5–5%. Targets were: peak particle velocities, and the values of dominant frequencies of each


where Fmax and Fmin are respectively the estimated maximum and minimum of the cost function. The temperature during the optimization is reduced by a damping function, that is,

T k ¼ T 0 ðaT Þk ;


where Tk is the value of temperature at the kth step of cooling process, aT is referred to as the temperature damping factor and falls in the range (0, 1). The magnitude of random walking, Dr, decreases with temperature by a damping function. The magnitude of random walking Drk at Tk is defined as:

Dr k ¼ Dr 0 ðar Þk1 ;


where ar is the damping factor of random walking in the range (0, 1). Dr0 is the initial value of the random walking, which is determined based on the properties of the cost function and the desired accuracy and resolution [45]. A pseudo-code for the SA algorithm used in this study is given in Fig. 6. The cost function (C), which was to be minimized, was the average of deviations of simulated peak particle velocities and dominant frequencies from their real values, in percent. The parameters of the SA algorithm used were set as T0 = 1000, aT = 0.9, Dr0 = 20% of range length for each unknown properties, and ar = 0.85. 4.2. Results The unknown geomechanical parameters in the numerical model were: elastic properties of the faults namely normal and shear stiffness coefficients, and plastic properties of faults namely cohesion and friction angle. A search range was considered for each of these parameters as mentioned in Table 5. The unknown parameter of numerical model was the additional damping ratio which

Fig. 7. Transversal component of particle velocity obtained by field measurements and by the 3DEC model: (a) in time domain and (b) in frequency domain.

Fig. 6. Pseudo-code for the SA algorithm used in this study.


N. Babanouri et al. / Computers and Geotechnics 49 (2013) 134–142

wave component at the monitoring spot which was located at the distance of 300 m, form the blast. Figs. 7–9 show blasting particle velocity components in time and frequency domains obtained by field measurements and by the 3DEC model at the monitoring point in the end of 20th iteration. As shown in parts ‘a’ of Figs. 7–9, calculated values of PPVs are in a good accordance with their real results (the maximum deviation is 15%). Values of predominant frequencies obtained in 3DEC model and site monitoring are also very close as can be seen in parts ‘b’ of Figs. 7–9 (the maximum deviation is 10%). Although vibration duration was not considered within targets (because it is of less importance compared to the amplitude and predominant frequency), but as can be seen in Figs. 7a, 8a, and 9a, durations simulated by 3DEC are very close to those recorded in the field. As a result of wave reflections, in both the real and simulated waveforms the low-amplitude waves are observed after 1 s when excitation was stopped. It is concluded that 3DEC well simulated the wave reflection during the passage of wave through the jointed media. When the search algorithm was proceeding, it was observed that variation of c and U of faults within their search ranges had no influence on 3DEC model responses. This is because the dynamic blasting loads were not strong enough to slide a large number of blocks along the faults so that they would exceed the elastic displacements. In this case, no values could be attributed to plastic parameters of faults. The optimized values of other parameters in the end of 20th iteration are presented in Table 6. As can be seen, the optimum value of additional damping is very low which indicates the material properties, the constitutive models, and the geometry of the 3DEC model were realistic enough to simulate the damping without a large amount of additional damping.

Fig. 9. Longitudinal component of particle velocity obtained by field measurements and by the 3DEC model: (a) in time domain and (b) in frequency domain.

Table 6 Optimized values of parameters in the 3DEC model. Parameter

Kn (GPa/m)

Ks (GPa/m)

n (%)

Optimum value




Values obtained for elastic properties of faults (Kn and Ks) are lower than those of in-filled clay material. Thus, assigning laboratory characteristics of in-filled clay to fault properties leads to over-estimation of stability of the Gol-e-Gohar iron ore mine walls where instabilities are mostly controlled by the faults.

5. Conclusions

Fig. 8. Vertical component of particle velocity obtained by field measurements and by the 3DEC model: (a) in time domain and (b) in frequency domain.

Numerical simulation of blast wave propagation through jointed rock masses is adversely affected by the factors such as difficulties of determining realistic dynamic input, weakness of existing constitutive models in reproducing the damping, and uncertainties associated with mechanical properties of rock masses. These factors could make numerical dynamic analyses absolutely inefficient and unrealistic. To overcome these shortcomings and enhance numerical dynamic analyses, three dimensional discrete element method was combined with other techniques in this paper. First, the linear superposition method was used to calculate dynamic input based on superposition of single-hole vibrations at the Gol-e-Gohar iron ore mine. Then, to any unknown parameter, instead of a fixed value, a range was assigned within which optimum value of the unknown parameter was back calculated by iterative comparison of simulated blast waveforms with the real seismograms. To find the


N. Babanouri et al. / Computers and Geotechnics 49 (2013) 134–142

optimum configuration of unknown parameters values efficiently and accurately, a simulated annealing algorithm was employed to supervise the iterations, in which targets were PPVs and predominant frequencies of waveforms monitored at a distance of 300 m form the blast. The unknown parameters were properties of in-filled faults and additional damping ratio, values of PPVs and predominant frequencies in the model and site monitoring were in a good accordance in the end of 20th iteration (maximum differences were 15% and 10%, respectively). Furthermore, durations of vibration simulated by model were very close to those recorded in the field that indicated the model well simulated the wave reflection. The results showed that the plastic parameters (c and U) of faults did not influence the model response, and no values could be, therefore, attributed to them in this case. On the other hand, values obtained for elastic properties of faults (Kn = 0.114 GPa/m and Ks = 0.025 GPa/m) were lower than those of in-filled clay material. Thus, application of laboratory characteristics of in-filled clay as fault properties leads to over-estimation of stability of the Gole-Gohar iron ore mine. The value of additional damping was very low (1.12%) which indicates that the material properties, constitutive models, and the model geometry were realistic enough to simulate the damping. Estimation of unknown properties of rock media could be made based on back analysis of the results obtained by blast monitoring. In fact, every blast could be regarded as a very large-scale loading test, which can provide valuable data. It seems that no significant attention has been yet paid to this side application of production blasting at mines. Although this study has dealt with the issue in a case study, the presented methodology could be used in order to study propagation of any blast wave provided that enough monitoring records are available. Acknowledgements The authors would like to thank Mr. H. R. Mohammadi, in rock mechanics division of Gol-e-Gohar iron ore mine, for providing the required site information. Special appreciation is also extended to Ms. S. Sarafrazi for her advices on the application of simulated annealing algorithm. References [1] Miller RK. The effects of boundary friction on the propagation of elastic waves. Bull Seismol Soc Am 1978;68:987–98. [2] Day SM. Test Problem for Plane Strain Block Motion Codes. S-Cubed Memorandum to Itasca; 1985. [3] Siskind DE, Stagg MS, Kopp JW, Dowding CH. Structure response and damage produced by ground vibration from surface mine blasting. US Bur Mines Rep Invest 8507, 1980. 74pp. [4] Ghosh A, Daemen Jack JK. A simple new blast vibration predictor (based on wave propagation laws). In: 24th US symposium on rock mechanics; 1983. [5] Taqieddin SA. Ground vibration levels: prediction and parameters. Min Sci Technol 1986;3:111–5. [6] Gupta RN, Roy PP, Bagchi A, Singh B. Dynamics effects in various rock mass and their predictions. J Mines, Metals Fuels 1987;35:455–62. [7] Wu YK, Hao H, Zhou YX, Chong K. Propagation characteristics of blast-induced shock waves in a jointed rock mass. Soil Dyn Earthquake Eng 1998;17:407–12. [8] Hakan A, Adnan K. The effect of discontinuity frequency on ground vibrations produced from bench blasting: a case study. Soil Dyn Earthquake Eng 2008;28:686–94. [9] Kuzu C. The importance of site-specific characters in prediction models for blast-induced ground vibrations. Soil Dyn Earthquake Eng 2008;28:405–14. [10] Dowding CH. Construction vibration. New Jersey: Prentice Hall, Upper Saddle River; 1996. [11] Hinzen KG. Modelling of blast vibrations. Int J Rock Mech Min Sci Geomech Abstr 1988;25:439–45. [12] Li Y. Underground blast induced ground shock and its modelling using artificial neural network. Comput Geotech 2005;32:164–78.

[13] Mohamed MT. Artificial neural network for prediction and control of blasting vibrations in Assiut (Egypt) limestone quarry. Int J Rock Mech Min Sci 2009;46:426–31. [14] Khandelwal M, Singh TN. Prediction of blast induced ground vibrations using artificial neural network. Int J Rock Mech Min Sci 2009;46:1214–22. [15] Khandelwal M. Evaluation and prediction of blast induced ground vibration using support vector machine. Int J Rock Mech Min Sci 2010;47:509–16. [16] Monjezi M, Ahmadi M, Sheikhan M, Bahrami A, Salimi AR. Predicting blastinduced ground vibration using various types of neural networks. Soil Dyn Earthquake Eng 2010;30:1233–6. [17] Sigre Y. Maitrise des tirs par modélisation sismique hybride à Val d’Azergues. Industrie minérale; 1994. [18] Mansouri H, Ebrahimi Farsangi MA. Blast vibration modeling in sar cheshmeh mine. In: 19th World mining congress; 2003. [19] Toraño J, Ramírez-Oyanguren P, Rodríguez R, Diego I. Analysis of the environment effects of ground vibrations produced by blasting in quarries. Int J Min Reclam Environ 2006;20(4):249–66. [20] Toraño J, Rodríguez R, Diego I, Rivas JM, Casal MD. FEM models including randomness and its application to the blasting vibrations prediction. Comput Geotech 2006;33(1):15–28. [21] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Geotechnique 1979;29:47–65. [22] Eberhardt E, Stead D. Mechanisms of slope instability in thinly bedded surface mine slopes. In: Proceedings of 8th congress of the international association for engineering geology and, environment; 1998. p. 3011–8. [23] Stead D, Eberhardt E, Coggan JS. Developments in the characterization of complex rock slope deformation and failure using numerical modeling techniques. Eng Geol 2006;83:217–35. [24] Bhasin R, Kaynia AM. Static and dynamic simulation of a 700 m high rock slope in western Norway. Eng Geol 2004;71:213–26. [25] Kveldsvik V et al. Dynamic distinct-element analysis of the 800 m high Aknes rock slope. Int J Rock Mech Min Sci 2009;46:686–98. [26] Liu YQ, Li HB, Zhao J, Li JR, Zhou QC. UDEC simulation for dynamic response of a rock slope subject to explosions. (SINOROCK2004 symposium). Int J Rock Mech Min Sci 2004;41:474. [27] Fan SC, Jiao YY, Zhao J. On modelling of incident boundary for wave propagation in jointed rock masses using discrete element method. Comput Geotech 2004;31:57–66. [28] Lei WD, Hefny AM, Yan S, Teng J. A numerical study on 2-D compressive wave propagation in rock masses with a set of joints along the radial direction normal to the joints. Comput Geotech 2007;34:508–23. [29] Zhao XB, Zhao J, Cai JG, Hefny AM. UDEC modelling on wave propagation across fractured rock masses. Comput Geotech 2008;35:97–104. [30] Chen SG, Zhao J. A study of UDEC modeling for blast propagation in jointed rock masses. Int J Rock Mech Min Sci 1998;35:93–9. [31] Chen SG, Cai JG, Zhao J. Discrete element modeling of an underground explosion in a jointed rock mass. Geotech Geol Eng 2000;18:59–78. [32] Zhao J, Cai JG. Transmission of elastic P-waves across single fractures with a nonlinear normal deformational behavior. Rock Mech Rock Eng 2001;34:3–22. [33] Kulatilake PHSW, Wang S, Stephansson O. Effect of finite size joints on the deformability of jointed rock in three dimensions. Int J Rock Mech Min Sci Geomech Abstr 1993;30:479–501. [34] Cundall PA. Formulation of a three-dimensional distinct element model-part I. A scheme to detect and represent contacts in a system composed of many polyhedral blocks. Int J Rock Mech Min Sci Geomech Abstr 1988;25:107–16. [35] Hart R, Cundall PA, Lemos J. Formulation of a three-dimensional distinct element model-part II. mechanical calculations for motion and interaction of a system composed of many polyhedral blocks. Int J Rock Mech Min Sci Geomech Abstr 1988;25:117–25. [36] Itasca Consulting Group Inc. 3DEC user’s guide, Minneapolis; 1999. [37] Lichun J, Liuqing H, Xiuying L. Investigation on the threshold control of safety blasting vibration velocity for the extraction of complicated orebody under railway. Min Sci Technol 2011;21:169–74. [38] Infanti N, Kanji MA. In situ shear strength, normal and shear stiffness determinations at Arua Vermelha project. In: Proceedings of third international congress of the international association of engineering geology, Madrid, vol. 2. 1978. p. 175–83. [39] Kuhlemeyer RL, Lysmer J. Finite element method accuracy for wave propagation problems. J Soil Mech Found Div ASCE 1973;99(SM5):421–7. [40] Blair DP, Minchinton A. Near-field blast vibration models. In: Eighth international symposium rock fragmentation by blasting, Santiago, Chile; 2006, pp. 152–9. [41] Stump BW, Reinke RE. Experimental confirmation of superposition from smallscale explosions. Bull Seismol Soc Am 1983;78:1059–73. [42] Itasca Consulting Group Inc. FLAC user’s guide, version 4.0, Minneapolis; 2002. [43] Shin JH, Moon HG, Chae SE. Effect of blast-induced vibration on existing tunnels in soft rocks. Tunn Undergr Space Technol 2011;26:51–61. [44] Kirkpatrick S, Gellato CD, Vecchi MP. Optimization by simulated annealing. Science 1983;220:671–80. [45] Chen Z, Wang J, Wang Y, Yin JH, Haberfield C. A three-dimensional slope stability analysis method using the upper bound theorem Part II: numerical approaches, applications and extensions. Int J Rock Mech Min Sci 2001;38:379–97.