# Blasting Vibration Safety Criterion Analysis with Equivalent Elastic Boundary: Based on Accurate Loading Model

Authors: Qingwen Li; Lan Qiao; Gautam Dasgupta; Siwei Ma; Liping Wang; Jianghui Dong
Journal: Shock and Vibration
 (2015)

---

## Abstract

In the tunnel and underground space engineering, the blasting wave will attenuate from shock wave to stress
wave to elastic seismic wave in the host rock. Also, the host rock will form crushed zone, fractured zone, and elastic seismic
zone under the blasting loading and waves. In this paper, an accurate mathematical dynamic loading model was built. And the
crushed zone as well as fractured zone was considered as the blasting vibration source thus deducting the partial energy for cutting host rock.
So this complicated dynamic problem of segmented differential blasting was regarded as an equivalent elastic boundary problem by taking advantage
of Saint-Venant’s Theorem. At last, a 3D model in finite element software FLAC3D accepted the constitutive parameters, uniformly distributed
mutative loading, and the cylindrical attenuation law to predict the velocity curves and effective tensile curves for calculating safety criterion formulas of
surrounding rock and tunnel liner after verifying well with the in situ monitoring data.


---

## Body

## 1. Introduction

With the rapid infrastructure development in China, the improvement of control level on the blasting technique, and the well adaptability in various engineering, the blasting construction are used more and more widely; also, the drill and blast method has been becoming the most primary method in tunneling and underground space construction; in China, there are more than 95% mountain tunnel constructions which adopt the drill and blast method [1]. Also, the blasting is producing huge help on the engineering yet generating many adverse effects on existing structures. Also, many scholars have done a lot of researches on adverse effects [2–4]. Hence, predicting the vibration safety criterion has become an important and effective method for analyzing the adverse effects. Jiang and Zhou used the safety criterion method to analyze tunnel liner structure [5]. Karadogan et al. gave a damage criteria norm for blast-induced ground vibrations in Turkey [6]. Yang et al. studied the safety distance for secondary shotcrete in Jinping-II deep-buried tunnels [7]. Though many works have been done by all the world scholars, the most attempts were based on the monitoring data to summarize the empirical equation [8–18], which is restricted for the real engineering blasting vibration. At the same time, some new methods are used to analyze the blasting problem, such as the soft computing method [19], artificial neural networks method [20], MEMS-based commutation module, [21] and RES-based model [22]. But it should be based on large number of in situ monitoring data, not considering the effect of vibration frequency and duration.The numerical simulation has become an important prediction method in the engineering, due to its high accuracy and less cost, as well as the fast development computer technique, to make many scholar use such ads ANSYS/LS-DYNA finite element software to solve the blasting vibration problems [23–29]. However, most models accepted the simplified triangle dynamic load, trapezoid dynamic load, or simplified data table which are not accurate enough for the blasting loading.In blasting, the gas pressure in hole reached 50 to 100 thousand times atmosphere. The high gas pressure could crush the host rock around blast hole to form the crushed zone. And because of higher dynamic compressive strength, the crushed zone will consume large partial energy to make the shock wave attenuate very fast. After propagating a distance, the pressure cannot crush the host rock any more, to make the shock wave attenuate as stress wave. Then the surrounding rock would suffer radial compression to generate a strong radial movement to lead the rock shell expansion to cause the circular tensile stress. Also, the tensile stress tends to be greater than critical dynamic tensile strength which is usually less than 1/10 compression strength of rock masses, to lead to the radial fractures formed. After those radial fractures connected with crushed zone, the blasting gas which looked like a wedge entered the radial fractures to make them enlarge. Hence, the fractured zone would be formed which also could consume some energy to make the vibration wave attenuate as elastic seismic wave whose impact region was the elastic seismic zone.In this paper, the accurate mathematical model of blasting loading was established on the analysis of the blasting pressure change, the blast hole volume expansion, the fracture development, and the blasting gas motion. All the blasting energy transformed to the kinetic energy of the host rock without considering energy loss, this energy was used to form the crushed zone and fractured zone, and the residual energy was used to produce the vibration in elastic seismic zone to propagate out. So the crushed zone and fractured zone were also regarded as the blasting vibration source by taking advantage of the unified mechanical and continuum damage mechanics thus deducting this partial energy for cutting host rock, and regarding the vibration wave on the elastic seismic zone as the cylindrical elastic seismic wave to attenuate along with radius direction [30].For the far field vibration analysis, the separate blasting loading of each hole could meet Saint-Venant’s theorem to be regarded as the uniformly distributed loading on its own segmentation cutting zone. Hence the complicated dynamic problem in segmented differential blasting turned into an initial equivalent elastic boundary problem alongside the uniformly distributed and the attenuation law. At last, a 3D model in finite element software FLAC3D accepted the parameters to predict the velocity curve and effective tensile curve for confirming safety criterion after verifying well with the in situ monitoring data.


## 2. Equivalent Boundary of Elastic Seismic Loading

Based on above analysis, for the far field vibration, after a distance propagating, the blasting wave would attenuate as the elastic seismic wave, and this elastic seismic wave followed a circular attenuation law along the radius direction. Also, the crushed zone and fractured zone were considered as blasting vibration source boundary by taking advantage of continuum damage mechanics theory and unified constitutive theory [30].


### 2.1. Seismic Elastic Velocity

During the blasting, assuming the host rock as non-compressible linear media, and all the blasting energy would transform to kinetic energy of host rock without any energy loss. So, according to dynamic gradient theory in the semi-infinite medium, the formula of the velocity can be shown as [31](1)V=ρ0q8ρrvsL-r-r-2+L-/22,where (2)vs=ln⁡L-+1+L-2-L-+1+L-2,L-=L02r0,r-=r2r0,where V is the vibration rounded velocity, ρ0 is the explosive density, q is the blasting energy, L0 is the charging length, r0 is radius of blast hole, and r is the propagated distance.For a single blast hole in semi-infinite medium, the mainly destructional form of the crushed zone with the radiusr1 and for the fractured zone with radius r2, so the formula can be expressed as [32] (3)r1=ρrcp25σc1/2Pσ*1/4r0,r2=υP1-υσt1/αr0,where υ is Poisson’s ratio, cp is the P-wave of rock mass, σc and σt are, respectively, the dynamic compressive strength and dynamic tensile strength of rock, σ* is the static compressive strength, P is the blasting loading on the blast hole wall, and α is the attenuation coefficient of blasting vibration wave.In above,(4)cp2=λ+2Gρr,λ=2Gυ1-2υ,where G is the shear modulus and ρr is the density of rock mass.Experimental study showed that the radius of crushed zone is 3 to 5 times blast hole radius and 10 to 15 times for fractured zone [33].


### 2.2. Seismic Elastic Loading on Equivalent Boundary

For the elastic seismic loading, the attenuation functions of shock wave and stress wave with the time and displacement, to make the blasting loadingP(x,t) worked on the blast hole wall was shown as [34](5)Pax,t=Px,tr0rα,where Pa(x,t) is the blasting loading after attenuating along with distance in radius and α is the attenuation coefficient of vibration wave; in above formula, for the shock wave (6)α=2+υ1-υ,andforstresswaveα=2-υ1-υ.For the underground space and tunneling blasting construction, there are, respectively, the initial cut hole, destruction rock hole, peripheral hole, and bottom hole arranged on the cutting surface, among, the initial cut holes were the inclined hole located in the center of cut surface, and the other holes were vertical hole surrounded the initial cut holes. Hence the blasting vibration must be a multiple and superimposed waves on cutting surface.(1) The Equivalent Blasting Loading of Initial Cut Zone. According to the above formula, for elastic dynamic loading of single blast hole, for the initial cut hole with ignoring the overlapping influence and energy loss, the equivalent elastic loading can be shown as(7)Pesx,t=Px,tr0r12+υ/1-υr1r22-υ/1-υ,where Pes(x,t) is dynamic load of single initial cut hole.Assuming the angle between initial cut hole and cut surface isθ, so the projection of blast hole on the cut surface is ellipse with long radius re and the short radius r0 (the radius of blast hole), as shown in Figure 1.Figure 1
Initial cut hole with cutting surface.From the above relationship, we can knowre=r0/sinθ. Usually, the initial cut hole was the first blasting segment, so the equivalent elastic dynamic load of initial cut section with considering the crushed zone and fractured zone as dynamic loading boundary was shown in Figure 2.Figure 2
Equivalent boundary of initial cut zone.The equivalent elastic loading of all the initial cut holes by taking advantage of Saint-Venant’s principle was shown as(8)Pex,t=k1SeSeqPx,tr0r12+υ/1-υr1r22-υ/1-υ,where Pe(x,t) is equivalent elastic loading of all the initial cut holes, k1 is the number of holes, Se is the area of a single elliptic fractured zone with taking the 10 to 15 times blast hole’s projection radius, and Seq is the area of whole initial cut section.(2) The Equivalent Blasting Loading of Other Segments. Except the initial cut holes, there was a free face due to last segmented blasting for other segments. And the free face make this segmented blasting could cut the whole rock mass between this segment and last segment, as shown in Figure 3. So taking r2=lk, the formula is carried out as below:(9)Pex,t=S0S-SeqPx,tr0r12+υ/(1-υ)r1lk2-υ/(1-υ),where S0 is the area of a segmented blasting, S is the area of whole cut surface, and lk is width of kth segmented blasting.Figure 3
Equivalent boundary of other segments.Based on the time interval of two segments, the total equivalent elastic loading of other segments was shown as(10)Ptx,t=∑m=1k=nS0S-SeqPx,t+τkr0r12+υ/(1-υ)r1lk2-υ/(1-υ),where n is the number of segments and τk is the time interval between two segments.


## 2.1. Seismic Elastic Velocity

During the blasting, assuming the host rock as non-compressible linear media, and all the blasting energy would transform to kinetic energy of host rock without any energy loss. So, according to dynamic gradient theory in the semi-infinite medium, the formula of the velocity can be shown as [31](1)V=ρ0q8ρrvsL-r-r-2+L-/22,where (2)vs=ln⁡L-+1+L-2-L-+1+L-2,L-=L02r0,r-=r2r0,where V is the vibration rounded velocity, ρ0 is the explosive density, q is the blasting energy, L0 is the charging length, r0 is radius of blast hole, and r is the propagated distance.For a single blast hole in semi-infinite medium, the mainly destructional form of the crushed zone with the radiusr1 and for the fractured zone with radius r2, so the formula can be expressed as [32] (3)r1=ρrcp25σc1/2Pσ*1/4r0,r2=υP1-υσt1/αr0,where υ is Poisson’s ratio, cp is the P-wave of rock mass, σc and σt are, respectively, the dynamic compressive strength and dynamic tensile strength of rock, σ* is the static compressive strength, P is the blasting loading on the blast hole wall, and α is the attenuation coefficient of blasting vibration wave.In above,(4)cp2=λ+2Gρr,λ=2Gυ1-2υ,where G is the shear modulus and ρr is the density of rock mass.Experimental study showed that the radius of crushed zone is 3 to 5 times blast hole radius and 10 to 15 times for fractured zone [33].


## 2.2. Seismic Elastic Loading on Equivalent Boundary

For the elastic seismic loading, the attenuation functions of shock wave and stress wave with the time and displacement, to make the blasting loadingP(x,t) worked on the blast hole wall was shown as [34](5)Pax,t=Px,tr0rα,where Pa(x,t) is the blasting loading after attenuating along with distance in radius and α is the attenuation coefficient of vibration wave; in above formula, for the shock wave (6)α=2+υ1-υ,andforstresswaveα=2-υ1-υ.For the underground space and tunneling blasting construction, there are, respectively, the initial cut hole, destruction rock hole, peripheral hole, and bottom hole arranged on the cutting surface, among, the initial cut holes were the inclined hole located in the center of cut surface, and the other holes were vertical hole surrounded the initial cut holes. Hence the blasting vibration must be a multiple and superimposed waves on cutting surface.(1) The Equivalent Blasting Loading of Initial Cut Zone. According to the above formula, for elastic dynamic loading of single blast hole, for the initial cut hole with ignoring the overlapping influence and energy loss, the equivalent elastic loading can be shown as(7)Pesx,t=Px,tr0r12+υ/1-υr1r22-υ/1-υ,where Pes(x,t) is dynamic load of single initial cut hole.Assuming the angle between initial cut hole and cut surface isθ, so the projection of blast hole on the cut surface is ellipse with long radius re and the short radius r0 (the radius of blast hole), as shown in Figure 1.Figure 1
Initial cut hole with cutting surface.From the above relationship, we can knowre=r0/sinθ. Usually, the initial cut hole was the first blasting segment, so the equivalent elastic dynamic load of initial cut section with considering the crushed zone and fractured zone as dynamic loading boundary was shown in Figure 2.Figure 2
Equivalent boundary of initial cut zone.The equivalent elastic loading of all the initial cut holes by taking advantage of Saint-Venant’s principle was shown as(8)Pex,t=k1SeSeqPx,tr0r12+υ/1-υr1r22-υ/1-υ,where Pe(x,t) is equivalent elastic loading of all the initial cut holes, k1 is the number of holes, Se is the area of a single elliptic fractured zone with taking the 10 to 15 times blast hole’s projection radius, and Seq is the area of whole initial cut section.(2) The Equivalent Blasting Loading of Other Segments. Except the initial cut holes, there was a free face due to last segmented blasting for other segments. And the free face make this segmented blasting could cut the whole rock mass between this segment and last segment, as shown in Figure 3. So taking r2=lk, the formula is carried out as below:(9)Pex,t=S0S-SeqPx,tr0r12+υ/(1-υ)r1lk2-υ/(1-υ),where S0 is the area of a segmented blasting, S is the area of whole cut surface, and lk is width of kth segmented blasting.Figure 3
Equivalent boundary of other segments.Based on the time interval of two segments, the total equivalent elastic loading of other segments was shown as(10)Ptx,t=∑m=1k=nS0S-SeqPx,t+τkr0r12+υ/(1-υ)r1lk2-υ/(1-υ),where n is the number of segments and τk is the time interval between two segments.


## 3. Continuum Damage Mechanics

For the host rock of tunnel, there was a damage zone due to internal stress redistribution and fracture extension. This damage made constitutive parameter become lower, which also generated impact on vibration wave propagation. The Canadian scholars measured that the main fracture extension depth is 0.5 m, and the maximum depth could reach 1.0 m with using in situ ultrasonic wave velocity testing technology [35]. So according to continuum damage mechanics theory, the effective stress under damage condition can be expressed as [36](11)σij*=σij′1-D,where σij* is the effective stress under damage condition, σij′ is effective stress of protolith, D is the scalar of damage condition, which can be carried out by Mazar damage model with two scalars being, respectively, tensile damage Dt and compression damage Dc, shown as [37](12a)Dcεc=1-e-αcεc-εc0/εc0,(12b)Dtεt=1-e-αtεt-εt0/εt0,where αc and αt are, respectively, compression and tensile coefficients, which depended on the material properties, for the rock, their value is usually 0.5, and εc0 and εt0 are, respectively, uniaxial compression and uniaxial tensile critical strain value. As shown below:(13a)εc=∑i=1,3εi-2,(13b)εt=∑i=1,3εi+2,εi-, εi+ are, respectively, the negative and positive strains.So, the damage scalar is the superposition of compression and tensile:(14)D=AcDc+AtDt,Ac+At=1,(15a)Ac=∑i=1,3Hiεi-εi-+εi+ε-i2,(15b)At=∑i=1,3Hiεi+εi-+εi+ε-i2,where ε-i2=∑i=1,3εi-+εi+2 is the effective strain, and when x<0, Hx=0, and when x>0, Hx=x.Based on above analysis, in the damage host rock, the formula of Young modulus between damage rock and original rock is(16)Ed=E01-D,where Ed is the Young modulus under damage rock and E0 is Young modulus under original rock.


## 4. Equivalent Process of Blasting Loading

After blasting, the dynamic loading on the blast hole wall can make the volume of blast hole enlarged and the fracture expanded; also, the gas pressure and dynamic loading must be reduced with volume enlargement. At last, the explosion gas rapidly overflowed and the applied force decayed to zero with fracture development to connect each other.Based on blasting mechanism, the process of blasting loading can be divided into four stages.(1)
The first stage: the dynamic loading will increase with time till reaching the peak intensity of blasting.(2)
The second stage: the blasting pressure will be reduced by the fracture expending, the fillings moving, and the volume increasing before the filling was ejected from the hole.(3)
The third stage: the explosive gas erupts quickly from the blast hole to lead to lower pressure after the filling was ejected.(4)
The fourth stage: the explosive gas rapidly overflowed and the applied force decayed to zero when fractures develop to connect together.


### 4.1. The First Stage of Blasting

The dynamic loading will increase with time till reaching the peak intensity of blasting when the detonation gas wave propagated to the bottom section of blast hole after exploding. Many research works show that the initial peak blasting loading must have relationship to detonation wave pressure, and according to the Chapman-Jouguet model, the detonation wave pressure in an explosion can be guided by the widely known equation [34](17)PD=1γ+1ρ0VD,where PD is the largest detonation wave pressure, ρ0 is the density of explosive, VD is detonation velocity, and γ is the ratio of the specific heats for the detonation gases; in this formula γ=3.0.The initial explosion pressure which was the explosion gas acted on the blast hole wall just after detonation is approximately the half of the detonation pressure for the coupled charge:(18)P1=ρ0VD2(γ+1).For decoupled charges, the initial explosion pressure also had the relationship to the proportion between the blast hole diameter and the charge diameter, so the formula is(19)P1=ρ0VD2(γ+1)ab2γ,where a is the charge diameter and b is the blast hole diameter.The rising time of loading was shown as(20)t1=LVD,where L is the length of blast hole.


### 4.2. The Second Stage of Blasting

The blasting pressure will reduce with the fracture expending, the filling moving, and the volume increasing of blast hole before the filling was ejected from the hole. So the volume of detonation gas with time can be expressed as(21)ΔVt=2πrtutLdt+2L∫0Laωηdη+14πrt2yt,where r(t) is the blast hole radius with time, u(t) is the expansion velocity of blast hole wall, ω(η) is the width of fracture, and y(t) is displacement of filling with time.According to gas law, gas pressure with volume changed in detonation cavity can be shown as [5] (22)P2=A1-ωR1Ve-R1V+B1-ωR2Ve-R2V+ωE0V.Put the formula (21) into formula (22):(23)P2t=A1-ωR1V0+ΔVte-R1V0+ΔVt+B1-ωR2V0+ΔVte-R2V0+ΔVt+ωE0V0+ΔVt,where V0 is the initial volume of blast hole, A, B, R1, R2, and ω are all the explosive material parameter, and E0 is initial energy of explosive.


### 4.3. The Third Stage of Blasting

The explosive gas erupts quickly out from the blast hole to lead to pressure getting lower after the filling was ejected or no filling blasting; according to gas dynamics theory, the blast hole in this period can be looked at as a bottle structure which follows the next assumptions, (a) because of the volume of blast hole was enlarged in last stage, which make the section of exit is smaller than inside, the shape of blast hole looks like a bottle; (b) because of very short time of reaction and better heat-insulating property of rock, the whole process are assumed as adiabatic process; (c) because of the very slow speed of volume expansion during last stage, the initial velocity of gas is zero at the beginning of this stage, sov0=0, as shown in Figure 4.Figure 4
The equivalent structure of third stage.In Figure4, V0, P0, ρ0, T0, and v0 are, respectively, the initial volume, pressure, density, temperature, and velocity. Also, Pe, ρe, Te, and ve are, respectively, the pressure, density, temperature, and velocity at the section of exit.For this bottle structure, the inner pressure change of blast hole had the relationship with gas erupted statement; whenPcr/Pa=1, the exit statement was critical condition, and the velocity of gas erupted was sound velocity; when Pcr/Pa<1, the exit statement was subcritical condition, and the pressure at the exit was equal to air pressure Pe=Pa; when Pcr/Pa>1, the exit statement was supercritical statement, and there was congested phenomenon at exit to make the velocity of gas also the sound velocity, Ma=vcr/vs=1, where Pcr, vcr are, respectively, critical statement of pressure and velocity.According to isentropic gas formula,Pcr/P0=Tcr/T0γ/(γ-1)=2/(γ+1)γ/(γ-1) to determine critical pressure value Pcr, where Tcr is the critical temperature, γ is the specific heat ratio of the gas, and R=297 J/(kg·K) is the gas constant, so (24)PcrP0=0.3536.Usually, the gas pressure in blast hole can reach p0=(50000~10000)pa, so the critical pressure pcr≫pa, which was supercritical statement to make the velocity of erupted gas be sound velocity.Based on above analysis, the gas flow formula at the supercritical statement was(25)qmcr=A2γ+1(γ+1)/2(γ-1)γP0ρ0.From time t to t+dt, based on the first law of thermodynamics with taking advantage of adiabatic process of gas in blast hole,(26)Pργ=const.So, (27)Pe(t)ρetγ=pe(t+dt)ρet+dtγ.Put the formula (25) into formula (27):(28)Pe(t)ρetγ=Pe(t+dt)ρe(t)-A/V02/(γ+1)(γ+1)/2γ-1γP0ρ0dtγ.So, (29)Pe(t+dt)=Petρet-A/V2/γ+1γ+1/2γ-1γP0ρ0dtρetγ.So, the formula of gas pressure change was(30)P3t+dt=Pe(t+dt)=Pet1-AV2γ+1γ+1/2γ-1γP0ρ0ρetdtγ.


### 4.4. The Fourth Stage of Blasting

The explosive gas rapidly overflows and the applied force decays to zero fast with fractures development to connect.


### 4.5. The Discipline of Loading

According the above separated analysis of blasting process, the blasting loading curve was shown as in Figure5.Figure 5
Blasting loading model.In Figure5, t0 is the beginning of blasting time, so t0=0; the initial dynamic loading P0=0; t1 is the maximum blasting loading time, P1 is the maximum blasting loading, t2 is the time of fillings ejected, P2 is the loading on the blast hole wall at that time, t3 is the time just before fractures connect each other, P3 is the loading on the blast hole wall at that time, t4 is the end time of blasting, and P4=0 is the gas pressure.


## 4.1. The First Stage of Blasting

The dynamic loading will increase with time till reaching the peak intensity of blasting when the detonation gas wave propagated to the bottom section of blast hole after exploding. Many research works show that the initial peak blasting loading must have relationship to detonation wave pressure, and according to the Chapman-Jouguet model, the detonation wave pressure in an explosion can be guided by the widely known equation [34](17)PD=1γ+1ρ0VD,where PD is the largest detonation wave pressure, ρ0 is the density of explosive, VD is detonation velocity, and γ is the ratio of the specific heats for the detonation gases; in this formula γ=3.0.The initial explosion pressure which was the explosion gas acted on the blast hole wall just after detonation is approximately the half of the detonation pressure for the coupled charge:(18)P1=ρ0VD2(γ+1).For decoupled charges, the initial explosion pressure also had the relationship to the proportion between the blast hole diameter and the charge diameter, so the formula is(19)P1=ρ0VD2(γ+1)ab2γ,where a is the charge diameter and b is the blast hole diameter.The rising time of loading was shown as(20)t1=LVD,where L is the length of blast hole.


## 4.2. The Second Stage of Blasting

The blasting pressure will reduce with the fracture expending, the filling moving, and the volume increasing of blast hole before the filling was ejected from the hole. So the volume of detonation gas with time can be expressed as(21)ΔVt=2πrtutLdt+2L∫0Laωηdη+14πrt2yt,where r(t) is the blast hole radius with time, u(t) is the expansion velocity of blast hole wall, ω(η) is the width of fracture, and y(t) is displacement of filling with time.According to gas law, gas pressure with volume changed in detonation cavity can be shown as [5] (22)P2=A1-ωR1Ve-R1V+B1-ωR2Ve-R2V+ωE0V.Put the formula (21) into formula (22):(23)P2t=A1-ωR1V0+ΔVte-R1V0+ΔVt+B1-ωR2V0+ΔVte-R2V0+ΔVt+ωE0V0+ΔVt,where V0 is the initial volume of blast hole, A, B, R1, R2, and ω are all the explosive material parameter, and E0 is initial energy of explosive.


## 4.3. The Third Stage of Blasting

The explosive gas erupts quickly out from the blast hole to lead to pressure getting lower after the filling was ejected or no filling blasting; according to gas dynamics theory, the blast hole in this period can be looked at as a bottle structure which follows the next assumptions, (a) because of the volume of blast hole was enlarged in last stage, which make the section of exit is smaller than inside, the shape of blast hole looks like a bottle; (b) because of very short time of reaction and better heat-insulating property of rock, the whole process are assumed as adiabatic process; (c) because of the very slow speed of volume expansion during last stage, the initial velocity of gas is zero at the beginning of this stage, sov0=0, as shown in Figure 4.Figure 4
The equivalent structure of third stage.In Figure4, V0, P0, ρ0, T0, and v0 are, respectively, the initial volume, pressure, density, temperature, and velocity. Also, Pe, ρe, Te, and ve are, respectively, the pressure, density, temperature, and velocity at the section of exit.For this bottle structure, the inner pressure change of blast hole had the relationship with gas erupted statement; whenPcr/Pa=1, the exit statement was critical condition, and the velocity of gas erupted was sound velocity; when Pcr/Pa<1, the exit statement was subcritical condition, and the pressure at the exit was equal to air pressure Pe=Pa; when Pcr/Pa>1, the exit statement was supercritical statement, and there was congested phenomenon at exit to make the velocity of gas also the sound velocity, Ma=vcr/vs=1, where Pcr, vcr are, respectively, critical statement of pressure and velocity.According to isentropic gas formula,Pcr/P0=Tcr/T0γ/(γ-1)=2/(γ+1)γ/(γ-1) to determine critical pressure value Pcr, where Tcr is the critical temperature, γ is the specific heat ratio of the gas, and R=297 J/(kg·K) is the gas constant, so (24)PcrP0=0.3536.Usually, the gas pressure in blast hole can reach p0=(50000~10000)pa, so the critical pressure pcr≫pa, which was supercritical statement to make the velocity of erupted gas be sound velocity.Based on above analysis, the gas flow formula at the supercritical statement was(25)qmcr=A2γ+1(γ+1)/2(γ-1)γP0ρ0.From time t to t+dt, based on the first law of thermodynamics with taking advantage of adiabatic process of gas in blast hole,(26)Pργ=const.So, (27)Pe(t)ρetγ=pe(t+dt)ρet+dtγ.Put the formula (25) into formula (27):(28)Pe(t)ρetγ=Pe(t+dt)ρe(t)-A/V02/(γ+1)(γ+1)/2γ-1γP0ρ0dtγ.So, (29)Pe(t+dt)=Petρet-A/V2/γ+1γ+1/2γ-1γP0ρ0dtρetγ.So, the formula of gas pressure change was(30)P3t+dt=Pe(t+dt)=Pet1-AV2γ+1γ+1/2γ-1γP0ρ0ρetdtγ.


## 4.4. The Fourth Stage of Blasting

The explosive gas rapidly overflows and the applied force decays to zero fast with fractures development to connect.


## 4.5. The Discipline of Loading

According the above separated analysis of blasting process, the blasting loading curve was shown as in Figure5.Figure 5
Blasting loading model.In Figure5, t0 is the beginning of blasting time, so t0=0; the initial dynamic loading P0=0; t1 is the maximum blasting loading time, P1 is the maximum blasting loading, t2 is the time of fillings ejected, P2 is the loading on the blast hole wall at that time, t3 is the time just before fractures connect each other, P3 is the loading on the blast hole wall at that time, t4 is the end time of blasting, and P4=0 is the gas pressure.


## 5. Project of Case Study

### 5.1. Engineering Introduction

The Chen-Chi West Connection Line Highway Tunnel, which was designed as a separate double tunnel each containing two lanes, was chosen as the case to be studied. The lengths of two tunnels are 1075 m and 1185 m, and the distance between them is between 18 m and 20 m. The cavern on the foothill region is 14 m wide and 12 m high. In this tunnel blasting, No. 2 rock emulsion explosive with diameter is 32 mm, length is 200 mm, and weight is 150 g was used. And positive multistage differential blasting method was taken by initial cut holes, destruction rock holes, peripheral holes, and bottom holes on the cutting surface. In working surface the distance between the holes was from 0.5 m to 0.7 m, the depth of the blast hole was 3.0 m, the charge length was 2.5 m without fillings, and the detonating cord across through the whole hole from opening to bottom. The blasting parameters were shown in Table1, and the form of segments was shown in Figure 6 and interval time was shown in Table 2.Table 1
Blasting parameters.
Charge diameter (mm)
Blast hole diameter (mm)
Holes number
Explosive density (kg/m3)
Detonation velocity (m/s)
32
40
80
1000
3200Table 2
The time interval of every multistage.
MS1
MS3
MS5
MS7
MS9
MS11
MS13
0 ms
50 ms
110 ms
200 ms
310 ms
460 ms
650 msFigure 6
Blasting segmented situation.


### 5.2. Constitutive Parameter

The formation lithology rock in the tunnel area is conglomerate, which belongs to the late Jurassic conglomerate, color is purple and structure is gravel with the diameter from 2 mm to 60 mm, and the gravels are cemented together with fillings to form the thick block layer. The compositions of the gravels are andesite, welded tuff, rhyolite, and siliceous rocks, and for the fillings they are debris of rock, fine sand, quartz, feldspar, and some other stable mineral sand. The constitutive parameters are shown in Table3.Table 3
Constitutive parameters of the rock.
Density (Kg/m³)
Compressive strength (MPa)
Tensile strength  (MPa)
Poisson’s ratio
Cohesion (MPa)
Friction angle (MPa)
Compressive hardening coefficient (MPa)
Tensile hardening coefficient (MPa)
2820
75.5
3.0
0.2
4.5
38
0.26
0.33


## 5.1. Engineering Introduction

The Chen-Chi West Connection Line Highway Tunnel, which was designed as a separate double tunnel each containing two lanes, was chosen as the case to be studied. The lengths of two tunnels are 1075 m and 1185 m, and the distance between them is between 18 m and 20 m. The cavern on the foothill region is 14 m wide and 12 m high. In this tunnel blasting, No. 2 rock emulsion explosive with diameter is 32 mm, length is 200 mm, and weight is 150 g was used. And positive multistage differential blasting method was taken by initial cut holes, destruction rock holes, peripheral holes, and bottom holes on the cutting surface. In working surface the distance between the holes was from 0.5 m to 0.7 m, the depth of the blast hole was 3.0 m, the charge length was 2.5 m without fillings, and the detonating cord across through the whole hole from opening to bottom. The blasting parameters were shown in Table1, and the form of segments was shown in Figure 6 and interval time was shown in Table 2.Table 1
Blasting parameters.
Charge diameter (mm)
Blast hole diameter (mm)
Holes number
Explosive density (kg/m3)
Detonation velocity (m/s)
32
40
80
1000
3200Table 2
The time interval of every multistage.
MS1
MS3
MS5
MS7
MS9
MS11
MS13
0 ms
50 ms
110 ms
200 ms
310 ms
460 ms
650 msFigure 6
Blasting segmented situation.


## 5.2. Constitutive Parameter

The formation lithology rock in the tunnel area is conglomerate, which belongs to the late Jurassic conglomerate, color is purple and structure is gravel with the diameter from 2 mm to 60 mm, and the gravels are cemented together with fillings to form the thick block layer. The compositions of the gravels are andesite, welded tuff, rhyolite, and siliceous rocks, and for the fillings they are debris of rock, fine sand, quartz, feldspar, and some other stable mineral sand. The constitutive parameters are shown in Table3.Table 3
Constitutive parameters of the rock.
Density (Kg/m³)
Compressive strength (MPa)
Tensile strength  (MPa)
Poisson’s ratio
Cohesion (MPa)
Friction angle (MPa)
Compressive hardening coefficient (MPa)
Tensile hardening coefficient (MPa)
2820
75.5
3.0
0.2
4.5
38
0.26
0.33


## 6. Numerical Simulation Analysis

### 6.1. Simulation Model

Based on the actual engineering geology and structures, a 3D model was built in the finite element software of FLAC3D to accept the constitutive parameters, blasting load, and wave propagation parameters, as shown in Figure7. The geometrical dimensions are length 228 m (X direction), width 203 m (Y direction), and height 159 m (Z direction). The finite element mesh consists of 386527 nodes and 372200 elements. And for the blasting zone, the group was between 100 m and 103 m (Y direction) after the cutting length 100 m, the initial cut hole zone was meshed as 192 elements (brownish zone) and other worked zone was meshed as 606 elements (orange zone), also, as shown in Figure 7.Figure 7
3D model in numerical simulation.


### 6.2. Dynamic Loading in Simulation

For the decoupled charge blasting, every stage pressure and arrival time were calculated by using above formulas. The detail values were shown in Table4.Table 4
The detail value of every stage pressure and arrival time.
P
D (MPa)
P
L (MPa)
P
0 (MPa)
P
1 (MPa)
P
2 (MPa)
P
3 (MPa)
P
4 (MPa)
Initial cut zone
2560
335.5
0
18.73
18.73
2.02
0
Other zones
2560
335.5
0
34.22
34.22
6.45
0
Arrival time (s)
t
0
=
0
t
1
=
0.01
t
2
=
0.01
t
3
=
0.122
t
4
=
0.122In the numerical simulation, it is difficult to determine the dumping of the rock mass directly, so, the field monitoring and calculation is necessary; in this paper, the dumping which was input into model was 0.015 based on the field monitoring and calculation repeatedly.


### 6.3. Verification of Simulation

In this simulation, the solve time was set to 2 s, and after 1284933 calculation steps, the curves of velocity were predicted from the model compared with 9 times field monitoring data at the same location. The simulation data and monitoring data were shown in Table5, and the typical velocity curves were shown in Figures 8 and 9.Table 5
The simulation and monitoring data of blasting vibration.
Date
Total charge (kg)
Distance (m)
Measured vertical PPV (cm/s)
Predicted vertical PPV (cm/s)
Error rate (%)
2013.5.18
240
64.3
1.263
1.293
+2.3
2013.5.22
240
58.8
1.519
1.644
+8.2
2013.5.24
240
57.2
1.705
1.771
+3.8
2013.5.27
240
54.8
1.999
1.998
−0.5
2013.6.02
240
50.0
2.505
2.510
+2
2013.6.03
240
53.2
2.248
2.176
−3.2
2013.6.11
240
57.2
2.304
1.771
−23.1
2013.6.12
240
58.8
1.582
1.641
+3.7
2013.6.14
240
62.3
1.747
1.403
−19.6Figure 8
Simulated velocity curves in surrounding rock.Figure 9
Monitored velocity curves in surrounding rock.As shown in Figure8, the predicted peak velocity was 2.505 cm/s compared with monitored peak velocity 2.51 cm/s which was obtained from Figure 9, so, the error rate is 2%; also, for those two curves, they had similar arrival time from 0.3 s to 0.7 s. From the error rate and similar arrival time, the simulation result in the surrounding rock was accurate enough for analysis. Also, for the other group data in Table 4, the error rates are from 0.5% to 23.1%, which all show that the simulation had the higher accurate results.


### 6.4. Predicted Results Analysis

After the simulation, the velocity curves and effective tensile stress curves with the distances 10 m, 20 m, 30 m, 50 m, and 80 m from the blasting point in surrounding rock and tunnel liner were obtained. The typical curves were shown in Figures10–13.Figure 10
Velocity curves in surrounding rock.Figure 11
Tensile stress curves in surrounding rock.Figure 12
Velocity curves in tunnel liner.Figure 13
Tensile stress curves in tunnel liner.For the surrounding rock, as shown in Figure10, the peak velocity is 7.099. And the maximum effective tensile stress is 0.4502, which was shown in Figure 11. Also, for the tunnel liner, it is easy to learn that the peak velocity is 2.719 and maximum effective tensile stress is 0.7434 from Figures 12 and 13.


### 6.5. Conforming with Safety Criterion

The peak velocity and maximum effective stress can be obtained from the simulated model, which were shown in Table6.Table 6
Safety criterion data.
Distance (m)
Surrounding rock
Tunnel liner
Peak velocity (cm/s)
Maximum effective tensile stress (MPa)
Peak velocity (cm/s)
Maximum effective tensile stress (MPa)
10
11.33
2.337
3.162
1.050
20
7.099
0.4502
2.719
0.7434
30
3.370
0.1795
1.918
0.5172
50
2.505
0.1233
1.057
0.2677
80
1.455
0.1111
0.7643
0.1784Summarizing data in Table6, the relationships between the effective tensile stress and peak velocity curves are, respectively, the safety criterion for the surrounding rock and tunnel liner, which were shown in Figures 14 and 15.Figure 14
Safety criterion in surrounding rock.Figure 15
Safety criterion in tunnel liner.As shown in Figure14, the safety criterion formula of the surrounding rock was (31)σt=0.0363PPV2-0.25PPV+0.4878,R2=0.991.In this case study, the rock dynamic tensile strength is 4 MPa, so after calculating, when the PPV in the surrounding rock at the junction of tunnel arch and wall reaches 13.865 cm/s, the effective tensile stress will approach the maximum tensile strength.As shown in Figure15, the safety criterion formula was (32)σt=0.0721PPV2+0.0593PPV+0.1089,R2=0.988.Also, for the tunnel liner, because there was a damage zone caused by excavated stress resilience and before blasts, the damage scale was D=0.2 which was calculated by Mazar damage model based on the field fracture surveyed, so the safety criterion formula of the damage tunnel liner is (33)σt*=0.0721PPV2+0.0593PPV+0.10891-D.The damage dynamic tensile strength was 3.2 MPa with the damage scale of D=0.2, so, when the PPV in the tunnel linear at the junction of tunnel arch and wall reaches 6.213 cm/s, the effective tensile stress will approach the maximum tensile strength.


## 6.1. Simulation Model

Based on the actual engineering geology and structures, a 3D model was built in the finite element software of FLAC3D to accept the constitutive parameters, blasting load, and wave propagation parameters, as shown in Figure7. The geometrical dimensions are length 228 m (X direction), width 203 m (Y direction), and height 159 m (Z direction). The finite element mesh consists of 386527 nodes and 372200 elements. And for the blasting zone, the group was between 100 m and 103 m (Y direction) after the cutting length 100 m, the initial cut hole zone was meshed as 192 elements (brownish zone) and other worked zone was meshed as 606 elements (orange zone), also, as shown in Figure 7.Figure 7
3D model in numerical simulation.


## 6.2. Dynamic Loading in Simulation

For the decoupled charge blasting, every stage pressure and arrival time were calculated by using above formulas. The detail values were shown in Table4.Table 4
The detail value of every stage pressure and arrival time.
P
D (MPa)
P
L (MPa)
P
0 (MPa)
P
1 (MPa)
P
2 (MPa)
P
3 (MPa)
P
4 (MPa)
Initial cut zone
2560
335.5
0
18.73
18.73
2.02
0
Other zones
2560
335.5
0
34.22
34.22
6.45
0
Arrival time (s)
t
0
=
0
t
1
=
0.01
t
2
=
0.01
t
3
=
0.122
t
4
=
0.122In the numerical simulation, it is difficult to determine the dumping of the rock mass directly, so, the field monitoring and calculation is necessary; in this paper, the dumping which was input into model was 0.015 based on the field monitoring and calculation repeatedly.


## 6.3. Verification of Simulation

In this simulation, the solve time was set to 2 s, and after 1284933 calculation steps, the curves of velocity were predicted from the model compared with 9 times field monitoring data at the same location. The simulation data and monitoring data were shown in Table5, and the typical velocity curves were shown in Figures 8 and 9.Table 5
The simulation and monitoring data of blasting vibration.
Date
Total charge (kg)
Distance (m)
Measured vertical PPV (cm/s)
Predicted vertical PPV (cm/s)
Error rate (%)
2013.5.18
240
64.3
1.263
1.293
+2.3
2013.5.22
240
58.8
1.519
1.644
+8.2
2013.5.24
240
57.2
1.705
1.771
+3.8
2013.5.27
240
54.8
1.999
1.998
−0.5
2013.6.02
240
50.0
2.505
2.510
+2
2013.6.03
240
53.2
2.248
2.176
−3.2
2013.6.11
240
57.2
2.304
1.771
−23.1
2013.6.12
240
58.8
1.582
1.641
+3.7
2013.6.14
240
62.3
1.747
1.403
−19.6Figure 8
Simulated velocity curves in surrounding rock.Figure 9
Monitored velocity curves in surrounding rock.As shown in Figure8, the predicted peak velocity was 2.505 cm/s compared with monitored peak velocity 2.51 cm/s which was obtained from Figure 9, so, the error rate is 2%; also, for those two curves, they had similar arrival time from 0.3 s to 0.7 s. From the error rate and similar arrival time, the simulation result in the surrounding rock was accurate enough for analysis. Also, for the other group data in Table 4, the error rates are from 0.5% to 23.1%, which all show that the simulation had the higher accurate results.


## 6.4. Predicted Results Analysis

After the simulation, the velocity curves and effective tensile stress curves with the distances 10 m, 20 m, 30 m, 50 m, and 80 m from the blasting point in surrounding rock and tunnel liner were obtained. The typical curves were shown in Figures10–13.Figure 10
Velocity curves in surrounding rock.Figure 11
Tensile stress curves in surrounding rock.Figure 12
Velocity curves in tunnel liner.Figure 13
Tensile stress curves in tunnel liner.For the surrounding rock, as shown in Figure10, the peak velocity is 7.099. And the maximum effective tensile stress is 0.4502, which was shown in Figure 11. Also, for the tunnel liner, it is easy to learn that the peak velocity is 2.719 and maximum effective tensile stress is 0.7434 from Figures 12 and 13.


## 6.5. Conforming with Safety Criterion

The peak velocity and maximum effective stress can be obtained from the simulated model, which were shown in Table6.Table 6
Safety criterion data.
Distance (m)
Surrounding rock
Tunnel liner
Peak velocity (cm/s)
Maximum effective tensile stress (MPa)
Peak velocity (cm/s)
Maximum effective tensile stress (MPa)
10
11.33
2.337
3.162
1.050
20
7.099
0.4502
2.719
0.7434
30
3.370
0.1795
1.918
0.5172
50
2.505
0.1233
1.057
0.2677
80
1.455
0.1111
0.7643
0.1784Summarizing data in Table6, the relationships between the effective tensile stress and peak velocity curves are, respectively, the safety criterion for the surrounding rock and tunnel liner, which were shown in Figures 14 and 15.Figure 14
Safety criterion in surrounding rock.Figure 15
Safety criterion in tunnel liner.As shown in Figure14, the safety criterion formula of the surrounding rock was (31)σt=0.0363PPV2-0.25PPV+0.4878,R2=0.991.In this case study, the rock dynamic tensile strength is 4 MPa, so after calculating, when the PPV in the surrounding rock at the junction of tunnel arch and wall reaches 13.865 cm/s, the effective tensile stress will approach the maximum tensile strength.As shown in Figure15, the safety criterion formula was (32)σt=0.0721PPV2+0.0593PPV+0.1089,R2=0.988.Also, for the tunnel liner, because there was a damage zone caused by excavated stress resilience and before blasts, the damage scale was D=0.2 which was calculated by Mazar damage model based on the field fracture surveyed, so the safety criterion formula of the damage tunnel liner is (33)σt*=0.0721PPV2+0.0593PPV+0.10891-D.The damage dynamic tensile strength was 3.2 MPa with the damage scale of D=0.2, so, when the PPV in the tunnel linear at the junction of tunnel arch and wall reaches 6.213 cm/s, the effective tensile stress will approach the maximum tensile strength.


## 7. Conclusion

For the far field vibration analysis, based on the unified theory, the crushed zone and fractured zone can be regarded as blasting vibration source. Also, the wave propagation was essentially elastic. It was proved to be very accurate after comparing between the simulated velocity curves and the field monitoring data.The accurate mathematic loading model was carried out by the analysis of the blasting pressure change, blast hole volume expansion, the fracture development, and the blasting gas motion; after well verification, this dynamic loading curve was more realistic and more accurate to the real loading, than other simplified loadings.For the far field vibration analysis, the separated blasting loadings of each hole were regarded as uniformly distributed loading on its own segmented zone by taking advantage of Saint-Venant theorem; it also met the requirement of analysis by comparing the monitoring data.Based on the predicted velocity curves and effective tensile stress curves in surrounding rock and tunnel linear, the safety criterion field can be carried out to analyze the safety and stability of every point in field. And through adopting the continuum damage mechanics theory, the safety criterion of tunnel linear was more corresponding to the real situation.



---

*Source: 604683-2015-05-20.xml*