#### Abstract

To study the dynamic response of shallow buried tunnel lining by drilling and blasting method, ABAQUS simulation software was used to establish a tunnel blasting finite element model based on the consideration of in situ stress. Dynamic static coupling numerical simulation was conducted to analyze the vibration response and stress response of tunnel lining. The simulation results were compared and analyzed with field monitoring data to obtain the dynamic response law of tunnel lining. The results show that in the same tunnel lining section, the vibration velocity response in the Y direction is the largest and the vibration velocity response in the Z direction is the smallest. The location of the peak particle velocity of the tunnel lining in three directions appears differently, and the location of the maximum MISES stress appears differently for different excavation sections. The arch shoulder is most affected by horizontal vibration in the X direction, and the vault is most affected by horizontal vibration in the Y direction. The dynamic response at the foot or arch shoulder position away from the lining section of the working face will show the “whip tip effect.” A sudden change in MISES stress occurs at the location of the footing in the 315° direction of the liner section, and the liner is not uniformly stressed in this range.

#### 1. Introduction

In the construction of tunnel projects, most of the rocky tunnels are constructed using the drilling and blasting method. The tunnel project is endowed with an initial stress field environment coupled with self-weight stress and tectonic stress, and shallow tunnels can only consider the initial stress field formed by gravity. The initial displacement field of the strata can generally be regarded as zero [1]. Jelusic et al. [2] used an adaptive network-based fuzzy inference system (ANFIS) to predict and evaluate blast-induced ground vibrations and frequencies. Tian et al. [3] explored the propagation law of blasting vibration velocity in strata and reduced the damage caused by blasting vibration to buildings by carrying out a series of blasting vibration tests. Zhao et al. and [4] Cao et al. [5] investigated the effect of blast-induced vibrations from an adjacent tunnel on the existing tunnel. Field monitoring experiments and a numerical method, i.e., finite element method (FEM), were adopted to study the subject. Yang et al. [6] studied the vibration characteristics on the tunnel surfaces and inside the surrounding rock based on the recorded field data. A three-dimensional dynamic finite-element model was used to verify the site survey results. Kim et al. [7]proposed that the large hole boring method using the state-of-the-art MSP (multisetting smart investigation of the ground and prelarge hole boring) machine (“MSP method”) can efficiently improve vibration reduction. These papers studied the blast vibration problem in tunnel engineering, and the results provide valuable references for further similar pieces of research. Lu et al., Tao et al., Xie et al., and Chen et al. [8–11] studied the effects of in situ stress on the generation and propagation of seismic waves, the mechanism of damage evolution in the surrounding rock, and the expansion pattern of blasting cracks in the surrounding rock. Ma et al., Yilmaz et al., and Omer et al. [12–14] studied the dynamic response, surrounding rock damage, and destabilization damage processes in deep tunnels under blasting load. Relevant studies have shown that in situ stress has a significant impact on rock crack development, surrounding rock damage, blasting stress wave generation, and its energy propagation. Hence, in situ stress has a significant impact on blasting dynamic response and must be taken seriously.

Xu et al. [15] theoretically investigated the dynamic responses of a multilayered half-space subjected to a spatially periodic harmonic moving load using the direct stiffness method. He et al. [16] presented a new semianalytical method to predict the three-dimensional dynamic response of a periodic jointed tunnel in the soil. The simulation is an effective and accurate method to study the vibration effects of blasting in tunnel drilling and blasting construction. However, limited to the own algorithm of the finite element software, it is not possible to calculate in the same analysis step in the in situ stress effect and explosive blast effect. In previous studies, only very few tunnel blasting construction simulations considered in situ stress, and most of the simulation studies were on deep buried tunnels, with few shallow buried tunnels. ABAQUS is a general software for finite element simulation that can complete both static and dynamic analysis. Its operation page is simple and visualized, and the difficulty of making a model is small. In this paper, the ABAQUS simulation software is used to study the blasting vibration response of shallow buried tunnel construction by the drilling and blasting method under the effect of in situ stress. Firstly, the force of the model under in situ stress is calculated, and in situ stress is balanced. Then, the static calculation result of the model is used as the initial condition for blasting analysis. Then, the blasting dynamic calculation is carried out to obtain the blasting vibration response of the model under the in situ stress condition. The simulation results are also compared with the field monitoring results to verify the reliability of the simulation. Then, the tunnel dynamic characteristics are further analyzed to obtain the blasting dynamic response law for shallow buried tunnel construction.

#### 2. Project Overview

Jinjing tunnel is located in Fujian Province, China, and it is a single-line railroad tunnel with a total length of 7292 m and a minimum burial depth of 10 m. The length of Class II surrounding rock area is 1540 m, and the length of Class III surrounding rock area is 3530 m, within which the full section blasting method is used. The length of each cycle is controlled at about 2.5 m. It has been shown that tunnel blasting excavation produces the highest blast loads in cutting holes [17]. In this project, the hole is filled with No.2 emulsified explosive, the diameter of the hole is 42 mm, and the charge of the cutting hole is 24 kg. The cutting hole is arranged in a compound wedge shape. The specific cross-sectional arrangement of the holes is shown in Figure 1. The tunnel dimension plot is shown in Figure 2.

#### 3. Numerical Model and Validation

A shallow buried section of the Jinjing tunnel was selected for numerical modeling, and the model was first calculated by implicit static analysis under the action of self-weight stress and overlying load, and in situ stress equilibrium was performed. Then, using the explicit-implicit coupling method, the static calculation results were used as the initial conditions for the blasting dynamic calculation, and then the explicit dynamic analysis was used to calculate the dynamic response of the model under the blasting action.

##### 3.1. Simulation Model

The simulation model is shown in Figure 3. The model size is 40 m × 36 m × 36 m. The surrounding rock element and the lining element use eight-node linear hexahedral C3D8R element, with 54,994 elements and 75,717 nodes. To avoid the reflection of stress waves caused by the finite boundary of the model, infinite element CIN3D8 was used around the model. The bottom of the model is fixed constrained, the displacement of the corresponding direction is constrained around, and the top of the model is a free surface. Too many holes will make meshing difficult and reduce computational convergence. Hence, set up two holes in the simulation and keep the total charge constant.

##### 3.2. In Situ Stress Load

A gravitational acceleration of 9.8 m·s^{−2} was applied in the vertical direction of the model, and a compressive load of 200 KPa was applied to the top of the shallow buried section because of the stacked load. Then, the in situ stress equilibrium was carried out, so that the geotechnical layer has stress without large deformation, and the results of in situ stress equilibrium are shown in Figure 4.

##### 3.3. Blast Load

In the simulation software, blast load simulation has two methods: explosives modeling and applied blast load profile. Dynamite blasting is completed in a very short time and accompanied by high energy release. In the existing literature [18, 19], the blasting load curve is much longer than the real blasting time, and there is no uniform standard for the peak blasting stress and blasting time. Hence, there are non-negligible disadvantages of using the blasting load curve to simulate the blasting effect. ABAQUS software can use CEL, SPH, and CONWEP methods for explosives modeling to simulate explosives blasting. The most suitable method should be selected according to the actual simulation case characteristics. For example, the CEL method is too expensive, and the CONWEP method is mostly used for TNT point explosion. The SPH method does not have the problem of large mesh deformation, and the software running load is much lower than the CEL method. Hence, the SPH method is used to simulate explosive explosions in this paper. Through the JWL equation of state to describe the relationship among the explosive detonation burst pressure, unit volume energy, and relative volume [20], as shown in equation (1).

In the equation, *P* is the pressure, *V* is the relative volume, which is dimensionless, *E* is the initial unit volume energy of the explosive, *A* and *B* are the explosive material-related parameters, and *R*_{1}, *R*_{2}, and *ω* are the explosive material-related constants, which are dimensionless. The parameters of 2# emulsified explosives are shown in Table 1 [21].

##### 3.4. Material Model and Parameters

We study the vibration of the lining and do not involve the damage study of the rock, and the rock material model is selected from the Mohr Coulomb model. According to the engineering geological data, the physical and mechanical parameters of the rock material are shown in Table 2.

The dynamic strength and static strength of tunnel lining differ greatly, and Wang [22] proposed a formula for calculating the static elastic modulus and dynamic elastic modulus of the perimeter pressure lining structure. The physical and mechanical parameters of the concrete lining are shown in Table 3.

In the equation, is the static modulus of elasticity, and is the dynamic modulus of elasticity.

##### 3.5. Simulation Results and Validation

Along the tunnel axis in the arch waist at the layout of five measurement points, measurement point 1 # is from the tunnel surface of 15 m, the measurement point spacing of 5 m, and the specific layout of measurement points is shown in Figure 5.

The peak particle velocity in the model is extracted from the points corresponding to the field measurement, and the field measurement data and simulation data are compared and analyzed, as shown in Tables 4 and 5.

The analysis of data shows that the simulation results with in situ stress are within 15% error, with the maximum error of 13.87% from the measured data in the field. The simulation results without in situ stress are within 30% error with the field measured data, with a maximum error of 27.58%. It can be seen that the numerical simulation results with in situ stress are more consistent with the field monitoring results. Hence, it is reasonable to consider in situ stress in the numerical simulation of blast vibration response in a shallow buried tunnel. As can be seen from Table 4, the peak particle velocity in the Y direction is greater than the peak particle velocity in the X and Z directions. Hence, the blast vibration velocity in the Y direction plays a controlling role in the blast vibration effect [23], and the vibration velocity in the Y direction is used to evaluate the blast vibration strength. The time course diagram of the simulation of vibration velocity in the Y direction at measurement point 3# is extracted, and the measured time course diagram is shown in Figures 6 and 7. Comparing the time course of vibration velocity plot, it can be seen that the simulated result of sustained higher vibration velocity is 36 ms and starts to decay to a very small value by 94 ms. The measured result is 21 ms for the sustained higher vibration velocity, which starts to decay to a very small value at 106 ms. The decay law of the simulated vibration velocity waveform with in situ stress is close to that of the measured vibration velocity waveform in the field. It further indicates that it is reasonable to consider in situ stress in the simulation of tunnel blasting vibration response.

#### 4. Analysis of Numerical Simulation Results

The peak particle velocity and peak MISES stress in the X, Y, and Z directions of the nodes of the section where the five measurement points are located are extracted, plotted in polar coordinates, and analyzed to obtain the distribution patterns of peak particle velocity and peak MISES stress in different sections of the tunnel. Because of the large amount of data from the simulation results, we wrote scripts to extract the data.

##### 4.1. Distribution Law of Lining Vibration Velocity

The nodal schematic of the tunnel lining section is shown in Figure 8, and the nodal vibration velocity peak polar coordinates of the measurement points are shown in Figure 9.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(k)**

**(l)**

**(m)**

**(n)**

**(o)**

The vibration velocity response of the lining presents a complex situation in tunnel blasting construction under the effect of in situ stress. An analysis of Figure 7 was performed. In the same tunnel lining section, the peak particle velocity in the Y direction is greater than the peak particle velocity in the X direction and the peak particle velocity in the Z direction. There are differences in the peak distribution of vibration velocity in the three directions at different locations of the same lining section. The peak particle velocity in the Z direction has the smallest value, and its distribution is the most uniform compared to the peak particle velocity in the X and Y directions. The absolute difference among the peak vibration velocities at different locations in the section is the smallest, which indicates that the blasting load has the smallest vibration effect on the longitudinal direction of the tunnel. Peak particle velocity in the X direction is not uniformly distributed along the lining section, and the maximum velocity appears in the arch waist and arch shoulder. The closer the lining section is to the tunnel face, the wider the distribution of the peak particle velocity in the section with a larger value. As the distance from the tunnel surface gradually increases, a large number of values of the X-vibration velocity peak distribution shifts to the arch shoulder, arch waist, and foot. The X-vibration velocity peak polar graph shows an asymmetric shape, with large abrupt changes in the four directions of 60°, 120°, 210°, and 330°, and the peak particle velocity curve becomes unsmooth. The shape of the tunnel section is obtained by the intersection of 3 circles, and these 4 directions coincide with the direction of the intersection of the 3 circles, which shows that the shape of the tunnel section has a large effect on the X-vibration velocity. As the distance from the section to the tunnel surface increases, the maximum value of X-vibration velocity does not decrease but increases in the arch shoulder and foot position of the lining section, and the phenomenon of “whip tip effect” appears. The peak Y-directional vibrational velocity is also unevenly distributed along the lining section, and the maximum values of vibrational velocity appear in the top, bottom, and foot of the arch. The Y-vibration velocity peak values are basically distributed in the top, bottom and foot of the arch regardless of the distance from the tunnel face, and the Y-vibration velocity peak polar diagram also shows an asymmetric shape. As the distance from the section to the tunnel face increases, the maximum values of Y-vibration velocity peaks do not differ greatly, however, they appear at different locations in the section. In the same position of the foot of the lining section, the peak Y-vibration velocity increases, and the phenomenon of “whip tip effect” appears.

For the reason that the vibration velocity peak polar plot shows asymmetry, it can be explained that although the model and constraints are symmetrical along the longitudinal axis of the tunnel, the tunnel section is not standard circular, the stress wave generated by the explosive explosion of the hollowing hole is not propagated to the tunnel section with the same intensity at the same time, and the blasting load on the tunnel is asymmetrical. Hence, the vibration velocity peak polar plot is asymmetrical.

For the “whip tip effect” of the peak horizontal vibration velocity in the X direction and horizontal vibration velocity in the Y direction, the vibration velocity time curve of the 2 nodes at the arch waist, which are far and near the distance to the tunnel face, is extracted for comparison and analysis. See Figures 10–13.

As shown in Figures 10–12 and 13, although the peak vibration velocity of the measurement point far from the tunnel face is larger than the peak vibration velocity of the near point, the decay speed of the vibration velocity of the far point is much faster than that of the near point, and the X- and Y-vibration velocity of the near point decays to a small value at 80 ms and 77 ms. The X- and Y-vibration velocity of the far-point decays to a small value at 30 ms and 29 ms. The reason for this phenomenon is that the farther the liner is from the palm surface, the lesser the energy generated by the blasting that spreads to the liner, and the liner vibration is not maintained by more energy. Hence, the decay will be accelerated naturally. We believe that the “whip tip effect” is related to the modal analysis of the lining. After propagating for a certain distance, the frequency of the wave generated by blasting may be close to the natural frequency of a certain modal of the lining, thus causing the lining vibration to become violent.

##### 4.2. Lining MISES Stress Distribution Law

The polar coordinates of the peak MISES stress at the section nodes are shown in Figure 14.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

The peak MISES stress distribution for tunnel blasting construction under in situ stress is different for different tunnel lining sections, each with its characteristics but also with the same points. Figure 13 is now analyzed. In the same tunnel lining section, the MISES stress peak polar plot shows an asymmetrical feature, and its maximum magnitude value appears in the arch shoulder. In the direction of 315° at the foot of the lining, the MISES stress peak curves undergo a large abrupt change, and the curves become unsmooth. It indicates that around the range of 315° in the lining section, the MISES stress value changes faster, its discontinuity increases, and the lining structure is not uniformly stressed in this small range. With the increase of the distance of the lining section from the tunnel face, the maximum value of the peak MISES stress will increase, however, the increase is limited, 1.46 MPa and 2.79 MPa, respectively. The reason for the increase of the maximum value of the peak MISES stress is closely related to the “whip tip effect” of the X- and Y-vibration velocity. The greater the vibration effect of the liner section, the greater the MISES stress. This phenomenon shows that it is not the case that the farther the lining section is from the working face, the smaller the force will be under the blasting action. As the distance from the working face increases, the smoothness of the MISES stress peak curve at the arch shoulder and arch waist increases, and the curve is fuller here, indicating that the stress at the arch shoulder and arch waist increases and is uniform as it moves away from the tunnel face.

#### 5. Conclusion and Discussion

(1)In the construction of tunnel engineering drilling and blasting method, the tunnel vibration response under blasting is not symmetrical along the longitudinal axis of the tunnel. The Y-vibration velocity response of the liner is the largest, Z vibration velocity response is the smallest, the arch shoulder is most affected by horizontal vibration in the X direction, and the vault is most affected by horizontal vibration in the Y direction. The shape of the tunnel liner section has a greater effect on the vibration response in the X direction.(2)At the foot or shoulder of the tunnel lining section, the peak X- and Y-vibration velocity of the lining far from the tunnel surface will be greater than the peak X- and Y-vibration velocity of the lining near the tunnel surface. The phenomenon of “whip tip effect” appears, however, the decay speed of the former will be significantly faster than that of the latter.(3)The maximum value of MISES stress in the tunnel lining occurs at the shoulder of the arch, and a sudden change in the peak MISES stress occurs in the 315° direction of the lining section, which should be noted in the design of the lining reinforcement within this range. Away from the tunnel face, the maximum value of MISES stress in the lining section will increase, and the MISES stress curve at the arch shoulder and arch waist will be smoother and fuller. The lining will be more uniformly stressed here.(4)Because of the limited computer capability, it is not possible to simulate the tunnel blasting construction of real projects in full model; however, the research method adopted in this paper has the advantages of simplicity and efficiency. The “whip tip effect” in this numerical simulation is only obtained within 40 m from the tunnel face, and the detailed explanation of this phenomenon is not explored in depth in this paper. We believe that the “whip tip effect” is related to the modal analysis of the lining, and the follow-up research direction is the relationship between blast vibration and modal analysis of the lining.#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was funded by the Natural Science Foundation of China (No. 52168055), the Natural Science Foundation of Jiangxi Province (20212ACB204001), and “Double Thousand Plan” Innovation Leading Talent Project of Jiangxi Province (jxsq2020101001).