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

*e 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. *e 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. *en, a detailed modelling process is presented and validated through a comparison with empirical formulas. *is 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. *e 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. *e 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.


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. Table 1 and Figure 1 describe several typical explosion accidents that occurred inside tunnels throughout China, thereby highlighting the importance of further researching such incidents. Some fruitful studies have been conducted on gas explosions and underground structures subjected to blast loading. For example, some scholars have performed scaleddown centrifuge modelling tests to describe the responses of tunnels under blast loading [1][2][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][5][6][7][8][9][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. is technology is also capable of displaying the details of the entire explosion process at each time step. Jayasinghe et al. [12][13][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][17][18][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. is paper presents a detailed background description of a tunnel project to quantitatively study gas explosions. Subsequently, a full-scale fluidstructure 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. e 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. e objective of this research, which is highly innovative and has scientific value, is to propose a feasible and effective approach for similar explosion events.

Project Background
e gas explosion that occurred inside the Luodaiguzhen tunnel is the only typical, recent example of an engineering case of such an incident worldwide. e tunnel, which is  2 Shock and Vibration 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. e 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 m 3 /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. is 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.
Figures 2(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). e region of construction and residential buildings near the tunnel entrance were completely destroyed, as shown in Figure 2(d).
e above descriptions suggest that the magnitude of the gas explosion within the Luodaiguzhen tunnel was tremendous, and the linings were badly destroyed. e 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 Figure 3. e 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. is 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.

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-todiameter 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. ese 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]. e 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. is 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 (CH 4 ) is the main combustible component of gas. Accordingly, the TNT-equivalence method is given by where E q is the ratio of the combustion heat of methane to that of TNT; Q G is the combustion heat of a unit mass of methane; J/kg; Q TNT is the combustion heat of a unit mass of TNT, 4.5 MJ/kg; M TNT is the equivalent weight of TNT; α is the methane concentration in the gas; V G is the volume of the accumulated gas, m 3 ; and ρ G is the density of methane, 0.716 kg/m 3 . is paper considers the most unfavourable situation by assuming that all methane has combusted completely, resulting in a maximum explosive energy. e combustion heat of methane is 55.64 MJ/kg according to the following expression [24]:    Shock and Vibration Consequently, the concentration of methane is 9.5%: where n 0 , 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. e amount of gas over ow detected in the ZK2 + 790∼ZK2 + 830 region during the construction of the tunnel was 1.08 m /min. Furthermore, the crosssectional area of the tunnel was 72.34 m 2 . 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).

e Numerical Model.
e 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. Figures 4(a) and 4(b) show the design and cross section profile drawings of the standard tunnel cross section and the numerical model, respectively. e units of r 1 to r 5 in Figure 4(a) are centimeters. e 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.
Considering the symmetry of the numerical model demonstrated in Figure 4(c), a quarter of the model is established to optimally utilize the computational resources, as shown in Figure 4(d).
e 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. e 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. Figure 4(d) depicts the five components of the FSI finite element model: the TNT source, air, first lining, secondary lining, and surrounding rock. e 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. e ALE multimaterial mesh model is achieved by using * CONSTRAINED_LAGRANGE_IN_SOLID in LS-DYNA [30]. e 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. is 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. rough refined meshing skills, the elements used to compose all of the components are eight-node solid hexahedron elements, which can improve the calculation efficiency. e 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.

TNT.
e 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: 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, R 1 , R 2 , and ω are empirical constants based on explosion experiments. e calculation parameters are listed in Table 2.

3.3.2.
Air. e air pressure, which is modelled using * MAT_NULL in LS-DYNA with a linear polynomial * EOS_LINEAR_POLYNOMIAL [30], is given by the following: where P is the air pressure; C 0 ∼C 6 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. e calculation parameters are listed in Table 3.

e Lining Structures.
e Riedel-Hiermaier-oma (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]. e RHT model takes ratedependent 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: where f c is the quasistatic uniaxial compressive strength of concrete, 25 MPa; p s � p/F rate (_ ε) is the quasistatic pressure; σ * TXC (p s ) is the equivalent stress strength of the quasistatic compressive meridian; R 3 (θ) is a scalar function of the load angle θ and the tensile-to-compressive meridian ratio; F rate (_ ε) is an influencing factor on the strain rate; F elastic is the scaling factor of the elastic deformation; F cap is the cap function; and A f and N f are the residual surface parameters. e compressive and tensile strain rate-dependent exponents β c and β t are 0.042 and 0.044, respectively, via the following [31]: 6 Shock and Vibration e accumulative damage factor D in the RHT model is defined as follows: where Δε p is the equivalent plastic strain increment; Δμ p is the plastic volumetric strain increment; D 1 and D 2 are the damage parameters with suggested values of 0.015 and 1.0, respectively; and ε min f is 0.0008 [34]. e 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. rough the following conversion formula, the equivalent tensile strength of the first lining R t 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: where f t is the tensile strength of concrete, 2.5 MPa; f y is the yield strength of steel, 300 MPa; and ρ is the steel content in the concrete, 1.2%.

e Surrounding
Rock. e 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. e constitutive relation of the surrounding rock can be described as follows: 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; E p is the plasticity hardening modulus; and E tan is the tangent modulus. e calculation parameters for the surrounding rock are shown in Table 5.

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: where ΔP f 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). e 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 Figure 5(a). Figure 5(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. ese findings also directly validate the modelling techniques

e Strength and Propagation Characteristics of the Blast
Shock Wave inside the Tunnel. Figure 6(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. e following can be discerned from Figure 6(b): (1) the strength of the shock wave attenuates signi cantly 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 di erent; (3) the radial curve exhibits clear uctuations, especially with regard to the abrupt, signi cant 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. Figure 7 illustrates the detailed propagation process of the shock wave inside the tunnel (the units in the diagrams are 10 2 GPa).      During the interaction between the shock wave and the concrete lining, the incident wave re ects o the lining due to the di erence in the acoustic impedance between the concrete and the air, leading to a signi cant increase in the strength of the shock wave. e strongest re ection e ect occurs in the corner of the sidewall, where the pressure increases to 47.89 MPa, as seen in Figures 7(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 re ected 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. ese phenomena can be clearly observed in Figures 7(f ) and 7(g).
Figures 7(h) and 7(i) demonstrate the eventual convergence of the re ected 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 ndings, the pressure eld is highly complicated, and the overpressure inside the tunnel remains high due to the arbitrary re ections of the shock waves, further indicating that the numerical method adopted in this paper can e ectively 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.
e E ective Stress Propagation in the Lining. Figure 8 depicts the contours of the e ective stress in the tunnel lining subjected to the explosion (the units in the diagrams are 10 2 GPa). e transition of the colours from blue to red represents an increase in the e ective stress. e 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 oor 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. e concentration of stress primarily occurs at the corner of the sidewall due to the "horn structure," which signi cantly strengthens the re ection e ect 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 su ers 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.

e 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. Figure 9 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 oor 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). e 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 10 clearly reveals that the peak values of the three indexes of all of the measurement points are relatively high in  e 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. is indicates that the arch of the lining in this region su ered signi cant deformation (i.e., outward expansion) upon being subjected to violent blast impact loading, as illustrated in Figure 10(b). e 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). e 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 re ection e ect 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 atten.     e variation curves of the dynamic responses of the sidewall that are clearly presented in Figure 11 demonstrate that the lower part of the sidewall clearly exhibited a greater response than the upper part. e 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. e 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. e velocity and displacement curves of point 9 clearly oscillate, as shown in Figure 12. e 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. e maximum peak displacement at point 9 is 28.4 cm downward, which is marginally lower than that at point 7. e 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.

e 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 Figure 13, 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. e 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. ese characteristics are in accordance with the previously analysed dynamic responses of the lining. Figure 14 shows the relationship curves of the peak principal stress against the distance to the detonation centre at all measurement points presented in Figure 9. e 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 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. e 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.

e 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. e 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]. Figure 15 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  e 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. Figure 16 depicts both the pressure-time curves of the shock wave that propagates through the air and strikes an obstacle and the simpli ed triangle impact loading curve. Both types of curves account for the enhanced strength of the shock wave due to re ection e ects, and the re ected pressure exhibits a certain relationship with the incident pressure in the diagram. e simpli ed loading curve can conserve considerably more computational resources than the typical re ected curve and has thus been adopted by many researchers [38][39][40].
e P-T curves, which are described by the fully coupled FSI model, of selected measurement points (Figure 17(a)) on the lining are presented in Figures 17(b) and 17(c). e curves in Figures 16 and 17 are dissimilar; furthermore, the P-T curves of the points at various locations on the lining shown in Figure 17 are also di erent 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 re ection e ect of the blast shock wave. Moreover, lines A and D and lines 3 and 4 also exhibit several peak pressures after 6 ms. e re ected 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 eld caused by the reections at di erent locations. e fully coupled blast loading model can represent all of the re ection e ects at various locations as the wave propagates inside the tunnel, while the simpli ed loading model cannot. e 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. is result is in accordance with the ndings 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.
e Damage Assessment of the Lining. e single degree of freedom (SDOF) model proposed by Fallah and Louca [41], which simpli es 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]. e SDOF method is used to assess the damage levels of an underground tunnel when undergoing several explosion cases [22]. In this section, this e ective method is used for reference to estimate the damage to the lining structure. e criterion is based on the maximum de ection at the mid-height of the structure, and di erent de ections correspond to di erent damage levels. e critical values of the de ections di erentiating various damage levels contained in the SDOF model are given in Table 6.
According to the SDOF approach, the positions of points 3, 6, and 9 in Figure 9(a) represent the arch, the sidewall, and the oor of the lining, respectively, as shown in Figure 18(a). e maximum de ection curves of these measured points are shown in Figure 18(b). Figure 18(b) reveals that the de ection curves of the arch and sidewall show decaying trends as the distance to the detonation site increases, while the de ection curve for the oor exhibits relative uctuations related to the propagation of the blast shock wave and the previously mentioned geometric construction of the tunnel oor. Table 7 describes the detailed damage characteristics for the arch, sidewall, and oor of the lining. e maximum de ections 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 oor 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.

Conclusions
(1) is paper reproduced the detailed propagation process of the blast shock wave inside the Luodaiguzhen tunnel. e wave expanded in the form of a spherical wave before propagating to the lining, and the energy attenuated quickly. e re ection of the incident wave from the lining signi cantly enhanced the strength of the wave. e strongest re ection e ect 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 eld was highly complicated due to the in nite number of wave re ections.     The typical incident wave P-T curve The typical reflected wave P-T curve The simplified triangle loading P-T curve 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) e 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 simpli ed blast loading model when describing the real responses and damage of the lining. e SDOF method was used to estimate the degrees of damage su ered 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.  Data Availability e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.   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   Shock and Vibration  19