Abstract
Underground tunnels are vulnerable to terrorists’ bombing attacks, which calls for studies on tunnel’s response to internal explosive loading. In this paper, the dynamic response of a cylindrical tunnel to an ideal centric point explosion was treated as an axisymmetric 2-dimensional problem, in which the tunnel was modeled with a continuous anisotropic shell, while the ground medium’s effect was accounted for with linear elastic Winkler springs and the explosive loading described by a temporal and spatial function. The governing equation of the motion is a fourth-order partial differential equation, for which a numerical method combining finite difference with the implicit Newmark-β method was adopted. This method avoided complicated integral transform and numerical inverse transformation, thus allowing efficient parameter study. The maximum radial displacement was found on the cricle of the center of explosive, where hoop stress is the maximum principal stress. The anisotropy showed little influence on maximum hoop stress. Within the range of ground medium’s modulus, minor influence on maximum hoop stress was incurred. This research may be helpful to hazard assessment and protective design for some critical subway tunnels.
1. Introduction
In recent years, increased terrorist bombing attacks on subway system have been observed [1–3]; the most recent one was the notorious 2017 Saint Petersburg subway bombing [4]. Besides threats from terrorism, accidental explosion may be another potential hazard. The aftermath of an internal blast may include casualties and structural damage [5]. A deformation beyond threshold may cause the collapse of the structure. Also many of these tunnels may be built across beneath water bodies, such as these in New York, London, Shanghai, and Paris. A blast may breach or tear apart the wall of the tunnel lining and cause a flooding, for which the aftermath may be intractable [6, 7]. Thus the resistance of a tunnel structure to an internal blast should be studied.
In the designing of tunnels, the effects from some natural hazards and man-made disturbances such as seismic events and fires are generally considered, whereas blast loads are rarely considered for civilian tunnels [8]. Therefore, under the background of increasing probability of terrorist bombing attacks, safety of municipal tunnels under explosions needs more attention.
So far experimental study on tunnel’s response to internal explosion can be scarcely found in literature. Preece et al. [9] and Liu and Nezili [10] conducted a series of scale model tests using centrifuge facility, where the tunnel was simulated by monolithic aluminum tube. Their results may not be applicable to reinforced concrete tunnel since mechanical behaviors of metals like aluminum are quite different from that of reinforced concrete. Bonalumi et al. [11] conducted internal blast tests on a concrete pipe embedded in soft soil purposed to study internal blast loading in cylindrical enclosure, whereas the response of structure is not accounted for. Zhao et al. [12] carried out a few full-scale tests on staggered segmental tunnel’s response to internal point-source explosions, concluding that damage mainly concentrated in joint areas; however, no theoretical study was made on the tests.
Tunnels’ response to blast loading was more studied with finite element modeling. Based on the results of numerical simulation with ABAQUS, Colombo et al. [13, 14] plotted pressure-impulse (P-I) diagrams, which is also called iso-damage diagram for segmental tunnels in case of an internal explosion. This provided a useful tool for damage assessment and antiexplosion design of segmental tunnel lining. However, their results were obtained on the presumption of an axisymmetric plane problem, that is, the blast loading keeps constant along the tunnel axis, which is unlikely the true case. Choi et al. [15, 16] did nonlinear finite element analyses on internal explosion of tunnels and evaluated the vulnerability of lining structure, only concluding with some descriptive remarks. Yu et al. [17] also simulated internal explosion effects on tunnel in soil, and the influences from the soil and charge location were analyzed; however, there were no clear conclusions. These FEM analyses allowed for reproduction of the process involving blast wave propagation and its interaction with the tunnel and the structure’s response. However, these researches rely much on commercial finite element code, and the validity of the results depends on the material models’ parameters they used. Until now, no validation against any test results can be found in literature. Also it is not efficient to do parameter studies when the model scale is big.
As for theoretical investigation, a tunnel’s response to internal blast loads is generally treated as an axisymmetric plane strain problem, where the blast loading is simplified as uniformly distributed loading which keeps constant along the axial direction, and the tunnel is modeled with continuous isotropic thin shell [18–20], but this can hardly be the real case since the explosion could not be a centric infinite line charge detonation. Among those theoretical studies, Gao et al. [21] studied a 3D problem of point-source explosion’s effects on a long straight cylindrical tunnel in soil. However, the method includes complicated integral transform and numerical inverse transforms, which makes their research inconvenient for engineering use. Another problem in their research is that the orthogonal anisotropy of tunnel structure, that is, the difference in mechanic behavior between the circumferential and axial directions, which might be important in its overall dynamic response, was seldom accounted for.
In this paper a 3-dimensional problem of a long cylindrical thin-walled tunnel under a point-source explosion was studied, where both the lining structure and the blast loading were described more realistically. To begin with, a description of the blast loading on the inner wall of the structure with simplified spatial and temporal function was elaborated. Then the tunnel structure was modeled as a long straight anisotropic cylindrical thin shell, of which the governing equation of its motion was established on the basis of classical theory of shell. Thereafter a solving method combining finite difference method and Newmark-β implicit integration method was developed and some numerical results of shell’s dynamic response and parameter study were displayed and discussed. At last some concluding remarks were given on the basis of the analyses and discussions. The method developed in this paper may aid in both blast hazard assessment and protective design for some critical part of subway tunnel systems.
2. The Description of Blast Loads due to a Centric Point-Source Explosion
Considering a point charge of explosive locating on the centre axis. After detonation, the blast wave’s propagation and its loading effects can be categorized into two phases. First is a period of spherical pressure wave propagation and its interaction with the structure in the near-field. After a distance of 4~6 times’ radii from the explosive centre, it evolves into second phase; that is, after complicated interactions in first phase the blast wave propagates along tunnel axis like a 1-dimensional plane wave [22, 23]. Therefore, the loading on the tunnel from inner explosion is quite complicated.
Blast wave with pressure of propagates at a velocity of along the inner surface of a cylindrical shell with a certain length; the overpressure of incident blast wave could be described by a temporal and spatial function [21] as follows:where , are attenuation parameters of blast waves with respect to time and space, respectively; the centric point source of explosion is chosen as the origin of the coordinate system; is overpressure of incident blast wave, which could be estimated with empirical equation given by M. A. Sadovskii [24]:where is scaled distance, with and denoting stand-off distance (m) to the centre of explosive charge and weight of the charge (kg).
Knowing incident blast wave loading in confined space of tunnel's enclosure, the effective loading acts on inner wall of the lining due to the wave's reflection, that is, reflected overpressure can be calculated with the following empirical formula for normal relfection [25]:where is incident overpressure and is ambient atmospheric pressure.
Pokrovsky [24] proved by experimental data that when incidence pressure is smaller than MPa, reflected pressure of an oblique reflection can be calculated by the formula for normal reflection.
3. Governing Motion Equation and Its Solution
Taking account of the stiffness difference between circumferential and axial directions, the circular tunnel lining structure could be treated as an anisotropic cylindrical shell, as graphically shown in Figure 1, and the ground resistance effects are modeled with radial linear elastic Winkler springs.
3.1. Governing Equation of Motion
As aforementioned, an explosive charge is detonated on centre axis of a circular tunnel with radius of . The lining undergoes radial deflection fluctuation under the drive of radial explosive loading. For the axisymmetric problem of a shell's radial vibration, imagine cutting a unit element for forces' analysis, and all the forces acting on this unit element are diagramed in Figure 1, where denotes axial moment acting on the shell surface, denotes circumferential moment, denotes circumferential force, and denotes axial force. Due to the axisymmetry, other forces and moments vanish.
For elastic soils the lining-ground interaction could be described by linear elastic Winkler spring; therefore the resistance force iswhere is dilatational stiffness of ground springs, for which Penzien [26] gave the method of calculation aswhere and are Young’s modulus and Poisson’s ratio of the ground medium. Soils’ damping in structure-soil interaction was neglected for the simplicity of the problem, also because this approach leads to bigger amplitude of response, thus making the assessment more conservative.
Considering the balance equation of all the forces acting on the unit element and projecting all the forces to radial direction we can getwhere denotes radial displacement of the shell.
Taking account of and and simplifying (5) we reach the following:Suppose , , , and are axial and circumferential elastic module and passion ratio of anisotropic shell material, respectively; is the material’s mass density.
According to classical theory of shell [27] we havewhere , , and .
Also maximum stress can be calculated as follows:
Substituting (8) into (6) we get forced vibration equation of a thin cylindrical shell:
At the beginning time , the structure is stationary before the loading takes effects; therefore, the system has initial conditions as
Supposing that the shell is straight and infinitely long, however it is unrealistic; we simply take a far enough boundary that within the period of computation keeps stationary. Here we choose the boundary condition at as follows:
Together with (11) the problem becomes a system of mixed initial initial conditions and boundary conditions.
3.2. Solution of the Problem
In this subsection, a numerical method for solving the problem is explained. With finite difference of 2-order accuracy, the partial derivatives in (10) can be written as follows.
For , with forward finite differenceFor ,For , with backward finite differencewhere is the index of spatial points.
The initial conditions areAnd the boundary condition isThe unconditionally stable implicit Newmark-β method [28] was incorporated into this problem to solve the time integration.
For a general dynamic problem, we have [M][U′′] + [C][U′] + [K][U] = [F], where [M], [U], [U′], [U′′], [K], [C], [F] denote matrix of masses, displacements, velocities, accelerations, stiffness, damping coefficients, and forces.
The finite difference approximations for the Newmark-β method can be written as follows:If and , then the Newmark-β method is implicit and unconditionally stable, with at least two-order accuracy.
4. Numerical Results and Parameter Studies
4.1. Basic Numerical Results
Concerning a numerical calculation case of a cylindrical tunnel shell of reinforced concrete, the parameters are
Blast wave parameters are Pa, Pa, , and .
The surrounding ground was at first treated as homogeneous soil, and its Young modulus was taken as MPa, . Later, more stiff soils will also be accounted for as surrounding media in this section.
The results of radial displacement of the shell at varied times is plotted in Figure 2; it can be seen that radial vibration of the shell is symmetrical with respect to the location of charge; that is, . The radial expansion of the centre point rises to maximum value very shortly after the acting of the explosive loading and then gradually decays. After it return to its initial location, there is a process of rebounding, that is, vibrating inwards, causing a radial contraction of the shell, which is evident from the curve of displacement at 40 msec. Therefore, the vibration of the shell includes two different phases, that is, radial expansion and contraction. Also it can be found that as the amplitude of vibration of centre point decays, the nearby points begin to vibrate, and this vibration propagates along the axis in two symmetric directions like a wave.
With the time history of radial displacement of all the points on the shell along the axial direction, the maximal axial stress and hoop stress on all points can be obtained, which are plotted in Figure 3. With both symmetric geometry of structure and the loading condition, it is intuitive to deduce that the stress outcomes are also symmetric. Like the displacement curves, the maximal stress was reached at centre point, which is most close to the explosion. Also the evolution profiles are similar to that of displacement. Notice that on all the points hoop stress values generally outweighed axial stress, no matter during expansion or contraction. This proved that hoop, i.e., circumferential direction of the shell, was more strained during this process.
For clarity, the centre point was given special attention and its axial stress and hoop stress components were compared in Figure 4, which shows clearly that hoop stress and axial stress components fluctuated like two waves synchronically, while the amplitude of the former is always nearly two times as big as the latter one.
4.2. Parameter Studies
Since it is already clear that the most intense response happens on a circle centred by the explosive charge, we only need to check the stress response on centre points.
First the influence of anisotropy of shell, which is accounted for by the coefficient of , is studied; the comparison of the hoop stress on centre point is plotted in Figure 5. It can be seen that, within the range of , the influence of anisotropy on hoop stress of centre point is not pronounced at the beginning; only a minor influence on response magnitude can be observed. The influence however grows more manifest in the later part of the motion, which can be found from the lags between the series of curves. As a convention, the maximum response is more concerned in engineering; thus the anisotropy might be neglected in analysis of limit response.
Another factor may be the stiffness of the surrounding ground materials like rock or soil where the cylindrical tunnel was embedded in. Since Young’s modulus of elastic soils can vary in a big range, a range of MPa was considered for homogeneous soils in this analysis, while a range of MPa was tested for more stiff rocks. The results are plotted in Figure 6. It can been seen that, within the ranges of values we checked, the influence from ground medium's elastic parameter is also negligible. As we see this, is too small to affect the response significantly.
5. Discussion and Conclusions
In this paper, an ideal centric point explosion in a cylindrical thin tunnel was simplified as an axisymmetric 2-dimensional problem, in which the tunnel was modeled with a continuous anisotropic shell, while the surrounding ground medium’s effect was accounted for with linear elastic Winkler springs. The damping of the ground medium in structure-ground interaction was neglected for the simplicity of the problem, also because this approach leads to bigger amplitude of response, thus making the assessment more conservative. The explosive loading produced by a centric point-source detonation was described by a temporal and spatial function. The governing equation of the motion was established and a numerical method combining finite difference with the implicit Newmark-β method was developed for solving the problem. Some discussions and concluding remarks could be made as follows.
() A long straight anisotropic shell model was adopted to model dynamic vibration of a cylindrical tunnel under a centric point-source explosion, and the governing equation of the motion was established.
() Combining finite difference method with Newmark-β method, a numerical procedure was developed for solving the governing equation of shell’s vibration problem.
() Numerical results were given to show some profiles of the vibration response of the anisotropic shell to centric point-source explosion, from which a vibration propagation and decay from charge centre to far end were shown.
() The shell’s vibration consists of an expansion phase and a contraction phase, during which the maximum hoop stress always outweigh axial stress, suggesting a possible failure mode of axial cracks due to a tensile strength criterion met by hoop stress, since, for engineering materials like concrete, the tensile strength is much lower than compressive strength.
() Parameter studies showed that the coefficient of anisotropy does not influence the response amplitude significantly when changing in the range of 0.75 to 1.25; however, it generally alters the period of the vibration.
() From our investigation, Young’s modulus of surrounding elastic medium shows limited influence on the shell’s dynamic response, even when we check the value from 15 MPa to 150 MPa.
This research may be helpful to blast hazard assessment and protective design for some critical subway tunnels.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research is financially supported by Natural Science Foundation of China (Projects no. 51478469 and no. 51021001); the support from CSC scholarship (no. 201503170208) is also acknowledged by the authors.