Dynamic Response Analysis of Underground Double-Line Tunnel under Surface Blasting

Blasting has been widely used for economical and rapid rock excavation in civil and mining engineering. In order to study the influence of relative horizontal distance and relative vertical distance between two tunnels on the dynamical response of the two tunnels, 10 numerical simulation cases are done by LS-DYNA 3D models under surface explosion by controlling the clear distance and height difference of double-line tunnel, and the ALE multimaterial fluid structure coupling algorithm is applied to analyze the dynamic response characteristics of double-line tunnel under different conditions. The numerical results show that the dynamic response characteristics of the tunnel lining are affected by the change of the clear distance and height difference of the tunnel. With the increase of the height difference between adjacent tunnels, the peak value of vibration velocity at the top of the lining on the blast face increases, which is due to the upward elevation of the right tunnel, which is more conducive to the reflection and superposition of stress waves. When the height difference of tunnel is 4–6 m, the vibration velocity and displacement of monitoring point C on the back blasting side will change abruptly, and the variation range of vibration velocity is about 25%, while the variation range of displacement is about 60%.


Introduction
Blasting is widely used in rock excavation, mining, and other hard rock engineering applications since it has the advantages of strong geological adaptability and economy. However, the blast-induced vibration may do harm to the surrounding buildings, tunnels, and other structures. At present, many studies have been done for blasting-induced damage or effects on surface buildings, while the research studies on underground structures are quite rare relatively. With the increase of underground space building, the underground tunnel engineering are booming and face lots of blasting problems, such as tunnel rock fall. erefore, it is very important to study how to effectively protect underground structures and minimize the loss under the strong explosion.
Many scholars have done research on the dynamic response characteristics of underground structures under explosive loading. Eitzenberger [1] observed through experimental and numerical investigations that the attenuation of the shock wave is controlled by the texture of the rock mass. Wu et al. [2] investigated the propagation of blastinduced shock wave in jointed rock mass using accelerometers and found that the attenuation of the shock wave is completely dependent on the distance from the charge weight, the angle between the strike of rock strata, and the direction of wave propagation. Yang et al. [3] investigated the response of blast-induced vibration on tunnel surfaces and inside surrounding rock using three-dimensional (3D) numerical analysis procedure. Numerical studies 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. Liang et al. [4] investigated the response of an existing tunnel subjected to blast-induced vibrations from a newly constructed tunnel placed adjacent to the existing tunnel. Mohammad and Rehan Sadique [5] considered an internal blast loading on a rock tunnel constructed in quartzite rock; the results show that the extent of damage in shallow depth tunnels is found to be more than that of the tunnels at higher depth of overburden. Liu et al. [6] investigated the explosion inside tunnel and a formula for the explosion blast wave overpressure at a certain distance from the detonation center point inside the tunnel was derived by using the dimensional analysis theory. Feldgun et al. [7] studied internal blast loading in a buried lined tunnel by the modified Godunov method, which considered all the stages of the process: detonation of the internal charge; the shock wave propagation in the internal gas and its following interaction with the cavity's shell lining including multiple reflections. Yang et al. [8] performed numerical modeling to assess the damage characteristics of an underwater tunnel subjected to blast loads and explore the potential mitigation measures based on coupled Lagrange and Euler (CLE) method.
e results show that the rigidity and load-carrying ability of the tunnel are significantly improved by bonding the CFRP cloth. e recommended thickness of the CFRP cloth is 0.5-0.835 mm. Koneshwaran et al. [9] investigated underground transport tunnels in blast loading; the results indicated that several bolts failed in the longitudinal direction due to redistribution of blast loading to adjacent tunnel rings, and the tunnel segments respond as arch mechanisms in the transverse direction and suffered damage mainly due to high bending stresses. e above researches are of great significance to understand the mechanism and process of underground space dynamic disaster. However, the researches on the dynamic response characteristics of underground double-line tunnel under the surface blasting are relatively rare. e research on the structures under explosive loading mainly adopts the methods of experiment, theoretical analysis, and numerical simulation. e explosion experiment is the most effective and direct method to study the dynamic response characteristics of the structures, but the destructive experimental conditions are harsh and costly. Under the blast impact load, the underground structure is not only affected by the stress waves from all directions, but also influenced by the nonperiodic transient action and considering the plastic strain of the material, which makes the problem become a highly nonlinear problem combining state nonlinearity and material nonlinearity [10][11][12][13][14]. erefore, it is difficult to realize the blast mechanics analysis of complex structures by theoretical means. For these reasons, the use of numerical method to study the complex progress between wave propagation and tunnel interaction is necessary and effective. In this paper, in order to study the influence of relative horizontal distance and relative vertical distance between two tunnels on the dynamical response of the two tunnels, 10 numerical simulation cases are done by LS-DYNA 3D models under surface explosion. e research results can provide a useful reference for the antiexplosion protection design of underground tunnel.

Finite Element Modeling of Rock and Lining.
In order to study the influence of relative horizontal distance and relative vertical distance between two tunnels on the dynamical response of the two adjacent tunnels, 10 numerical simulation cases are done by LS-DYNA 3D models under surface explosion. e model is constructed in cm-g-μs unit system. Figure 1 shows the geometric model of the double-line tunnel. Figures 2 and 3 are the isometric side view and front view of the finite model of the double-line tunnel, respectively. e size of whole model is 80 m × 10 m × 42 m. e outer diameter of the lining is 10 m, the wall thickness of the lining is 0.3 m, and the longitudinal length of the lining is 10 m. e charge of the TNT explosive is 40 kg. e distance between the side and top of the left lining and the center of TNT explosive is 30 m and 20 m, respectively, as shown in Figure 1. ree degrees of freedom (UX, UY, and UZ) of the bottom of the finite element model are constrained. In order to simulate the infinite region and eliminate the influence of the reflected stress wave on the simulation results, except for the top free surface, the other boundary segment of the model is controlled by keyword * Boundary_Non_Reflecting, which can absorb the expansion wave and shear wave passing through the interface. Considering the influence of lining clear distance and height difference on dynamic response characteristics of double line tunnel under surface blasting, the following parametric analysis scheme is formulated. When analyzing the influence of lining clear distance variable on dynamic characteristics, ensure that the height difference of two linings is 0 m, the position of left lining relative to TNT explosive remains unchanged at the same time, only change the horizontal clear distance of right lining, and the horizontal clear distances of lining are 5 m, 7 m, 9 m, 11 m, and 13 m, respectively. In order to simplify the subsequent analysis, the corresponding conditions are marked as A1, A2, A3, A4, and A5, respectively, as shown in Figure 4 and Table 1. When analyzing the influence of lining height difference variables on dynamic characteristics, the horizontal clear distance between two linings is 9 m, and the position of left lining relative to TNT explosive remains unchanged. e height differences between two linings are 0 m, 2 m, 4 m, 6 m, and 8 m, respectively. In order to simplify the subsequent analysis, the corresponding conditions are marked as A3, B3, C3, D3, and E3, respectively, as shown in Figure 4 and Table 1.

ALE Multimaterial Fluid Structure Coupling Algorithm.
e Arbitrary Lagrange-Euler (ALE) coupling algorithm is used in LS-DYNA 3D to solve the fluid structure coupling problem in this paper, which is the same as the description of Euler algorithm. It can be understood that there are two layers of grids overlapped together, but the difference is that the grid in ALE algorithm is not fixed and can move arbitrarily in space. e ALE algorithm first performs several Lagrange time step calculations, in which the element mesh deforms with the material flow and then performs ALE time step calculation. (1) e boundary conditions of the modified object remain unchanged, and the internal element is meshed to keep the topological relationship of the mesh unchanged. is step is called smooth step. (2) e element parameters (density, energy, stress tensor, etc.) and node velocity vector in the deformed mesh are transferred to the new mesh, which is called expectation step. ALE algorithm on the one hand retains the advantages of Lagrange algorithm [10,15,16]; that is, it can accurately detect the boundary of the mesh, but also inherits the main advantages of Euler algorithm, which can well solve the problem of element distortion, make up for the shortcomings of the two algorithms, and is very suitable for large deformation analysis.
In this simulation, the rock mass and lining are described by the Lagrange method; the explosive and air are described by the Euler method. is method couples fluid and solid together to realize the nonlinear coupling between fluid medium and rock mass model by keyword * Constrained_Lagrange_In_Solid. e fluid structure coupling algorithm is to couple the structure and fluid together through certain constraint method to realize the transmission of mechanical parameters. e main constraint methods [17] are velocity constraint, acceleration constraint, and penalty function constraint. e advantage of this algorithm is that the fluid element and structural element on the coupling surface do not need to be corresponded one by one, which greatly reduces the workload of mesh generation. e calculation steps of velocity and acceleration constraints are as follows: e fluid element with structural nodes is searched, and the structural node parameters (mass, momentum, and nodal force) are assigned to the fluid element nodes. (1) Calculate the new acceleration of fluid node (velocity): e acceleration (velocity) of constrained structure node is as follows: where m n and m o represent the nodal mass of fluid element before and after distribution, respectively; M and F are the momentum and nodal force, respectively; a and v are the acceleration and velocity of the node; h is the number of nodes contained in a single fluid element; and f and s are symbols of fluid and solid element.

Constitutive Model of Air.
e air model [19] is described by keyword * MAT_NULL combined with multilinear equation of state * EOS_LINEAR_POLYNOMIAL in this study.
For the convenience of calculation, air is regarded as an ideal gas, in which the parameters are as follows: For air in cm-g-μs unit system, the parameters of the null material model are as follows: density ρ � 0.0012 g/cm 3 and dynamic viscosity coefficient Mu � 0.001.

Constitutive Model of Concrete.
e concrete material adopts the model which adapts to high pressure and high strain rate by means of keyword * MAT_JOHNSON_HOLMQUIST_CONCRETE (HJC). e material model of HJC consists of state equation, yield equation, and damage equation. e equation of state can be divided into elastic stage, plastic stage (internal porosity compression, porosity reduction), and fully dense stage (internal porosity compression, damage, and fine grain cracks). e characteristic of the HJC model is that it can reflect the dynamic response of brittle materials such as concrete under large strain, high strain rate and high pressure, and material damage effect. It is especially suitable for the study of dynamic response of concrete structure under explosion load. e equations of state of HJC material model are as follows.
Elastic loading and unloading (p < p c ) is given by Loading of plastic transition zone (p c ≤ p ≤ p 1 ) is given by Unloading of plastic transition zone (p c ≤ p ≤ p 1 ) is given by Fully compacted loading (p > p 1 ) is given by In the fully compacted unloading (p > p 1 )stage, the material is completely destroyed: e yield equation of HJC material model is as follows: e damage equation of HJC material model is as follows: where k e � p c /μ c is the bulk modulus, which is the ratio of crushing volume pressure p c and crushing volume μ c strain in uniaxial compression test, p 1 represents compaction pressure of concrete material, and μ 1 is the compaction volume strain. K 1 is the plastic bulk modulus; p max and μ max are the maximum volume pressure and volume strain before unloading. In this stage, the porosity in the material is compressed, the material is damaged, and the crushing crack begins to appear; σ * is the standardized equivalent stress, D is the damage value (0 ≤ D ≤ 1.0), p * is the standardized hydrostatic pressure, ε * � (ε/ε 0 ) is the dimensionless strain rate, A is the normalized cohesive strength; B is the normalized pressure hardening coefficient; N is the pressure where p is hydrostatic stress, q represents von Mises stress, and K is a scalar parameter that determines the shape of the yield surface and maintains the convexity of the yield surface in the deviatoric (p) plane. r is the third invariant of the deviatoric stress tensor: where J 2 and J 3 represent the second and the third deviatoric stress invariants, respectively.
where φ is the angle of friction and C is the dilation angle. e material properties used for quartzite rock mass are given in Table 3, obtained from triaxial test [21]. Figure 5 shows the pressure nephogram of rock mass under A1 condition. At the moment of TNT explosive explosion, the explosive volume expands rapidly in a very short time and rapidly changes from solid state to high-pressure gas state. e high-pressure gas acts on the rock and produces stress wave (t � 0.5 ms) on the rock mass, with the peak pressure of about 228 MPa. Due to the nonreflective boundary condition, the shock wave will not be reflected at the rock mass boundary. With the further propagation and diffusion of the shock wave in the rock mass, the amplitude of the stress wave is greatly weakened. When t � 6.3 ms, the peak value of the pressure decreases to about 4.32 MPa. Combined with Figure 5 (t � 10 ms, t � 17 ms, and the velocity time history curve of lining monitoring points in Figure 6), at t � 10 ms, the stress wave propagates to the left lining, and at t � 17 ms, the stress wave propagates to the right lining. With the further attenuation of the shock wave, the stress wave gradually becomes elastic wave, which propagates at the elastic wave speed which does not disturb the physical state of the rock mass. It can also be seen from Figure 7 that the peak velocity of monitoring point A is about 7.5 cm/s at t � 27 ms, and the time from stress wave propagation to monitoring point A to reaching the peak velocity of monitoring point A is about 17 ms. e peak velocity of point B is about 6.4 cm/s at t � 46 ms, and the time from stress wave propagation to point B to the peak velocity is about 29 ms. It can be seen that the time course of reaching the peak speed is different, and the time difference is about 12 ms. is is because the stress wave produces a series of reflection and diffraction between two tunnels, so its dynamic response process is different from that of a single tunnel. In view of this, this paper starts the follow-up research, that is, by changing the horizontal spacing and height difference of the tunnel to investigate the dynamic response characteristics of the double-line tunnel under the surface explosion load. Figures 8 and 9 show the equivalent stress nephogram of concrete lining under different conditions (considering horizontal clear distance and height difference) at the same time (t � 20 ms), respectively. For the tunnel on the left, the stress distribution on the left side of YOZ symmetry plane of the lining is basically the same. For the right side of the YOZ symmetry plane of the lining, the stress distribution shows obvious differences. According to the color corresponding to the cloud figure, the darker the color of the element is, the smaller the equivalent stress is. It can be seen from Figure 8 that with the increase of the clear distance of the lining (5 m, 7 m, 9 m, 11 m, and 13 m), the stress on the right side of the symmetry plane of YOZ of the lining is not regularly decreasing. It can be seen that under A3 and A4 conditions, the overall value of the equivalent stress in this area is obviously smaller than the other three conditions (A1, A2, and A5).

Distribution Law of Equivalent Stress of Lining.
is shows that the influence of the tunnel spacing on the dynamic response characteristics of the lining is not linear. Similarly, it can be seen in Figure 9 that the dynamic response characteristics of the lining also show a nonlinear relationship with the variation of the lining height difference (0 m, 2 m, 4 m, 6 m, and 8 m). Figure 10 is the layout of monitoring points for the left lining. Monitoring points A, B, C, and D are all located on the left lining. Figure 11 shows the relationship curve between the velocity of the monitoring point and the horizontal clear distance of the tunnel. It can be seen from Figure 11 that the velocity response of the monitoring point at the top of the lining blast facing surface is lower than that of the side wall under the action of explosion load. Taking the clear distance of the lining as an example, the peak vibration velocity of the Shock and Vibration monitoring point B at the top of the lining blast facing surface is about 7.5 cm/s, and the peak vibration velocity of the monitoring point A at the side wall of the lining blast facing surface is about 6.5 cm/s. e vibration velocity is about 13% lower than that of B. At the same time, it can be seen that the peak velocity of monitoring points A and B almost does not change with the increase of the horizontal spacing of the lining, which indicates that the influence of the change of the tunnel spacing on the vibration velocity of the top of the lining face and the side wall can be ignored.

1.000e -06
Effective stress (v -m) 9.500e -07 9.000e -07 8.500e -07 8.000e -07 7.500e -07 7.000e -07 6.500e -07 6.000e -07 5.500e -07 5.000e -07 4.500e -07 4.000e -07 3.500e -07 3.000e -07 2.500e -07 2.000e -07 1.500e -07 1.000e -07 5.000e -08 0.000e + 00     0.27 MPa, and the peak value of equivalent stress decreases by about 10%. For point D, the peak value of equivalent stress decreases from 0.28 MPa to 0.24 MPa, and the peak value of equivalent stress decreases by about 14%. As for Figures 13 and 14, it can be seen from the analysis that no matter for the X-direction stress or Y-direction stress of the monitoring point, the peak stress of the monitoring point A has almost no obvious change. is is because although the refraction and diffraction effect of the rock stress wave between the two linings is produced, whose influence is mainly on the lining on the back blasting side, the influence on the monitoring point of the lining side wall on the blast facing side can be ignored. For point C on the blasted side, it can be seen from the analysis of Figures 13 and 14 that the stress components of X and Y are in a downward trend on the whole, while for the monitoring point D on the bottom of the lining, it can be seen that the stress component in the Y-direction has a slight increment with the increase of the net distance of the tunnel, and the stress peak value in the Ydirection increases from 0.15 MPa to 0.18 MPa, an increase of about 20%. Figure 15 shows the relation between the velocity of the monitoring point and the height difference of the tunnel. It can be seen from Figure 15 that under the action of surface blasting, similar to the previous analysis, the peak value of vibration velocity at monitoring point A on the blast facing surface of lining side wall does not change significantly. For the monitoring point B on the blasting face of lining, with the increase of tunnel height difference, the peak value of vibration velocity increases. is is because as the right tunnel rises upward, which is more conducive to the reflection and superposition of stress wave, so the vibration velocity of measurement point B increases. It can be seen that the vibration velocity of monitoring point B increases from 7.3 cm/s to 7.7 cm/s, with an increase of about 5.5%. When the height difference of lining is between 0 and 4 m, the sensitivity of vibration velocity of side wall monitoring point C is low, and there is almost no significant change in vibration velocity. However, when the height difference is between 4 and 6 m, the vibration velocity changes from 5 cm/s to 5.5 cm/s, and the increase of vibration velocity in this interval is about 10%. For point D at the bottom of the lining, the vibration velocity changes significantly with the increase of tunnel height difference. It can be seen that the vibration velocity of point D increases from 5.7 cm/s to 7.1 cm/s, and the peak velocity increases by about 25%. Figures 16-18 are the curves of the relationship between the equivalent stress of the monitoring point and the height difference of the tunnel, the relationship between the Xdirection stress of the monitoring point and the height difference of the tunnel, and the relationship between the Ydirection stress of the monitoring point and the height difference of the tunnel, respectively. It can be seen from Figure 16 that, with the increase of the tunnel height difference, the equivalent stress of points A and B on the blasting face remains stable, and the peak values of equivalent stress are about 0.36 MPa and 0.43 MPa, respectively. For the back burst surface lining test point C, the peak value of equivalent stress increases slightly with the increase of tunnel height difference, and the peak value of equivalent stress increases from 0.24 MPa to 0.26 MPa, with an increase of about 8.3%. For the back blasting face lining monitoring point D, with the increase of tunnel height difference, the peak value of equivalent stress of point D presents a small downward trend, and the peak value of equivalent stress decreases from 0.29 MPa to 0.27 MPa, with a decrease of about 6.9%. As for the X-direction stress component of the monitoring points, it can be seen from the analysis of Figure 17 that except for point A on the blasting face, the Xdirection stress component of the other monitoring points increases first when the height difference of the lining is 0∼4 m and then decreases when the height difference is about 4∼8 m. It shows that the reflected shock wave has the greatest influence on the X-direction stress component of the monitoring point on the back blasting surface of the left tunnel when the height difference of the adjacent tunnel is 4 m and the horizontal clear distance is constant. It can be seen from Figure 18 that the peak value of the Y-direction stress at the monitoring point B on the blasting face of the lining decreases slowly with the increase of the tunnel height difference, while for the monitoring point C on the blasting face of the lining, the component of the Y-direction stress decreases slightly with the increase of the tunnel height difference. Figure 19 is the layout of the monitoring point path of the back burst surface of the left lining. Figures 20 and 21 show the relationship between the displacement along the path (along the negative direction of Z) of the horizontal net distance and the tunnel height difference (t � 31 ms), respectively. It can be seen from the figure that the displacement of the lining monitoring point is decreasing along the path direction as a whole. It can also be seen that the displacement of the monitoring points of condition 1, A1 and A5 is larger than A2, A3, and A4, which shows that the variation of displacement of the monitoring points is not very significant due to the lining spacing. For condition 2, it can be seen from Figure 21 that with the increase of the tunnel height difference, the displacement of the monitoring point increases with the increase when the path length is less than 4 m. In addition, the displacement increases obviously in the interval with the tunnel height difference of 4-6 m.   Taking the initial point of the path (distance � 0) as an example, the displacement of the monitoring point under C3 condition (height difference � 4 m) is about 0.0025 cm. e displacement of the monitoring point is about 0.004 cm under D3 condition (height difference � 6 m), which is about 60% higher than that of the other. Figure 22 is the path layout of rock monitoring points. e monitoring points are arranged clockwise along the zaxis. It can be seen that the variation law of the vibration velocity along the path of the monitoring points the same no matter in condition 1 or condition 2. e vibration velocity of the tunnel top reaches the peak at about 100 cm along the path, rather than just above the top of the tunnel, as shown in Figures 23 and 24. It can also be seen that the lowest vibration velocity is at the bottom of the tunnel. It can be seen from the dispersion degree of the vibration velocity value of two curves that the dispersion degree of condition 1 will be higher, which indicates that the change of the horizontal spacing of the tunnel will have more significant influence on the vibration velocity of each monitoring point than the change of the tunnel height difference.

Conclusions
Based on LS-DYNA 3D nonlinear finite element software, a full coupling model of TNT explosive-surrounding rocklining structure-air is established. Lagrange algorithm is used for lining and rock mass materials. Euler algorithm is used for air and explosive materials. e nonlinear coupling between Euler fluid domain and Lagrange structure domain is realized by ALE multimaterial fluid structure coupling algorithm. By controlling the clear distance and height difference of double-line tunnel, a variety of numerical simulation cases are formulated. e dynamic response characteristics of double-line tunnel under different conditions and surface explosion load are systematically studied. e main conclusions are as follows: Under the action of surface blasting, the influence of the tunnel clear distance and height difference on the dynamic response characteristics of the lining presents a nonlinear variation law except for the monitoring point A on the blasting face. e peak velocity of monitoring points A and B on the blast facing surface hardly changes with the increase of the horizontal spacing of the lining, indicating that the influence of the change of the tunnel spacing on the vibration velocity of the top and side walls of the blast facing surface of the lining can be ignored. For point C on the side wall of back blasting face, the peak velocity of point C decreases about 26.5% with the increase of horizontal clear distance of lining. For point D at the bottom of back blasting face, the peak velocity of point D decreases about 20.9% with the increase of horizontal clear distance of lining. With the increase of the horizontal distance of the tunnel, the equivalent stress of the measuring points decreases on the whole.
With the increase of the height difference between adjacent tunnels, the peak value of vibration velocity of lining monitoring point B increases, which is more conducive to  the reflection and superposition of stress waves due to the upward elevation of the right tunnel. When the height difference between adjacent tunnels is 4 m, the reflected shock wave has the greatest influence on the X-direction stress component of the measuring points on the back blasting surface of the left tunnel. When the tunnel height difference is about 4∼6 meters, the vibration velocity and displacement of monitoring point C on the back blasting side will change abruptly, and the variation range of vibration velocity is about 25%, while the variation range of displacement is about 60%.
Data Availability e data are available and explained in this article; readers can access the data supporting the conclusions of this study. Disclosure e authors would like to declare that the work described herein is original research and has not been previously published elsewhere.

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding the publication of this paper.