Applied Science, Faculty of
Mining Engineering, Keevil Institute of
DSpace
UBCV
Mitelman, Amichai
Master of Applied Science - MASc
University of British Columbia
This thesis presents the application of the finite-discrete element method for simulation of the impact of external blast loads on tunnels in rock. An extensive database of field tests of underground explosions above tunnels is used for calibrating and validating the proposed numerical method. The numerical results are shown to be in good agreement with published data for large-scale physical experiments. 1D and 2D model results are compared to analytical spalling equations and to the field test findings. It was found that only the 2D models are suitable for support design. The influence of rock strength on tunnel durability to withstand blast loads is investigated. It was found that higher rock strength will increase the tunnel resistance to the load on one hand, but decrease attenuation on the other hand. Thus, under certain conditions, results for weak and strong rock masses are similar. Finally, a discussion on tunnel support design to withstand blasting is presented. A distinction between heavy spalling and light rockfall is made based on an estimation of the ratio of peak stress of the arriving wave to the rock tensile strength. Accordingly, different design approaches are suggested: for heavy spalling a low impedance isolating layer between the tunnel liner and surrounding rock is recommended. For light rockfall, a simplified static FEM analysis procedure is presented.
NUMERICAL ANALYSIS OF THE EFFECTS OF EXTETNAL BLASTS ON TUNNELS by Amichai Mitelman A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Mining Engineering) THE UNIVERISTY OF BRITISH COLUMBIA (Vancouver) June, 2015 © Amichai Mitelman, 2015 ii Abstract This thesis presents the application of the finite-discrete element method for simulation of the impact of external blast loads on tunnels in rock. An extensive database of field tests of underground explosions above tunnels is used for calibrating and validating the proposed numerical method. The numerical results are shown to be in good agreement with published data for large-scale physical experiments. 1D and 2D model results are compared to analytical spalling equations and to the field test findings. It was found that only the 2D models are suitable for support design. The influence of rock strength on tunnel durability to withstand blast loads is investigated. It was found that higher rock strength will increase the tunnel resistance to the load on one hand, but decrease attenuation on the other hand. Thus, under certain conditions, results for weak and strong rock masses are similar. Finally, a discussion on tunnel support design to withstand blasting is presented. A distinction between heavy spalling and light rockfall is made based on an estimation of the ratio of peak stress of the arriving wave to the rock tensile strength. Accordingly, different design approaches are suggested: for heavy spalling a low impedance isolating layer between the tunnel liner and surrounding rock is recommended. For light rockfall, a simplified static FEM analysis procedure is presented. iii Preface This thesis is original and independent work done by the author. The author was the lead author in the two published manuscripts based on the work for this thesis and the co-author was the thesis supervisor. The paper “Modelling of blast-induced damage in tunnels using a hybrid finite-discrete numerical approach”, published in the Journal of Rock Mechanics and Geotechnical Engineering (2014) is based on Chapters 3, 7 and 8 in this thesis. The paper “Analysis of Tunnel Support Design to Withstand Spalling Induced by Blasting,” submitted Feb. 2015 and currently in the revision process, is based on Chapters 5, 6 and 9. iv Table of Contents Abstract ...................................................................................................................................... ii Preface ....................................................................................................................................... iii Table of Contents ...................................................................................................................... iv List of Tables............................................................................................................................. vi List of Figures .......................................................................................................................... vii List of Abbreviations................................................................................................................. ix Acknowledgements .................................................................................................................... x Dedication ................................................................................................................................. xi 1 Introduction ......................................................................................................................... 1 1.1 Problem Statement ...................................................................................................... 1 1.2 Research Objectives .................................................................................................... 1 1.3 Thesis Organization ..................................................................................................... 2 2 Literature Review ............................................................................................................... 3 2.1 Tunnels under Blast Loads .......................................................................................... 3 2.2 The ERA Underground Explosion Tests ................................................................... 10 3 General Methodology and Program Overview ................................................................. 14 3.1 General Methodology ................................................................................................ 14 3.2 Numerical Approach ................................................................................................. 16 4 Assessment and Characterization of Blast Load............................................................... 20 4.1 Methodology ............................................................................................................. 20 4.2 Model Setup .............................................................................................................. 22 4.3 Analysis of Results and Discussion ........................................................................... 27 5 One-Dimensional Analysis of Blast Load Using ELFEN ................................................ 31 5.1 Introduction ............................................................................................................... 31 5.2 Model Setup .............................................................................................................. 34 5.3 Analysis of Results and Discussion ........................................................................... 36 6 Two-Dimensional Analysis .............................................................................................. 39 6.1 Methodology ............................................................................................................. 39 6.2 Model Setup .............................................................................................................. 40 6.3 Analysis of Results and Discussion ........................................................................... 41 v 6.4 Influence of Joints on Spalling .................................................................................. 46 6.4.1 Methodology and Model Setup .......................................................................... 46 6.4.2 Analysis of Results and Discussion of Model with Joints ................................. 48 6.4.3 Analysis of Results and Discussion of Model with Material Heterogeneity ..... 50 7 Modeling of the ERA Tests using a FEM/DEM Approach .............................................. 52 7.1 Methodology ............................................................................................................. 52 7.2 Calibration Model Setup ........................................................................................... 53 7.3 Models Setup for Additional Damage Zones ............................................................ 56 7.4 Analysis of Results and Discussion ........................................................................... 57 8 Comparison of Different Rock Types ............................................................................... 65 8.1 Methodology ............................................................................................................. 65 8.2 Analysis of Results and Discussion ........................................................................... 68 9 Support Design to Withstand Blasting ............................................................................. 71 9.1 Introduction and Methodology .................................................................................. 71 9.2 Support Design for Zone 3 Type Damage ................................................................. 72 9.2.1 Model Setup ....................................................................................................... 74 9.2.2 Analysis of Results and Discussion ................................................................... 75 9.3 Support Design for Zone 4 Type Damage ................................................................. 76 9.3.1 Analysis of Results and Discussion ................................................................... 78 9.4 Proposed Flowchart for Support Selection ................................................................ 80 10 Conclusions and Recommendations ............................................................................. 82 10.1 Research Conclusions ................................................................................................ 82 10.2 Recommendation for Further Studies ........................................................................ 83 References ................................................................................................................................ 85 Appendix 1- Matlab Formulation for Attenuation Influence on 1D Spall ............................... 90 vi List of Tables Table 1 Maximum scaled distance of damage zone (modified from the COE Manual TM-5-857) .......................................................................................................................................... 12 Table 2 Percentage of damaged area measured in the ERA tests (modified from EM 1110-345-432) ................................................................................................................................... 14 Table 3 TNT properties ............................................................................................................ 24 Table 4 Sandstone averaged intact properties from the ERA tests. ......................................... 24 Table 5 Mohr-Coulomb fit to Hoek-Brown failure criterion for the sandstone material......... 25 Table 6 PPV and stress results from AUTODYN model. ........................................................ 27 Table 7 Elastic (dynamic) constants used in AUTODYN model ............................................ 29 Table 8 Properties for the 1D bar models in ELFEN ............................................................... 34 Table 9 Number of spalls in the ELFEN models compared to spall theory ............................ 38 Table 10 Material properties for the sandstone model in ELFEN ........................................... 53 Table 11 Distances of models corresponding to damage zones ............................................... 57 Table 12 Comparison of results for damage zones from the ERA tests and ELFEN models.. 64 Table 13 Properties for the different rock types....................................................................... 67 Table 14 PPV attenuation coefficients for the different rock types ......................................... 70 Table 15 Properties of the fiber reinforced concrete ................................................................ 75 Table 16 Input parameters for static analysis ........................................................................... 78 vii List of Figures Figure 1 Thesis structure, main research chapters ..................................................................... 2 Figure 2 Damage zones defined in the ERA tests .................................................................... 12 Figure 3 Yield surface and softening curve for the Rankine rotating crack in ELFEN (modified from ELFEN user’s manual, Rockfield 2005). ....................................................... 17 Figure 4 Bar model with joint to investigate transmission and reflection through joints in ELFEN ..................................................................................................................................... 19 Figure 5 Procedure for blast load estimation from AUTODYN to ELFEN model ................. 21 Figure 6 Gauge layout for AUTODYN model ........................................................................ 23 Figure 7 3D view of the axisymmetric model in AUTODYN ................................................. 23 Figure 8 Comparison of AUTODYN results with results from ERA tests and Hendron’s tests. .................................................................................................................................................. 28 Figure 9 Stress time-history for Gauge 1. ................................................................................ 28 Figure 10 ELFEN normalized face load time curve. ............................................................... 30 Figure 11 Spalling process due to a compressive pulse reflecting at a free face ..................... 32 Figure 12 1D Model setup and load curve ............................................................................... 35 Figure 13 Mesh for bar models. ............................................................................................... 35 Figure 14 Spalling in Models 1-4............................................................................................. 36 Figure 15 Graphs of analytical equations with ELFEN models for (a) spall thickness to tensile strength and (b) spall velocity to tensile strength..................................................................... 38 Figure 16 Definition of angle of fracturing and fracture thickness .......................................... 40 Figure 17 2D model setup. ....................................................................................................... 41 Figure 18 Results of fracturing for: (a) Model 1, 𝜎𝑡 = 2 𝑀𝑃𝑎 (b) Model 2, 𝜎𝑡 = 4 𝑀𝑃𝑎 (c) Model 3, 𝜎𝑡 = 6 𝑀𝑃𝑎 (d) Model 4, 𝜎𝑡 = 8 𝑀𝑃𝑎 ................................................................... 42 Figure 19 ELFEN Graphs of 2D analysis results for (a) Angle of fracturing vs. tensile strength (b) Fracture thickness vs. tensile strength ................................................................................ 43 Figure 20 Graphs for first spall appearance to tensile strength ................................................ 44 Figure 21 Principal stress trajectories for the wave-tunnel interaction from Model #1 after: (a) 2.11 ms (b) 2.73 ms (c) 3 ms .................................................................................................... 45 Figure 22 Model with explicit joints ........................................................................................ 48 Figure 23 Spalling above tunnel roof in (a) model with joints (b) identical model with no joints ......................................................................................................................................... 49 Figure 24 Comparison of (a) model with heterogeneity and (b) regular model ...................... 50 Figure 25 Wave propagation in (a) model without material heterogeneity (b) model with material heterogeneity .............................................................................................................. 51 Figure 26 Model setup: (a) model constraints (b) model geometry and damping regions ...... 55 viii Figure 27 PPV attenuation results of calibrated ELFEN model in comparison to the mean attenuation established by Ambraseys and Hendron (1968) and from the ERA tests in sandstone. ................................................................................................................................. 56 Figure 28 Fracture patterns from (a) Rock Mechanics for Underground Mining, (modified from Brady and Brown, 2006) (b) ELFEN Model #1 .............................................................. 58 Figure 29 Principal stresses time-history for radial and tangential stresses (in Pa units). ....... 60 Figure 30 Fracturing results for Models 1-3 at t=10 ms corresponding to (a) Zone 2 (b) Zone 3 (c) Zone 4 and (d) Zone with no damage. ................................................................................ 61 Figure 31 Spalling in Zone 3, full model and close-up ............................................................ 62 Figure 32 Fracturing results for the 10 m diameter tunnel in different rock types: (a) limestone (b) sandstone (c) basalt (d) granite ........................................................................................... 68 Figure 33 Stress attenuation vs. scaled distance as measured for the simulated rock types .... 69 Figure 34 Simulated results of liners with 65 cm thickness for a) model with gap between liner and rock b) model with adhesive bond between liner and rock....................................... 75 Figure 35 Phase2 Model for light damage assessment ............................................................ 77 Figure 36 Capacity diagrams for a 15 and 25 cm liner from Phase2 models. ......................... 79 Figure 37 Flowchart for design of support for tunnel under blasting ...................................... 80 ix List of Abbreviations 𝜎𝑡 : Tensile strength 𝜎𝑚 : Amplitude of blast wave 𝑐 : Wave velocity 𝐺𝑓 : Fracture energy E : Young’s Modulus G : Shear modulus PPV : Peak particle velocity 𝑓 : Frequency 𝜌 : Density λ : Wavelength x Acknowledgements I would like to express my sincere gratitude to my supervisor, Dr. Davide Elmo, for his continuous support and guidance, and for sharing his broad knowledge in rock mechanics and vast experience in numerical modeling. I could not have imagined a better supervisor and teacher! I would also like to thank my thesis committee members, Profs. Eberhardt and Scoble, for taking the time to read my thesis and add their invaluable comments. Special thanks to Dr. Melanie Armstrong and the Rockfield support team for their help and useful advice for modeling in ELFEN. Finally, I would like to thank my wife and daughters for joining me on this journey, and for being an everlasting source of joy and inspiration in my life. xi Dedication “Give me shelter.” The Rolling Stones 1 Introduction 1.1 Problem Statement Tunnels can be subjected to external blast loads from nearby blasting such as production blasting in mining, civil works or weapons attack. For the engineer given the task of designing tunnel support to withstand blast loads or to evaluate the performance of an existing tunnel, very few published guidelines can be found. The most extensive tests of damage to tunnels from blasting took place between the years 1948-1952, conducted by Engineering Research Associates together with the US Army Corps of Engineers (COE) and the US Bureau of Mines. In these tests measurements of the amount of broken rock and the velocity of the rock ejection into the tunnel were made so that tunnel support could be designed to withstand these. However, establishing tunnel support from the test findings is not a straight forward process; as the tests serve as a valuable database, there are still many questions that require answers in order to properly address the design problem. At the time of the tests, both the fields of rock mechanics and numerical modeling were, for the most, at their infancy. Taking advantage of both the valuable knowledge gained in rock engineering and the high computational capabilities of today, many new insights can be reached. 1.2 Research Objectives The primary research objectives of this thesis are: 1. To validate the suitability of the FEM/DEM code ELFEN to model tunnels under blast loads. 2. To determine the influence of rock mass properties on tunnel durability to blasting. 3. To develop a method for tunnel support design to withstand blasting. In order to achieve these objectives several secondary objectives are addressed: 1. Establishing the proper failure criterion acceptable for modeling the subject of blast loads impacting tunnels. 2 2. Investigation of a conceptual 1D model in ELFEN and comparison to an analytical 1D solution to evaluate the suitability of ELFEN for simulating the spalling process. 3. Comparison of the above 1D model to 2D models to conclude if the 1D simplification is suitable for the subject in question. 4. Investigation of the spalling evolution in all models so that acceptable values for the problem of the impact on the liner can be established. 1.3 Thesis Organization The thesis is divided into 10 chapters; following this introductory chapter are a chapter summarizing the literature review and a chapter discussing the general methodology undertaken. The final chapter summarizes the conclusions and recommendations. The 6 research chapters in between, follow the structure displayed in Figure 1 below: Figure 1 Thesis structure, main research chapters 3 2 Literature Review 2.1 Tunnels under Blast Loads The vast majority of studies and publications on the subject of support design for tunnels deal with static loads. In addition to the static loads that are inevitably imposed by the in-situ stresses and the redistribution of stresses due to the tunnel excavation, tunnels can also be subjected to dynamic loads resulting from earthquakes, rockbursts and nearby blasting. Analysis of rock dynamics necessitate time marching solution schemes; they are further complicated by the fact that the increased loading rate induces a change in the mechanical behavior of the rock materials and rock masses. There are no existing standard methods for rock dynamic testing; in rock engineering practice, guidelines and design methodologies for dealing with dynamic problems are generally lacking (Zhou and Zhao, 2011). A tunnel can be subjected to blasting from drill and blast excavation in an adjacent tunnel, mine production blasting, accidental explosions or weapons attack. In comparison to surface structures, tunnels are considered highly resistant to dynamic loading. A study by Rozen (1978) on 71 cases of tunnels subjected to earthquake shaking showed that little to no damage occurred in tunnels situated in rock. This can be explained by the nature of tunnels being constrained to the movement of the surrounding ground while surface structures have greater freedom to develop large displacements. Wang (1993) summarizes guidelines for seismic design of tunnels including types of deformations of the tunnel liner. Unlike earthquakes, blasting operations are characterized by high frequencies. The relationship between frequency f and wavelength λ is: 𝝀 =𝒄𝒇 [1] where c is the wave velocity. As a result, the shorter wavelengths apparent in blasting will not cause the same types of deformations as in earthquakes. The typical wavelength resultant from blasting will be smaller than the structure and therefore cannot cause shaking of the structure in a coherent manner (Dowding, 1996). In addition, as attenuation of the wave is related to the number of cycles, waves with higher frequencies will attenuate much more rapidly. Due to these differences, blasting can cause heavy damage to tunnels at very close 4 ranges, whereas only localized types of damage will occur at larger distances. In both cases, damage is in the form of detachment of rocks that fly into the tunnel space with an initial velocity. As the distance of tunnel to blast increases, the intensity of cracking and the velocity of the detached fragments decrease. Due to the complexity of full analysis of rock masses subjected to blasting the traditional practical approach applied in the field of blast design in construction and mining projects is largely based on simplified empirical approaches. For instance, damage thresholds for different structures were established and correlated to measured Peak Particle Velocity (PPV) values that caused cracking damage. These empirical findings have been incorporated into codes and regulations for mining and construction. Of the different variables that could be measured, PPV was found to be the most effective in correlating structural damage in residential structures close to blasting operations (Dowding, 1996). Other codes consist of a combined PPV and frequency damage threshold in which higher frequencies permit higher PPVs (Cameron, Forsyth and Kleine, 2014). PPV values are commonly plotted as a function of scaled distances. For cylindrical shaped charges, the scaled distance is the distance divided by the square root of the explosive weight. For spherically propagating charges, where the length of the explosive is small compared to its other dimensions the scaled distances is the distance divided by the cube root of the explosive weight. The use of scaled distances is efficient, as results are not limited to the single weight of the explosive tested in the field or in the numerical model, and conclusions can readily be made for all explosive weights. In practice the attenuation parameters 𝐾 and 𝑛 are found by fitting curves of PPV attenuation versus the scaled distance using trial blasts conducted at the site according to: 𝑷𝑷𝑽 = 𝑲 ∗ (𝑹𝑾𝒃)−𝒏 [2] where 𝑅 is the radial distance from the blast, 𝑊 is the weight of the explosive, 𝑏 refers to either the cube root or the square root according to the shape of the charge, and 𝐾 and 𝑛 are constants that are dependent on the local geology. 5 Most commonly, PPV measurements are made in the context of avoiding damage to surface structures where the waves of concern are surface waves. In this thesis the subject of interest is damage to underground tunnels, where the waves that reach the tunnel from blasting are body waves. In reference to body waves propagating spherically in rocks, different physical tests have been carried out to establish mean values of 𝐾and 𝑛. For example, Ambraseys and Hendron (1968) conducted a series of tests and found that 𝑛=2.6 and 𝐾=11.45 for scaled distances smaller than 4.14 m/kg13⁄ , whereas for larger scaled distances values of 𝑛=1.6 and 𝐾=2.08 were derived . Henrych (1979) conducted tests in granite and found 𝑛 and 𝐾 to be 0.7 and 2 respectively. Johnson, Rizzo and Rozen (1988) conducted tests in limestone and chalk and obtained results of 𝑛=2.75 and 𝐾=11.43-12. In the ERA tests (ERA, 1953) a large number of underground explosions were detonated above tunnels, and damage zones as a function of scaled distance from the charge were empirically defined. The work carried out in this thesis will rely on results and observations from the ERA tests for numerical model calibration and analysis, therefore a separate section in this thesis is dedicated to summarize the main findings of these tests. It is noted that the ERA tests and previous references cited, deal with explosions that are fully coupled with rock. When the charge hole has greater volume than the explosive volume (i.e. uncoupled state), the PPV generated by the same quantity of explosives for a given distance is substantially reduced (Zhou, 2011). The relation between peak particle velocity and stress in an elastic wave stems from Newton's second law of motion and is: 𝝈𝒎𝒂𝒙 = 𝝆𝒄 ∗ 𝐏𝐏𝐕 [3] where the product 𝜌𝑐 is termed “ acoustic impedance”. Since the acoustic impedance in a rock medium is close to constant, the PPV attenuation is proportional to the stress attenuation. It is noted that unlike vibration of structures where the velocity of the ground is the direct cause of structural damage, in tunnels the peak stress of the impact wave is the main cause of dynamic 6 failure. According to Equation [3] different rock materials will respond differently to equal PPVs according to their acoustic impedance. In the context of damage to tunnels from blasting, once the attenuation has been established for the specific site, the safe distance and/or magnitude of blasts that could cause damage can be determined. From the ERA results it was found that the minimum PPV that causes damage in an unlined tunnel is 460 mm/sec. Langefors and Kihlstrom (1976) obtained similar results in tests they conducted. Kendorski et al. (1973) tested bolted and lined tunnels in heavily jointed schist using ANFO as an explosive and found that hairline cracks formed in the liners under PPVs of 900 mm/sec and that displacement of the cracks was associated with PPVs of 1220 mm/sec. With the PPV value of 50 mm/sec commonly used as a damage threshold for surface structures, the numbers above demonstrate once again the high resistance of tunnels to dynamic loads. Due to the scatter in PPV results there is much room for engineering judgment in the determination of the damage threshold. Jiang and Zhou (2011) determined the safe PPV for a lined tunnel subjected to nearby blasting from an adjacent tunnel to be as low as 110 mm/sec. Chen (2006) concluded that the safety of a tunnel liner is governed by the weak link of the adhesion in the rock to liner interface, with the adhesion most commonly being lower than the tensile strength of the rock. With the emergence of many powerful numerical modeling programs and the advance and availability in computational capabilities, many attempts have been made to harness these for modeling of blasts in rock for both mining and construction purposes. A preliminary issue that must be addressed prior to conducting numerical simulations is to determine what input parameters should be used for numerical analysis of dynamic problems in rock. It is clear from many works that there is some dependency of the loading rate on the rock’s strength parameters. For example, Kubota et al. (2007) used an explosive to load a water filled pipe that impacted on a sandstone specimen. The water pipe length was adjusted in order to measure the response to different strain rates. It was found that the tensile strength of the sandstone increases by a power of one third of the strain rate. Zhao (2000) conducted dynamic tests on granite and found that the Hoek-Brown failure criterion is applicable for 7 dynamic problems. The same author concluded that the dynamic tensile and shear strengths increase with loading rate with semi-log relations and that the parameter 𝑚 is left unaffected by the strain rate. Nevertheless, it can be seen in different publications of numerical studies that static parameters are commonly used for dynamic problems. This is due to the fact that there are currently no existing standard methods for rock dynamic testing and no accepted empirical formulas based on an extensive database that relate static to dynamic strength (Zhou & Zhao, 2011). For the far field impact of blasting which is the focus of this thesis, different methods of numerical modeling successfully predicted outcomes measured in the field. Wei and Zhao (2008) modeled effects of a blast on a rock cavern in massive granite with LS-DYNA (Livemore software, 2006). Ma et al. (1998) conducted a field explosion test in limestone and used AUTODYN (Century Dynamics, 2006) to model the wave propagation and resultant damage. They used a piecewise Drucker-Prager strength criterion and a rate-dependent constitutive relationship to model the energy dissipation caused by damage and permanent plastic deformation. Both LS-DYNA and AUTODYN are continuum based programs, and do not have the capability of explicitly modeling material discontinuities. Discontinuities in rock masses, commonly termed as joints, have a highly significant influence on rock mass strength and generally govern rock engineering; it can be said that the greatest effort in static design in rock is related to the effect of joints. In numerical modeling two main approaches in dealing with rock joints exist, (a) an equivalent continuum approach based on the GSI rating system (Hoek et al. 1995), where a reduction in the overall strength of the rock mass relative to the intact rock specimen are applied on the basis of the joint characteristics, and (b) explicit modeling of the joints where the intact rock and joints are each assigned their own strength properties. Elmo (2006) summarizes and discusses the different implicit and explicit numerical methods that are used to model brittle failure in rock, as well as their advantages and limitations. In the context of dynamic design of tunnels, it is assumed that joints have a significant influence as well. Pyrak-Nolte (1988) derived analytical solutions for P and S waves arriving at different incident angles to a joint. Her model, termed as the displacement-discontinuity model, defines the joint as a non-welded contact represented by two essential boundary 8 conditions: (1) the stress is continuous across both sides of the joint; (2) the displacements across the two sides of the joint are discontinuous, and the difference in displacement between the two sides divided by the stress is equal to the stiffness of the joint. Pyrak-Nolte's model was confirmed in laboratory tests. Different studies have been carried out utilizing numerical programs with explicit joint modeling capabilities to investigate the effects of these on wave propagation. Wang et al. (2005) used 3DEC (Itasca, 2005) and found that transmission and reflection through a single joint and two parallel joints are in good agreement to the analytical solutions derived by Pyrak-Nolte. The study of the effect of the joint orientation relative to the wave propagation was studied as well by Wang using a coupled method combining LS-DYNA and UDEC (Itasca, 2005) and found that maximum reflection is obtained when the joint dips in an angle of 45°. Chen et al. (1998 and 2000) modeled a similar problem using a two staged process; in the first stage the code AUTODYN was used to simulate the blast load, and in the second stage these loads were recorded into UDEC to investigate the blast wave propagation. Note that UDEC is a Discrete Element Method (DEM) code that explicitly models the joints of the rock mass, but cannot model the explosion process in its original formulation. The numerical model results obtained by Chen were in good agreement to a field test, and the results demonstrated the attenuation effect of joints on wave propagation A key limitation of Pyrak-Nolte’s equations, as noted in her Doctoral thesis, is that the equations are not applicable when the wavelength is larger than the joint spacing. In this case multiple reflections and transmissions interacting with each other would complicate the process. As joint spacing in rock masses is most commonly in the scale of decimeters, and wavelengths in the scale of meters, this limitation almost always prevails. Later experiments conducted by Pyrak-Nolte et al. (1990) all showed that when joints are closely spaced the final transmitted wave is larger than the solution obtained using the analytical transmission coefficient for a single joint powered by the number of joints. Cai and Zhao (2000) developed theoretical methods to solve the problem of transmission and reflection across multiple joints. They defined the transmission and reflection as a function of 9 the ratio of joint spacing to wavelength and found that from a certain ratio the final transmission coefficient remains constant. Lei et al. (2007) explored the effect of multiple joints on blast wave propagation as well and obtained results from numerical simulations that were in agreement with a field test. In light of the above findings, it is clear that joints have some role in increasing attenuation, and therefore are favorable for the tunnel durability to dynamic loads, as explicitly mentioned by many of the above authors. It is debated whether the effect of joints on reduction of tunnel strength in rock dynamics is as pronounced as it is in static rock engineering. Johnson, Rizzo and Rozen (1988) proposed an empirical correction factor to the empirical guidelines established in the ERA tests based on the Rock Quality Index (RQD, proposed by Deere et al., 1966) of the rock mass. They found that a lower RQD adversely affects the tunnel’s ability to withstand blast loads. As there are empirical equations relating the RQD to the joint frequency (Palmstrom, 1982), it can be concluded from their findings that the overall influence of the joints on tunnel durability is negative. However, they note that their findings are based on a limited database and that additional and extensive field work must be carried out to properly consider the influence of rock mass characteristics. It is noted that RQD is not related to the properties of the joints (such as roughness, aperture, etc.) and therefore could not serve as a single index to assess this problem. There is very little published work that directly addresses the problem of design of support systems for tunnels under blast loads. Further to the ERA tests, the COE manual (EM 1110-345-432, 1962) discusses in detail the issue of support design for tunnels under blasting but in a qualitative manner. Amongst the issues discussed are: 1. The effect of impedance mismatch on the tunnel liner. 2. Comparison to one-dimensional spalling phenomenon in metal bars. 3. The recommended reduction of impact energy of the spalling on the liner as a function of distance from the blast. Zhao et al. (2010) proposed a simplified method of dynamic liner design based on one-dimensional spalling similar to the manner it is discussed in the COE manual to estimate the load on the tunnel liner. This load is used as an input to a Single Degree of Freedom (SDOF) 10 model to compute the displacements and stresses in the liner. Insights and limitations of the one-dimensional analysis will be discussed in Chapters 5and 6 in detail. It is clear that under certain conditions, tunnel support such as shotcrete liners and bolts that are usually calculated to sustain the static loading may not be sufficient in the case of additional dynamic loads. Malmgrenand and Nordlund (2006) studied the behavior of roof wedges supported by shotcrete and subjected to blast-induced vibrations in the Kiirunavaara mine and found that wedges can be ejected by a dynamic load even if the static safety factor was greater than 10. Rockbursts are another form of dynamic loading that may be imposed on a tunnel. Rock bursting into the tunnel occurs as a result of high stresses in deep mines or tunnels. It is similar to spalling damage from blast loading as pieces of fragmented rock fly into the tunnel with an initial velocity. The velocities measured in rock bursting are in the range of 3 m/sec, and can exceed up to 10 m/sec (Kaiser, 2012). Similar velocities were measured under blast loads in the ERA tests. The phenomenon of smaller rocks ejected first at higher velocities followed by bigger fragments ejected at lower velocities approaching the measured PPV was observed in both rockbursting and in blasting damage in the ERA tests. Kaiser (2012) lays out general guidelines for dynamic support to withstand rockbursts, and mentions that a design tool called BurstSupport is currently being developed. An important difference between a support system in a rockburst prone tunnel to a support system designed to withstand blasting is that in the former, the addition of the confining pressures imposed by the support system act to reduce the probability of bursting to occur, whereas in the case of blasting the support can merely be designed to absorb the impacts induced by the blast. Nonetheless, it is the author’s view that work related to rockburst support and blast load support can shed light on each other, as the analysis of tunnel support subjected to a series of spalling impacts may be similar. Currently, both subjects are far from being fully developed. 2.2 The ERA Underground Explosion Tests The most extensive field tests conducted in order to determine the durability of tunnels to withstand external explosions took place between the years 1948-1952, conducted by 11 Engineering Research Associates (ERA) and ordered by the US Army Corps of Engineers (COE) and the US Bureau of Mines. In these tests a large number of explosions were detonated above tunnels, and damage zones as a function of the scaled distance from the charge were empirically defined. These tests will be used in the current research for comparison and calibration of numerical models. Three different sources that refer to the ERA tests are included in this thesis: 1. Underground Explosion Test Program Final Report (1953) - a detailed description of the testing procedure, intact rock properties, and results. The author was able to find only the report of the tests undertaken in sandstone. 2. TM-5-857 (1961)- A COE manual that summarizes the characteristics of damage for the four zones of failure for granite and sandstone, with recommendations for practical design. 3. EM 1110-345-432 (1962)- A latter version of TM-5-857, with a more detailed discussion of basic wave theory and possible failure modes. The ERA tests were conducted mainly in granite (12 tests) and sandstone (17 tests) with other tests in limestone. Single delay charges in the range of 145 to 145,000 kg of TNT were detonated above the tunnels. All charges were buried charges, detonated at scaled depths that are assumed to be fully coupled with the ground. All tunnels were unlined, with diameters of 2, 5 and 10 m. The tests result scatter was 25%. Based on observation of the results, four zones of damage from total collapse to light damage were defined, as shown in Figure 2. The same terminology of damage zones is adopted in this research. 12 Figure 2 Damage zones defined in the ERA tests End of zone of damage [scaled distance] Zone 4 Light damage Zone 3 Moderate damage Zone 2 Heavy damage Zone 1 Total collapse Material 3.02 1.87 1.43 0.71 Granite 2.94 1.83 1.39 0.69 Sandstone Table 1 Maximum scaled distance of damage zone (modified from the COE Manual TM-5-857) Table 1 shows the scaled distances of the damage zones measured from the blast to the end of the damage zones. Damage zones are independent of the tunnel diameter. It is noted in the 13 COE Manual TM-5-857 that the damage zones do not have to appear all at once; for example, a blast detonated in granite at a scaled distance that corresponds to Zone 3 will cause the damage associated with zone 3 directly below the blast, and Zone 4 type damage in the adjacent area from a scaled distance of 1.87 to 3.02. As part of the ERA tests spall velocities were measured with high speed cameras. However, the corresponding results are presumed to be inaccurate. The maximum velocities are 30, 18, and 6-9 m/sec for the first three zones respectively. In Zone 4 there is no information pertaining spall velocity, and it is presumed that in this region rocks merely drops to the floor. It was observed that the first spalls carried the maximum velocities and were of smaller size than the subsequent spalls that carried lower velocities. At close distances, pieces of rock fly inwards in the radial direction and as the distance increases the spall ejects only from the roof downwards. Results of volume and area of damage for each zone were measured. The area of damage was defined as the difference between the area of the tunnel cross section before and after the blast divided by the area of the cross section before the blast. It was found that for larger tunnels the volume of damage decreases. In the COE Manual there are charts for correction factors that need to be applied when correlating the damage area to the scaled distance. The correction factor is a function of the scaled tunnel diameter which is defined as the diameter of the tunnel divided by the cube root of the explosive weight. Table 2 presents the area of damage for tunnels with scaled diameters smaller than 0.18 and for a tunnel with a scaled diameter of 1.1, where a correction factor of 0.2 is applied accordingly. The scaled diameter of 1.1 corresponds to a diameter of 10 m under a load of 500 kg of TNT, which will be used for the numerical models presented in this thesis. 14 Zones 1 2 3 4 % of damaged area to tunnel area Tunnels with a scaled diameter less than 0.18 Complete closure 30-80 30-5 0-5 Tunnels with a scaled diameter of 1.1 20-35 6-16 6-1 0-1 Table 2 Percentage of damaged area measured in the ERA tests (modified from EM 1110-345-432) Hendron (1977) summarized additional results of the ERA tests pertaining PPV values measured in the free face of the tunnel. The mean values were found to be 12, 4, and 0.8-1.9 m/sec for damage Zones 2-4 respectively. Since the PPV in the free field is the contribution of the incident and reflected wave, then it is approximately double the PPV value that would be measured in the same distance from the blast not at a free face. Accordingly, in the current analysis comparison of PPV attenuation to the ERA tests for calibration purposes was obtained by dividing by two the values of PPV in the free field given by Hendron with respect to the mean distance to the damage zone. 3 General Methodology and Program Overview 3.1 General Methodology For designing a tunnel support system to withstand blasting it is essential to specify the different stages of the problem under observation, which are: 1. Occurrence of a detonation of an explosive fully coupled in rock. 2. The chemical reaction induces very high gas pressures that act on the surrounding rock. 3. The immediate vicinity of the blast undergoes major damage, and a shock wave is formed. 15 4. The wave propagates and attenuates until it reaches the tunnel. 5. According to the magnitude of the wave, damage to the tunnel may occur and is in the form of detachment of rocks (termed as spalling) that fly into the tunnel. When considering the methodology for analysis of blasting impact on tunnels in rock, it is important to appreciate that both blasting processes and analysis of rock mass behavior are subjects that include several inherent inaccuracies: 1. Rock masses are heterogeneous and anisotropic, and can vary within the same site in both type and strength. 2. The discontinuities in the rock mass vary in orientation, spacing, roughness, persistence, etc. 3. Rock mass strength is typically estimated using empirical methods. 4. Wave attenuation and PPV exhibit large scatter in results, and can be established only by best fits to trial blasts measurements. 5. The impact of multiple joints on wave transmission and reflection is not yet properly understood. 6. Rock fragmentation is a complex phenomenon and is generally assessed by using empirical models. 7. The knowledge and study of rock dynamic parameters is generally lacking. 8. Explosives do not behave in an ideal manner, i.e. they do not reach their theoretical thermo-hydrodynamic value. 9. Pressure time-histories of blasting exhibit a great deal of irregularity. In order to address both the complexity of the problem under consideration and the inaccuracies that may arise to the points mentioned above, the analysis is divided into stages. Each stage is compared to experimental data to ensure that modelled results are within a realistic range. Models are also set up in a way that they best address the questions that the models are supposed to answer. Elmo (2006) defines an effective numerical model as one that captures the key features and attempts to simplify reality rather than replicating it. An effort is made throughout the model 16 setup and calibration processes that both materiel parameters and the chosen failure criterion are kept as simple as possible. It is argued that modeling methods that require a large amount of input parameters that are hard to estimate and/or require complicated and expensive laboratory tests would be of little practical use. The input parameters follow values commonly used for static analysis. As discussed in Section 2.1, the response of rock to dynamic loads is affected by the loading rate. However the knowledge in this field is still lacking. As discussed in Chapter 2, different authors have argued that there is some proportionality between the static and dynamic strength of the rock and that the dynamic parameters do not vary greatly from the static parameters. Therefore it is the author’s opinion that using values known from static analysis should not fail to capture the general trends when comparing different rock types. 3.2 Numerical Approach To the author’s knowledge, there is no software program that can model the problem as a whole, and it is argued that the different problem steps require using different numerical programs. An integrated approach is used in this current analysis, based on the numerical codes AUTODYN (Century Dynamics, 2014) and ELFEN (Rockfield, 2014). For the first step of estimating the load generated by the explosive detonation, AUTODYN is used. AUTODYN is a finite-difference software, specially designed to solve a wide variety of non-linear problems including the explosion of condensed high explosive material. AUTODYN is capable of generating pressure-time history of explosives using the Jones-Wilkins-Lee (JWL) equations (AUTODYN User Manual, 2014). ELFEN is a combined finite/discrete-element program (FEM/DEM) for 2D and 3D modeling of jointed rock subjected to quasi-static or dynamic loading conditions. Its ability to reproduce quasi-brittle fracture has been demonstrated through application to a large number of materials including ceramics, concrete and rock. ELFEN uses fracture mechanics principles integrated within the FEM/DEM formulation to realistically reproduce the behavior of fractured rock materials (Rockfield, 2005). In this thesis ELFEN is used for 2D the modeling of wave propagation and its interaction with a tunnel. The Rankine Rotating Crack is used as the failure criterion. The parameters 17 required for this criterion are the tensile strength of the rock and the fracture energy. The Rankine rotating crack failure criterion is based on the concept of Mode 1 fracturing studied in fracture mechanics and is well suitable for representing brittle materials. Once the maximum principal stress reaches the tensile strength the elastic modulus is degraded in the direction of the major principal stress invariant (Rockfield, 2005). The yield surface and softening mechanism are depicted in Figure 3 below. Figure 3 Yield surface and softening curve for the Rankine rotating crack in ELFEN (modified from ELFEN user’s manual, Rockfield 2005). It is noted that ELFEN enables selecting a minimum element size for the formation of additional fracturing that is independent of the initial element size. Throughout all models undertaken in this work the minimum element size for fracturing was chosen to be equal to the initial element size. This was done to avoid long computational times, as the time step in explicit modeling is dependent on the minimum element size. The fracture energy is the integral of the stress over the strain in the softening portion of the stress-strain curve shown in Figure 3. The fracture energy 𝐺𝐹 is related to the critical stress intensity factor 𝐾𝐼𝐶 and elastic modulus E (Rockfield, 2005): 𝑮𝑭 =𝑲𝑰𝑪𝟐𝑬 [4] 18 This equation implies that rocks with a higher elastic modulus will have a lower fracture energy and will fail in a more brittle manner. Fracture toughness can be measured in a laboratory. Alternatively, fracture toughness can also be estimated using an empirical relationship to the rock’s tensile strength (Zhang, 2002): 𝝈𝒕 = 𝟔. 𝟖𝟖𝑲𝑰𝑪 [5] Therefore it is possible to define the fracture energy using the two equations above based on knowledge of the elastic modulus and the tensile strength. This approach is used for estimation of the fracture energy used in the ELFEN models presented in this thesis. A major advantage of simulating damage to rock with the fracturing feature in ELFEN is that fracturing results enable an explicit and immediate depiction of damage. In other works carried out (and mentioned in the literature review) where continuum based codes were used to assess tunnel damage from blasting, the analysis of results included observation of the plastic zone and/or measurement of the PPV. It is clear that fracturing results give a better indication of failure, as the plastic zone does not necessarily imply failure. As for using PPV as an indication of failure, the author believes that for simulating the problem of tunnels under blast loads, stress is a better damage indicator than PPV, as rock fails due to stress. As can be seen in Equation [3], stress is proportionally related to PPV through the acoustic impedance of the rock. Although rocks with higher acoustic impedance typically have a higher limiting strength, the actual limiting factor is stress rather than PPV. Moreover, it is the author’s view that the proper methodology should use PPV attenuation to calibrate the numerical model according to data collected from the site rather than using it is a failure indicator. It is noted that in continuum based codes fracturing can be simulated using erosion controls where elements are deleted once they exceed a defined limit (e.g. maximum stress or strain). The setback of this method is that conservation of energy is not maintained as parts of the body disappear, whereas in the fracture simulation in ELFEN the total mass remains constant throughout the solution process. Other codes that use different variations of DEM (discrete-element-modeling) for the simulation of brittle fracturing are available. Amongst them are the Particle Flow Code (PFC) 19 and the Voronoi tessellation generator in UDEC (Itasca, 2012). Kazerani and Zhao (2010) and Lisjak and Grasselli (2014) review these methods and discussed some of the disadvantages associated with the more complex calibration processes required by those models. A unique capability of ELFEN is its ability to explicitly model joints without limitations of block creation (Elmo, 2006). Preliminary work carried out by the other showed that models as the one depicted in the figure below, yielded results that were in agreement with the analytical equations derived by Pyrak-Nolte (1988). However, they did not conform to the concept of conservation of energy as defined by Pyrak-Nolte: 𝑻𝟐 + 𝑹𝟐 = 𝟏 [6] where 𝑇 and 𝑅 are the wave transmission and reflection coefficients respectively. Figure 4 Bar model with joint to investigate transmission and reflection through joints in ELFEN When adding multiple joints to the model, the transmission and reflection results were not consistent with works carried out by others (Cai and Zhao, 2000; Lei et al. 2007; Zhou et al. 2012). With lack of substantial knowledge of the effect of multiple joints on wave propagation it was decided that for the models in this thesis the rock mass would be modeled as an equivalent continuum and that wave attenuation will be achieved by means of calibrating the damping percentage to attenuation measured in the ERA tests. 20 4 Assessment and Characterization of Blast Load 4.1 Methodology The first step of solving the problem of a tunnel under blast loading is the estimation of the load created by the release of explosive energy. The numerical code AUTODYN (Century Dynamics, 2007) is used for this stage. AUTODYN consists of coupled Euler-Lagrange solvers. Both Euler and Lagrange solvers run simultaneously and interact with each other. The material around the explosive is meshed using an Euler solver where the mesh is fixed in space and the material flows through the mesh, thus numerical errors and instabilities caused by large deformations are avoided. Conversely, the rock material is meshed using a Lagrange mesh which is more suitable for modelling the rock farther away from the blast where smaller deformations are expected. The modeled region around the detonation point undergoes large deformations where the volumetric strains have a greater impact than the deviatoric strains. The relationship between pressure and volumetric strain is defined through the Equation of State (EOS). Knowledge of rock input parameters for the EOS is hard to obtain as laboratory tests for establishing these are expensive and few published experiments are available. To overcome this issue, numerical gauges are setup in the model so that the PPV attenuation is recorded and compared with the ERA test results. The bulk and shear moduli for the EOS are initially inserted as the average static moduli as documented in the ERA tests, and calibrated through a trial and error process. Once the AUTODYN model is calibrated to the ERA tests the corresponding pressure time history can be used as an input in the ELFEN models. In the region close to the blast the wave is in the form of a shock wave and travels at sonic speed, conforming to the Rankine–Hugoniot relations which are calculated by AUTODYN. The integrated AUTODYN-ELFEN modelling approach is designed to take advantage of the strengths of AUTODYN for blast modelling. The large deformations that the sonic region undergoes due to the extreme pressures generated by the blast are accommodated by the AUTODYN coupled Euler-Lagrange solver. In order to avoid errors stemming from large deformations the stress input for the ELFEN models are defined using the history points at a 21 scaled distance of 0.3 away (which is 2.4 m for 500 kg of TNT) from the blast location in AUTODYN. The procedure is schematically illustrated in Figure 5. Figure 5 Procedure for blast load estimation from AUTODYN to ELFEN model Due to the limitations on accurate design of blasting impact on tunnels, the main objective of this modeling exercise is to assess the peak pressure and the dominant frequency in order to establish the general magnitude of the impact. Since this thesis will focus on the far field impact of blasts, an attempt to accurately model the closer vicinity of the blast would be superfluous. Extensive tests of ground shock effects in soils carried out by Drake (1983) found that the pressure time histories consist of a compressive pulse with a short rise time followed by a negligible negative tensile pulse. Dominant frequency has been found to be a sufficient parameter for assessment of blasting impact (Dowding, 1985). 22 4.2 Model Setup The research in this thesis deals with a spherically shaped and fully confined blast from TNT explosive. The model is setup in 2D axisymmetric conditions in order that the 3D nature of spherically propagating blast waves is achieved while still maintaining computational efficiency. Taking advantage of symmetry, only a quarter of the full problem can be modeled. The model dimensions are 25x25 m with the TNT material located at the origin as shown in Figure 6 and Figure 7. An Euler mesh is assigned starting from the origin and extending to 5x5 m. Both the boundaries of the Euler and Lagrange domains are assigned to be impedance boundaries so that wave reflection effects are minimized. The outer boundaries are fixed. Numerical gauges are set according to scaled distances corresponding to measurements from the ERA tests. Kuhlmeyer and Lysmer (1973) found that in order to obtain realistic results in dynamic simulations the mesh size should be no larger than one-eighth of the minimal wavelength. Because of the fast computing times required by the proposed AUTODYN model, an even finer mesh size was used in the models. It is recommended by the AUTODYN User’s Manual to use Lagrange type elements that are about twice the size of the Euler elements for flow constraints to function properly. Accordingly, in the current models the mesh size for the Lagrange type rock elements was 9 cm, while a mesh of 4.5 cm was used for the TNT and Air Euler elements. 23 Figure 6 Gauge layout for AUTODYN model Figure 7 3D view of the axisymmetric model in AUTODYN 24 In order to conduct a dynamic analysis in AUTODYN, each material model requires a strength model, a failure model, and an equation of state (EOS) to be defined. The properties of TNT were obtained from AUTODYN’s material library and are listed in Table 3. The EOS for TNT is assigned according to the empirical JWL equations. Property Value Density [kg/m³] 1630 CJ detonation velocity [m/sec] 6930 CJ Pressure [GPa] 21 Table 3 TNT properties The input parameters for the air material were also obtained from AUTODYN’s material library. The EOS for air is the ideal gas equation. The internal energy corresponding to the atmospheric pressure was assigned for the air material as an initial condition. The mean values of the sandstone material properties from the different sites of the ERA tests are listed in Table 4. Value Property 2300 Density [kg/m³] 17.2 Young modulus [GPa] 86 Compressive strength [MPa] 1.8 Tensile strength [MPa] 0.2 Poisson Ratio Table 4 Sandstone averaged intact properties from the ERA tests. Note that the sandstone properties from the ERA tests refer to intact rock material. No rating of the rock mass with respect to the discontinuities is available from the ERA tests since the tests were conducted before the development of the rock classification systems used today. It is known that the presence of joints causes a reduction in strength. However the ERA tests do not explicitly take this important aspect under consideration. 25 It is assumed that the rock mass was relatively massive as no support was installed in the tunnels where the ERA tests were conducted. In weak rock masses excavation without support for long periods of time is not possible. Referring to the stand-up time chart developed by Bienawski (1989) since the ERA tunnels did not require support, it is inferred that the sandstone rock was characterized by a relatively high RMR rating. RMR (Rock Mass Rating, proposed by Bienawski 1989) is a commonly used rock mass classification system used in rock mechanics (Brady & Brown, 2004). Similarly, GSI (Geological Strength Index, proposed by Hoek et al., 1995) is another rock mass classification system that has the advantage of being linked to the Hoek-Brown failure criterion system (Hoek et al., 2002). Assuming massive rock mass conditions, the author believes that a GSI rating of 80 could be assumed for the ERA sandstone material. Using the Hoek-Brown failure criterion, the material properties used in the model to represent the sandstone material are listed in Table 5. The equivalent Mohr-Coulomb parameters are also indicated. 80 GSI Hoek-Brown parameters 0.2 Poisson ratio 17 mi 0 D 64 Friction angle (degrees) Equivalent Mohr-Coulomb parameters (for sig3max=0.75MPa) 3.3 Cohesion [Mpa] 15,200 Erm [Mpa] Table 5 Mohr-Coulomb fit to Hoek-Brown failure criterion for the sandstone material. The first four rows in Table 5 are used as inputs for the Hoek-Brown equations and the three bottom rows refer to the parameters used as inputs for modeling the sandstone in AUTODYN. The parameter mi can be found based on the rock type where 17 is the mean value for sandstone (RocScience, 2002). D is the blasting damage factor for tunnels where poor quality drill and blast through their excavation process was used, and is assumed to be zero. 26 Although AUTODYN does not have a Mohr-Coulomb failure model, the linear Drucker-Prager strength model can be used for simulation of the sandstone since the Drucker-Prager failure criterion is a generalization of the Mohr-Coulomb model (AUTODYN user manual, 2014). The Drucker-Prager model is suitable for modeling pressure sensitive materials such as rock, and is defined by a friction angle and cohesion strength parameters. It differs from the Mohr-Coulomb model in that its 3D surface is circular, rendering computer implementation easier. For AUTODYN, values for the yielding stress at zero pressure and the slope of the failure surface must be assigned; these are assigned using the equivalent Mohr-Coulomb parameters listed in Table 5. Failure of the rock is defined as a function of plastic strain. In general, rock is a relatively brittle material and cannot undergo a large amount of strain. Therefore a plastic strain limit of 0.0025 was assumed in the current models. Other control parameters (e.g. element erosion, damping) were kept as default. EOS (equation of state) for a material defines the relationship between the hydrostatic pressure, density, and internal energy (or temperature). In modeling of blasting in rock it is assumed that temperature can be neglected (Chen et al., 2000), therefore the EOS is a function of pressure to volume alone. This relationship is affected by the strain rate of the load applied, implying that the bulk modulus increases with increasing strain rate. Data regarding the impact of strain rate is hard to obtain and requires laboratory testing methods that are both expensive and not yet well established, therefore a simple linear equation of state was assumed. The measured elastic constants for the intact rock (Table 4) were first reduced using the Hoek-Brown GSI system (Table 5). The initial bulk and shear moduli are gradually increased in order that PPV results are consistent with the ERA tests. In order to minimize scatter in modeling results two models are undertaken, one with an explosive weight of 200 kg, and one with an explosive weight of 500 kg. Gauges are set in the scaled distances corresponding to the four zones of failure in the ERA tests as shown in Figure 6. Gauges are set in the x and y axis, and along a line inclined 45°, and the results are compared with the ERA test results, as discussed in the following section. 27 4.3 Analysis of Results and Discussion Table 6 shows the results for the numerical gauges presented in Figure 6. The results for gauges 5-8 are obtained by summing the vector components in the x and y directions. The results of all gauges at equal distances are averaged for the results of both the 200 kg and the 500 kg models. Generally, the modeling results show a good agreement with the scaling law that allows comparing different explosive weights by using a scaled distance. 200 kg of TNT 500 kg of TNT Average PPV [m/sec] Average Stress [MPa] gauge PPV [m/sec] Stress [MPa] PPV [m/sec] Stress [MPa] 1 80.8 360 33 250 5 36.6 310 58.2 374 9 53.0 281 54 347 Average 56.8 300 48.4 320 48.4 320 2 6.3 50 5.9 41.5 6 2.2 13 5.2 32.4 10 2.4 17 5.6 34.6 Average 3.6 27 5.6 36 5.6 33 3 2.0 18 4.1 30 7 1.6 9.1 3.1 19.1 11 1.6 12 2.6 21.9 Average 1.7 13 3.3 24 3.3 18.5 4 0.9 6.8 1.9 16 8 0.8 4.9 1.6 8.35 12 0.8 5.9 0.9 5.13 Average 0.8 5.88 1.5 9.82 1.5 7.8 Table 6 PPV and stress results from AUTODYN model. 28 Figure 8 Comparison of AUTODYN results with results from ERA tests and Hendron’s tests. Figure 9 Stress time-history for Gauge 1. 29 Figure 8 shows the modeling results in comparison with results of the ERA tests and Ambraseys and Hendron’s PPV attenuation formula (1968). In the intermediate scaled distances of 1.4 𝑚/𝑘𝑔1/3 results are in good agreement, while the AUTODYN model results are seen to attenuate less rapidly with increasing distance. However, this is not a concern as the objective is to obtain the magnitude of the load at a distance closer to the blast. Table 7 lists the final dynamic elastic constants used in the models after the trial and error process. With comparison to the initial static constants an increase of 25% is shown. Dynamic Static 0.2 0.2 Poisson ratio 18.9 15.1 E [GPa] 7.9 6.3 G [GPa] 10.5 8.4 Bulk Modulus [GPa] Table 7 Elastic (dynamic) constants used in AUTODYN model Pressure-time histories are observed in order to assess the dominant frequency of the blasting loads. There are different methods of assessing the dominant frequency, most notably FFT (Fast Fourier Transformations) are used for this purpose. In this work a simplified method was used: the time difference between the highest peak and following peak is used to give a good estimation of the time period and dominant frequency (Dowding, 1996). The results show that the time periods at a distance close to the blast range from 0.5 to1 milliseconds. The load inserted into the ELFEN models calibrated to the ERA tests will therefore be 320 MPa, which is the average stress recorded in a scaled distance of 0.3 m/kg13⁄ (the average of gauges 1, 5, and 9 shown in Table 6). The negative pulse that follows the positive compressive pulse is neglected, and the rise time is assumed to be a small portion of the total duration (Drake, 1983). The time period is chosen to be 8 milliseconds, which is the approximate average of the time periods recorded at this distance .The load function shown in Figure 10 below is used for all ELFEN models in this thesis. 30 Figure 10 ELFEN normalized face load time curve. 00.20.40.60.81-0.1 0.1 0.3 0.5 0.7 0.9Normalized Load Time [ms] 31 5 One-Dimensional Analysis of Blast Load Using ELFEN 5.1 Introduction As discussed by Jing (2003), the most important step in numerical modeling is not running the calculations, but the conceptualization of the problem regarding the dominant processes, properties, and significant parameters. In this context, the use of one-dimensional simplification before performing more complex simulations is helpful in gaining a better understanding of the failure mechanism. Subsequently, the 2D results will be compared to the 1D solutions to determine whether the 1D simplification is acceptable with respect to the general problem. This study will focus on examining the phenomenon termed spalling in which the reflecting tensile wave generates fragments of rock that detach from the tunnel boundaries and eject into the tunnel space with an initial velocity. The 1D simplification is discussed in detail in the COE Manual (EM 1110-345-432, 1962) and used by Zhao et al. (2010) as a basis for developing a simplified approach to dynamic liner design. The following discussion is based mainly on these two references. The COE Manual illustrates the spalling mechanism using a 1D example of a compressive pulse that reaches and reflects back from a free surface, as depicted in Figure 11 below. The theory is backed by laboratory experiments on steel bars subjected to explosive pulses. 32 Figure 11 Spalling process due to a compressive pulse reflecting at a free face In Figure 11 shows a compressive stress pulse travelling to the right with the tensile reflection illustrated at the bottom and travelling to the left. The stress at the free face will always equal zero because the transmission and reflection of the wave cancel each other. According to the notation used in Figure 11, the first spall will form when: 𝝈𝒎 − 𝝈𝒙 = 𝝈𝒕 [7] where 𝜎𝑚 is the peak stress of the wave, 𝜎𝑡 is the tensile strength of the rock material, and 𝜎𝑥 is the compressive stress of the incident wave at the point of the first spall. Using geometry and given the wave velocity 𝑐, the above equation can be rearranged to find the thickness of the first spall ℎ as: 𝒉 =𝒄𝒕𝟐𝝈𝒕𝟐𝝈𝒎 [8] where 𝑡2 is the time of the descending part of the wave. It is noted that the rise time of the wave does not affect the results. 33 The PPV of a tensile wave is opposed to the direction of its propagation and therefore the velocity of the first spall is the PPV contribution from both the incident wave and the reflected wave. Using the relationship between stress and PPV given in Equation [3], the velocity of the first spall 𝑉𝑜 is given by: 𝑽𝒐 =(𝝈𝒎+𝝈𝒙)𝝆𝒄=𝟐𝝈𝒎𝝆𝒄(𝟏 −∆𝒕𝒕𝟐) [9] where ∆𝑡 is the time interval from the time when the peak stress reaches the free surface to the time of the first spall. As a new free face is formed, the remainder of the compressive pulse, now with a peak stress of 𝜎𝑥 will once again reflect from the newly formed free face. The second spall will occur if the tensile stress will once again exceed the strength of the rock. According to the COE Manual the number of the spalled layers will be given by the integer number smaller than the ratio 𝜎𝑚/𝜎𝑡. Equations [8] and [9] imply that an increase in the tensile strength of the rock for a given blast load will result in a thicker spall, but with a lower spall velocity; as the ratio 𝜎𝑚/𝜎𝑡 increases, there will be a larger number of spalls with greater velocity but smaller thickness. Zhao et al. (2010) used this 1D concept for establishing a simplified method for design of a tunnel support system consisting of shotcrete and bolts, assuming the shotcrete layer between each pair of bolts behaves as a beam subjected to a dynamic impact. The magnitude of the impact is assumed to be the thickness of the spall with an initial calculated velocity according to Equations [8] and [9]. Zhao's 1D simplification neglects the effect of wave attenuation. To overcome this problem, an enhanced 1D model incorporating the effect of attenuation is formulated in this research using Matlab and is presented in Appendix 1. A parametric study of values has found that attenuation of the wave will cause the first spall thickness to decrease by 5-15% with the deviation dependent on the attenuation, the wavelength and the ratio 𝜎𝑚/𝜎𝑡. The geometry of the 1D simplification differs from the case of a spherical wave interacting with an arched tunnel roof. The effect of lateral resistance cannot be captured in the 1D 34 analysis. The methodology undertaken is therefore to begin with the conceptual 1D analysis and then carry out 2D analyses to investigate the influence of the geometrical difference. For this purpose, a model of a bar subjected to a 1D compressive pulse is undertaken using ELFEN. A parametric study is then used to compare ELFEN results of spalling thickness and velocity with the above Equations [8] and [9]. 5.2 Model Setup Four models where set up with only the tensile strength 𝜎𝑡 varying while all other parameters are kept constant. The parameters are presented in Table 8 below. Note that the Rankine rotating crack failure criterion is used in all ELFEN models in this thesis. Preliminary parametric studies have shown that the fracture energy 𝐺𝑓 has negligible impact on the modeling results. Therefore fracture energy is kept constant without the need to adjust 𝐺𝑓 to 𝜎𝑡 as discussed in Chapter 3. Model # 1 2 3 4 𝜎𝑡 [MPa] 2 4 6 8 𝜎𝑚 [MPa] 10 E [GPa] 15.1 𝜌 [kg/m³] 2300 c [m/sec] 2562 𝑡2 [ms] 0.7 𝐺𝑓 [J/m²] 22 Table 8 Properties for the 1D bar models in ELFEN The model geometry consists of a bar in the dimensions of 0.5x5.0 m. Preliminary modeling has shown that the thickness of the bar has no impact on the results. The right end of the bar is constrained in both vertical and horizontal directions, and the horizontal boundaries of the bar are constrained in the vertical direction to maintain the one-dimensionality nature of the problem. The pulse is applied to the left free end as a face load with a load time curve identical to that shown in Figure 10. No damping is used in the models, as the purpose of 35 these models is to compare results to the analytical equations, which do not account for wave attenuation. Figure 12 1D Model setup and load curve As the wave travels to the right and reaches the constrained end, it reflects back as a compressive wave. When it travels back and reaches the free end on the opposite side it will reflect once again as a tensile wave. Spalling will therefore occur on the left side of the model. The mesh size is chosen to be 5 cm as shown in Figure 13. The fracture formation is constrained to occur only along the boundaries of the meshed elements. In ELFEN it is currently not possible to use a rectangular structured mesh when adopting a fracturing criterion. The possible influence of the mesh topology on the model results will be discussed in the following section. Figure 13 Mesh for bar models. 36 5.3 Analysis of Results and Discussion Spalling occurred in all four models as shown in Figure 14 below. As the wave reflects from the free end the first cracks appear. Eventually, the spalled pieces of rock detach and begin to travel to the right at a constant velocity. Figure 14 Spalling in Models 1-4 Figure 15 and Table 9 show the results for the 1D ELFEN analysis model results compared to the theoretical solutions including spall thickness, spall velocity and number of spalls. The modeled results for spall thickness are in agreement with the analytical equations. Results of the first spall velocity are relatively lower in ELFEN, but follow the analytical trend. The numbers of spalls are higher in the ELFEN models. It is recognized that the differences in velocity and number of spalls may be related to mesh effects. As the wave reflects, the cracks do not form all across the model to create an immediate spall, as is the underlying assumption for derivation of the analytical one-dimensional equations. Cracks form at one orientation of the mesh and continue to elongate after the pulse has already continued to travel by. This occurrence is assumed to be related to 37 the chosen mesh topology. As stated in the previous section, in ELFEN it is not possible to set a rectangular structured mesh in fracture modeling. In the finite element solution process the stresses are calculated at the nodes, therefore it is possible that in the current models failure occurs at specific nodes due to their orientation relative to the wave advancement. It is reasoned that the fact that the fracture does not form all at once affects the results of the number of spalls and spall velocity. Some portion of the wave that would have been trapped in the ejected portion continues to travel within the bar and creates more spalls than the analytical solutions. In turn, the same effect causes spall velocities to be lower, due to the lower amount of trapped energy. It can be argued that irregularity of the mesh serves as a good representation of reality as rocks tend to exhibit a great deal of heterogeneity at both the laboratory and field scales. It is assumed that formation of a clear straight cut through rock will seldom occur, as the fracturing process is affected by the randomness of the rock grain size and distribution, pre-existing micro-cracks, and strength anisotropy. To the author’s knowledge, no publications of laboratory experiments addressing the issue of spalling of rock specimens are available to help validate or refute this issue. Hence, it is possible that the heterogeneous nature of the rock may cause the spalling process to be more similar to the ELFEN results. It is noted that maximum velocities of spalled rocks recorded in the ERA tests were found to be even higher than the analytical predictions. This issue will be further examined in the 2D analyses. 38 Figure 15 Graphs of analytical equations with ELFEN models for (a) spall thickness to tensile strength and (b) spall velocity to tensile strength Model # 1 2 3 4 𝜎𝑡 [MPa] 2 4 6 8 Number of Spalls Analytical solution 5 2 1 1 ELFEN 5 4 2 2 Table 9 Number of spalls in the ELFEN models compared to spall theory 39 6 Two-Dimensional Analysis 6.1 Methodology In the previous chapter one-dimensional spalling was simulated with ELFEN. Results of spall thickness, spall velocity and number of spalls from ELFEN were compared with the analytical solution presented in the COE manual (EM 1110-345-432, 1962). The purpose of this chapter is to compare results of 2D models of circular tunnels with both analytical solutions and the 1D model results from Chapter 5. The suitability of the 1D simplification is questioned as blast waves propagate spherically. Moreover, tunnel roofs are typically arched rendering the 1D approximation even less suitable. The influence of the model geometry on the simulated results will be investigated. While static design of tunnels can usually be reduced to a 2D plane strain problem, blasting is a 3D problem. The blast wave propagates spherically, arriving at different locations of the tunnel boundary at different magnitudes. The 2D simplification is probably somewhat conservative as the plane strain simplification assumes that the maximum blast load acts infinitely; on the other hand the effect of out of plane shear waves that can impact on the tunnel are not accounted for. Works done by others (Zhao et al. 2012; Wang, Konietzky and Shen, 2009; Chen et al. 2000) showed that 2D models were successful in predicting blasting outcomes tested in the field. In the 2D problem the final damage extent can be measured with respect to the tunnel shape. Angle of fracturing is herein defined as the angle measured from the tunnel center to the most distant cracks in the lateral direction, as shown in Figure 16. Fracture thickness is herein defined as the thickness of fracturing above the tunnel roof and is also illustrated in Figure 16. 40 Figure 16 Definition of angle of fracturing and fracture thickness 6.2 Model Setup Four models were set up using the same properties listed in Table 8 for the 1D modeling. As there is some attenuation of the wave due to the geometrical spreading, the load assigned to act upon the arc was 16 MPa so that it would arrive at approximately 10 MPa at the tunnel roof, thus 𝜎𝑚 is the same as in the 1D models. No additional damping is assigned; as mentioned, the effect of damping on fracturing can be neglected as investigated and briefly discussed in the previous chapter. The loading function is identical to that used in Chapter 5 and shown in Figure 10. The load is applied at a distance of 7 m which is larger than two times the wavelength and assumed to be sufficient to avoid effects of additional reflections from the free face of the arc and to limit other edge effects. The tunnel diameter is 10 m. The 2D model setup is shown in Figure 17Figure 17 2D model setup., including the different mesh regions. 41 Figure 17 2D model setup. Kuhlmeyer and Lysmer (1973) found that in order to obtain realistic results in dynamic simulations the mesh size must be no larger than one-eighth the minimal wavelength. Since results could be obtained in a reasonable amount of time for this problem, an even finer mesh size was defined with the mesh size set to 12 cm in the zone of interest. In the perimeter area which does not affect the tunnel a coarse mesh size of 3 m is assigned, with a size transitioning zone between the two regions. 6.3 Analysis of Results and Discussion In contrast to the 1D models, spalling occurs only in Model #1 where the ratio 𝜎𝑚/𝜎𝑡 is 5. Results of the four models are presented in Figure 18 and Figure 19. It is evident that as the rock material’s tensile strength increases the angle of fracturing and the fracturing thickness both decrease. This is consistent with the ERA test observations where in the zones of greater damage spalling occurred over a greater portion of the tunnel boundary with rock ejecting into the tunnel space in the radial direction, whereas under less intense loading damage was limited to the tunnel roof. 42 It is noted that applying damping to simulate the effect of realistic wave attenuation would cause fracturing to spread around the tunnel circumference in a more parallel manner than as can be seen in the results below. Figure 18 Results of fracturing for: (a) Model 1, 𝝈𝒕 = 𝟐 𝑴𝑷𝒂 (b) Model 2, 𝝈𝒕 = 𝟒 𝑴𝑷𝒂 (c) Model 3, 𝝈𝒕 = 𝟔 𝑴𝑷𝒂 (d) Model 4, 𝝈𝒕 = 𝟖 𝑴𝑷𝒂 43 (a) (b) Figure 19 ELFEN Graphs of 2D analysis results for (a) Angle of fracturing vs. tensile strength (b) Fracture thickness vs. tensile strength The first cracks occur at a distance almost equal to the calculated distance using the analytical equation for spall thickness derived from the 1D simplification. Results of distance of the first crack from the tunnel boundary in the 2D analysis and the distance from the bar free face in the 1D analytical equations are plotted in the graph in Figure 20. 01020304050607080901002 3 4 5 6 7 8Angle of fracturing [degrees] Tensile strength [MPa] 01020304050607080901002 3 4 5 6 7 8Fracture thickness [m] Tensile strength [MPa] 44 Figure 20 Graphs for first spall appearance to tensile strength However, once the first cracks are formed, the fracturing evolution in the 2D models differs from the 1D models. Whilst in the fracturing process in 1D is shown to occur more as a detachment of distinct spalls, the 2D results show that final fracturing is in the form of a network of coalescing cracks. In the 2D models the fracturing may propagate downwards and to the sides as well as upwards from the first cracks. The difference between the 1D and 2D model results can most likely be attributed to the geometry of the wave and free face. Figure 21 below displays contours of principal stresses with the stress trajectories representing the vector direction of the principal stresses. Figure 21a shows the incident wave before reaching the tunnel and Figure 21b shows the wave reflecting and causing a complex stress redistribution governed by the problem geometry. The different points of the wave arrive at the tunnel roof at different time steps, reflect at different angles and interact with each other. Once cracking is initiated, the failure process is further complicated as portions of the waves interact with these new surfaces. 00.10.20.30.40.50.60.70.80 2 4 6 8 10Spall thickness [m] Tensile strength [MPa] 2D modelsAnalytical solutions45 Figure 21 Principal stress trajectories for the wave-tunnel interaction from Model #1 after: (a) 2.11 ms (b) 2.73 ms (c) 3 ms In the 2D models only in Model #1, with the weakest tensile strength, fracturing causes the rock to fully detach and spall. In this model the spalls have different velocities: the first layer of spalls consists of vertical velocities in the ranges 0.8-3.7 m/sec and the upper layers consist of velocities in the range of 0.5-1 m/sec. The spalls more distant from the center have horizontal components as well as vertical components in the ranges of 0.4-0.9 m/sec. 46 The spall with maximum velocity of 3.7 m/sec is higher than the velocity of the first spall in the 1D Model #1 and the velocity given by the analytical equation, which are 2.3 m/sec and 3 m/sec respectively. This is consistent with observations from the ERA tests where it was reported that the first spalled layer ejected at velocities greater than the recorded PPVs. Although no spalling occurs in the models with lower ratios of 𝜎𝑚/𝜎𝑡 this does not necessarily imply that these will not occur in reality. Rock masses are jointed and this presumably increases the probability of spalling. This aspect will be considered in the following section which deals with spalling of jointed rock masses. Results show that using the 1D simplification may lead to erroneous conclusions regarding liner design. As mentioned, the 2D model results demonstrate that fracturing may propagate downwards and eventually the liner itself would spall as well as the rock. This point will further be discussed and demonstrated in Chapter 9. In summary, it is seen that aside from results of the first fracture appearance, the damage results in the 2D models greatly differ from the 1D simplification. The 2D results are found to be more consistent with the findings of the ERA tests in terms of fracture propagation and spalling velocities. The difference between the 1D and 2D models with respect to spall occurrence requires further investigation, and can be determined through comparison to the ERA test results, as will be further discussed in Chapter 7. 6.4 Influence of Joints on Spalling 6.4.1 Methodology and Model Setup It is assumed that the joints in the rock mass have a significant influence on spalling. It was reported in the ERA tests that in some instances the shape of the tunnel roof after blasting was irregular, with overbreaking that occurred most likely due to existing joint sets. As mentioned, preliminary models for this work of wave propagation through joints yielded results that were not consistent with works carried out by other authors. The transmission and reflection through multiple joints has a significant impact on wave attenuation. It is noted that 47 when modeling transmission through single joints results were in agreement with analytical solutions. In order to avoid erroneous conclusions related to the presence of joints, two methods are considered as part of this research: 1. Explicitly modeling the joints but placing them vertically and only in the close area of the tunnel roof so that their interference with wave propagation is minimal. The joints are assigned distinct strength properties utilizing the Distinct Element Method capability of ELFEN. 2. Applying heterogeneity to the material strength properties. In ELFEN it is possible to specify a statistical distribution to each elastic and plastic parameter. A random number generator is used to create values for each element according to the upper and lower variation limit and type of distribution (Rockfield, 2005). As the joints in the rock mass have a weakening effect, the initiation of failure is not limited to a certain location as is the case in the models in the previous section where all elements are of equal strength properties. The effect of the material heterogeneity is therefore assumed to be similar to that of the joints. A major disadvantage of the material heterogeneity method in this case is that it is not possible to assign a statistical distribution to the tensile strength in ELFEN. As mentioned in the previous chapters, the tensile strength governs spalling failure. The models used for both methods are identical to Model #3 presented in Section 6.2. This model was chosen because no spalling occurred in it, thus it can be examined if the addition of the joints will cause spalling to occur. The geometry of the model containing the joints is presented in the Figure 22 below. The length of the 3 joints is assigned to be 1 m, so that their upper end is slightly above the area that fractured in the model without the joints. The joints are assigned a cohesion value of zero and friction angle of 30°. The joint stiffness is assigned a high value of 12 GPa/m, equal almost to the stiffness of the material; this serves as a means of minimizing wave reflections. 48 Figure 22 Model with explicit joints For the model using the material heterogeneity method normal distributions are assigned to all elastic constants (Young’s modulus, Poisson ratio, shear modulus, and density) and to the fracture energy. All parameters are given a lower bound of 35% of the mean and an upper bound of 5% in order to increase the amount of weakness zones. 6.4.2 Analysis of Results and Discussion of Model with Joints As shown in Figure 23a, spalling occurs in the model including the pre-inserted joints. Results for the model with no joints are shown in Figure 23b for comparison. The pieces of rock between the joints detach with a velocity of approximately 1 m/sec. This velocity is lower than the velocity obtained from the 1D model and from the analytical equation (1.8 and 2.4 m/sec respectively). Possibly, this reduction in velocity can be attributed to loss of energy from friction that occurs along the joints. The fracturing initiates at a manner similar to the model without joints but also propagates outwards and causes additional fracturing and spalling. 49 (a) (b) Figure 23 Spalling above tunnel roof in (a) model with joints (b) identical model with no joints With relation to disruption to the form of the wave front it can be seen that as the incident wave arrives at the tunnel roof there is no influence of the joints as the principal stresses are directed nearly parallel to the joints. Once the wave begins reflecting and the interaction with the tunnel roof causes the principal stresses to change their orientation throughout time, the influence of the joints becomes evident as the stress distribution is discontinuous. 50 To summarize, the contribution of the joints to failure is evident in the explicit joint simulation which appears to be a valid method for investigation of the impacts of joints on tunnel spalling under blast load. 6.4.3 Analysis of Results and Discussion of Model with Material Heterogeneity In the model with material heterogeneity no spalling occurs. In terms of the amount and spreading of the fractures, results are similar to the model without joints carried out in Section 6.2 with the material heterogeneity function making only a minor contribution to the formation of cracks, as can be seen in Figure 24.. Figure 24 Comparison of (a) model with heterogeneity and (b) regular model With regards to wave propagation, some minor disruption to the wave form is evident in the material heterogeneity model. In Figure 25 the contours of the principal stresses of the waves in the model with and without material heterogeneity are displayed and the discontinuous form of the wave due to simulation of material heterogeneity can be observed. To summarize, the equivalent continuum approach with material heterogeneity method does not appear to be a valid method of investigating the effect of joints on tunnel spalling under blast loads. It is possible however that enabling randomization of the tensile strength could yield improved results. 51 Figure 25 Wave propagation in (a) model without material heterogeneity (b) model with material heterogeneity 52 7 Modeling of the ERA Tests using a FEM/DEM Approach 7.1 Methodology The objective of this section is to extend the 2D investigation presented in 6 by calibrating and comparing results to the ERA underground explosion tests. The calibration process will focus on finding the appropriate damping percentage by comparing PPV attenuation to the ERA tests. The following procedure is implemented in the ELFEN models in this chapter: (1) An initial model with a scaled distance of blast to tunnel corresponding to damage zone 4 in the ERA tests is created in ELFEN. (2) The damping percentage in the model above is calibrated so that PPV results are consistent with the PPV attenuation measured in the ERA tests. (3) Other models (simulating Zones 2, Zone 3, and model of no damage) are then assigned the same damping percentage as the calibrated model. (4) Results of damage area and velocity of the spalled rock in all models are then compared to those of the ERA tests. It is noted that simulating Zone 1 with the methodology used would be meaningless; as mentioned in Chapter 4 the load recorded from AUTODYN and applied into ELFEN was at some distance from the actual blast in order that the advantages of AUTODYN in modeling large deformations are utilized. Only a thin layer of less than 3 m of rock would remain between the tunnel and arc; this would have a significant and artificial weakening effect on results. It is assumed that for Zone 2 where the distance of the tunnel to the arc is 10 m, this approximation is still within reason. Deng et al. (2014) found the initial stresses to have little influence on results in depths of up to 200 m. The ERA tests were conducted at shallow depths and therefore the existing in-situ stresses are not accounted for in the following models, and gravity loading is not applied. 53 7.2 Calibration Model Setup The material properties for the sandstone material are based on the ERA tests and are listed in Table 10 below. The same methodology used for the AUTODYN model in 4 was used for this model, namely assuming a GSI of 80 to account for strength reduction due to scaling effects. The method of obtaining the value of fracture energy 𝐺𝑓 is discussed in Chapter 3. Property Value E [GPa] 15.1 𝜌 [kg/m³] 2300 Poisson ratio 0.25 𝜎𝑡 [MPa] 4 𝐺𝑓 [J/m²] 22 Table 10 Material properties for the sandstone model in ELFEN The tensile strength used for the model exceeds the tensile strength of the intact rock specimen measured in the ERA tests, with the latter being only 1.8 MPa. The tensile strength of rock is commonly equal to approximately one tenth of the compressive strength (Sheorey, 1997). The mean compressive uniaxial strength for the sandstone in the ERA tests measured 80 MPa. Hence, a low value of 1.8 MPa although possible, does not seem likely. This low result could be related to the method of measurement, (e.g. when applying direct tensile tests often the results obtained are meaningless due to specimen preparation problems). Furthermore, it is known that the dynamic loading effect will cause the tensile strength to be higher than the static strength, as discussed in the literature review. Hence, a tensile strength of 4 MPa is used. The load is a face load applied upon the upper arc as shown in Figure 26. The shape and magnitude of the pulse is in accordance with the one established using AUTODYN and displayed in Figure 10 in Chapter 4. The radius of the arc equals 2.4 m which with respect to the explosive weight of 500 kg simulated in AUTODYN corresponds to a scaled distance of 0.3 kg/𝑚13, as shown in Figure 5. 54 Element size is chosen to be 12 cm in the high density mesh region, and the mesh outside the area of interest was set up to be as large as 2 m in order that computational efficiency is increased, similar to the models presented in Chapter 6. Two damping regions are defined as displayed in Figure 26. In the first region the damping percentage is set to a lower value because high attenuation is already caused by the heavy fracturing in the loading area. The distance from the tunnel to the blast is set as 23 m corresponding to Zone 4 type damage from the ERA tests. The size of the model is 40x30 m. Non-reflecting boundaries are an efficient means of avoiding undesirable wave reflections from the model boundaries. However, those are not available for 2D modeling in ELFEN. The model dimensions were therefore set up so that such reflections could not influence the solution and the wave completely attenuates before reaching the boundaries. (a) 55 (b) Figure 26 Model setup: (a) model constraints (b) model geometry and damping regions For the same reason symmetry is not utilized, i.e. the problem is not cut into half through the center of the tunnel’s vertical axis, as reflections from the symmetry axis could greatly disrupt the simulation. Only half of the blast is modeled as the blast would propagate upward as well. For this reason the lines above the arc are constrained in the vertical direction. After using a trial and error process, damping percentages of 1 and 6.5 were selected for regions 1 and 2 respectively. Results of PPV corresponding to these damping percentages with comparison to PPV results from the ERA tests and from Ambraseys and Hendron’s empirical PPV attenuation formula (1968) are displayed in Figure 27. 56 Figure 27 PPV attenuation results of calibrated ELFEN model in comparison to the mean attenuation established by Ambraseys and Hendron (1968) and from the ERA tests in sandstone. 7.3 Models Setup for Additional Damage Zones In the previous section calibration of the PPV attenuation was established with a model where the distance from the tunnel roof to the blast load corresponded to Zone 4 type damage. In this section models with different distances of blast to tunnel corresponding to the other zones of damage are considered. The models’ setup is identical to the one described in the previous section where only the distance from the tunnel roof to the blast, namely the center of the loaded arc, is altered. These distances and the corresponding damage zones are listed in Table 11. Model #3 is set to a distance where no damage should occur. 57 Model 1 2 3 Zone 2 3 - Scaled distance from blast [kg/m13] 1.26 1.76 3.28 Distance from blast [m] 10 14 27 Table 11 Distances of models corresponding to damage zones 7.4 Analysis of Results and Discussion A qualitative description of explosive interaction with rock is summarized by Brady and Brown (2006). Subsequently to blasting, the region surrounding the blasthole experiences expansion and dense fracturing, while farther away from the blast radial cracks appear. Finally, when the wave is reflected from a free face (in this case the tunnel boundary) additional fractures may appear. It can be seen in Figure 28 that the fracture pattern generated by the ELFEN model matches the description given by Brady and Brown. 58 Figure 28 Fracture patterns from (a) Rock Mechanics for Underground Mining, (modified from Brady and Brown, 2006) (b) ELFEN Model #1 10mRadial CracksDensely fractured zoneTensile reflection fracturingTunnelBlast LoadRadial CracksDensely Fractured zoneBlasta)b)59 According to the explanations given by Hagan (1973) and Bauer (1978) crushing occurs around a blasthole wall when the pressure in the detonation front exceeds the dynamic compressive strength of the rock. Fracturing in ELFEN using the Rankine failure criterion is caused only by tensile failure. Although the loading pulse is purely compressive, tensile fracturing still occurs around the blasthole. The modeling results are explained in the context of: 1. Both compressive stress pulses in the radial and tangential directions are followed by a small tensile tail. This can possibly be attributed to the unloading effect. 2. In the tangential direction the expansion of the blast hole can cause tensile hoop stresses. 3. Once a fracture is created due to one of the above reasons the stresses redistribute; hence cracks may continue to propagate. Figure 29 below shows the stress history for a node below the loaded arc in the radial and tangential directions in black and green, respectively. The results appear to corroborate the first process responsible for fracture propagation. The second reason mentioned above is not apparent in the results: tensile stresses in the tangential directions should have occurred concurrent to the compressive stresses in the radial direction and this does not appear in the stress history. 60 Figure 29 Principal stresses time-history for radial and tangential stresses (in Pa units). 61 Figure 30 Fracturing results for Models 1-3 at t=10 ms corresponding to (a) Zone 2 (b) Zone 3 (c) Zone 4 and (d) Zone with no damage. 62 Figure 31 Spalling in Zone 3, full model and close-up Results of fracturing for models 1-3 and for the Zone 4 model are displayed in Figure 30. It is noted that in order to display the spalling in a vivid manner, such as shown in Figure 31, long solving durations are required so that the spalls gain noticeable distance from their initial location. It is assumed that the areas that are densely fractured around the tunnel boundary will eventually fully detach and fall into the tunnel, as established in Chapter 6. These densely fractured areas are measured for Models 1 and 2 to compare to the damage areas recorded in 63 the ERA tests (listed in Table 2). The damage areas for the models of Zone 2 and Zone 3 are measured as 7% and 4% respectively. These numbers fall within the ranges of the ERA tests. The velocities of the detaching spalls in the vertical direction in Zone 2 range from 5-21 m/sec, with the velocity decreasing for each subsequent layer. Comparing to the ERA tests ELFEN results are within reason as the maximum velocities measured in the ERA tests where 18 m/sec. The results for damage area in the ERA tests spread over a significantly wider range for Zone 2. Observation of fracturing results of the Zone 2 model may provide an explanation for this: at closer distances the tangential cracking extending from the area of the blast may coincide with the fracturing damage. These fractures may coalesce and result in greater damage. Taking into account the random nature of fracturing, this may contribute to the greater variability of results. For Zone 3 the velocity of the first layer of spalls is within the range of 7-8 m/sec which is within the range of velocities recorded at the ERA tests (6-9 m/sec). The following layers of rock that spall have velocities in the range of 0.8-5 m/sec. The phenomenon observed in the ERA tests where the maximum velocities belong to small pieces of rock that eject first is reproduced in both the Zone 2 and Zone 3 models and can be observed in Figure 31. For the Zone 4 Model the angle of fracturing is measured to be approximately 40º and the fracturing thickness is approximately 30 cm. Results from the ERA tests characterized Zone 4 type damage as light damage with the damaged area ranging from 0-1% of the tunnel area and the spall velocity to be very low or almost zero. Spalling in Zone 4 was assumed to result due to the presence of joints. It is not possible to conclude from the simulated fracturing results what the damage area would be. It was demonstrated in Section 6.4 that the joints and heterogeneity of the rock have a significant impact on the spalling process. It can be concluded that fracturing results are in general agreement with the damage description of the ERA tests as the angle of fracturing is low, and falling of rocks would occur only due to the contribution of joints. In Model #3 (model of the zone with no damage) no fracturing above the tunnel occurred, as anticipated. 64 Maximum spall velocity [m/sec] Damaged area [%] ERA ELFEN ERA ELFEN Zone 2 18 21 6-16 7 Zone 3 6-9 8 1-6 4 Zone 4 0 0 0-1 NA Table 12 Comparison of results for damage zones from the ERA tests and ELFEN models As noted, no gravity was applied to the models and therefore the measured spall velocities are purely a function of their initial spall ejection. In the context of support design (that will further be discussed in Chapter 9), a support system would absorb only the initial velocity, and therefore velocity increase from gravitational acceleration would have no effect on the loading of tunnel support. To summarize, after the calibration process of the wave attenuation, the FEM/DEM approach is shown to yield tunnel response in agreement with the physical tests, as summarized in Table 12. In addition, in the area of the blast the fracturing pattern matches the qualitative results encountered in the field and described in the literature. 65 8 Comparison of Different Rock Types 8.1 Methodology This chapter presents numerical modeling results for additional rock types in an attempt to examine the influence of rock mass strength on tunnel response to blast loads. In the previous chapter it has been established that results of numerical FEM/DEM model results are consistent with underground explosion field tests. The following chapter attempts to study and compare the response of a tunnel in different rock types. RocLab and the GSI system (RocScience, 2002) are used for obtaining typical values of different rock types, as will be discussed. For a given distance of a tunnel to the detonation point, the blast wave attenuation has a major impact on failure. It has been found in the previous models that the ratio 𝜎𝑚/𝜎𝑡 determines the tunnel’s response to the blast load, where 𝜎𝑚 is defined as the peak stress of the wave arriving at the tunnel free face. For a given load, the wave attenuation in different rock masses will vary and 𝜎𝑚 will vary accordingly. In Chapter 6 the effect of varying 𝜎𝑡 while keeping 𝜎𝑚 constant was studied. In this chapter 𝜎𝑡 varies according to the rock type and 𝜎𝑚 varies according to the resultant wave attenuation. In an elastic analysis, the wave will attenuate only due to geometrical spreading, with no dependency on the material properties of the medium. When the Rankine fracturing failure criterion is used, different rock types will respond differently to the load and the resultant attenuation will vary due to energy dissipated through fracturing. It is known that for massive rocks the attenuation will be less rapid than in weak and fractured rocks. Data from Johnson, Rozen and Rizzo (1988) and Zhou and Ong (1996) shows that the higher the acoustic impedance of the rock mass, the lower the attenuation coefficient. Different empirical attenuation equations have been proposed and are summarized by Zhou (2011). Wave attenuation is considered to be site specific. Furthermore, even in tests undertaken in the same site a large scatter in results is often encountered (Dowding, 1996). Joints have a significant influence on wave attenuation and are not simulated explicitly in this study as discussed previously. 66 Due to these limitations it is assumed that true wave attenuation rates cannot be predicted by means of numerical simulation alone. The purpose of the study in this chapter is to examine if the general trend is consistent with typical values encountered in the field and mentioned in published literature. Typical values of the attenuation coefficient 𝑛 are 1.5 for hard and strong rocks and greater than 2 for weak and soft rocks (Zhou, 2011). The hypothesis being tested is to confirm whether an initial comparison between rock types based on their strength properties can be made using ELFEN and the GSI rating system. It is customary that sensitivity analyses are conducted by varying one parameter while keeping the others constant so that the impact of each parameter is studied. However, using this methodology to study behavior of rock could yield superficial results. For example, increasing the tensile strength to high levels while the Young’s modulus is kept constant would result in a deformable rock with high strength, which would not represent true material conditions. As there is assumed to be some covariance between rock properties, the methodology used for this section is to attempt to anticipate the behavior of additional rock types using typical values for these rock types. In addition to the sandstone material modeled in the previous chapter, models with properties typical for limestone, basalt and granite are considered in the analysis. These types of rock were chosen as their strength spans from weak to strong; this way results can reflect the influence of the overall rock strength on the durability of the tunnel to blast loading. Consistent with this trend, GSI values from low to high were chosen for the rock masses. RocLab was used to establish the rock mass Young’s modulus and the Hoek-Brown parameter 𝑚𝑖. RocLab is a free software program manufactured by RocScience, used for determining rock mass strength parameters, based on the latest version of the generalized Hoek-Brown failure criterion (Rocscience Inc., 2002). The compressive strength of each rock type was chosen according to typical values encountered in the field; the Hoek-Brown parameter 𝑀𝑅 (modulus ratio) is set according to the mean values of each rock type. The Young’s modulus was estimated as well on the basis of the rock type using the empirical method developed by Hoek using the Modulus Ratio 𝑀𝑅 (Hoek, 2006). 67 All rocks were assumed to have a Poisson ratio with the typical value of 0.25. The shear modulus was calculated according to the Poisson ratio and Young’s modulus using the elastic relationship. The tensile strength is assumed to be one tenth of the compressive strength (Sheory, 1997). The fracture energy was calculated based on the tensile strength and Young’s modulus using Equations [4] and [5]. The properties of the rock types are listed in Table 13. Limestone Sandstone Basalt Granite GSI 50 80 70 90 MR 900 - 350 425 sigci [MPa] 30 80 100 150 Tensile strength [MPa] 3 4 10 15 Gf [J/m²] 23 22 83 78 Density [kg/m³] 2600 2300 2700 2700 Erm [GPa] 8.3 15.1 25.6 61.1 Shear Modulus [GPa] 3.3 6 10.2 24.4 Poisson ratio 0.25 0.25 0.25 0.25 Table 13 Properties for the different rock types The model setup (i.e. mesh size, damping, constraints, load functions) is identical to the one presented in Section 7.2. The distance from the blast to tunnel is 14 m corresponding to Zone 3 type damage. The ERA results for granite found that damage results in granite were found to be only slightly favorable compared to the sandstone (see Table 1). A limited number of tests were conducted in limestone and were found to be in the same range of results as well. The author was not able to obtain the detailed report of the granite field tests and the elastic properties of the intact granite specimens. However, granite properties encountered in the field are consistently in the strong end of the rock strength spectrum. The minor difference in damage results between the weak sandstone and the presumably strong granite appear as counter intuitive at first sight. The following section will address this phenomenon. 68 8.2 Analysis of Results and Discussion Results are shown in Figure 32. Observation of fracturing results in the area above the tunnel roof reveals that there is only minor decrease of the fracture thickness and the angle of fracturing as the rock mass is stronger. Contrastingly, in the area close to the blast the densely fractured area is greater for the weaker rock types. Figure 32 Fracturing results for the 10 m diameter tunnel in different rock types: (a) limestone (b) sandstone (c) basalt (d) granite These results are in agreement with the ERA tests for sandstone and granite, where the extent of damage was found to be only slightly more favorable for the granite. Chapter 6 showed that when equal pulses act upon the tunnel roof, damage is directly related to the tensile strength of the rock. The discrepancy between the difference in tensile strength of the granite and 69 sandstone to the similar amount of damage is attributed to the difference in attenuation. For the stronger rocks attenuation is less rapid, and the arriving pulse is larger. Figure 33 shows the peak stress as a function of the scaled distance for the different simulated rock types. It can be seen that the stronger the rock, the stronger the arriving peak stress. Figure 33 Stress attenuation vs. scaled distance as measured for the simulated rock types As the mesh is identical for all models, four nodes with varying vertical distance from the load where used to record the PPV. Using a curve fitting tool in Matlab, the PPV attenuation coefficients 𝑛 and K are found for each of the different rock types. The results are presented in Table 14 below. 70 Attenuation coefficients K n Stress [MPa] Limestone 15.8 2.5 Sandstone 16.4 2.6 Basalt 15.5 2.4 Granite 9.9 1.9 Table 14 PPV attenuation coefficients for the different rock types Results show that the FEM/DEM simulation and the Rankine fracturing failure criterion yield attenuation rates in agreement with the aforementioned values encountered in field measurements: around 1.5 for strong rocks, and above 2 for weak rocks. The numerical explanation for this is that the larger extent of fracturing that occurs around the area of the blast for the weaker rock dissipates more energy. As stated, this is not to say that numerical modeling is capable of predicting actual attenuation rates, as these are highly site specific. However, it is believed that the ELFEN models capture the general trend of the different rock types and can be used for gaining an initial estimation and understanding of the tunnel response to blasting in different rock types. In-situ measurements can later be used to modify the damping percentage so the attenuation rate is better fitted to the actual field behavior. In terms of tunnel response to blasting, the two contrasting effects of the rock tensile strength, i.e. wave attenuation and rock resistance to fracturing, even out each other. It is apparent that tunnel durability to blast load cannot be estimated in a straightforward manner; "the stronger the better" does not apply in this case, and both effects should be weighed carefully. Additional work is recommended to be carried out to better implement numerical codes as a means of attenuation prediction. Field data of wave attenuation can be correlated to the various rock and joint properties. The use of implicit and explicit modeling of jointed rock masses and different constitutive strength models can be considered and compared to establish such capabilities. 71 9 Support Design to Withstand Blasting 9.1 Introduction and Methodology This chapter will discuss the final design stage, the design of support that can prevent spalling into the tunnel. Contrary to the preceding chapters, this chapter does not rely on field tests; rather it discusses numerical modeling methods to establish the support required to withstand the spalling load. The underlying assumption for this chapter is that the demands from a support system under blasting are: a) preventing spalling and rock fall into the tunnel, b) to remain generally intact and c) to carry the dead loads subsequent to the blasting. Different damage criteria may require some modifications of the proposed modeling procedure. Rock bolts and/or shotcrete liners are widely used as support systems for ensuring the stability of underground excavations. In order to prevent spalling into the tunnel, a support system consisting of bolts alone may not be effective; even when considering a closely spaced support system, fracturing and spalling could occur between bolts. Note that from the ERA tests no qualitative conclusions regarding spall size distribution were made and the size of spalls had been found to vary between the different sites. A combination of shotcrete and bolts is arguably a more efficient means of support. As stated by Kaiser (2012) in his work related to liner design to withstand rockbursts, current design methods fail to qualitatively consider the full function of the bolts. Proper modeling of the bolts would have to consider not only the bolts themselves, but would also have to account for the complicated 3D dynamic interaction between the bolts, the grout, the liner and the rock. It is noted that bolt simulation in 2D poses an additional problem related to impedance mismatch: the simulated bolt acts as an infinite plain and would thus cause wave reflections that would not occur in reality. The addition of wired mesh or fiber reinforcement to the tunnel liner increases the liner’s bending capacity and imparts ductility to an otherwise brittle material. Similar to bolts, the modeling of wire mesh in 2D could pose technical modeling difficulties: the wire mesh would require very small elements, thus significantly decreasing computational efficiency. Additionally, the wire mesh has an important lateral resistance component. 72 In order to avoid the complexities involved with the simulation of bolts and wire mesh reinforced concrete, the support considered for this thesis is a homogenous fiber reinforced concrete liner. Note that the method of application of the liner (sprayed shotcrete, cast-in-place concrete, etc.) is assumed to have minor influence on support behavior in the blast related context. Using the terminology of the ERA tests, design of the liner for the intensive damage in Zones 1 and 2 would be impractical (EM 1110-345-432, 1962). Therefore the focus of investigations in this chapter will be on Zone 3 and Zone 4 type damage. The different damage characteristics encountered in Zones 3 and 4 were discussed and demonstrated in the previous chapters, and support design for these damage types is therefore discussed separately. 9.2 Support Design for Zone 3 Type Damage Under spalling conditions the significance of applying a layer of low density material that separates the dynamic liner from the surrounding rock is discussed in the COE Manual (EM 1110-345-432, 1962). When a wave reaches an interface between two different materials the transmission and reflection of the incident wave stress 𝜎𝐼 depends on the acoustic impedance 𝜌1𝑐1 and 𝜌2𝑐2 of the two materials. The equations for the transmitted and reflected stress derived from basic elastic wave theory are: 𝝈𝑻 =𝟐𝝈𝑰(𝝆𝟐𝒄𝟐)𝝆𝟐𝒄𝟐+𝝆𝟏𝒄𝟏 Equation [10] 𝝈𝑹 =𝝈𝑰(𝝆𝟐𝒄𝟐−𝝆𝟏𝒄𝟏)𝝆𝟐𝒄𝟐+𝝆𝟏𝒄𝟏 Equation [11] The elastic wave velocity c of a material is related to the elastic modulus and density: 𝒄 = √𝑬𝝆 Equation [12] 73 In the case of the wave travelling through the rock-liner interface, as concrete is commonly similar to rock in both density and elastic modulus, the transmitted wave will be close to the incident wave and the reflected wave will be of negligible magnitude. As demonstrated in the 2D models presented in 6, under Zone 3 type damage fractures propagate downwards from the first crack initiation. In this case a concrete liner would spall in a similar manner to the rock and be an ineffective solution. Counter intuitively, applying a high strength concrete would result in an increased transmitted wave and this could exacerbate the liner response. The solution is therefore to separate the dynamic liner from the rock by means of an intermediate layer with an acoustic impedance that is as low as possible. For example, consider a layer of styrofoam with a specific gravity of 0.6 and an elastic modulus of 0.3 GPa; sandstone with material properties from the ERA tests used in the previous chapters; and a concrete liner with an elastic modulus of 30 GPa and a density of 2500 kg; using Equations [10] and [11] will show that 96% of the wave will reflect back to the sandstone and only a negligible 8% will eventually be transmitted into the concrete liner. This shows that such a material is perfectly capable of isolating the liner from the spalling process and allowing it to function as a protective layer. Installing shotcrete upon a styrofoam layer would require some practical considerations during application, such as installing an initial thin layer with small aggregates and under low pressure, to allow satisfactory adhesion between the styrofoam and shotcrete. A clear disadvantage of such a supporting system is that the isolated liner can only serve as a protective layer for dynamic loading and therefore a separate static support system according to the rock mass properties and the static loading would have to be applied prior to the dynamic support. It is noted that steel rebar in the concrete liner could possibly decrease spalling of the concrete depending on the spacing of the bars and size of spalls. Possibly, under moderate blasting loads such a support system would suffice to prevent spalling without the need of an isolating layer. The sufficient thickness for the isolating layer must be larger than the displacement of the wave at the free face. In the model of Zone 3 type damage from Chapter 7 the maximum displacement measured at the tunnel roof before spalling begins is less than 2 mm. This in not 74 to say that a layer of over 2 mm thickness is sufficient; the engineer must account for static displacements as well. Kavvadas (2005) discusses long term creep deformations of tunnels and states that the majority of the deformation occurs immediately after excavation, but the displacement slowly increases for months after the excavation. If the protective liner is applied prior to the final deformation, the isolating layer thickness would decrease and its density would increase accordingly. 9.2.1 Model Setup To compare the performance of a tunnel liner with and without an isolation layer, two models are considered. The models are identical with the exception of the rock-liner boundary condition. For the first model, the liner is assumed to be bonded to the rock, and the adhesive strength for the bond between the liner and rock is assumed to be 1 MPa (Ahmed & Ansell, 2012). For the second model, a gap of 3 cm is left between the rock and the liner. The properties of the fiber reinforced concrete are listed in Table 15. The failure criterion for the concrete liner is a combined Mohr-Coulomb and Rankine model to account for both tensile and shear failure. The properties of the rock are the sandstone properties from the ERA tests used in the previous chapters. The ratio of the peak stress of the blast wave to the tensile strength is assigned a high value of 5 to allow intensive spalling to occur, as was observed in Chapter 6. The model geometry is identical to the 2D analyses presented in 6, with only the addition of the liner. The liner is shaped as a semicircle and constrained at its bottom. No damping is assigned to the liner as in impact problems damping has a negligible influence (Chopra, 2011). The liner thickness is gradually increased by increments of 5 cm until failure does not occur. 75 Properties of the Fiber Reinforced Concrete Density [kg/m³] 2500 E [GPa] 35 Shear Modulus [GPa] 13.4 Poisson Ratio 0.2 Tensile Strength [MPa] 6.5 Gf [J/m²] 50 Cohesion [MPa] 5 Friction Angle 52 Table 15 Properties of the fiber reinforced concrete 9.2.2 Analysis of Results and Discussion For the model with the liner bonded to the rock, under low thickness values of 20-25 cm the fracturing due to the wave reflection is limited to the rock, and the liner fails as a result of spalling impact. As the thickness is increased fracturing propagates from the rock downward to the liner, and failure is combined of both inner spalling of the concrete and rock impact. Figure 34 Simulated results of liners with 65 cm thickness for a) model with gap between liner and rock b) model with adhesive bond between liner and rock 76 Simulated results of fracturing are shown in Figure 34 where it can be seen that the bonded liner is fully fractured and damaged, and the isolated liner remains intact. The modeling results demonstrate the ineffectiveness of the bonded fiber reinforced liner to resist intensive blast damage. However, even when isolation is applied, a considerable thickness of 65 cm is required to prevent damage, which suggests that a more detailed investigation to obtain more efficient support systems is essential. A 3D analysis could account for the additional out-of-plane resistance of the tunnel liner, and would also be a better means of modeling additional support systems such as bolts, mesh reinforced concrete and steel sets embedded in concrete. However, a full 3D simulation of spalling behavior with explicit fracturing would be highly expensive from a computational perspective, and is beyond the scope of this thesis. 9.3 Support Design for Zone 4 Type Damage In models corresponding to Zone 4 type damage it was demonstrated that fracturing occurs at some distance from the tunnel roof and only joints can cause spalling to occur. In this case a single concrete liner could serve as an effective solution and should be designed to withstand the anticipated impact, and an isolation layer is not required. This assumption is valid as long the liner thickness is less than the distance of first crack appearance (which can be obtained using Equation [8]). For a thicker liner, additional consideration may be required. In the ERA tests the area of damaged rock relative to the tunnel area was for tunnels of 10 m diameter was approximately 1% and the broken pieces of rock were observed to have a velocity of 0-1 m/s. It is therefore proposed that for these cases a simple static analysis would be a sufficient means of support design. Low velocity spalls that may eject due to pre-existing joints can be accounted for by some increase of the static load. As an example, an analysis is carried out using the FEM code Phase2 (RocScience, 2011), as shown in Figure 35. A 10 m diameter tunnel is supported by liner elements. An area above the tunnel roof which corresponds to 2% of the tunnel area is assigned extremely weak strength properties to account for the 1% of rock falls measured in the ERA tests and additional movements and fractures that may occur from the blast load. On its lateral boundaries, the 77 degraded material is bound by joints, to enable increased vertical displacements. Analyses of behavior of shotcrete supported rock wedges subjected to blast-induced vibrations conducted by Malmgren and Nordlund (2006) found that for a given mass, impact loading can result in stresses and deformations even greater than 10 times the static load imposed by gravity. Thus, a load of 10 kN is assigned to account for additional dynamic loads from the blast vibration. The rock and liner properties are listed in Table 16 where the liner is assumed to be reinforced with a 30x30 cm steel wire mesh with a diameter of 6 mm. Figure 35 Phase2 Model for light damage assessment 78 Unit weight [MN/m³] GSI Young's Modulus [GPa] Rock 0.027 70 14 Degraded Rock 0.027 10 0.6 Liner 0.025 - 35 Table 16 Input parameters for static analysis Compared to heavy spalling, the stresses that will develop in the liner under small PPVs are smaller and therefore the initial state of stresses on the liner should not be neglected. Hence, the problem is staged as follows: 1. Initial stage to allow stress distribution of in-situ stresses. 2. Softening of the inclusion to account for deformations that occur prior to support installation, according to the core replacement technique proposed by Hoek et al. (2008). 3. Installation of liner. 4. Degradation of the rock material above the tunnel, activation of joints at its boundaries, and activation of vertical load. 9.3.1 Analysis of Results and Discussion Support capacity diagrams of bending moment vs. thrust are presented in Figure 36. It can be seen that for the 15 cm liner at one point the FoS (Factor of Safety) was less than one, and at 2 points less than 1.2. With an increase of liner thickness to 25 cm all points of the liner possess a FoS of greater than 1.2. For both cases, shear forces in the liner were all above a FoS of 1.2. The advantage of the proposed method is that Phase2 analyses can be set up and run at very short times, enabling optimization and/or evaluation of an array of varying conditions. 79 Figure 36 Capacity diagrams for a 15 and 25 cm liner from Phase2 models. 80 9.4 Proposed Flowchart for Support Selection To summarize the discussion in this chapter, a design flowchart is presented in Figure 37. Figure 37 Flowchart for design of support for tunnel under blasting As stated in 8, PPV attenuation is a site specific characteristic, and requires assessment using field tests, or empirical data. Once the attenuation is established, an estimation of the peak stress arriving at the tunnel boundary 𝜎𝑚 can be made. As demonstrated in Chapter 6, spalling 81 occurred only under a ratio of 𝜎𝑚/𝜎𝑡 of 5. This finding applies to the specific characteristics based on the ERA tests; different cases (e.g. different tunnel dimensions, strength parameters, loading functions, etc.) may show that spalling can occur at lower ratios. It is therefore suggested as an initial and conservative estimate to use a lower ratio of 3 to assess spalling occurrence. In the case that spalling is assumed to occur, it is recommended that an isolation layer be applied between the rock and liner. If further investigations show that a very thick liner is required, advanced 3D modeling and/or possibly simplification of the problem to a single equivalent wedge may be carried out. It is noted that such investigations were not carried out as part of this thesis. Under a low stress to strength ratio where spalling is not assumed to occur, a simple static analysis is proposed and outlined in Section 9.3. If the rock is jointed and a loose wedge above the tunnel is assumed to exist, modeling should account for the additional loads due to wedge ejection. These can be accounted within the static analysis by using equivalent static forces. It is emphasized that the design methodology is subject to various assumptions and modeling limitations mentioned throughout the preceding chapters. Thus, the flowchart is a basic outline of the design process based on judgment, and should not be seen as a comprehensive summary of the problem. The subject of support design for blasting induced damage is far from being fully developed and additional recommended work is discussed in Chapter 10. 82 10 Conclusions and Recommendations 10.1 Research Conclusions The hybrid FEM/DEM approach was used to model tunnels in rock subjected to external blast loading. Comparison to field tests verified this method as a suitable means of estimating the extent of damage to the tunnel. FEM-DEM codes combine the aspects of both finite elements and discrete elements, and incorporate fracture-mechanics principles to allow for the realistic simulation of the transition of brittle materials from a continuum to a discontinuum state. In the context of spalling damage to tunnels, this method provides important advantages compared to non-hybrid methods (FEM and DEM): a direct comparison to field test results can be made (i.e. spall velocity and fractured area measurement). Moreover, the method can readily be extended to investigate means of support, by simulating the dynamic impact of the spalled fragments with the tunnel support system. Chapter 3 lays out the general methodology and numerical approaches used for this thesis. The final results that show consistency with the ERA field tests suggest that the proposed approach, including its assumptions and simplifications, is valid. In Chapter 4 assessment of the blast load function was made using the numerical code AUTODYN. In this chapter it was assumed that a simplified triangular pulse load based on estimation of the peak stress and dominant frequency is sufficient for the subsequent modeling stage. In Chapter 5 one-dimensional models of bars subjected to dynamic pulses were presented. Results of the 1D models conducted in ELFEN were found to be in agreement with analytical solutions. In Chapter 6 2D models of a spherically propagating wave interacting with a circular tunnel were presented. Comparison to the 1D models showed that there are significant differences between the two modeling approaches and that 2D modelling is better at simulating the damage observed in physical tests in terms of spall velocity, spall size, and occurrence of spalling. In Chapter 7 the 2D models were further developed to replicate the ERA tests. An initial calibration process was undertaken in order that the wave attenuates in a manner similar to the 83 attenuation measured in the tests. Altering the blast to tunnel distance to the corresponding distances of the damage zones defined in the ERA tests showed that results of damaged area and maximum spall velocity are in agreement with the test results. Moreover, fracturing pattern in the close field of the blast was shown to be consistent with the pattern encountered in the field and described in the literature. The tensile based Rankine rotating crack failure criterion was found to be sufficient for damage simulation without the need for including a shear failure criterion, which is commonly used in FEM/DEM simulations in the field of rock mechanics. In Chapter 8 a comparison of different rock type parameters spanning from weak to strong was made. Results showed that due to the contrasting effects of the rock strength (i.e. tensile strength and wave attenuation) the tunnel response to a given blast load in weak and strong rocks is similar. Additionally, results showed that after an initial calibration process, results of attenuation rates of weak to strong rock are consistent with typical measurements encountered in the field. In Chapter 9 a discussion regarding the selection of tunnel support to withstand blasting impact is presented. Modeling results demonstrated that under heavy spalling an isolating layer between the rock and liner is critical for preventing fracturing of the tunnel liner. For light damage conditions, a simpler static FEM analysis procedure was proposed. A decision flowchart that summarizes the proposed process for support design is presented. 10.2 Recommendation for Further Studies To the author’s opinion, numerical modeling for the purpose of simulating blasting in rock is a subject far from being fully developed. Future studies that can advance the knowledge and design methods of rock blasting in general, and tunnel design to withstand blasting in particular, include: 1. Comparison of blasting field measurements with numerical simulations to investigate the ability of different numerical techniques to simulate wave attenuation. In this context, it is required to conduct further study of the effects of discontinuities in the rock mass on wave propagation. 84 2. Field experiments to further advance the understanding of the influence of rock parameters on tunnel response, and the performance of different support systems under blast induced damage. 3. Advanced 3D numerical modeling of rock spalling, to investigate 3D effects, and the performance of the addition of rock bolts, wire mesh, and/or steel sets to the liner support. 85 References 1. Ambraseys N.R., Hendron A.J. (1968) “Dynamic behavior of rock masses”, Rock mechanics and Engineering Practices, London: Wiley; 1968.p. 203–7 2. AUTODYN (2014) User Manual, revision 3.0. Century Dynamics, San Ramon, California. 3. Bauer A., Glynn G., Heater R. and Katsabanis P., (1984), "A laboratory comparative study of slurries, emulsions, and heavy ANFO explosives." Proceedings of the 10th Conference on Explosives and Blasting Technique, Society of Explosives Engineers. 4. Bieniawski, Z.T. (1989),” Engineering rock mass classifications”, New York: Wiley 5. Brady, BHG & Brown, ET (2006). Rock mechanics for underground mining (3rd Edition). 6. Cai J. G., Zhao J. (2000),” Effects of multiple parallel fractures on apparent attenuation of stress waves in rock masses”, International Journal of Rock Mechanics and Mining Sciences, 37(4): pp 661-682. 7. Cameron A., Forsyth B.,Kleine T.(2014) “Blast Design and Assessment for Surface Mines and Quarries”, online course: www.edumine.com 8. Chen M., Wenbo L., Changping Y., (2006) “Blasting vibration criterion for a rock-anchored beam in an underground powerhouse”, Tunnelling and Underground Space Technology Vol. 22 pp 69–79. 9. Chen S.G., J. G. Cai, J. Zhao, Y. X. Zhou, (2000), "Discrete element modeling of an underground explosion in a jointed rock mass" Geotechnical and Geological Engineering Vol. 18 pp 59-78. 10. Chen, S. G. and Zhao, J. (1998) "A study on UDEC modeling of blast wave propagation in jointed rock masses." Int. J. Rock Mech. & Min. Sci. 35(1), 93-99. 11. Chopra A.(2011), “Dynamics of structures” (4th Edition) (Prentice-Hall International Series in Civil Engineering and Engineering Mechanics) 12. Deere D.U., Miller R.P., (1966), “Engineering classification and index properties for intact rock”, Tech Rept. No. AFWL-TR-65-116, Air Force Weapons Lab, Kirkland Air Force Base, New Mexico, p. 300 86 13. Deng X.F., Zhu B. J., Cheb S. G., Zhou Y. X., Zhao J. (2014), “Numerical study on tunnel damage subjected to blast-induced shock wave in jointed rock masses”, Tunnelling and Underground Space Technology 43 (2014) 88–100 14. Dowding, C. H.(1996), “Construction Vibrations”, Prentice-Hall, Englewood Cliffs 15. Dowding, C. H., and Rozen, A. (1978), “Damage to rock tunnels from earthquake shaking,” Journal of the Geotechnical Engineering Division, ASCE, Vol. 104. 16. Drake, J. L. and Little, C.D. (1983), “Ground shock from penetrating conventional weapons,” Proceedings of the Symposium of the Interaction of Non-nuclear Munitions with Structures, US Air Force Academy, Colorado, Springs, CO. 17. Elmo D., (2006), “Evaluation of a hybrid FEM/DEM approach for determination of rock mass strength using a combination of discontinuity mapping and fracture mechanics modelling, with particular emphasis on modelling of jointed pillars” , PhD thesis, University of Exeter, Cornwall 18. EM 1110-345-432 US Corps of Engineers Manual (1962) "Engineering and Design of Underground Installations in Rock- Tunnels and Linings.” Armed Services Technical Information Agency, Arlington, Virginia. 19. Engineering Research Associates, (1953) “Underground explosion test program final report”, Vol. 2 Sandstone Rock, Sacramento District. 20. Hagan T.N., (1973). "Rock breakage by explosives." National Symp. on Rock Fragmentation, Australian Geomechanics Society. 21. Hendron, A. J. (1977) "Engineering of Rock Blasting in Civil Projects," Structural and geotechnical mechanics : a volume honoring Professor Nathan M. Newmark,(editor, W. J. Hall), Prentice Hall, Englewood Cliffs, Nj, pp.242-277 22. Hendrych, J. (1979) “The dynamics of explosions and its use”, Elsevier, Czechoslovakia. 23. Hoek E.(2014)," Rock mass properties",www.rocscience.com/education/hoeks_corner 24. Hoek, E., (2006) “Empirical estimation of rock mass modulus”, International journal of rock mechanics and mining science, Vol. 43, Issue 2, pp 203-215 25. Hoek-Brown failure criterion. www.rocscience.com, Toronto, Ontario, Canada 26. Itasca, (2010). Slope Model description of formulation with verification and example problems. Revision 2. Itasca Consulting Group Inc., Minneapolis, United States 87 27. Jiang N., Zhou C., (2012) “Blasting vibration safety criterion for a tunnel liner structure”, Tunneling and Underground Space Technology 32 (2012) 52–57. 28. Jing L. (2003), “A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering”, Int. J. Rock Mech. Min. Sci. Vol. 40. pp. 283-353. 29. Johnson, W.J, Rozen A. and Rizzo P.C, (1988),“Explosions in soils: the effects of soil properties on shock attenuation”. In: Proceedings, 23rd Department of Defense Explosives Safety Seminar, Atlanta, GA, USA, pp.1831–1841. 30. Kaiser P. K, Cai., M., (2012) “Design of rock support system under rockburst condition”, Journal of Rock Mechanics and Geotechnical Engineering. 2012, 4 (3) pp 215–227 31. Kavvadas M. J.(2005), “Monitoring ground deformation in tunnelling: Current practice in transportation tunnels”, Engineering Geology 79 (2005) pp. 93 – 113 32. Kazerani, T., Zhao, J., (2010). “Micromechanical parameters in bonded particle method for modelling of brittle material failure”, Int. J. Numer. Anal. Methods Geomech. 34 (18), 1877–1895. 33. Kendorski F.S., Jude C.V., Duncan W.M., (1973), “Effect of blasting on shotcrete lined tunnels”, Mining Engineering, Vol. 25, No. 12, pp. 38-41 34. Kubota S., Ogataa Y., Wadaa Y., Simangunsonga G., Hideki S, Kikuo M.,(2007) “Estimation of dynamic tensile strength of sandstone”, International Journal of Rock Mechanics & Mining Sciences 45 (2008) 397–406 35. Kuhlmeyer, R. L. and Lysmer, J. (1973), Finite element method accuracy for wave propagation problems. J. Soil Mech. Foundations Div., 99, 421-427. 36. Langefors, U., Kihlstroem, B., (1978), “The modern technique of rock blasting”, GeoRef, Copyright American Geological Institute. 37. Lei W.D., HefnyA.M.. Yan S, TengJ., (2007), "A numerical study on 2-D compressive wave propagation in rock masses with a set of joints along the radial direction normal to the joints", Computers and Geotechnics 34 pp 508–523 38. Lisjak A., Grasselli G., (2014), “A review of discrete modeling techniques for fracturing processes in discontinuous rock masses”, Journal of Rock Mechanics and Geotechnical Engineering, V 6, N 4, pp 301-314, 88 39. LS-DYNA. (2006) Keyword user’s manual. Version 970, Livermore: Livermore software technology corporation (LSTC). 40. Ma G.W., Hao H., Zhou Y.X., (1998) “Modeling of wave propagation induced by underground explosion”, Computers and Geotechnics, Vol. 22, No. 3/4, pp. 283±303, 41. Malmgren L.,Nordlund E., (2006) “Behavior of shotcrete supported rock wedges subjected to blast-induced vibrations”, International Journal of Rock Mechanics & Mining Sciences 43 (2006) 593–615 42. Palmstrom, A. (1982), “The volumetric joint count –a useful and simple measure of the degree of jointing”, IVth Int. Congress IAEG, New Delhi, pp.v221-v228. 43. Pyrak-Nolte, L. J. (1988) “Seismic Visibility of Fractures”, Ph.D. thesis, University of California, Berkeley, California. 44. Pyrak-Nolte L.J., Myer, L.R., Cook, N.G.W., (1990), “Anisotropy in seismic velocities and amplitudes from multiple parallel fractures”, Journal of Geophysical Research, v 95, n B7, pp 11,345-11,358 45. Rockfield. (2005). Rockfield Software Ltd. Technium, Kings Road, Prince of Wales Dock, Swansea, SA1 8PH, UK. 46. Rocscience Inc. (2002), RocLab Version 1.0 - Rock Mass Strength Analysis using the Generalized Hoek-Brown failure criterion 47. Rocscience Inc. (2011), Phase2 Version 8.0 - Finite Element Analysis for Excavations and Slopes, www.rocscience.com Toronto, Ontario, Canada. 48. Sheorey, P. R.(1997), “Empirical rock failure criteria”, Rotterdam: A.A. Balkema, 176p. 49. TM-5-857" (1961), "Design of Underground Installations in Rock, U.S.A.E Part 4 "Penetration and Explosion Effects". 50. Wang W., Li, Zou Y., ZhangY (2005), "3DEC modeling on effect of joints and interlayer on wave propagation", Trans. Nonferrous Met. Soc. China Vol. 16 pp 728-734 51. Wang Z.L., Konietzky H., R.F. Shen (2009) "Coupled finite element and discrete element method for underground blast in faulted rock masses", Soil Dynamics and Earthquake Engineering 29. 89 52. Wang, J, N., (1993), “Sesimic design of tunnels: A state-of-the-art approach”, Parsons, Brinckerhoff, Quade and Douglas Inc, New York. 53. Wei X., Zhao Z., (2008) "Response characteristics of underground rock cavern subjected to blast load", World Tunnel Congress, Underground Facilities for Better Environment and Safety – India 54. Wu Y.K., H. Hao, Y.X. Zhou & K. Chong (1998) “Propagation characteristics of blast-induced shock waves in a jointed rock mass”, Soil Dynamics and Earthquake Engineering, Volume 17, Issue 6, August 1998, Pages 407–412 55. Zhang Z.X. (2002) “An empirical relation between mode I fracture toughness and the tensile strength of rock”, International Journal of Rock Mechanics & Mining Sciences Vol. 39 pp 401–406 56. Zhao J., (2000), “Applicability of Mohr-Coulomb and Hoek-Brown strength criteria to the dynamic strength of brittle rock”, International Journal of Rock Mechanics and Mining Sciences 37. Pp 1115-1121 57. Zhao P J, Lok T, Yin Z, Zhou Z., (2010) “Simplified design of rock cavern concrete lining to resist shock loading”, Journal of Central South University of Technology Volume 17, Issue 5 , pp 1087-1094 58. Zhou Y, Zhao J., (2011) “Advances in Rock Dynamics and Applications-Introduction”, Advances in Rock Dynamics and Applications, CRC Press 59. Zhou Y., (2011), “Explosion loading and tunnel response”, Advances in Rock Dynamics and Applications, CRC Press, Chapter 17. 60. Zhou, Y. and Ong, Y.H., (1996),”Ground shock prediction methods – a critical appraisal”, Proceedings of the 1st Asia-Pacific Conference on Shock and Impact Loads on Structures, pp.477–483. 90 Appendix 1- Matlab Formulation for Attenuation Influence on 1D Spall % parametric study of attenuation effect on 1D spalling clear all W=500^(1/3); % scaled weight [kg^0.33] FF=15*1000; % distance from blast to free face[mm] T=0.001; % Time Period [sec] tensile=-4e6; % tensile strength of the rock [Pa] E=15.1e9 ; % Young's modulus [Pa] ro=2300; % density [kg] c=sqrt(E/ro); % P-wave velocity [m/sec] wl=1000*round(c*T); % wl is wavelength [mm] wave=zeros(1,wl); % incident wave R1=zeros(1,FF+wl); % incident in space R2=zeros(1,FF+wl); % reflection in space Total=zeros(1,FF+wl); % sum of incident and reflected waves reflect=zeros(1,wl); % reflected wave wave1=zeros(1,wl); % incident wave % variables denoted 1 refer to wave without attenuation R11=zeros(1,FF+wl); % R21=zeros(1,FF+wl); % Total1=zeros(1,FF+wl); % reflect1=zeros(1,wl); % peak1=ro*c*11.45*(0.001*(FF-wl)/W)^-2.8; % peak at arrival to free face [m] h=-c*0.8*T*tensile*1000/(2*peak1) % analytical solution for verification for j=FF-wl:FF; % formulation for attenuating wave peak=ro*c*11.45*(0.001*j/W)^-2.8; % Ambraseys & Hendron attenuation wave(1,1)=0; wave(1,0.8*wl)=peak; % m1=peak/(0.8*wl); % slope for decay for i=1:0.8*wl-1; % loop for decay portion wave(1,i+1)=m1*(i+1); end m2=peak/(0.2*wl); % slope for rise for i=1+0.8*wl:wl; % loop for rise portion wave(1,i)=wave(1,0.8*wl)-m2*(i-(1+0.8*wl)); end R1=zeros(1,FF+wl); R1(j:(j+wl-1))=wave; 91 reflect=-wave; reflect=fliplr(reflect); % mirror image of incident wave R2((FF-k):(FF-k+wl-1))=reflect; Total(1:FF)=R1(1:FF)+R2(1:FF); max_tensile=min(Total); k=k+1; if max_tensile<tensile % tension-negative convention ff=find(Total==max_tensile); % location of maximum tensile b=ff; failure_location=(FF-b) % [mm] location of failure from free face break end end k=0; % comparison for constant wave wave1(1,1)=0; % constant wave formulation wave1(1,0.8*wl)=peak1; % m1=peak1/(0.8*wl); % slope for decay for i=1:0.8*wl-1; % loop for decay portion wave1(1,i+1)=m1*(i+1); end m2=peak1/(0.2*wl); % slope for rise for i=1+0.8*wl:wl; % loop for rise portion wave1(1,i)=wave1(1,0.8*wl)-m2*(i-(1+0.8*wl)); end for j=FF-wl:FF; R11=zeros(1,FF+wl); R11(j:(j+wl-1))=wave1; reflect1=-wave1; reflect1=fliplr(reflect1); R21((FF-k):(FF-k+wl-1))=reflect1; Total1(1:FF)=R11(1:FF)+R21(1:FF); max_tensile1=min(Total1); k=k+1; if max_tensile1<tensile % tension-negative convention ff=find(Total1==max_tensile1); b1=ff; failure_location1=(FF-b1) % [cm] location of failure from free face break 92 end end delta=(failure_location1-failure_location)/failure_location % [%] difference between spalling with and without attenuation
Thesis/Dissertation
eng
Mining Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivs 2.5 Canada
Graduate
Numerical analysis of the effects of external blasts on tunnels
Text