Abstract
Tunnel blast-induced vibration probably causes damage to the rock mass surrounding the tunnel surface and also to the rock mass of the slope at the tunnel entrance. It is important to simultaneously monitor the vibration on the tunnel surface and on the tunnel entrance slope face, especially when the blasting work face is close to tunnel entrance. During blasting excavation of the traffic tunnel at Baihetan hydropower station, vibration monitors were installed both on the tunnel surface and on the tunnel entrance slope face. Based on the monitoring data, a comparative study is conducted on the peak particle velocity (PPV) and frequency characteristics of the vibrations at these two locations. A three-dimensional FEM simulation of the tunnel blast is then performed to verify the field test results. The field monitoring and the numerical simulation show that there is significant difference between the vibration on the tunnel surface and that on the tunnel entrance slope face. The vibration on the tunnel surface has greater PPV and faster attenuation, while the tunnel entrance slope face has higher frequency and faster reduction rate in the center frequency. These differences are attributed to the different wave types and wave propagation paths. The tunnel surface is mainly surface waves transmitted through the damaged rock mass around the tunnel profile, while the tunnel entrance slope face originates mainly from the body waves transmitted via the undamaged rock mass inside the mountain. The comparisons of the monitored vibrations with the velocity limits specified in the Chinese standard show that the most dangerous vibration that may exceed the limit occurs on the tunnel surface. Therefore, the maximum charge weight used in the tunnel blast is determined by the vibration on the tunnel surface. Under different control standards, the allowable maximum charge weight per delay is further discussed.
1. Introduction
Construction of hydropower stations in the areas of high mountains and deep canyons involves large-scale excavation of underground caverns, such as powerhouses, diversion tunnels, tailrace tunnels, and traffic tunnels. The drilling and blasting method is the most commonly adopted for underground cavern excavation [1]. When explosives in blastholes are detonated, a part of the explosion energy is used for rock fragmentation, and the rest of the energy is dissipated in the form of seismic waves (blasting vibration), air blasts, flying rocks, and noise. Among these negative effects, the blasting vibration is a major concern for designers and constructors, because it will cause the maximum damage to the surrounding structures [2]. With regard to tunnel blasts in mountains, the blasting vibration will not only endanger the stability of the tunnel itself, but also adversely affect the stability of the slops at the tunnel entrances.
Blasting vibration velocity or acceleration monitoring is one of the most direct and effective methods to evaluate the negative effects of blasting vibration. Ramulu et al. [3] assessed the rock damage due to repeated blasting in a railway tunnel through borehole imaging observation and vibration velocity monitoring. Yu et al. [4] studied the safety criterion of blasting vibration on the tunnel wall based on vibration velocity and acceleration monitoring as well as 3D numerical simulations. During the blasting vibration monitoring for tunnel blasts, the monitors are typically arranged on the tunnel surface to measure the vibration responses of the surrounding rock masses or linings. The vibration monitoring for the slopes at the tunnel entrances is less involved. In fact, under blasting vibration, the rock slopes are also easily subjected to instable failure, especially in the presence of rock discontinuities [5, 6]. In some blasting vibration standards, the allowable peak particle velocity (PPV) for rock slopes is even smaller than that for traffic and mine tunnels [7]. However, because the vibration monitoring on the tunnel entrance slop face attracts much less attention, the vibration characteristics and dynamic responses of the tunnel entrance slopes under tunnel blasts are not very clear at present.
During tunnel blasts, the vibration monitors arranged on the tunnel surface mainly receive the surface waves coming through the tunnel wall [8], as shown in Figure 1. On the other hand, for the monitors on the tunnel entrance slope face, the received waves include not only the surface waves spreading via the tunnel surface, but also the body waves passing through the mountain and their reflected waves on the free surface. The body waves contain compressional P-waves and transverse S-waves. The surface waves include Rayleigh waves (R-waves) and Love waves (Q-waves). R-waves are generated whenever a free surface exists, but Q-waves are observed only when a soft superficial layer covers a stiffer medium [9]. Gutowski and Dym [10] studied the propagation velocity and energy transmission of the different types of waves (excluding the Q-wave). In the case of a point source, the P-wave spreads fastest but carries the least amount of the seismic energy, only accounting for about 7%. The S-wave transmits slower than the P-wave but faster than the R-wave and generally carries about 26% of the total seismic energy. The R-wave is the slowest among the three types of waves but carries the most amount of energy, about 67% [10]. The attenuation of the Rayleigh surface wave with an increase in distance is slower than that of the body wave [11]. With regard to the body wave radiation from a cylindrical blasthole, the analytical study shows that a significant amount of energy is transported in the shear-vertical (SV) mode [12].
Due to the different types of waves, blast-induced vibration inside the rock medium differs from that on the rock surface. In this regard, Wu et al. [13] and Wu and Hao [14] experimentally and numerically studied the characteristics of the vibration inside the geological medium and the vibration on the ground surface caused by underground blasts. They found that because of the wave refection on the ground surface, the PPV and principal frequency of the vibration on the ground surface are higher than those of the vibration inside the geological medium at the same distance. Based on the field monitoring, Shi et al. [1] investigated the difference between the vibration on the ground surface and the vibration in the underground tunnel from an open-pit bench blast. The results also show that the surface PPV is higher than the underground value at the same distance. Yang et al. [8] made a comparative study of tunnel blast-induced vibration on tunnel surfaces and inside surrounding rock. The field and numerical results show that, compared with the inside vibration, the tunnel surface vibration has a higher, more readily attenuated PPV and a lower frequency with a slower rate of decline in the dominant frequency. During the tunnel excavation as illustrated in Figure 1, the blast-induced vibrations on the tunnel surface and on the tunnel entrance slope face have different types of waves. Therefore, similar to the above cases, the vibration characteristics on these two sites are definitely different. However, such a comparative study has not been reported in the open literature. It is well known that, in blasting design, the maximum charge weight per delay is determined by the most dangerous vibration that may exceed the vibration limit specified in the standard. With regard to tunnel blasts near the entrance, the most dangerous vibration may occur on the tunnel surface or on the tunnel entrance slope face. Therefore, comparing the blasting vibration on the tunnel surface with that on the tunnel entrance slope face is of practical significance for the blasting design of tunnel excavation.
The paper is organized as follows: Section 2 introduces the blasting design and vibration monitoring of a traffic tunnel blast at the Baihetan hydropower station in China. In this project, vibration monitors were arranged both on the tunnel surface and on the tunnel entrance slope face. Section 3 presents the PPV and frequency characteristics of the tunnel blast-induced vibrations on the tunnel surface and the tunnel entrance slope face based on the field monitoring data analysis. A three-dimension finite element simulation of the tunnel blast-induced rock vibrations is presented in Section 4 to verify the field results. Section 5 focuses on discussing the maximum charge weight per delay of the tunnel blast design according to the vibrations on the tunnel surface and the tunnel entrance slope face. Section 6 concludes and summarizes the manuscript.
2. Blasting Design and Vibration Monitoring
2.1. General Overview
Baihetan hydropower station is the largest hydropower station currently under construction in the world. It is located on the lower reach of Jinsha River in China. This hydropower station mainly consists of a double-curvature arch dam, power houses, headrace structures, and flood discharge structures. Due to the limited site conditions in the high mountain and deep canyon area, the power houses, diversion tunnels, tailrace tunnels, and flood discharge tunnels are all built underground in the mountain. To enter these houses and tunnels for construction and operation, some traffic tunnels also need to be excavated in the mountain. The field test investigated in this study was carried out in the No. 2 traffic tunnel on the right bank, as shown in Figure 2. The No. 2 traffic tunnel is located at the elevation of 850 m and passes through the relatively complete cryptocrystalline basalt. According to the laboratory and field tests, the rock density is 2700−2800 kg/m3, the rock unconfined compressive strength is 70−95 MPa, the Young modulus of the rock mass is 15−18 GPa, and the acoustic P-wave velocity of the rock mass is in the range of 4200−4700 m/s. During the field test, the distance from the blasting work face to the tunnel entrance is approximately 100 m. The slope at the No. 2 traffic tunnel entrance has an angle of 28−63°.
2.2. Blasting Parameters and Vibration Monitoring
The drilling and blasting method is used in the excavation of the No. 2 traffic tunnel. The maximum excavation height is 7 m, and the excavation width is 8 m. In order to reduce the dynamic disturbance to the surrounding structures, the full-face millisecond delay blasting is adopted in the tunnel excavation. The cutting blastholes are first detonated, followed in sequence by the stoping blastholes and contour blastholes. These blastholes are divided into six delays with time intervals of 90–150 ms, as shown in Figure 3. The used blastholes are 50 mm in diameter and 3.0 m in length. Decoupled charges are used in the blasts. The explosive columns filled in the cutting blastholes and stoping blastholes have a diameter of 32 mm, and the explosive columns in the contour blastholes are 25 mm in diameter. The used emulsion-type explosive has a density of 1080–1400 kg/m3 and a velocity of detonation in the range of 4500–5500 m/s.
During the blasting excavation in the No. 2 traffic tunnel, velocity sensors are arranged on the tunnel surface to monitor the tunnel vibration and also placed on the tunnel entrance slope face to monitor the slope vibration, as shown in Figure 4. Inside the tunnel, there are seven monitoring points arranged on the floor along the tunnel axis. The horizontal distance between the monitoring points and the tunnel wall is about 1.5 m, and the distance to the blasting work face varies from 30 to 88 m. There are five monitoring points arranged on the tunnel entrance slope face at different elevations. These external monitoring points are located directly above the tunnel. Their horizontal distance to the blasting work face varies from 32 to 100 m, and the vertical distance to the tunnel roof changes from 48 to 20 m. At each monitoring point, a triaxial vibration velocity sensor is installed to monitor the longitudinal, transverse, and vertical velocity histories. The monitored signals are then transmitted to the data logger for extraction and analysis. Through the data logger, the sampling rate of the blasting vibration monitoring can be set. In order to retain high-frequency waves, the sampling rate is set to 8 kHz in this monitoring. The vibration velocity histories measured in two consecutive blasting tests are collected and analyzed in this study.
The cutting blastholes are detonated under the condition with only one free surface, i.e., the blasting work face. The detonation of the cutting holes creates new free surfaces for the detonation of the stoping holes. Similarly, the stoping hole blast also generates new free surfaces for the contour hole blast. The stoping holes and contour holes are both detonated under the condition with multiple free surfaces. Due to the lack of sufficient free surfaces, the cutting hole blast often produces the maximum vibration velocity among all of these delays. The vibration velocity histories monitored at this site also demonstrate it. Because of this, only the vibration caused by the cutting hole blast in the first delay is investigated. In fact, it makes little sense to incorporate the vibrations produced in different delays together for statistical analysis because the blastholes in different delays have different arrangements, charge structures, and free surface conditions.
3. Comparisons Based on the Field Monitoring
3.1. Velocity Histories and PPV
Figure 5(a) presents the vibration velocity histories measured at the No. 6 monitoring point on the tunnel surface. The distance from the No. 6 monitoring point to the blasting work face is approximately 80 m. On the tunnel floor, the vertical vibration velocity is the maximum velocity component. Figure 5(b) gives the vibration velocity histories recorded at the No. 9 monitoring point on the tunnel entrance slope face. The slant distance between this point and the blasting center is about 76 m. On the tunnel entrance slope face, the maximum vibration velocity component also occurs in the vertical direction. It agrees with the experimental and numerical results concluded by other researchers that the vibration velocity perpendicular to the free surface is often the maximum component [13, 14]. The vertical PPV at the No. 6 monitoring point on the tunnel surface is 0.8 cm/s, which is greater than the vertical PPV of 0.4 cm/s at the No. 9 monitoring point on the tunnel entrance slope face. The other two velocity components at the No. 6 monitoring point are also greater compared with the No. 9 monitoring point. The distance from the No. 6 monitoring point to the blasting work face is almost the same as that of the No. 9 monitoring point, or even farther. According to the comparison at these two points, it is seen that the blasting vibration velocity on the tunnel surface is greater than that on the tunnel entrance slope face at the same distance.
(a)
(b)
In order to further demonstrate the above result, the PPV variation with distance for the blasting vibrations on the tunnel surface and on the slope face is investigated based on the field monitoring data. Compared with the PPV of the maximum vibration velocity component, the PPV of the vector velocity provides a more accurate measurement for the real vibration level [15]. In most blasting vibrations, the vector PPV occurs at the same time as the maximum component PPV. However, the vector PPV is generally greater in magnitude due to the contribution of the other two velocity components. In this study, the vector PPV is used to represent the vibration magnitude. In the vector summation, the resultant velocity history is computed as follows:where is the resultant velocity history and , , and are the longitudinal, transverse, and vertical velocity histories.
The blasting vibration magnitude is related to many factors, such as blasthole arrangements, explosive types, charge structures, geological conditions, rock properties, and wave propagation distance [4]. Therefore, it is very difficult to accurately predict the blasting vibration magnitude, even for a single blasthole detonated in a homogeneous medium. Alternatively, empirical charge weight scaling laws are usually employed to describe the PPV attenuation relation with increasing distance. In this study, it is considered that the vector PPV follows the most typical charge weight scaling law:where d is the distance between the monitoring point and the blasting source (m), W is the charge weight per delay (kg), and K and α are site constants. For the cutting hole blast in the field test, the charge weight W is 24 kg. The ratio of the charge length (3 m) to the charge diameter (0.032 m) is about 93, which is a long cylindrical charge. Therefore, the distance is scaled with the square root of the charge weight [16].
It should be noted that the scaling law in (2) is not a fundamental equation for blasting vibration prediction. It only provides a pragmatic assessment of the blasting vibration magnitude for a given charge weight and distance. Since the scaling law has no real physical meaning, the slant distance to the blasting source is adopted in (2) for the monitoring points on the tunnel entrance slope face. Under this scenario, the PPV attenuation along the horizontal distance and the vertical distance is not considered separately, but incorporated into the attenuation along the slant distance. According to the field monitoring data, the vector PPV versus scaled distance on a log-log plot is given in Figure 6. The linear regression analysis is then conducted on these scattered data to yield the site constants K and α. The intercept of the best-fitting line at is the value of logK, and the slope is the value of α. After the site constants K and α are determined, the vector PPV variation with scaled distance can be obtained as follows. For the blasting vibration on the tunnel surface,
For the blasting vibration on the tunnel entrance slope face,
According to the charge weight scaling equations in (3) and (4), the PPV attenuation relations with increasing distance at W = 24 kg are shown in Figure 7. It is clearly seen that, within the distance of 90 m, the PPV of the blasting vibration on the tunnel surface is greater than that on the tunnel entrance slope face. Beyond 90 m, the PPV on the tunnel surface is smaller because it decays faster with an increscent distance. Actually, the blasting vibration in the far field beyond 90 m is not a concern of interest in this study because the PPV at this location is much smaller than the velocity limits of rock damage to the tunnel or slope. Therefore, in the region of our attention, the blasting vibration on the tunnel surface exceeds that on the tunnel entrance slope face. The above charge weight scaling laws in (3) and (4) are derived based on the 10–14 sets of scattered data during the two blasting tests. Therefore, an appropriate statistical test needs to be conducted on the scattered data to show that there is a statistical difference between the PPV on the tunnel surface and the PPV on the slope face. To achieve this, hypothesis tests are performed on the scattered data shown in Figure 6 by using the statistical package SPSS. The hypothesis tests show that, at the significance level of 0.05, the slopes and intercepts of the two fitted lines shown in Figure 6 are significantly different from each other. This means that there is a significant difference between the two fitted lines of the scattered data. Therefore, in the statistical sense, the PPV of the blasting vibration on the tunnel surface is significantly different and furthermore greater than the PPV on the tunnel entrance slope face.
As mentioned earlier, it is the different types of waves that cause the difference between the vibration on the tunnel surface and the vibration on the tunnel entrance slope face. The vibration waveforms monitored on the tunnel surface are mainly the surface waves transmitted from the blasting source via the tunnel floor, while the vibration waveforms recorded on the tunnel entrance slope face are mainly composed of the body waves spreading directly through the mountain and their reflection waves generated on the slope face. For the body waves at nonvertical incidence on the slope face, mode-converted body waves and surface waves are generated [12]. It is well known that when waves propagate to a free surface, the magnitude of the waves tends to be amplified due to the wave reflection [13, 14]. Even in the presence of this amplification effect, the vibration on the tunnel entrance slope face is still smaller than that on the tunnel surface. This is because the body waves carry only a small amount of explosion seismic wave energy from the blasting source, and the majority of the energy is transmitted outside in the form of surface waves [9, 10]. For the blasts of the traffic tunnel, the cylindrical explosives are detonated at an extremely shallow depth, within 3.0 m beneath the blasting work face. In this situation, even more energy is transmitted through the tunnel surface in surface waves.
It is seen from Figure 7 that the PPV of the blasting vibration on the tunnel surface decays faster with distance in comparison with the PPV on the tunnel entrance slope face. The attenuation of blasting vibration is related not only to the wave types, but also to the geological conditions and rock properties on the wave propagation path. For the point source in a homogeneous medium, the body wave is attenuated with an increase in distance (D) by D−1, and the Rayleigh surface wave is attenuated by D−1/2 [11]. In the simplest model, the attenuation of the surface wave is slower than that of the body wave. It is contrary to the field monitoring result that the surface waves measured on the tunnel surface decay faster. This indicates that, in the field test, the attenuation of blasting vibration is dominated by the geological conditions and rock properties on the wave propagation path. For the blasting vibration on the tunnel surface, the wave propagation path is the rock mass immediately surrounding the tunnel wall and floor from the blasting work face. On the other hand, for the blasting vibration on the tunnel entrance slope face, the waves reach the slope face through the rock mass inside the mountain. The No. 2 traffic tunnel is excavated in relatively complete basalt, and the geological conditions around the tunnel and inside the mountain are not much different. However, the properties of the rock mass around the tunnel profile are different significantly from those of the rock mass inside the mountain due to the blasting excavation. Inside the mountain, the rock mass is undamaged, whereas the stress waves and high-pressure gases produced by blasts cause cracks to extend into the rock mass surrounding the tunnel profile. These cracks can persist into a depth of 2-3 m beyond the tunnel surface [17]. Due to existence of these blast-induced cracks, the integrity of the rock mass around the tunnel surface is reduced compared with the undamaged rock mass. Some acoustic testing results have shown that, during tunneling blasts, the P-wave velocity of the damaged rock mass around the tunnel profile is reduced by 40%–50% [8]. It is well known that the damaged rock mass with lower integrity transmits waves at faster attenuation. Therefore, the blasting vibration on the tunnel surface that is transmitted through the damaged rock mass decays faster.
3.2. Amplitude-Frequency Spectra and Center Frequency
The structure damage under blasting vibration is dependent not only on the PPV, but also on the vibration frequency. When the vibration frequency is approaching or equal to the natural frequency of the structure, the vibration magnitude on the structure will be amplified. The natural frequencies of most structures are lower. Therefore, under the same PPV, the low-frequency vibration has greater potential to cause structural damage than does the high-frequency vibration. Because of this, the frequency content is considered in most of the current blasting vibration standards [18]. In these standards, the structure damage criterion is given as a PPV limit based upon predominant frequency. With regard to tunnels and rock slopes, the Chinese standard from GB6722-2014 is tabulated in Table 1. It is seen that, for the blasting vibration at higher frequency, a greater PPV is allowable, and the allowable PPV is reduced as the frequency decreases. Therefore, extra care should be taken for the low-frequency vibration in the far field, where it probably exceeds the PPV limit at low frequency.
By performing the Fourier transformation on the velocity histories plotted in Figure 5, the amplitude-frequency spectra of the vibrations on the tunnel surface and slope face are obtained, as shown in Figure 8. To facilitate comparison, the normalized amplitude A/Amax is presented on the ordinate axis. For the vibration at the No. 9 monitoring point on the tunnel entrance slope face, the frequency is distributed in the band of 0−600 Hz. For the vibration at the No. 6 monitoring point on the tunnel surface, the amplitude of the spectra has almost dropped to zero at about 300 Hz. In the frequency band of 100−300 Hz, the amplitude of the vibration on the slope face is significantly greater than that on the tunnel surface. The direct comparison of the amplitude-frequency spectra between the two points shows that the vibration on the tunnel entrance slope face has a wider frequency band and more high-frequency content.
(a)
(b)
(c)
In order to quantify conveniently, several average parameters, such as dominant frequency, center frequency, and spectral bandwidth are used to characterize the amplitude-frequency spectrum. In this study, the center frequency is adopted as the characteristic frequency of the amplitude-frequency spectrum. It is defined aswhere fc is the center frequency, fi is the individual frequency in the spectrum, and A (fi) is the amplitude associated with each frequency fi.
For the blasting vibration at the No. 6 monitoring point on the tunnel surface, the center frequencies of the longitudinal, transverse, and vertical velocity components are 119 Hz, 124.8 Hz, and 122.8 Hz, respectively. With regard to the vibration at the No. 9 monitoring point on the slope face, the center frequencies are correspondingly 152.4 Hz, 142.7 Hz, and 136.7 Hz. For these two monitoring points at almost equal distance, the frequency of the vibration on the tunnel entrance slope face is higher than that of the vibration on the tunnel surface. It is assumed that the center frequency fc also follows the charge weight scaling law [8, 13]. For the vibrations monitored at different distances, the center frequency versus scaled distance is shown in Figure 9. The center frequency presented in this figure is the average of the longitudinal, transverse, and vertical components. Similarly, the linear regression is performed on the scattered data to give the center frequency variation with scaled distance. For the blasting vibration on the tunnel surface,
With regard to the vibration on the tunnel entrance slope face,
The hypothesis tests on the scattered data show that there is a significant difference between the two fitted lines presented in Figure 9. This indicates that the center frequency of the vibration on the tunnel surface is statistically different from that on the slope face. According to the above scaling laws in (6) and (7), the center frequency of the vibration at varying distance is shown in Figure 10. Within the distance of 100 m, the center frequency of the vibration on the tunnel entrance slope face is significantly higher than that on the tunnel surface. However, the center frequency of the vibration on the slope face declines faster as the distance increases. Exceeding the 100 m distance, it gradually approaches the center frequency of the vibration on the tunnel surface. However, anyway, in the distance of interest, the blasting vibration on the tunnel entrance slope face has higher frequency. This is because the vibration on the slope face originates mainly from the body waves that travel in the undamaged rock mass inside the mountain, while the vibration on the tunnel surface is the surface waves that spread in the damaged rock mass surrounding the tunnel profile. Under a given input source, a relatively competent transmission medium attenuates the amplitude at high frequency more slowly compared with a fissured medium and thus corresponds to a higher center frequency.
During explosion seismic wave propagation in a viscoelastic medium, the seismic Q attenuates the amplitude of the spectrum as the propagation distance increases. However, the amplitudes at different frequencies are not equally attenuated. The amplitude at high frequency is attenuated faster than that at low frequency with an increase in the distance [19]. Consequently, the relative frequency content of the amplitude-frequency spectrum is changed and the center frequency is reduced as the distance increases. According to this mechanism, the center frequency declines faster with distance for the waves carrying more frequency content. Therefore, the vibration on the tunnel entrance slope face that has higher frequency corresponds to a faster reduction in the center frequency.
The above results are obtained based on a small amount of field monitoring data, and furthermore there is a scatter for the field data. This kind of scatter is probably caused by some complicated and uncertain factors, such as rock discontinuities on the vibration propagation paths, rock looseness at the monitoring locations, and instrument differences at each monitoring point. To compensate for this deficiency, numerical simulation of the tunnel blast is conducted to demonstrate the field monitoring results. In the numerical model, these uncertain factors that influence the blasting vibration can be avoided.
4. Comparisons Based on Numerical Simulation
4.1. Numerical Model
Rock blasting is an important class of fluid-structure interaction (FSI) problems. One difficulty in handling numerically the FSI stems from the fact that the structural equations are usually formulated with material (Lagrangian) coordinates, while the fluid equations are typically written by using spatial (Eulerian) coordinates. A straightforward approach to the solution of the coupled fluid-structure dynamic equations requires moving at each time step, at least the portions of the fluid grid that are close to the moving and flexing structure grid. This approach is appropriate for small deformation but may lead to severe grid distortions when the structure undergoes large deformation [20]. Some approaches have emerged as an alternative to handle the FSI computations. Donea et al. [21] presented arbitrary Lagrangian–Eulerian (ALE) finite element methods for the prediction of the nonlinear responses of fluid-structure systems exposed to transient dynamic loading. In the ALE formulation, the fluid domain is formulated by an ALE kinematical description in which the grid points can be moved with the fluid in a normal Lagrangian fashion, be held fixed in an Eulerian manner, or be moved in some arbitrarily specified ways to give a continuous rezoning capability [21]. Recently, smoothed finite element methods (S-FEMs) based on adaptive mesh are also applied to the FSI computations to reduce overall cost [22]. Among these S-FEMs, the edge-based S-FEM (ES-FEM) is found to be the most computationally efficient [23, 24]. With regard to meshfree methods, Rabczuk and Eibl [25] and Fakhimi and Lanari [26] described the fluid interactions with large motion of structures by using a smooth particle hydrodynamics (SPH) method. Rabczuk and Belytschko [27] and Rabczuk et al. [28] also proposed a cracking-particle approach for modeling discrete cracks of fracturing structures under impulsive loads. This meshfree method is simple and does not require any modifications when structures fail and allow fluid to flow through cracks.
The numerical modeling of the tunnel blast-induced rock vibration in this study is implemented by using the commercial finite element software ANSYS/LS-DYNA. It is well known that the dynamic program LS-DYNA is one of the few codes that can simulate the explosive detonations and the interactions between detonation products (fluids) and structures by using an ALE algorithm. In this program, the simulation of explosive detonations is implemented by using the Jones–Wilkins–Lee (JWL) equation of state (EOS) [27, 28]. Accuracy of the ALE modeling is highly dependent on mesh refinement. Therefore, it is very computationally expensive to use the ALE algorithm to simulate the explosive-rock interactions. Because of this, the ALE approach is usually used to handle the near-field rock fragmentation in simple cases, such as two-dimensional or single-hole problems. The far-field rock responses under practical three-dimensional and multihole blasts are rarely simulated by using the ALE approach due to low computational efficiency. The interest of this numerical study is the far-field rock vibration caused by practical tunnel blasts involving multiple blastholes. Therefore, in order to improve computational efficiency, the FSI is not handled in the numerical modeling. Instead, the blast loading process is treated as a pressure history acting on the blast-created free surfaces (blasting excavation boundaries).
Consistent with the above field data analysis, only the cutting hole blast is simulated in this study. Based on the above assumptions, a numerical model is developed in the Lagrangian coordinate, in which the cutting blast-created free surfaces are the blast loading boundaries. The model measures 150 m long, 50 m wide, and 75 m high, as shown in Figure 11. The traffic tunnel is located 20 m above the bottom of the model. In this tunnel, the rock mass in the first 100 m has been excavated, and the distance from the tunnel entrance to the blasting work face is 100 m. The slope at the tunnel entrance is simplified to one with four levels according to the field geography. From the top down, each level has a height of 8 m, 10 m, 10 m, and 20 m, and the slope angles are 28°, 34°, 63°, and 39°, respectively.
Ground vibration caused by blasting is attributed to elastic seismic wave propagation in media. The elastic seismic waves cannot directly cause damage to the rock medium. The interest of this study is the far-field vibration that is transmitted in the form of elastic seismic waves. Therefore, in this numerical study, the rock medium is considered to be a linear elastic material. It should be noted that the rock mass surrounding the tunnel profile is damaged by blasting during excavation. The properties of the damaged rock mass are inferior to those of the undamaged rock mass inside the mountain. In order to estimate the rock mass strength and deformation modulus under excavation conditions, a parameter called the disturbance factor is introduced into the Hoek–Brown failure criterion to deal with the strength and deformation modulus reduction due to blast damage and stress relaxation. With regard to the rock mass deformation modulus Em (GPa), it is given by [29]:where σci is the unconfined compressive strength of the intact rock material, GSI is the geological strength index that depends on rock structures and joint conditions, and D is the disturbance factor that represents the degree of disturbance or damage to the rock mass subjected to blasting and stress relaxation. The value of D varies from 0 for undamaged rock masses to 1.0 for severely damaged rock masses. For undamaged rock masses at D = 0, the deformation modulus Em0 is
Based on (8) and (9), the ratio of the deformation modulus between damaged rock masses and undamaged rock masses is then obtained as follows:
When the undamaged rock mass deformation modulus Em0 is known, the damaged deformation modulus Em can be estimated by
Guidelines for estimating the disturbance factor D due to blast damage and stress relaxation are given by Hoek and Brown [29]. For the severe blast damage extending 2.0–3.0 m into the surrounding rock mass of a tunnel, the disturbance factor D = 1.0 at the tunnel surface with a linear decrease to D = 0 at 2.0 m into the surrounding rock mass is assigned. For the linear elastic material, the deformation modulus Em is equal to Young’s modulus E. According to the above selection of the disturbance factor D, Young’s modulus of the damaged rock mass around the tunnel profile linearly increases from E0/2 on the tunnel surface to E0 at 2.0 m depth, where E0 is Young’s modulus of the undamaged rock mass. In addition to the rock mass surrounding the tunnel profile, the rock mass immediately behind the slope face is also damaged by blasting and stress relaxation. The field acoustic tests in combination with the Hoek–Brown criterion show that the damaged zone extends 2.0–3.2 m behind the slope face, and the most severe disturbance factor D = 1.0 is reached on the slope face [30]. Similarly, a linearly increasing Young’s modulus from E0/2 to E0 is assigned to the damaged rock mass within 3.2 m behind the slope face. The density and Poisson’s ratio are considered unchanged for the damaged rock mass. The energy dissipation in the dynamic system is implemented through setting damping. In this numerical simulation, the damping is reproduced by the classical Rayleigh approach. The damping matrix in the Rayleigh damping is defined as a linear combination of the mass matrix and the stiffness matrix. The mass and stiffness proportional damping constants are set to 6.0 and 0.05 in the present study.
During the tunnel blast, in addition to the dynamic blasting load, the rock mass is also subjected to static in situ stress due to gravity. Implementation of the numerical modeling of the tunnel blast-induced rock vibration involves two steps, static stress initialization and dynamic blast loading. This systematic process can be implemented in ANSYS/LS-DYNA by using its implicit and explicit solutions in sequence, as shown in Figure 12. The bottom and side boundaries are set to be fixed, and the gravity is first preloaded on the entire model. The implicit solver in ANSYS is then performed to compute the initial stress and deformation of the rock elements. This information is written into a database file and delivered to the explicit program LS-DYNA to update the stress and deformation of the explicit elements. After that, a blasting pressure history is applied on the explicit elements of the blasting excavation boundaries. Finally, the explicit solver in LS-DYNA is operated to computer the blast-induced rock vibration. In order to prevent wave reflection on the artificial finite boundaries, nonreflecting boundary conditions are exerted on the bottom and sides of the model.
4.2. Blast Loading
As mentioned above, a pressure history is used in this numerical study to represent the blast loading process. There are various functions attempting to approximate the pulsing pressure history on the blasthole wall, such as Gaussian functions, triangular load functions, and pressure decay functions. In this study, the decay function is employed. It is expressed as follows [31]:where is the pressure history applied on the blasthole wall, is the peak of the pressure history, β is a constant that depends on the rise time of the pressure history, and t is time.
For a decoupled cylindrical charge, the peak pressure on the blasthole wall can be estimated bywhere ρe is the explosive density, VOD is the velocity of detonation, dc is the charge diameter, db is the blasthole diameter, and γ is the adiabatic exponent of detonation gases. For the explosives commonly used in rock blasting, the adiabatic exponent γ approximates 1.5.
For the pressure history in (12), the peak pressure is attained at the time . Then, the constant β is given by
In (12), the angular frequency of the pressure wave ω and the constant β have the following relationship [31]:
Then, the frequency f in Hz can be written in terms of the constant β as
In this numerical study, the rise time tr = 110 μs is considered. According to (14)−(16), the rise time of 110 μs can yield the maximum frequency of 1000 Hz that is required. The pressure decay function in (12) is the blasting load acting on the blasthole wall. It needs to be equivalently transferred to the blasting excavation boundary to deal with simulation of the multihole blasts. When a row of blastholes is detonated precisely at the same time, blast-induced cracks grow preferentially along the connecting line between the holes. Therefore, for the cutting hole blast in Figure 3, the connecting lines between the cutting holes in the same delay can be considered as the blasting excavation boundary. The blasthole diameter is much smaller than the spacing between the adjacent blastholes, and hence the equivalent transference from the pressure on the blasthole wall to the pressure on the excavation boundary follows the relationshipwhere Pe (t) is the equivalent pressure history on the blasting excavation boundary and S is the spacing between the adjacent blastholes in the same delay.
For the cylindrical charges with a limited VOD, the detonation propagation within explosive columns has an important effect on seismic wave radiation [12]. For most explosives and rock types, S-Mach waves are always generated due to the detonation propagation. When the velocity of detonation propagation is high enough, P-Mach waves are also generated. If the velocity of detonation within explosive columns exceeds the P-wave velocity in rock media, the seismic wave radiation from a cylindrical charge is completely dominated by P-Mach and S-Mach waves. These Mach waves propagate out in cone shapes; i.e., they do not spread out in a radial direction. In order to mimic the detonation propagation, a traveling pressure at the velocity of detonation is applied on the blasting excavation boundary over the entire blasthole length of 3.0 m, as shown in Figure 11. Along the 3.0 m length, there are 12 elements with a uniform size of 0.25 m. Consequently, the delay interval for the start of pressure between each element is set to 59 μs, which is equal to the duration of the detonation wave propagation through an element.
4.3. Validation of the Numerical Modeling
In the numerical modeling of blast-induced rock vibration, the output of the velocity histories is highly dependent on the input parameters of rock mass properties and blast loading conditions. Our qualitative sensitivity analyses show that the PPV and frequency of the vibration velocity are mostly affected by Young’s modulus of the rock mass, and the peak pressure and rise time of the blasting load. As Young’s modulus of the rock mass increases, the PPV decreases but the frequency content increases. With an increase in the peak pressure of the blasting load, the PPV is significantly increased. When the rise time of blast loading becomes longer, the vibration frequency is significantly reduced. According to the test data of the rock mass properties that are given in the tunnel design document, the mean values are assigned to the undamaged rock mass inside the mountain as follows: rock mass density ρ = 2750 kg/m3, Young’s modulus E0 = 16.5 GPa, and Poisson’s ratio ν = 0.25. The mean explosive density ρe = 1240 kg/m3 and velocity of detonation VOD = 5000 m/s are adopted as provided by the explosive manufacturer. The other input parameters about the blasting geometry remain the same as the site survey. These parameters yield a peak pressure of 101.6 MPa for the blasting load on the blasting excavation boundary, as shown in Figure 11. Note that there are many uncertainties in the input parameters, particularly in the geomaterials like rock. In order to quantify the influence of the uncertain input parameters on the simulation results, it is necessary to carry out an uncertainty analysis. However, this content is beyond the scope of this paper. Nevertheless, a comparison between the simulated vibration waves and the field monitoring data is made in the following text to validate the numerical modeling.
The geometric model is discretized into a mesh of solid elements. The element size is determined by the wavelength of the wave passing through the element. Blair [12] recommended 6–12 elements per wavelength to avoid wave distortion. For the rock mass used in this numerical modeling, the Rayleigh wave velocity that is the lowest wave velocity is 1727 m/s. From the amplitude-frequency spectra presented in Figure 8, the highest frequency does not exceed 1000 Hz. Consequently, the shortest wavelength of interest is approximately 1.7 m. According to Blair’s recommendation [12], the minimum element size in this numerical modeling should be shorter than 0.28 m. As the distance to the blast loading boundary increases, the amplitudes at high frequencies decline. At far distance, most of the vibration amplitudes are distributed in the frequency bands lower than 500 Hz. Therefore, a gradient mesh that varies from 0.25 m near the blast loading boundary to a larger size on the borders is employed to adapt to the decreasing vibration frequency with distance. Under this restriction, structured discretizations with different numbers of hexahedral elements from 1,250,000 to 6,740,000 are first considered. The obtained vertical velocity histories at No. 6 observation point on the tunnel surface are shown in Figure 13. When the number of the structured discretization elements is increased from 1,250,000 to 6,740,000, the simulated waveforms are almost the same with only a slight difference of less than 6% in the PPV and center frequency. Then, an unstructured discretization arrangement with 6,270,000 tetrahedral elements is considered. The simulated result almost coincides with the results of the structured discretizations, as shown in Figure 13. However, the computational cost of the unstructured discretization is higher than that of the structured discretization at the approximate number of elements. Under the above-mentioned element size restriction, the vibration velocity-time curves for the different discretizations do not show any mesh dependence.
Another measure for checking the quality of the numerical modeling is to calculate the energy balance through [28]where Eint is the internal energy, Ekin is the kinetic energy, Edis is the dissipated energy due to the damping property of the rock material, ΔWTB is the external work done at NB boundary nodes at NT time steps, N is the number of nodes, t0 is the initial time, and tn is the duration of computation.
Figure 14(a) presents the energy-time histories for the structured discretization of 6,740,000 elements. When the blasting load acts on the excavation boundaries, the external work done by the blasting load is transformed into the internal energy and kinetic energy of the rock mass. As the rock mass rebounds in the falling stage of the blasting load and the vibration decays with time, most of the energy is dissipated in the form of heat. Only a small portion is stored in the rock mass in the form of internal energy. Figure 14(b) shows the relative errors in the energy balance for the different discretizations. The relative errors for the structured 6,740,000-element discretization and the unstructured 6,270,000-element discretization are less than 5%. Even for the coarser structured discretizations with 3,130,000 elements and 1,250,000 elements, the relative errors still do not exceed 10%. An acceptable energy balance is achieved for the different discretizations in the numerical modeling. From the errors in the energy balance, the most accurate results are obtained for the structured discretization of 6,740,000 elements. Therefore, this mesh discretization is finally adopted to simulate the tunnel blast-induced rock vibration. Figure 15 shows the simulated vertical velocity histories at No. 6 observation point on the tunnel surface and No. 9 observation point on the tunnel entrance slope face by using this mesh. The simulated vibration waves agree well with the field monitoring data. This indicates that the selected rock mass and explosive parameters and the used mesh discretization are appropriate.
(a)
(b)
(a)
(b)
The purpose of this numerical study is not to pursue the accurate agreement between the simulated vibration waves and the field monitoring data, but to verify the comparison relationship between the vibration on the tunnel surface and the vibration on the slope face. Therefore, some assumptions and simplifications are made in the numerical modeling to improve the computational efficiency. The treatment of the explosive detonation as a pressure history probably causes some deviations in the simulated results, particularly in the near field of the blasting source. Fortunately, many studies have shown that such a treatment can achieve satisfactory simulation accuracy for the far-field rock vibration [8, 32]. Therefore, the deviations due to the assumptions and simplifications of the numerical modeling do not change the relative relationship between the vibration on the tunnel surface and the vibration on the slope face.
4.4. Numerical Results
From the blasting work face to the tunnel entrance, 30 observation points are arranged on the tunnel floor along the tunnel axis to observe the blasting vibration on the tunnel surface. The longitudinal distance between the adjacent observation points is 2−5 m. The distance from the observation points to the side wall of the tunnel is also set to 1.5 m to keep it consistent with the field monitoring. There are also 30 observation points selected on the ground to record the vibration on the tunnel entrance slope face. These external observation points are all located exactly above the observation points inside the tunnel. According to the simulated vibration waveforms at various distances, the vector PPV and the center frequency against scaled distance are shown in Figures 16 and 17. By solving the slope and intercept of the best-fit lines on these discrete points, the charge weight scaling laws of the vector PPV and the center frequency are obtained as follows. For the vibration on the tunnel surface,
For the vibration on the tunnel entrance slope face,
It is seen from Figures 16 and 17 that most of the discrete points are clustered immediately near the fitted lines. Compared with the field monitoring data in Figures 6 and 9, a better fit is achieved on the simulated data because those uncertain factors are eliminated in the numerical simulation. The hypothesis tests on the simulated data in Figures 16 and 17 show that, at the significance level of 0.05, there is a significant difference between the two fitted lines. The numerical simulation results also prove that the vibration on the tunnel surface is obviously different from that on the tunnel entrance slope face in magnitude and frequency. From the comparison between (19) and (21), although the PPV on the tunnel surface decays faster with an increase in distance, it still exceeds the PPV on the slope face at the same distance within 100 m. The comparison between (20) and (22) shows that the vibration on the slope face has a higher and more readily reduced center frequency. All of these results based on the numerical simulation are consistent with the findings of the field monitoring. This verifies the vibration characteristics on the tunnel surface and tunnel entrance slope face obtained in this study.
The rock vibration response under blasting is a very complicated topic because it is influenced by many factors, involving explosive properties, blasthole arrangements, detonation sequences, detonator accuracy, and lithological and geological conditions. Among these factors, some are fully considered, some are simplified, and some minor ones are even completely ignored. Because of this, there is some deviation between the simulated vibration waveforms and the monitored waves. For example, the vector PPV and center frequency of the simulated vibration are higher. In fact, it is impossible to exactly match the simulated waves with the monitored ones. Notwithstanding the limitations of the numerical simulation, it still demonstrates the comparison between the vibration on the tunnel surface and the vibration on the tunnel entrance slope face.
5. Discussion
It is reasonable to expect that the blasting vibration magnitude is increased as the charge weight W increases. It is also reasonable to expect that the vibration magnitude is decreased as the distance d between the blasting source and the monitor increases. This is the basis of the charge weight scaling law in (2). According to the scaling laws in (3) and (4), the vector PPV on the tunnel surface and tunnel entrance slope face at various charge weights and distances can be estimated. Figure 18 shows the vector PPV variation with an increase in the charge weight for different blast-to-monitor distances. During the excavation of the No. 2 traffic tunnel, the vibration caused by blasting may endanger the tunnel itself or the slope at the tunnel entrance. It is well known that the structure damage under blasting vibration depends not only on the PPV but also on the frequency. The amplitude-frequency spectrum analysis of the monitored vibration waves shows that within 100 m distance, the center frequencies exceed 50 Hz, whether on the tunnel surface or on the slope face. According to the Chinese standard from GB6722-2014 in Table 1, the maximum allowable vibration velocities for the traffic tunnel and slope are 20 cm/s and 15 cm/s, respectively. The vibration velocity limits are also drawn in Figure 18. It is seen that the PPV of the vibration on the tunnel entrance slope face is situated far below the standard of 15 cm/s when the charge weight varies from 10 to 50 kg and the distance varies from 10 to 50 m. Under the same variation in the charge weight and distance, the PPV of the vibration on the tunnel surface is mostly below the control standard of 20 cm/s. However, at the distance d = 10 m, the PPV on the tunnel surface exceeds the velocity limit specified in the standard when the charge weight is more than 35 kg. This indicates that the rock mass surrounding the tunnel surface is more adversely affected by the tunnel blast vibration. Then, during the blasting excavation of the No. 2 traffic tunnel, the criterion of blasting vibration control is that the PPV on the tunnel surface must not exceed 20 cm/s.
The factors that influence the blasting vibration magnitude can be summarized into two categories. One is the blasting source factors, and the other one is the propagation path factors. Accordingly, blasting vibration control is generally implemented through two approaches. One is to reduce vibration intensity at the source, and the other one is to accelerate vibration attenuation on the wave propagation path. In practice, for example, presplitting cracks and manually excavated grooves are set on the propagation path to accelerate blasting vibration attenuation. However, setting additional presplitting cracks and grooves increases costs. Consequently, reducing the vibration intensity through adjusting source factors becomes the main means of blasting vibration control. The adjustable blasting source factors mainly include total charge weight, charge weight per delay, blasthole arrangements, charge structures, detonation sequences, and delay intervals. Among these factors, the charge weight per delay is the most related to the vibration magnitude and furthermore is a factor that is easier to adjust. Therefore, the change of the charge weight per delay is preferentially selected to control blasting vibration. When the charge weight scaling law for the PPV and the PPV limit are known, the allowable charge weight per delay can be estimated bywhere Wc is the allowable charge weight per delay and PPVc is the velocity limit specified in the standard.
As mentioned earlier, the allowable charge weight per delay is determined by the most dangerous vibration that may exceed the velocity limit specified in the standard. With regard to the blasting excavation of the No. 2 traffic tunnel, the most dangerous vibration occurs on the tunnel surface, and hence PPVc = 20 cm/s, K = 42.39, and α = 1.477 are considered. Substituting these parameters into (23) yields the allowable charge weight per delay for the protection targets at different distances. In order to protect the rock mass surrounding the tunnel surface at a distance of 10 m, the charge weight per delay must not exceed 36 kg. If a stricter control standard, for example, PPVc = 15 cm/s, is adopted, the allowable charge weight per delay should be limited to 25 kg. During the blasting excavation of the No. 2 traffic tunnel, the used maximum charge weight per delay is 24 kg. This indicates that the blasting design used in the traffic tunnel excavation is appropriate.
The dynamic stability of the tunnel entrance slope under the tunnel blast-induced vibration is analyzed by using the equivalent acceleration of blasting vibration and the Sarma method of the rigid body limit equilibrium analysis [33]. Because the direction of the blasting vibration inertial force changes with time due to the phase change of the vibration waves, the safety factor of the dynamic stability fluctuates around the static stability safety factor under gravity. It finally tends to the static value due to the attenuation of blasting vibration with time. Compared with the static value, the dynamic stability safety factor of the tunnel entrance slope is increased or decreased by less than 3.5% under the tunnel blast-induced vibration. The minimum value of the dynamic stability safety factor is 2.57, which is greater than the value specified in the slope design standard. Therefore, it is considered that the tunnel entrance slope is stable under the tunnel blast-induced vibration. The monitored vibration velocities on the tunnel entrance slope face, which do not exceed the PPV limits specified in the standard, also demonstrate the dynamic stability of the slope.
6. Conclusions
Based on the blasting vibration monitoring in the tunnel blast test, a comparative study was carried out on the vibration on the tunnel surface and the vibration on the slope face at the tunnel entrance. The statistical tests on the field monitoring data show that there is a significant difference between the vibrations at these two locations. The vibration on the tunnel surface is mainly the surface waves transmitted through the rock mass surrounding the tunnel profile. On the other hand, the vibration on the tunnel entrance slope face is mainly composed of the body waves transmitted through the rock mass inside the mountain and the reflection waves due to the body wave incidence on the free slope face. Because the majority of explosion seismic wave energy is transmitted outside in the form of surface waves, the PPV on the tunnel surface is greater than that on the slope face at the same distance. The rock mass around the tunnel profile is damaged by blasting, and it results in a faster PPV attenuation with distance for the vibration on the tunnel surface. The vibration on the slope face has a wider frequency band and a higher center frequency because the body waves travel in the undamaged rock mass inside the mountain. With an increase in the traveling distance, the amplitude at high frequency is attenuated faster than that at low frequency. Consequently, the vibration on the slope face that carries more high-frequency content corresponds to a faster reduction in the center frequency as the distance increases. The above conclusions obtained from the field monitoring data have been verified by a three-dimensional dynamic FEM simulation.
According to the charge weight scaling laws derived from the field monitoring data, the PPV on the tunnel surface and that on the tunnel entrance slope face are estimated and compared with the PPV limits in the Chinese standard from GB6722-2014. It is found that when the charge weight varies from 10 to 50 kg and the distance varies from 10 to 100 m, the PPV on the slope face is much lower than the limit, but the PPV on the tunnel surface has the possibility of exceeding the limit. Therefore, the maximum charge weight per delay used in the blasting design of the tunnel blast is determined by the vibration on the tunnel surface. Under the control standard of 15 cm/s at a distance of 10 m, the allowable charge weight per delay should be limited within 25−36 kg.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant nos. 52179102, 51969015, and U1765207) and the Natural Science Foundation of Jiangxi Province (Grant nos. 20204BCJ23002 and 20192ACB21019). The authors wish to express their thanks for all the support they received.