Propagation of shock wave and structure dynamic response of explosion in a subway station: a case study of Wuhan subway station

2.1. Numerical modeling

A numerical model is established based on the Mafangshan station of Wuhan metro line 8 project. Mafangshan station is a typical two-layer three-span island subway station with reinforced concrete box frame structure. The width of the standard section of the station is 23.3 m, the length is about 120 m, and the total construction area of the station is approximately 8196 m², the overlying soil thickness of station is about 3.0 m. The dimensions of the main structure of the station are shown in Fig. 2. The size of the pillar in Fig. 2 is 800 mm×1100 mm, and the distance between pillars is 7.2 m, both vertically and horizontally.

Fig. 2The cross-section of Mafangshan station (unit: mm)

The numerical model of Mafangshan station is established by ANSYS/ LS-DYNA. The established finite element model adopts the cm-g-μs unit system. Due to the location of accidental explosions in subway stations is unpredictable. This paper is unable to conduct numerical simulation research on setting explosives in every place of subway station. In order to facilitate modeling calculation and research analysis, the explosion point is selected in the center of the station platform layer, and the height of the platform is as 40 cm. The explosive source is TNT square explosive, and different explosive charges of 5 kg, 10 kg, 20 kg, 30 kg, and 40 kg are considered to analyze the dynamic responses of the structure under different explosion loadings. The actual construction of the station is simplified in the numerical model, according to the symmetry of the station, a quarterly model of the subway station is built. The depth of the station is 3 m, the sidewall and lower part of the station are also made with a 3 m thick soil layer. The rest of the dimensions are determined by the actual size of the subway station. Finally, we obtained a numerical model of 6100×2180×1465 (cm), as shown in Fig. 3(a). The $X$-axis is the longitudinal direction of the station, the $Y$-axis is the vertical direction, and the $Z$-axis is the horizontal direction. The gap in the middle floor is the station stairs and elevator shaft. The finite element model adopts Soild164 element. Material models include explosive, air, reinforced concrete and soil. Explosive, air and reinforced concrete are modeled using Euler-mesh, and the soil layer is modeled using Lagrange-mesh, as shown in Fig. 3(b). Fluid structure interaction (FSI) algorithm is adopted to calculate the concrete, explosive and the air due to the interaction of them [17].

Fig. 3The finite element model of the subway station

a)

b)

The soil around the station is a semi-infinite space. Due to the limitation of modeling, it is unable to build the whole soil around the station. Therefore, an only limited model with station structure and its surrounding soil is established. However, there will be reflection superposition of the stress waves at the model boundary without boundary conditions treatments, which will lead to imprecise results of the simulation. Lysmer proposed a non-reflecting boundary algorithm [18]. When non-reflecting boundary algorithm is set on the boundary, the boundary can be simulated as an infinite field. As the wave propagates to the boundary, it can be transmitted to an infinite distance. It is achieved by applying shear and viscous normal stresses on the boundary, as shown in the following formula:

where $S_n$ is normal stress, $S_{\tau }$ is Tangential stress, $\alpha$ is the damping coefficient, $\rho$ is medium density, $C_n$ and $C_{\tau }$ are the velocities of longitudinal wave and transverse wave, respectively, $V_n$ and $V_{\tau }$ are particles on the boundary of the normal velocity and tangential velocity, respectively.

In this paper, the negative $Z$ direction boundary surface, the negative $Y$ direction boundary surface of the model and the right boundary of the model in the $X$ direction are set with the non-reflective boundary algorithm. In addition, Both the negative $X$ direction boundary surface and the positive $Z$ direction boundary surface impose symmetry constraints, the positive $Y$ direction boundary surface is free, as shown in Fig. 3(b).

2.2. The material model for explosive

In this paper, the high-energy explosive material model *MAT_HIGH_EXPLOSIVE_BURN [19] and the state equation of Jones-Wilkins-Lee [20] are used to simulate the TNT explosive. The high-energy explosive material model used to specify the density, detonation velocity and detonation pressure of explosive. The JWL equation of state is used to simulate the relationship between specific volume and pressure during explosive detonation. The equation defines the pressure as:

where $P$ is pressure, $A$, $B$, $R_1$, $R_2$, and $\omega$ are material constants, $V$ is the relative volume, $E_0$ is the initial specific internal energy.

The explosive, which is used in the numerical simulation, is TNT explosive, which relevant parameters refer to relevant literature [21]. The explosive density is 1630 kg/m³, the detonation pressure is 21 GPa, and the detonation velocity is 6930 m/s. The parameters of the explosive used in the Eq. (1) are shown in Table 1.

2.3. The material model for air

The air in station is defined as an ideal gas, which is simulated by *MAT_NULL material model and the LINEAR_POLYNOMIA equation of state. The equation of state defines the pressure as:

where $C_0$-$C_6$ is the material constant, $E_0$ is the specific internal energy, $\mu$ is the specific energy.

The air density is taken as 1.29 kg/m³. The state equation of ideal air is referred to the literature [22], and its parameters are shown in Table 2.

2.4. The material model for soil

The soil around the subway station is simulated by *MAT_DRUCKER_PRAGER material model [23], which is used to define the parameters of the yield surface as the friction angle $\phi$ and cohesion $c$ of the soil. This material model uses the improved Drucker Prager yield criterion, so that the shape of the yield surface can be distorted into a more realistic soil model. The expression of yield surface is as follows:

where $T$ is the shear stress intensity, ${\sigma }_m$ is the average stress, $\beta$, ${\sigma }_y$ are the model parameters, which can be determined by the friction angle $\phi$ and cohesion $c$ of soil:

where $\phi$ is the friction angle, $c$ is the cohesion.

The main parameters of the soil model around the station are shown in Table 3.

2.5. The material model for reinforced concrete

The primary structural material of the station is reinforced concrete. Reinforced concrete consists of concrete and reinforcement, if the two materials are treated as different elements, the calculation results are accurate, but the modeling is complex and challenging to achieve. Therefore, we through the method of equal modulus to define the reinforced concrete as an equivalent concrete material, which approximately simulates the reinforced concrete structure [24]. The material adopts *MAT_PLASTIC_KINEMATIC material model, which is suitable for isotropic plastic follow-up strengthening material including strain rate effect. This model can be used to simulate the stress-strain behavior of concrete material. The yield condition is as follows:

where ${\xi }_{ij}$, ${\sigma }_y$can be determined by the following formula:

where $S_{ij}$ is Cauchy stress tensor, $p$ and $c$ are the input constant, $\beta$ is hardening parameter, ${\sigma }_0$ is the yield stress, $\epsilon$ and ${\epsilon }_{eff}^p$ are strain rate and effective plastic strain respectively, $E_P$ is the plastic hardening modulus.

The parameters of reinforced concrete model are shown in Table 4.

4.1. Propagation of shock wave

Different from the propagation law of shock wave in infinite space, the propagation law of shock wave in subway station is complicated due to the characteristics of its semi-closed structure [29]. In this paper, the shock wave propagation process of the 20 kg explosive model is analyzed by the post-processing software LS-PrePost, as shown in Fig. 6.

The detonation source is set at a low height, after reaching the platform layer, the shock wave reflected from the platform layer will overlap with the original shock wave before reflection. Therefore, it can be seen from Fig. 6 that the initial propagation of the shock wave is in the shape of a hemisphere. Then, at about 6ms, the shock wave propagates to the roof of the platform floor and reflection will occur after the contact with floor. It can be found that the superimposed pressure of the reflected wave and the incident wave is much higher than the pressure of the incident wave.

As the shock wave continues to propagate, it will spread to the track travel area and the sidewall. As the propagation distance increases, the energy of the shock wave will gradually decrease, and the pressure value at the sidewall will obviously decrease. However, due to the small space in the lower part of the track area, multiple reflections will occur between the track and the platform structure, so the pressure in this area will increase significantly. In addition, there is a bulge on the wavefront at the head of the stairs, indicating that the head of the stairs plays a role of explosion relief, transmitting part of the pressure to the station hall floor.

Then, at about 90 ms, from the overall overhead view of the station shown in Fig. 7, it can be found that the wavefront of the shock wave gradually becomes a plane from the original spherical or irregular convex surface and propagates to both ends of the station.

Fig. 6Nephogram of shock wave at different times (charge of TNT: 20 kg)

a)$t=$1.8 ms

b)$t=$6.0 ms

c)$t=$15.0 ms

d)$t=$32.8 ms

e)$t=$43.6 ms

f)$t=$70.4 ms

4.2. Structural displacement response of explosion

The numerical model with the explosive charge of 20 kg is shown and analyzed in LS-PrePost, and it showed that the plate, pillar, wall of the station had an absolute displacement after the explosive explosion. However, the charge of explosives is not very large, the range of structural displacement response is limited, and the impact of shock wave on the structure becomes relatively small after the fourth pillar. Therefore, only the structure within a specific range of the explosion is analyzed. Considering that the impact of shock wave on the pillar is mainly concentrated on the pillar close to the explosive, the structure within the range of the two pillars close to the explosive in the 1/4 model is analyzed. The upper and lower pillars are named as $C_1$, $C_2$, $C_3$ and $C_4$, respectively (Fig. 3(a)).

Fig. 7Nephogram of shock wave propagation at 90ms (charge of TNT: 20 kg)

We first analyze the platform plate in the station structure. Starting from the platform plate at the bottom of the explosive, a node is taken every 5 m and a total of 4 points are considered to obtain the displacement and time history curve of its $Y$-axis (plate displacement is mainly $y$-direction displacement), as shown in Fig. 8. It can be seen from Fig. 8 that the platform plate under the explosives has the maximum displacement, about 35.2 mm. As the distance between the nodes and the explosive center increases, the energy of shock wave decreases, and the peak displacement of the platform plate decreases. In addition, under the action of structural damping, the peak value of displacement will also gradually decrease with time [30].

Fig. 8Platform plate nodes at different distances from the explosion center displacement-time curve (charge of TNT: 20 kg)

Table 6 shows the maximum displacement response of the station structure under different explosive charges. The plate displacement in Table 6 is taken from the position on the same vertical axis as the explosive, and the plate structure displacement is mainly the vertical displacement ($Y$-axis); The two ends of the pillar are constrained, so the displacement of the pillar in Table 6 is primarily the horizontal displacement. (that is, the $XZ$ plane displacement of the model), and the maximum displacement is concentrated in the middle of each pillar. As can be seen from Table 6, in different explosives, due to the closest distance to the explosion source, the displacement response under the platform board explosion center is the most significant, reaching 13.3 mm, 22.8 mm, 35.2 mm, 45.6 mm and 68.8 mm respectively. The second is the middle floor directly above the explosion source, and the displacement of each pillar is mostly a few millimeters. It is reasonable that the displacement response of the structure increases with the charge of explosive.

4.3. Structural stress response of explosion

Fig. 9 shows the v-s stress nephogram of a 20 kg explosive model after the explosion. According to the effective stress in Fig. 9, the stresses in the platform, middle floor and the pillar near the explosive in the station structure are obvious after the explosion. Because the pillar closest to the explosive is under more significant stress, the effective stress in pillar $C_1$ and $C_2$ under 20 kg of explosive was studied. Fig. 10 is an effective stress nephogram on the pillar at different times. It can be seen in Fig. 10 that stress first appeared at the bottom of the pillar about 9 ms after the explosion. Then the stress propagated to both ends of the pillar, the obvious stress appeared at the junction of the pillar and the board. Besides, it can be found that there is a regular vibration of the stress value in the pillar, which is particularly obvious at the junction of the pillar and the middle floor. Then the effective stress in the pillar gradually weakens and dissipates. The results obtained with other explosives were similar.

Fig. 9V-s effective stress nephogram of the subway station at 50 ms

Because the pillar plays the role of supporting the overall structure in the subway station, it is the most important force-bearing component. The pillar C₁ located at the lower of the station is the pillar with the most stress in the station, so to further obtain the stress of the structure in the station, the pillar $C_1$ is taken for further analysis. The highest point, lowest point and middle point of the centerline of each face of pillar $C_1$ are selected as monitoring points. A total of 12 monitoring points on the pillar, as shown in Fig. 11. After the explosion, the pillar is exposed to the shock wave of the explosion, so tensile and compressive stress will be generated in the pillar. The pillar material is a reinforced concrete material, and its tensile performance is far less than its compressive performance. Therefore, the tensile stress of each monitoring point is used as a criterion for whether the pillar is damaged.

Fig. 10Effective stress nephogram of the pillar under explosion at different times (charge of TNT: 20 kg)

a)$t=$9.2 ms

b)$t=$22.4 ms

c)$t=$31.0 ms

d)$t=$70.4 ms

e)$t=$93.0 ms

f)$t=$126.8 ms

Table 7 is a statistical table of the maximum tensile stress at each monitoring point in the pillar $C_1$. It can be seen from Table 7 that under the same explosive charge, the tensile stress at each monitoring point of the pillar is different and the tensile stress distribution is uneven. In addition, by comparing the tensile stress of each monitoring point of the pillar under different explosive charges, it can be found that the maximum tensile stress in the pillar is always at the bottom of the pillar, followed by the top of the pillar, and the tensile stress in the middle of the pillar is the smallest. Further analysis of the tensile stress at the bottom of the pillar, in the case of different explosive charges, the tensile stress at the monitoring point b₃ is the largest among the four monitoring points at the bottom of the pillar. That is, the tensile stress at the bottom of the frontal surface perpendicular to the $Z$ axis is the largest. In the analysis of the explosion limit of reinforced concrete structures, according to the theory of equal strength, the structural ultimate tensile stress strength of concrete materials equivalent to reinforced concrete can be estimated to be 10.0 MPa [31]. In the data in Table 7, under the action of the explosive shock wave of 5 kg TNT explosive, the tensile stress values in pillar $C_1$ did not exceed the structural ultimate tensile stress, and the structure was safe. Under the action of 10 kg TNT explosive, most areas in the pillar Destruction; under the effect of TNT explosives of more than 20 kg, the structure in the pillar has been destroyed.

Fig. 11Selection of monitoring points

According to the data in Table 7, it can be obtained that the value of the tensile stress generated at $b_3$ in the pillar is the largest under the impact of the explosive shock wave. Therefore, the tensile stress at this point is used as the criterion for determining whether the pillar undergoes is damaged. When the maximum tensile stress at point $b_3$ is greater than or equal to the structural ultimate tensile stress of the material, the structure fails. According to the data in Table 7, the explosive charge and the tensile stress value at $b_3$ are linearly fitted, and the functional relationship between the two is approximately determined from a statistical perspective. The relationship between the tensile stress at the pillar $b_3$ and the charge of explosive can be obtained, as shown in Fig. 12, and the formula is as follows:

where $\sigma$ is the tensile stress at $b_3$, and $Q$ is the charge of TNT explosive.

Fig. 12The fitting formula of tensile stress value and explosive quantity at b3

According to the ultimate tensile stress of the station pillar, combined with the fitting formula, it can be roughly deduced that when the charge of TNT explosive is 6.77 kg, the tensile stress at the bottom of the pillar is equal to the ultimate tensile stress of the structure. At this time, cracks will occur at the bottom of the pillar, and the structure will be damaged.

The explosive model used in this paper is TNT, which is a simple substance explosive and has a wide range of uses in the military. At present, with the development of explosive manufacturing technology, the types of explosives used by terrorists have shown an increasing trend. In addition to TNT explosives, such as RDX, rock emulsion explosives, ammonium nitrate explosives, and various homemade explosives have also been used in terror. Therefore, in consideration of other types of explosives, the TNT equivalent method is used to convert the charge of TNT explosives that cause cracks at the bottom of the station pillar in this paper to other types of explosives. According to the relevant literature [32], the formula for calculating the TNT explosive equivalent of different types of explosives is:

where $W_S$ is the actual explosive charge, $W$ is the explosive charge of the TNT explosive, $Q_S$ is the explosive heat of the actual explosive, and $Q$ is the explosive heat of the TNT explosive.

According to the explosive heat of different types of explosives, combined with the formula, we can get the charges of explosives required to crack the bottom of pillar $C_1$ under the action of different kinds of explosives, which are shown in Table 8.