Studying the Effect of Burial Depth on the Dynamic Response of Tunnels under the Blast Loading

Studying the Effect of Burial Depth on the Dynamic Response of Tunnels under the Blast Loading

WANG Yuan1, Mahdi Shadab Far2

1Professor in Geotechnical Engineering, Hohai University, Nanjing, China

2Ph.D. Student in Bridge and Tunnel Engineering, Hohai University, Nanjing, China

Abstract

Today, underground structures were constructed and exploited considering conditions of cities expansion. Despite mechanized drilling methods, it is preferable to use blasting and drilling in excavation of tunnels particularly in hard rocks for different reasons such as high flexibility of this method. For this reason, it is important to study dynamic behavior of tunnel and rock under blast loading. In this paper, effect of burial depth on dynamic response of tunnel was studied using a parametric study and attempt has been made to give a simplified formula for estimating radial displacement time-history around tunnel in the desirable depth. Results of this study were not complex and can be easily used for real projects.

Key Words: depth of tunnel, radial displacement, blast loading, dynamic response of tunnel

1. Introduction

   There are two major reasons for tendency of human to use underground spaces:

   • Use of underground space for transportation, storage and transfer of material

   • Extraction of material from mines

     which in both cases drilling of tunnels and access roads are the main parts of operation [9].

     Tunnel drilling method can be classified into two groups considering material of the ground. The first group includes tunnel drilling in soft soil and rocks. For drilling in the grounds, traditional mechanical methods (discontinuous drilling such as drilling with shovel) or mechanized methods (continuous methods such as Tunnel Boring Machines (TBM)) are used. The second group includes tunnel boring in hard rocks in which it is preferred to use drilling and blasting which are regarded as traditional and discontinuous methods [7,16].

     Drilling and blasting method was first applied by

     Kaspar Weindl in 1627 in a silver mine in Slovakia [5].

     This method is suitable for hard rocks such as granite, gneiss, basalt and quartz and for soft rocks such as marl, lime, clay and plaster. Although this method is applied for the rock with different conditions, drilling and blasting method is preferred in the following cases:

   • Very hard rocks

   • Noncircular sections

     In general, parameters affecting blasting can be classified as two groups of controllable parameters and non-controllable parameters. Controllable parameters are those which can be corrected or changed by trial and error like the applicable explosives or the desired blasting methods. Non-controllable parameters are those which blasting expert cannot control such as geomechanical characteristics of the rock which should be blasted.

     Selection of the quantity and type of the suitable explosive for blasting is one of the important aspects in blasting design. Some main parameters relating to explosives include velocity of detonation (VOD), density, impedance characteristics and type of explosive. Parameters relating to geomechanical characteristics of rock mass and geology of the region have considerable effect on suitable crushing, vibration control, fly rock and safety aspects. The most important characteristics of rock which are effective on result of blasting include:

   • Tensile and compressive strength of rock: Most rocks have compressive strength of 8 to 10 times as much as tensile strength. This problem is an important factor in rock blasting.

   • Specific weight: The rocks with high specific weight are harder than the rocks with lower specific weight for blasting.

   • Speed of seismic waves: Speed of seismic waves of all kinds of rock varies between 1500 and 6000 m/s. Hard rocks with high speed of seismic waves are easily detonated particularly when the explosives with high detonation rate are used.

   • Rock mass structure: One of the design stages includes study of rock structure and its other characteristics. Rock structure includes joints and other weak parts. In case discontinuities are vertical to tunnel axis, the progression will be good but if discontinuities are parallel to tunnel axis, the progression will be weak and tunnel breast work form will be uneven [14].

     Of other parameters which can be very effective on dynamic behavior of rock during detonation are burial depth and tunnel installation depth. Changes in depth of tunnel considerably affect vertical stiffness of rock [11]. For this reason, it will be effective on displacement around the tunnel. In this paper, a parametric study investigates effect of burial depth parameter on dynamic response of tunnel and at the end a simplified method is presented for estimating response of radial displacement of tunnel crest.

2. Review of literature

   Considering growing application of detonation operation in mineral and civil industries, it is very important to control results of detonation operation such as ground vibration, back break, undesirable displacement, formation and expansion of crashed areola and generally dynamic behavior of the environment around the tunnel. For this reason, different researchers tried to study such subjects since early 50s by supplying mathematical, experimental and analytical models.

   In this regard, Ash (1963) suggested an experimental model to control undesired factors and improve results of detonation. In these models, radius of the damaged zone is calculated as a coefficient of borehole radius using relative density of the rock mass and relative strength of the explosive [1]. N. M. Gay et al. (1971) in a report entitled destruction of model earthen tunnels by internal explosive detonation studied failure mechanism in different models of tunnel considerable behavioral characteristics for the surrounding environment and different geometrical characteristics for tunnel. They also studied effect of tunnel cover depth (burial depth) on failure mechanism of tunnel [4]. Djordjevic (1999) presented a two-component model of cracking around the detonated zone given compressive- shear failure and tensile failure [3]. S. Y. Rohani et al. (2010) modeled tunnel in discontinuous rock using discrete element method and studied effect of detonation and discontinuities inclination angle on optimal depth of tunnel [14]. Anirban De et al. (2010) studied effect of detonation on underground tunnels using some physical centrifuge models. They studied wave-tunnel-structure interaction by making some

   prototypes in different burial depths and different soils. Results of their work can be used for designing underground structures, improving the available structures and also for calibrating numerical methods [2]. Yang et al. (2010) evaluated shallow tunnels of metro affected by detonation. They studied stability and safety of tunnels lining using explicit dynamic nonlinear finite element method. Results of their work show that the upper part of tunnel lining is more susceptible than other parts. They conclude that tunnels with depths of higher than 7 m are safe for detonations on the ground surface with an explosive containing less than 500 g of TNT [17]. Rezvanifar et al. (2012) modeled and analyzed a tunnel structure in weak rock mass under bomb explosion load containing 10 tons of TNT using Finite Element Method (FEM). They studied pressure waves velocity caused by detonation and its effects on stability of tunnel in different depths and evaluated optimal depth of tunel to be used for passive defense structures [12]. Panje et al. (2013) studied behavior of shallow tunnels using Boundary Element Method based on full-plane elastostatics fundamental solutions given linear behavior of rock mass under both gravity loads and tractions and studied effect of depth on stresses caused by surface charge [10].

   In this paper, after mentioning main hypotheses of modeling and analysis, attempt has been made to study the behavior of rock around tunnel during detonation. After performing 6 full dynamic analysis on the tunnels with different depths, effect of burial depth on dynamic response of tunnel was evaluated using a parametric study and then a proposed method was presented for estimating radial displacement time-history above crest of the tunnel.

3. Problem solving methodology

   1. Analysis method

      Two common methods of dynamic analysis include equivalent static method and full dynamic method. In equivalent static method which was introduced by Idriss and Seed, linear analysis is first done by assigning initial values to damping ratio and shear modulus in different zones of the model. Then, maximum shear strain cycle is recorded for each element and damping value and shear modulus relating to the next stage are determined using curves extracted from laboratory tests. New values of damping ratio and shear modulus are used again in another numerical analysis. This process continues until the values of numerical analysis and laboratory tests match. At the end, values of computer simulation can be confirmed. In the full dynamic method, modeling analysis process

      is done only in one phase provided that suitable relation

      2 2

      is used in modeling and damping and shear modulus

      t ( 1

      )

      (4)

      are automatically modeled on the strain surface considering that relation [15].

      Although each one of the two equivalent static and full dynamic methods has its own weaknesses and strengths, full dynamic method has been used in this paper due to intensive vibrations caused by blast

      max

      Where is corresponding angular speed with the maximum vibration frequency of the system and is a fraction of critical damping. Values of these parameters are calculated as follows:

      2

      loading. Main characteristics of full dynamic method can be summarized as follows:

      max

      td

      (5)

      • Using this method, different frequencies are naturally simulated in the model.

      • Irreversible displacements and other permanent changes are automatically modeled.

      • In this method, increase of plastic strains can be modeled according to plasticity theory, if

        0.45

        td

        min

        min

        (6)

        (7)

        necessary.

      • Shear and pressure waves are applied in a simulation, therefore, materials respond to their common and complementary effects which are very important in intensive vibrations.

   2. Problem formulation

      For dynamic modeling of the problem, finite difference formulation was used in the FLAC software. In this regard, stiffness of each studied zone is calculated as follows [18]:

      Where min is damping ratio and min is angular frequency for Reyleigh damping state.

   3. Important points in dynamical modeling

      In dynamical modeling, three important issues were considered which are studied in this section.

      • Dynamic loading and boundary conditions of the model

      • Mechanical damping of the model

      • Conditions for wave transmission inside the model

        kz (K

   1. L2

   G) d T

   3 Az

   (1)

   1. Blast loading

      Dynamic loading in the software can be applied on

      Where K is Young's modulus, G is shear modulus, Ld is length of element diameter, Az is area of element and T is thickness of plane which is considered equal to 1 for plane strain. Solution time steps are determined with the following relation.

      A

      boundaries of the model or internal node points. In this paper, since source of blast wave is the tunnel, the second method was used. Dynamic loading in the software can be applied on the model in one of the four forms of acceleration time-history, velocity time- history, stress or pressure time-history and force time-

      cr x

      cr x

      t min

      Cp max

      (2)

      history. In this paper, stress time-history was used to model the dynamic blast load.

      Where Cp is velocity of pressure wave (p-wave), A is area of the triangular zone and xmax is maximum length of the studied zone. In case there is no damping, a Factor of Safety was also used to estimate critical time interval. Thus, solution time step is proposed as follows:

      t

      Blast loading is a kind of impact load which is propagated in soil and rock around the tunnel. Large part of blasting energy is dissipated during crushing process of the environment around charged bore hole and only small part passes the ground mass as elastic waves and reaches the surrounding environment. Vibrations transfer resulting from blast waves is

      td cr

      2

      (3)

      affected by different factors such as geological structure of the environment, mechanical characteristics

      In cases which damping is also used in equations, time interval steps are reduced and its value is calculated as follows according to Belytschko relation [18].

      of rock and type and rate of the explosive. For this reason, accurate modeling of blasting load is very complex. In this paper, the relation introduced by International Society of Explosive Engineers has been used for determining input compressive blast load (Equation 8).

      d e

      d e

      P 2.5D2 106 (kbar)

      3

      (8)

      for analyzing the model. To assign damping in the studied materials, two kinds of damping can be used:

      Where is explosive density (gr/cm ) and De is

      detonation velocity (m/s).

      Dynamic pressure caused by detonation should be defined as a time function. Blast loading time history is obtained using decay function introduced by Stareld

      • Rayleigh Damping

      • Local Damping

      In the first method, a damping matrix like C is related to component of mass matrix M and stiffness matrix K using coefficients and (Equation 11).

      and Pugliese (1968). Blast compressive load is

      obtained as follows:

      C M K

      (11)

      D B e

      D B e

      P 4P xp 1 Bt exp( 2Bt) (9)

      2

      Where PB is decoupled detonation pressure which is calculated by relation 10, B is experimental load coefficient (=16.338) and t represents time.

      Where coefficients and are mass and stiffness proportional coefficients, respectively.

      The following figure shows damping ratio against frequency. Curve (a) represents that only mass damping component has been considered and curve (b) represents that only stiffness damping component has

      been considered. Curve (c) shows general state and

      P P ( dc )3

      (10)

      outcome of these two.

      d

      d

      B d

      h

      Where dc is charge diameter and dh is borehole diameter.

      Regarding the boundary conditions of the model, it should be noted that restrained boundaries cause reflection of vibrational waves in dynamic analyses and energy dissipation doesnt occur really. In other words, in case the restrained boundaries are not corrected in dynamic state, they are reflected inside the environment when waves reach these boundaries and make the problem face error (in real estate, the mentioned waves pass boundary of the system).

      Use of larger models can compensate for this problem because damping of material in large geometries will cause suitable energy dissipation and waves reflection will be negligible enough after colliding with the restrained boundaries. Hwever, this method increases volume of calculations and consequently increases analysis time.

      Better method for preventing undesirable refection of waves inside the model is to use quiet boundaries and energy absorbent. For this purpose, independent

      dampers are used in vertical and horizontal directions.

      Figure 1 Damping ratio against Frequency

      As observed above, mass damping effect is higher in low frequencies while stiffness damping effect is higher in high frequencies. General curve which includes combination of two mentioned components reaches minimum point with the following value.

      (12)

      min

      When quiet boundaries are used, special attention

      min

      (13)

      should be paid to condition of static loadings and make the mode reach static stability if possible.

      On the other hand, central frequency in which total mass and stiffness damping will be equal to half of total damping is defined as follows:

   2. Damping

      Damping has lower effect on control of structures

      responses to impact loads than periodic and harmonic

      f min min 2

      (14)

      loads because response of structure to an impact load is formed shortly before damping forces can absorb high energy from the structure. For this reason, un-damped response is used for blasting analyses in many cases. In this paper, although damping effect does not seem considerable, damped full dynamic analysis was used

      In this paper, 5% of Rayleigh Damping was introduced for the software using parameters fmin and min with command SET dy_damp rayleigh f.

      Local Damping was first used for quasi-static analyses but it was later used in dynamic analyses due to its specific characteristics. This damping acts by adding or deleting one mass to and from a node point in

      a fluctuation period. When mark of speed in fluctuation changes, mass is added to the model and when the speed reaches minimum or maximum point, it will be deleted. For this reason, increase of kinetic energy is deleted twice in each fluctuation period. In case D is a part of critical damping, local damping coefficient will be equal to:

      To do dynamic analysis of blasting, a tunnel with horseshoe section was modeled in depth of 26 m under the ground surface in a rock. Geometrical specifications of the model are shown in Figure 2 and Table 1.

      Ground Surface

      L D

      (15)

      Because vibration frequency is not considered in calculation of local damping coefficient, it is easier to use local damping coefficient than Rayleigh Damping. However, local damping was not used due to complexity of wave propagation caused by blast in the model.

   3. Wave transmission conditions

   During the dynamic analysis, when wave propagates

   Point (a)

   W

   H

   Center Line

   L

   in soil environment, numerical turbulence may occur under undesired conditions. Frequency applied on the environment and wave velocity are effective on numerical accuracy of wave transmission conditions.

   Figure 2 Schematic plan of tunnel and the surrounding rock

   Table 1 Geometrical specifications of tunnel

   For this purpose, it is necessary to satisfy condition

   Depth of tunnel

   Width of tunnel

   Height of tunnel

   recommended by Kuhlemeyer and Lysmer [13]. These two researchers suggested that size of the zones in the studied continuous medium is smaller than value of the

   H (m) W (m) L (m) 26 11 9

   following relation:

   l

   100

   (16)

   The studied rock was modeled given linear elastic behavior with density of 2 kg/m3, bulk modulus of 54e3 N/m2 and shear modulus of 45e3 N/m2 and in

   Where l is size of the largest element and is wavelength when the highest frequency occurs in the medium. On the other hand, in case C is wave speed of p or s wave in the medium and T is wave period, will be equal to:

   dimensions of 70 m*100 m. Blast loading was also applied on the model according to the method described in Section 3.3.1 and given the characteristics mentioned in Table 2.

   Table 2 Specifications of blast loading

   CT

   (17)

   Parameter Amount

4. Problem solving and study of results

   Before dynamic analysis of the model, there should be static equilibrium in it. For this reason, the model was first statically analyzed and the following cases

   raised time (s) duration (s)

   specific gravity of emulsion blast explosion velocity (m/s)

   0.00006

   0.0006

   1.2

   5900 m/s

   were studied for dynamic analysis after ensuring equilibrium of the system.

   • It was assured that model meshing and dimensions of the available zones satisfy wave transmission conditions.

   • Desired mechanical damping was defined considering material and frequency of the material.

   • Boundary conditions and dynamic loading were applied on the model.

   • Dynamic analysis was done using history command.

   After analyzing the model, in the post-processing stage, the field output of dynamic displacement to static displacement ratio which is a dimensionless parameter and is called dimensionless displacement in this paper was calculated and the dimensionless displacement contour lines were drawn (Figure 3). As Figure 3 shows, dynamic displacement to static displacement ratio increases by ascending order when it reaches the tunnel.

   0.00 10.00 20.00 30.00 40.00 50.00 60.00 70.00 80.00 90.00 100.00 meter

   0.00

   -10.00

   -20.00

   -30.00

   3.0

   4.0

   5.0

   2.0

   1.0

   -40.00

   -50.00

   -60.00

   -70.00

   5.0

   4.0

   3.0

   2.0

   1.0

   meter

   Contour Lines of the Dimensionless Displacement

   Figure 3 Contour lines of the dimensionless displacement in the surrounding rock of the tunnel caused by blast loading

   To study the effect of depth on dynamic response of tunnel, we considered point (a) on top of the tunnel crest (Figure 2) as reference point. The reason for selection of point (a) is that radial displacement of this node is highly affected by changes of the ground surface. Then, five other dynamic analyses were done on the tunnel by selecting 6 different burial depths (Table 3).

   Table 3 Six different depths of tunnel

   H (m)

   Depth 1 6

   Depth 2 16

   Depth 3 26

   Depth 4 36

   Depth 5 46

   Depth 6 56

   Time histories of radial dynamic displacement to static displacement ratio on node (a) for all six analyses are shown in Figure 4.

   Figure 4 Radial dynamic displacement to static displacement ratio on node (a) for different depths of tunnel

   With reducing the tunnel depth, major change is made in vertical stiffness of the environment. For this reason, this change caused maximum dimensionless radial displacement of node (a) in depth of 6 m to be 4 times as much as the corresponding value in depth of 56 m. To study how this trend is performed and to present a stepwise method for estimating displacement in desired depth, we conducted a parametric study.

   We first should have calculated the governing equation of dimensionless radial displacement time- history curves obtained from 6 analyses. But due to complexity of the curves, it was not possible to fit a

   function. For this reason, we classified all time-history curves into two groups: one group includes curve peak section (from the start of the curve to t=0.012 s) and another group includes descending section of the curve

   an exponential function (y=cedx) with their second part (Figures 5.a to 5.i). In this regard, equation of each time-history curve turned into the following equation:

   y ax2 bx x 0.12

   (from t=0.012 s to end of the curve). Then, we fitted a

   dx

   (17)

   quadratic polynomial function with the intercept of 0 (y=ax2+bx ) with the first part of each curve and fitted

   y ce x 0.12

   (a) (b)

   (c) (d)

   (e) (f)

   (g) (h)

   (i) (j)

   (k) (l)

   Figure 5 Dimensionless radial displacement time-history fitted curves for 6 different depths

   Equations of the dimensionless radial displacement time-history curves for different depths of the tunnel are summarized in Table 4.

   Table 4 Governing equations for dimensionless radial displacement time-history curves in 6 different depths

   Depth 1 y = -86253×2 + 1162.2x y = 3.1393e-1.231x	Depth 1	6	-86253 1162.2 3.1393 -1.231
   Depth 2 y = -59897×2 + 853.16x y = 2.4899e0.1682x	Depth 2	16	-59897 853.16 2.4899 0.1682
   Depth 3 y = -53571×2 + 742.1x y = 1.9751e0.2184x	Depth 3	26	-53571 742.1 1.9751 0.2184
   Depth 4 y = -31991×2 + 475.58x y = 1.5184e-0.506x	Depth 4	36	-31991 475.58 1.5184 -0.506
   Depth 5 y = -25734×2 + 375.85x y = 1.1788e-0.169x	Depth 5	46	-25734 375.85 1.1788 -0.169

   Depth 1 y = -86253×2 + 1162.2x y = 3.1393e-1.231x	Depth 1	6	-86253 1162.2 3.1393 -1.231
   Depth 2 y = -59897×2 + 853.16x y = 2.4899e0.1682x	Depth 2	16	-59897 853.16 2.4899 0.1682
   Depth 3 y = -53571×2 + 742.1x y = 1.9751e0.2184x	Depth 3	26	-53571 742.1 1.9751 0.2184
   Depth 4 y = -31991×2 + 475.58x y = 1.5184e-0.506x	Depth 4	36	-31991 475.58 1.5184 -0.506
   Depth 5 y = -25734×2 + 375.85x y = 1.1788e-0.169x	Depth 5	46	-25734 375.85 1.1788 -0.169

   Part 1 Part 2

   As Table 5 shows, we ordered coefficients of the fitted equation against corresponding depth of tunnel.

   Table 5 Coefficients of the dimensionless radial displacement time-history curves equations vs. corresponding depth of tunnel

   H (m) a b c d

   Depth 6

   y = -27326×2 + 372.08x

   y = 1.0122e-0.611x

   Depth 6 56 -27326 372.08 1.0122 -0.611

   Then, using the least square method, we calculated relations between coefficients a to d and depth of tunnel (Figures 6.4 to 6.d). By this way, the following equations were obtained.

   was dynamically analyzed in 6 different depths under blast loading. Then, time-history curve of the node above tunnel crest was extracted in all case loads. The equation of each time-history curve was determined as

   a function comprising of a quadratic polynomial

   a 28242ln(h) 138118 b 375.1ln(h) 1867.4 c 0.98ln(h) 5.0318

   (19.a)

   (19.b)

   (19.c)

   function and an exponential function. Then, relation between coefficients of these equations and depth of tunnel was calculated using sum of least squares method and presented as 4 formulas.

   d 8105 p 0.0086p 0.2772h 2.5331

   (19.d)

   Using the above equations, coefficients a to d can be calculated only with depth of tunnel and then dimensionless radial displacement time-history curves can be predicted using Equation 18. The predicted curve relates to the point above tunnel crest which is highly affected by change of tunnel depth.

5. Conclusion

   In this paper, we studied effect of depth on dynamic response of tunnel in rock after describing factors affecting dynamic behavior of tunnel affected by blast load. For this purpose, a tunnel with horseshoe section

   To use the formulas presented in this paper, depth of tunnel should be substituted in the equations and coefficients a to d should be calculated. Then, dimensionless radial displacement time-history can be estimated above crest of tunnel by substituting coefficients calculated in Relation 18. To use the method introduced in this paper, it should be noted that the proposed equations are applied for shallow and semi-deep tunnels and using them for deep tunnels (above 70 m) has high error and is not provable. On the other hand, the proposed equations were extracted according to the problem hypotheses for tunnel in sound rock. In case tunnel is available in jointed rocks, use of the presented formulas may not give accurate results.

   (a) (b)

   (c) (d)

   Figure 6 Relations between coefficients a to d and depth of tunnel

   Effect of mechanical specifications of rock, the presence of joint with geometry and different behavioral characteristics, geometrical shape of tunnel and finally type and rate of explosive on dynamic behavior of tunnel affected by blast load can be studied using the general method proposed in this paper as new research subjects in future.

6. References

1. Ash, R. L. The Mechanics of Rock Breakage, Standards for Blasting Design. Vol. 56. Pit and Quarry. Print.

2. De, Anirban, Thomas F. Zimmie, Tarek Abdoun, and Anthony Tessari. "Physical Modeling of Explosive Effects on Tunnels." Proc. of Fourth International Symposium on Tunnel Safety and Security, Germany, Frankfurt Am Main. 2010. 159-

   68. Print.

3. Djordjevic, N. "A Two-component Model of Blast Fragmentation." Proceedings of Sixth International Symposium on Rock Fragmentation by Blasting. Johannesburg: SAIMM, 1999. 213-19. Print. Ser. 21.

4. Gay, Wallace M. Destruction of Model Earthen Tunnels by Internal Explosive Detonation. Vicksburg, MS: Waterways Experiment Station, 1971. Print.

5. Hollister, Graham. "The Introduction of Powder." Bulletin of the Peak District Mines Historical Society 12.3 (1994). Print.

6. Johansen, John, and C. F. Mathiesen. Modern Trends in Tunnelling and Blast Design. Rotterdam, Netherlands: A.A. Balkema, 2000. Print.

7. Kolymbas, Dimitrios. Tunelling and Tunnel Mechanics: A Rational Approach to Tunnelling. Berlin: Springer, 2005. Print.

8. Liu, Chih-Hsing, and Kok-Meng Lee. "Dynamic Modeling of Damping Effects in Highly Damped Compliant Fingers for Applications Involving Contacts." Journal of Dynamic Systems, Measurement, and Control 134.011005 (2012): 1-

   9. Print.

9. Olofsson, Stig O. Applied Explosives Technology for Construction and Mining. Arla, Sweden: APPLEX, 1990. Print.

10. Panje, Mahdi, Jafar A. Marnani, Hamid A. Elahi, Hamid Kohsari, and Mohammad A. Pira. "Evaluation of Effective Parameters on Stress Behavior of 2D Shallow Tunnels Using Boundary Elements Method." Transportation Research Journal (TRJ) 10.1 (2013): 17-28. Web.

11. Rahimian, Mohammad, Babak Omidvar, Benyamin Kiomarsi, and Alireza Sanayi Ha. "Dynamic Response of Circular Tunnels without Lining under the Influence of Internal Explosion."

    Journal of Civil Engineering and Surveying 44th ser. 1 (2009): 51-60. Print.

12. Rezvanifar, Sajad, Hamed Pourbaferani, and Naser Madani. "Analyzing Tunnel Stability under Explosive Loads Using Dynamic Modeling, Case Study: Basij Twin Tunnels." Proc. of Fourth Iranian Mining Engineering Conference, Tehran University, Tehran, Iran. 2012. Print.

13. RL, Kuhlemeyer, and J. Lysmer. "Finite Element Method Accuracy for Wave Propagation Problems." Soil Mechanics & Foundations 99 (1973): 421-27. Print.

14. Rohani, Seyed Yahya, Fereidoon Khosravi, and Gholam Reza Cheraghi Aliabadi. "Effects of Explosion and Discontinuity Dip Angles on the Optimum Depth of Tunnel." Passive Defense Quarterly 5th ser. 1.1 (2011). Web.

15. Seed, H. Bolton, Robert T. Wong, I. M. Idriss, and K. Tokimatsu. "Closure to Moduli and Damping Factors for Dynamic Analyses of Cohesionless Soils by H. Bolto Seed, Robert T. Wong, I. M. Idriss, and K. Tokimatsu (November, 1986, Vol. 112, No. 11)." Journal of Geotechnical Engineering 114.8 (1988): 958. Print.

16. Sinha, R. S. Underground Structures: Design and Instrumentation. Amsterdam: Elsevier, 1989. Print.

17. Yang, Yubing, Xiongyao Xie, and Rulu Wang. "Numerical Simulation of Dynamic Response of Operating Metro Tunnel Induced by Ground Explosion." Journal of Rock Mechanics and Geotechnical Engineering 4th ser. 2 (2010): 373-

    84. Print.

18. Zamiran, Siavash. Modeling and Analysis of Soil and Rock Structures in FLAC Software. 1st ed. Vol. 1. Tehran: Noavar, 2012. Print. ISBN: 978- 600-168-055-7