# Numerical Investigation of the Dynamic Responses and Damage of Linings Subjected to Violent Gas Explosions inside Highway Tunnels

Authors: Zhipeng Li; Shunchuan Wu; Ziqiao Cheng; Yibo Jiang
Journal: Shock and Vibration
 (2018)

---

## Abstract

The linings of structures suffer severe damage when subjected to internal explosions, which cause numerous casualties and incalculable economic losses. In this paper, a violent gas explosion that occurred inside a highway tunnel in the city of Chengdu, China, is studied through numerical simulations. The evaluated energy of the gas explosion was equivalent to 2428.9 kg of TNT. A fully coupled numerical model consisting of five parts is established with dimensions consistent with the real prototype dimensions and by considering fluid-structure interaction (FSI) effects. Then, a detailed modelling process is presented and validated through a comparison with empirical formulas. This paper investigates the strength and propagation characteristics of a blast shock wave inside the tunnel, and both the effective stresses and dynamic responses of the lining are analysed under the blast impact loading. The damage mechanism is studied, and the evolution of the lining damage is reproduced, the results of which show good agreement with the actual conditions. Moreover, in terms of the responses and damage of the lining, the fully coupled blast loading model has obvious advantages in comparison with the simplified blast loading model. Furthermore, the damage assessment of the lining conducted using the single degree of freedom (SDOF) method agrees well with the results of the numerical simulation and site investigations. The comprehensive numerical simulation technique used in the present paper and its results could represent valuable references for future research on violent explosions within tunnels or very large underground structures and provide relevant information for the blast-resistant design of such structures.


---

## Body

## 1. Introduction

Tunnel engineering plays an important role in improving transportation networks. However, a security incident within the relatively confined space of a tunnel could lead to human casualties, substantial property and economic losses, and adverse social consequences. Consequently, an explosion accident constitutes one of the most disastrous types of events. Table1 and Figure 1 describe several typical explosion accidents that occurred inside tunnels throughout China, thereby highlighting the importance of further researching such incidents.Table 1
Explosion accidents within tunnels.
Time and location
Tunnel
Burst source
Casualties
December 22 2005, Sichuan province
Dongjiashan tunnel
Gas explosion
44 people dead
April 8, 2011, Gansu province
Qidaoliang tunnel
Petrol transport vehicles
4 people dead
December 25, 2012, Shanxi province
Nanlvliangshan tunnel
Emulsion explosive
8 people dead
March 1, 2014, Shanxi province
Yanhou tunnel
Methanol transport vehicles
40 people dead
May 5, 2017, Guizhou province
Qishanyan tunnel
Gas explosion
12 people dead
May 23, 2017, Hebei province
Futuyu tunnel
LNG transport vehicles
12 people deadFigure 1
Tunnel explosion accidents. (a) Explosion accident inside the Nanlvliangshan tunnel. (b) Explosion accident inside the Yanhou tunnel.
(a)
(b)Some fruitful studies have been conducted on gas explosions and underground structures subjected to blast loading. For example, some scholars have performed scaled-down centrifuge modelling tests to describe the responses of tunnels under blast loading [1–3]. Research experiments/cases have revealed the different explosion mechanisms/locations and consequently proposed approaches to predict the destructive power of related explosions/seismic sources [4–10]. Li et al. [11] investigated the dynamic response of a masonry structure subjected to an internal gas explosion through 16 full-scale in situ tests and carried out intensive parametric numerical studies on the resistance of such structures to explosions. However, few full-scale explosion experiments have been conducted within a highway tunnel because it is both dangerous and expensive; furthermore, it is nearly impossible to obtain such a large amount of explosives. With the rapid development of computational explosion mechanics and improvements to modern computational performance standards, numerical simulation technologies are now capable of addressing these challenging engineering problems in a timely and cost-effective manner. This technology is also capable of displaying the details of the entire explosion process at each time step. Jayasinghe et al. [12–14] carried out a series of numerical simulations of the dynamic responses of pile structures under blasting loads using ANSYS/LS-DYNA, nonlinear explicit finite element commercial software. Xia et al. [15] analysed the ground vibrations of a water tunnel induced by excavation blasting through experiments and numerical methods. Additionally, Hao et al. [16–19] performed numerous numerical studies on the dynamic responses and damage mechanisms of concrete structures under explosions using the arbitrary Lagrangian–Eulerian (ALE) technique in LS-DYNA and compared the results with field experiments to demonstrate the feasibility of the method. Yang et al. [20] also utilized LS-DYNA to study the dynamic responses of operating metro tunnels subjected to ground explosions. Furthermore, Mobaraki and Vaghefi [21] researched the responses of various buried tunnels under surface explosions from the detonation of 1000 kg of TNT and examined the antiexplosion characteristics of different cross section tunnel shapes. Mussa et al. [22] investigated the dynamic responses of buried rectangular tunnels under the explosions of different magnitudes of TNT and additionally studied the damage degree and resistance of those tunnels to blast loading. Accordingly, the abovementioned studies are important references with respect to the research methods, constitutive models, and calculated parameters of the different materials involved in the numerical simulations. Consequently, those studies provide quantitative data for the responses of structures subjected to blast loading.Tunnel lining plays an important role in the overall stability of a tunnel; thus, it is crucial to study the dynamic responses and damage of lining upon experiencing a violent gas explosion within a tunnel. This paper presents a detailed background description of a tunnel project to quantitatively study gas explosions. Subsequently, a full-scale fluid-structure interaction (FSI) numerical model is established using the ALE technique in LS-DYNA, and proper constitutive models and key parameters are chosen and modified for the different materials modelled. The results of the numerical investigation reveal the propagation of the blast shock wave inside the tunnel and the detailed dynamic responses, damage mechanism and damage processes of the lining, and provide a damage assessment. Finally, the numerical simulation results are validated through a comparison with the findings of field investigations, and empirical methods are presented. The objective of this research, which is highly innovative and has scientific value, is to propose a feasible and effective approach for similar explosion events.


## 2. Project Background

The gas explosion that occurred inside the Luodaiguzhen tunnel is the only typical, recent example of an engineering case of such an incident worldwide. The tunnel, which is located in Chengdu, Sichuan province, China, represents the key component of an important infrastructure project constructed in a public-private partnership for the development of Chengdu and its surrounding districts. The tunnel spans a total length of 2915 m across Longquan Mountain, the strata of which are rich in gas. During a previous phase of surveying and engineering design, the maximal gas overflow reached 0.52 m3/min, classifying this tunnel as a high-gas tunnel [23].A lack of ventilation over a span of 10 days led to the voluminous accumulation of gases. Consequently, on February 24, 2015, a violent explosion was triggered inside the Luodaiguzhen tunnel during the Chinese Lunar New Year holiday when several workers created open flames due to poor management practices and an absence of safety awareness. This major accident killed 8 people and injured approximately 200 more; meanwhile, the linings of the tunnel and nearby residential buildings also suffered damage to different degrees, generating direct economic losses reaching 136 million dollars.Figures2(a) and 2(b) depict the “mushroom cloud” that was produced by the explosion and a scene showing the site conditions and the deployment of emergency rescue operations following the accident. Damaged equipment and construction tools were blasted out of the tunnel; in particular, a loader was fragmented into pieces and blown more than 120 m away, as seen in Figure 2(c). The region of construction and residential buildings near the tunnel entrance were completely destroyed, as shown in Figure 2(d). The above descriptions suggest that the magnitude of the gas explosion within the Luodaiguzhen tunnel was tremendous, and the linings were badly destroyed.Figure 2
Site conditions of the Luodaiguzhen tunnel gas explosion.
(a)
(b)
(c)
(d)The results of a survey conducted at the site of the Luodaiguzhen explosion revealed 3 detonations: one was located at ZK2 + 810 in the left hole, and the other two were located at K2 + 300 and K2 + 587 in the right hole, as shown in Figure3. The linings were damaged as follows. In the region within 5 m of the detonation, the secondary lining completely failed and collapsed; furthermore, the arch and sidewall areas of the first lining were fractured into blocks. In the region between 5 and 10 m of the detonation, the secondary lining fractured into blocks and lost its overall stability, but it did not collapse completely; meanwhile, the first lining in the sidewall area exhibited serious failure. Finally, at distances greater than 10 m from the detonation site, numerous cracks were created throughout the secondary lining dominated by longitudinal cracks with local circumferential cracks. According to these damage traits, the lining can be divided into three zones characterized by complete failure, severe failure, and general failure. This paper performs a numerical investigation of the damaged section (ZK2 + 790∼ZK2 + 830) in the left hole of the tunnel, where the strength grade of the concrete used in the tunnel lining is 25 MPa.Figure 3
Schematic diagrams of the tunnel explosion locations. (a) The left hole of the tunnel. (b) The right hole of the tunnel.
(a)
(b)


## 3. Finite Element Model

### 3.1. Equivalent Gas Explosion Loading

A gas explosion is essentially the rapid ignition of combustible gases when the accumulated concentration of gases reaches a critical value and encounters an open fire. Gas deflagration occurs when it experiences constrained conditions and encounters an obstruction during the spreading of a combustion wavefront. If the constrained conditions are strengthened and the obstructions are greatly increased, then the conversion of gas combustion from deflagration to detonation happens; as a consequence, the constraining structure will suffer enormous damage because the detonation generates a tremendous overpressure [24]. Compared with the overpressure from chemical high explosives, the overpressure from gas is lower but lasts longer [11]. However, the value of the length-to-diameter ratio of the Luodaiguzhen tunnel is high, the construction equipment acting as obstacles are numerous, as seen in Figures 2(c) and 2(d), and accumulated gas inside the tunnel is significantly abundant. These conditions would cause the energy generated by the gas explosion to increase greatly [24]. In general, the gas explosion in the Luodaiguzhen tunnel had a violent detonation.At present, the explosion model of the combustible gas needs to be further studied, and TNT could be used to accurately describe the detonation state of explosives [25, 26]. The TNT-equivalence method suggested by references [27, 28] is employed to evaluate the energies released by various explosions and equate those energies to a certain weight of TNT. This approach can obtain a satisfactory approximation and has been widely used [26]; for example, studies conducted on the explosion of dangerous chemicals at the Port of Tianjin, China, and on terrorist bombing attacks have both used the TNT-equivalence method [29]. Although there are differences between the gas explosion and the TNT-equivalency method, the peak pressure is often overestimated, this is also an effective method to some extent [11]. Consequently, the TNT material is modelled within LS-DYNA; the equivalent TNT weight of the accumulated gases in the Luodaiguzhen tunnel can be calculated as follows.Methane (CH4) is the main combustible component of gas. Accordingly, the TNT-equivalence method is given by(1)Eq=QGQTNT,(2)MTNT=α∗Eq∗VG∗ρG,where Eq is the ratio of the combustion heat of methane to that of TNT; QG is the combustion heat of a unit mass of methane; J/kg; QTNT is the combustion heat of a unit mass of TNT, 4.5 MJ/kg; MTNT is the equivalent weight of TNT; α is the methane concentration in the gas; VG is the volume of the accumulated gas, m3; and ρG is the density of methane, 0.716 kg/m3.This paper considers the most unfavourable situation by assuming that all methane has combusted completely, resulting in a maximum explosive energy. The combustion heat of methane is 55.64 MJ/kg according to the following expression [24]:(3)CH4+2O2⇒CO2+2H2O+890.3kJConsequently, the concentration of methane is 9.5%:(4)C0=20.950.2095+n0%,where n0, the value of which is 2 according to Formula (3), is the number of moles of oxygen required for 1 mol of methane to combust completely.The amount of gas overflow detected in the ZK2 + 790∼ZK2 + 830 region during the construction of the tunnel was 1.08 m³/min. Furthermore, the cross-sectional area of the tunnel was 72.34 m2. Hence, the total volume of the accumulated gas was 2893.6 m³; the equivalent weight of TNT is 2433.6 kg via Formulas (1) and (2).


### 3.2. The Numerical Model

The geometrical shape and size of the numerical model created in this section is entirely consistent with the geometry of the Luodaiguzhen tunnel. Moreover, the explosive charge that formed the blast shock wave is modelled within the fully coupled FSI model; hence, the explosive detonation, the blast shock wave propagation within the tunnel, and the interactions of the linings with the wave are all simulated accordingly. Figures4(a) and 4(b) show the design and cross section profile drawings of the standard tunnel cross section and the numerical model, respectively. The units of r1 to r5 in Figure 4(a) are centimeters. The position of the detonation (i.e., the centre of the explosion) is located at the centre of the tunnel, where the distance to the arch vault and floor are both 3.6 m.Figure 4
The fully coupled FSI numerical model. (a) The design diagram of the standard tunnel cross section. (b) The cross section of the numerical model and the detonation position. (c) Schematic diagram of the numerical research model. (d) Components of the ALE model.
(a)
(b)
(c)
(d)Considering the symmetry of the numerical model demonstrated in Figure4(c), a quarter of the model is established to optimally utilize the computational resources, as shown in Figure 4(d). The translational displacements of the nodes normal to the planes of symmetry are constrained, and nonreflecting boundaries are applied to the outer surfaces (except for the planes of symmetry) to simulate an infinite domain.The preprocessing scheme for the numerical simulation developed in this paper includes the generation of the model and the setting of parameters. First, a 3D geometrical model is established and then meshed to create all of the different components and boundary conditions; this construction process and some other basic steps are carried out in the ANSYS software, after which the key file is output. Second, some parameters (e.g., the constitutive models, equations of state (EOSs) of different material components, and FSI settings) incorporated in the key file are modified in LS-DYNA.Figure4(d) depicts the five components of the FSI finite element model: the TNT source, air, first lining, secondary lining, and surrounding rock. The TNT and air domains are modelled in an Eulerian mesh to remove problems associated with mesh distortions induced by large-scale deformation, while the rest of the components are modelled in a Lagrangian mesh to better understand the responses of the structures when subjected to an internal explosion. The ALE multimaterial mesh model is achieved by using ∗CONSTRAINED_LAGRANGE_IN_SOLID in LS-DYNA [30]. The TNT component is modelled by using ∗INITIAL_VOLUME_FRACTION_GEOMETRY in the air domain, and the equivalent TNT weight is obtained simply by establishing the filled position (i.e., the coordinates of the detonation) and the radius. This method can construct the geometry of the TNT source in the form of a concentrated spherical explosive while simplifying the mesh generation to a large extent. Through refined meshing skills, the elements used to compose all of the components are eight-node solid hexahedron elements, which can improve the calculation efficiency. The mesh size of the elements near the detonation is set to 7 cm to ensure a higher accuracy, and the size increases gradually with the distance away from the detonation. Moreover, elements belonging to different components are aligned at the same positions to significantly improve the coupling effect and achieve substantially better interactions between the blast shock wave and the structures.


### 3.3. The Material Model

#### 3.3.1. TNT

The pressure generated by the equivalent TNT weight, which is modelled using∗MAT_HIGH_EXPLOSIVE_BURN and ∗EOS_JWL in LS-DYNA [30], is given by the following:(5)P=A1−ωR1Ve−R1V+B1−ωR2Ve−R2V+ωEV,where P is the pressure; V is the ratio of the current volume to the initial volume; E is the initial energy; and A, B, R1, R2, and ω are empirical constants based on explosion experiments. The calculation parameters are listed in Table 2.Table 2
Calculation parameters for the TNT component.
ρ (kg/m3)
V
D (m/s)
P
C
−
J (GPa)
A (GPa)
1630
6930
21
3.737 × 102
B (GPa)
R
1
R
2
ω
3.747
4.15
0.9
0.35


#### 3.3.2. Air

The air pressure, which is modelled using∗MAT_NULL in LS-DYNA with a linear polynomial ∗EOS_LINEAR_POLYNOMIAL [30], is given by the following:(6)P=C0+C1μ+C2μ2+C3μ2+C4+C5μ+C6μ2E,where P is the air pressure; C0∼C6 are constants; μ=ρ/ρ0−1, where ρ/ρ0 indicates the ratio of the current density to the initial density; and E is the initial energy per unit volume. The calculation parameters are listed in Table 3.Table 3
Calculation parameters for the air component.
ρ (kg/m3)
C0∼C3
C4∼C5
C6
E
0 (J/m3)
V
0
1.290
0
0.4
0
2.5E5
1


#### 3.3.3. The Lining Structures

The Riedel–Hiermaier–Thoma (RHT) constitutive model, which was proposed by Riedel et al. and given a set of standard parameters [31], has obvious advantages when describing the dynamic characteristics of concrete materials, and it has been validated in the literature [32, 33]. The RHT model takes rate-dependent mechanical behaviours into account and is capable of describing the elastic deformation phase, the linear hardening phase, and the damage-softening phase of concrete under blast impact loading, and these phases are represented through the elastic limit surface, the maximum failure surface, and the residual strength surface, respectively, which are described as follows:(7)σelasticp,θ,ε˙=fc⋅σTXC∗ps,el⋅R3θ⋅Frateε˙⋅Felastic⋅Fcap,σfailp,θ,ε˙=fc⋅σTXC∗ps⋅R3θ⋅Frateε˙,σresidual=Af×p∗Nf,where fc is the quasistatic uniaxial compressive strength of concrete, 25 MPa; ps=p/Frateε˙ is the quasistatic pressure; σTXC∗ps is the equivalent stress strength of the quasistatic compressive meridian; R3θ is a scalar function of the load angle θ and the tensile-to-compressive meridian ratio; Frateε˙ is an influencing factor on the strain rate; Felastic is the scaling factor of the elastic deformation; Fcap is the cap function; and Af and Nf are the residual surface parameters.The compressive and tensile strain rate-dependent exponentsβc and βt are 0.042 and 0.044, respectively, via the following [31]:(8)βc=420+3fc,βt=220+fc.The accumulative damage factorD in the RHT model is defined as follows:(9)D=∑Δεp+ΔμpD1p∗+pspall∗D2,D1p∗+pspall∗D2≥εfmin,where Δεp is the equivalent plastic strain increment; Δμp is the plastic volumetric strain increment; D1 and D2 are the damage parameters with suggested values of 0.015 and 1.0, respectively; and εfmin is 0.0008 [34].The steel ribs and reinforced nets in the first lining effectively improve the tensile strength of the concrete. However, if these reinforcing elements are built into the numerical model, the computational efficiency of the numerical model will significantly decrease. In this paper, the reinforcement provided by the steel ribs and reinforced nets are equivalent to increase the tensile strength of the concrete. Through the following conversion formula, the equivalent tensile strength of the first liningRt is 6.1 MPa [35], while the other calculation parameters of the RHT model for the first lining are consistent with those for the secondary lining, as shown in Table 4:(10)Rt=ft1+fyft⋅ρ,where ft is the tensile strength of concrete, 2.5 MPa; fy is the yield strength of steel, 300 MPa; and ρ is the steel content in the concrete, 1.2%.Table 4
Calculation parameters for the secondary lining.
ρ (kg/m3)
G (GPa)
E
P
S
F
B
0
B
1
T
1 (GPa)
A
N
F
C (MPa)
2314
16.7
2
1.22
1.22
35.27
1.6
0.61
25
f
c
∗
f
t
∗
Q
0
B
T
2
e
0
c
e
0
t
e
c
e
t
0.18
0.1
0.6805
0.0105
0
3E-11
3E-12
3E22
3E22
β
c
β
t
P
T
F
g
c
∗
g
t
∗
X
I
D
1
D
2
E
f
min
4.2E-2
4.4E-2
0.001
0.53
0.7
0.5
0.015
1.0
0.0008
A
f
N
f
A
1 (GPa)
A
2 (GPa)
A
3 (GPa)
P
el (GPa)
P
comp (GPa)
N
α
0
0.7
0.8
35.27
39.58
9.04
0.023
6
3.0
1.188


#### 3.3.4. The Surrounding Rock

The surrounding rock is modelled using∗MAT_PLASTIC_KINEMATIC [30], which can more effectively represent the effect of a high strain rate on the surrounding rock under blast impact loading. The constitutive relation of the surrounding rock can be described as follows:(11)σy=1+ε˙/C1/pσ0+βEpεp,eff,(12)Ep=EtanE0/E0−Etan,where σ0 is the initial yield strength of the rock; σy is the hardening yield strength; β is the hardening parameter, 0≤β≤1; ε˙ is the strain rate; C and p are the strain rate parameters for the Cowper–Symonds strain rate model; εp,eff is the effective plastic strain; Ep is the plasticity hardening modulus; and Etan is the tangent modulus. The calculation parameters for the surrounding rock are shown in Table 5.Table 5
Calculation parameters for the surrounding rock.
ρ (kg/m3)
E (GPa)
P
r
σ
y (MPa)
E
tan (GPa)
2700
23
0.22
5.6
7.0
β
C
P
F
s
V
p
1.0
2.5E-6
4.0
0.01
0


### 3.4. Validation of the Numerical Model

Since the case studied in this paper is related to a violent explosion that occurred within a highway tunnel and that was caused by the detonation of 2428.9 kg of TNT, it is almost impossible to validate the model by comparing it with in situ tests or by adopting existing research findings. However, extensive research has been conducted to propose equations that can predict the overpressure generated by the propagation of a blast shock wave in a free-air field. Accordingly, the representative equations proposed by Henrych and Major [36] and Equations (11) and (12) were constructed based upon numerous experiments:(13)ΔPf=1.40717Z+0.55397Z2−0.03572Z3+0.000625Z4,0.05≤Z≤0.3,0.61938Z−0.03262Z2+0.21324Z3,0.3≤Z≤1,0.0662Z+0.405Z2+0.3288Z3,1≤Z≤10,Z=RW3.where ΔPf is the peak overpressure of a blast shock wave, MPa; Z is the scaled distance, m·kg−1/3; R is the distance from the measurement points to the detonation, m; and W is the weight of the TNT explosive, kg.Hence, the modelling method presented in this paper is validated in this section by comparing the results calculated using these empirical prediction formulas in the form of peak overpressures against the scaled distance (SD). The validation model, which consists of the detonation of 100 kg of TNT and a square free-air field with a width of 4.5 m, is built using the abovementioned modelling techniques with the arrangement of measurement points shown in Figure5(a).Figure 5
Validation model and the results of a comparison. (a) Numerical model and the distribution of the measurement points. (b) Contrasting results of the strength of the shock wave.
(a)
(b)Figure5(b) shows the comparison results, from which the approximate trends of the two curves clearly demonstrate the decrease in the peak overpressure with an increase in the SD, indicating a reliable correlation between the present numerical simulation results and the output from the empirical prediction formulas. These findings also directly validate the modelling techniques used in this paper; this approach can be applied to future research.


## 3.1. Equivalent Gas Explosion Loading

A gas explosion is essentially the rapid ignition of combustible gases when the accumulated concentration of gases reaches a critical value and encounters an open fire. Gas deflagration occurs when it experiences constrained conditions and encounters an obstruction during the spreading of a combustion wavefront. If the constrained conditions are strengthened and the obstructions are greatly increased, then the conversion of gas combustion from deflagration to detonation happens; as a consequence, the constraining structure will suffer enormous damage because the detonation generates a tremendous overpressure [24]. Compared with the overpressure from chemical high explosives, the overpressure from gas is lower but lasts longer [11]. However, the value of the length-to-diameter ratio of the Luodaiguzhen tunnel is high, the construction equipment acting as obstacles are numerous, as seen in Figures 2(c) and 2(d), and accumulated gas inside the tunnel is significantly abundant. These conditions would cause the energy generated by the gas explosion to increase greatly [24]. In general, the gas explosion in the Luodaiguzhen tunnel had a violent detonation.At present, the explosion model of the combustible gas needs to be further studied, and TNT could be used to accurately describe the detonation state of explosives [25, 26]. The TNT-equivalence method suggested by references [27, 28] is employed to evaluate the energies released by various explosions and equate those energies to a certain weight of TNT. This approach can obtain a satisfactory approximation and has been widely used [26]; for example, studies conducted on the explosion of dangerous chemicals at the Port of Tianjin, China, and on terrorist bombing attacks have both used the TNT-equivalence method [29]. Although there are differences between the gas explosion and the TNT-equivalency method, the peak pressure is often overestimated, this is also an effective method to some extent [11]. Consequently, the TNT material is modelled within LS-DYNA; the equivalent TNT weight of the accumulated gases in the Luodaiguzhen tunnel can be calculated as follows.Methane (CH4) is the main combustible component of gas. Accordingly, the TNT-equivalence method is given by(1)Eq=QGQTNT,(2)MTNT=α∗Eq∗VG∗ρG,where Eq is the ratio of the combustion heat of methane to that of TNT; QG is the combustion heat of a unit mass of methane; J/kg; QTNT is the combustion heat of a unit mass of TNT, 4.5 MJ/kg; MTNT is the equivalent weight of TNT; α is the methane concentration in the gas; VG is the volume of the accumulated gas, m3; and ρG is the density of methane, 0.716 kg/m3.This paper considers the most unfavourable situation by assuming that all methane has combusted completely, resulting in a maximum explosive energy. The combustion heat of methane is 55.64 MJ/kg according to the following expression [24]:(3)CH4+2O2⇒CO2+2H2O+890.3kJConsequently, the concentration of methane is 9.5%:(4)C0=20.950.2095+n0%,where n0, the value of which is 2 according to Formula (3), is the number of moles of oxygen required for 1 mol of methane to combust completely.The amount of gas overflow detected in the ZK2 + 790∼ZK2 + 830 region during the construction of the tunnel was 1.08 m³/min. Furthermore, the cross-sectional area of the tunnel was 72.34 m2. Hence, the total volume of the accumulated gas was 2893.6 m³; the equivalent weight of TNT is 2433.6 kg via Formulas (1) and (2).


## 3.2. The Numerical Model

The geometrical shape and size of the numerical model created in this section is entirely consistent with the geometry of the Luodaiguzhen tunnel. Moreover, the explosive charge that formed the blast shock wave is modelled within the fully coupled FSI model; hence, the explosive detonation, the blast shock wave propagation within the tunnel, and the interactions of the linings with the wave are all simulated accordingly. Figures4(a) and 4(b) show the design and cross section profile drawings of the standard tunnel cross section and the numerical model, respectively. The units of r1 to r5 in Figure 4(a) are centimeters. The position of the detonation (i.e., the centre of the explosion) is located at the centre of the tunnel, where the distance to the arch vault and floor are both 3.6 m.Figure 4
The fully coupled FSI numerical model. (a) The design diagram of the standard tunnel cross section. (b) The cross section of the numerical model and the detonation position. (c) Schematic diagram of the numerical research model. (d) Components of the ALE model.
(a)
(b)
(c)
(d)Considering the symmetry of the numerical model demonstrated in Figure4(c), a quarter of the model is established to optimally utilize the computational resources, as shown in Figure 4(d). The translational displacements of the nodes normal to the planes of symmetry are constrained, and nonreflecting boundaries are applied to the outer surfaces (except for the planes of symmetry) to simulate an infinite domain.The preprocessing scheme for the numerical simulation developed in this paper includes the generation of the model and the setting of parameters. First, a 3D geometrical model is established and then meshed to create all of the different components and boundary conditions; this construction process and some other basic steps are carried out in the ANSYS software, after which the key file is output. Second, some parameters (e.g., the constitutive models, equations of state (EOSs) of different material components, and FSI settings) incorporated in the key file are modified in LS-DYNA.Figure4(d) depicts the five components of the FSI finite element model: the TNT source, air, first lining, secondary lining, and surrounding rock. The TNT and air domains are modelled in an Eulerian mesh to remove problems associated with mesh distortions induced by large-scale deformation, while the rest of the components are modelled in a Lagrangian mesh to better understand the responses of the structures when subjected to an internal explosion. The ALE multimaterial mesh model is achieved by using ∗CONSTRAINED_LAGRANGE_IN_SOLID in LS-DYNA [30]. The TNT component is modelled by using ∗INITIAL_VOLUME_FRACTION_GEOMETRY in the air domain, and the equivalent TNT weight is obtained simply by establishing the filled position (i.e., the coordinates of the detonation) and the radius. This method can construct the geometry of the TNT source in the form of a concentrated spherical explosive while simplifying the mesh generation to a large extent. Through refined meshing skills, the elements used to compose all of the components are eight-node solid hexahedron elements, which can improve the calculation efficiency. The mesh size of the elements near the detonation is set to 7 cm to ensure a higher accuracy, and the size increases gradually with the distance away from the detonation. Moreover, elements belonging to different components are aligned at the same positions to significantly improve the coupling effect and achieve substantially better interactions between the blast shock wave and the structures.


## 3.3. The Material Model

### 3.3.1. TNT

The pressure generated by the equivalent TNT weight, which is modelled using∗MAT_HIGH_EXPLOSIVE_BURN and ∗EOS_JWL in LS-DYNA [30], is given by the following:(5)P=A1−ωR1Ve−R1V+B1−ωR2Ve−R2V+ωEV,where P is the pressure; V is the ratio of the current volume to the initial volume; E is the initial energy; and A, B, R1, R2, and ω are empirical constants based on explosion experiments. The calculation parameters are listed in Table 2.Table 2
Calculation parameters for the TNT component.
ρ (kg/m3)
V
D (m/s)
P
C
−
J (GPa)
A (GPa)
1630
6930
21
3.737 × 102
B (GPa)
R
1
R
2
ω
3.747
4.15
0.9
0.35


### 3.3.2. Air

The air pressure, which is modelled using∗MAT_NULL in LS-DYNA with a linear polynomial ∗EOS_LINEAR_POLYNOMIAL [30], is given by the following:(6)P=C0+C1μ+C2μ2+C3μ2+C4+C5μ+C6μ2E,where P is the air pressure; C0∼C6 are constants; μ=ρ/ρ0−1, where ρ/ρ0 indicates the ratio of the current density to the initial density; and E is the initial energy per unit volume. The calculation parameters are listed in Table 3.Table 3
Calculation parameters for the air component.
ρ (kg/m3)
C0∼C3
C4∼C5
C6
E
0 (J/m3)
V
0
1.290
0
0.4
0
2.5E5
1


### 3.3.3. The Lining Structures

The Riedel–Hiermaier–Thoma (RHT) constitutive model, which was proposed by Riedel et al. and given a set of standard parameters [31], has obvious advantages when describing the dynamic characteristics of concrete materials, and it has been validated in the literature [32, 33]. The RHT model takes rate-dependent mechanical behaviours into account and is capable of describing the elastic deformation phase, the linear hardening phase, and the damage-softening phase of concrete under blast impact loading, and these phases are represented through the elastic limit surface, the maximum failure surface, and the residual strength surface, respectively, which are described as follows:(7)σelasticp,θ,ε˙=fc⋅σTXC∗ps,el⋅R3θ⋅Frateε˙⋅Felastic⋅Fcap,σfailp,θ,ε˙=fc⋅σTXC∗ps⋅R3θ⋅Frateε˙,σresidual=Af×p∗Nf,where fc is the quasistatic uniaxial compressive strength of concrete, 25 MPa; ps=p/Frateε˙ is the quasistatic pressure; σTXC∗ps is the equivalent stress strength of the quasistatic compressive meridian; R3θ is a scalar function of the load angle θ and the tensile-to-compressive meridian ratio; Frateε˙ is an influencing factor on the strain rate; Felastic is the scaling factor of the elastic deformation; Fcap is the cap function; and Af and Nf are the residual surface parameters.The compressive and tensile strain rate-dependent exponentsβc and βt are 0.042 and 0.044, respectively, via the following [31]:(8)βc=420+3fc,βt=220+fc.The accumulative damage factorD in the RHT model is defined as follows:(9)D=∑Δεp+ΔμpD1p∗+pspall∗D2,D1p∗+pspall∗D2≥εfmin,where Δεp is the equivalent plastic strain increment; Δμp is the plastic volumetric strain increment; D1 and D2 are the damage parameters with suggested values of 0.015 and 1.0, respectively; and εfmin is 0.0008 [34].The steel ribs and reinforced nets in the first lining effectively improve the tensile strength of the concrete. However, if these reinforcing elements are built into the numerical model, the computational efficiency of the numerical model will significantly decrease. In this paper, the reinforcement provided by the steel ribs and reinforced nets are equivalent to increase the tensile strength of the concrete. Through the following conversion formula, the equivalent tensile strength of the first liningRt is 6.1 MPa [35], while the other calculation parameters of the RHT model for the first lining are consistent with those for the secondary lining, as shown in Table 4:(10)Rt=ft1+fyft⋅ρ,where ft is the tensile strength of concrete, 2.5 MPa; fy is the yield strength of steel, 300 MPa; and ρ is the steel content in the concrete, 1.2%.Table 4
Calculation parameters for the secondary lining.
ρ (kg/m3)
G (GPa)
E
P
S
F
B
0
B
1
T
1 (GPa)
A
N
F
C (MPa)
2314
16.7
2
1.22
1.22
35.27
1.6
0.61
25
f
c
∗
f
t
∗
Q
0
B
T
2
e
0
c
e
0
t
e
c
e
t
0.18
0.1
0.6805
0.0105
0
3E-11
3E-12
3E22
3E22
β
c
β
t
P
T
F
g
c
∗
g
t
∗
X
I
D
1
D
2
E
f
min
4.2E-2
4.4E-2
0.001
0.53
0.7
0.5
0.015
1.0
0.0008
A
f
N
f
A
1 (GPa)
A
2 (GPa)
A
3 (GPa)
P
el (GPa)
P
comp (GPa)
N
α
0
0.7
0.8
35.27
39.58
9.04
0.023
6
3.0
1.188


### 3.3.4. The Surrounding Rock

The surrounding rock is modelled using∗MAT_PLASTIC_KINEMATIC [30], which can more effectively represent the effect of a high strain rate on the surrounding rock under blast impact loading. The constitutive relation of the surrounding rock can be described as follows:(11)σy=1+ε˙/C1/pσ0+βEpεp,eff,(12)Ep=EtanE0/E0−Etan,where σ0 is the initial yield strength of the rock; σy is the hardening yield strength; β is the hardening parameter, 0≤β≤1; ε˙ is the strain rate; C and p are the strain rate parameters for the Cowper–Symonds strain rate model; εp,eff is the effective plastic strain; Ep is the plasticity hardening modulus; and Etan is the tangent modulus. The calculation parameters for the surrounding rock are shown in Table 5.Table 5
Calculation parameters for the surrounding rock.
ρ (kg/m3)
E (GPa)
P
r
σ
y (MPa)
E
tan (GPa)
2700
23
0.22
5.6
7.0
β
C
P
F
s
V
p
1.0
2.5E-6
4.0
0.01
0


## 3.3.1. TNT

The pressure generated by the equivalent TNT weight, which is modelled using∗MAT_HIGH_EXPLOSIVE_BURN and ∗EOS_JWL in LS-DYNA [30], is given by the following:(5)P=A1−ωR1Ve−R1V+B1−ωR2Ve−R2V+ωEV,where P is the pressure; V is the ratio of the current volume to the initial volume; E is the initial energy; and A, B, R1, R2, and ω are empirical constants based on explosion experiments. The calculation parameters are listed in Table 2.Table 2
Calculation parameters for the TNT component.
ρ (kg/m3)
V
D (m/s)
P
C
−
J (GPa)
A (GPa)
1630
6930
21
3.737 × 102
B (GPa)
R
1
R
2
ω
3.747
4.15
0.9
0.35


## 3.3.2. Air

The air pressure, which is modelled using∗MAT_NULL in LS-DYNA with a linear polynomial ∗EOS_LINEAR_POLYNOMIAL [30], is given by the following:(6)P=C0+C1μ+C2μ2+C3μ2+C4+C5μ+C6μ2E,where P is the air pressure; C0∼C6 are constants; μ=ρ/ρ0−1, where ρ/ρ0 indicates the ratio of the current density to the initial density; and E is the initial energy per unit volume. The calculation parameters are listed in Table 3.Table 3
Calculation parameters for the air component.
ρ (kg/m3)
C0∼C3
C4∼C5
C6
E
0 (J/m3)
V
0
1.290
0
0.4
0
2.5E5
1


## 3.3.3. The Lining Structures

The Riedel–Hiermaier–Thoma (RHT) constitutive model, which was proposed by Riedel et al. and given a set of standard parameters [31], has obvious advantages when describing the dynamic characteristics of concrete materials, and it has been validated in the literature [32, 33]. The RHT model takes rate-dependent mechanical behaviours into account and is capable of describing the elastic deformation phase, the linear hardening phase, and the damage-softening phase of concrete under blast impact loading, and these phases are represented through the elastic limit surface, the maximum failure surface, and the residual strength surface, respectively, which are described as follows:(7)σelasticp,θ,ε˙=fc⋅σTXC∗ps,el⋅R3θ⋅Frateε˙⋅Felastic⋅Fcap,σfailp,θ,ε˙=fc⋅σTXC∗ps⋅R3θ⋅Frateε˙,σresidual=Af×p∗Nf,where fc is the quasistatic uniaxial compressive strength of concrete, 25 MPa; ps=p/Frateε˙ is the quasistatic pressure; σTXC∗ps is the equivalent stress strength of the quasistatic compressive meridian; R3θ is a scalar function of the load angle θ and the tensile-to-compressive meridian ratio; Frateε˙ is an influencing factor on the strain rate; Felastic is the scaling factor of the elastic deformation; Fcap is the cap function; and Af and Nf are the residual surface parameters.The compressive and tensile strain rate-dependent exponentsβc and βt are 0.042 and 0.044, respectively, via the following [31]:(8)βc=420+3fc,βt=220+fc.The accumulative damage factorD in the RHT model is defined as follows:(9)D=∑Δεp+ΔμpD1p∗+pspall∗D2,D1p∗+pspall∗D2≥εfmin,where Δεp is the equivalent plastic strain increment; Δμp is the plastic volumetric strain increment; D1 and D2 are the damage parameters with suggested values of 0.015 and 1.0, respectively; and εfmin is 0.0008 [34].The steel ribs and reinforced nets in the first lining effectively improve the tensile strength of the concrete. However, if these reinforcing elements are built into the numerical model, the computational efficiency of the numerical model will significantly decrease. In this paper, the reinforcement provided by the steel ribs and reinforced nets are equivalent to increase the tensile strength of the concrete. Through the following conversion formula, the equivalent tensile strength of the first liningRt is 6.1 MPa [35], while the other calculation parameters of the RHT model for the first lining are consistent with those for the secondary lining, as shown in Table 4:(10)Rt=ft1+fyft⋅ρ,where ft is the tensile strength of concrete, 2.5 MPa; fy is the yield strength of steel, 300 MPa; and ρ is the steel content in the concrete, 1.2%.Table 4
Calculation parameters for the secondary lining.
ρ (kg/m3)
G (GPa)
E
P
S
F
B
0
B
1
T
1 (GPa)
A
N
F
C (MPa)
2314
16.7
2
1.22
1.22
35.27
1.6
0.61
25
f
c
∗
f
t
∗
Q
0
B
T
2
e
0
c
e
0
t
e
c
e
t
0.18
0.1
0.6805
0.0105
0
3E-11
3E-12
3E22
3E22
β
c
β
t
P
T
F
g
c
∗
g
t
∗
X
I
D
1
D
2
E
f
min
4.2E-2
4.4E-2
0.001
0.53
0.7
0.5
0.015
1.0
0.0008
A
f
N
f
A
1 (GPa)
A
2 (GPa)
A
3 (GPa)
P
el (GPa)
P
comp (GPa)
N
α
0
0.7
0.8
35.27
39.58
9.04
0.023
6
3.0
1.188


## 3.3.4. The Surrounding Rock

The surrounding rock is modelled using∗MAT_PLASTIC_KINEMATIC [30], which can more effectively represent the effect of a high strain rate on the surrounding rock under blast impact loading. The constitutive relation of the surrounding rock can be described as follows:(11)σy=1+ε˙/C1/pσ0+βEpεp,eff,(12)Ep=EtanE0/E0−Etan,where σ0 is the initial yield strength of the rock; σy is the hardening yield strength; β is the hardening parameter, 0≤β≤1; ε˙ is the strain rate; C and p are the strain rate parameters for the Cowper–Symonds strain rate model; εp,eff is the effective plastic strain; Ep is the plasticity hardening modulus; and Etan is the tangent modulus. The calculation parameters for the surrounding rock are shown in Table 5.Table 5
Calculation parameters for the surrounding rock.
ρ (kg/m3)
E (GPa)
P
r
σ
y (MPa)
E
tan (GPa)
2700
23
0.22
5.6
7.0
β
C
P
F
s
V
p
1.0
2.5E-6
4.0
0.01
0


## 3.4. Validation of the Numerical Model

Since the case studied in this paper is related to a violent explosion that occurred within a highway tunnel and that was caused by the detonation of 2428.9 kg of TNT, it is almost impossible to validate the model by comparing it with in situ tests or by adopting existing research findings. However, extensive research has been conducted to propose equations that can predict the overpressure generated by the propagation of a blast shock wave in a free-air field. Accordingly, the representative equations proposed by Henrych and Major [36] and Equations (11) and (12) were constructed based upon numerous experiments:(13)ΔPf=1.40717Z+0.55397Z2−0.03572Z3+0.000625Z4,0.05≤Z≤0.3,0.61938Z−0.03262Z2+0.21324Z3,0.3≤Z≤1,0.0662Z+0.405Z2+0.3288Z3,1≤Z≤10,Z=RW3.where ΔPf is the peak overpressure of a blast shock wave, MPa; Z is the scaled distance, m·kg−1/3; R is the distance from the measurement points to the detonation, m; and W is the weight of the TNT explosive, kg.Hence, the modelling method presented in this paper is validated in this section by comparing the results calculated using these empirical prediction formulas in the form of peak overpressures against the scaled distance (SD). The validation model, which consists of the detonation of 100 kg of TNT and a square free-air field with a width of 4.5 m, is built using the abovementioned modelling techniques with the arrangement of measurement points shown in Figure5(a).Figure 5
Validation model and the results of a comparison. (a) Numerical model and the distribution of the measurement points. (b) Contrasting results of the strength of the shock wave.
(a)
(b)Figure5(b) shows the comparison results, from which the approximate trends of the two curves clearly demonstrate the decrease in the peak overpressure with an increase in the SD, indicating a reliable correlation between the present numerical simulation results and the output from the empirical prediction formulas. These findings also directly validate the modelling techniques used in this paper; this approach can be applied to future research.


## 4. Results and Discussion

### 4.1. The Strength and Propagation Characteristics of the Blast Shock Wave inside the Tunnel

Figure6(a) depicts the arrangement of measurement points along the radial and longitudinal directions of the tunnel. Six points are equidistantly distributed within 6 m of the detonation in both directions, while the rest of the points are arranged along the longitudinal direction at a 2 m interval. Figure 6(b) shows the maximum peak overpressure curves of the points in the form of the overpressure against the SD.Figure 6
Arrangement of the measurement points and the maximum peak overpressure curves. (a) Arrangement of the measurement points. (b) Peak overpressure curves of the measurement points.
(a)
(b)The following can be discerned from Figure6(b): (1) the strength of the shock wave attenuates significantly in the vicinity of the detonation; (2) the maximum peak overpressures at the points distributed along both directions are basically equal near the detonation site and decrease gradually with an increase in the distance from the detonation, but the variation trends of the curves are different; (3) the radial curve exhibits clear fluctuations, especially with regard to the abrupt, significant increase in the peak values at the points close to the corner of the sidewall, while the longitudinal curve is relatively stable as the SD increases, although the peak values of the points are still high (approximately 3∼5 MPa). Accordingly, it is necessary to investigate the propagation characteristics of the blast shock wave inside the tunnel.Figure7 illustrates the detailed propagation process of the shock wave inside the tunnel (the units in the diagrams are 102 GPa).Figure 7
Propagation process of the blast shock wave inside the tunnel. (a) 0 ms. (b) 0.5 ms. (c) 0.98 ms. (d) 1.3 ms. (e) 2.2 ms. (f) 5 ms. (g) 7 ms. (h) 9.3 ms. (i) 12 ms.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Note that the pressure inside the tunnel (0.1 MPa) is equal to the standard atmospheric pressure at 0 ms. The shock wave expands in the form of a spherical wave during the initial state of the explosion, but the overpressure decreases quickly and decays from 153.7 MPa to 12.9 MPa with an increase in the area of the wavefront, as shown in Figures7(a) and 7(c).During the interaction between the shock wave and the concrete lining, the incident wave reflects off the lining due to the difference in the acoustic impedance between the concrete and the air, leading to a significant increase in the strength of the shock wave. The strongest reflection effect occurs in the corner of the sidewall, where the pressure increases to 47.89 MPa, as seen in Figures7(d) and 7(e).During the explosion, a negative pressure region is formed in the vicinity of the detonation site as the air burns out. Meanwhile, the reflected shock wave propagates towards the detonation site and propagates along the longitudinal direction of the tunnel coincident with the original incident shock wave accompanied by a gradual decrease in the negative pressure region. These phenomena can be clearly observed in Figures7(f) and 7(g).Figures7(h) and 7(i) demonstrate the eventual convergence of the reflected shock wave along the longitudinal axis of the tunnel and the disappearance of the negative pressure region; after these events, the shock waves continue to propagate back and forth along their own trajectories.According to the above numerical findings, the pressure field is highly complicated, and the overpressure inside the tunnel remains high due to the arbitrary reflections of the shock waves, further indicating that the numerical method adopted in this paper can effectively reveal the details of the explosion event inside the tunnel. Meanwhile, it is very important to investigate the responses of the lining subjected to severe blast impact loading.


### 4.2. The Effective Stress Propagation in the Lining

Figure8 depicts the contours of the effective stress in the tunnel lining subjected to the explosion (the units in the diagrams are 102 GPa). The transition of the colours from blue to red represents an increase in the effective stress. The blast impact loading, which is transferred from the air in the form of a spherical shock wave as shown in Figure 8(a), simultaneously acts on the floor and vault of the lining at 0.98 ms. During the interaction between the blast shock wave and the lining, the stress propagates along the circumferential and longitudinal directions of the tunnel lining until the entire lining is subjected to a high blast impact loading. The concentration of stress primarily occurs at the corner of the sidewall due to the “horn structure,” which significantly strengthens the reflection effect and consequently increases the strength of the blast shock wave. Figure 8 also demonstrates that the stress acting on the lining does not decrease as the duration of the impact loading increases, attributed to the complicated propagation features of the blast shock wave inside the tunnel, as previously mentioned. Moreover, the lining in the vicinity of the detonation site suffers from considerable permanent deformation, which is evident in Figure 8(f). Hereafter, this paper will focus on the dynamic responses of the lining under high blast loading.Figure 8
Effective stress cloud charts of the tunnel lining. (a) 0.98 ms. (b) 1.3 ms. (c) 2.2 ms. (d) 3.6 ms. (e) 8.2 ms. (f) 15.6 ms.
(a)
(b)
(c)
(d)
(e)
(f)


### 4.3. The Dynamic Responses of the Lining

In this section, some measurement points are selected on the lining along the circumferential and longitudinal directions to investigate the dynamic responses of the lining when subjected to a blast impact loading. Figure9 illustrates the schematic arrangement of the measurement points. Along the circumferential direction (Figure 9(a)), points 1∼5 divide the arch of the lining into 4 equal parts (point 3 is located at the middle of the arch), points 5∼7 divide the sidewall of the lining into 2 equal parts (point 6 is located at the middle of the sidewall), and points 7∼9 divide the floor of the lining into 2 equal parts. Along the longitudinal direction, the points within 10 m of the detonation site are arranged at a 1 m interval, and the points at a distance greater than 10 m away from the explosion centre are arranged at an interval of 2 m, as shown in Figure 9(b). The dynamic responses of the lining are analysed in terms of the velocity, displacement, and acceleration of the measurement points, the numerical results of which are presented hereafter.Figure 9
Schematics of the measurement points on the lining. (a) Cross section view. (b) Longitudinal cross section view.
(a)
(b)Figure10 clearly reveals that the peak values of the three indexes of all of the measurement points are relatively high in the region within 5 m of the detonation site because the blast impact loading is the strongest immediately adjacent to the detonation site. The maximum values of the resultant velocity, displacement, and acceleration at point 1 (located at the vault of the arch) resulting from the blast shock wave that arrived first at the vault are 57.5 m/s, 39.1 cm, and 18482.7 g, respectively.Figure 10
Dynamic responses of the lining arch. (a) Velocity responses. (b) Displacement responses. (c) Acceleration responses.
(a)
(b)
(c)The trends of the velocity curves of the points are all consistent with those of the displacement curves, and the peak displacements of all of the measurement points are greater than 6.5 cm in the region within 5 m of the detonation site. This indicates that the arch of the lining in this region suffered significant deformation (i.e., outward expansion) upon being subjected to violent blast impact loading, as illustrated in Figure10(b). The peak velocities and displacements at point 5 are relatively higher than those at the other points at distances within 5∼14 m of the explosion centre, as shown in Figures 10(a) and 10(b). The peak accelerations at point 3 are larger than those at the other points at distances within 3∼10 m of the detonation, indicating that the middle part of the arch experienced a high blast impact loading, as seen in Figure 10(c), which was caused by the complicated reflection effect of the blast shock wave in this region. Due to the decrease in the blast loading with an increase in the distance from the detonation centre, all of the peak values of these factors decay rapidly; moreover, the general trends of the peak curves presented in Figure 10 tend to flatten.The variation curves of the dynamic responses of the sidewall that are clearly presented in Figure11 demonstrate that the lower part of the sidewall clearly exhibited a greater response than the upper part. The maximum peak velocity and acceleration at point 7 (Figure 9), where the “horn structure” at the corner of the sidewall significantly increases the reflection effect and consequently strengthens the blast shock wave, are 32.4 m/s and 8648 g, respectively. The maximum displacement at point 6 is 38.2 cm, which is 95.5% of the thickness of the structure, thereby demonstrating that the middle of the sidewall suffers a severe deflection, leading to complete failure.Figure 11
Dynamic responses of the lining sidewall. (a) Velocity responses. (b) Displacement responses. (c) Acceleration responses.
(a)
(b)
(c)The velocity and displacement curves of point 9 clearly oscillate, as shown in Figure12. The reasons for this oscillation may be twofold: (1) the blast impact loading first acts on point 9, which is the projected position of the detonation centre onto the tunnel floor, during the propagation of the blast shock wave identical to point 1(Figure 9), and (2) the tunnel drain channel is located beneath point 9, and this hollow structure complicates both the shock wave and the impact loading at the position of point 9. The maximum peak displacement at point 9 is 28.4 cm downward, which is marginally lower than that at point 7. The peak acceleration at point 8 is larger than those of the other points at distances within 10 m of the detonation site, indicating that the impact loading acting on this point is larger.Figure 12
Dynamic responses of the lining floor. (a) Velocity responses. (b) Displacement responses. (c) Acceleration responses.
(a)
(b)
(c)


### 4.4. The Damage to the Lining

#### 4.4.1. The Damage Mechanism of the Lining

Cloud charts of the damage traits of the lining after being subjected to the blast impact loading are clearly presented in Figure13, in which a damage level of 1 indicates that the structure has been completely damaged, and a level of 0 denotes a fully intact structure. The lining in the vicinity of the detonation clearly experienced large-scale damage; in addition, the arch and sidewall expanded outward substantially, while the floor moved downward perceptibly. Far away from the detonation, the failure traits of the lining are primarily similar to those of the lining in the local damage area, particularly in the corner of the sidewall and the floor, and dense cracks are outspread over the arch. These characteristics are in accordance with the previously analysed dynamic responses of the lining.Figure 13
Overall degree of damage of the tunnel lining. (a) Overall damage situation. (b) Details of the cross section. (c) Details of the tunnel floor.
(a)
(b)
(c)Figure14 shows the relationship curves of the peak principal stress against the distance to the detonation centre at all measurement points presented in Figure 9. The positive and negative values of the peak principal stresses illustrated in Figure 14 represent the tensile and compressive stresses acting on the lining, respectively.Figure 14
Stress states of the lining.Figure14 demonstrates that the tensile stresses acting on the lining are high (approximately 5 MPa) and that the trends of the curves are stable (i.e., straight lines). However, the compressive stresses experienced by the lining show a wide range of changes; the values in the vicinity of the detonation site are the highest, and the stresses tend to gradually stabilize with an increase in the distance.The damage mechanism of the lining can be concluded as follows. (1) A local failure zone was formed in the lining at distances within 5.1 m from the detonation site when it was subjected to the significant stress exerted by the blast shock wave; then, under the alternating actions of the strongest compressive and tensile stresses, the failure zone rapidly expanded and eventually developed into a zone of complete failure. (2) In the region between 5.1 and 10.9 m from the explosion centre, the dense cracks within the lining propagated and became interconnected under higher compressive and tensile stresses, and the intact lining was divided into blocky structures, thereby leading to the loss of its overall stability and the formation of a severely damaged zone. (3) In the region beyond 10.9 m from the detonation site, a generally damaged zone was formed, in which only cracks formed throughout the lining.


#### 4.4.2. The Damage Process of the Lining

To understand the entire damage process of the lining in detail, this paper utilizes the erosion algorithm for components in LS-DYNA; this algorithm has the ability to simulate the physical failure process of the lining by deleting the components that satisfy specific erosion criteria. The peak tensile strain of concrete is approximately 2E-4 under quasistatic loading, while the strain rate and the corresponding strength amplification factor can exceed 10∼100 s−1 and 5.0 under blast impact loading, respectively. However, the enhancement factor of the failure strain is less than that of the strength factor; thus, in this section, a tensile failure strain of 1E-3 is adopted for the concrete as the eroding plastic strain, meaning that a component of the lining fails and is deleted upon reaching this threshold value [37].Figure15 illustrates the numerical simulation of the full damage process of the lining and provides comparisons with the results of field investigations. A locally damaged zone first appeared in the arch and the middle of the lining floor; then, a failure zone formed in the corner of the sidewall due to the “horn structure” that significantly strengthened the blast impact loading. Hence, the locally damaged zone expanded along the longitudinal and circumferential directions of the lining before finally developing into a completely damaged zone, as clearly shown in Figures 15(a)–15(d). The strength of the blast loading decreased during the propagation of the shock wave, but the complete failure of the lining in the vicinity of the detonation aggravated the progressive failure of the lining far away from the detonation, resulting in the formation of a severely damaged zone between 5.1 and 10.9 m from the detonation site, at which distances the intact lining was cut into blocks of various sizes as demonstrated in Figures 15(d)–15(e). For the zone situated more than 10.9 m from the explosion centre, cracks were produced that propagated across the arch under high tensile stresses, as shown in Figures 15(e)–15(f); hence, a generally damaged zone developed. Comparisons of the damage features between the numerical and field investigations are illustrated in Figure 15(g), from which it is evident that the results of the numerical simulation agree well with the actual situation and that the research methods used in this paper are both feasible and effective.Figure 15
Simulations of the damage and comparison with the actual lining situation. (a) 1.5 ms. (b) 2 ms. (c) 2.5 ms. (d) 4 ms. (e) 8 ms. (f) 20 ms. (g) Damage details of the numerical and site investigations.
(a)
(b)
(c)
(d)
(e)
(f)
(g)


### 4.5. Discussion

#### 4.5.1. The Blast Impact Loading

In this paper, the dynamic responses and damage characteristics of the lining are directly related to the mode of the blast impact loading. The abovementioned research findings are based on the FSI numerical model. The keys of this model are that the modelled explosive charge forms the blast shock wave and that the lining is subjected to a fully coupled blast loading. Figure16 depicts both the pressure-time curves of the shock wave that propagates through the air and strikes an obstacle and the simplified triangle impact loading curve. Both types of curves account for the enhanced strength of the shock wave due to reflection effects, and the reflected pressure exhibits a certain relationship with the incident pressure in the diagram. The simplified loading curve can conserve considerably more computational resources than the typical reflected curve and has thus been adopted by many researchers [38–40].Figure 16
P-T curves of the blast loading.The P-T curves, which are described by the fully coupled FSI model, of selected measurement points (Figure17(a)) on the lining are presented in Figures 17(b) and 17(c).Figure 17
P-T curves of different positions on the lining described by the numerical simulation. (a) Locations of the measurement points on the lining. (b) P-T curves of the measurement points on the same cross section. (c) P-T curves of the measurement points at the same positions.
(a)
(b)
(c)The curves in Figures16 and 17 are dissimilar; furthermore, the P-T curves of the points at various locations on the lining shown in Figure 17 are also different both with regard to the trends of the curves and the peak pressures. All of the curves clearly experience multiple peak values during 0∼4 ms, which can be explained by the reflection effect of the blast shock wave. Moreover, lines A and D and lines 3 and 4 also exhibit several peak pressures after 6 ms. The reflected peak pressures of the measurement points are hardly describable using a given relationship with the incident pressure; one of the main reasons for this might be the variations in the shock wave field caused by the reflections at different locations. The fully coupled blast loading model can represent all of the reflection effects at various locations as the wave propagates inside the tunnel, while the simplified loading model cannot. The wave transformed from a compressive wave to a tensile wave; under the tensile stress, cracks emerged throughout the lining, and the damaged lining broke into blocks of concrete, forming zones of severe and general damage. This result is in accordance with the findings of Section 4.4. Based on the analysis in this section, use of the fully coupled blast loading model has obvious advantages for the research object of this paper.


#### 4.5.2. The Damage Assessment of the Lining

The single degree of freedom (SDOF) model proposed by Fallah and Louca [41], which simplifies the structure into an equivalent perfect elastic-plastic model, is applied to quantitatively study and assess the damage of all types of structures, especially for the column structure [42, 43]. The SDOF method is used to assess the damage levels of an underground tunnel when undergoing several explosion cases [22]. In this section, this effective method is used for reference to estimate the damage to the lining structure. The criterion is based on the maximum deflection at the mid-height of the structure, and different deflections correspond to different damage levels. The critical values of the deflections differentiating various damage levels contained in the SDOF model are given in Table 6.Table 6
Damage assessment criterion of the SDOF approach.
Damage level
Maximum displacement (cm)
Low damage (LD)
<2
Moderate damage (MD)
2∼4
High damage (HD)
4∼8
Full damage (FD)
>8According to the SDOF approach, the positions of points 3, 6, and 9 in Figure9(a) represent the arch, the sidewall, and the floor of the lining, respectively, as shown in Figure 18(a). The maximum deflection curves of these measured points are shown in Figure 18(b).Figure 18
Evaluation of the lining damage based on the SDOF method. (a) Arrangement of the measurement points. (b) Maximum deflection curves.
(a)
(b)Figure18(b) reveals that the deflection curves of the arch and sidewall show decaying trends as the distance to the detonation site increases, while the deflection curve for the floor exhibits relative fluctuations related to the propagation of the blast shock wave and the previously mentioned geometric construction of the tunnel floor. Table 7 describes the detailed damage characteristics for the arch, sidewall, and floor of the lining. The maximum deflections exceed 8 cm in the arch and sidewall corresponding to full damage (FD) in the region within 5 m of the detonation site, while the degree of damage ranges from FD to moderate damage (MD) in the region between 6 and 10 m from the explosion centre, and low damage (LD) occurs when the distance is greater than 10 m. Meanwhile, the general trend of the floor curve is attenuated. In summary, it is clearly evident that the degrees of the damage of the lining evaluated using the SDOF method agree well with the site conditions and the numerical results.Table 7
Damage characteristics of the lining.
Damage degree
Parts of the lining
Distance to the detonation (m)
0
1
2
3
4
5
6
7
8
9
10
12
14
16
18
20
Arch
FD
FD
FD
FD
FD
FD
HD
HD
MD
MD
MD
LD
LD
LD
LD
LD
Sidewall
FD
FD
FD
FD
FD
FD
FD
HD
HD
HD
MD
LD
LD
LD
LD
LD
Floor
FD
FD
HD
MD
HD
HD
FD
FD
FD
HD
MD
LD
MD
MD
LD
LD


## 4.1. The Strength and Propagation Characteristics of the Blast Shock Wave inside the Tunnel

Figure6(a) depicts the arrangement of measurement points along the radial and longitudinal directions of the tunnel. Six points are equidistantly distributed within 6 m of the detonation in both directions, while the rest of the points are arranged along the longitudinal direction at a 2 m interval. Figure 6(b) shows the maximum peak overpressure curves of the points in the form of the overpressure against the SD.Figure 6
Arrangement of the measurement points and the maximum peak overpressure curves. (a) Arrangement of the measurement points. (b) Peak overpressure curves of the measurement points.
(a)
(b)The following can be discerned from Figure6(b): (1) the strength of the shock wave attenuates significantly in the vicinity of the detonation; (2) the maximum peak overpressures at the points distributed along both directions are basically equal near the detonation site and decrease gradually with an increase in the distance from the detonation, but the variation trends of the curves are different; (3) the radial curve exhibits clear fluctuations, especially with regard to the abrupt, significant increase in the peak values at the points close to the corner of the sidewall, while the longitudinal curve is relatively stable as the SD increases, although the peak values of the points are still high (approximately 3∼5 MPa). Accordingly, it is necessary to investigate the propagation characteristics of the blast shock wave inside the tunnel.Figure7 illustrates the detailed propagation process of the shock wave inside the tunnel (the units in the diagrams are 102 GPa).Figure 7
Propagation process of the blast shock wave inside the tunnel. (a) 0 ms. (b) 0.5 ms. (c) 0.98 ms. (d) 1.3 ms. (e) 2.2 ms. (f) 5 ms. (g) 7 ms. (h) 9.3 ms. (i) 12 ms.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Note that the pressure inside the tunnel (0.1 MPa) is equal to the standard atmospheric pressure at 0 ms. The shock wave expands in the form of a spherical wave during the initial state of the explosion, but the overpressure decreases quickly and decays from 153.7 MPa to 12.9 MPa with an increase in the area of the wavefront, as shown in Figures7(a) and 7(c).During the interaction between the shock wave and the concrete lining, the incident wave reflects off the lining due to the difference in the acoustic impedance between the concrete and the air, leading to a significant increase in the strength of the shock wave. The strongest reflection effect occurs in the corner of the sidewall, where the pressure increases to 47.89 MPa, as seen in Figures7(d) and 7(e).During the explosion, a negative pressure region is formed in the vicinity of the detonation site as the air burns out. Meanwhile, the reflected shock wave propagates towards the detonation site and propagates along the longitudinal direction of the tunnel coincident with the original incident shock wave accompanied by a gradual decrease in the negative pressure region. These phenomena can be clearly observed in Figures7(f) and 7(g).Figures7(h) and 7(i) demonstrate the eventual convergence of the reflected shock wave along the longitudinal axis of the tunnel and the disappearance of the negative pressure region; after these events, the shock waves continue to propagate back and forth along their own trajectories.According to the above numerical findings, the pressure field is highly complicated, and the overpressure inside the tunnel remains high due to the arbitrary reflections of the shock waves, further indicating that the numerical method adopted in this paper can effectively reveal the details of the explosion event inside the tunnel. Meanwhile, it is very important to investigate the responses of the lining subjected to severe blast impact loading.


## 4.2. The Effective Stress Propagation in the Lining

Figure8 depicts the contours of the effective stress in the tunnel lining subjected to the explosion (the units in the diagrams are 102 GPa). The transition of the colours from blue to red represents an increase in the effective stress. The blast impact loading, which is transferred from the air in the form of a spherical shock wave as shown in Figure 8(a), simultaneously acts on the floor and vault of the lining at 0.98 ms. During the interaction between the blast shock wave and the lining, the stress propagates along the circumferential and longitudinal directions of the tunnel lining until the entire lining is subjected to a high blast impact loading. The concentration of stress primarily occurs at the corner of the sidewall due to the “horn structure,” which significantly strengthens the reflection effect and consequently increases the strength of the blast shock wave. Figure 8 also demonstrates that the stress acting on the lining does not decrease as the duration of the impact loading increases, attributed to the complicated propagation features of the blast shock wave inside the tunnel, as previously mentioned. Moreover, the lining in the vicinity of the detonation site suffers from considerable permanent deformation, which is evident in Figure 8(f). Hereafter, this paper will focus on the dynamic responses of the lining under high blast loading.Figure 8
Effective stress cloud charts of the tunnel lining. (a) 0.98 ms. (b) 1.3 ms. (c) 2.2 ms. (d) 3.6 ms. (e) 8.2 ms. (f) 15.6 ms.
(a)
(b)
(c)
(d)
(e)
(f)


## 4.3. The Dynamic Responses of the Lining

In this section, some measurement points are selected on the lining along the circumferential and longitudinal directions to investigate the dynamic responses of the lining when subjected to a blast impact loading. Figure9 illustrates the schematic arrangement of the measurement points. Along the circumferential direction (Figure 9(a)), points 1∼5 divide the arch of the lining into 4 equal parts (point 3 is located at the middle of the arch), points 5∼7 divide the sidewall of the lining into 2 equal parts (point 6 is located at the middle of the sidewall), and points 7∼9 divide the floor of the lining into 2 equal parts. Along the longitudinal direction, the points within 10 m of the detonation site are arranged at a 1 m interval, and the points at a distance greater than 10 m away from the explosion centre are arranged at an interval of 2 m, as shown in Figure 9(b). The dynamic responses of the lining are analysed in terms of the velocity, displacement, and acceleration of the measurement points, the numerical results of which are presented hereafter.Figure 9
Schematics of the measurement points on the lining. (a) Cross section view. (b) Longitudinal cross section view.
(a)
(b)Figure10 clearly reveals that the peak values of the three indexes of all of the measurement points are relatively high in the region within 5 m of the detonation site because the blast impact loading is the strongest immediately adjacent to the detonation site. The maximum values of the resultant velocity, displacement, and acceleration at point 1 (located at the vault of the arch) resulting from the blast shock wave that arrived first at the vault are 57.5 m/s, 39.1 cm, and 18482.7 g, respectively.Figure 10
Dynamic responses of the lining arch. (a) Velocity responses. (b) Displacement responses. (c) Acceleration responses.
(a)
(b)
(c)The trends of the velocity curves of the points are all consistent with those of the displacement curves, and the peak displacements of all of the measurement points are greater than 6.5 cm in the region within 5 m of the detonation site. This indicates that the arch of the lining in this region suffered significant deformation (i.e., outward expansion) upon being subjected to violent blast impact loading, as illustrated in Figure10(b). The peak velocities and displacements at point 5 are relatively higher than those at the other points at distances within 5∼14 m of the explosion centre, as shown in Figures 10(a) and 10(b). The peak accelerations at point 3 are larger than those at the other points at distances within 3∼10 m of the detonation, indicating that the middle part of the arch experienced a high blast impact loading, as seen in Figure 10(c), which was caused by the complicated reflection effect of the blast shock wave in this region. Due to the decrease in the blast loading with an increase in the distance from the detonation centre, all of the peak values of these factors decay rapidly; moreover, the general trends of the peak curves presented in Figure 10 tend to flatten.The variation curves of the dynamic responses of the sidewall that are clearly presented in Figure11 demonstrate that the lower part of the sidewall clearly exhibited a greater response than the upper part. The maximum peak velocity and acceleration at point 7 (Figure 9), where the “horn structure” at the corner of the sidewall significantly increases the reflection effect and consequently strengthens the blast shock wave, are 32.4 m/s and 8648 g, respectively. The maximum displacement at point 6 is 38.2 cm, which is 95.5% of the thickness of the structure, thereby demonstrating that the middle of the sidewall suffers a severe deflection, leading to complete failure.Figure 11
Dynamic responses of the lining sidewall. (a) Velocity responses. (b) Displacement responses. (c) Acceleration responses.
(a)
(b)
(c)The velocity and displacement curves of point 9 clearly oscillate, as shown in Figure12. The reasons for this oscillation may be twofold: (1) the blast impact loading first acts on point 9, which is the projected position of the detonation centre onto the tunnel floor, during the propagation of the blast shock wave identical to point 1(Figure 9), and (2) the tunnel drain channel is located beneath point 9, and this hollow structure complicates both the shock wave and the impact loading at the position of point 9. The maximum peak displacement at point 9 is 28.4 cm downward, which is marginally lower than that at point 7. The peak acceleration at point 8 is larger than those of the other points at distances within 10 m of the detonation site, indicating that the impact loading acting on this point is larger.Figure 12
Dynamic responses of the lining floor. (a) Velocity responses. (b) Displacement responses. (c) Acceleration responses.
(a)
(b)
(c)


## 4.4. The Damage to the Lining

### 4.4.1. The Damage Mechanism of the Lining

Cloud charts of the damage traits of the lining after being subjected to the blast impact loading are clearly presented in Figure13, in which a damage level of 1 indicates that the structure has been completely damaged, and a level of 0 denotes a fully intact structure. The lining in the vicinity of the detonation clearly experienced large-scale damage; in addition, the arch and sidewall expanded outward substantially, while the floor moved downward perceptibly. Far away from the detonation, the failure traits of the lining are primarily similar to those of the lining in the local damage area, particularly in the corner of the sidewall and the floor, and dense cracks are outspread over the arch. These characteristics are in accordance with the previously analysed dynamic responses of the lining.Figure 13
Overall degree of damage of the tunnel lining. (a) Overall damage situation. (b) Details of the cross section. (c) Details of the tunnel floor.
(a)
(b)
(c)Figure14 shows the relationship curves of the peak principal stress against the distance to the detonation centre at all measurement points presented in Figure 9. The positive and negative values of the peak principal stresses illustrated in Figure 14 represent the tensile and compressive stresses acting on the lining, respectively.Figure 14
Stress states of the lining.Figure14 demonstrates that the tensile stresses acting on the lining are high (approximately 5 MPa) and that the trends of the curves are stable (i.e., straight lines). However, the compressive stresses experienced by the lining show a wide range of changes; the values in the vicinity of the detonation site are the highest, and the stresses tend to gradually stabilize with an increase in the distance.The damage mechanism of the lining can be concluded as follows. (1) A local failure zone was formed in the lining at distances within 5.1 m from the detonation site when it was subjected to the significant stress exerted by the blast shock wave; then, under the alternating actions of the strongest compressive and tensile stresses, the failure zone rapidly expanded and eventually developed into a zone of complete failure. (2) In the region between 5.1 and 10.9 m from the explosion centre, the dense cracks within the lining propagated and became interconnected under higher compressive and tensile stresses, and the intact lining was divided into blocky structures, thereby leading to the loss of its overall stability and the formation of a severely damaged zone. (3) In the region beyond 10.9 m from the detonation site, a generally damaged zone was formed, in which only cracks formed throughout the lining.


### 4.4.2. The Damage Process of the Lining

To understand the entire damage process of the lining in detail, this paper utilizes the erosion algorithm for components in LS-DYNA; this algorithm has the ability to simulate the physical failure process of the lining by deleting the components that satisfy specific erosion criteria. The peak tensile strain of concrete is approximately 2E-4 under quasistatic loading, while the strain rate and the corresponding strength amplification factor can exceed 10∼100 s−1 and 5.0 under blast impact loading, respectively. However, the enhancement factor of the failure strain is less than that of the strength factor; thus, in this section, a tensile failure strain of 1E-3 is adopted for the concrete as the eroding plastic strain, meaning that a component of the lining fails and is deleted upon reaching this threshold value [37].Figure15 illustrates the numerical simulation of the full damage process of the lining and provides comparisons with the results of field investigations. A locally damaged zone first appeared in the arch and the middle of the lining floor; then, a failure zone formed in the corner of the sidewall due to the “horn structure” that significantly strengthened the blast impact loading. Hence, the locally damaged zone expanded along the longitudinal and circumferential directions of the lining before finally developing into a completely damaged zone, as clearly shown in Figures 15(a)–15(d). The strength of the blast loading decreased during the propagation of the shock wave, but the complete failure of the lining in the vicinity of the detonation aggravated the progressive failure of the lining far away from the detonation, resulting in the formation of a severely damaged zone between 5.1 and 10.9 m from the detonation site, at which distances the intact lining was cut into blocks of various sizes as demonstrated in Figures 15(d)–15(e). For the zone situated more than 10.9 m from the explosion centre, cracks were produced that propagated across the arch under high tensile stresses, as shown in Figures 15(e)–15(f); hence, a generally damaged zone developed. Comparisons of the damage features between the numerical and field investigations are illustrated in Figure 15(g), from which it is evident that the results of the numerical simulation agree well with the actual situation and that the research methods used in this paper are both feasible and effective.Figure 15
Simulations of the damage and comparison with the actual lining situation. (a) 1.5 ms. (b) 2 ms. (c) 2.5 ms. (d) 4 ms. (e) 8 ms. (f) 20 ms. (g) Damage details of the numerical and site investigations.
(a)
(b)
(c)
(d)
(e)
(f)
(g)


## 4.4.1. The Damage Mechanism of the Lining

Cloud charts of the damage traits of the lining after being subjected to the blast impact loading are clearly presented in Figure13, in which a damage level of 1 indicates that the structure has been completely damaged, and a level of 0 denotes a fully intact structure. The lining in the vicinity of the detonation clearly experienced large-scale damage; in addition, the arch and sidewall expanded outward substantially, while the floor moved downward perceptibly. Far away from the detonation, the failure traits of the lining are primarily similar to those of the lining in the local damage area, particularly in the corner of the sidewall and the floor, and dense cracks are outspread over the arch. These characteristics are in accordance with the previously analysed dynamic responses of the lining.Figure 13
Overall degree of damage of the tunnel lining. (a) Overall damage situation. (b) Details of the cross section. (c) Details of the tunnel floor.
(a)
(b)
(c)Figure14 shows the relationship curves of the peak principal stress against the distance to the detonation centre at all measurement points presented in Figure 9. The positive and negative values of the peak principal stresses illustrated in Figure 14 represent the tensile and compressive stresses acting on the lining, respectively.Figure 14
Stress states of the lining.Figure14 demonstrates that the tensile stresses acting on the lining are high (approximately 5 MPa) and that the trends of the curves are stable (i.e., straight lines). However, the compressive stresses experienced by the lining show a wide range of changes; the values in the vicinity of the detonation site are the highest, and the stresses tend to gradually stabilize with an increase in the distance.The damage mechanism of the lining can be concluded as follows. (1) A local failure zone was formed in the lining at distances within 5.1 m from the detonation site when it was subjected to the significant stress exerted by the blast shock wave; then, under the alternating actions of the strongest compressive and tensile stresses, the failure zone rapidly expanded and eventually developed into a zone of complete failure. (2) In the region between 5.1 and 10.9 m from the explosion centre, the dense cracks within the lining propagated and became interconnected under higher compressive and tensile stresses, and the intact lining was divided into blocky structures, thereby leading to the loss of its overall stability and the formation of a severely damaged zone. (3) In the region beyond 10.9 m from the detonation site, a generally damaged zone was formed, in which only cracks formed throughout the lining.


## 4.4.2. The Damage Process of the Lining

To understand the entire damage process of the lining in detail, this paper utilizes the erosion algorithm for components in LS-DYNA; this algorithm has the ability to simulate the physical failure process of the lining by deleting the components that satisfy specific erosion criteria. The peak tensile strain of concrete is approximately 2E-4 under quasistatic loading, while the strain rate and the corresponding strength amplification factor can exceed 10∼100 s−1 and 5.0 under blast impact loading, respectively. However, the enhancement factor of the failure strain is less than that of the strength factor; thus, in this section, a tensile failure strain of 1E-3 is adopted for the concrete as the eroding plastic strain, meaning that a component of the lining fails and is deleted upon reaching this threshold value [37].Figure15 illustrates the numerical simulation of the full damage process of the lining and provides comparisons with the results of field investigations. A locally damaged zone first appeared in the arch and the middle of the lining floor; then, a failure zone formed in the corner of the sidewall due to the “horn structure” that significantly strengthened the blast impact loading. Hence, the locally damaged zone expanded along the longitudinal and circumferential directions of the lining before finally developing into a completely damaged zone, as clearly shown in Figures 15(a)–15(d). The strength of the blast loading decreased during the propagation of the shock wave, but the complete failure of the lining in the vicinity of the detonation aggravated the progressive failure of the lining far away from the detonation, resulting in the formation of a severely damaged zone between 5.1 and 10.9 m from the detonation site, at which distances the intact lining was cut into blocks of various sizes as demonstrated in Figures 15(d)–15(e). For the zone situated more than 10.9 m from the explosion centre, cracks were produced that propagated across the arch under high tensile stresses, as shown in Figures 15(e)–15(f); hence, a generally damaged zone developed. Comparisons of the damage features between the numerical and field investigations are illustrated in Figure 15(g), from which it is evident that the results of the numerical simulation agree well with the actual situation and that the research methods used in this paper are both feasible and effective.Figure 15
Simulations of the damage and comparison with the actual lining situation. (a) 1.5 ms. (b) 2 ms. (c) 2.5 ms. (d) 4 ms. (e) 8 ms. (f) 20 ms. (g) Damage details of the numerical and site investigations.
(a)
(b)
(c)
(d)
(e)
(f)
(g)


## 4.5. Discussion

### 4.5.1. The Blast Impact Loading

In this paper, the dynamic responses and damage characteristics of the lining are directly related to the mode of the blast impact loading. The abovementioned research findings are based on the FSI numerical model. The keys of this model are that the modelled explosive charge forms the blast shock wave and that the lining is subjected to a fully coupled blast loading. Figure16 depicts both the pressure-time curves of the shock wave that propagates through the air and strikes an obstacle and the simplified triangle impact loading curve. Both types of curves account for the enhanced strength of the shock wave due to reflection effects, and the reflected pressure exhibits a certain relationship with the incident pressure in the diagram. The simplified loading curve can conserve considerably more computational resources than the typical reflected curve and has thus been adopted by many researchers [38–40].Figure 16
P-T curves of the blast loading.The P-T curves, which are described by the fully coupled FSI model, of selected measurement points (Figure17(a)) on the lining are presented in Figures 17(b) and 17(c).Figure 17
P-T curves of different positions on the lining described by the numerical simulation. (a) Locations of the measurement points on the lining. (b) P-T curves of the measurement points on the same cross section. (c) P-T curves of the measurement points at the same positions.
(a)
(b)
(c)The curves in Figures16 and 17 are dissimilar; furthermore, the P-T curves of the points at various locations on the lining shown in Figure 17 are also different both with regard to the trends of the curves and the peak pressures. All of the curves clearly experience multiple peak values during 0∼4 ms, which can be explained by the reflection effect of the blast shock wave. Moreover, lines A and D and lines 3 and 4 also exhibit several peak pressures after 6 ms. The reflected peak pressures of the measurement points are hardly describable using a given relationship with the incident pressure; one of the main reasons for this might be the variations in the shock wave field caused by the reflections at different locations. The fully coupled blast loading model can represent all of the reflection effects at various locations as the wave propagates inside the tunnel, while the simplified loading model cannot. The wave transformed from a compressive wave to a tensile wave; under the tensile stress, cracks emerged throughout the lining, and the damaged lining broke into blocks of concrete, forming zones of severe and general damage. This result is in accordance with the findings of Section 4.4. Based on the analysis in this section, use of the fully coupled blast loading model has obvious advantages for the research object of this paper.


### 4.5.2. The Damage Assessment of the Lining

The single degree of freedom (SDOF) model proposed by Fallah and Louca [41], which simplifies the structure into an equivalent perfect elastic-plastic model, is applied to quantitatively study and assess the damage of all types of structures, especially for the column structure [42, 43]. The SDOF method is used to assess the damage levels of an underground tunnel when undergoing several explosion cases [22]. In this section, this effective method is used for reference to estimate the damage to the lining structure. The criterion is based on the maximum deflection at the mid-height of the structure, and different deflections correspond to different damage levels. The critical values of the deflections differentiating various damage levels contained in the SDOF model are given in Table 6.Table 6
Damage assessment criterion of the SDOF approach.
Damage level
Maximum displacement (cm)
Low damage (LD)
<2
Moderate damage (MD)
2∼4
High damage (HD)
4∼8
Full damage (FD)
>8According to the SDOF approach, the positions of points 3, 6, and 9 in Figure9(a) represent the arch, the sidewall, and the floor of the lining, respectively, as shown in Figure 18(a). The maximum deflection curves of these measured points are shown in Figure 18(b).Figure 18
Evaluation of the lining damage based on the SDOF method. (a) Arrangement of the measurement points. (b) Maximum deflection curves.
(a)
(b)Figure18(b) reveals that the deflection curves of the arch and sidewall show decaying trends as the distance to the detonation site increases, while the deflection curve for the floor exhibits relative fluctuations related to the propagation of the blast shock wave and the previously mentioned geometric construction of the tunnel floor. Table 7 describes the detailed damage characteristics for the arch, sidewall, and floor of the lining. The maximum deflections exceed 8 cm in the arch and sidewall corresponding to full damage (FD) in the region within 5 m of the detonation site, while the degree of damage ranges from FD to moderate damage (MD) in the region between 6 and 10 m from the explosion centre, and low damage (LD) occurs when the distance is greater than 10 m. Meanwhile, the general trend of the floor curve is attenuated. In summary, it is clearly evident that the degrees of the damage of the lining evaluated using the SDOF method agree well with the site conditions and the numerical results.Table 7
Damage characteristics of the lining.
Damage degree
Parts of the lining
Distance to the detonation (m)
0
1
2
3
4
5
6
7
8
9
10
12
14
16
18
20
Arch
FD
FD
FD
FD
FD
FD
HD
HD
MD
MD
MD
LD
LD
LD
LD
LD
Sidewall
FD
FD
FD
FD
FD
FD
FD
HD
HD
HD
MD
LD
LD
LD
LD
LD
Floor
FD
FD
HD
MD
HD
HD
FD
FD
FD
HD
MD
LD
MD
MD
LD
LD


## 4.5.1. The Blast Impact Loading

In this paper, the dynamic responses and damage characteristics of the lining are directly related to the mode of the blast impact loading. The abovementioned research findings are based on the FSI numerical model. The keys of this model are that the modelled explosive charge forms the blast shock wave and that the lining is subjected to a fully coupled blast loading. Figure16 depicts both the pressure-time curves of the shock wave that propagates through the air and strikes an obstacle and the simplified triangle impact loading curve. Both types of curves account for the enhanced strength of the shock wave due to reflection effects, and the reflected pressure exhibits a certain relationship with the incident pressure in the diagram. The simplified loading curve can conserve considerably more computational resources than the typical reflected curve and has thus been adopted by many researchers [38–40].Figure 16
P-T curves of the blast loading.The P-T curves, which are described by the fully coupled FSI model, of selected measurement points (Figure17(a)) on the lining are presented in Figures 17(b) and 17(c).Figure 17
P-T curves of different positions on the lining described by the numerical simulation. (a) Locations of the measurement points on the lining. (b) P-T curves of the measurement points on the same cross section. (c) P-T curves of the measurement points at the same positions.
(a)
(b)
(c)The curves in Figures16 and 17 are dissimilar; furthermore, the P-T curves of the points at various locations on the lining shown in Figure 17 are also different both with regard to the trends of the curves and the peak pressures. All of the curves clearly experience multiple peak values during 0∼4 ms, which can be explained by the reflection effect of the blast shock wave. Moreover, lines A and D and lines 3 and 4 also exhibit several peak pressures after 6 ms. The reflected peak pressures of the measurement points are hardly describable using a given relationship with the incident pressure; one of the main reasons for this might be the variations in the shock wave field caused by the reflections at different locations. The fully coupled blast loading model can represent all of the reflection effects at various locations as the wave propagates inside the tunnel, while the simplified loading model cannot. The wave transformed from a compressive wave to a tensile wave; under the tensile stress, cracks emerged throughout the lining, and the damaged lining broke into blocks of concrete, forming zones of severe and general damage. This result is in accordance with the findings of Section 4.4. Based on the analysis in this section, use of the fully coupled blast loading model has obvious advantages for the research object of this paper.


## 4.5.2. The Damage Assessment of the Lining

The single degree of freedom (SDOF) model proposed by Fallah and Louca [41], which simplifies the structure into an equivalent perfect elastic-plastic model, is applied to quantitatively study and assess the damage of all types of structures, especially for the column structure [42, 43]. The SDOF method is used to assess the damage levels of an underground tunnel when undergoing several explosion cases [22]. In this section, this effective method is used for reference to estimate the damage to the lining structure. The criterion is based on the maximum deflection at the mid-height of the structure, and different deflections correspond to different damage levels. The critical values of the deflections differentiating various damage levels contained in the SDOF model are given in Table 6.Table 6
Damage assessment criterion of the SDOF approach.
Damage level
Maximum displacement (cm)
Low damage (LD)
<2
Moderate damage (MD)
2∼4
High damage (HD)
4∼8
Full damage (FD)
>8According to the SDOF approach, the positions of points 3, 6, and 9 in Figure9(a) represent the arch, the sidewall, and the floor of the lining, respectively, as shown in Figure 18(a). The maximum deflection curves of these measured points are shown in Figure 18(b).Figure 18
Evaluation of the lining damage based on the SDOF method. (a) Arrangement of the measurement points. (b) Maximum deflection curves.
(a)
(b)Figure18(b) reveals that the deflection curves of the arch and sidewall show decaying trends as the distance to the detonation site increases, while the deflection curve for the floor exhibits relative fluctuations related to the propagation of the blast shock wave and the previously mentioned geometric construction of the tunnel floor. Table 7 describes the detailed damage characteristics for the arch, sidewall, and floor of the lining. The maximum deflections exceed 8 cm in the arch and sidewall corresponding to full damage (FD) in the region within 5 m of the detonation site, while the degree of damage ranges from FD to moderate damage (MD) in the region between 6 and 10 m from the explosion centre, and low damage (LD) occurs when the distance is greater than 10 m. Meanwhile, the general trend of the floor curve is attenuated. In summary, it is clearly evident that the degrees of the damage of the lining evaluated using the SDOF method agree well with the site conditions and the numerical results.Table 7
Damage characteristics of the lining.
Damage degree
Parts of the lining
Distance to the detonation (m)
0
1
2
3
4
5
6
7
8
9
10
12
14
16
18
20
Arch
FD
FD
FD
FD
FD
FD
HD
HD
MD
MD
MD
LD
LD
LD
LD
LD
Sidewall
FD
FD
FD
FD
FD
FD
FD
HD
HD
HD
MD
LD
LD
LD
LD
LD
Floor
FD
FD
HD
MD
HD
HD
FD
FD
FD
HD
MD
LD
MD
MD
LD
LD


## 5. Conclusions

(1)
This paper reproduced the detailed propagation process of the blast shock wave inside the Luodaiguzhen tunnel. The wave expanded in the form of a spherical wave before propagating to the lining, and the energy attenuated quickly. The reflection of the incident wave from the lining significantly enhanced the strength of the wave. The strongest reflection effect was located in the corner of the sidewall because of its “horn structure”. Consequently, the internal pressure within the tunnel was incredibly high, and the wave field was highly complicated due to the infinite number of wave reflections.
(2)
The blast impact loading first acted on the arch vault and floor of the lining and then propagated along the radial and longitudinal directions until the entire lining was subjected to the blast impact loading. The corner of the sidewall evidently exhibited a concentration of stress. Moreover, the lining in the vicinity of the explosive centre incurred considerable outward expansion and deformation.
(3)
The maximum peak values of the velocity, displacement, and acceleration at the arch vault of the lining reached 57.5 m/s, 39.1 cm, and 18482.7 g, respectively. The corresponding values at the corner of the sidewall were 32.4 m/s, 29.6 cm, and 8648 g. These two locations on the lining were clearly affected the most by the blast impact loading.
(4)
The analysis of the peak principal stresses in the lining showed that the lining was subjected to high tensile stresses (approximately 5 MPa) during the explosion, while the compressive stress varied from 10 MPa to 100 MPa. Combined with the damage results calculated during the numerical simulation, the region of the lining within 5.1 m of the detonation site was completely destroyed and consequently collapsed under significant stresses. Meanwhile, the region of the lining between 5.1 and 10.9 m from the explosion centre was fragmented into blocks and lost its overall stability under higher compressive and tensile stresses, thereby developing into a zone of severe damage. Finally, the region of the lining beyond 10.9 m from the detonation was generally damaged by tensile stresses.
(5)
The damage process of the lining was reproduced through the erosion algorithm for the components of the lining. In comparison, the simulated damage traits of the lining agreed well with the site conditions. Moreover, the fully coupled blast loading model was superior to the simplified blast loading model when describing the real responses and damage of the lining. The SDOF method was used to estimate the degrees of damage suffered by the lining, and the results were consistent with the actual conditions. Meanwhile, the numerical simulations conducted in this paper were validated, and they are thus reasonable and feasible for the study of large structures subjected to violent internal explosions.



---

*Source: 2792043-2018-10-08.xml*