Modeling Mitigation Effects of Watershield on Shock Waves

The object of this analysis is to investigate mitigation effects of watershield on air blast waves. To examine the water mitigation concept, features of the free-field detonation process are studied from a series of one-dimensional simulations using a multimaterial Eulerian finite element technique. Five different shock Hugoniots for water are compared, and the most accurate data are suggested. To verify the numerical procedure, results are compared with available experimental data for UNDEX problem and analytical predictions for air shocks. For the case of contact watershield, the magnitude of peak pressure generally decreases and the shock arrival time increases with increasing thickness of watershield. The total pressure impulse is reduced significantly at near field. Non-contact watershield was also examined, and was found to provide a better design criterion based on the further decay of peak pressure.

The object of this analysis is to investigate mitigation effects of watershield on air blast waves. To examine the water mitigation concept, features of the free-field detonation process are studied from a series of one-dimensional simulations using a multimaterial Eulerian finite element technique. Five different shock Hugoniots for water are compared, and the most accurate data are suggested. To verify the numerical procedure, results are compared with available experimental data for UNDEX problem and analytical predictions for air shocks. For the case of contact watershield, the magnitude of peak pressure generally decreases and the shock arrival time increases with increasing thickness of watershield. The total pressure impulse is reduced significantly at near field. Non-contact watershield was also examined, and was found to provide a better design criterion based on the further decay of peak pressure.

Introduction
Blast waves in air and underwater explosion and the consequent damage to structures have been of great interest in both defense and academic areas. Both stateof-the-art computational and experimental studies in these areas have been provided by many researchers.
However, the purpose of the current analysis is to investigate the influence of water placed in the near proximity of explosives, not necessarily in contact, on the generation and propagation of air blast waves. Water placed near the explosives can reduce the peak shock pressure and total pressure loading to nearby structures by absorbing the detonation energy. The storage cell test results demonstrated that water can reduce the peak gas pressure and total impulse by as much as 90%, at least for the range of test parameters [6]. The addition of 13.5 lb water reduces the peak gas pressure of 51.3 to 5.8 psi for the detonation of 4.67 lb TNT.
The high pressure shock waves produced by detonation propagate outward in all directions and aerosolize the water placed near the center of detonation. For this to occur, water medium delays the shock front and reduces initial peak pressure. Furthermore, the aerosolized water expanding with detonation products absorbs detonation energy of the explosive. Typical heats of detonation for TNT is 980 cal/g. On the other hand, 539 calories are required to vaporize 1 gram water. Hence, water absorbing the significant detonation energy can actually reduce the total pressure impulse.
Recent large-scale experiments performed in the KLOTZ-Club tunnel in Ålvdalen [5], however, did not show any effect of the water mitigation due to other variations in the experiment setup. In order to exploit the potential advantages of water concept, onedimensional spherically infinite computational domain has been considered. Basically, this is not an underwater explosion (UNDEX) problem, nor an acoustic water-air interface problem, but a non-linear wave propagation problem. To this end, the polynomial equation of state (EOS) of water which is based on experimental shock Hugoniot data of published works are obtained and implemented into a non-linear Eulerian finite element program which solves the effects of shock physics. For the case of contact watershield, the ratio of the watershield thickness (w) to initial explosive radius (t) is increased from values (w/t) of 0 to 1. w/t = 0 corresponds to air shock without watershield. For the case of non-contact watershield, ratios of air-gap (a) to watershield thickness (w) are 1 2 , 1 and 3 2 for w/t = 2 5 , and 1 2 and 3 2 for w/t = 3 5 . Examined parameters are the peak pressure, impulse and shock arrival time.
The equations of state for the materials considered in the study are described first. Numerical approach and verification is then addressed. The following two sections deal, respectively, with problem description and calculated results. The final section is a summary of this study and suggests some directions for further research.

Equation of state
The pressure-volume-energy behavior of the detonation product gases of the TNT is modeled with the standard Jones-Wilkins-Lee (JWL) EOS and detonation velocity of 6930 m/s [3]. The equation for pressure P is: where the parameters for TNT are A = 3.712 Mbar, B = 0.03231 Mbar, ω = 0.30, R 1 = 4.15, R 2 = 0.95, η = ρ/ρ 0 , ρ is the overall material density, ρ 0 = 1630 kg/m 3 (reference density), E = 4.29 × 10 3 kJ/kg (specific internal energy per unit mass). The Chapman-Jouguet pressure, P CJ , is 0.21 Mbar. Air is modeled as an ideal gas which uses a gamma law equation of state: where γ (=1.4) is the ratio of specific heats. The initial density of air, ρ 0 , is 1 kg/m 3 . To satisfy standard atmosphere pressure of 1 bar, the initial internal energy, E, is 2.5 bar, which is determined from Eq. (2) and the value of γ at unperturbed stage. The unit of E are the unit of pressure. The linear polynomial EOS is usually used to model gas with the gamma law EOS. In the simulations, this is achieved by setting higher order terms zero in the polynomial EOS. Water is one of the most extensively studied materials under shock loading, i.e., at high pressures, densities, and temperatures [8,9]. The shock compression data on water from experiments are analyzed using the Rankine-Hugoniot conservation equations which relate kinematic parameters to pressure, density, and internal energy in the shock state. There are two techniques for shock experiments on water. The first method is that the shock generated by plate impact compresses water sample. The second method is similar except for the explosive-driven shock. The range of shocked pressure in water is extensive. High explosives have been used to shock water to 100 GPa [12]. A two-stage light-gas gun has been used to achieve pressures up to 230 GPa [9]. Shock compression is also accompanied by an increase in temperature from 100 • C to several 10000 • C for shocked water. Once shock Hugoniot data are obtained from experiments, those can be implemented into a finite element code. The shock Hugoniot data is correlated by the cubic shock velocity-particle velocity equation, where, C 0 is the sound speed at undisturbed state, S 1 , S 2 , and S 3 are the coefficients of the slope of the U s − u p curve. U s and u p are the shock wave and fluid particle velocities, respectively. Most calculations showed that the water Hugoniot was the most important component in detonation study. Hence, five approximations from shock Hugoniot experimental data were compared due to uncertainties in the data for water. The first model was provided by Fusheng et al. [4] where two straight U s − u p curves were fitted from experiments. They observed a phase transition at 25 GPa. This pressure corresponds to the particle velocity of 3.5 km/s. The phase transition was accompanied by a volume contraction. Mitchell and Nellis [9], however, observed no indication of an energy absorption mechanism, like ionization, along the Hugoniot curve. This does not mean that such a process is not occurring, but only that it is not "turning on" strongly or reaching a threshold in this range. If an energy absorption process becomes significant, a change in slope should be observed. Shock wave and thermodynamic data pertaining to the Hugoniot centered 20 • C and 1 bar are tabulated in Table 1 [12]. There is a linear data from the HULL code. A linear curve fit of data determined by Sandia National Laboratories [15] is also available. A curve fit of data using the nonlinear U s − u p provided by Steinberg [13] is also widely used. Finally, the fifth model is the curve fit data provided by Lalle [8]. Lalle performed shock loading experiments on water in the 0-30 GPa pressure Table 1 Shock-wave and thermodynamic data for water Hugoniot [12] Pressure (GPa) V (cm 3  range to determine shock Hugoniot. He obtained three linear U s − u p curves with two kicks. The temperature for the first kick is 100 • C. The parameters for these models are summarized in Table 2. Figure 1 shows the comparison of the non-dimensionalized values of shock wave velocity and particle fluid velocity for above parameters. It is observed that U s − u p curve can not be described by a straight line. Steinberg shock Hugoniot model which fits well to the experimental data of Lalle [8] will be used for the current analysis. Based on the previous data, water can be modeled as a compressible fluid with a Mie-Gruneisen EOS which uses cubic shock velocity-particle velocity to define pressure for compressed and expanded materials. In compression (µ > 0), the pressure is given by, Fig. 1. Non-dimensionalized shock wave velocity vs fluid particle velocity.
and in tension (µ < 0), the pressure is given by, where, γ 0 is the Gruneisen gamma, and a is the first order volume correction to γ 0 . The density of water is 1000 kg/m 3 . The initial specific internal energy is 2.026 bar; it represents the specific internal energy necessary to give the water an initial pressure equal to standard atmospheric pressure. This is determined from Eq. (4) and 0.4934 of the value of γ 0 [13]. The polynomial equation of state can be also used to model water. In this expression, the pressure is related Table 2 Shock Hugoniot parameters for water to the relative volume and specific internal energy by a cubic equation. In compression, the pressure is given by, and in tension, the pressure is given by, where, µ = η − 1 and a 1 , a 2 , a 3 , b 0 , b 1 , and b 2 are constants for water. Constants for these equations were determined by fitting Mie-Gruneisen EOS [13] to the polynomial EOS. That is, the constants in Eqs (6) and (7) were determined, for both compression and expansion states, by matching terms in Eqs (4) and (5) with Steinberg shock Hugoniot data. The procedure was provided by Chisum and Shin [2]. The values of the atmosphere condition EOS parameters determined by this procedure, appropriate for condensation values on the order of µ < 0.8, are as follows: The initial specific internal energy is 205.9 J/kg which is determined from Eq. (6) and the above parameters; it also represents the specific internal energy necessary to give the water an initial pressure equal to standard atmospheric pressure.

Numerical modeling and verification
The multimaterial Eulerian finite element program MSC/DYTRAN [10] was applied to model the current problem including the explosion event, water shock, and blast waves propagation. The multimaterial Eulerian scheme in this code allows up to nine different Eulerian materials to be present in a given problem, and is likely suitable for the current study. The Eulerian scheme uses the basic conservation equations, in conjunction with constitutive equation and equation of state, to compute the solution in space. The description of the numerical approach is provided by Chisum and Shin [2].
To verify both the adequacy of the state equations and the numerical procedure for the simplest possible case, a comparison to experimental data provided by Weingarten [15] was conducted for a one-dimensional spherically symmetric problem. Four grams of Pento-lite including detonator were placed near the center of pressurized tank (1150 psi), and a pressure signature was measured at a distance of 13.97 cm from the explosive center.
Non-rectangular hexahedron elements, whose top and bottom were parallel but whose sides had a slope of 0.1, were used to model the TNT and fluid medium [2]. The mesh included 999 cells to a radius of 6 m. The first 0.00825 m to be filled with four gram TNT contains 25 equal length cells. The next 0.00825 m also contains 25 equal length cells. This region will be filled with water for the next problem. The cell size was then gradually increased by x 0 /x i = 40 until the mesh reached the outer boundary. x 0 and x i are the largest and smallest cell thicknesses, respectively. A flow boundary condition was specified in this outer boundary to satisfy the free-field condition. Pentolite is modeled with the standard JWL equation of state. The density of Pentolite is 1700 kg/m 3 . The values used for JWL Pentolite were obtained from the handbook by Dobratz [3].  Figure 2 shows the pressure signatures from the numerical results and the corresponding experimental data measured at a distance of 13.97 cm. The first peak is the result of the initial shock wave, while the second pulse comes from bubble collapse. Although the numerical result of the first peak pressure is 30% higher than the experimental data, the impulse is almost identical. The good agreement in pressure signature shows that the multimaterial Eulerian scheme is capable of calculating the UNDEX problem. Figure 3 illustrates the one-dimensional free-field model setup. The same geometry and computational elements with the previous case were used to study the effects of water mitigation on blast waves. To do this, four different cases were analyzed. The features of shock waves in both air and water medium only are first examined to provide a reference data and confidence in numerical procedure. Numerical results are verified using analytical equations for blast waves and empirical equations for UNDEX. The third problem studied consists of the TNT charge of equal mass and contact watershield with ratio of the watershield thickness, w, to initial explosive radius, t, varying from values of 0 to 1. w/t = 0 corresponds to air shocks without watershield.

Problem description
The final case studied consisted of the TNT charge of equal mass and non-contact watershield with various airgap thickness. The peak pressure, shock arrival time, and total pressure loadings are calculated at three   Table 3.

Shock waves in water
To explore the water behavior, the detonation process in one atmosphere water medium was considered in the section. Even though this is not a underwater explosion problem deep under the sea, it can provide a base line data and confidence in numerical modeling by comparing the results to the empirical equations. Pressure signatures calculated at three locations (A, B, and C in Fig. 3) from the numerical calculations show the typical exponential decay waveforms. At arrival time, t a , the pressure rises quickly to a peak value, P max , and gradually returns to the ambient pressure. Since there exists no high hydrostatic pressure, no secondary bubble collapse waves which comes from the pressure difference between the hydrostatic and bubble pressures are observed. In this case, the oscillation of this bubble is, if ever, not violent.
Empirical equations for predicting the UNDEX process have been developed for several explosives, such as peak shock pressure, bubble pulse period, and maximum bubble radius. In this case, the empirical equations for peak shock pressure can still be used, since it is a function of the weight of the explosive and standoff distance to the measuring locations. The approximation for peak shock pressure [14] for TNT is, where, A 1 = 1.18, K 1 = 22505, R is stand-off distance (ft), and W is charge weight (lb). Table 4 summarizes the results from the numerical calculations and empirical equation for three locations. The numerical results show a good agreement at close location (A), and less agreement in peak pressure at far location (C) when compared to the empirical results.

Air shocks
Blast waves formed in air using the same computational domain is considered in this section. The pressure signatures calculated at A, B, C for w/t = 0 are displayed in Fig. 4, showing the typical blast waveforms. The pressure rises quickly to a peak pressure, P max at t a , and returns to the ambient pressure. The pressure then drops to a partial vacuum state. At close location (A), a secondary shock is also observed in this time period. Numerous repeated shocks of small amplitude at various times after the first shock have indeed been observed [1]. These are caused by the successive implosion of rarefaction waves from the contact surface between explosion products and air.
Analytic equations to predict characteristics of air shocks have been provided by Kinney and Graham [7]. For computer calculations, it is convenient to have the overpressure-distance data in analytic form and given by, where P 0 /P a is the ratio of overpressure to ambient pressure and Z (m/kg 1/3 ) is the scaled distance away from the center with an energy release of one kilogram TNT in the standard atmosphere, where P 0 = P max − P a . Since shock velocity is uniquely related to overpressure ratio, the time required for that shock front to travel out to various distances can be also given by, where C 0 is the speed of sound in the undisturbed atmosphere and x c is the charge radius. Table 5 summarizes the results for the magnitude of peak pressure and shock arrival time (t a /kg 1/3 ) from the numerical calculations and the above analytic equations. Good Table 5 Comparison of calculations and analytic results of overpressure in air The magnitudes of peak overpressure in each location are underestimated compared to the results from the analytic equations. This can be due to the spatial discretization scheme to solve the discontinuities at shock front, such that the shock front is no longer discontinuity.

Contact watershield
The previous examples dealing with the shock waves in water and air blast waves provide confidence in numerical modeling. In this section, the effect of watershield in contact with the explosives on the generation of blast waves was examined. In Fig. 4, the calculated pressure signatures at A, B, C locations are displayed for ratio of the watershield thickness to initial explosive radius of 0, 2 5 , and 1. It is observed, as was anticipated, that the peak pressure continues to decrease with increasing the thickness of watershield and the shock arrival time is also delayed. This trend is realistic if we suppose the limiting case when the ratio of water thickness to explosive radius goes to infinite. Then, this becomes the underwater explosion situation, and the pressure in air should drop to atmosphere pressure. It is also interesting to note that the blast wave in air initially travels faster than in water due to the large Mach number, and high compressibility near the explosion, but it decelerates quickly.
The typical exponential pressure decay waveforms were observed up to w/t = 1 5 . For the case of w/t > 2 5 , it has the interesting and unexpected property that pressure cutoff is observed at close location (A). In order to interpret the pressure cutoff behavior, the interface between the watershield (including explosion gas products) and air was calculated and the prediction of interface growth is displayed in Fig. 5. In the figure, A is the closest measurement distance from the charge. Some gas bubble oscillations, which are not violent, can be seen. The oscillation period and maximum ra-    dius increases with increasing watershield thickness. It can be observed that when this interface reaches the observation locations, the cutoff starts to occur. This phenomenon seems to be attributed to the water state which absorbs the detonation energy of the explosive. The initial shock waves are able to change the water quickly into a mist of water droplets suspended in the atmosphere. By changing the water mist from a liquid state to a vapor state, huge amounts of energy released from detonation is dissipated. Hence, this phenomenon results in a reduced total pressure impulse. The test results demonstrated that water can reduce the peak gas pressure and total impulse by as much as 90% [6]. Fig. 9. Calculated peak pressures of non-contact watershield normalized by the peak pressure of no watershield, w/t = 2 5 , w/t = 3 5 . Figure 6 shows the results for the magnitude of peak pressures normalized by the peak pressure without watershield. Several things are evident from the result. Most of all, it can be observed that the influence of watershield is most significant at close location (A), where the magnitude of peak pressure quickly starts to decrease and it decays more than 30% for w/t = 1. On the other hand, at far location (C), the magnitude remains almost constant and decays only 2% for w/t = 1. The peak pressure at B location shows the intermediate characteristics. The normalized shock arrival time is displayed in Fig. 7. The influence of watershield is most significant at near field location, where shock front is delayed more than 70% for w/t = 1 case due to the non-linear wave propagation speed at early stage. The shock arrival time seems to be a lin- ear function of the ratio of watershield to explosive radius. Also plotted in Fig. 8 is the results for the total pressure impulses calculated at three locations. Although the impulse is significantly reduced at A, it seems to be difficult to separate any mitigation effect of the water at B and C locations. This is because at far locations neither the magnitude of the peak pressure is reduced nor does the water absorb detonation energy. According to this, significant reduced impulse is observed only in the region where water mist is suspended in the atmosphere, which is near field.

Non-contact watershield
Finally, airgap is maintained between the high explosive and watershield in order to investigate the per- Fig. 11. Impulses calculated at A, B, and C for w/t = 2 5 , a/w = 0, formance of non-contact watershield compared to the previous case. The ratio of the watershield thickness to initial explosive radius (w/t) are 2 5 , 3 5 . The ratio of airgap to watershield thickness (a/w) are 1 2 , 1 and 3 2 for w/t = 2 5 , and 1 3 and 3 2 for w/t = 3 5 . Figure 9 shows the results for the magnitude of peak pressures normalized by the magnitude of no watershield case. Obviously, it can be observed that the noncontact watershield produces less peak pressure than the contact watershield does. In the cases of a/w = 1 and 3 2 for w/t = 2 5 , peak pressures are less than the peak pressure would have been in the case of contact watershield for w/t = 2 5 (see Fig. 6). This is also true for w/t = 3 5 . That is, non-contact watershield provides better water mass efficiency. Again, this effect is most significant at near field location. The normalized shock arrival time is displayed in Fig. 10. The trend is similar to the contact watershield case. For w/t = 3 5 , the wave front is delayed further compared to the case of w/t = 2 5 . Also plotted in Fig. 11 is the results for the total pressure impulses calculated at three locations. At close location (A), the impulse is further reduced compared to contact watershield case. However, it is not reduced at B and C locations. This is due to the same reason which is discussed before. It can be concluded that this geometry provides a better design criterion.

Conclusions
The multimaterial Eulerian finite element approach has been used to model underwater explosion and air shocks. It is then used to model the blast waves with contact/non-contact watershield. The computational simulation included modeling the formation and propagation of water/air-shock waves.
It has been shown that water medium delays the shock front and reduces the magnitude of initial peak shock pressure. At a later time period, the water mist expanding with detonation products absorbs detonation energy of the explosive, and this can actually reduce the total pressure impulse. The peak pressure generally decreases with increasing thickness of watershield. The total pressure impulse is reduced significantly at near field, not far field. Finally, non-contact watershield was examined and it was found to provide a better design criterion based on the further decay of peak pressure and less total impulse. Two-dimensional calculations in either fully or partially confined boundary such as ammunition storage facilities was not considered here for the purpose of understanding the physics involved. Not only water but other potential materials need to study, although water is most cost effective. These area should be studied in the future.