# Complex Modeling of the Effects of Blasting on the Stability of Surrounding Rocks and Embankment in Water-Conveyance Tunnels

Authors: Xian-qi Zhou; Jin Yu; Jin-bi Ye; Shi-yu Liu; Ren-guo Liao; Xiu-wen Li
Journal: Complexity
 (2018)

---

## Abstract

Blasting in water-conveyance tunnels that cross rivers is vital for the safety and stability of embankments. In this work, a tunnel project that crosses the Yellow River in the north district of the first-phase Eastern Line of the South-to-North Water Diversion Project was selected as the research object. A complex modeling and numerical simulation on embankment stability with regard to the blasting power of the tunnel was conducted using the professional finite difference software FLAC3D to disclose the relationships between the blasting seismic waves with vibration velocity and embankment displacement under different excavation steps. Calculation results demonstrated that displacement generally attenuated from the tunnel wall to the internal structure of rocks under the effect of blasting seismic waves. The tunnel wall was in tension, and tensile stress gradually transformed into compressive stress with increased depth into the rocks. The curtain-grouting zone was mainly concentrated in the compressive zones. For different excavation steps, the vibration velocity at different feature points was high at the beginning of blasting and then gradually decreased. The resultant displacement was relatively small in the early excavation period and slowly increased as blasting progressed. The effects of different excavation steps on the safety of surrounding rocks and embankment under blasting seismic waves were simulated. We found that the blasting-induced vibration velocity was within the safe range of the code and that the calculated displacement was within the allowed range. Numerical simulation was feasible to assess the safety and stability of engineering projects.


---

## Body

## 1. Introduction

Currently, borehole-blasting method is widely applied in tunnel construction. The advantages of this method are simple operation and low cost. However, blasting seismic waves influence the building structures on the Earth’s surface to a certain extent [1]. The effects of blasting seismic waves on the safety of surrounding rocks and embankment should be considered when a tunnel runs through large rivers because the embankment safety of rivers is important. To solve these problems, abundant theoretical and experimental studies have been reported. On the basis of elastic stress wave theory, Blair and Jiang [2] investigated the propagation law on blasting seismic waves during vertical column charge explosion under elasto-plastic conditions. They concluded that the surface vibration intensity in regions far from the explosive source increased with the increase of detonation velocity, whereas the surface vibration intensity in regions close to the explosive source increased to a critical value first and gradually decreased. Moreover, the surface vibration velocity was negatively related with the damping of the propagation medium. Several scholars have discussed the explosive load dynamic effect of architectural structures and the propagation law and characteristics of blasting seismic waves. Research conclusions have been widely applied in engineering practices [3–5]. Ozer investigated the influences of blasting construction in subway tunnels on ground vibration. The acquired data were analyzed through a field experiment, and the maximum vibration velocity and frequency during blasting were estimated [6]. Kuzu evaluated the effects of blasting in underground tunnels on overground architectures in urban areas and tested the degree of injury of these overground architectures under different blasting conditions through a field test [7]. Dowding et al. discussed the effects of excavation blasting of rocks on the strain of urban structures through a case study. In addition, 8 blasting tests at 10 points of the structure and foundation rocks of the structure were observed. The velocity at middle point of rocks during blasting was significantly higher than the normal value, and the strain of external wall was similar or slightly lower than the deformation of the masonry structure or materials for wall [8]. Minaev assessed the relationship between the blasting in dense sandstone and the stability of upper structure. He provided a theoretical proof and applied it in actual engineering [9]. Apart from theoretical analysis and experimental test, Valliappan and Ang analyzed the effects of blasting excavation in underground tunnels on lower architectural and overground structures. They concluded that the generated high-intensity seismic waves from blasting excavation of underground tunnels generated outstanding amplitudes when they were propagated on the surface of rocks. To analyze the effects of underground blasting on the Earth’s surface, the influences of underground tunnel blasting on other structures were simulated by finite element software [10]. Eslami and Goshtasbi investigated the blasting principle and simulated the influences of specific parameters on surface structures during blasting by using 3DEC. Therefore, a research method for the accurate prediction of surface structural damages caused by blasting in underground tunnels was proposed [11]. Nguyen et al. indicated that blasting in road constructions caused the damages of adjacent structures. In addition, he simulated blasting in tunnel engineering by using the FEM software MIDAS GTS NX and compared it with empirical observation values. On the basis of the analysis results, the propagation velocity of surface blasting vibration was related with the radius from the measuring point to the blasting source [12]. Lu et al. conducted a field test on blasting of a tunnel that ran through the EN SHI Airport. On the basis of monitoring data analysis, the natural vibration frequency of the runway was significantly different from blasting frequency. The influence mechanism was analyzed by the dynamic finite element software LS-DYNA, and the final safety threshold was defined [13]. Despite the analysis on the influences of blasting vibration in tunnels based on existing software, several researchers have developed several specific analysis tools. Xu and Deng believed that special tools were required to investigate the unsteady behavior caused by blasting seismic waves in tunnel construction. Hence, an analysis tool was developed based on wavelet theory and was used to analyze the vibration response of one large dam in China caused by blasting of adjacent tunnel [14]. Swoboda et al. assumed that a static analysis in tunnel blasting ignored the dynamic effect. The author introduced a numerical model that learned the changes of rocks caused by blasting vibrations [15]. Ning et al. analyzed the blasting vibration in jointed rocks by discontinuous deformation analysis (DDA). In DDA, the dynamic parameter adjustment and nonreflection boundary conditions were considered and the subblock DDA was applied in the software [16]. The characteristics of surrounding tunnel rocks under blasting vibrations have been investigated [17–20]. However, these studies did not consider the influences of blasting vibrations in underwater tunnels on surrounding rocks. Relevant studies on embankment safety under the effect of blasting vibrations in a tunnel that ran through the Yellow River have mainly concentrated in the safety evaluation [21–23] and characteristic test of soil mass [24, 25]. Current studies mainly focus on propagation laws of blasting seismic waves in road tunnels, subways, and effects of blasting excavation on surrounding rocks of underground tunnels by numerical simulation, field monitoring, and so on. However, only few studies have discussed the influences of blasting vibrations on the stability of surrounding rocks and embankment during blasting of the tunnel. In this study, a numerical simulation based on the blasting in a tunnel that runs through the Yellow River was conducted to investigate the influences of blasting on the stability of surrounding rocks and embankment.


## 2. Rock Constitutive Model Based on the FLAC Software

### 2.1. Mohr–Coulomb Model of Rocks

The failure envelope of the Mohr–Coulomb model in FLAC is composed of the Mohr–Coulomb criterion (shearing yield function) and critical tensile stress (yield function of tensile stress). The flow rule of tensile stress is the flow rule of association, whereas the flow rule of shearing stress is the flow rule of nonassociation.


#### 2.1.1. Elasticity Rule Expressed by Increments

The FLAC3D calculates the principal stresses (σ1, σ2, and σ3) of the rock unit when using the Mohr–Coulomb model. The value and application directions of principal stresses are calculated from the tensor of stress and are ordered based on their values (pressure stress is negative).
(1)σ1≤σ2≤σ3.The corresponding principal strain increment can be expressed as(2)Δei=Δeie+Δeip,i=1−3,where superscripts e and p are the elasticity and plasticity parts of the strain. The plastic strain component is nonzero when the material is in the plastic flow state. The Hooker theorem expressed by the principal stress and principal strain increment can be written as follows:
(3)Δσ1=α1Δe1e+α2Δe2e+Δe3e,Δσ2=α1Δe2e+α2Δe1e+Δe3e,Δσ3=α1Δe3e+α2Δe2e+Δe1e,where α1=K+4G/3 and α2=K−2G/3. K and G are the bulk and shear moduli of rocks, respectively.


#### 2.1.2. Material Yield and Potential Function

On the basis of the ranking criteria of principal stress in (1), the failure criteria of materials can be determined by the envelope line on plane (σ1,σ3) (Figure 1).Figure 1
Mohr–Coulomb failure criteria.The failure envelope line of material AB section, which is defined by the Mohr–Coulomb yield function, is expressed as(4)fs=σ1−σ3Nϕ+2cNϕ.The BC section is defined by the tensile stress yield function.(5)ft=σt−σ3,where φ is the frictional angle, c is the cohesive force, and σt is the tensile strength of bulk materials.
(6)Nϕ=1+sinϕ1−sinϕ.The shearing stress yield criterion only considers the influences of the maximum and minimum principal stresses and assumes that the intermediate principal stress is insignificant. For materials with nonzero frictional angle, the tensile strength should be not higher than a specific valueσmaxt.
(7)σmaxt=ctanϕ.The potential function of shearing stressgs corresponds to the flow rule of nonassociation. This function is defined as follows:
(8)gs=σ1−σ3Nψ,where ψ is the dilation angle and
(9)Nψ=1+sinψ1−sinψ.Considering that the failure criterion of tensile stress is the flow rule of association, the following formula expresses that(10)gt=−σ3.For the unit in the 3D stress space under the influence of tensile shear stress, the flow rules of tensile and shearing stresses in the Mohr–Coulomb model on the side bearing the combined stress can be defined by one functionhσ1,σ3. hσ1,σ3 represents the bisector of shearing stress yield function fs=0 and the tensile stress yield function ft=0 on plane σ1,σ3 (Figure 2). Function h is defined as follows:
(11)h=σ3−σt+αPσ1−σP,where αP and σp are constants.
(12)αP=1+Nϕ2+Nϕ,σp=σtNϕ−2cNϕ.Figure 2
Plastic flow rule in the Mohr–Coulomb model.The failure state of one point in planeσ1,σ3 can be determined by function h based on the stress state at this point. As shown in Figure 2, zones 1 and 2 are the negative and positive zones of function h, respectively. If the stress state of one point is in zone 1, this point experiences shearing failure and the stress state of this point returns to shear yield line fs based on the flow rule regulated by potential function gs. If the stress state of one point is in zone 2, this point experiences tension and compression failures and the stress state of this point returns to tensile yield line ft=0 based on the flow rule regulated by potential function gt.On the basis of stress orders defined in (1), the points in the tensile-shear critical state automatically operate above the flow rules. The following method is easy to be realized in the program due to the small strain increments. In each calculation time step, only one plastic flow rule and corresponding stress correction are considered. Specifically, stress is alternatively corrected by two criteria when the stress state of one point lies on the boundary of the tensile-shear yield zone. In this process, the yield criteria can reach considerable calculation accuracy. This process depends on strain increment.


#### 2.1.3. Stress Correction under the Plastic State

Shear yield should be considered. The flow rule can be written as(13)Δeip=λ3∂g3∂σi,i=1,3,where λ3 is the unknown variable that should be determined. On the basis of the definition of function gs, the calculation of partial derivative expresses that
(14)Δe1p=λ3,Δe2p=0,Δe3p=−λ3NΨ.As shown in (2), the elastic strain increment is expressed as the total strain increment minus the plastic strain increment. On the basis of the flow rule described in (14), the law of elasticity in (3) is changed into
(15)Δσ1=σ1Δe1+α2Δe2+Δe3−λ3α1−α2Nψ,Δσ2=σ1Δe2+α2Δe1+Δe3−λ3α2−α2Nψ,Δσ3=σ1Δe3+α2Δe1+Δe2−λ3α2−α1Nψ.SuperscriptsN and O are the current and previous time step stress states. They are defined as
(16)σtN=σ10+Δσt,t=1.3.By combining (14) and (13), we obtain
(17)σtN=σ1∫=λ2a1−a2Nψ,σ2N=σ2∫=λ2a1−a2Nψ,σ3N=σ3∫=λ2a2−a1Nψ.SuperscriptI is the stress of the previous time step that is predicted based on the elasticity rule plus the stress increment. Therefore,
(18)σ1I=σ10+α1Δe1+α2Δe2+Δe3,σ2I=σ20+α1Δe2+α2Δe2+Δe3,σ3I=σ30+α1Δe3+α2Δe1+Δe2.The unknown parameterλf can be calculated by positioning the current stress state on the shearing yield plane. In the yield function f3=0, σ1 and σ3 are replaced by σ1N and σ3N. The transformation shows that
(19)λ3=f3σ1I⋅σ1Ia1−a2Nψ−a2−a1NψNϕ.In the following text, the tensile yield is considered. The flow rule of tension is expressed as follows:(20)Δe1p=λt∂gt∂σ1,t=1.3,where λ3 is an unknown variable. On the basis of the definition of gs, the calculation on partial derivative shows that
(21)Δe1p=0,Δe2p=0,Δe3p=−λt.Similar to the above deduction process of shear stress correction, we obtain(22)σ1N=σ1∫+λtα2,σ2N=σ2∫+λtα2,σ3N=σ3∫+λtα1,where
(23)λt=fψσ3tα1=αt−σ3tα1.


### 2.2. Computer Implementation of Plastic Flow Rule

To calculate the stress state of materials based on Mohr–Coulomb model, FLAC3D calculates the stress increment by the Hooker theorem based on the total strain increment of the current time step. Elastic budget stress (σy∫) is acquired by adding the calculated stress increment on the stress of the previous time step. The principle budget stresses (σ1∫, σ1t, and σ1I) and the corresponding principal directions are ordered based on the ranking rule regulated in (1). The elastic budget stress is corrected, and a new stress state is determined if the stress state satisfies the composite yield condition. Under the premise of meeting the yield criteria, two forms of stress state exist, namely, hσ1∫⋅σ3I≤0 and hσ1∫⋅σ3I≥0 ((11)). In the first situation, the occurrence of shear failure at this point is determined by the program. Stress is calculated again based on (17) in which λt is determined by (19). In the second situation, the materials develop tensile failure and the stress state is determined again based on (22) and (23). Suppose the main direction of the principal stress remains unchanged during the plastic stress correction. Therefore, the stress tensor component that corresponds to the system reference coordinates can be calculated based on the principal stress value.In FLAC3D, the tensile strength of materials is defaulted to zero. If the given tensile strength is higher than σmaxt, the tensile strength in the program is σmaxt. If the maximum principal stress that is calculated by one unit (σ3) is higher than tensile strength σt, this unit loses the tensile strength (0). In this way, the tensile softening effect of materials can be simulated.


## 2.1. Mohr–Coulomb Model of Rocks

The failure envelope of the Mohr–Coulomb model in FLAC is composed of the Mohr–Coulomb criterion (shearing yield function) and critical tensile stress (yield function of tensile stress). The flow rule of tensile stress is the flow rule of association, whereas the flow rule of shearing stress is the flow rule of nonassociation.


### 2.1.1. Elasticity Rule Expressed by Increments

The FLAC3D calculates the principal stresses (σ1, σ2, and σ3) of the rock unit when using the Mohr–Coulomb model. The value and application directions of principal stresses are calculated from the tensor of stress and are ordered based on their values (pressure stress is negative).
(1)σ1≤σ2≤σ3.The corresponding principal strain increment can be expressed as(2)Δei=Δeie+Δeip,i=1−3,where superscripts e and p are the elasticity and plasticity parts of the strain. The plastic strain component is nonzero when the material is in the plastic flow state. The Hooker theorem expressed by the principal stress and principal strain increment can be written as follows:
(3)Δσ1=α1Δe1e+α2Δe2e+Δe3e,Δσ2=α1Δe2e+α2Δe1e+Δe3e,Δσ3=α1Δe3e+α2Δe2e+Δe1e,where α1=K+4G/3 and α2=K−2G/3. K and G are the bulk and shear moduli of rocks, respectively.


### 2.1.2. Material Yield and Potential Function

On the basis of the ranking criteria of principal stress in (1), the failure criteria of materials can be determined by the envelope line on plane (σ1,σ3) (Figure 1).Figure 1
Mohr–Coulomb failure criteria.The failure envelope line of material AB section, which is defined by the Mohr–Coulomb yield function, is expressed as(4)fs=σ1−σ3Nϕ+2cNϕ.The BC section is defined by the tensile stress yield function.(5)ft=σt−σ3,where φ is the frictional angle, c is the cohesive force, and σt is the tensile strength of bulk materials.
(6)Nϕ=1+sinϕ1−sinϕ.The shearing stress yield criterion only considers the influences of the maximum and minimum principal stresses and assumes that the intermediate principal stress is insignificant. For materials with nonzero frictional angle, the tensile strength should be not higher than a specific valueσmaxt.
(7)σmaxt=ctanϕ.The potential function of shearing stressgs corresponds to the flow rule of nonassociation. This function is defined as follows:
(8)gs=σ1−σ3Nψ,where ψ is the dilation angle and
(9)Nψ=1+sinψ1−sinψ.Considering that the failure criterion of tensile stress is the flow rule of association, the following formula expresses that(10)gt=−σ3.For the unit in the 3D stress space under the influence of tensile shear stress, the flow rules of tensile and shearing stresses in the Mohr–Coulomb model on the side bearing the combined stress can be defined by one functionhσ1,σ3. hσ1,σ3 represents the bisector of shearing stress yield function fs=0 and the tensile stress yield function ft=0 on plane σ1,σ3 (Figure 2). Function h is defined as follows:
(11)h=σ3−σt+αPσ1−σP,where αP and σp are constants.
(12)αP=1+Nϕ2+Nϕ,σp=σtNϕ−2cNϕ.Figure 2
Plastic flow rule in the Mohr–Coulomb model.The failure state of one point in planeσ1,σ3 can be determined by function h based on the stress state at this point. As shown in Figure 2, zones 1 and 2 are the negative and positive zones of function h, respectively. If the stress state of one point is in zone 1, this point experiences shearing failure and the stress state of this point returns to shear yield line fs based on the flow rule regulated by potential function gs. If the stress state of one point is in zone 2, this point experiences tension and compression failures and the stress state of this point returns to tensile yield line ft=0 based on the flow rule regulated by potential function gt.On the basis of stress orders defined in (1), the points in the tensile-shear critical state automatically operate above the flow rules. The following method is easy to be realized in the program due to the small strain increments. In each calculation time step, only one plastic flow rule and corresponding stress correction are considered. Specifically, stress is alternatively corrected by two criteria when the stress state of one point lies on the boundary of the tensile-shear yield zone. In this process, the yield criteria can reach considerable calculation accuracy. This process depends on strain increment.


### 2.1.3. Stress Correction under the Plastic State

Shear yield should be considered. The flow rule can be written as(13)Δeip=λ3∂g3∂σi,i=1,3,where λ3 is the unknown variable that should be determined. On the basis of the definition of function gs, the calculation of partial derivative expresses that
(14)Δe1p=λ3,Δe2p=0,Δe3p=−λ3NΨ.As shown in (2), the elastic strain increment is expressed as the total strain increment minus the plastic strain increment. On the basis of the flow rule described in (14), the law of elasticity in (3) is changed into
(15)Δσ1=σ1Δe1+α2Δe2+Δe3−λ3α1−α2Nψ,Δσ2=σ1Δe2+α2Δe1+Δe3−λ3α2−α2Nψ,Δσ3=σ1Δe3+α2Δe1+Δe2−λ3α2−α1Nψ.SuperscriptsN and O are the current and previous time step stress states. They are defined as
(16)σtN=σ10+Δσt,t=1.3.By combining (14) and (13), we obtain
(17)σtN=σ1∫=λ2a1−a2Nψ,σ2N=σ2∫=λ2a1−a2Nψ,σ3N=σ3∫=λ2a2−a1Nψ.SuperscriptI is the stress of the previous time step that is predicted based on the elasticity rule plus the stress increment. Therefore,
(18)σ1I=σ10+α1Δe1+α2Δe2+Δe3,σ2I=σ20+α1Δe2+α2Δe2+Δe3,σ3I=σ30+α1Δe3+α2Δe1+Δe2.The unknown parameterλf can be calculated by positioning the current stress state on the shearing yield plane. In the yield function f3=0, σ1 and σ3 are replaced by σ1N and σ3N. The transformation shows that
(19)λ3=f3σ1I⋅σ1Ia1−a2Nψ−a2−a1NψNϕ.In the following text, the tensile yield is considered. The flow rule of tension is expressed as follows:(20)Δe1p=λt∂gt∂σ1,t=1.3,where λ3 is an unknown variable. On the basis of the definition of gs, the calculation on partial derivative shows that
(21)Δe1p=0,Δe2p=0,Δe3p=−λt.Similar to the above deduction process of shear stress correction, we obtain(22)σ1N=σ1∫+λtα2,σ2N=σ2∫+λtα2,σ3N=σ3∫+λtα1,where
(23)λt=fψσ3tα1=αt−σ3tα1.


## 2.1.1. Elasticity Rule Expressed by Increments

The FLAC3D calculates the principal stresses (σ1, σ2, and σ3) of the rock unit when using the Mohr–Coulomb model. The value and application directions of principal stresses are calculated from the tensor of stress and are ordered based on their values (pressure stress is negative).
(1)σ1≤σ2≤σ3.The corresponding principal strain increment can be expressed as(2)Δei=Δeie+Δeip,i=1−3,where superscripts e and p are the elasticity and plasticity parts of the strain. The plastic strain component is nonzero when the material is in the plastic flow state. The Hooker theorem expressed by the principal stress and principal strain increment can be written as follows:
(3)Δσ1=α1Δe1e+α2Δe2e+Δe3e,Δσ2=α1Δe2e+α2Δe1e+Δe3e,Δσ3=α1Δe3e+α2Δe2e+Δe1e,where α1=K+4G/3 and α2=K−2G/3. K and G are the bulk and shear moduli of rocks, respectively.


## 2.1.2. Material Yield and Potential Function

On the basis of the ranking criteria of principal stress in (1), the failure criteria of materials can be determined by the envelope line on plane (σ1,σ3) (Figure 1).Figure 1
Mohr–Coulomb failure criteria.The failure envelope line of material AB section, which is defined by the Mohr–Coulomb yield function, is expressed as(4)fs=σ1−σ3Nϕ+2cNϕ.The BC section is defined by the tensile stress yield function.(5)ft=σt−σ3,where φ is the frictional angle, c is the cohesive force, and σt is the tensile strength of bulk materials.
(6)Nϕ=1+sinϕ1−sinϕ.The shearing stress yield criterion only considers the influences of the maximum and minimum principal stresses and assumes that the intermediate principal stress is insignificant. For materials with nonzero frictional angle, the tensile strength should be not higher than a specific valueσmaxt.
(7)σmaxt=ctanϕ.The potential function of shearing stressgs corresponds to the flow rule of nonassociation. This function is defined as follows:
(8)gs=σ1−σ3Nψ,where ψ is the dilation angle and
(9)Nψ=1+sinψ1−sinψ.Considering that the failure criterion of tensile stress is the flow rule of association, the following formula expresses that(10)gt=−σ3.For the unit in the 3D stress space under the influence of tensile shear stress, the flow rules of tensile and shearing stresses in the Mohr–Coulomb model on the side bearing the combined stress can be defined by one functionhσ1,σ3. hσ1,σ3 represents the bisector of shearing stress yield function fs=0 and the tensile stress yield function ft=0 on plane σ1,σ3 (Figure 2). Function h is defined as follows:
(11)h=σ3−σt+αPσ1−σP,where αP and σp are constants.
(12)αP=1+Nϕ2+Nϕ,σp=σtNϕ−2cNϕ.Figure 2
Plastic flow rule in the Mohr–Coulomb model.The failure state of one point in planeσ1,σ3 can be determined by function h based on the stress state at this point. As shown in Figure 2, zones 1 and 2 are the negative and positive zones of function h, respectively. If the stress state of one point is in zone 1, this point experiences shearing failure and the stress state of this point returns to shear yield line fs based on the flow rule regulated by potential function gs. If the stress state of one point is in zone 2, this point experiences tension and compression failures and the stress state of this point returns to tensile yield line ft=0 based on the flow rule regulated by potential function gt.On the basis of stress orders defined in (1), the points in the tensile-shear critical state automatically operate above the flow rules. The following method is easy to be realized in the program due to the small strain increments. In each calculation time step, only one plastic flow rule and corresponding stress correction are considered. Specifically, stress is alternatively corrected by two criteria when the stress state of one point lies on the boundary of the tensile-shear yield zone. In this process, the yield criteria can reach considerable calculation accuracy. This process depends on strain increment.


## 2.1.3. Stress Correction under the Plastic State

Shear yield should be considered. The flow rule can be written as(13)Δeip=λ3∂g3∂σi,i=1,3,where λ3 is the unknown variable that should be determined. On the basis of the definition of function gs, the calculation of partial derivative expresses that
(14)Δe1p=λ3,Δe2p=0,Δe3p=−λ3NΨ.As shown in (2), the elastic strain increment is expressed as the total strain increment minus the plastic strain increment. On the basis of the flow rule described in (14), the law of elasticity in (3) is changed into
(15)Δσ1=σ1Δe1+α2Δe2+Δe3−λ3α1−α2Nψ,Δσ2=σ1Δe2+α2Δe1+Δe3−λ3α2−α2Nψ,Δσ3=σ1Δe3+α2Δe1+Δe2−λ3α2−α1Nψ.SuperscriptsN and O are the current and previous time step stress states. They are defined as
(16)σtN=σ10+Δσt,t=1.3.By combining (14) and (13), we obtain
(17)σtN=σ1∫=λ2a1−a2Nψ,σ2N=σ2∫=λ2a1−a2Nψ,σ3N=σ3∫=λ2a2−a1Nψ.SuperscriptI is the stress of the previous time step that is predicted based on the elasticity rule plus the stress increment. Therefore,
(18)σ1I=σ10+α1Δe1+α2Δe2+Δe3,σ2I=σ20+α1Δe2+α2Δe2+Δe3,σ3I=σ30+α1Δe3+α2Δe1+Δe2.The unknown parameterλf can be calculated by positioning the current stress state on the shearing yield plane. In the yield function f3=0, σ1 and σ3 are replaced by σ1N and σ3N. The transformation shows that
(19)λ3=f3σ1I⋅σ1Ia1−a2Nψ−a2−a1NψNϕ.In the following text, the tensile yield is considered. The flow rule of tension is expressed as follows:(20)Δe1p=λt∂gt∂σ1,t=1.3,where λ3 is an unknown variable. On the basis of the definition of gs, the calculation on partial derivative shows that
(21)Δe1p=0,Δe2p=0,Δe3p=−λt.Similar to the above deduction process of shear stress correction, we obtain(22)σ1N=σ1∫+λtα2,σ2N=σ2∫+λtα2,σ3N=σ3∫+λtα1,where
(23)λt=fψσ3tα1=αt−σ3tα1.


## 2.2. Computer Implementation of Plastic Flow Rule

To calculate the stress state of materials based on Mohr–Coulomb model, FLAC3D calculates the stress increment by the Hooker theorem based on the total strain increment of the current time step. Elastic budget stress (σy∫) is acquired by adding the calculated stress increment on the stress of the previous time step. The principle budget stresses (σ1∫, σ1t, and σ1I) and the corresponding principal directions are ordered based on the ranking rule regulated in (1). The elastic budget stress is corrected, and a new stress state is determined if the stress state satisfies the composite yield condition. Under the premise of meeting the yield criteria, two forms of stress state exist, namely, hσ1∫⋅σ3I≤0 and hσ1∫⋅σ3I≥0 ((11)). In the first situation, the occurrence of shear failure at this point is determined by the program. Stress is calculated again based on (17) in which λt is determined by (19). In the second situation, the materials develop tensile failure and the stress state is determined again based on (22) and (23). Suppose the main direction of the principal stress remains unchanged during the plastic stress correction. Therefore, the stress tensor component that corresponds to the system reference coordinates can be calculated based on the principal stress value.In FLAC3D, the tensile strength of materials is defaulted to zero. If the given tensile strength is higher than σmaxt, the tensile strength in the program is σmaxt. If the maximum principal stress that is calculated by one unit (σ3) is higher than tensile strength σt, this unit loses the tensile strength (0). In this way, the tensile softening effect of materials can be simulated.


## 3. Calculation Model and Calculation Concept

### 3.1. Calculation Model

To analyze the effect of blasting in tunnels on the embankment of the Yellow River, the following parts are selected for the 3D finite element calculation modeling (Figure3). This modeling zone includes different geological strata and three faults (f11, f12, and f14) that cover a total length of 120 m.Figure 3
Position sketch of the 3D calculation model.On the basis of the above regions with considerations to the grouting curtain area, a method similar to the plane model is adopted. The general large 3D finite element model is shown in Figure4. The grouting curtain area is shown in Figure 5, and the water-conveyance tunnel is shown in Figure 6. The three faults are shown in Figure 7.Figure 4
3D computation mode.Figure 5
Grouting curtain.Figure 6
Water-conveyance tunnel.Figure 7
Faults.


### 3.2. Calculation Loads

Similar to plane computation, three major types of loads are used:(1)
Gravity loads
(2)
Water loads
(3)
Blasting seismic wavesThe gravity load is applied as the physical power. Water load is applied at top of the model as the uniformly distributed loads. The blasting seismic wave is applied on the chamber surrounding. On the basis of the provided calculation data, the velocity wave is selected as the applied seismic wave (Figure8).Figure 8
Blasting velocity wave.


### 3.3. Calculation Parameters

On the basis of the provided data, the following parameters are selected as material parameters in this calculation (Table.1).Table 1
Material parameters.
Strata
Density
Elastic modulus
Poisson ratio
Friction angle
Cohesion
Tensile strength
Kg/m3
Pa
°
Pa
Pa
Modern alluvial siltQ43al
2680
0.8e9
0.42
30
0.1e6
0.0e6
C3g2
2700
5.0e9
0.40
35
0.6e6
2.0e6
C3g1
2720
12.0e9
0.35
40
2.0e6
3.0e6
C2z
2750
18.0e9
0.30
45
2.50e6
7.0e6
Fault
2680
0.8e9
0.42
30
0.1e6
0.0e6
Grout curtain
3025
19.8e9
0.33
49.5
2.75e6
7.70e6


### 3.4. Calculation Concept

On the basis of the “advanced pregrouting construction scheme,” that is, 50 m of pregrouting and 40 m of excavation requirements, the water-conveyance tunnel in this model is divided into four sections (Figure9). In the calculation process, the blasting seismic wave is applied on the face and tunnel wall at each excavation. The selected region on the tunnel wall is 10 m away from the face (Figure 10). The blasting seismic wave is applied on the tunnel wall in the final excavation, that is, the fourth excavation. The application range is also 10 m away from the exit (Figure 11).Figure 9
Segmented excavation and zoning of blasting seismic wave application.Figure 10
First excavation and blasting seismic wave application.Figure 11
Fourth excavation and blasting seismic wave application.


## 3.1. Calculation Model

To analyze the effect of blasting in tunnels on the embankment of the Yellow River, the following parts are selected for the 3D finite element calculation modeling (Figure3). This modeling zone includes different geological strata and three faults (f11, f12, and f14) that cover a total length of 120 m.Figure 3
Position sketch of the 3D calculation model.On the basis of the above regions with considerations to the grouting curtain area, a method similar to the plane model is adopted. The general large 3D finite element model is shown in Figure4. The grouting curtain area is shown in Figure 5, and the water-conveyance tunnel is shown in Figure 6. The three faults are shown in Figure 7.Figure 4
3D computation mode.Figure 5
Grouting curtain.Figure 6
Water-conveyance tunnel.Figure 7
Faults.


## 3.2. Calculation Loads

Similar to plane computation, three major types of loads are used:(1)
Gravity loads
(2)
Water loads
(3)
Blasting seismic wavesThe gravity load is applied as the physical power. Water load is applied at top of the model as the uniformly distributed loads. The blasting seismic wave is applied on the chamber surrounding. On the basis of the provided calculation data, the velocity wave is selected as the applied seismic wave (Figure8).Figure 8
Blasting velocity wave.


## 3.3. Calculation Parameters

On the basis of the provided data, the following parameters are selected as material parameters in this calculation (Table.1).Table 1
Material parameters.
Strata
Density
Elastic modulus
Poisson ratio
Friction angle
Cohesion
Tensile strength
Kg/m3
Pa
°
Pa
Pa
Modern alluvial siltQ43al
2680
0.8e9
0.42
30
0.1e6
0.0e6
C3g2
2700
5.0e9
0.40
35
0.6e6
2.0e6
C3g1
2720
12.0e9
0.35
40
2.0e6
3.0e6
C2z
2750
18.0e9
0.30
45
2.50e6
7.0e6
Fault
2680
0.8e9
0.42
30
0.1e6
0.0e6
Grout curtain
3025
19.8e9
0.33
49.5
2.75e6
7.70e6


## 3.4. Calculation Concept

On the basis of the “advanced pregrouting construction scheme,” that is, 50 m of pregrouting and 40 m of excavation requirements, the water-conveyance tunnel in this model is divided into four sections (Figure9). In the calculation process, the blasting seismic wave is applied on the face and tunnel wall at each excavation. The selected region on the tunnel wall is 10 m away from the face (Figure 10). The blasting seismic wave is applied on the tunnel wall in the final excavation, that is, the fourth excavation. The application range is also 10 m away from the exit (Figure 11).Figure 9
Segmented excavation and zoning of blasting seismic wave application.Figure 10
First excavation and blasting seismic wave application.Figure 11
Fourth excavation and blasting seismic wave application.


## 4. Calculation Result Analysis

### 4.1. Selection of the Analyzed Fracture and Characteristic Points

To analyze the influences of excavation steps on rock mass, grouting curtain, and embankment, four sections (A-A, B-B, C-C, and D-D) are selected from the model zone for analysis (Figure12). To analyze the influences of excavation steps on embankment, 6 characteristic points are selected from the embankment surface for calculation analysis. The 6 characteristic points are located in the center in which one side of the computation zone has 3 points each. The locations of different characteristic points are shown in Figure 13.Figure 12
Excavation fractures.Figure 13
Key points on the embankment surface.


### 4.2. Analysis on the Profile Calculation Results

For an accurate analysis on the maximum effect of blasting seismic wave on rock mass, grouting curtain, and embankment, only the section where the explosion source lies is analyzed. In other words, only section A-A is analyzed in the first explosion, and the stress displacements on the remaining three sections are not analyzed. The rest can be performed in the same manner.


#### 4.2.1. Stress and Displacement Analysis on Section A-A during the First Excavation

As shown in Figures14–16 under the effect of blasting seismic wave, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass. The maximum horizontal and maximum vertical displacements of the tunnel wall are 60.00 and 60.00 mm, respectively. The maximum resultant displacement is 65.00 mm. In the grouting curtain zone, the horizontal displacement mainly ranges between 0.00 and 30.00 mm, the vertical displacement concentrates in 0.00–30.00 mm, and the resultant displacement is mainly in the range of 20.00~40.00 mm.Figure 14
Displacement of section A-A along X.Figure 15
Displacement of section A-A along Z.Figure 16
Resultant displacement of section A-A.The distribution law of the major principal stress is shown in Figure17. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to the compression state. The grouting curtain zone mainly concentrates in the compressive zones. The maximum tensile stress on the tunnel wall is 1.20 MPa, and the compressive stress in the grouting curtain is −0.20~−0.40 MPa.Figure 17
Major principal stress on section A-A.


#### 4.2.2. Stress and Displacement Analysis on Section B-B during the Second Excavation

As shown in Figures18–20 under the effect of blasting seismic wave, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass. The maximum horizontal and maximum vertical displacements of the tunnel wall are 43.00 and 43.20 mm, respectively. The maximum resultant displacement is 44.00 mm. In the grouting curtain zone, the horizontal displacement mainly ranges between 2.10 and 20.00 mm, the vertical displacement concentrates in 2.27~15.60 mm, and the resultant displacement is mainly in the range of 17.00~27.00 mm.Figure 18
Displacement of section B-B along X.Figure 19
Displacement of section B-B along Z.Figure 20
Resultant displacement of section B-B.The distribution law of the major principal stress is shown in Figure21. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to compression state. The grouting curtain zone mainly concentrates in the compressive zones. The maximum tensile stress on the tunnel wall is 1.65 MPa, and the compressive stress in the grouting curtain is −0.50~−1.00 MPa.Figure 21
Major principal stress on section B-B.


#### 4.2.3. Stress and Displacement Analysis on Section C-C during the Third Excavation

As shown in Figures22–24 under the effect of blasting seismic wave, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass. The maximum horizontal and maximum vertical displacements of the tunnel wall are 50.00 and 58.00 mm, respectively. The maximum resultant displacement is 58.00 mm. In the grouting curtain zone, the horizontal displacement mainly ranges between 2.63 and 28.95 mm, the vertical displacement concentrates between 0.00 and 33.00 mm, and the resultant displacement is mainly in the range of 25.80~40.44 mm.Figure 22
Displacement of section C-C along X.Figure 23
Displacement of section C-C along Z.Figure 24
Resultant displacement of section C-C.The distribution law of the major principal stress is shown in Figure25. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to compression state. The grouting curtain zone mainly concentrates in the compressive zones. The maximum tensile stress on the tunnel wall is 1.35 MPa, and the compressive stress in the grouting curtain is −0.20~−0.72 MPa.Figure 25
Major principal stress on section C-C.


#### 4.2.4. Stress and Displacement Analysis on Section D-D during the Fourth Excavation

As shown in Figures26–28 under the effect of blasting seismic wave, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass. The maximum horizontal and maximum vertical displacements of the tunnel wall are 35.00 and 35.00 mm, respectively. The maximum resultant displacement is 38.00 mm. In the grouting curtain zone, the horizontal displacement mainly ranges between 0.00 and 20.00 mm, the vertical displacement concentrates between 5.00 and 20.00 mm, and the resultant displacement is mainly in the range of 16.00~24.00 mm.Figure 26
Displacement of section D-D along X.Figure 27
Displacement of section D-D along Z.Figure 28
Resultant displacement of section D-D.The distribution law of the major principal stress is shown in Figure29. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to compression state. The grouting curtain zone mainly concentrates in the compressive zones. The maximum tensile stress on the tunnel wall is 1.40 MPa, and the compressive stress in the grouting curtain is −0.04~−0.40 MPa.Figure 29
Major principal stress on section D-D.On the basis of the above displacement and stress maps, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass under the effect of blasting seismic wave. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to compression state. The curtain-grouting zone mainly concentrates in the compressive zones.


### 4.3. Time-History Response Analysis on Characteristic Points

#### 4.3.1. Time-History Response Analysis on Velocity at the Characteristic Points on the Embankment under Different Excavation Steps

As shown in Figure30, the maximum vibration velocities at characteristic points P1, P2, P3, P4, P5, and P6 are 2.47, 2.13, 1.36, 1.88, 2.24, and 0.95 cm/s, respectively.Figure 30
Time-history response curve of velocity at characteristic points on the embankment in the first excavation.As shown in Figure31, the maximum vibration velocities at P1, P2, P3, P4, P5, and P6 are 3.28, 3.94, 3.43, 3.00, 4.33, and 3.54 cm/s, respectively.Figure 31
Time-history response curve of velocity at characteristic points on the embankment in the second excavation.Figure32 shows that the maximum vibration velocities at P1, P2, P3, P4, P5, and P6 are 1.17, 2.38, 3.00, 0.79, 1.22, and 1.59 cm/s, respectively.Figure 32
Time-history response curve of velocity at characteristic points on the embankment in the third excavation.As shown in Figure33, the maximum vibration velocities at P1, P2, P3, P4, P5, and P6 are 2.35, 2.46, 2.65, 1.81, 2.10, and 2.19 cm/s, respectively.Figure 33
Time-history response curve of velocity at characteristic points on the embankment in the fourth excavation.As shown in Figures30–33 in each excavation step, the vibration velocity is high in the beginning of blasting and gradually decreases until it becomes stable. The vibration velocities at all characteristic points at each excavation step are generally small that range between 0.79 and 4.33 cm/s.


#### 4.3.2. Time-History Response Analysis of Resultant Displacement at Characteristic Points under Different Excavation Steps

As shown in Figure34, the maximum vibration displacements at P1, P2, P3, P4, P5, and P6 are 2.69, 3.09, 5.33, 2.27, 2.86, and 5.22 mm, respectively.Figure 34
Time-history response curve of resultant displacement at characteristic points on the embankment in the first excavation.As shown in Figure35, the maximum vibration displacements at P1, P2, P3, P4, P5, and P6 are 13.31, 17.29, 19.20, 12.51, 16.29, and 18.56 mm, respectively.Figure 35
Time-history response curve of resultant displacement at characteristic points on the embankment in the second excavation.As shown in Figure36, the maximum vibration displacements at P1, P2, P3, P4, P5, and P6 are 16.60, 13.86, 15.28, 16.23, 13.71, and 14.57 mm, respectively.Figure 36
Time-history response curve of resultant displacement at characteristic points on the embankment in the third excavation.As shown in Figure37, the maximum vibration displacements at P1, P2, P3, P4, P5, and P6 are 16.17, 18.24, 19.93, 14.73, 17.29, and 19.19 mm, respectively.Figure 37
Time-history response curve of resultant displacement at characteristic points on the embankment in the fourth excavation.As shown in Figures34–37 in each excavation step, the resultant displacement at different characteristic points caused by blasting excavation continuously increases. This condition demonstrates that tunnel blasting causes a continuously increasing resultant deformation of the earth surface. The resultant displacement caused by the first step of blasting excavation is relatively small and gradually increases in the following excavation steps. This condition is because of the tunnel that begins to be connected slowly that results in considerable cavities. Hence, the resultant displacement gradually increases. Displacements of the representative points P1, P2, and P3 are much greater than those of the representative points P4, P5, and P6 due to different positions of these points, to specific, the representative points P1, P2, and P3 are above the excavation tunnel and the representative points P4, P5, and P6 are on one side of the excavation tunnel.


## 4.1. Selection of the Analyzed Fracture and Characteristic Points

To analyze the influences of excavation steps on rock mass, grouting curtain, and embankment, four sections (A-A, B-B, C-C, and D-D) are selected from the model zone for analysis (Figure12). To analyze the influences of excavation steps on embankment, 6 characteristic points are selected from the embankment surface for calculation analysis. The 6 characteristic points are located in the center in which one side of the computation zone has 3 points each. The locations of different characteristic points are shown in Figure 13.Figure 12
Excavation fractures.Figure 13
Key points on the embankment surface.


## 4.2. Analysis on the Profile Calculation Results

For an accurate analysis on the maximum effect of blasting seismic wave on rock mass, grouting curtain, and embankment, only the section where the explosion source lies is analyzed. In other words, only section A-A is analyzed in the first explosion, and the stress displacements on the remaining three sections are not analyzed. The rest can be performed in the same manner.


### 4.2.1. Stress and Displacement Analysis on Section A-A during the First Excavation

As shown in Figures14–16 under the effect of blasting seismic wave, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass. The maximum horizontal and maximum vertical displacements of the tunnel wall are 60.00 and 60.00 mm, respectively. The maximum resultant displacement is 65.00 mm. In the grouting curtain zone, the horizontal displacement mainly ranges between 0.00 and 30.00 mm, the vertical displacement concentrates in 0.00–30.00 mm, and the resultant displacement is mainly in the range of 20.00~40.00 mm.Figure 14
Displacement of section A-A along X.Figure 15
Displacement of section A-A along Z.Figure 16
Resultant displacement of section A-A.The distribution law of the major principal stress is shown in Figure17. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to the compression state. The grouting curtain zone mainly concentrates in the compressive zones. The maximum tensile stress on the tunnel wall is 1.20 MPa, and the compressive stress in the grouting curtain is −0.20~−0.40 MPa.Figure 17
Major principal stress on section A-A.


### 4.2.2. Stress and Displacement Analysis on Section B-B during the Second Excavation

As shown in Figures18–20 under the effect of blasting seismic wave, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass. The maximum horizontal and maximum vertical displacements of the tunnel wall are 43.00 and 43.20 mm, respectively. The maximum resultant displacement is 44.00 mm. In the grouting curtain zone, the horizontal displacement mainly ranges between 2.10 and 20.00 mm, the vertical displacement concentrates in 2.27~15.60 mm, and the resultant displacement is mainly in the range of 17.00~27.00 mm.Figure 18
Displacement of section B-B along X.Figure 19
Displacement of section B-B along Z.Figure 20
Resultant displacement of section B-B.The distribution law of the major principal stress is shown in Figure21. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to compression state. The grouting curtain zone mainly concentrates in the compressive zones. The maximum tensile stress on the tunnel wall is 1.65 MPa, and the compressive stress in the grouting curtain is −0.50~−1.00 MPa.Figure 21
Major principal stress on section B-B.


### 4.2.3. Stress and Displacement Analysis on Section C-C during the Third Excavation

As shown in Figures22–24 under the effect of blasting seismic wave, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass. The maximum horizontal and maximum vertical displacements of the tunnel wall are 50.00 and 58.00 mm, respectively. The maximum resultant displacement is 58.00 mm. In the grouting curtain zone, the horizontal displacement mainly ranges between 2.63 and 28.95 mm, the vertical displacement concentrates between 0.00 and 33.00 mm, and the resultant displacement is mainly in the range of 25.80~40.44 mm.Figure 22
Displacement of section C-C along X.Figure 23
Displacement of section C-C along Z.Figure 24
Resultant displacement of section C-C.The distribution law of the major principal stress is shown in Figure25. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to compression state. The grouting curtain zone mainly concentrates in the compressive zones. The maximum tensile stress on the tunnel wall is 1.35 MPa, and the compressive stress in the grouting curtain is −0.20~−0.72 MPa.Figure 25
Major principal stress on section C-C.


### 4.2.4. Stress and Displacement Analysis on Section D-D during the Fourth Excavation

As shown in Figures26–28 under the effect of blasting seismic wave, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass. The maximum horizontal and maximum vertical displacements of the tunnel wall are 35.00 and 35.00 mm, respectively. The maximum resultant displacement is 38.00 mm. In the grouting curtain zone, the horizontal displacement mainly ranges between 0.00 and 20.00 mm, the vertical displacement concentrates between 5.00 and 20.00 mm, and the resultant displacement is mainly in the range of 16.00~24.00 mm.Figure 26
Displacement of section D-D along X.Figure 27
Displacement of section D-D along Z.Figure 28
Resultant displacement of section D-D.The distribution law of the major principal stress is shown in Figure29. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to compression state. The grouting curtain zone mainly concentrates in the compressive zones. The maximum tensile stress on the tunnel wall is 1.40 MPa, and the compressive stress in the grouting curtain is −0.04~−0.40 MPa.Figure 29
Major principal stress on section D-D.On the basis of the above displacement and stress maps, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass under the effect of blasting seismic wave. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to compression state. The curtain-grouting zone mainly concentrates in the compressive zones.


## 4.2.1. Stress and Displacement Analysis on Section A-A during the First Excavation

As shown in Figures14–16 under the effect of blasting seismic wave, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass. The maximum horizontal and maximum vertical displacements of the tunnel wall are 60.00 and 60.00 mm, respectively. The maximum resultant displacement is 65.00 mm. In the grouting curtain zone, the horizontal displacement mainly ranges between 0.00 and 30.00 mm, the vertical displacement concentrates in 0.00–30.00 mm, and the resultant displacement is mainly in the range of 20.00~40.00 mm.Figure 14
Displacement of section A-A along X.Figure 15
Displacement of section A-A along Z.Figure 16
Resultant displacement of section A-A.The distribution law of the major principal stress is shown in Figure17. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to the compression state. The grouting curtain zone mainly concentrates in the compressive zones. The maximum tensile stress on the tunnel wall is 1.20 MPa, and the compressive stress in the grouting curtain is −0.20~−0.40 MPa.Figure 17
Major principal stress on section A-A.


## 4.2.2. Stress and Displacement Analysis on Section B-B during the Second Excavation

As shown in Figures18–20 under the effect of blasting seismic wave, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass. The maximum horizontal and maximum vertical displacements of the tunnel wall are 43.00 and 43.20 mm, respectively. The maximum resultant displacement is 44.00 mm. In the grouting curtain zone, the horizontal displacement mainly ranges between 2.10 and 20.00 mm, the vertical displacement concentrates in 2.27~15.60 mm, and the resultant displacement is mainly in the range of 17.00~27.00 mm.Figure 18
Displacement of section B-B along X.Figure 19
Displacement of section B-B along Z.Figure 20
Resultant displacement of section B-B.The distribution law of the major principal stress is shown in Figure21. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to compression state. The grouting curtain zone mainly concentrates in the compressive zones. The maximum tensile stress on the tunnel wall is 1.65 MPa, and the compressive stress in the grouting curtain is −0.50~−1.00 MPa.Figure 21
Major principal stress on section B-B.


## 4.2.3. Stress and Displacement Analysis on Section C-C during the Third Excavation

As shown in Figures22–24 under the effect of blasting seismic wave, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass. The maximum horizontal and maximum vertical displacements of the tunnel wall are 50.00 and 58.00 mm, respectively. The maximum resultant displacement is 58.00 mm. In the grouting curtain zone, the horizontal displacement mainly ranges between 2.63 and 28.95 mm, the vertical displacement concentrates between 0.00 and 33.00 mm, and the resultant displacement is mainly in the range of 25.80~40.44 mm.Figure 22
Displacement of section C-C along X.Figure 23
Displacement of section C-C along Z.Figure 24
Resultant displacement of section C-C.The distribution law of the major principal stress is shown in Figure25. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to compression state. The grouting curtain zone mainly concentrates in the compressive zones. The maximum tensile stress on the tunnel wall is 1.35 MPa, and the compressive stress in the grouting curtain is −0.20~−0.72 MPa.Figure 25
Major principal stress on section C-C.


## 4.2.4. Stress and Displacement Analysis on Section D-D during the Fourth Excavation

As shown in Figures26–28 under the effect of blasting seismic wave, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass. The maximum horizontal and maximum vertical displacements of the tunnel wall are 35.00 and 35.00 mm, respectively. The maximum resultant displacement is 38.00 mm. In the grouting curtain zone, the horizontal displacement mainly ranges between 0.00 and 20.00 mm, the vertical displacement concentrates between 5.00 and 20.00 mm, and the resultant displacement is mainly in the range of 16.00~24.00 mm.Figure 26
Displacement of section D-D along X.Figure 27
Displacement of section D-D along Z.Figure 28
Resultant displacement of section D-D.The distribution law of the major principal stress is shown in Figure29. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to compression state. The grouting curtain zone mainly concentrates in the compressive zones. The maximum tensile stress on the tunnel wall is 1.40 MPa, and the compressive stress in the grouting curtain is −0.04~−0.40 MPa.Figure 29
Major principal stress on section D-D.On the basis of the above displacement and stress maps, the displacement generally attenuates from the tunnel wall to the internal structure of rock mass under the effect of blasting seismic wave. The tunnel wall is mainly in the tension state. The tensile stress gradually weakens with its deep concentration on the rock mass. The rock mass turns to compression state. The curtain-grouting zone mainly concentrates in the compressive zones.


## 4.3. Time-History Response Analysis on Characteristic Points

### 4.3.1. Time-History Response Analysis on Velocity at the Characteristic Points on the Embankment under Different Excavation Steps

As shown in Figure30, the maximum vibration velocities at characteristic points P1, P2, P3, P4, P5, and P6 are 2.47, 2.13, 1.36, 1.88, 2.24, and 0.95 cm/s, respectively.Figure 30
Time-history response curve of velocity at characteristic points on the embankment in the first excavation.As shown in Figure31, the maximum vibration velocities at P1, P2, P3, P4, P5, and P6 are 3.28, 3.94, 3.43, 3.00, 4.33, and 3.54 cm/s, respectively.Figure 31
Time-history response curve of velocity at characteristic points on the embankment in the second excavation.Figure32 shows that the maximum vibration velocities at P1, P2, P3, P4, P5, and P6 are 1.17, 2.38, 3.00, 0.79, 1.22, and 1.59 cm/s, respectively.Figure 32
Time-history response curve of velocity at characteristic points on the embankment in the third excavation.As shown in Figure33, the maximum vibration velocities at P1, P2, P3, P4, P5, and P6 are 2.35, 2.46, 2.65, 1.81, 2.10, and 2.19 cm/s, respectively.Figure 33
Time-history response curve of velocity at characteristic points on the embankment in the fourth excavation.As shown in Figures30–33 in each excavation step, the vibration velocity is high in the beginning of blasting and gradually decreases until it becomes stable. The vibration velocities at all characteristic points at each excavation step are generally small that range between 0.79 and 4.33 cm/s.


### 4.3.2. Time-History Response Analysis of Resultant Displacement at Characteristic Points under Different Excavation Steps

As shown in Figure34, the maximum vibration displacements at P1, P2, P3, P4, P5, and P6 are 2.69, 3.09, 5.33, 2.27, 2.86, and 5.22 mm, respectively.Figure 34
Time-history response curve of resultant displacement at characteristic points on the embankment in the first excavation.As shown in Figure35, the maximum vibration displacements at P1, P2, P3, P4, P5, and P6 are 13.31, 17.29, 19.20, 12.51, 16.29, and 18.56 mm, respectively.Figure 35
Time-history response curve of resultant displacement at characteristic points on the embankment in the second excavation.As shown in Figure36, the maximum vibration displacements at P1, P2, P3, P4, P5, and P6 are 16.60, 13.86, 15.28, 16.23, 13.71, and 14.57 mm, respectively.Figure 36
Time-history response curve of resultant displacement at characteristic points on the embankment in the third excavation.As shown in Figure37, the maximum vibration displacements at P1, P2, P3, P4, P5, and P6 are 16.17, 18.24, 19.93, 14.73, 17.29, and 19.19 mm, respectively.Figure 37
Time-history response curve of resultant displacement at characteristic points on the embankment in the fourth excavation.As shown in Figures34–37 in each excavation step, the resultant displacement at different characteristic points caused by blasting excavation continuously increases. This condition demonstrates that tunnel blasting causes a continuously increasing resultant deformation of the earth surface. The resultant displacement caused by the first step of blasting excavation is relatively small and gradually increases in the following excavation steps. This condition is because of the tunnel that begins to be connected slowly that results in considerable cavities. Hence, the resultant displacement gradually increases. Displacements of the representative points P1, P2, and P3 are much greater than those of the representative points P4, P5, and P6 due to different positions of these points, to specific, the representative points P1, P2, and P3 are above the excavation tunnel and the representative points P4, P5, and P6 are on one side of the excavation tunnel.


## 4.3.1. Time-History Response Analysis on Velocity at the Characteristic Points on the Embankment under Different Excavation Steps

As shown in Figure30, the maximum vibration velocities at characteristic points P1, P2, P3, P4, P5, and P6 are 2.47, 2.13, 1.36, 1.88, 2.24, and 0.95 cm/s, respectively.Figure 30
Time-history response curve of velocity at characteristic points on the embankment in the first excavation.As shown in Figure31, the maximum vibration velocities at P1, P2, P3, P4, P5, and P6 are 3.28, 3.94, 3.43, 3.00, 4.33, and 3.54 cm/s, respectively.Figure 31
Time-history response curve of velocity at characteristic points on the embankment in the second excavation.Figure32 shows that the maximum vibration velocities at P1, P2, P3, P4, P5, and P6 are 1.17, 2.38, 3.00, 0.79, 1.22, and 1.59 cm/s, respectively.Figure 32
Time-history response curve of velocity at characteristic points on the embankment in the third excavation.As shown in Figure33, the maximum vibration velocities at P1, P2, P3, P4, P5, and P6 are 2.35, 2.46, 2.65, 1.81, 2.10, and 2.19 cm/s, respectively.Figure 33
Time-history response curve of velocity at characteristic points on the embankment in the fourth excavation.As shown in Figures30–33 in each excavation step, the vibration velocity is high in the beginning of blasting and gradually decreases until it becomes stable. The vibration velocities at all characteristic points at each excavation step are generally small that range between 0.79 and 4.33 cm/s.


## 4.3.2. Time-History Response Analysis of Resultant Displacement at Characteristic Points under Different Excavation Steps

As shown in Figure34, the maximum vibration displacements at P1, P2, P3, P4, P5, and P6 are 2.69, 3.09, 5.33, 2.27, 2.86, and 5.22 mm, respectively.Figure 34
Time-history response curve of resultant displacement at characteristic points on the embankment in the first excavation.As shown in Figure35, the maximum vibration displacements at P1, P2, P3, P4, P5, and P6 are 13.31, 17.29, 19.20, 12.51, 16.29, and 18.56 mm, respectively.Figure 35
Time-history response curve of resultant displacement at characteristic points on the embankment in the second excavation.As shown in Figure36, the maximum vibration displacements at P1, P2, P3, P4, P5, and P6 are 16.60, 13.86, 15.28, 16.23, 13.71, and 14.57 mm, respectively.Figure 36
Time-history response curve of resultant displacement at characteristic points on the embankment in the third excavation.As shown in Figure37, the maximum vibration displacements at P1, P2, P3, P4, P5, and P6 are 16.17, 18.24, 19.93, 14.73, 17.29, and 19.19 mm, respectively.Figure 37
Time-history response curve of resultant displacement at characteristic points on the embankment in the fourth excavation.As shown in Figures34–37 in each excavation step, the resultant displacement at different characteristic points caused by blasting excavation continuously increases. This condition demonstrates that tunnel blasting causes a continuously increasing resultant deformation of the earth surface. The resultant displacement caused by the first step of blasting excavation is relatively small and gradually increases in the following excavation steps. This condition is because of the tunnel that begins to be connected slowly that results in considerable cavities. Hence, the resultant displacement gradually increases. Displacements of the representative points P1, P2, and P3 are much greater than those of the representative points P4, P5, and P6 due to different positions of these points, to specific, the representative points P1, P2, and P3 are above the excavation tunnel and the representative points P4, P5, and P6 are on one side of the excavation tunnel.


## 5. Conclusions

On the basis of the 3D calculation on tunnel blasting, several major conclusions can be drawn which are as follows:(1)
The displacement generally attenuates from the tunnel wall to the internal structure of rock mass under the effect of blasting seismic wave. The tunnel wall is mainly in the tension state. The tensile stress is gradually converted into compressive stress with its deep concentration on the rock mass. The curtain-grouting zone mainly concentrates in the compressive zones.
(2)
On the basis of the calculation of characteristic points at the top of the embankment, the maximum vibration velocity at different characteristic points ranges between 0.79 and 4.33 cm/s under different blasting excavation steps. The maximum displacement along X changes between −13.11 and 13.28 mm, the maximum displacement along Y fluctuates within −12.39~14.97 mm, and the maximum displacement along Z ranges between −19.12 and 5.23 mm. The variation range of the maximum resultant displacement is in the range of 2.27~19.93 mm. All displacements are in the allowed range.
(3)
Considerable tunnel holes occur and the resultant displacement gradually increases with the continuation of blasting excavation.
(4)
The blasting vibration wave velocity is lower than the safety value in the code based on the vibration wave velocity at different characteristic points on the embankment. Therefore, the embankment is stable and safe during tunnel blasting.
(5)
A finite element numerical simulation on complicated problems in large actual engineering projects is conducted. The subsequent construction is guided based on the calculation results that achieves ideal outcomes. The proposed numerical simulation is significantly fast, convenient, and accurate and can be promoted in actual engineering.



---

*Source: 4654315-2018-08-12.xml*