Abstract
Studies of shock waves from explosions involve the dynamic analysis of geometric, material, and contact nonlinearities. The action of the shock wave on soil (pressurization and decompression) is completed within milliseconds. The energy transfer process of shock waves has a significant impact on engineering structures; an analyst evaluating the effectiveness of safety measures must first understand the dynamic characteristics of this energy transfer. Thus, examining the movement patterns of the media through which shock waves travel is of primary concern. Construction designs must consider the propagation characteristics of shock waves in soil. In this study, the finite element method was employed to analyze the dynamic and mechanical behavior of soil by applying the multimaterial arbitrary LagrangianEulerian algorithm to a 3D solid element with eight nodes and six faces. A structural model of the fluidsolid coupling was constructed to examine the shock wave propagation process resulting from shallow underground explosions. In addition, ground surface acceleration values obtained from field tests were used to validate results of the numerical analysis of time and spacedependent shock wave propagation and attenuation in soil. The objective of this study was to perform a rational analysis of the seismic effects of shallow underground explosions to provide a reference for earthquake engineering.
1. Introduction
Explosions rapidly release energy, causing a sudden increase in temperature and pressure in the surrounding medium, resulting in a phenomenon known as overpressure. Blast waves cause surface particles to vibrate, and the inertia affects material stability. Seismic waves in the ground are primarily produced by shock waves in the air. A seismic wave first propagates in the form of a shock wave, which is later converted to an elastic seismic wave, causing surface particles to vibrate [1, 2]. Surface vibrations from a shallow underground explosion form a spherical shock wave in the soil. Shock waves near the blast source exhibit high velocity and magnitude, resulting in higher velocity and magnitude in the corresponding soil compression wave [3]. With the same charge weight, the error of margin of the peak blast pressure is 15 % compared with experimental values [4]. In shallow explosions, the blast pressure compresses the soil, forming a crater and creating shock waves in both the air and soil. If a charge is buried deep underground, as is the case in drilling and blasting, the blast effect is insufficient to displace the soil. Thus, the detonation cannot directly create a shock wave in the air, and a shock wave forms only in the soil [5]. The mechanical behavior of an explosion is complex; the load time is short and the vibration frequency is high. For surface or shallow underground explosions, the magnitude of destruction is analyzed on the basis of the magnitude of surface vibrations.
The seismic effect of an explosion refers to the seismic phenomenon generated in a defined area as a result of the detonation of a charge within a geological environment. Crandle [6] analyzed the propagation of blast waves and found that when seismic waves produced by an explosion propagate through a medium, the surrounding particles vibrate in response to the process. When the vibration reaches a specific magnitude, structures within a defined area are affected to varying degrees, even to the point of destruction. Koga and Matsuo [7] examined the detonation of charges in surface soil and found that a portion of the energy is released as waves traveling outward through the medium. In the propagation process, a portion of the energy is converted into seismic waves that cause the particles in the medium to vibrate, resulting in vibration on the ground surface. Pullen et al. [8] constructed a simple soil composition model in DYNA2D, a hydrodynamic finite element code, to simulate the propagation behaviors of a complex stress wave. Analysis results of the initial pressure, strain, and impulse roughly corresponded to the experimental results.
The seismic effect of an explosion causes surface particles to vibrate. The motion characteristics of the medium are key parameters affecting the blast load. Particle displacement, velocity, and acceleration can be applied to evaluate the basic parameters of the vibration effect. Because an underground explosion causes surface vibrations, geological structures are subjected to timedependent effects such as displacement, velocity, acceleration, stress, and deformation. Inertia, which is produced by resistance to the acceleration of a structure with mass, affects the stability of the structures. In studies that have analyzed explosions or earthquakes, the vibrational velocity, acceleration, and displacement of particles have been used as physical parameters to assess the magnitude of blast vibrations [912]. In the present study, a numerical analysis model was constructed and validated by using ground surface acceleration readings obtained through explosion field experiments. The energy of the shock waves caused by shallow underground explosions and the area of effect were also analyzed. The results may serve as a reference for antiblast safety engineering.
2. Shallow underground explosion field experiment
Fig. 1 shows the configuration of the field site for the explosion field experiment. The experiment used 0.25 lb (113.389 g) of TNT buried 30 cm beneath the ground surface. Triaxial accelerometers were placed 350, 500, and 650 cm from the center of the charge to measure the ground surface acceleration in the $X$, $Y$, and $Z$ directions for a quantitative analysis of the propagation and attenuation of the ground surface motion caused by the underground explosions. Oscilloscopes, signal conditioning processors, power supplies, and data acquisition systems were also used in the field experiments. The signals produced must pass through signal conditioning processors linked to an oscilloscope for signal transmission and storage.
Fig. 1Configuration of the explosives in the field experiment
3. Finite element analysis methods
LSDYNA, a nonlinear finite element analysis code based on continuum mechanics theories, was employed to analyze the shallow underground explosions and the resulting propagation of blast waves in soil. LSDYNA primarily uses an explicit timeintegration method to solve dynamic timeintegration problems and is a conditionally stable calculation method. Explicit timeintegration methods are based on the d’Alembert principle. The single degreeoffreedom elastic damping system is expressed in Eq. (1), a nonlinear system is expressed in Eq. (2), and the matrix equation for the motion with multiple degrees of freedom is expressed in Eq. (3). This code is advantageous because it can be used to analyze nonlinear and large deformation behaviors in 3D structures, and is particularly suitable for solving nonlinear dynamic problems such as highvelocity collisions and explosions [13]:
where $m$ represents mass; $c$ denotes damping; $k$ refers to stiffness; $\ddot{u}$ represents acceleration; $\dot{u}$ is velocity; $u$ represents displacement; $p\left(t\right)$ denotes external forces matrix; $\ddot{U}$ represents the acceleration matrix for the node; $\dot{U}$ is the velocity matrix for the node; $U$ denotes the displacement matrix for the node; $P\left({t}_{n}\right)$ represents the external forces matrix; $M$ is the mass matrix; $C$ refers to the damping matrix; and $K$ is the stiffness matrix.
Explicit timeintegration based on the central difference scheme is suitable for solving nonlinear problems. The method assumes the relationships among acceleration, velocity and shifting are expressed in Eqs. (4) and (5), respectively. The variance between the above two equations is $\mathrm{\Delta}{t}^{2}$. To find the answer for shifting at the time $t+\mathrm{\Delta}t$, consider the dynamic equation for time ($t$), as shown in Eq. (6). Apply procedures Eqs. (4) and (5) in (6), as shown in Eq. (7). Procedure Eq. (8) reveals that ${U}_{t+\mathrm{\Delta}t}$ is apparently derived from the shifting of the previous time; only the balance conditions at time ($t$) need to be taken into consideration. This is the benefit of explicit time integration:
where ${\ddot{U}}_{t}$ represents the acceleration matrix for a node at time $t$; ${\dot{U}}_{t}$ is the velocity matrix for a node at time $t$; and ${U}_{t}$ is the displacement matrix for a node at time $t$.
Explicit timeintegration requires a timestep $\mathrm{\Delta}t$ to be smaller than the critical time interval; this is the conditional stability of explicit timeintegration. However, $\mathrm{\Delta}t$ is calculated differently for different element types. When calculating time $T$, LSDYNA automatically divides $T$ into $T/\mathrm{\Delta}t$ cycles; thus, the time required for code calculations depends on $\mathrm{\Delta}t$. In this study, a solid element comprising eight nodes was used; the conditions for stability are expressed in Eqs. (8) and (9), and the characteristic length is given in Eq. (10). The propagation of velocity in an elastic material with a constant bulk modulus is shown in Eq. (11). The LSDYNA user’s manual suggests using a timestep control coefficient smaller than 0.67 when analyzing blast effects [14]:
where $Q$ is the function of volume viscosity coefficients ${C}_{a}$ and ${C}_{b}$; ${L}_{e}$ is the characteristic length; ${\dot{\epsilon}}_{kk}$ is the strain rate tensor; ${v}_{e}$ is the element volume; ${{A}_{e}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the area at the longest side; ${c}_{v}$ is the sonic velocity in materials; $\rho $ is the density; $E$ is Young’s modulus; $G$ is shear modulus and $\upsilon $ is Poisson’s ratio.
3.1. Numerical computational model
In this study, the multi–material arbitrary LagrangianEulerian (MMALE) algorithm was employed to analyze the time history of shock wave dynamic reactions produced by the shallow underground explosion of TNT. To construct a numerical analysis model of the fluidsolid coupling, the explosives and air were defined as a Eulerian grid and the soil as a Lagragian grid. The numerical analysis involved a 3D solid element with eight nodes and six faces, as well as explosives, air, and soil defined as MMALE grids. The fluid and solid grids were independent of each other. The fluid–solid coupling was constructed by overlapping the grids. The setting CONSTRAINED_LAGRANGE_IN_SOLID was applied to construct a numerical analysis model of the fluidsolid coupling based on the MMALE algorithm. This model enabled the analysis of the time and spacedependent propagation and attenuation of shock waves in the soil.
Fig. 2 shows a quartersymmetry model used for numerical analysis. The numerical analysis model was constructed according to the explosion field test conditions by using the following units: $cm$, $g$, $\mu s$. To exploit the symmetry of the model, a quartersymmetry model with nonreflecting boundary conditions was used for the analysis. For the field experiment, the rectangular TNT charges weighed 0.25 lb (113.389 g) and had a density of 163 g/cm^{3}. In the quartersymmetry model, the dimensions of the charges were 1.64×1.64×9.3 cm. The charges were buried at a depth of 30 cm and the detonation point was located at the center of the charge. Air was defined as an ideal gas. The ideal gas law was applied with the dimensions of 190×660×120 cm in the quartersymmetry model. The dimensions of the soil in the model were 180×660×120 cm. Based on the convergence analysis of the numerical model, the grid density was defined as the smallest width of the TNT; the Eulerian grid was defined as 0.5fold the smallest side of the charge; and the timestep control coefficient was defined as 0.3 [13, 14].
Fig. 2Quartersymmetry model used in the numerical analysis
Table 1 lists the material parameters for the explosives, air, and soil. To model the air, the material identified as #9 (MAT_NULL), in the LSDYNA code was combined with the EOS_LINEAR_POLYNOMIAL EOS to describe the characteristics of the material, as shown in Eq. (12) [14]. The parameters used were density (RO), pressure cutoff (PC), dynamic viscosity coefficient (MU; $\mu =1/V1$), relative volume (TEROD and CEROD), Young’s modulus (YM), Poisson’s ratio (PR), initial energy per unit volume (${E}_{0}$), and ${C}_{1}$, ${C}_{2}$, ${C}_{3}$, ${C}_{4}$, ${C}_{5}$, ${C}_{6}$ are constants. For ideal gases, the coefficients in the LINEAR_POLYNOMIAL EOS are ${C}_{1}$, ${C}_{2}$, ${C}_{3}$, ${C}_{6}=0$ and ${C}_{4}={C}_{5}=\gamma 1$. The EOS is expressed in Eq. (13). The parameters used were $\gamma $ represents the rate of change to the specific air temperature, $\rho $ represents the current air density, ${\rho}_{0}$ is the initial air density value and $\rho /{\rho}_{0}$ denotes the relative density:
Table 1Material parameters of soil, air, and TNT explosive
Element  Material and equation of state parameters (unit system: g, cm, μsecond)  
Soil  MAT_SOIL_AND_FOAM  
RO  $G$  BULK  ${A}_{0}$  ${A}_{1}$  ${A}_{2}$  PC  
2.6  0.000147  0.00729  7.105E07  4.198E08  6.215E10  0.0  
Air  MAT_NULL  
RO  PC  MU  TEROD  CEROD  YM  PR  
0.00129  0.0  0.0  0.0  0.0  0.0  0.0  
EOS_LINEAR_POLYNOMIAL  
${C}_{0}$  ${C}_{1}$, ${C}_{2}$, ${C}_{3}$, ${C}_{6}$  ${C}_{4}$  ${C}_{5}$  ${V}_{0}$  ${E}_{0}$  
–1.07E06  0.0  0.4  0.4  1.0  2.53E06  
TNT  MAT_HIGH_EXPLOSIVE_BURN  
RO  $D$  PCJ  BETA  $K$  $G$  SIGY  
1.63  0.693  0.21  0.0  0.0  0.0  0.0  
EOS_JWL  
$A$  $B$  ${R}_{1}$  ${R}_{2}$  OMEGA  ${E}_{0}$  ${V}_{0}$  
3.712  0.03231  4.15  0.95  0.3  0.07  1.0 
To model the explosives, the material identified as #8 (MAT_HIGH_EXPLOSIVE_BURN), in LSDYNA was combined with the JonesWilkinsLee (JWL) EOS to describe the characteristics of the high explosives, as shown in Eq. (14) [15]. Relevant parameters taken from the Lawrence Livermore National Laboratory Explosives Handbook [16] were the density (RO) of the material, detonation velocity ($D$), ChapmanJouget pressure (PCJ), beta burn flag (BETA), bulk modulus ($K$), shear modulus ($G$), and yield stress (SIGY). The parameters for the EOS were relative volume ($V$), initial energy per unit of volume (${E}_{0}$), and internal energy of the material (${E}_{m}$); $A$, $B$, ${R}_{1}$, ${R}_{2}$, and $\omega $ were constants representing the characteristics of the explosives:
Soil is a compressible porous material. The Krieg yield criterion is based on the theory of isotropic elasticity, which states that the behavior of materials at the yield point can be differentiated into hydrostatic pressure and shear stress. This model is suitable for analyzing material failure under stress. The yield function ($f$) is expressed in Eq. (15). LSDYNA was used to calculates the shear failure criterion according to the relationship between the average stress and extent of damage [17]. The material model MAT_SOIL_AND_FOAM in LSDYNA was used to analyze the dynamic conditions of soil after being subjected to blast effects [14, 15]. The soil parameters were density (RO), shear modulus ($G$), bulk modulus ($K$), pressure cutoff ($C$) and ${A}_{0}$, ${A}_{1}$, ${A}_{2}$ is the yield function constant for plastic yield function below. Indoor experiments of the soil showed that its water content was 11.8 %, $c=\mathrm{}$8.43 kPa, and $\varphi =\mathrm{}$14°. Under the Unified Soil Classification System, this soil was classified as SPSM, which is poorly graded sand containing silt:
with ${A}_{0}={c}^{2}$, ${A}_{1}=2c\mathrm{t}\mathrm{a}\mathrm{n}\varphi $, ${A}_{2}=\mathrm{t}\mathrm{a}{\mathrm{n}}^{2}\varphi $. Where ${J}_{2}$ is the second invariant of the deviatoric stress; $p$ is the pressure; and ${A}_{0}$, ${A}_{1}$, ${A}_{2}$ are the parameters related to the shear on the yield surface.
4. Results and discussion
4.1. Analysis of vibrations recorded in the explosion field test
This study was conducted to analyze the propagation characteristics of shock waves and the dynamic reactions of the ground surface caused by shallow underground explosions, as well as to use the quantitative readings of ground surface acceleration to verify the results of the numerical analysis. To determine the blastwave propagation characteristics, the ground surface acceleration was measured at 350, 500, and 650 cm horizontally from the surface of the detonation. Table 2 shows the peak ground surface acceleration at points along the $X$, $Y$, and $Z$ axes. Figs. 3, 4, and 5 show the acceleration time history of the ground surface acceleration at 350, 500, and 650 cm from the blast source, respectively. The results show that given the same charge weight, the vertical ground surface acceleration is faster than the horizontal surface acceleration. In addition, the acceleration decreases with distance, and the shock wave was substantially weaker at 500 cm from the blast source. Attenuation occurred more rapidly in the vertical direction than in the horizontal direction. The overall trend corresponded to the attenuation characteristics of blast wave energy.
Table 2Peak ground surface acceleration along the X, Y, and Z axes
Distance (cm)  Field experiment  Numerical analysis  Relative of error  
Horizontal  Vertical  Horizontal  Vertical  Horizontal  Vertical  
$Z$ axis  $X$ axis  $Y$ axis  $Z$ axis  $X$ axis  $Y$ axis  $Z$ axis  $X$ axis  $Y$ axis  
(gal)  (gal)  (%)  
350  0.699  0.671  0.802  0.628  0.601  0.737  –10.16  –10.43  –8.11 
500  0.680  0.552  0.783  0.616  0.504  0.725  –9.41  –8.70  –7.41 
650  0.246  0.158  0.293  0.220  0.143  0.266  –10.57  –9.49  –9.22 
Fig. 3Time history of the surface acceleration at 350 cm from the blast source
Fig. 4Time history of the surface acceleration at 500 cm from the blast source
Fig. 5Time history of the surface acceleration at 650 cm from the blast source
Fig. 6Time history of the surface acceleration at 350 cm from blast source (numerical analysis)
Fig. 7Time history of the surface acceleration at 500 cm from the blast source (numerical analysis)
Fig. 8Time history of the surface acceleration at 650 cm from the blast source (numerical analysis)
4.2. Validation of the numerical model
Fig. 6, 7, and 8 show the acceleration time history obtained from the numerical analysis of the ground surface acceleration at 350, 500, and 650 cm from the blast source, respectively. The vertical ground surface acceleration was used to validate the results of the numerical analysis. The relative margin of error was calculated by subtracting the experimental results from the numerical analysis results, divided by the experimental reading, and multiplied by 100. The vertical ground surface acceleration relative margins of error at 350, 500, and 650 cm from the blast source were –8.11 %, –7.41 %, and –9.22 %, respectively. All margins of error were within 15 %. According to previous research [4], these results validate the reliability of the numerical analysis model and confirm that the MMALE algorithm can accurately analyze the propagation characteristics of shock waves from explosions.
4.3. Dynamic characteristics of the blast wave
An explosion is a phenomenon involving a rapid release of energy. The high temperatures and pressures quickly decompress the surrounding medium, resulting in shock waves. The blast wave affects the internal stability of materials. Structural designs must include improved antishockwave protection measures to ensure safety. Ground surface vibrations caused by explosions are the result of energy transfer and conversion. The propagation of blast waves through the air is termed airinduced ground motion, whereas that through surface materials is referred to as a directly transmitted ground shock. An explosion produces a spherical shock wave that causes vibrations in the ground surface. Surface and underground explosions cause shock waves in the soil. Shallow underground explosions cause shock waves to form in both the air and soil. Drilling and blasting underground explosives cause shock waves directly in the soil and motion in the ground surface, both of which are affected by the distance from the explosion.
Fig. 9Time history of changes in shock wave pressure in a soil medium
a)$t=$ 100 μs
b)$t=$ 2,000 μs
c)$t=$ 4,000 μs
d)$t=$ 6,000 μs
e)$t=$ 8,000 μs
f)$t=$ 10,000 μs
Fig. 9 shows the changes in shock wave pressure in a soil medium over time. The energy released during an explosion forms a dynamic pressure wave that produces a shock wave in air or soil. The blast transmits a stress wave through soil in the form of a compression wave. The blast pressure pushes outward as a sphere, but the shock wave energy decays within milliseconds. Explosions cause stress in the medium to propagate outward, forming a cavity at the center of the explosion. Surrounding the cavity, from innermost to outermost layers are the fracture, elasticplastic, and elastic zones. In the propagation process, the shock wave first forms an elastic shock wave and then a plastic shock wave. Plastic shock waves are affected by unloading disturbances. The pressure and wave velocities decay continuously, and the impedance changes with the pressure. When the pressure wave reaches a free surface and is reflected in the form of a tensile wave, the resulting action of the expansion wave, pressure wave, and pressure of the explosive gases generates tensile waves and shear waves at the ground surface. The analysis results show that the shock waves near the blast source were fast and strong. These results correspond to the compression waves in the soil, which are key parameters influencing ground surface vibrations.
5. Conclusions
In this study, field tests were conducted and LSDYNA was employed to examine the formation of highstress shock waves during explosions and their characteristics and actions in the soil by exploring the dynamic reactions of the ground surface under identical energy conditions during an explosion. The results obtained on the propagation of shock waves in soil and the characteristics of the dynamic reactions at the ground surface may serve as a reference for future studies on planning and design of protective constructions. The results show that using the finite element method to simulate the dynamic characteristics of soil after an explosion provides accurate simulation results. Characteristics of the stress wave propagation were determined by examining the dynamic behavior of soil after being subjected to a blast. The largest and primary stress wave weakened with increasing distance from the blast point. The overall trend observed by numerical analysis matched the attenuation characteristics of the blast wave energy, indicating that the MMALE algorithm can accurately analyze the dynamic reactions of shock waves in soil. An understanding of the phenomena involved in the propagation of vibrational waves through the ground surface from the perspective of controlling vibrations caused by explosions is necessary for blastresistant planning and design in construction projects.
References

Structures to Resist the Effects of Accidental Explosions. Technical Manual TM51300, US Army Engineers Waterways Experimental Station, Washington, DC, 1990.

Wang J. Simulation of Landmine Explosion Using LSDYNA 3D Software: Benchmark Work of Simulation in Soil and Air. Aeronautical and Maritime Research Laboratory DSTOTR1168, 2001.

Wang Z. Q., Y. Lu. Numerical analysis on dynamic deformation mechanism of soils under blast loading. Soil Dynamics and Earthquake Engineering, Vol. 23, 2003, p. 705714.

Kivity Y., Shafri D., BenDor G., Sadot O., Anteby I. The blast wave resulting from an accidental explosion in an ammunition magazine. International Symposium on Military Aspects of Blast and Shock Conference, Canada, 2006.

Fundamental of Protection Design for Conventional Weapons. Technical Manual TM 58551, US Army Engineers Waterways Experimental Station, Washington, DC, 1998.

Crandle F. J. Ground vibration due to blasting and its effect upon structures. Boston Society of Civil Engineering, Vol. 36, 1949, p. 222245.

Koga Y., Matsuo O. Shaking table tests of embankments resting on liquefiable sandy ground. Soils and Found, Vol. 30, Issue 4, 1990, p. 162174.

Pullen A. D., Ianucci L., Newman J. B., Perry S. H. Experimental Investigation of the dynamic loading of soil and the use of finite element simulation for analysis. International Conference – Structures Under Shock and Impact, 1989, p. 443454.

Stamatopoulos C., Petridis P., Bassanou M., Allkja S., Loukatos N., Small A. Improvement of dynamic soil properties induced by preloading verified by a field test. Engineering Geology, Vol. 163, 2013, p. 101112.

Dowding C. H. Blast Vibration Monitoring and Control. PrenticeHall, Englewood Cliffs, NJ, 1985.

Ak H., Iphar M., Yavuz M., Konuk A. Evaluation of ground vibration effect of blasting operations in a magnesite mine. Soil Dynamics and Earthquake Engineering, Vol. 29, Issue 4, 2009, p. 669676.

Wang I. T. Numerical analysis of the dynamic response of slope under the blasting vibration effect. International Journal of Computer and Electrical Engineering, Vol. 6, Issue 4, 2014, p. 351357.

LSDYNA Theoretical Manual. Livermore Software Technology Corporation, 1998.

LSDYNA Keyword User’s Manual. Livermore Software Technology Corporation, 2009.

LSDYNA Theory Manual. Livermore Software Technology Corporation, 2006.

Dobratz B. M. LLNL Explosive Handbook Properties of Chemical Explosives and Explosive Simulants. Lawrence Livemore National Laboratory, 1981.

Len S. Geomaterial Modeling with LSDYNA. Livermore Software Technology Corporation, 2001.
About this article
This study was partially supported by the Ministry of Science and Technology, Taiwan, R.O.C. through Grant MOST 1042221E145004.