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.
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 (, , and ) 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).
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 and . and 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 () (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, is the cohesive force, and 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 .
The potential function of shearing stress 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 . represents the bisector of shearing stress yield function and the tensile stress yield function on plane (Figure 2). Function is defined as follows: where and are constants.
The failure state of one point in plane can be determined by function 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 , 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 based on the flow rule regulated by potential function . 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 based on the flow rule regulated by potential function .
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 where is the unknown variable that should be determined. On the basis of the definition of function , 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 and are the current and previous time step stress states. They are defined as
By combining (14) and (13), we obtain
Superscript is the stress of the previous time step that is predicted based on the elasticity rule plus the stress increment. Therefore,
The unknown parameter can be calculated by positioning the current stress state on the shearing yield plane. In the yield function , and are replaced by and . The transformation shows that
In the following text, the tensile yield is considered. The flow rule of tension is expressed as follows: where is an unknown variable. On the basis of the definition of , the calculation on partial derivative shows that
Similar to the above deduction process of shear stress correction, we obtain where
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 () is acquired by adding the calculated stress increment on the stress of the previous time step. The principle budget stresses (, , and ) 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, and ((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 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 , the tensile strength in the program is . If the maximum principal stress that is calculated by one unit () is higher than tensile strength , 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 (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.
3.2. Calculation Loads
Similar to plane computation, three major types of loads are used: (1)Gravity loads(2)Water loads(3)Blasting seismic waves
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 data, the velocity wave is selected as the applied seismic wave (Figure 8).
3.3. 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).
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 (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.
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 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.00 mm, and the resultant displacement is mainly in the range of 20.00~40.00 mm.
The distribution law of the major principal stress is shown in Figure 17. 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.
4.2.2. Stress and Displacement Analysis on Section B-B during the Second Excavation
As shown in Figures 18–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.
The distribution law of the major principal stress is shown in Figure 21. 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.
4.2.3. Stress and Displacement Analysis on Section C-C during the Third Excavation
As shown in Figures 22–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.
The distribution law of the major principal stress is shown in Figure 25. 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.
4.2.4. Stress and Displacement Analysis on Section D-D during the Fourth Excavation
As shown in Figures 26–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.
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.40 MPa.
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 Figure 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 31, 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 32 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.
As shown in Figure 33, 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.
4.3.2. Time-History Response Analysis of Resultant Displacement at Characteristic Points under Different Excavation Steps
As shown in Figure 34, 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.
As shown in Figure 35, 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.
As shown in Figure 36, 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.
As shown in Figure 37, 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.
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 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.
Data Availability
The data used to support the findings of this study are currently under embargo while the research findings are commercialized. Requests for data, 12 months after publication of this article, will be considered by the corresponding author.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the Scientific Research Funds of Huaqiao University, the National Natural Science Foundation of China (51679093), the Open Research Fund of State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, China (SKLGDUEK1701), and the Promotion Program for Young and Middle-aged Teacher in Science and Technology Research of Huaqiao University, China (ZQN-PY112).