- Email: [email protected]
On modelling of incident boundary for wave propagation in jointed rock masses using discrete element method
Computers and Geotechnics 31 (2004) 57–66 www.elsevier.com/locate/compgeo
On modelling of incident boundary for wave propagation in jointed rock masses using discrete element method S.C. Fan *, Y.Y. Jiao, J. Zhao Protective Technology Research Centre, School of Civil and Environmental Engineering, Nanyang Technological University, Nanyang Avenue 639798, Singapore Received 7 January 2003; received in revised form 17 October 2003; accepted 20 November 2003
Abstract Numerical modelling is an effective alternative for studying the mechanical behaviour of rock masses. Very often, in solving the same problem with the same code, different models yield remarkably different results. There are many reasons. Amongst them, the boundary conditions are often overlooked but it has significant influence on the stress wave propagation in jointed rock masses. Usually, in solving dynamic problems, the incident boundaries are described either by stress-history input (SHI) or velocity-history input (VHI). This paper presents the study on the effect of these two approaches using discrete element method (DEM) on wave propagation and wave attenuation in jointed rock masses. Two cases were designed to illustrate the points. The first one shows onedimensional P-wave propagation along a rock bar having an internal cross joint. The second one shows the blasting wave propagating through jointed rock masses. The blasting wave is induced by an underground explosion. The rock bar/masses are discretized into elements and the analyses adopt the DEM, which is available in the code ÔUDECÕ. The results reveal that the VHI approach can yield results in accord with field-test data, while the SHI approach may lead to unreasonable results. The two approaches may yield the same results for homogeneous rock, but definitely not the same for jointed rock masses. In addition to the numerical illustrations, this paper also presents some insights of the reason behind it. 2003 Elsevier Ltd. All rights reserved. Keywords: Discrete element method; Wave propagation; Jointed rock; Boundary effects
1. Introduction In rock engineering problems, handling the discontinuities in rock masses is a challenging task. The discontinuities appear in the form of faults, joints or bedding planes. The presence of the discontinuities in rock masses has great influence on the responses to either static or dynamic loads [1], and renders the numerical simulations more complicated. It demands special treatment for the discontinuities, namely joints. There are many numerical tools available. The popular ones include the finite element method (FEM), boundary element method (BEM), finite difference method (FDM), and discrete element method (DEM), just name a few. In FEM, joints are modelled using a kind of
*
Corresponding author. Tel.: +65-6792-1620; fax: +65-6792-1650. E-mail address: [email protected] (S.C. Fan).
0266-352X/$ - see front matter 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2003.11.002
special ‘‘joint elements’’, which simulate the interface discontinuity [2,3]. In FDM, joints are simulated using slide-lines [4]. Very often, BEM is employed for the far field across the joint while FEM is used for near field, and their boundary interface is specified according to the properties of joint [5]. The above-mentioned treatments are efficient when there are only a few joints. For heavily jointed rock masses, an alternative is to treat it as a continuum having equivalent properties [6]. In fact, the equivalent continuum may behave like a jointed rock masses under some circumstances, particularly under static loads. Generally, they fail to represent the dynamic responses and no equivalent properties could be found. On the other hand, the discrete element method (DEM), which is developed to model discontinuous problems, has been recognized to be a better alternative [7]. In DEM, a rock mass is treated as an assemblage of discrete blocks by fictitious or real joints, of which some are the contact interfaces between distinct bodies.
58
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
Individual block can be defined either as rigid or deformable by specifying the material parameters. At the contact interfaces, constitutive laws governing the force– displacement relationship are specified. The solution scheme is identical to that used by the explicit FDM for continuum analysis. The dynamic behaviour is represented numerically by a temporal stepwise algorithm. By using the explicit dynamic approach, DEM is able to simulate large displacements across the discontinuities without encountering much numerical stability problems. The DEM algorithm was coded and available commercially (code-named ÔUDECÕ) [8]. Among the numerical analysts, saying has it that ‘‘garbage in and garbage out’’. Very often, in solving the same engineering problem with the same code, different models yield remarkably different results. There are many reasons: Putting aside those due to inexperienced users, there are many physical phenomena needed to be addressed. Zukas and Scheffler [9,10] investigated the influence of mesh configuration and material interface on numerical results; Johnson [11] studied the effects of introducing an artificial viscosity when using the method called smooth particle hydrodynamics (SPH); Chapman et al. [12] reported their parametric studies on blast wave simulation. Usually, when solving dynamic problems, there can be three kinds of incident boundary (which is defined as boundary that accepts dynamic load), namely displacement-history input (DHI), stress-history input (SHI), and velocity-history input (VHI). Zhao and Valliappan [13–16] proposed a DHI-based wave input procedure in the finite and infinite coupled model to simulate wave-scattering problems in infinite media. They showed that the DHI approach was equivalent to the then prevailing wave-input procedures and had the merit of higher computational efficiency; Lemos [17] employed DEM to investigate the effect of discontinuity in a jointed rock mass subjected to dynamic load generated from a line source. LemosÕ model has only one single discontinuity, and the effects of the two approaches (VHI and SHI) on the resulting slip along the discontinuity were reported. Effects due to other factors such as damping, in situ stress, and material model, etc. were recorded in other literatures [8,17,18]. Amongst the above-mentioned factors, the incident boundary conditions are often overlooked, but it has significant influence on the stress wave propagation in jointed rock masses. This paper presents the study on the effect of the two popular approaches (VHI and SHI) on wave propagation and wave attenuation in jointed rock masses. Two cases were designed to illustrate the points. The first one shows the one-dimensional P-wave propagation along a rock bar having an internal cross joint. The second one shows the blasting wave propagating through jointed rock masses. The blasting wave is induced by an underground explosion.
2. Modelling of boundary conditions In reality, engineering problems always involve infinite or very large domain. Many numerical techniques have been developed to handle those situations, such as special infinite elements in FEM and infinite domain in BEM. Nevertheless, in solving dynamic problems, be it for a finite or infinite domain, the boundary conditions need to be defined. It could be reflective or non-reflective. When the DEM is employed, the domain of interest needs to be limited to a finite size, and some kinds of artificial boundaries are introduced. The artificial boundaries are usually modelled as non-reflective boundaries. Besides, there are natural boundaries which are reflective such as ground surface, excavation surfaces, construction surfaces, etc. In general, boundary conditions can be classified under three categories: prescribed force/stress (PS) boundary; prescribed displacement/velocity (PV) boundary; or boundary with both prescribed stress and velocity. The reflective or non-reflective boundary falls in one of these categories. What follows are treatments for various boundary condition. 2.1. Incident boundary In the dynamic analysis, there are two prevailing approaches of inputting the dynamic excitation, namely VHI and SHI. Correspondingly, the incident boundaries are PV and PS boundary, respectively. In the simulation using DEM, the temporal incremental equation of motion is ðtþDt=2Þ
u_ i
ðtDt=2Þ
¼ u_ i þ
nX
X o Dt ðtÞ ðtÞ ðtDt=2Þ ; F i þ a Þ Fi signðu_ i m ð1Þ
ðtþDtÞ
ðtDtÞ
in which u_ i and u_ i denote the velocity at gridpoint i at time t þ Dt and t Dt, respectively, Fit is the resultant force at grid-point i at time t, a is the damping coefficient, and m is the nodal mass. The primary variables of Eq. (1) are velocities. Therefore, when SHI approach is adopted, the incident boundary is PS type and its implementation is not straightforward. Firstly, the equivalent force at grid-point has to be evaluated, i.e., Fib ¼ rbij nj Ds;
ð2Þ
where nj is the outward normal vector of the boundary segment/surface and Ds is the length/area of the b boundary segment/surface on which the P stress-history rij is imposed. Then the total force Fi acted on the boundary grid-point i is the sum of all applied forces and those derived from internal or inter-block stresses, i.e., Fi ¼ Fib þ Fic þ Fiz þ Fig ;
ð3Þ
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
in which Fic is the resultant contact forces between blocks and exists only at grid-points along the block boundary; Fiz denotes the contribution from internal stresses in the adjacent zone around the grid-point. Fiz can be calculated through Z Fiz ¼ rij nj ds; ð4Þ C
where rij is the zone stress tensor; nj is the unit outward normal to the contour C [8]; Fig is the body force due to gravity, i.e., Fig ¼ mg gi ;
ð5Þ
where P mg is the lumped mass at the grid-point. The value of Fi is then applied in Eq. (1) to calculate resultant grid-point velocity. On the other hand, when VHI approach is adopted, it becomes very straightforward and it is not necessary to invoke Eq. (1) because the primary variables are velocities. The velocity at the boundary grid-point in each subsequent time step can be updated directly form the input data, i.e., ðtþDt=2Þ
u_ i
ðtþDt=2Þ
¼ U_ i
;
ð6Þ
ðtþDt=2Þ
is the velocity-history input, described as where U_ i a function of time. 2.2. Other boundaries Other boundaries include natural boundaries and artificial boundaries. A natural boundary could be ground surface, excavation surface or other construction surface, which is often a reflective boundary and treated in the same way as the PS boundary mentioned above. For the artificial boundaries, they are often non-reflective so that they can silent the outward propagating wave. In the current DEM simulation, two types of artificial boundaries are employed: viscous (non-reflective) boundaries and free-field boundaries [8]. At the viscous boundaries, independent linear viscous dashpots are put in place in the normal and shear directions to absorb any waves coming at incident angles greater than 30. At the free-field boundaries, special boundary elements representing the extended infinite medium are put in place so that all outgoing waves are completely absorbed. It is worth noting that the non-reflection boundaries are PS boundaries where stresses are prescribed. However, the stresses are usually not known in advance. An approximate technique is to derive from the velocities as follows: rn ¼ 2ðqCp vn Þ;
ð7Þ
rs ¼ 2ðqCs vs Þ;
ð8Þ
59
in which rn and rs are the normal and shear stress components, respectively; q is the mass density; Cp and Cs are, respectively, the speed of P-wave and S-wave propagating through the medium; and vn and vs are the normal and transverse particle velocity, respectively. Note that a factor of two is applied to double the stresses in order to overcome the effect of the viscous boundary.
3. Consequence of VHI and SHI approaches The effect of VHI and SHI approaches adopted in DEM simulation of wave propagation and wave attenuation in jointed rock masses are investigated through two case-studies. The analyses were carried out using the UDEC code. 3.1. One-dimensional P-wave propagation along a rock bar To fend off disturbance arising from other factors, the P-wave propagation through a simple one-dimensional rock bar is computed to illustrate the consequence due to different approaches (VHI and SHI). Two scenarios are considered. One is a perfectly elastic rock bar without any joint inside. The other is the same rock bar but having a transverse joint at the mid-length. The length and height of the rock bar are 10.0 and 1.0 m, respectively. The properties of the parent rock bar are: density q ¼ 2650 kg/m3 , bulk modulus K ¼ 44:1 GPa, and shear modulus G ¼ 33:8 GPa. The properties of joint are: normal stiffness kn ¼ 100 GPa/m, shear stiffness ks ¼ 100 GPa/m, cohesion C ¼ 2 MPa, friction angle / ¼ 35, and no tensile strength is considered. No damping is taken into account (DAMP 0 0). In order to simulate the one-dimensional P-wave propagation, viscous boundaries are specified at the upper, lower, and right boundaries. A train of sinusoidal P-velocity waves for VHI approach and corresponding sinusoidal P-stress waves for SHI approach are applied to the left boundary of the rock bar, respectively. The amplitude of the incident velocity wave is 0.065 m/s; the frequency is 200 Hz; and the whole duration is 0.02 s. For the VHI approach, the input is straightforward. But for the SHI approach, the equivalent stress-history is derived from the velocity-history, or vice versa: rn ¼ qCp vn ;
ð9Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where the speed of P-wave Cp ¼ ðK þ 4G=3Þ=q, and thus the amplitude of the stress wave is 1.0 MPa. The resulting velocity-histories at 10 points along the bar are computed. The computational model is shown in Fig. 1. The five points to the left of the joint are denoted by Ai ði ¼ 1; . . . ; 5Þ and the five points to the right of the joint are denoted by Bi ði ¼ 1; . . . ; 5Þ.
60
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
Fig. 1. One-dimensional P-wave propagating along a rock bar.
Fig. 3. Resulting velocity-histories by VHI approach for scenario 2: having a joint. (a) At five points to the left of joint. (b) At five points to the right of joint.
Fig. 2. Resulting velocity-histories at point A1 and B5 for scenario 1: no joint. (a) Results by SHI approach. (b) Results by VHI approach.
• Scenario 1. In scenario 1, no joint exists. Regardless of the VHI or SHI approach being used, the resultant velocity waveform arriving at the far end (point B5 ) was found similar to the input waveform. Little distortion was observed (see Fig. 2(a) and (b)). • Scenario 2. In Scenario 2, there is a transverse joint at the mid-length. First, VHI approach is adopted and the resulting velocity-histories at the 10 points are traced. Fig. 3(a) shows the velocity-histories at the five points (A1 ; . . . ; A5 ) to the left of the joint and Fig. 3(b) for five points (B1 ; . . . ; B5 ) to the right of the joint. It can be seen from Fig. 3(b) that the negative parts of the waveform almost vanished immediately after travelling across the joint. On the other hand, positive parts of the waveform cross the joint with only half waveforms survived. It is worth noting that the incident waves are reflected back from the joint, the reflected waves superimpose with the incident waves (see Fig. 3(a)), and it results in magnified
waveforms at points A2 ; . . . ; A5 after initial 3/4 of a period. In addition, the farther the point is away from the joint, the weaker the magnification effect is. For easy visualization of the consequence, the velocityhistories at point A1 and B5 are plotted together in Fig. 4. Note should be taken that the waveform at
Fig. 4. Resulting velocity-histories by VHI approach for scenario 2: having a joint.
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
point A1 does not change because Eq. (1) is not invoked at the incident boundary. Subsequently, the SHI approach is employed. Results are different from those by VHI approach. Fig. 5(a) portrays the velocity-histories at the five points to the left of the joint. It can be seen that the waveforms are totally distorted due to superimposition of the reflected waves and the incident waves. In addition, the mean velocities all drift away from the initial zero value. It is also observed that the resulting velocity waves at A1 ; . . . ; A5 are all magnified in the process of propagation. The magnification at point A1 is the least, and velocities at A1 become negative at late time. Fig. 5(b) plots the resulting velocity-histories at the five points to the right of the joint. It can be seen that tensile parts of the waveforms are also filtered off by the joint, whereas the compressive parts of waveforms only have one half waveform survived after travelling through the joint. For easy visualization, Fig. 6 shows the contrast of velocity-histories between point A1 and B5 . Fig. 7 shows the maximum axial displacements obtained from the two approaches. It can be seen that SHI approach leads to relatively larger backward displacements than the VHI approach at locations to the left of the joint. Furthermore, SHI approach yields the maxi-
61
Fig. 6. Resulting velocity-histories by SHI approach for scenario 2: having a joint.
Fig. 7. Maximum axial displacements obtained from the two approaches.
mum backward displacement at the incident boundary, while the VHI approach results in steady zero displacement at the incident boundary. Figs. 8 and 9 plot the stress-histories of rx for locations to the left of the joint for the SHI and VHI approaches, respectively. As can be seen, waveforms in
Fig. 5. Resulting velocity-histories by SHI approach for scenario 2: having a joint. (a) At five points to the left of joint. (b) At five points to the right of joint.
Fig. 8. Resulting rx -histories by SHI approach for scenario 2: having a joint.
62
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
through the joint efficiently. To the contrary, in the SHI approach, the wave reflection from the joint allows the backward rigid-body motion and leads to opening of the joint. As such, most of the energy returns back to the left part of the bar. The mechanism will be discussed in Section 4. 3.2. Blasting wave propagation through jointed rock masses
Fig. 9. Resulting rx -histories by VHI approach for scenario 2: having a joint.
both are distorted beyond the initial 3/4 period, and the stress amplitude decreases when it is further away from the incident boundary. In the SHI approach, the resultant stress wave is identical to the original incident wave, and the waveforms at different locations appear similar to each other. In the VHI approach, the resultant stress waves are magnified and it remains so at the incident boundary even after 3/4 period; and the waveforms appear much more disturbed. These phenomena are reasonable and apprehensible. At the joint, the arrival waves are reflected and then the reflected waves superimposes with the incoming waves. In the SHI approach, the incident boundary cannot be sufficiently constrained so that a backward rigid-body movement becomes possible and it leads to stress relaxation at locations to the left of the joint. On the other hand in the VHI approach, the movement of the incident boundary follows exactly the prescribed input waveform and therefore the unwanted backward rigid-body movements are eliminated. As a result, the disturbance and magnification of the waveforms can be revealed. The two scenarios above clearly demonstrate the effect of the way in defining the incident boundary conditions when joints exist in the fractured rock masses. In responding to different incident boundary conditions, the manner of wave propagation is different due to the wave reflection from the joint, and it results in significantly disturbed wave form. Since DEM is a velocitybased method, the numerical results will be sensitive to the prescribed incident boundary conditions and the behaviour at the joint as well. As a result, the outcome is input-dependent. In the VHI approach in scenario 2, the incident boundary is constrained to the prescribed position with respect to time (in accordance of the input waveform), the backward rigid-body motion of the left part of the bar is not allowed to take place when the wave reflects back from the joint, and thus the joint doss not open up. Consequently, the wave is able to transmit
A small-scale field test of ground blasting was reported in literature [19]. The test results are used as benchmark against the results obtained by the current study using VHI and SHI approaches. The explosion chamber, having the dimensions of 8 4 2 m, is located at 115 m below the ground surface. The effective TNT charge weight is W ¼ 606 kg with a loading density of 10 kg/m3 . A DEM model was constructed accordingly as shown in Fig. 10. The analyses are carried out using the UDEC code. In the DEM model, there are two sets of joints: one set dipping at 45 with a spacing of 2 m and the other dipping at 45 with a spacing of 2 m. The parent rock material is assumed linear elastic, while all joints satisfy the Coulomb-slip constitutive model. The properties of parent rock are: density q ¼ 2650 kg/m3 , bulk modulus K ¼ 41 GPa, and shear modulus G ¼ 31 GPa. The properties of joints are: normal stiffness kn ¼ 100 GPa/ m, the shear stiffness ks ¼ 100 GPa/m, cohesion C ¼ 2 MPa, friction angle / ¼ 25, and no tensile strength is considered. At the artificial boundaries (bottom, left, and right), ÔquietÕ-type viscous boundaries are put in place to absorb the outgoing waves. The ceiling and the
Fig. 10. DEM computational model for the ground blasting test.
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
63
walls around the chamber are the incident boundaries. It should be noted that in dynamic modelling, zero or small damping is usually adopted [8]. To study the effect, two different damping ratios (zero and a relatively small ratio of 0.01) are used in the simulations for both approaches (VHI and SHI). Totally, four (2 2) combinations are considered. The blast load, represented by a triangular overpressure-history on the chamber walls, is obtained by the following formulae [20]: 14:0717 5:5397 0:3572 þ 2 3 R R R 0:00625 6 0:3; þ 4 ½kg=cm2 0:05 6 R R 6:1938 0:3262 DpU ¼ 2 R R 2:1324 6 1; þ 3 ½kg=cm2 0:3 6 R R 0:662 4:05 3:288 6 10; DpU ¼ þ 2 þ 3 ½kg=cm2 1 6 R R R R
DpU ¼
ð10Þ DpUr ¼ 2DpU þ
6DpU2 ; DpU þ 7p0
s þ 0:264R 2 0:129R 3 ffiffiffiffiffi ¼ 103 ð0:107 þ 0:444R p 3 W 4 Þ 0:05 6 R 6 3 ½s=kg1=3 ; þ 0:0335R
ð11Þ
ð12Þ
where Dp/ is the maximum overpressure at the blasting wave front; Dp/r is the reflected wave overpressure at the chamber walls; s (measured in second) is the duration of the overpressure; R (measured in metre) is the distance pffiffiffiffiffi ¼ R= 3 W of a point away from the charge centre; R (unit ¼ m/kg1=3 ) is the scaled distance; W (measured in kg) is the charge weight; and p0 (measured in kg/cm2 ) is the atmospheric pressure. From the test configuration, with the loading non-uniformity taken into account, Eq. (11) yields a maximum overpressure equal to 30.23 MPa, and Eq. (12) yields the duration equal to 2.5 ms. This triangular form of stress-history is used directly in the SHI approach. But for the VHI approach, the equivalent velocity-history needs to be computed. The approximation formula (Eq. (9)) is not adopted here. Instead, a more accurate method is used. The velocityhistory at the incident boundary is obtained through a separate DEM analysis on the same model without joints. It yields an equivalent velocity-history of similar triangular form having a maximum velocity of 2.03 m/s over duration of 2.5 ms. Fig. 11(a) plots the resulting velocity-histories obtained by VHI approach against field test data at 8 m above the detonation (using damping ratio of 0.01), while Fig. 11(b) shows the results obtained from SHI approach (using damping ratio of 0.01). We first focus
Fig. 11. Resulting velocity-history at a location 8 m above detonation. (a) Results obtained from VHI approach (damping ratio ¼ 0.01). (b) Results obtained from SHI approach (damping ratio ¼ 0.01).
on the peak velocity, which is considered a key factor controlling material damage. The VHI approach yields a peak velocity of 0.989 m/s, which is close to the field test result. But the SHI approach yields a peak velocity of 0.723 m/s, which is nearly 24% lower than the field test result. When zero damping is considered in the simulation, similar phenomena (as shown in Fig. 12) can be observed, i.e., the peak velocity is 1.013 m/s from the VHI approach and 0.792 m/s from the SHI approach, compared to the peak value of 0.95 m/s obtained from the field test. It shows that, despite the damping applied, the VHI approach is superior to the SHI approach in simulating blasting wave transmission through jointed rock masses. From the results of the four combinations, one can also see that the simulated waveforms agree preferably well with those registered in the field test, only with minor differences in other attributes. It is noteworthy that two major different consequences can be observed between the damped and zero-damping assumption. First, the amplitude obtained from the damped assumption is smaller than that from the zerodamping when same input approach is adopted [cf.
64
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
Fig. 13. Attenuation of peak particle velocity (PPV) (damping ratio ¼ 0.01).
R PPV ¼ 1:5472 Q1=3
2:4699 ðfor VHI approachÞ; ð15Þ
Fig. 12. Resulting velocity-history at a location 8 m above detonation. (a) Results obtained from VHI approach (zero damping). (b) Results obtained from SHI approach (zero damping).
1.013 (for zero damping and VHI); 0.989 (for damped and VHI); 0.792 (for zero damping and SHI); and 0.723 (for damped and SHI)). Second, fibrillation towards the end is observed when zero damping is considered. It indicates that a broader frequency band of the blasting wave can be generated form the zero-damping case. Fig. 13 compares the attenuation of peak particle velocity (PPV) obtained from the two approaches against the field test data (for the cases with damping ratio of 0.01). The PPV trend lines for both approaches are also plotted. Is can be seen that the PPV trend line obtained from the VHI approach is in better accord with the field test data. We note that literature [19] gives the trend line for measured PPV data as follows: 2:5 R PPV ¼ 1:8 ðfor field test dataÞ: ð13Þ Q1=3 The best-fit PPV trend lines plotted in Fig. 13 are: 2:6426 R PPV ¼ 0:9619 ðfor SHI approachÞ; Q1=3 ð14Þ
in which R is the distance (measured in metre) away from the charge center and Q is the TNT charge weight (measured in kilogram). Once again, by comparing Eqs. (14) and (15) against Eq. (13), we notice that the PPV trend line obtained by the VHI approach is closer to the field test data. The results of the numerical simulations show that the VHI approach is better than the SHI in simulating the wave transmission across joints. It concurs to their performance observed in the one-dimensional wave propagation problem in previous section. The same reason applies: the VHI approach is able to provide sufficient constraint to the incident boundary and thus prevent the joint from opening. In fact in a real blast, due to the expanding gas in the blasting cavity, the displacement-history at the blasting wall is always constrained. This reality can be reproduced if the VHI approach is adopted, but not the SHI approach. In addition, the damping effect observed from the numerical results also demonstrates an interesting dilemma: an artificial damping may yield ÔbeautifulÕ waveforms agreeable with test data but leads to smaller amplitude; whereas zero-damping assumption can reproduce no only reasonable amplitude but also the strong fibrillation towards the end. Against this background, a step-wise damping or a very small damping ratio may be adopted in a blast-wave simulating.
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
4. Discussions
65
approach is derived from its capability of reproducing this constraint.
4.1. Equivalence of VHI and SHI approaches It has been well recognized that the VHI and SHI approaches are equivalent for homogenous medium without any interior discontinuities. The stress and velocity are related by Eq. (9). However, it is often misinterpreted and applied incorrectly to non-homogenous medium such as jointed rock masses. As can be seen from Eq. (1), which is invoked when SHI approach is applied, that the velocity is a function of material properties (m) and the damping characteristics (a). Applying Eq. (1) will result in a velocity-history not equivalent to the actual one at the incident boundary. It was clearly demonstrated in Case I (in Section 3.1) that when the waves transmitted through the homogenous rock bar they exhibited little distortion. But in rock masses having joints, even in the case of perfect elasticity, it is difficult to obtain a SHI equivalent to the VHI. 4.2. Consequence at the incident boundary In VHI approach, the velocity-history imposed at the incident boundary is actually a kind of displacement constraint, and thus restraints the boundary to the prescribed position. As long as the velocity conditions are determined [6,17], it will yield correct boundary movements. However, in the SHI approach, Eq. (1) is invoked to derive the equivalent velocity-history (VH) at the incident boundary. Nevertheless, the derived VH is not really equivalent because the effects of the interior joints are not taken into account. In fact, in jointed rock masses, the forwarding waves are reflected from joint planes near the incident boundary, and then superimpose with the incident waves, resulting in false backward displacement of the incident boundary, which is not constrained in the SHI approach. It can be seen clearly from Fig. 7. 4.3. Consequence at the interior joint planes The false backward displacement at the incident boundary mentioned in the previous section generates a series of domino effects. It may lead to opening of the adjacent joint planes. As a result, wave cannot be transmitted efficiently. On the other hand, in the VHI approach, the incident boundary is restrained, and thus movement of adjacent joint planes are limited to a less extent than that in the SHI approach. The separation at joints is lesser and it leads to better transmission of waves through the jointed rock masses. Therefore, the VHI approach exhibits better simulation quality, as demonstrated in the Case II study (in Section 3.2). In reality, due to the existence of expanding gas in the blasting cavity, the displacement-history at the blasting wall is always constrained. The effectiveness of VHI
5. Conclusions This paper focuses on the modelling of dynamic inputs and their effects on wave propagation and wave attenuation in jointed rock masses. The problem lies in the presence of discontinuities. In homogenous rock masses having no joints, the VHI approach and SHI approach are equivalent. However, their non-equivalence in jointed rock masses has been often overlooked. The consequences are illustrated in two case-studies: one-dimensional P-wave propagation along a rock bar and blasting wave propagation in jointed rock masses induced by an underground explosion. The VHI approach exhibits superior simulation capability, while the SHI approach yields inferior results. As a kind of displacement constraint, the VHI approach can reflect correct boundary behaviour, and hence ensure good wave transmission through joints. In contrast, the SHI approach is unable to provide adequate displacement constraint, and thus leads to false backward movement at the incident boundary and subsequent wider opening of adjacent joints, which diminishes the wave transmissibility.
Acknowledgements The second author thank the leave of absence granted by Institute of Rock and Soil Mechanics of the Chinese Academy of Sciences.
References [1] Goodman RE. Methods of geological engineering in discontinuous rocks. St. Paul: West Publishing; 1976. [2] Ghaboussi J, Wilson EL, Isenberg J. Finite elements for rock joints and interfaces. J Soil Mech Found Div, ASCE 1973;99:833– 48. [3] Goodman RE, Taylor RL, Brekke T. A model for the mechanics of jointed rock. J Soil Mech Found Div, ASCE 1968;94:637–59. [4] Schwer LE, Lindberg HE. Application brief: a finite element slideline approach for calculating tunnel response in jointed rock. Int J Numer Anal Methods Geomech 1992;16:529–40. [5] Beer G. Implementation of combined boundary element-finite element analysis with application in geomechanics. In: Banerjee PK, Watson JO, editors. Developments in boundary element methods. London: Applied Science; 1986. p. 191–226. [6] Chen SG, Zhao J. A study of UDEC modeling for blast wave propagation in jointed rock masses. Int J Rock Mech Min Sci 1998;35:93–9. [7] Cundall PA. A computer model for simulating progressive large scale movements in blocky rock systems. In: Proc. Symp. Int. Soc. Rock Mech., vol. 1, Nancy, France; 1971 [paper II-8]. [8] Itasca Consulting Group. UDEC reference manual, Version 3.0. USA: Minneapolis, MN; 1996.
66
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
[9] Zukas JA, Scheffer DR. Practical aspects of numerical simulation of dynamic events: effects of meshing. Int J Impact Engng 2000;24:925–45. [10] Scheffer DR, Zukas JA. Practical aspects of numerical simulation of dynamic events: material interfaces. Int J Impact Engng 2000;24:821–42. [11] Johnson GR. Artificial viscosity effects for SPH impact computations. Int J Impact Engng 1996;18(5):477–88. [12] Chapman TC, Rose TA, Smith PD. Blast wave simulation using AUTODYN2D: a parametric study. Int J Impact Engng 1995;16(5/6):777–87. [13] Zhao C, Valliappan S. A numerical model for wave scattering problems in infinite media due to P and SV wave incidences. Int J Numer Methods Engng 1992;33:1661–82. [14] Zhao C, Valliappan S. An efficient wave input procedure for infinite media. Commun Numer Methods Engng 1993;9:407– 15.
[15] Zhao C, Valliappan S. Incident P and SV wave scattering effects under different canyon topographic and geological conditions. Int J Numer Anal Methods Geomech 1993;17:73–94. [16] Zhao C, Valliappan S. Seismic wave scattering effects under different canyon topographic and geological conditions. Soil Dyn Earthquake Engng 1993;12:129–43. [17] Lemos JA. Distinct element model for dynamic analysis of jointed rock with application to dam foundations and fault motion. Ph.D. thesis, University of Minnesota, USA; 1987. [18] Chen SG. Discrete element modeling of jointed rock masses under dynamic loading. Ph.D. thesis, Nanyang Technological University, Singapore; 1999. [19] Ma GW, Hao H, Zhou YX. Modeling of wave propagation induced by underground explosion. Comput Geotech 1998;22(3/ 4):283–303. [20] Henrych J. The dynamics of explosion. New York: Elsevier Scientific; 1979.
On modelling of incident boundary for wave propagation in jointed rock masses using discrete element method S.C. Fan *, Y.Y. Jiao, J. Zhao Protective Technology Research Centre, School of Civil and Environmental Engineering, Nanyang Technological University, Nanyang Avenue 639798, Singapore Received 7 January 2003; received in revised form 17 October 2003; accepted 20 November 2003
Abstract Numerical modelling is an effective alternative for studying the mechanical behaviour of rock masses. Very often, in solving the same problem with the same code, different models yield remarkably different results. There are many reasons. Amongst them, the boundary conditions are often overlooked but it has significant influence on the stress wave propagation in jointed rock masses. Usually, in solving dynamic problems, the incident boundaries are described either by stress-history input (SHI) or velocity-history input (VHI). This paper presents the study on the effect of these two approaches using discrete element method (DEM) on wave propagation and wave attenuation in jointed rock masses. Two cases were designed to illustrate the points. The first one shows onedimensional P-wave propagation along a rock bar having an internal cross joint. The second one shows the blasting wave propagating through jointed rock masses. The blasting wave is induced by an underground explosion. The rock bar/masses are discretized into elements and the analyses adopt the DEM, which is available in the code ÔUDECÕ. The results reveal that the VHI approach can yield results in accord with field-test data, while the SHI approach may lead to unreasonable results. The two approaches may yield the same results for homogeneous rock, but definitely not the same for jointed rock masses. In addition to the numerical illustrations, this paper also presents some insights of the reason behind it. 2003 Elsevier Ltd. All rights reserved. Keywords: Discrete element method; Wave propagation; Jointed rock; Boundary effects
1. Introduction In rock engineering problems, handling the discontinuities in rock masses is a challenging task. The discontinuities appear in the form of faults, joints or bedding planes. The presence of the discontinuities in rock masses has great influence on the responses to either static or dynamic loads [1], and renders the numerical simulations more complicated. It demands special treatment for the discontinuities, namely joints. There are many numerical tools available. The popular ones include the finite element method (FEM), boundary element method (BEM), finite difference method (FDM), and discrete element method (DEM), just name a few. In FEM, joints are modelled using a kind of
*
Corresponding author. Tel.: +65-6792-1620; fax: +65-6792-1650. E-mail address: [email protected] (S.C. Fan).
0266-352X/$ - see front matter 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2003.11.002
special ‘‘joint elements’’, which simulate the interface discontinuity [2,3]. In FDM, joints are simulated using slide-lines [4]. Very often, BEM is employed for the far field across the joint while FEM is used for near field, and their boundary interface is specified according to the properties of joint [5]. The above-mentioned treatments are efficient when there are only a few joints. For heavily jointed rock masses, an alternative is to treat it as a continuum having equivalent properties [6]. In fact, the equivalent continuum may behave like a jointed rock masses under some circumstances, particularly under static loads. Generally, they fail to represent the dynamic responses and no equivalent properties could be found. On the other hand, the discrete element method (DEM), which is developed to model discontinuous problems, has been recognized to be a better alternative [7]. In DEM, a rock mass is treated as an assemblage of discrete blocks by fictitious or real joints, of which some are the contact interfaces between distinct bodies.
58
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
Individual block can be defined either as rigid or deformable by specifying the material parameters. At the contact interfaces, constitutive laws governing the force– displacement relationship are specified. The solution scheme is identical to that used by the explicit FDM for continuum analysis. The dynamic behaviour is represented numerically by a temporal stepwise algorithm. By using the explicit dynamic approach, DEM is able to simulate large displacements across the discontinuities without encountering much numerical stability problems. The DEM algorithm was coded and available commercially (code-named ÔUDECÕ) [8]. Among the numerical analysts, saying has it that ‘‘garbage in and garbage out’’. Very often, in solving the same engineering problem with the same code, different models yield remarkably different results. There are many reasons: Putting aside those due to inexperienced users, there are many physical phenomena needed to be addressed. Zukas and Scheffler [9,10] investigated the influence of mesh configuration and material interface on numerical results; Johnson [11] studied the effects of introducing an artificial viscosity when using the method called smooth particle hydrodynamics (SPH); Chapman et al. [12] reported their parametric studies on blast wave simulation. Usually, when solving dynamic problems, there can be three kinds of incident boundary (which is defined as boundary that accepts dynamic load), namely displacement-history input (DHI), stress-history input (SHI), and velocity-history input (VHI). Zhao and Valliappan [13–16] proposed a DHI-based wave input procedure in the finite and infinite coupled model to simulate wave-scattering problems in infinite media. They showed that the DHI approach was equivalent to the then prevailing wave-input procedures and had the merit of higher computational efficiency; Lemos [17] employed DEM to investigate the effect of discontinuity in a jointed rock mass subjected to dynamic load generated from a line source. LemosÕ model has only one single discontinuity, and the effects of the two approaches (VHI and SHI) on the resulting slip along the discontinuity were reported. Effects due to other factors such as damping, in situ stress, and material model, etc. were recorded in other literatures [8,17,18]. Amongst the above-mentioned factors, the incident boundary conditions are often overlooked, but it has significant influence on the stress wave propagation in jointed rock masses. This paper presents the study on the effect of the two popular approaches (VHI and SHI) on wave propagation and wave attenuation in jointed rock masses. Two cases were designed to illustrate the points. The first one shows the one-dimensional P-wave propagation along a rock bar having an internal cross joint. The second one shows the blasting wave propagating through jointed rock masses. The blasting wave is induced by an underground explosion.
2. Modelling of boundary conditions In reality, engineering problems always involve infinite or very large domain. Many numerical techniques have been developed to handle those situations, such as special infinite elements in FEM and infinite domain in BEM. Nevertheless, in solving dynamic problems, be it for a finite or infinite domain, the boundary conditions need to be defined. It could be reflective or non-reflective. When the DEM is employed, the domain of interest needs to be limited to a finite size, and some kinds of artificial boundaries are introduced. The artificial boundaries are usually modelled as non-reflective boundaries. Besides, there are natural boundaries which are reflective such as ground surface, excavation surfaces, construction surfaces, etc. In general, boundary conditions can be classified under three categories: prescribed force/stress (PS) boundary; prescribed displacement/velocity (PV) boundary; or boundary with both prescribed stress and velocity. The reflective or non-reflective boundary falls in one of these categories. What follows are treatments for various boundary condition. 2.1. Incident boundary In the dynamic analysis, there are two prevailing approaches of inputting the dynamic excitation, namely VHI and SHI. Correspondingly, the incident boundaries are PV and PS boundary, respectively. In the simulation using DEM, the temporal incremental equation of motion is ðtþDt=2Þ
u_ i
ðtDt=2Þ
¼ u_ i þ
nX
X o Dt ðtÞ ðtÞ ðtDt=2Þ ; F i þ a Þ Fi signðu_ i m ð1Þ
ðtþDtÞ
ðtDtÞ
in which u_ i and u_ i denote the velocity at gridpoint i at time t þ Dt and t Dt, respectively, Fit is the resultant force at grid-point i at time t, a is the damping coefficient, and m is the nodal mass. The primary variables of Eq. (1) are velocities. Therefore, when SHI approach is adopted, the incident boundary is PS type and its implementation is not straightforward. Firstly, the equivalent force at grid-point has to be evaluated, i.e., Fib ¼ rbij nj Ds;
ð2Þ
where nj is the outward normal vector of the boundary segment/surface and Ds is the length/area of the b boundary segment/surface on which the P stress-history rij is imposed. Then the total force Fi acted on the boundary grid-point i is the sum of all applied forces and those derived from internal or inter-block stresses, i.e., Fi ¼ Fib þ Fic þ Fiz þ Fig ;
ð3Þ
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
in which Fic is the resultant contact forces between blocks and exists only at grid-points along the block boundary; Fiz denotes the contribution from internal stresses in the adjacent zone around the grid-point. Fiz can be calculated through Z Fiz ¼ rij nj ds; ð4Þ C
where rij is the zone stress tensor; nj is the unit outward normal to the contour C [8]; Fig is the body force due to gravity, i.e., Fig ¼ mg gi ;
ð5Þ
where P mg is the lumped mass at the grid-point. The value of Fi is then applied in Eq. (1) to calculate resultant grid-point velocity. On the other hand, when VHI approach is adopted, it becomes very straightforward and it is not necessary to invoke Eq. (1) because the primary variables are velocities. The velocity at the boundary grid-point in each subsequent time step can be updated directly form the input data, i.e., ðtþDt=2Þ
u_ i
ðtþDt=2Þ
¼ U_ i
;
ð6Þ
ðtþDt=2Þ
is the velocity-history input, described as where U_ i a function of time. 2.2. Other boundaries Other boundaries include natural boundaries and artificial boundaries. A natural boundary could be ground surface, excavation surface or other construction surface, which is often a reflective boundary and treated in the same way as the PS boundary mentioned above. For the artificial boundaries, they are often non-reflective so that they can silent the outward propagating wave. In the current DEM simulation, two types of artificial boundaries are employed: viscous (non-reflective) boundaries and free-field boundaries [8]. At the viscous boundaries, independent linear viscous dashpots are put in place in the normal and shear directions to absorb any waves coming at incident angles greater than 30. At the free-field boundaries, special boundary elements representing the extended infinite medium are put in place so that all outgoing waves are completely absorbed. It is worth noting that the non-reflection boundaries are PS boundaries where stresses are prescribed. However, the stresses are usually not known in advance. An approximate technique is to derive from the velocities as follows: rn ¼ 2ðqCp vn Þ;
ð7Þ
rs ¼ 2ðqCs vs Þ;
ð8Þ
59
in which rn and rs are the normal and shear stress components, respectively; q is the mass density; Cp and Cs are, respectively, the speed of P-wave and S-wave propagating through the medium; and vn and vs are the normal and transverse particle velocity, respectively. Note that a factor of two is applied to double the stresses in order to overcome the effect of the viscous boundary.
3. Consequence of VHI and SHI approaches The effect of VHI and SHI approaches adopted in DEM simulation of wave propagation and wave attenuation in jointed rock masses are investigated through two case-studies. The analyses were carried out using the UDEC code. 3.1. One-dimensional P-wave propagation along a rock bar To fend off disturbance arising from other factors, the P-wave propagation through a simple one-dimensional rock bar is computed to illustrate the consequence due to different approaches (VHI and SHI). Two scenarios are considered. One is a perfectly elastic rock bar without any joint inside. The other is the same rock bar but having a transverse joint at the mid-length. The length and height of the rock bar are 10.0 and 1.0 m, respectively. The properties of the parent rock bar are: density q ¼ 2650 kg/m3 , bulk modulus K ¼ 44:1 GPa, and shear modulus G ¼ 33:8 GPa. The properties of joint are: normal stiffness kn ¼ 100 GPa/m, shear stiffness ks ¼ 100 GPa/m, cohesion C ¼ 2 MPa, friction angle / ¼ 35, and no tensile strength is considered. No damping is taken into account (DAMP 0 0). In order to simulate the one-dimensional P-wave propagation, viscous boundaries are specified at the upper, lower, and right boundaries. A train of sinusoidal P-velocity waves for VHI approach and corresponding sinusoidal P-stress waves for SHI approach are applied to the left boundary of the rock bar, respectively. The amplitude of the incident velocity wave is 0.065 m/s; the frequency is 200 Hz; and the whole duration is 0.02 s. For the VHI approach, the input is straightforward. But for the SHI approach, the equivalent stress-history is derived from the velocity-history, or vice versa: rn ¼ qCp vn ;
ð9Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where the speed of P-wave Cp ¼ ðK þ 4G=3Þ=q, and thus the amplitude of the stress wave is 1.0 MPa. The resulting velocity-histories at 10 points along the bar are computed. The computational model is shown in Fig. 1. The five points to the left of the joint are denoted by Ai ði ¼ 1; . . . ; 5Þ and the five points to the right of the joint are denoted by Bi ði ¼ 1; . . . ; 5Þ.
60
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
Fig. 1. One-dimensional P-wave propagating along a rock bar.
Fig. 3. Resulting velocity-histories by VHI approach for scenario 2: having a joint. (a) At five points to the left of joint. (b) At five points to the right of joint.
Fig. 2. Resulting velocity-histories at point A1 and B5 for scenario 1: no joint. (a) Results by SHI approach. (b) Results by VHI approach.
• Scenario 1. In scenario 1, no joint exists. Regardless of the VHI or SHI approach being used, the resultant velocity waveform arriving at the far end (point B5 ) was found similar to the input waveform. Little distortion was observed (see Fig. 2(a) and (b)). • Scenario 2. In Scenario 2, there is a transverse joint at the mid-length. First, VHI approach is adopted and the resulting velocity-histories at the 10 points are traced. Fig. 3(a) shows the velocity-histories at the five points (A1 ; . . . ; A5 ) to the left of the joint and Fig. 3(b) for five points (B1 ; . . . ; B5 ) to the right of the joint. It can be seen from Fig. 3(b) that the negative parts of the waveform almost vanished immediately after travelling across the joint. On the other hand, positive parts of the waveform cross the joint with only half waveforms survived. It is worth noting that the incident waves are reflected back from the joint, the reflected waves superimpose with the incident waves (see Fig. 3(a)), and it results in magnified
waveforms at points A2 ; . . . ; A5 after initial 3/4 of a period. In addition, the farther the point is away from the joint, the weaker the magnification effect is. For easy visualization of the consequence, the velocityhistories at point A1 and B5 are plotted together in Fig. 4. Note should be taken that the waveform at
Fig. 4. Resulting velocity-histories by VHI approach for scenario 2: having a joint.
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
point A1 does not change because Eq. (1) is not invoked at the incident boundary. Subsequently, the SHI approach is employed. Results are different from those by VHI approach. Fig. 5(a) portrays the velocity-histories at the five points to the left of the joint. It can be seen that the waveforms are totally distorted due to superimposition of the reflected waves and the incident waves. In addition, the mean velocities all drift away from the initial zero value. It is also observed that the resulting velocity waves at A1 ; . . . ; A5 are all magnified in the process of propagation. The magnification at point A1 is the least, and velocities at A1 become negative at late time. Fig. 5(b) plots the resulting velocity-histories at the five points to the right of the joint. It can be seen that tensile parts of the waveforms are also filtered off by the joint, whereas the compressive parts of waveforms only have one half waveform survived after travelling through the joint. For easy visualization, Fig. 6 shows the contrast of velocity-histories between point A1 and B5 . Fig. 7 shows the maximum axial displacements obtained from the two approaches. It can be seen that SHI approach leads to relatively larger backward displacements than the VHI approach at locations to the left of the joint. Furthermore, SHI approach yields the maxi-
61
Fig. 6. Resulting velocity-histories by SHI approach for scenario 2: having a joint.
Fig. 7. Maximum axial displacements obtained from the two approaches.
mum backward displacement at the incident boundary, while the VHI approach results in steady zero displacement at the incident boundary. Figs. 8 and 9 plot the stress-histories of rx for locations to the left of the joint for the SHI and VHI approaches, respectively. As can be seen, waveforms in
Fig. 5. Resulting velocity-histories by SHI approach for scenario 2: having a joint. (a) At five points to the left of joint. (b) At five points to the right of joint.
Fig. 8. Resulting rx -histories by SHI approach for scenario 2: having a joint.
62
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
through the joint efficiently. To the contrary, in the SHI approach, the wave reflection from the joint allows the backward rigid-body motion and leads to opening of the joint. As such, most of the energy returns back to the left part of the bar. The mechanism will be discussed in Section 4. 3.2. Blasting wave propagation through jointed rock masses
Fig. 9. Resulting rx -histories by VHI approach for scenario 2: having a joint.
both are distorted beyond the initial 3/4 period, and the stress amplitude decreases when it is further away from the incident boundary. In the SHI approach, the resultant stress wave is identical to the original incident wave, and the waveforms at different locations appear similar to each other. In the VHI approach, the resultant stress waves are magnified and it remains so at the incident boundary even after 3/4 period; and the waveforms appear much more disturbed. These phenomena are reasonable and apprehensible. At the joint, the arrival waves are reflected and then the reflected waves superimposes with the incoming waves. In the SHI approach, the incident boundary cannot be sufficiently constrained so that a backward rigid-body movement becomes possible and it leads to stress relaxation at locations to the left of the joint. On the other hand in the VHI approach, the movement of the incident boundary follows exactly the prescribed input waveform and therefore the unwanted backward rigid-body movements are eliminated. As a result, the disturbance and magnification of the waveforms can be revealed. The two scenarios above clearly demonstrate the effect of the way in defining the incident boundary conditions when joints exist in the fractured rock masses. In responding to different incident boundary conditions, the manner of wave propagation is different due to the wave reflection from the joint, and it results in significantly disturbed wave form. Since DEM is a velocitybased method, the numerical results will be sensitive to the prescribed incident boundary conditions and the behaviour at the joint as well. As a result, the outcome is input-dependent. In the VHI approach in scenario 2, the incident boundary is constrained to the prescribed position with respect to time (in accordance of the input waveform), the backward rigid-body motion of the left part of the bar is not allowed to take place when the wave reflects back from the joint, and thus the joint doss not open up. Consequently, the wave is able to transmit
A small-scale field test of ground blasting was reported in literature [19]. The test results are used as benchmark against the results obtained by the current study using VHI and SHI approaches. The explosion chamber, having the dimensions of 8 4 2 m, is located at 115 m below the ground surface. The effective TNT charge weight is W ¼ 606 kg with a loading density of 10 kg/m3 . A DEM model was constructed accordingly as shown in Fig. 10. The analyses are carried out using the UDEC code. In the DEM model, there are two sets of joints: one set dipping at 45 with a spacing of 2 m and the other dipping at 45 with a spacing of 2 m. The parent rock material is assumed linear elastic, while all joints satisfy the Coulomb-slip constitutive model. The properties of parent rock are: density q ¼ 2650 kg/m3 , bulk modulus K ¼ 41 GPa, and shear modulus G ¼ 31 GPa. The properties of joints are: normal stiffness kn ¼ 100 GPa/ m, the shear stiffness ks ¼ 100 GPa/m, cohesion C ¼ 2 MPa, friction angle / ¼ 25, and no tensile strength is considered. At the artificial boundaries (bottom, left, and right), ÔquietÕ-type viscous boundaries are put in place to absorb the outgoing waves. The ceiling and the
Fig. 10. DEM computational model for the ground blasting test.
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
63
walls around the chamber are the incident boundaries. It should be noted that in dynamic modelling, zero or small damping is usually adopted [8]. To study the effect, two different damping ratios (zero and a relatively small ratio of 0.01) are used in the simulations for both approaches (VHI and SHI). Totally, four (2 2) combinations are considered. The blast load, represented by a triangular overpressure-history on the chamber walls, is obtained by the following formulae [20]: 14:0717 5:5397 0:3572 þ 2 3 R R R 0:00625 6 0:3; þ 4 ½kg=cm2 0:05 6 R R 6:1938 0:3262 DpU ¼ 2 R R 2:1324 6 1; þ 3 ½kg=cm2 0:3 6 R R 0:662 4:05 3:288 6 10; DpU ¼ þ 2 þ 3 ½kg=cm2 1 6 R R R R
DpU ¼
ð10Þ DpUr ¼ 2DpU þ
6DpU2 ; DpU þ 7p0
s þ 0:264R 2 0:129R 3 ffiffiffiffiffi ¼ 103 ð0:107 þ 0:444R p 3 W 4 Þ 0:05 6 R 6 3 ½s=kg1=3 ; þ 0:0335R
ð11Þ
ð12Þ
where Dp/ is the maximum overpressure at the blasting wave front; Dp/r is the reflected wave overpressure at the chamber walls; s (measured in second) is the duration of the overpressure; R (measured in metre) is the distance pffiffiffiffiffi ¼ R= 3 W of a point away from the charge centre; R (unit ¼ m/kg1=3 ) is the scaled distance; W (measured in kg) is the charge weight; and p0 (measured in kg/cm2 ) is the atmospheric pressure. From the test configuration, with the loading non-uniformity taken into account, Eq. (11) yields a maximum overpressure equal to 30.23 MPa, and Eq. (12) yields the duration equal to 2.5 ms. This triangular form of stress-history is used directly in the SHI approach. But for the VHI approach, the equivalent velocity-history needs to be computed. The approximation formula (Eq. (9)) is not adopted here. Instead, a more accurate method is used. The velocityhistory at the incident boundary is obtained through a separate DEM analysis on the same model without joints. It yields an equivalent velocity-history of similar triangular form having a maximum velocity of 2.03 m/s over duration of 2.5 ms. Fig. 11(a) plots the resulting velocity-histories obtained by VHI approach against field test data at 8 m above the detonation (using damping ratio of 0.01), while Fig. 11(b) shows the results obtained from SHI approach (using damping ratio of 0.01). We first focus
Fig. 11. Resulting velocity-history at a location 8 m above detonation. (a) Results obtained from VHI approach (damping ratio ¼ 0.01). (b) Results obtained from SHI approach (damping ratio ¼ 0.01).
on the peak velocity, which is considered a key factor controlling material damage. The VHI approach yields a peak velocity of 0.989 m/s, which is close to the field test result. But the SHI approach yields a peak velocity of 0.723 m/s, which is nearly 24% lower than the field test result. When zero damping is considered in the simulation, similar phenomena (as shown in Fig. 12) can be observed, i.e., the peak velocity is 1.013 m/s from the VHI approach and 0.792 m/s from the SHI approach, compared to the peak value of 0.95 m/s obtained from the field test. It shows that, despite the damping applied, the VHI approach is superior to the SHI approach in simulating blasting wave transmission through jointed rock masses. From the results of the four combinations, one can also see that the simulated waveforms agree preferably well with those registered in the field test, only with minor differences in other attributes. It is noteworthy that two major different consequences can be observed between the damped and zero-damping assumption. First, the amplitude obtained from the damped assumption is smaller than that from the zerodamping when same input approach is adopted [cf.
64
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
Fig. 13. Attenuation of peak particle velocity (PPV) (damping ratio ¼ 0.01).
R PPV ¼ 1:5472 Q1=3
2:4699 ðfor VHI approachÞ; ð15Þ
Fig. 12. Resulting velocity-history at a location 8 m above detonation. (a) Results obtained from VHI approach (zero damping). (b) Results obtained from SHI approach (zero damping).
1.013 (for zero damping and VHI); 0.989 (for damped and VHI); 0.792 (for zero damping and SHI); and 0.723 (for damped and SHI)). Second, fibrillation towards the end is observed when zero damping is considered. It indicates that a broader frequency band of the blasting wave can be generated form the zero-damping case. Fig. 13 compares the attenuation of peak particle velocity (PPV) obtained from the two approaches against the field test data (for the cases with damping ratio of 0.01). The PPV trend lines for both approaches are also plotted. Is can be seen that the PPV trend line obtained from the VHI approach is in better accord with the field test data. We note that literature [19] gives the trend line for measured PPV data as follows: 2:5 R PPV ¼ 1:8 ðfor field test dataÞ: ð13Þ Q1=3 The best-fit PPV trend lines plotted in Fig. 13 are: 2:6426 R PPV ¼ 0:9619 ðfor SHI approachÞ; Q1=3 ð14Þ
in which R is the distance (measured in metre) away from the charge center and Q is the TNT charge weight (measured in kilogram). Once again, by comparing Eqs. (14) and (15) against Eq. (13), we notice that the PPV trend line obtained by the VHI approach is closer to the field test data. The results of the numerical simulations show that the VHI approach is better than the SHI in simulating the wave transmission across joints. It concurs to their performance observed in the one-dimensional wave propagation problem in previous section. The same reason applies: the VHI approach is able to provide sufficient constraint to the incident boundary and thus prevent the joint from opening. In fact in a real blast, due to the expanding gas in the blasting cavity, the displacement-history at the blasting wall is always constrained. This reality can be reproduced if the VHI approach is adopted, but not the SHI approach. In addition, the damping effect observed from the numerical results also demonstrates an interesting dilemma: an artificial damping may yield ÔbeautifulÕ waveforms agreeable with test data but leads to smaller amplitude; whereas zero-damping assumption can reproduce no only reasonable amplitude but also the strong fibrillation towards the end. Against this background, a step-wise damping or a very small damping ratio may be adopted in a blast-wave simulating.
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
4. Discussions
65
approach is derived from its capability of reproducing this constraint.
4.1. Equivalence of VHI and SHI approaches It has been well recognized that the VHI and SHI approaches are equivalent for homogenous medium without any interior discontinuities. The stress and velocity are related by Eq. (9). However, it is often misinterpreted and applied incorrectly to non-homogenous medium such as jointed rock masses. As can be seen from Eq. (1), which is invoked when SHI approach is applied, that the velocity is a function of material properties (m) and the damping characteristics (a). Applying Eq. (1) will result in a velocity-history not equivalent to the actual one at the incident boundary. It was clearly demonstrated in Case I (in Section 3.1) that when the waves transmitted through the homogenous rock bar they exhibited little distortion. But in rock masses having joints, even in the case of perfect elasticity, it is difficult to obtain a SHI equivalent to the VHI. 4.2. Consequence at the incident boundary In VHI approach, the velocity-history imposed at the incident boundary is actually a kind of displacement constraint, and thus restraints the boundary to the prescribed position. As long as the velocity conditions are determined [6,17], it will yield correct boundary movements. However, in the SHI approach, Eq. (1) is invoked to derive the equivalent velocity-history (VH) at the incident boundary. Nevertheless, the derived VH is not really equivalent because the effects of the interior joints are not taken into account. In fact, in jointed rock masses, the forwarding waves are reflected from joint planes near the incident boundary, and then superimpose with the incident waves, resulting in false backward displacement of the incident boundary, which is not constrained in the SHI approach. It can be seen clearly from Fig. 7. 4.3. Consequence at the interior joint planes The false backward displacement at the incident boundary mentioned in the previous section generates a series of domino effects. It may lead to opening of the adjacent joint planes. As a result, wave cannot be transmitted efficiently. On the other hand, in the VHI approach, the incident boundary is restrained, and thus movement of adjacent joint planes are limited to a less extent than that in the SHI approach. The separation at joints is lesser and it leads to better transmission of waves through the jointed rock masses. Therefore, the VHI approach exhibits better simulation quality, as demonstrated in the Case II study (in Section 3.2). In reality, due to the existence of expanding gas in the blasting cavity, the displacement-history at the blasting wall is always constrained. The effectiveness of VHI
5. Conclusions This paper focuses on the modelling of dynamic inputs and their effects on wave propagation and wave attenuation in jointed rock masses. The problem lies in the presence of discontinuities. In homogenous rock masses having no joints, the VHI approach and SHI approach are equivalent. However, their non-equivalence in jointed rock masses has been often overlooked. The consequences are illustrated in two case-studies: one-dimensional P-wave propagation along a rock bar and blasting wave propagation in jointed rock masses induced by an underground explosion. The VHI approach exhibits superior simulation capability, while the SHI approach yields inferior results. As a kind of displacement constraint, the VHI approach can reflect correct boundary behaviour, and hence ensure good wave transmission through joints. In contrast, the SHI approach is unable to provide adequate displacement constraint, and thus leads to false backward movement at the incident boundary and subsequent wider opening of adjacent joints, which diminishes the wave transmissibility.
Acknowledgements The second author thank the leave of absence granted by Institute of Rock and Soil Mechanics of the Chinese Academy of Sciences.
References [1] Goodman RE. Methods of geological engineering in discontinuous rocks. St. Paul: West Publishing; 1976. [2] Ghaboussi J, Wilson EL, Isenberg J. Finite elements for rock joints and interfaces. J Soil Mech Found Div, ASCE 1973;99:833– 48. [3] Goodman RE, Taylor RL, Brekke T. A model for the mechanics of jointed rock. J Soil Mech Found Div, ASCE 1968;94:637–59. [4] Schwer LE, Lindberg HE. Application brief: a finite element slideline approach for calculating tunnel response in jointed rock. Int J Numer Anal Methods Geomech 1992;16:529–40. [5] Beer G. Implementation of combined boundary element-finite element analysis with application in geomechanics. In: Banerjee PK, Watson JO, editors. Developments in boundary element methods. London: Applied Science; 1986. p. 191–226. [6] Chen SG, Zhao J. A study of UDEC modeling for blast wave propagation in jointed rock masses. Int J Rock Mech Min Sci 1998;35:93–9. [7] Cundall PA. A computer model for simulating progressive large scale movements in blocky rock systems. In: Proc. Symp. Int. Soc. Rock Mech., vol. 1, Nancy, France; 1971 [paper II-8]. [8] Itasca Consulting Group. UDEC reference manual, Version 3.0. USA: Minneapolis, MN; 1996.
66
S.C. Fan et al. / Computers and Geotechnics 31 (2004) 57–66
[9] Zukas JA, Scheffer DR. Practical aspects of numerical simulation of dynamic events: effects of meshing. Int J Impact Engng 2000;24:925–45. [10] Scheffer DR, Zukas JA. Practical aspects of numerical simulation of dynamic events: material interfaces. Int J Impact Engng 2000;24:821–42. [11] Johnson GR. Artificial viscosity effects for SPH impact computations. Int J Impact Engng 1996;18(5):477–88. [12] Chapman TC, Rose TA, Smith PD. Blast wave simulation using AUTODYN2D: a parametric study. Int J Impact Engng 1995;16(5/6):777–87. [13] Zhao C, Valliappan S. A numerical model for wave scattering problems in infinite media due to P and SV wave incidences. Int J Numer Methods Engng 1992;33:1661–82. [14] Zhao C, Valliappan S. An efficient wave input procedure for infinite media. Commun Numer Methods Engng 1993;9:407– 15.
[15] Zhao C, Valliappan S. Incident P and SV wave scattering effects under different canyon topographic and geological conditions. Int J Numer Anal Methods Geomech 1993;17:73–94. [16] Zhao C, Valliappan S. Seismic wave scattering effects under different canyon topographic and geological conditions. Soil Dyn Earthquake Engng 1993;12:129–43. [17] Lemos JA. Distinct element model for dynamic analysis of jointed rock with application to dam foundations and fault motion. Ph.D. thesis, University of Minnesota, USA; 1987. [18] Chen SG. Discrete element modeling of jointed rock masses under dynamic loading. Ph.D. thesis, Nanyang Technological University, Singapore; 1999. [19] Ma GW, Hao H, Zhou YX. Modeling of wave propagation induced by underground explosion. Comput Geotech 1998;22(3/ 4):283–303. [20] Henrych J. The dynamics of explosion. New York: Elsevier Scientific; 1979.