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

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 FLAC 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.


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][4][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][18][19][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][22][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.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).

Rock Constitutive Model
The corresponding principal strain increment can be expressed as 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: 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).The failure envelope line of material AB section, which is defined by the Mohr-Coulomb yield function, is expressed as The BC section is defined by the tensile stress yield function.
where φ is the frictional angle, c is the cohesive force, and σ t is the tensile strength of bulk materials.
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 σ t max .
The potential function of shearing stress g s corresponds to the flow rule of nonassociation.This function is defined as follows: where ψ is the dilation angle and Considering that the failure criterion of tensile stress is the flow rule of association, the following formula expresses that 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 function h σ 1 , σ 3 .h σ 1 , σ 3 represents the bisector of shearing stress yield function f s = 0 and the tensile stress yield function f t = 0 on plane σ 1 , σ 3 (Figure 2).Function h is defined as follows: where α P and σ p are constants.
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 f s based on the flow rule regulated by potential function g s .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 f t = 0 based on the flow rule regulated by potential function g t .
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.

Stress Correction under the Plastic State. Shear yield should be considered. The flow rule can be written as
where λ 3 is the unknown variable that should be determined.
On the basis of the definition of function g s , the calculation of partial derivative expresses that 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 Superscripts N and O are the current and previous time step stress states.They are defined as 3 Complexity By combining ( 14) and ( 13), we obtain

17
Superscript I is the stress of the previous time step that is predicted based on the elasticity rule plus the stress increment.Therefore, The unknown parameter λ f can be calculated by positioning the current stress state on the shearing yield plane.In the yield function f 3 = 0, σ 1 and σ 3 are replaced by σ N 1 and σ N 3 .The transformation shows that In the following text, the tensile yield is considered.The flow rule of tension is expressed as follows: where λ 3 is an unknown variable.On the basis of the definition of g s , the calculation on partial derivative shows that Δe

Complexity
The Yellow River       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 FLAC 3D , the tensile strength of materials is defaulted to zero.If the given tensile strength is higher than σ t max , the tensile strength in the program is σ t max .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.

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 (Figure 3).This modeling zone includes different geological strata and three faults (f11, f12, and f14) that cover a total length of 120 m.
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 Figure 4.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.The 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  6 Complexity data, the velocity wave is selected as the applied seismic wave (Figure 8).

Calculation Parameters.
On the basis of the provided data, the following parameters are selected as material parameters in this calculation (Table.1).
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 (Figure 9).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).

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.

Stress and Displacement Analysis on Section A-A
during the First Excavation.As shown in Figures 14-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.00mm, and the resultant displacement is mainly in the range of 20.00~40.00mm.                      13 Complexity 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.00mm.
The distribution law of the major principal stress is shown in Figure 29.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.40MPa.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.30, 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.
As shown in Figure 33 As shown in Figures 34-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   Figure 31: Time-history response curve of velocity at characteristic points on the embankment in the second excavation.16 Complexity are above the excavation tunnel and the representative points P4, P5, and P6 are on one side of the excavation tunnel.

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.97mm, and the maximum displacement along Z ranges between −19.12 and 5.23 mm.The    17 Complexity variation range of the maximum resultant displacement is in the range of 2.27~19.93mm.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.
above deduction process of shear stress correction, we obtain

Figure 2 :
Figure 2: Plastic flow rule in the Mohr-Coulomb model.

Figure 3 :
Figure 3: Position sketch of the 3D calculation model.

Figure 9 :
Figure 9: Segmented excavation and zoning of blasting seismic wave application.

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 (Figure 12).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 13 :Figure 12 :
Figure 13: Key points on the embankment surface.

Figure 10 :
Figure 10: First excavation and blasting seismic wave application.Figure11: Fourth excavation and blasting seismic wave application.

Figure 11 :Figure 14 :
Figure 10: First excavation and blasting seismic wave application.Figure11: Fourth excavation and blasting seismic wave application.

Figure 15 :Figure 16 :
Figure 15: Displacement of section A-A along Z.

Figure 17 :
Figure 17: Major principal stress on section A-A.

Figure 18 :
Figure 18: Displacement of section B-B along X.

Figure 19 :
Figure 19: Displacement of section B-B along Z.

Figure 20 :
Figure 20: Resultant displacement of section B-B.

Figure 21 :
Figure 21: Major principal stress on section B-B.

Figure 22 :
Figure 22: Displacement of section C-C along X.

Figure 23 :
Figure 23: Displacement of section C-C along Z.

Figure 24 :
Figure 24: Resultant displacement of section C-C.

Figure 26 :
Figure 26: Displacement of section D-D along X.

Figure 27 :
Figure 27: Displacement of section D-D along Z.

Figure 28 :
Figure 28: Resultant displacement of section D-D.

Figure 29 :
Figure 29: Major principal stress on section D-D.

Figure 30 :
Figure 30: Time-history response curve of velocity at characteristic points on the embankment in the first excavation.

Figure 32 :
Figure 32: Time-history response curve of velocity at characteristic points on the embankment in the third excavation.

Figure 33 :
Figure 33: Time-history response curve of velocity at characteristic points on the embankment in the fourth excavation.

Figure 34 :
Figure 34: Time-history response curve of resultant displacement at characteristic points on the embankment in the first excavation.

Figure 35 :
Figure 35: Time-history response curve of resultant displacement at characteristic points on the embankment in the second excavation.

Figure 36 :
Figure 36: Time-history response curve of resultant displacement at characteristic points on the embankment in the third excavation.

Figure 37 :
Figure 37: Time-history response curve of resultant displacement at characteristic points on the embankment in the fourth excavation.
, 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.As shown in Figures 30-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.