Numerical analysis of underground tunnels subjected to surface blast loads

The increased terrorist and vandalism attacks on important public structures and utilities have raised the vital necessity for the investigation of performance of structures under blast loads to improve the design and enhance the behavior of structures subjected to such threats. In this study, 3D finite element analysis is used to study the effect of surface explosions on the response of RC bored tunnels. The soil behavior is modelled using Drucker-Prager Cap model. Two types of soil are investigated, and the blast load is considered through various weights of TNT explosive charges at heights of 0.50 m and 1.0 m from ground surface. To study the effect of horizontal standoff distance, six different horizontal distances are considered. The results show that the soil type has a significant effect on tunnel response due to surface blasts. Also the weight and the location of charge have a great effect on the safety of the tunnel. Finally, a parametric study is established to define the borders of the restricted area around the tunnel location to be safe.


INTRODUCTION
lasts have a devastating effect on structures, as well as, human beings. Blasts may affect humans by injury, health problems and loss of lives. When an explosive charge is near the ground surface, the super and shallow structures, foundations and shallow tunnels, are subjected to ground shock. When the charge is exploded near the ground surface, there will be two types of transmitted energy that cause ground shock. A part of energy goes directly through the soil causing a direct ground shock, while additional part is transferred through the air and compresses the ground surface and sends pressure into the ground under layers known as air-induced-ground shock [1]. If the explosion is directly above the surface, a conical-shaped crater is formed in the ground [2,3].

B
The pressure induced on structures by air, surface and buried charges can be calculated using many of empirical equations that is determined by experimental and analytical studies [4][5][6][7][8]. Most of these equations are used to calculate different parameters at a certain position such as maximum pressure for different types of soil, pressure history and crater dimensions. The destructive effects of blast loads have directed many investigators to study how the explosions affect the response of structures and the methods to dilute the risks using both tests and simulation models. The damage mode of tunnels in both dry and saturated soil were studied through experimental and numerical studies of blast loads [9][10][11][12][13][14][15]. The explosions tests are very limited due to their complexity and high cost. To reduce the number of blast tests, verified numerical methods are usually used to complement tests to provide parametric studies. According to researchers, computer packages such as ABAQUS, LS-DYNA, AUTODYNA and ANSYS are used to perform finite element analysis of structures through explosions [9,11,16,17]. The CONWEP algorithm available in the ABAQUS/Explicit is used to calculate the dynamic, time-dependent air-blast loading. The CONWEP algorithm implemented in ABAQUS/Explicit is based on the program developed by the U.S. Army to calculate conventional weapons effects. CONWEP implemented in software is used to model the blast load. The destruction of surface dynamic loads on underground structures depends on many factors such as the weight of TNT, the scaled distance, soil properties, position of structure from the TNT and the material of the structure [11,13,15,[18][19][20][21]. The layer thickness and moisture content also ground stiffness have significant influence on underground structures behavior such as lining stress and damage due to surface blasting, [14]. The extension of soil liquefaction is smaller if used large amount of explosive because the main part of the blast propagated energy transmits and destructs the sub-structure with small part fades into the surrounded soil. For saturated soils subjected to blast loads, the pressure and the particle velocity are higher than that of bulk soils, [22][23][24][25]. In current study, detailed Finite Element Analysis (FEA) using ABAQUS 6.14 is used to study the response of both of soil and Reinforced Concrete (RC) tunnel due to the surface blast. The soil effect was studied through two types of soils with 100 kg of TNT at 0.50 m of ground surface. The weight of TNT charges effect is considered through four separated weights of TNT charges of 100, 200, 300 and 400 kg of TNT at a height of 1.0 m above ground surface. Six different horizontal distances from the tunnel crown of 0.0, 1.   The model in the FEA is simulated using lagrangian three-dimensional solid continuum elements. The air and the TNT charge are not modelled in the research as the soil behavior and tunnel performance are the main parameters that are taken in consideration. The pressure resulted from the blast is only modelled using CONWEP implemented in ABAQUS, [23,26,27]. The soil behavior is simulated using Drucker-Prager Cap model in ABAQUS. This model presents perfect plasticity as well as isotropic hardening, [28][29][30]. Two types of soil (A and B) are used in the recent study. Soil A is selected to represent a firm low plasticity clay (USCS classification is CL) and the parameters were extracted from the laboratory test results in [28,31,32]. Soil B is selected to represent a very stiff clay (CH) and the parameters were extracted from [28]. Both of the two soils are used to represent and compare the behavior of the tunnel under different types of soil. The soil properties are obtained from [26,28], and are represented in Tab The soil and tunnel are modeled using C3D8R elements. To model the RC tunnel, the Concrete Damage Plasticity model is used. The concrete properties and the damaged plasticity model are summarized in Tab. 2 according to [26,33].
In the recent research, the global damping was neglected, [34][35][36]. In ABAQUS/Explicit a small amount of numerical damping is introduced by default in the form of bulk viscosity to control high frequency oscillations and to improve the modeling of high-speed dynamic events, [28][29][30]. ABAQUS/Explicit contains two forms of bulk viscosity, linear (default=.06), and quadratic (default=1.2), which can be defined for the whole model at each step of the analysis. In addition, Infinite elements and impedance conditions also add damping to a model.  The steel bars in RC tunnel are modelled as elastic-perfectly plastic materials and T3D2 elements are used for steel bars. The contact between bars and RC tunnel is simulated using embedded region contact available in ABAQUS. The steel properties is shown in Tab. 3. Fig. 2

MODEL VERIFICATION
or verification, the FE model is first produced without the tunnel structure, to show the propagation of shock waves inside the soil caused by surface blasts. A soil block with dimensions of 30 × 30 × 50 mm and 5m of non-reflected boundaries, in both X and Y directions, is used to model soil A (clayey soil), as shown in Fig. 4. The bottom side of the soil block is fixed in Z direction. A spherical charge of a weight 100 kg is used, and a scaled distance R= 0.5 m is set from the ground surface to the center of the charge.

Soil peak pressure
To confirm that the CONWEP will work well with the target surface X-Y surface, two finite element analyses are performed; an elastic-plastic analysis (FE Plastic) using the Drucker-Prager Cap model and a simple elastic analysis (FE Elastic). Eqn.
(1) by Nagy et al. [26] is used to evaluate soil peak pressure (PP). Eqn. (1) is based upon the empirical Eqn. 2, as presented in the technical manual TM5-855-1 [6] published by the Department of Defense of the US Army.
where R is the distance from the element to the charge center, W is charge weight with the TNT equivalent weight, and C and n are empirical constants that depend on the soil type.  P is peak free-field shock stress, f is a coupling factor, (ρc) is the acoustic impedance (ft/s), n\ is the attenuation coefficient, c is the seismic velocity (ft/s), ρ is mass density. The factors C and n have limits depending on soil type, where the constant C is a characteristic value from the TM5-855-1 design manual and depends on the charge material and soil properties. The constant n is the attenuation factor which influenced by the soil generally. In this case there is upper and lower limits for C and n because the soil parameters are not certain. The limits of C and n are [(1.12-0.65) and (2.75-2.50)] respectively. For the case of Nagy, et al. [26], Eqn. (1), the coupling factor is taken equal to 0.40, as the center of the 100-kg TNT spherical charge is in direct contact with the soil surface; Therefore, the bottom hemisphere of the spherical charge is buried into the soil. However, in this analysis, and in order to use the CONWEP option, the center of the charge must be located at a distance above the soil surface so any part of the charge is located above the surface. Therefore, the coupling factor should be reduced to 0.14 as in this case the center of the charge would be in air [6]. The comparison of the empirical equations of TM5-855 and finite element results is presented in Fig. 5. To get the pressure values directly without any interpolation, depths are chosen to correspond to the location of element centroids. Six points at distances of (1.00, 1.40, 2.60, 3.00, 4.00 and 5.40) m from the center of the TNT charge to the elements centroid are selected to represent the pressure time histories. The pressure-time histories in the soil at selected depth locations are shown in Fig. 6. All values of pressures at different scaled distances in elastic analysis are within range of upper or lower limits. Pressure in plastic analysis at scaled distances between 1 and 10 m/kg 1/3 , are within range with upper and lower limits. The soil pressure at all the selected depths has reduced lower than 0.5 MPa after 0.025 sec and that is a significant drop in the soil pressure.
Values of pressure with scaled distances less than 1 m/kg 1/3 didn't show a fair agreement with the empirical equations because of the formation of soil crater and the heaving of the soil.

Soil crater
When a blast wave hits the ground surface, the ground surface starts to scour and detonation gases penetrate the surface to form ejecta, which is thrown into the air [8]. The apparent crater radius and depth are measured relative to the original surface as shown in Fig. 7. At 0.03 second, the finite element analysis predicts the apparent crater radius, Ra, of 1.69 m and the apparent crater depth, D a , of 1.39 m as shown in Fig. 8.  To get the crater depth a D empirically, Eqn. 3 provided by [5] is used.
Where a V is The apparent crater volume, 3 m and depends on the cratering efficiency 0 (V ) of the explosive, 3 m /t , TNTequivalent charge weight (w), kg. and the height of the explosive (H) above the surface, m. Several values for the cratering efficiency of the explosive are considered as the water content of the soil in the finite element model is not known. The values are (312, 47 and 34 m3 /t) for wet clay, dry clay and dry alluvium correspondingly. To get the crater radius a R , Eqns.s 4 and 5 provided by [7,8]   is the density of the unreacted explosive, D is the detonation velocity, HOB is the ratio of the charge height above the surface to the charge radius. Tab. 5 illustrates the comparison of the empirical equations provided by [7,8,27]. The developed model is able to predict the apparent crater dimensions with a good level of agreement with the empirical equations.  or each study, two paths are chosen as shown in Fig. 9. A circular path called path 1 and a longitudinal path called path 2. Path 1 is around the perimeter of the tunnel and passes through its crown starts from node 4 to node 3 to crown at node 1 and then passes through node 2 and finally to node 4 again. Path 2 is 30m along the tunnel length and passes through the tunnel crown.

Effect of Soil Type
When the TNT charge is exploded, the wave is transmitted into the soil and hits the tunnel. Negative and positive periods of pressure are shaped along the tunnel body in all directions F Fig. 10 shows the kinetic energy (KE) measured for the whole volume of RC tunnel elements to show the effect of soil type on the behavior of the tunnel. The tunnel starts to response when the soil delivers the pressure at (t =11.25ms and 8.5ms) for soil (A) and (B) respectively. The response reached its crowning through 4.01ms and 1.76ms for soil (A) and (B) respectively. The tunnel dissipated 66.29% of its KE within 38.75ms of initial response period for soil (A) while for soil (B), the loss of energy is more and about 91.5% at 41.5ms.     From Fig. 11a-b and 12, it can be seen that in case of soil B the maximum pressure reached the tunnel is reduced by 58.9%.    As shown in Fig. 13.a-b and 14, the displacement is reduced by about 87%. The main reason for this reduction is the high stiffness of soil (B) regarding to soil (A).

Effect of Weight of TNT
In the second study, for cases 3, 4, 5 and 6, four different TNT charges(100, 200, 300 and 400 kg) were applied at constant scaled distance from ground surface, 1.0 m. Soil (A) is used in these cases. Fig. 15a and Fig. 15-b show the kinetic energy and plastic dissipation energy for the tunnel respectively. Fig. 15a-f show the pressure displacement along Paths 1 and 2 at the four charges at different time histories. Fig. 17a-b show the displacement along Paths 1 and 2 for the four cases.  The kinetic energy difference due to the TNT change is significant; the kinetic energy is increased by 64% when the charge weight increased from 100kg to 400 kg of TNT, (Fig. 15a). The same as the kinetic energy, the plastic dissipation energy increased as the charge weight increased to reach a difference of 65.81%between TNT= 100 and TNT= 400, (Fig. 15b). As shown in Fig. 16a-f, as the weight of TNT increases the pressure reaches the tunnel increases. The maximum pressure that reached the crown of the tunnel due to the four cases is 17.6 MPa, 21.96 MPa, 24.4 MPa, and 26.9 MPa. The pressure of the last case of 400 TNT is still bigger despite of running time, at 0.05 sec is still 13.43 MPa, despite of case 3, TNT=100, it is decreased to 4.80MPa and this is because of the enormous energy released from the weight of 400 TNT kg. From Fig.  17a-b, the maximum vertical displacement occurred in case 6 (TNT=400) is 8.69 cm and decreases as the TNT weight decreases to reach 3.7 cm for case 3 (TNT=100). Tab. 6 shows the maximum Von-Mises stress with different charges. As shown in Tab. 6, when charge reached 300 kg the stress reached the tunnel is more than the ultimate compressive strength of concrete 50 MPa. At great charges a protection like layers of foam or RC slab should be created for the tunnel to reduce energy conveyed to it, [16,37]. TNT

Effect of Horizontal Distance of TNT Charge
The last study, for cases (  As shown in Fig. 18a, the kinetic energy is reduced by 6.1% 13.3% 21.9% and 44.9% at R, 1.50, 3.0, 5.0 and 10.0 respectively. Also, the plastic dissipation energy is decreased to 12.64%, 40.79%, 62.68%, and 84.47% at R, 1.50, 3.0, 5.0 and 10.0m respectively, (Fig. 18b). As shown in Fig. 19a-c, the maximum pressure occurred at the crown of the tunnel in the six cases (7, 8, 9, 10, 11 and 12) are 15 Fig. 20a-b the maximum displacement at tunnel crown is 5.39 cm at R=0, while the displacement decreases to 2.61 cm at R=10.0m.
According to the previous results in this study, the scaled distance from the tunnel body shows a major effect on the tunnel response exposed to blast loads. Tab. 7 shows the the maximum Von-Mises stress with the six different horizontal distances. The stress reduces when the distance increases, the maximum stress is 40.39 MPa and occurs when the charge is directly above the tunnel when R=0.0 while the stress reaches 12.55 Mpa when R=10.0. In case of 100 kg of TNT the tunnel body will be safe without any additional precautions.

Restricted Area
A Vehicle Borne Improvised Explosive Device (VBIED) is usually used as a weapon of vandalism to destruct structures. VBIED can be carried in different types of vehicles to carry a huge amount of explosive with no suspicion. Air-blast pressure as well as ground shock are generated from the VBIED explosion because the VBIED is relatively nearby the ground surface. Tab. 8 shows the vehicles with the extreme load of TNT that can be supported without any suspicion [38]. The previous results from third parametric study can be used to investigate the efficient restricted area from the crown of the tunnel to be safe against any vandalism attack.
To study safe restricted area for tunnel, different types of vehicle (compact sedan (227 kg of TNT), sedan (454 kg of TNT) and cargo van (1814 kg of TNT)) are studied according to ultimate compressive Von-Mises stress of 50 MPa. The maximum Von-Mises for each case at different horizontal distances are shown in Tab. 9 to check the safe restricted area. The tunnel is safe for compact sedan with no need for restricted area but not for sedan or cargo van. The safe restricted zone for the sedan is 2.50 m from the crown of the tunnel (left and right the tunnel). For the cargo van the safe restricted zone is 5.0 m from the crown of the tunnel (left and right the tunnel).

CONCLUSIONS
finite element model is developed to study the performance of underground tunnel under surface blast load. The tunnel lining is reinforced concrete. Two types of soils are considered with different TNT weights and different horizontal distances. the following conclusions are obtained: 1) The developed model is able to predict the pressure propagating into the soil with a fair agreement to the empirical equations.
2) According to the kinetic energy (KE) of the tunnel, energy dissipation is more than in the weak soil (A) and the pressure reaching the tunnel is significantly affected by the type of soil. Stiff soils work as a protective shield and reduce both pressure and displacement of the tunnel.
3) The charge weight has a significant effect on tunnel response, as in high charges, the tunnel can't withstand the TNT energy and may be destructed so layers of protection materials like foam or RC slab should be used. 4) Based on the parametric study, the borders of the restricted area around the tunnel location to be safe is suggested. the tunnel is safe for compacted sedan but the safe restricted zone for the tunnel for sedan and cargo van is 2.50 m and 5.0 m respectively, from the crown of the tunnel (left and right the tunnel). 5) In future studies, the effect of moisture content in soil on the reaction of the tunnel under surface blast loads should be included and also field testing can be attained to develop a suggested design criterion for predicting the design of different types of tunnels under explosions.