Comparation of Analytical Models and Review of Numerical Simulation Method for Blast Wave Overpressure Estimation after the Explosion

A R T I C L E I N F O A B S T R A C T Article history: Received: 14 December, 2020 Accepted: 19 January, 2021 Online: 05 February, 2021 A comparative analysis of formulas for blast wave overpressure is presented in the paper, and models were compared with available experimental data. The Kinney and Shin models show the best agreement with experimental data (Kingery-Bulmash) for free airburst, while for surface burst, Swisdak, Vanuci, and Jeon models predict test data most accurately. One of the novelties in the paper is introduction of new exponential and power functions for blast overpressure estimation, giving good agreement with experimental data. Also, several numerical simulations of free airburst explosions were performed to introduce methodology, and compare the data obtained with experimental data. A detailed description of the procedure for these simulations was provided – a contribution to numerical modeling of blast wave phenomena.


Introduction
There has been an increase in the research of explosive blast effects since numerous accidents (ie Beirut, 2020) happened, and terrorist attacks were carried out on civilians, with a large number of casualties during last decades. These include i.e. the bombing of US embassies (Ankara, 1958, Beirut, 1984, World Trade Center attack (New York, 2001), Mumbai attacks (2008), Mariott Hotel attack (Islamabad, 2008), Baghdad bombings (2016), Atatürk Airport attack in Turkey (2016), Nairobi hotel attack (Kenya, 2019), Jolo Cathedral bombings (Philippines, 2019), many attacks in Pakistan (Quetta, 2020) and also in larger area of Afganistan (lately Kabul, Kuz-Kunar, 2020) and many more.
An explosion presents violent release of energy. The detonation presents a rapid chemical reaction proceeding through the explosive charge at supersonic speed. The detonation wave converts explosive into a hot, dense, high-pressure gas, which is the source of powerful blast waves around explosion site ( Figure  1).
Pressures at the Chapman-Jouguet point can go up to 400 kbar. Acceleration of the air particles (when the detonation is complete) beyond the face of the explosive is in order of 10 11 g. Close to the face of the charge, at near-field scaled distance (where Z=R/m 1/3 ): Z = 0,054 m/kg 1/3 , the maximal temperature is around 200 000 K. The temperature rises almost instantaneously with the shock front arrival but the maximal temperature occurs at the expanding detonation products front [1].
In this paper, a comparation of analytical methods for blast wave overpressure is presented and compared to experimental data in Section 3. Also, numerical simulations of free airburst explosions were done in Ansys Autodyn for 1 kg, 10 kg, 10 T of TNT, and obtained data were compared with available empirical data, in Section 4. Section 5 presents conclusions.

Explosion types and blast wave parameters
There are two main groups of blast loads, regarding explosive confinement: unconfined (free air blast explosions), an air burst close to the ground, and surface burst) and confined (fully vented, partially confined, and fully confined). In free air spherical explosions, shock wave moves from the center of the explosion and hits the target with no amplification. Hemispherical explosions are surface explosions where the detonation is close to the ground, and where the incident shockwave is strengthened due to reflections from the ground. As the blast wave propagate, Mach front ( Figure 2) can be formed by superposition of the incident and reflected waves [2]. Confined explosion relates to detonation inside structures. Figure 2: Free-air burst, air burst near ground and surface burst [2] The pressure wave function in the air is schematicaly shown in Figure 3. There is a rapid rise from ambient pressure to a peak incident overpressure Pso which decays to the ambient value in time to (positive phase duration). The front arrives at a location in time tA. Negative phase (duration to-) is characterized by underpressure Pso-. The incident impulse of the blast wave is the integrated area under the pressure-time curve (marked as is for the positive and is-for the negative phase) [2].  [2] In technical manuals [2], the overpressure exponential curve is often replaced by a linear curve. The negative phase is often ignored in structural calculations.
When the blast wave hits the surface, the reflected pressure value can be more than 20 times larger than the incident pressure value [3].

Prediction of blast wavefront parameters
A large number of studies were conducted in order to understand better blast effects after the explosion and the response of structures to blast loads. Empirical, semi-empirical, and numerical methods are mostly used for the prediction of blast effects.
Empirical methods are correlations with test data. The accuracy of empirical formula is usually lower for near-field explosions.
Semi-empirical methods are partly based on physical models. These methods rely also on experimental data and their accuracy is usually better than purely empirical methods. They are used in certain blast codes (programs) Numerical (computational fluid dynamics, CFD) methods are usually the most accurate ones and are based on equations describing basic laws of continuum mechanics (conservation of mass, momentum, and energy). The physical behavior of materials is generaly described by constitutive equations [4].
Some of these programs (ie AUTODYN) can also be used to model the near-field explosion effects, in order to make new formulas for near-field air-blast estimations, and to update old formulas (ie. for reflection coefficients).  Reference [5] includes formulas for air bursts and surface bursts, for estimation of the values for incident and reflected pressures (and for other parameters as well). The positive phase pressures, impulses, durations, and other parameters of shock environment for free air and surface burst (for TNT) are given in Figs. 5 and 6 versus scaled distance Z (from Z=0,05 m/kg 1/3 to Z=40 m/kg 1/3 ).
According to Hopkinson-Cranz law, a dimensional scaled distance Z = R/W 1/3 , where R is distance from the detonation point, and W is the mass of explosive charge (equivalent to TNT mass). The values of TNT equivalent mass factors can be found in papers and technical manuals.
Using diagrams, to obtain absolute value of particular parameter, their scaled value is multiplied by a W 1/3 to take the size of the explosive charge into calculation. Pressure and velocity are not scaled [5].
Calculations for surface burst are more straightforward than for airburst near the ground because of complex wave reflections in the latter case.
High-order polynomial equations functions made using regression analysis of test data [5] are used in the CONWEP.
In CONWEP (utilizes Kingery-Bulmash polynomials), following formula is used for defining P-t curve: Here A is decay coefficient, which can be calculated by curve fitting of experimental P(t) curve in positive phase [4]. Besides Eq (1), formulas of Flynn and Etheridge can be used for prediction of positive blast pressure profile.
Automated formulas for surface blast parameters can also be found in programs BECV4 and A.T. -Blast. Table 1 below shows peak reflected overpressures Pr for surface burst (ie. in Figure 6) with different W-R (mass-distance) combinations. These values can be obtained using polynomial equations, or programs. Experiments have shown that human blast tolerance depends on magnitude of blast wave pressure and the shock duration. Tests showed that the lungs is the critical organ in injuries related to blast waves. The air bubbles release from damaged lung alveoli into the body vascular system accounts for most fatalities. The severe lunghaemorrhage occurs at pressure levels above 5,5 bar, while the lethality due to lung damage is for pressures 6,9 to 8,3 bar, Near 100 per cent deaths are confirmed for pressure levels 13,8 -17,2 bar. Besides these primary effects of blast on humans, there are secondary (fragments from surrounding structures and buildings collapse), tertiary (blast wave and winds that can sweep people) and quaternary effects (other explosion-related injuries) [6].
Blast effects may damage structures by two main type of loading: diffraction and drag loading. Diffraction loading relates to loading the structure from all sides, where the blast overpressure is applied to sides of the object nearly simultaneously (i.e. on a building oriented towards explosion blast wave would arrive on the front sides and roof at nearly the same time). Ductile targets (i.e. made of metal) can be crushed. Brittle targets (i.e. made of concrete) will most likely shatter. The loading can be estimated from the peak overpressure [7].
Drag loading, related to dynamic pressure, is the force which acts on surfaces perpendicular to the pressure wave. The drag load is less than the diffraction loading, and it is applied for a longer time period. The drag load reverses direction, which can tear objects apart. Flexible targets are not damaged significantly by diffraction loading but can be prone to drag loading injuries.
Objects not fixed can be thrown several meters away. Humans are quite vulnerable to this type of loading. Lightly protected equipment can also be damaged by drag loading [7].
Prediction formulas for blast loads on structures can be found in technical manuals. Table 2 shows expected damage on structures loaded with the blast wave. Minor damage in some buildings 7,6 -12,4 Metal panels deformed 12,4 -20 Concrete walls damaged over 34, 5 Wooden buildings demolished 27,6 -48, 3 Major damage to steel objects 41,4 -62, 1 Heavy damage to reinforced buildings 69 -82,7 Probable demolition of most buildings Blast overpressure is the most important parameter in blast effects modelling and many analytical relations can be found in the literature, most of them based on test data at different scaled distances and charge sizes.
The formulas for a free spherical airburst explosion can also be used for surface explosions by increasing the charge mass using the multiplication with a coefficient (reflection factor) 2η, where η takes into account energy used for deformation of the base material, with values ranging from 0,55 for water to 1 for steel.
The formulas for spherical airbursts and hemispherical surface burst blast wave pressure PS, presented by several authors, are given in Tables 3 and 4.  Tables 3 and 4 because of space, but they can be found in their original paper [5]. Constants for formula by Jeon are also omitted from Table 4, and can be found in [9]. .. ..
..   for airburst, obtained using formulas from different models, compared to experimental [5] curve (Kingery-Bulmash polynomials) In this paper we made comparison of blast overpressures (0,15 < Z < 10) for airburst and surface (hemispherical) burst, using different models (presented in Table 3), and results are shown in Figs. 7 and 8. The experimental data [5] are included for reference as the widely approach for blast parameters prediction. Figure 7 shows most curves deviate from the test data for small scaled distances (ie for Z < 0,4). This may be since some of these equations (Brode, Mills) were developed for nuclear blast.

Numerical simulations
Numerical simulations of explosions were done in Ansys AUTODYN (version 2019), engineering software package that use finite difference, finite volume, and finite element techniques to solve a wide variety of non-linear problems in solid, fluid, and gas dynamics. This type of program is sometimes referred to as a hydrocode. The phenomena to be studied with such a program can be characterized as time-dependent with geometric nonlinearities (large strains and deformations) and material non-linearities (plasticity, failure, strain-hardening and softening, multiphase equations of state).
When simulating explosions, material properties can be selected from AUTODYN library (Table 5). Air uses Ideal gas equation of state, where the pressure P is related to the energy e (with adiabatic exponent ) as: This form of an equation is useful for its simplicity and computation ease, where only the value of  needs to be supplied.
The values of the constants A, R1, B, R2 and  from Eq. (40) for some explosives have been estimated from experiments (cylinder test -expansion of a hollow metal cylinder filled with explosive), and are available in AUTODYN material library. AUTODYN (1D, 2D, 3D) has been verified for use in blast effects estimation in many research papers.
In this research, we give an example of numerical simulation procedure of blast wave formation in the air after the detonation of 1 kg, 10 kg, and 10 T of TNT. Default values for the TNT and air equations of state were used, as well as other parameters specified in AUTODYN library (Table 5).  Zones filled with the air were given initial energy of 2,06810 5 mJ/mm 3 [13] to provide an ambient pressure of 101,3 kPa. A detonation point is located at the center of explosive (0,0,0) to start the explosion at time zero.
Gauge points (for pressure reading, Figure 10  In AUTODYN, the programs' Euler solver was used in a wedge-shaped grid (1D model) at the apex of which the explosive charge (TNT) was located. For a mass of 1 kg, the explosive radius was 52,7 mm, for 10 kg TNT this radius was 113,56 mm, for 10 T it was 1135,6 mm. Outflow transmitting boundary condition was used to eliminate the wave reflection [13]. In reference [13], cell size values of 3 mm were recommended as a compromise between simulation duration and accuracy, but the results are dependent on the mesh so we used a cell size of 1 mm for comparison with the test data. Reference [1] suggests cell size equal to 0,002 times the distance to the monitoring location on the reflecting surface provided results within 10% of the converged value for 0,0553 ≤ Z < 40 m/kg 1/3 . The reference [13] suggests quadratic viscosity values of 0,1 which we adopted in our simulations.
As a part of the mesh size sensibility study, we first numerically simulated the blast wave formation after the explosion, with three different mesh sizes (1 mm, 10 mm, and 20 mm, quadrilateral cell size) for 10 kg TNT charge, as an example. Results can be seen in Table 6, and in Fig 11. Table 6 shows that for coarser grid (20 mm) results deviate substantially from experimental data comparing to finer mesh (1 mm).  Figure 11 shows mesh size influence on shock wave form at a distance of 1m, after detonation of 10 kg TNT. The coarser the grid, the slower the rate of rise of pressure front and the flatter the curve shape. Peak overpressure values are different for different cell sizes (smaller cell sizes giving larger pressure values). Figure 11: Mesh size influence on shock wave form at a distance of 1m, for 10 kg TNT explosive charge Figure 12 shows the motion of a pressure wave along the mesh during simulation, for the case of 1 kg TNT.
During the simulation, AUTODYN shows rarefaction wave going toward detonation point (explosive charge). It reduces the pressures and density behind the front of the expanding detonation products.
This rarefaction wave, moving inward, is then reflected outward from the centre of the explosive forming, so called, secondary shock wave. Formation of a secondary shock wave, after the detonation of 10 kg explosive charge (TNT) can be seen in Figure 13.   Figure 15 illustrates the wave form P(t), obtained in simulations, at different distances from an explosion, for 1 kg and 10 kg TNT charges.
The second peak on the first (black) curve is likely caused by the blast wave reflection at the material interface (TNT/Air) because of different acoustic impedances of these materials (it should be noted that this second peak is of little practical significance [13]).
We can see from Figure 15 that blast wave pressures drop significantly with distances in air.
Generally, when the blast wave moves through the atmosphere, pressures are rapidly decreased and have a brief existence span, (measured in milliseconds; as shown clearly in Figure 15). The results from numerical simulations for 1 kg, 10 kg, and 10 T of TNT explosive charge were compared with experimental data [5], and relative differences presented in Table 7.
Relative differences observed between simulation results and test data were smaller than 5,9 % for all cases for 1 kg TNT, smaller than 4 % for 10 kg TNT, and less than 2,1 % for 10 T charge.
AUTODYN User's Manual reports of errors around 4% for peak incident overpresure between numerical simulations and test data. As can be seen from the results, AUTODYN can successfully model blast wave formation in the air after the detonation of high explosive. Results from these simulations can be remapped into 3D urban scenario simulation of explosion, which can be follow-up research.

Conclusions
In this paper we made a comparation of available formulas for blast wave overpressure. These formulas were compared with available empirical data [5]. The Kinney and Shin models showed the best agreement with experimental data for free airburst, while for surface burst, Swisdak, Vanuci, and Jeon models predicted test data most accurately.
Also, we performed a numerical simulation of airburst detonation after the explosion of TNT charge in Ansys AUTODYN, for the cases of 1 kg, 10 kg, and 10 T of TNT explosive charge, with a fine mesh (1mm), where obtained data were compared with experimental data (Kingery & Bulmash). Relative differences observed between simulation results and test data were smaller than 5,9 % for all cases for 1 kg TNT, smaller than 4 % for 10 kg TNT, and less than 2,1 % for 10 T charge. In this part, we also gave a detailed description of the procedure for these simulations as a valuable tool for blast wave phenomena investigation.
The novelty in the paper is that we introduced new exponential and power functions for surface blast overpressure estimation, with small relative differences compared to experimental data [5]. These formulas, with acceptable accuracy, are somewhat simpler than many of the formulas used for blast overpressure estimation, and can be implemented faster.
The following work in this field can be pointed towards complex geometry urban scenario blast effect modeling, where surface burst effects dominate, with complex reflections of shock wave present.