International Journal of Rock Mechanics and Mining Sciences 104 (2018) 131–146
Contents lists available at ScienceDirect
International Journal of
Rock Mechanics and Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
Dynamic stress concentration and energy evolution of deep-buried tunnels under blasting loads
Xibing Lia,b, Chongjin Lia,b,⁎, Wenzhuo Caoc, Ming Taoa,b
a School of Resources and Safety Engineering, Central South University, Changsha, Hunan, China
b Hunan Key Laboratory of Resources Exploitation and Hazard Control for Deep Metal Mines, Changsha, Hunan, China c Department of Earth Science and Engineering, Imperial College, London SW7 2AZ, United Kingdom
T
ARTICLE INFO
Keywords:
Underground tunnel Dynamic stress concentration Energy evolution
Blasting load
Numerical simulation
ABSTRACT
1. Introduction
In recent years, the exhaustion of mineral resources in shallow depths, and the rapid development of tunneling and hydropower en- gineering, have considerably motivated the tunnel excavations to ex- tend to depth. However, due to the complicated geological environment in which deep excavations are carried out, a large number of un- conventional rock failure phenomena such as spalling,1,2 zonal disin- tegration phenomenon3,4 and rockburst hazards5,6 have been observed during underground excavations. These accidents or hazards will bring about damages to equipment and delays of excavation operation, and even pose great threats to the safety of construction personnel. There- fore, it is an urgent issue to figure out the mechanism of the engineering disasters occurring in deep excavations.
In practice, underground rocks and ores are naturally stressed by gravitational and tectonic stress. When an underground tunnel is ex- cavated, the previous stress states existing in rock mass are disturbed, with the radial principal stress being released and tangential principal
stress concentrating in the periphery of the tunnel.7–9 In this process, the strain energy releases at some locations while accumulating at other locations, which leads to different mechanical responses of under- ground tunnels under dynamic disturbance.10,11 In addition, during the underground excavation process, the excavation damaged zone (EDZ) is formed in the proximity of the excavated tunnel. To date, considerable research efforts were devoted to investigating the formation of EDZ and the fracture mechanisms of surrounding rock during underground ex- cavations.12–17 For instance, a series of studies have been carried out at the Underground Research Laboratory (URL) since 1983 to study ex- cavation responses when underground openings were excavated.13,14 Findings of these works showed that various factors such as the near- field stress history, geological variability, excavation method, tunnel geometry, and confining pressure are responsible for the excavation damage and instability of underground openings. The presence of the EDZ around an underground opening in turn has a great influence on the mechanical, hydraulic, and thermal characteristics of surrounding rock masses. However, previous research works on the instability of
A theoretical formulation was first established to evaluate the dynamic stress concentration factor (DSCF) around a circular opening under conditions of blasting stress wave incidence. A two-dimensional numerical model was then constructed by the particle flow code (PFC) in order to simulate the dynamic responses around an underground tunnel subjected to blasting load. In the simulation, a series of horizontal blasting stress waves were applied to an underground tunnel under various in situ stress states, and then the dynamic responses around the tunnel were analyzed from the viewpoint of the dynamic stress concentration and energy evolution. The results of theoretical analysis indicated that obvious dynamic effects occur at tunnel boundary during blasting stress wave incidence, and the DSCF at the roof and floor of the tunnel is much larger than that at two sidewalls when blasting stress wave was applied to left model boundary. The numerical results showed that high static compressive stress concentration around the underground tunnel results in the accumulation of substantial strain energy at the same location. The roof and floor of the tunnel are more prone to dynamic failures during the blasting loading process. In addition, the analysis of energy dissipation indicated that the strain energy reduction and the residual kinetic energy are positively related to the lateral pressure coefficient and the burial depth of the tunnel, and the residual kinetic energy is much larger than the strain energy reduction under the same condition. Furthermore, for an underground tunnel subjected to high in situ stress, the blasting stress wave with lower amplitude is sufficient to trigger severe dynamic failures.
⁎ Corresponding author at: School of Resources and Safety Engineering, Central South University, Changsha, Hunan, China. E-mail address: lcj2015@csu.edu.cn (C. Li).
https://doi.org/10.1016/j.ijrmms.2018.02.018
Received 12 April 2017; Received in revised form 30 November 2017; Accepted 7 February 2018
Available online 28 February 2018
1365-1609/ © 2018 Elsevier Ltd. All rights reserved.
X. Li et al.
International Journal of Rock Mechanics and Mining Sciences 104 (2018) 131–146
underground openings focus on static and quasi-static conditions, and few reports have considered the effect of dynamic disturbance. Many evidences showed that, during underground excavations, dynamic disturbance such as explosion-induced vibrations from adjacent tunnel and stress impact from neighboring rockbursts have a significant in- fluence on existing tunnels.18–20 Therefore, the dynamic disturbance is an important factor to be considered when studying the stability of deep-buried tunnel.
The drill and blast (D&B) method is extensively used in mining and tunneling engineering, because it is still an economical and efficient excavation approach for rock fracture and fragmentation.21 When the drill and blast method is used in underground excavations, the blasting vibration is generated and propagates to the deep of surrounding rock mass in the form of stress waves, which may cause damage to not only the surrounding rock mass but also nearby structures.22,23 Therefore, many researchers have conducted a lot of studies on the dynamic re- sponses of structures subjected to blast-induced stress waves, aiming at putting forward more reasonable and effective support schemes. For instance, Malmgren and Nordlund24 analyzed the dynamic behaviors of shotcrete supported rock wedges subjected to blast-induced vibrations based on field measurement data in the Kiirunavaara mine, and in- dicated that a wedge can be ejected by a dynamic load even if the static safety factor is larger than 10. Therefore, the support system was sug- gested to be able to consume energy in order to support the rock wedges subjected to blasting loads. In addition, Dhakal and Pan25 carried out numerical parametric analyses to investigate the response character- istics of structures subjected to blasting-induced ground motion char- acterized by short duration, large amplitude and high frequency. Chen et al.26 theoretically investigated the dynamic responses of under- ground arch structures subjected to conventional blasting loads, with considering the influence of the curvature of structure surface and the arrival time difference. Their results indicated that the protective structures are better to be constructed in a site with smaller acoustic impedance and larger attenuation factor. Moreover, Mitelman and Elmo27 pointed out that when the blast-induced stress waves arrive at the tunnel boundary, they are reflected and converted into tensile stress waves, which can cause rock fragments to fly into the tunnel (i.e. spall failure). Based on the modeling results, the authors proposed a new approach for tunnel support designs to withstand spalling induced by blasting loads.
With rapid development of computer technology, numerical simu- lation techniques have become economical and powerful tools for modeling rock mechanics and rock engineering.28,29 Using numerical analysis methods, many researchers have carried out various studies on the dynamic response of rock mass and underground structures under dynamic disturbance. The boundary element method (BEM) was used by Stamos and Beskos30 to determine the dynamic response of large three-dimensional underground structures subjected to dynamic loads or seismic waves. Wang et al.31 analyzed dynamic fracture behaviors of rock in tension due to blast loading using a finite element method (FEM) code LS-DYNA. Ning et al.32,33 implemented the discontinuous deformation analysis (DDA) to simulate the rock mass failures by the blast-induced high pressure expansion. In their numerical model, the whole process of the blast chamber expansion, explosion gas penetra- tion, rock mass failure and cast, and the formation of the final blasting pile can be wholly reproduced. In addition, the finite difference method (FDM) based program FLAC3D was used by Wang et al.34 to study the dynamic response of underground gas storage salt cavern under seismic loads. As for dynamic responses of underground tunnels under dynamic disturbance, Zhu et al.35 used a finite element code RFPA-Dynamics to simulate the rockburst of underground opening triggered by dynamic disturbance; Li and Weng36 investigated the fracturing behaviors of deep-buried opening subjected to dynamic disturbance using LS-DYNA; Wang and Cai37 used the SPECFEM2D, a software package based on the spectral element method (SEM), to study the effect of the wavelength- to-excavation span ratio on ground motion near excavation boundaries
induced by seismic waves. By making full use of the advantages of various numerical methods, the hybrid finite-discrete element method (FEM-DEM) becomes an alternative numerical method to model blast- induced crack evolution and stress wave propagation in rock mass.38,39 Recently, the discrete element method (DEM) code PFC2D was used to study the dynamic features of stress wave propagating through rock joints40,41 and to simulate the excavation unloading process of under- ground tunnel in high stress rock mass.42,43 These works validated the feasibility and accuracy of PFC2D to simulate the dynamic process of rock. However, few numerical models were established by PFC2D to simulate the dynamic failure characteristics of underground tunnels subjected to blasting load.
This paper reports on an investigation into the dynamic stress concentration and energy evolution of an underground tunnel during blasting loading process. A two-dimensional mathematical physical model with a circular hole was first established to determine the dy- namic stress concentration factor (DSCF) around tunnel boundary under blasting stress wave incidence. Using the theoretical formula- tions, DSCF and dynamic effects at tunnel boundary was obtained. Then a two-dimensional numerical simulation model established by PFC2D was introduced to verify against the theoretical solution. Based on the numerical model, the distributions of tangential stress and strain energy around a circular tunnel under different in situ stress states were first discussed, and then parametric analyses were carried out to investigate the evolution and dissipation of strain energy around an underground tunnel under different in situ stress environments and various wave- forms of blasting stress wave. Findings of the present study indicated that dynamic disturbance and high in situ stress are two important factors to trigger dynamic failure around an underground tunnel. This paper provides an insight into the mechanism of rockbursts in the periphery of underground tunnels, as well as guidance for the design and support of deep-buried tunnels.
2. Theoretical formulation of the dynamic response of a circular hole
2.1. Dynamic response behaviors of circular hole under harmonic wave incidence
In theoretical analysis, it is assumed that a circular tunnel is ex- cavated along the direction parallel to the principal stress, so the pro- blem can be approximately regarded as a plane strain case. For an underground tunnel subjected to dynamic stress waves, according to the superposition principle,44 the stress, displacement and velocity components of the rock mass around the tunnel can be obtained by superimposing the static component induced by in situ stress with the dynamic component induced by incident plane wave under unstressed condition. However, due to the stress and deformation induced by in situ stress are time-independent, the dynamic stress wave is only con- sidered when theoretically investigating the dynamic response of un- derground tunnel. In the view of wave mechanics, the problem of the interactions between stress wave and underground opening can be re- gard as the initial-boundary value problem of wave equation. In this section, we focus on the dynamic responses of underground tunnel subjected to blasting load, which can be simplified as an analysis of circular hole subjected to a plane P wave as shown in Fig. 1, where x and y are the Cartesian coordinate system, θ and r are the Polar co- ordinate system, and a is the radius of tunnel. As the transient response induced by any form of transient loading can be determined by su- perposing harmonic waves of all frequencies, it is necessary to first determine a theoretical formulation under harmonic wave excitation, which was described in detail by Mow and Pao.45
As shown in Fig. 1, a harmonically time-varying incident plane P wave propagates along the positive direction of axis x, and the incident wave can be expressed as:
132
X. Li et al.
International Journal of Rock Mechanics and Mining Sciences 104 (2018) 131–146
(11)
2μ ∞
σ = ∑ (ε inφ ε(1) + A ε(3) + B ε(3))sin(nθ)e−iωt
rθ r2 n 041 n41 n42 n=0
2μ ∞
σ = ∑ (ε inφ ε(1) + A ε(3) + B ε(3))cos(nθ)e−iωt
θθ r2 n 021 n21 n22 n=0
(12) where ε(1), ε(3), ε(1)⋯ etc. are defined as a part of the contributions to
11 11 21
Fig. 1. Simplified model of interactions between incident P wave and underground tunnel.
circular frequency, and cp is the P wave velocity.
In terms of the wave function expansion method, the incident wave
function can be expanded as:
the stresses due to various waves, and the superscripts represent the type of Bessel function.45
The boundary condition at r=a is σrr r=a = 0, σrθ r=a = 0, thus An and Bn can be obtained from the boundary condition:
(i) i(αx−ωt) φ =φ0e
E (3) 42
(1) where α = ω/cp is the P wave number, φ0 is the amplitude, ω is the
11
E (3) 41
(3) E11
E (1) 11
E (1)
A =−εinφ 41
n n0E(3)E(3)
E (3) 12
12 E (3)
42 (13) (1)
∞
φ(i) = φ ∑ ε inJ (αr)cos(nθ)e−iωt
E (1) n n0E(3)E(3)
E (3) B=−εinφ 41 41
E11
0nn
n=0 (2)
11
E (3) 41
The stress field around the tunnel can be determined once the coefficients An and Bn are known. In this paper, we are interested in the dynamic responses at tunnel boundary, in which the tangential stress only exists. By substituting Eqs. (13) and (14) into Eq. (12), and letting r = a, we obtain the tangential stress at tunnel boundary:
⎧1 n=0 where Jn (x) is the first type of Bessel function, and εn = ⎨⎩2 n ≥ 1.
When an incident plane P-wave propagates through a circular hole, a compressional wave (P wave) and a shear wave (SV wave) arise from the circular hole boundary, because the reflecting surface is not per- pendicular to the direction of P wave incidence. The SH wave is not generated because it causes rock particles to oscillate perpendicular to the analyzed plane.44 P and SV waves can be expressed as:
12 E (3)
42 (14) where E (3)⋯ is the value of ε (3)⋯ evaluated at r = a.
∞
φ(r) = ∑ A H(1) (αr)cos(nθ)e−iωt
41∞
σθθ |r=a = μβ2φ ⎛1 − ⎞ ∑ εn in+1sn cos(nθ)e−iωt
11 11
nn n=0
(3)
π 0 ⎝ κ2 ⎠ n=0 (15) whereκistheratioofPtoSwavevelocity,κ= cp = β = λ+2μ,and
∞
ψ(r) = ∑ B H(1) (βr)sin(nθ)e−iωt
cs α μ
Dn (16)
nn
n=0 (4)
Nn
N = (n2 − 1)βaH(1) (βa) − (n3 − n + 1β2a2)H(1) (βa)
D = αaH(1) (αa)N − H(1) (αa)⎡(n3 − n + 1β2a2)βaH(1) (βa) n n−1 n n ⎢⎣ 2 n−1
−(n2 + n − 1β2a2)H(1) (βa)⎤ 4 n ⎥⎦
where φ(r) and ψ(r) are the reflected P wave and the reflected S wave, respectively, which represent waves diverging from the origin, β = ω/cs is the S wave number, cs is the S wave velocity, Hn(1) (x) is the first type of Hankel function, An and Bn are coefficients of the expressions that can be determined from the appropriate boundary conditions.
The total wave in rock mass can be obtained by adding the reflected wave to the incident wave:
(5)
Sn=
n n−1 2n (17)
∞
φ = φ(i) + φ(r) = ∑ [φ ε inJ (αr) + A H(1) (αr)cos(nθ)]e−iωt
(18)
0nn nn n=0
When an incident P-wave propagates in the intact rock medium, the stress intensity of the P wave in the propagation direction is given by σ0 = μβ2φ0. This can serve as the normalized factor, and then the dy- namic stress concentration factor (DSCF) at tunnel boundary can be defined as the ratio of σθθ to σ0:
∞
ψ = ψ(r) = ∑ B H(1) (βr)sin(nθ)e−iωt
n=0
nn
Radial stress σrr, tangential stress σθθ and shear stress σrθ can be
described in terms of displacement potential:
σrr=λ∇2φ+2μ⎡∂2φ+ ∂⎛1∂ψ⎞⎤ ⎢⎣ ∂ r 2 ∂ r ⎝ r ∂ θ ⎠ ⎥⎦
∑ n n+1 n 0 n=0
−iωt
(6)
σ41∞ σ|=σ=π⎛⎝1−κ2⎞⎠εiscos(nθ)e
θθ r=a θθ
2.2. Dynamic responses of a circular hole under transient wave incidence
The steady-state responses of a circular tunnel under a harmonic P- wave have been obtained in Section 2.1. In practice, we are more in- terested in the transient responses of the tunnel under an aperiodic disturbance such as blasting load. In order to obtain the transient re- sponses of the tunnel induced by blasting load, it is necessary to first determine the blast load variation applied to the model boundary.
In blasting process, the explosion-induced load variation is ex- tremely complex in time domain, especially if several deferred-time detonation segments are adopted in the excavations. Therefore, it is necessary to seek a relatively simple equivalent blast loading curve in
(19)
2 ⎡1∂φ 1∂2ψ ⎛1∂ψ σθθ=λ∇φ+2μ⎢r∂r+r2∂θ2 +⎜r2∂θ−r2∂r∂θ⎟⎥
(7)
(8)
(9)
1 ∂2ψ⎞⎤ ⎣ ⎝ ⎠⎦
σ = μ⎧2⎡1 ∂2φ − 1 ∂φ⎤ + ⎡ 1 ∂2ψ − r ∂ ⎛1 ∂ψ⎞⎤⎫ r θ ⎨⎩ ⎢⎣ r ∂ r ∂ θ r 2 ∂ θ ⎥⎦ ⎢⎣ r 2 ∂ θ 2 ∂ r ⎝ r ∂ r ⎠ ⎥⎦ ⎬⎭
where λ and μ are Lamé constants.
By substituting Eqs. (5) and (6) into Eqs. (7)–(9), we obtain:
2μ ∞
σ = ∑ (ε inφ ε(1) + A ε(3) + B ε(3))cos(nθ)e−iωt
rr r2 n 0 11 n 11 n 12 (10) n=0
133
X. Li et al.
International Journal of Rock Mechanics and Mining Sciences 104 (2018) 131–146
theoretical and numerical analysis. Based on previous publications re- lating to blasting procedures,11,36,46 the blasting load can be simplified as a triangular load. The entire blasting processes can be reduced to linear loading and linear unloading process, which can be expressed as:
alternatively in terms of sine or cosine transforms. When the input is an impulse function, the impulse response can be expressed by:
2∞
gδ(xi,t)= ∫ R(ω)cosωtdω
⎧0 ,t<0
π0 (26) Then the response due to a Heaviside step input can be obtained
⎪ t Pbm , 0 ≤ t < tr Pb(t)= tr
from the integral of the impulse response:
⎨ts−tPbm ,tr≤t
(23) where Re ζ = αa. By substituting Eqs. (19) and (23) into Eq. (22), we
⎪0 ⎩
, t ≥ ts
When 0 ≤ t < tr, we have:
t 1 2 ∞ R(ω)sinω(t−τ) g(xi, t) = ∫ t dτπ ∫ ω
0r0
= 2 ∫∞ R(ω)(1− cosωt)dω
πtr 0 ω2 Whentr ≤t