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

Authors: Pei Zhang; Jianhua Cai; Feng Zong; Yanpeng He; Qiong Wang
Journal: Shock and Vibration
 (2021)

---

## Abstract

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


---

## Body

## 1. 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. Therefore, 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 blast-induced 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. The results show that the rigidity and load-carrying ability of the tunnel are significantly improved by bonding the CFRP cloth. The 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.The 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. The research on the structures under explosive loading mainly adopts the methods of experiment, theoretical analysis, and numerical simulation. The 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–14]. Therefore, 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. The research results can provide a useful reference for the antiexplosion protection design of underground tunnel.


## 2. Three-Dimensional Finite Element Modeling

### 2.1. 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. The 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. The size of whole model is 80 m × 10 m × 42 m. The 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. The charge of the TNT explosive is 40 kg. The 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. Three 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.Figure 1
Isometric side view of geometric model.Figure 2
Equiaxed side view of finite element model.Figure 3
Front view of finite element model.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 Figure4 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. The 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.Figure 4
Details of numerical simulation cases.Table 1
Details of numerical simulation cases.
Numerical simulation case numberHorizontal clearance of lining (m)Numerical simulation case numberHeight difference of lining (m)A15A30A27B32A39C34A411D36A513E38


### 2.2. ALE Multimaterial Fluid Structure Coupling Algorithm

The 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. The 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) The boundary conditions of the modified object remain unchanged, and the internal element is meshed to keep the topological relationship of the mesh unchanged. This step is called smooth step. (2) The 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. This method couples fluid and solid together to realize the nonlinear coupling between fluid medium and rock mass model by keyword∗Constrained_Lagrange_In_Solid. The fluid structure coupling algorithm is to couple the structure and fluid together through certain constraint method to realize the transmission of mechanical parameters. The main constraint methods [17] are velocity constraint, acceleration constraint, and penalty function constraint. The 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. The calculation steps of velocity and acceleration constraints are as follows:The 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)mnM,Ff,i=m0M,Ff,i+himM,fs.Calculate the new acceleration of fluid node (velocity):(2)avf,i=FMf,imfn,i.The acceleration (velocity) of constrained structure node is as follows:(3)avs=∑1ihiavf,i,where mn and mo 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.


### 2.3. Constitutive Model of Materials

#### 2.3.1. Constitutive Model of TNT Explosive

Explosive material model is described by keyword∗MAT_HIGH_EXPLOSIVE_BURN combined with JWL equation of state provided by LS-DYNA. The P-V relation of the JWL equation of state is as follows:(4)P=A1−ωR1Ve−R1V+B1−ωR2Ve−R2V+ωEV,where A, B, R1, R2, and ω are material constants for TNT explosive, parameters A and B represent the magnitudes of pressure, ρ is the density of the explosive, and E is the specific internal energy at atmospheric pressure.For TNT  explosive [18], ρ = 1630 kg/m3; detonation wave speed D = 6930 m/s; Chapman–Jouget pressure PCJ=0.255×1011Pa; A = 373800 MPa; B = 3747 MPa; R1 = 4.15; R2 = 0.9; and ω=0.35.


#### 2.3.2. Constitutive Model of Air

The air model [19] is described by keyword ∗MAT_NULL combined with multilinear equation of state ∗EOS_LINEAR_POLYNOMIAL in this study.(5)P=C0+C1μ+C2μ2+C3μ2+C4+C5μ+C6μ2e.For the convenience of calculation, air is regarded as an ideal gas, in which the parameters are as follows:C0 = C1 = C2 = C3 = C6 = 0, C4 = C5 = 0.4. For air in cm-g-μs unit system, the parameters of the null material model are as follows: density ρ=0.0012g/cm3 and dynamic viscosity coefficient Mu = 0.001.


#### 2.3.3. Constitutive Model of Concrete

The concrete material adopts the model which adapts to high pressure and high strain rate by means of keyword∗MAT_JOHNSON_HOLMQUIST_CONCRETE (HJC). The material model of HJC consists of state equation, yield equation, and damage equation. The 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). The 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. The equations of state of HJC material model are as follows.Elastic loading and unloading (p<pc) is given by(6)p=keμ.Loading of plastic transition zone (pc≤p≤p1) is given by(7)p=p1+p1−pcμ−μ1μ−μ1.Unloading of plastic transition zone (pc≤p≤p1) is given by(8)p−pmax=1−μmax−μcμ1−μcKe+μmax−μcμ1−μcK1μ−μmax.Fully compacted loading (p>p1) is given by (9)p=k1μ−μ11+μ1+k2μ−μ11+μ12+k3μ−μ11+μ13.In the fully compacted unloading (p>p1)stage, the material is completely destroyed:(10)p−pmax=k1μ−μ11+μ1−μ−μ11+μ1max.The yield equation of HJC material model is as follows:(11)σ∗=A1−D+BP∗N1+Clnε∗.The damage equation of HJC material model is as follows:(12)D=∑Δεp+ΔμpD1P∗+T∗D2,where ke=pc/μc is the bulk modulus, which is the ratio of crushing volume pressure pc and crushing volume μc strain in uniaxial compression test, p1 represents compaction pressure of concrete material, and μ1 is the compaction volume strain. K1 is the plastic bulk modulus; pmax 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 hardening index; C represents the strain rate coefficient; Δεp and Δμp are the equivalent plastic strain increment and plastic volume strain increment of an element in a calculation cycle, T∗ is the standardized maximum tensile stress, and D1 and D2 are the damage constant of the concrete. The parameters of concrete in this study are shown in Table 2.Table 2
Concrete properties [20].
MidRo (g/cm3)G (Mbar)ABCNFc (Mbar)2.40.1140.791.60.0070.613.00E − 04T (Mbar)Eps0EfminSfmaxPc (Mbar)Uc (Mbar)Pl (Mbar)ul4.00E − 051.00E − 060.0171.60E − 040.0010.0080.1D1D2K1 (Mbar)K2 (Mbar)K3 (Mbar)fs0.0410.85−1.712.080.3


#### 2.3.4. Constitutive Model of Rock

The rock is simulated by means of the Drucker–Prager model. The yield criteria of the Drucker–Prager model are given by(13)t−p⋅tanβ−d=0,(14)t=q2a−brq3,(15)a=1+1K,b=1−1K,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:(16)r=27J321/3,where J2 and J3 represent the second and the third deviatoric stress invariants, respectively.(17)tanβ=6sinφ3−sinφ,(18)d=6Ccosφ3−sinφ.where φ is the angle of friction and C is the dilation angle.The material properties used for quartzite rock mass are given in Table3, obtained from triaxial test [21].Table 3
Quartzite rock properties [21].
DescriptionQuartzite rockWeathering conditionSlightly to moderatelySpecific gravity (G)2.65Density (ρ) (kg/m3)2550Elastic modulus (E) (GPa)28Poisson’s ratio (ν)0.25Angle of internal friction (φ)45°In situ stress ratio (K0)0.5Dilation angle (ψ)5°Cohesion (c) (MPa)2.3σc (MPa)40RQD range75–80RMR47


## 2.1. 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. The 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. The size of whole model is 80 m × 10 m × 42 m. The 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. The charge of the TNT explosive is 40 kg. The 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. Three 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.Figure 1
Isometric side view of geometric model.Figure 2
Equiaxed side view of finite element model.Figure 3
Front view of finite element model.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 Figure4 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. The 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.Figure 4
Details of numerical simulation cases.Table 1
Details of numerical simulation cases.
Numerical simulation case numberHorizontal clearance of lining (m)Numerical simulation case numberHeight difference of lining (m)A15A30A27B32A39C34A411D36A513E38


## 2.2. ALE Multimaterial Fluid Structure Coupling Algorithm

The 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. The 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) The boundary conditions of the modified object remain unchanged, and the internal element is meshed to keep the topological relationship of the mesh unchanged. This step is called smooth step. (2) The 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. This method couples fluid and solid together to realize the nonlinear coupling between fluid medium and rock mass model by keyword∗Constrained_Lagrange_In_Solid. The fluid structure coupling algorithm is to couple the structure and fluid together through certain constraint method to realize the transmission of mechanical parameters. The main constraint methods [17] are velocity constraint, acceleration constraint, and penalty function constraint. The 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. The calculation steps of velocity and acceleration constraints are as follows:The 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)mnM,Ff,i=m0M,Ff,i+himM,fs.Calculate the new acceleration of fluid node (velocity):(2)avf,i=FMf,imfn,i.The acceleration (velocity) of constrained structure node is as follows:(3)avs=∑1ihiavf,i,where mn and mo 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.


## 2.3. Constitutive Model of Materials

### 2.3.1. Constitutive Model of TNT Explosive

Explosive material model is described by keyword∗MAT_HIGH_EXPLOSIVE_BURN combined with JWL equation of state provided by LS-DYNA. The P-V relation of the JWL equation of state is as follows:(4)P=A1−ωR1Ve−R1V+B1−ωR2Ve−R2V+ωEV,where A, B, R1, R2, and ω are material constants for TNT explosive, parameters A and B represent the magnitudes of pressure, ρ is the density of the explosive, and E is the specific internal energy at atmospheric pressure.For TNT  explosive [18], ρ = 1630 kg/m3; detonation wave speed D = 6930 m/s; Chapman–Jouget pressure PCJ=0.255×1011Pa; A = 373800 MPa; B = 3747 MPa; R1 = 4.15; R2 = 0.9; and ω=0.35.


### 2.3.2. Constitutive Model of Air

The air model [19] is described by keyword ∗MAT_NULL combined with multilinear equation of state ∗EOS_LINEAR_POLYNOMIAL in this study.(5)P=C0+C1μ+C2μ2+C3μ2+C4+C5μ+C6μ2e.For the convenience of calculation, air is regarded as an ideal gas, in which the parameters are as follows:C0 = C1 = C2 = C3 = C6 = 0, C4 = C5 = 0.4. For air in cm-g-μs unit system, the parameters of the null material model are as follows: density ρ=0.0012g/cm3 and dynamic viscosity coefficient Mu = 0.001.


### 2.3.3. Constitutive Model of Concrete

The concrete material adopts the model which adapts to high pressure and high strain rate by means of keyword∗MAT_JOHNSON_HOLMQUIST_CONCRETE (HJC). The material model of HJC consists of state equation, yield equation, and damage equation. The 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). The 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. The equations of state of HJC material model are as follows.Elastic loading and unloading (p<pc) is given by(6)p=keμ.Loading of plastic transition zone (pc≤p≤p1) is given by(7)p=p1+p1−pcμ−μ1μ−μ1.Unloading of plastic transition zone (pc≤p≤p1) is given by(8)p−pmax=1−μmax−μcμ1−μcKe+μmax−μcμ1−μcK1μ−μmax.Fully compacted loading (p>p1) is given by (9)p=k1μ−μ11+μ1+k2μ−μ11+μ12+k3μ−μ11+μ13.In the fully compacted unloading (p>p1)stage, the material is completely destroyed:(10)p−pmax=k1μ−μ11+μ1−μ−μ11+μ1max.The yield equation of HJC material model is as follows:(11)σ∗=A1−D+BP∗N1+Clnε∗.The damage equation of HJC material model is as follows:(12)D=∑Δεp+ΔμpD1P∗+T∗D2,where ke=pc/μc is the bulk modulus, which is the ratio of crushing volume pressure pc and crushing volume μc strain in uniaxial compression test, p1 represents compaction pressure of concrete material, and μ1 is the compaction volume strain. K1 is the plastic bulk modulus; pmax 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 hardening index; C represents the strain rate coefficient; Δεp and Δμp are the equivalent plastic strain increment and plastic volume strain increment of an element in a calculation cycle, T∗ is the standardized maximum tensile stress, and D1 and D2 are the damage constant of the concrete. The parameters of concrete in this study are shown in Table 2.Table 2
Concrete properties [20].
MidRo (g/cm3)G (Mbar)ABCNFc (Mbar)2.40.1140.791.60.0070.613.00E − 04T (Mbar)Eps0EfminSfmaxPc (Mbar)Uc (Mbar)Pl (Mbar)ul4.00E − 051.00E − 060.0171.60E − 040.0010.0080.1D1D2K1 (Mbar)K2 (Mbar)K3 (Mbar)fs0.0410.85−1.712.080.3


### 2.3.4. Constitutive Model of Rock

The rock is simulated by means of the Drucker–Prager model. The yield criteria of the Drucker–Prager model are given by(13)t−p⋅tanβ−d=0,(14)t=q2a−brq3,(15)a=1+1K,b=1−1K,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:(16)r=27J321/3,where J2 and J3 represent the second and the third deviatoric stress invariants, respectively.(17)tanβ=6sinφ3−sinφ,(18)d=6Ccosφ3−sinφ.where φ is the angle of friction and C is the dilation angle.The material properties used for quartzite rock mass are given in Table3, obtained from triaxial test [21].Table 3
Quartzite rock properties [21].
DescriptionQuartzite rockWeathering conditionSlightly to moderatelySpecific gravity (G)2.65Density (ρ) (kg/m3)2550Elastic modulus (E) (GPa)28Poisson’s ratio (ν)0.25Angle of internal friction (φ)45°In situ stress ratio (K0)0.5Dilation angle (ψ)5°Cohesion (c) (MPa)2.3σc (MPa)40RQD range75–80RMR47


## 2.3.1. Constitutive Model of TNT Explosive

Explosive material model is described by keyword∗MAT_HIGH_EXPLOSIVE_BURN combined with JWL equation of state provided by LS-DYNA. The P-V relation of the JWL equation of state is as follows:(4)P=A1−ωR1Ve−R1V+B1−ωR2Ve−R2V+ωEV,where A, B, R1, R2, and ω are material constants for TNT explosive, parameters A and B represent the magnitudes of pressure, ρ is the density of the explosive, and E is the specific internal energy at atmospheric pressure.For TNT  explosive [18], ρ = 1630 kg/m3; detonation wave speed D = 6930 m/s; Chapman–Jouget pressure PCJ=0.255×1011Pa; A = 373800 MPa; B = 3747 MPa; R1 = 4.15; R2 = 0.9; and ω=0.35.


## 2.3.2. Constitutive Model of Air

The air model [19] is described by keyword ∗MAT_NULL combined with multilinear equation of state ∗EOS_LINEAR_POLYNOMIAL in this study.(5)P=C0+C1μ+C2μ2+C3μ2+C4+C5μ+C6μ2e.For the convenience of calculation, air is regarded as an ideal gas, in which the parameters are as follows:C0 = C1 = C2 = C3 = C6 = 0, C4 = C5 = 0.4. For air in cm-g-μs unit system, the parameters of the null material model are as follows: density ρ=0.0012g/cm3 and dynamic viscosity coefficient Mu = 0.001.


## 2.3.3. Constitutive Model of Concrete

The concrete material adopts the model which adapts to high pressure and high strain rate by means of keyword∗MAT_JOHNSON_HOLMQUIST_CONCRETE (HJC). The material model of HJC consists of state equation, yield equation, and damage equation. The 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). The 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. The equations of state of HJC material model are as follows.Elastic loading and unloading (p<pc) is given by(6)p=keμ.Loading of plastic transition zone (pc≤p≤p1) is given by(7)p=p1+p1−pcμ−μ1μ−μ1.Unloading of plastic transition zone (pc≤p≤p1) is given by(8)p−pmax=1−μmax−μcμ1−μcKe+μmax−μcμ1−μcK1μ−μmax.Fully compacted loading (p>p1) is given by (9)p=k1μ−μ11+μ1+k2μ−μ11+μ12+k3μ−μ11+μ13.In the fully compacted unloading (p>p1)stage, the material is completely destroyed:(10)p−pmax=k1μ−μ11+μ1−μ−μ11+μ1max.The yield equation of HJC material model is as follows:(11)σ∗=A1−D+BP∗N1+Clnε∗.The damage equation of HJC material model is as follows:(12)D=∑Δεp+ΔμpD1P∗+T∗D2,where ke=pc/μc is the bulk modulus, which is the ratio of crushing volume pressure pc and crushing volume μc strain in uniaxial compression test, p1 represents compaction pressure of concrete material, and μ1 is the compaction volume strain. K1 is the plastic bulk modulus; pmax 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 hardening index; C represents the strain rate coefficient; Δεp and Δμp are the equivalent plastic strain increment and plastic volume strain increment of an element in a calculation cycle, T∗ is the standardized maximum tensile stress, and D1 and D2 are the damage constant of the concrete. The parameters of concrete in this study are shown in Table 2.Table 2
Concrete properties [20].
MidRo (g/cm3)G (Mbar)ABCNFc (Mbar)2.40.1140.791.60.0070.613.00E − 04T (Mbar)Eps0EfminSfmaxPc (Mbar)Uc (Mbar)Pl (Mbar)ul4.00E − 051.00E − 060.0171.60E − 040.0010.0080.1D1D2K1 (Mbar)K2 (Mbar)K3 (Mbar)fs0.0410.85−1.712.080.3


## 2.3.4. Constitutive Model of Rock

The rock is simulated by means of the Drucker–Prager model. The yield criteria of the Drucker–Prager model are given by(13)t−p⋅tanβ−d=0,(14)t=q2a−brq3,(15)a=1+1K,b=1−1K,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:(16)r=27J321/3,where J2 and J3 represent the second and the third deviatoric stress invariants, respectively.(17)tanβ=6sinφ3−sinφ,(18)d=6Ccosφ3−sinφ.where φ is the angle of friction and C is the dilation angle.The material properties used for quartzite rock mass are given in Table3, obtained from triaxial test [21].Table 3
Quartzite rock properties [21].
DescriptionQuartzite rockWeathering conditionSlightly to moderatelySpecific gravity (G)2.65Density (ρ) (kg/m3)2550Elastic modulus (E) (GPa)28Poisson’s ratio (ν)0.25Angle of internal friction (φ)45°In situ stress ratio (K0)0.5Dilation angle (ψ)5°Cohesion (c) (MPa)2.3σc (MPa)40RQD range75–80RMR47


## 3. Dynamic Response Characteristics of Double-Line Tunnel during Surface Blasting

### 3.1. Propagation Law of Stress Wave in Double-Line Tunnel

Figure5 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. The 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.Figure 5
Nephogram of geotechnical pressure at typical time (Mbar). (a)T = 0.5 ms. (b) T = 6.3 ms. (c) T = 10 ms. (d) T = 17 ms.
(a)(b)(c)(d)Figure 6
Layout of monitoring points.It can also be seen from Figure7 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. The 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. This 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.Figure 7
Velocity time-history curve of monitoring points.


### 3.2. Distribution Law of Equivalent Stress of Lining

Figures8 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). This 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 8
Equivalent stress nephogram of lining under condition 1 att = 20 ms (Mbar). (a) A1. (b) A2. (c) A3. (d) A4. (e) A5.
(a)(b)(c)(d)(e)Figure 9
Equivalent stress nephogram of lining under condition 2 att = 20 ms (Mbar). (a) A3. (b) B3. (c) C3. (d) D3. (e) E3.
(a)(b)(c)(d)(e)


### 3.3. Parameter Analysis of Monitoring Points

Figure10 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 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. The 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.Figure 10
Layout of monitoring points.Figure 11
Relationship curve between monitoring points speed and horizontal clear distance of tunnel.As for the monitoring point C on the lining side wall of the back blasting face and the monitoring point D at the bottom of the lining, it can be seen from Figure11 that the vibration velocity of the monitoring point C is lower than that of the monitoring point D, and with the increase of the horizontal clear distance of the tunnel, the vibration velocity of the monitoring points C and D decreases significantly. For the monitoring point C on the back blasting face, with the increase of the horizontal clear distance of the lining, the vibration velocity of the monitoring point C decreases from 6.4 cm/s to 4.7 cm/s, and the peak velocity decreases by about 26.5%. For point D at the bottom of back blasting surface, the vibration velocity of point D decreases from 6.7 cm/s to 5.3 cm/s with the increase of horizontal clear distance of lining, and the peak velocity decreases by about 20.9%.Figures12–14 are the curves of the relationship between the equivalent stress of the monitoring point and the horizontal clear distance of the tunnel, the relationship between the stress of the monitoring point in the X-direction and the horizontal clear distance of the tunnel, and the relationship between the stress of the monitoring point in the Y-direction and the horizontal clear distance of the tunnel, respectively. It can be seen from Figure 12 that, similar to the variation law of velocity monitoring point, the peak value of equivalent stress of lining side wall monitoring point A and lining top point B on the blast facing face does not change with the increase of horizontal spacing of lining. The peak value of equivalent stress of monitoring point B is about 0.42 MPa, and that of monitoring point A is about 0.36 MPa. The peak value of equivalent stress at point A is reduced by about 14%.Figure 12
Relationship curve between equivalent stress of monitoring point and horizontal clear distance of tunnel.Figure 13
Relationship curve betweenX-direction stress of monitoring points and horizontal clear distance of tunnel.Figure 14
Relationship curve betweenY-direction stress of monitoring points and horizontal clear distance of tunnel.For monitoring point C on lining side wall of back blasting face and monitoring point D on lining bottom, from the overall analysis of Figure12, it can be seen that with the increase of horizontal spacing of tunnel, the equivalent stress of monitoring point decreases on the whole. For point C, the peak value of equivalent stress decreases from 0.3 MPa to 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 Figures13 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. This 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 Y-direction increases from 0.15 MPa to 0.18 MPa, an increase of about 20%.Figure15 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. This 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%.Figure 15
Relationship curve between monitoring points speed and horizontal height difference of tunnel.Figures16–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 X-direction stress of the monitoring point and the height difference of the tunnel, and the relationship between the Y-direction 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 X-direction 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 16
Relationship curve between equivalent stress of monitoring point and horizontal height difference of tunnel.Figure 17
Relationship curve between theX-direction stress of monitoring points and tunnel height difference.Figure 18
Relationship curve between theY-direction stress of monitoring points and tunnel height difference.Figure19 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. The 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 19
Path setting of lining monitoring points.Figure 20
Variation of displacement of lining monitoring points along the path in case 1.Figure 21
Variation of displacement of lining monitoring points along the path in case 2.Figure22 is the path layout of rock monitoring points. The monitoring points are arranged clockwise along the z-axis. 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. The 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.Figure 22
Path setting of tunnel rock mass monitoring points.Figure 23
Variation velocity of lining monitoring points along the path in case 1.Figure 24
Variation velocity of lining monitoring points along the path in case 2.


## 3.1. Propagation Law of Stress Wave in Double-Line Tunnel

Figure5 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. The 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.Figure 5
Nephogram of geotechnical pressure at typical time (Mbar). (a)T = 0.5 ms. (b) T = 6.3 ms. (c) T = 10 ms. (d) T = 17 ms.
(a)(b)(c)(d)Figure 6
Layout of monitoring points.It can also be seen from Figure7 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. The 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. This 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.Figure 7
Velocity time-history curve of monitoring points.


## 3.2. Distribution Law of Equivalent Stress of Lining

Figures8 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). This 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 8
Equivalent stress nephogram of lining under condition 1 att = 20 ms (Mbar). (a) A1. (b) A2. (c) A3. (d) A4. (e) A5.
(a)(b)(c)(d)(e)Figure 9
Equivalent stress nephogram of lining under condition 2 att = 20 ms (Mbar). (a) A3. (b) B3. (c) C3. (d) D3. (e) E3.
(a)(b)(c)(d)(e)


## 3.3. Parameter Analysis of Monitoring Points

Figure10 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 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. The 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.Figure 10
Layout of monitoring points.Figure 11
Relationship curve between monitoring points speed and horizontal clear distance of tunnel.As for the monitoring point C on the lining side wall of the back blasting face and the monitoring point D at the bottom of the lining, it can be seen from Figure11 that the vibration velocity of the monitoring point C is lower than that of the monitoring point D, and with the increase of the horizontal clear distance of the tunnel, the vibration velocity of the monitoring points C and D decreases significantly. For the monitoring point C on the back blasting face, with the increase of the horizontal clear distance of the lining, the vibration velocity of the monitoring point C decreases from 6.4 cm/s to 4.7 cm/s, and the peak velocity decreases by about 26.5%. For point D at the bottom of back blasting surface, the vibration velocity of point D decreases from 6.7 cm/s to 5.3 cm/s with the increase of horizontal clear distance of lining, and the peak velocity decreases by about 20.9%.Figures12–14 are the curves of the relationship between the equivalent stress of the monitoring point and the horizontal clear distance of the tunnel, the relationship between the stress of the monitoring point in the X-direction and the horizontal clear distance of the tunnel, and the relationship between the stress of the monitoring point in the Y-direction and the horizontal clear distance of the tunnel, respectively. It can be seen from Figure 12 that, similar to the variation law of velocity monitoring point, the peak value of equivalent stress of lining side wall monitoring point A and lining top point B on the blast facing face does not change with the increase of horizontal spacing of lining. The peak value of equivalent stress of monitoring point B is about 0.42 MPa, and that of monitoring point A is about 0.36 MPa. The peak value of equivalent stress at point A is reduced by about 14%.Figure 12
Relationship curve between equivalent stress of monitoring point and horizontal clear distance of tunnel.Figure 13
Relationship curve betweenX-direction stress of monitoring points and horizontal clear distance of tunnel.Figure 14
Relationship curve betweenY-direction stress of monitoring points and horizontal clear distance of tunnel.For monitoring point C on lining side wall of back blasting face and monitoring point D on lining bottom, from the overall analysis of Figure12, it can be seen that with the increase of horizontal spacing of tunnel, the equivalent stress of monitoring point decreases on the whole. For point C, the peak value of equivalent stress decreases from 0.3 MPa to 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 Figures13 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. This 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 Y-direction increases from 0.15 MPa to 0.18 MPa, an increase of about 20%.Figure15 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. This 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%.Figure 15
Relationship curve between monitoring points speed and horizontal height difference of tunnel.Figures16–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 X-direction stress of the monitoring point and the height difference of the tunnel, and the relationship between the Y-direction 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 X-direction 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 16
Relationship curve between equivalent stress of monitoring point and horizontal height difference of tunnel.Figure 17
Relationship curve between theX-direction stress of monitoring points and tunnel height difference.Figure 18
Relationship curve between theY-direction stress of monitoring points and tunnel height difference.Figure19 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. The 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 19
Path setting of lining monitoring points.Figure 20
Variation of displacement of lining monitoring points along the path in case 1.Figure 21
Variation of displacement of lining monitoring points along the path in case 2.Figure22 is the path layout of rock monitoring points. The monitoring points are arranged clockwise along the z-axis. 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. The 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.Figure 22
Path setting of tunnel rock mass monitoring points.Figure 23
Variation velocity of lining monitoring points along the path in case 1.Figure 24
Variation velocity of lining monitoring points along the path in case 2.


## 4. Conclusions

Based on LS-DYNA 3D nonlinear finite element software, a full coupling model of TNT explosive-surrounding rock-lining structure-air is established. Lagrange algorithm is used for lining and rock mass materials. Euler algorithm is used for air and explosive materials. The 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. The dynamic response characteristics of double-line tunnel under different conditions and surface explosion load are systematically studied. The 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.The 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 theX-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%.



---

*Source: 9226615-2021-07-30.xml*