Numerical Method for Predicting the Blast Wave in Partially Confined Chamber

Blast waves generated by cylindrical TNT explosives in partially confined chamber were studied numerically and experimentally. Based on the classical fifth-order weighted essentially nonoscillatory finite difference schemes (fifth-order WENO schemes), the 1D, 2D, and 3D codes for predicting the evolution of shock waves were developed. A variety of benchmark-test problems, including shock tube problem, interacting blast wave, shock entropy wave interaction, and double Mach reflection, were studied. Experimental tests of explosion events in a partially confined chamber were conducted.Then, the 3D code was employed to predict the overpressure-time histories of certain points of chamber walls. Through comparing, a good agreement between numerical prediction and experimental results was achieved. The studies in this paper provide a reliable means to predict the blast load in confined space.


Introduction
A blast wave from an explosion is the rapid release of energy from energetic materials in a very short time and provokes a sharp peak of pressure.A supersonic wave is subsequently generated and propagates away from the detonation origin.The consequence of blast accidents usually means severe damage of structures and loss of life.And even worse, explosions which occurred in confined space may be extremely severe and cause more damage than a similar external free-air explosion due to the aggravation effect.Researches of internal explosion events usually aim for predicting both the blast loading acting on the inner space boundary and the dynamic responses of structural elements, which are essential to the design of protective facilities [1][2][3].
The main feature of a blast wave is the presence of discontinuities, such as shocks and contact discontinuities.Actually, the evolution of a blast wave is a convection dominated problem, which can be described by the hyperbolic conservation laws.For decades, researchers have been contriving to describe the profile of blast wave numerically.Feldgun et al. [4,5] used AUTODYN commercial program implementing the Eulerian multimaterial approach to investigate characteristics of an interior explosion with or without venting.Baum et al. [6] developed a coupled computational fluid dynamics (CFD) and computational structural dynamics (CSD) methodology to simulate the blast waves generated by explosives in a test facility with rigid and deformable walls.Benselama et al. [7] used a secondorder accurate finite difference scheme to simulate the blast waves generated by the TNT explosives in free air using the axial symmetry Euler equations.The predicted results were in good agreement with the experimental data.Alpman [8] developed an in-house computational fluid dynamics code using Euler equation and adaptive grids to simulate the one-dimensional spherical blast wave in the open space, and the explosive was modeled as a high pressure sphere with JWL equation of state.Alpman [9] performed onedimensional spherical blast wave simulations in the open space using Runge-Kutta Discontinuous Galerkin Method and the explosive was modeled as a high pressure sphere with JWL equation of state.The simulations of blast wave evolution were conducted by solving Euler equations using a finite volume method, in which second-order accuracy was achieved by extrapolating density, velocity, and pressure at the cell interfaces.
In order to explore the actual characteristics of confined explosion and validate the accuracy of numerical methods, experimental tests were carried out by researchers [10][11][12][13][14].
In the present work, the 1D, 2D, and 3D codes of predicting the evolution of shock waves were developed by employing the fifth-order weighted essentially nonoscillatory finite difference (WENO) schemes.The planar shock tube and double Mach refection were simulated to preliminary verify the reliability of the codes.Besides, experimental tests of partially confined explosions were carried out to validate the 3D code.This paper is organized as follows.In Section 2, the three-dimensional Euler equations and dimension split method are presented.The basic theory of the WENO schemes and numerical discretization of the Euler equations are described in Section 3 in detail.One-dimensional planar shock tube problems and two-dimensional double Mach reflection example are presented to validate the precision and reliability of the numerical method implemented in the codes in Section 4. Experimental designs are presented in Section 5.In 3D numerical simulations, comparisons with experiments are presented in Section 6.Finally, the conclusions are drawn in Section 7.

Governing Equations and Dimension Split Method
2.1.Governing Equations.The blast wave is the consequence of detonated explosive.In the numerical model, the detonation products of explosives are usually treated as statically pressurized gas, which is specified within certain space occupying the volume of the explosive or slightly larger.The pressure of the gas is set to be the instantaneous detonation pressure and the specific internal energy is set to be the detonation heat of the charge.Afterwards, the expansion of the pressurized gas produces the blast waves traveling outwards with an inward rarefaction counterpart.Euler equations are widely used to model many interesting physical phenomena related to propagation of blast wave.The general formulation of conservative Navier-Stokes equations for the model described above reads as where  is the conservative variables vector and (), (), and () are the fluxes in , , and  directions, respectively. is the source term vector that may include the viscous and thermal diffusion energy terms and can be defined as follows: where   ,   are mass and energy source terms and   ,  V , and   are momentum source terms in , , and  directions, respectively.They may include physical dissipative terms, such as viscous and thermal effects.But with regard to the issue of blast wave propagation, the physical dissipative terms are fairly small compared to convective terms which are dominant, and they are not taken into consideration in the present calculation.Thus, (2) is simplified as threedimensional Euler equations: where  is the fluid density and , V, and  are fluid velocity in , , and  directions, respectively. is the static pressure,  is total energy, and  is total energy per unit mass.The ideal gas law is employed as equation of state for calorically perfect gas; it is given as follow: where  is the ratio of specific heat and constant value of 1.4 is defined in the present calculation.

Dimension Split
In each direction, (7) or (8) or (9) can be considered as one-dimensional hyperbolic equations, respectively, as follows:

Numerical Methods
3.1.Fifth-Order WENO Schemes.Consider an uniform grid defined by the points   = Δ,  = 0, . . ., , which are also called cell centers, with cell boundaries given by  +1/2 =   + Δ/2, where Δ is the uniform grid spacing.The semidiscrete form of (10) by the method of lines yields a system of ordinary differential equations [15], where   () is a numerical approximation to the point value (  , ).A conservative finite difference formulation for hyperbolic conservation laws requires high-order consistent numerical fluxes at the cell boundaries in order to form the flux differences across the uniformly spaced cells.The conservative property of the spatial discretization is obtained by implicitly defining the numerical flux function ℎ() as such that the spatial derivative in (11) is exactly approximated by a conservative finite difference formula at the cell boundaries, where ℎ ±1/2 = ℎ( ±1/2 ).
High-order polynomial interpolations to ℎ ±1/2 are computed by using known grid values of  = (  ).The classical fifth-order WENO scheme uses a 5-point stencil is the third-degree polynomial below, defined in each one of the stencils. where are Lagrangian interpolation coefficients.The expanded form of ( 15) is as follows: The weights   are defined as The coefficients  0 = 1/10,  1 = 3/5, and  2 = 3/10 are the ideal weights since they generate the central upstream fifth-order scheme for the 5-point stencil  5 .  is referred to as the nonlinear weights.The parameter  is used to avoid the division by zero in the denominator and  = 2 is chosen to increase the difference of scales of distinct weights at nonsmooth parts of the solution.
The smoothness indicators   measure the regularity of the polynomial approximation   (  ) at the stencil   and are given by The expressions of   in terms of the known grid values in each 3-point stencil are given by

Lax-Friedrichs Flux Splitting Method.
For the sake of stability, the finite difference procedure described above is applied to  + () and  − () separately, where  ± () correspond to a flux splitting, with In order to satisfy the up-winding regulation, a biased stencil with one more point to the left and to the right of  +1/2 was used to reconstruct  + () and  − (), respectively.It is further required that  ± () are smooth functions of .The most commonly used flux splitting is the Lax-Friedrichs splitting [16], with From ( 22), it could be obtained that 3.3.Characteristic Matrix Splitting.In order to obtain better results, the characteristic matrix is split to consider the upwind characteristics of the Euler equations.Here, onedimensional Euler equations are taken as an example.It is defined that  = ()/; then the left and right eigenvector matrices,  and , are given in ( 25) and (26), respectively.
where ( Then one-dimensional Euler equations are presented as follows: The diagonal matrix Λ can be split into two parts Λ = Λ + +Λ − ; here the Lax-Friedrichs flux splitting method is used as described in Section 3.2.Then (28) can be derived After the convection flux was determined by using the fifth-order WENO schemes presented in Section 3.1 and the results are recorded in  ± , then the results should be projected to the original physical space through  =  ± .Finally, the whole one-dimensional Euler equations are updated in each time step.

Numerical Discretization in Time.
In the developed codes, the ODEs such as (11) resulting from the semidiscrete form of the PDEs such as (10) are evolved in time by the thirdorder total variation diminishing Runge-Kutta scheme (RK-TVD): Δ ( (2) ) . (30)

Numerical Examples
In order to verify the developed codes in the present work, three classical Euler problems and double Mach reflection problem are employed to demonstrate the resolution of the code in one-and two-dimensional space.computational domain is set to be 1 and the initial condition is defined as follows: (1, 0, 1) , 0 ≤  < 0.5 (0.125, 0, 0.1) , 0.5 ≤  ≤ 1.
The initial contact surface is located at 0.5, and the boundary condition is set to be outflow.The total grid number is 200 and the final time is 0.18.The simulation results are compared with the exact solution.The density curve and the pressure curve are shown in Figures 1(a) and 1(b), respectively.The calculated results indicate that the code based on the fifth-order WENO schemes is capable of capturing the shock wave, rarefaction wave, and the contact wave.

Interacting Blast
Wave.This problem was firstly discussed by Woodward and Colella [18].It refers to the interaction between the shock wave, contact wave, and rarefaction wave.The length of computational domain is set to be 1 and the initial condition is defined as follows: (1, 0, 0.01) , 0.1 ≤  ≤ 0.9 (1, 0, 100) , 0.9 ≤  ≤ 1. ( The boundary condition is set to be reflection at both ends of the computational domain.The total grid number is 400 and the final time is 0.038.The simulation results are compared with the exact solution.The density curve and the pressure curve are shown in Figures 2(a) and 2(b), respectively.The simulation results indicate that the code based on the fifth-order WENO schemes can be used to describe the interaction of the complicated waves.

Shock Entropy Wave Interaction.
This problem [19] describes the interaction between shock wave and flow field with small density fluctuations.The length of computational domain is set to be 1 and the initial condition is defined as follows: (, , ) = { (3.857143, 2.629369, 10.33333) , −5 ≤  < −4 (1 + 0.2 sin 5, 0, 1) , −4 ≤  ≤ 5. ( The boundary condition is set to be outflow.The total number of grids is 400 and the final time is 1.8.The simulation results are compared with the exact solution.The density and the pressure distributions are shown in Figures 3(a) and 3(b), respectively.The simulation results indicate that the code based on the fifth-order WENO schemes is capable of describing the interactions between the shock wave and entropy wave.

Double Mach Reflection.
This problem was firstly proposed and studied in detail by Woodward and Colella [18].The computational domain is 4 × 1 as shown in Figure 4, and the reflecting wall lies at the bottom of the computational domain for 1/6 <  < 4. A right-moving shock with Mach number of 10 is initially positioned at  = 1/6,  = 0 and makes a 60 ∘ angle with the -axis, which is represented by a black line in Figure 4.The propagation direction of the shock wave is described by the arrow.The exact postshock condition is imposed for the bottom part from  = 0 to  = 1/6, and reflective boundary condition is used for the rest.At the top boundary, the flow values are set to be the exact solution of the oblique shock.Inflow and outflow boundary conditions are defined for the left and right boundaries.The unshocked   We have performed a numerical study on the resolution of the fifth-order WENO finite difference methods implemented in the codes developed in the present work on problems containing both discontinuities and complex solution structures.The simulation results show that the method proposed in this paper is capable of predicting the evolution of shock waves.The work in Section 4 encourages us to implement this method for the research of blast waves in partially confined chamber.number 1 and number 2 were installed on wall A, as shown in Figure 8(a).The PCB pressure transducer number 3 was installed on wall B, as shown in Figure 8(b).Two PCB pressure transducers number 4 and number 5 were installed on wall C, as shown in Figure 8(c).The HBM data acquisition system is used to acquire the transient overpressure data with the sampling frequency of 10 6 /s.

Test Cases. Three different doses of cylindrical TNT
explosives with the density of 1630 kg/m 3 were used in experimental tests, as shown in Table 1.Four typical explosion cases were conducted in the chamber shown in Table 2, and each case was repeated at least three times.

3D WENO Simulation and Validation
Through the experiments, the overpressure histories of the gauging points on the walls of the chamber were obtained.The experimental data set would enable the validation of the 3D WENO code developed in this paper.

Explosive Modeling.
In the numerical calculation, a simplified model was used for introducing explosion effects to the flow field.Firstly, it was assumed that the TNT explosive detonated instantaneously, and the explosive was modeled as an isobaric high pressure cylinder.Initial parameters of high pressure detonation gas and air are listed in Table 3.Secondly, the volume of the cylindrical TNT explosives occupied few grids, which will cause serious numerical oscillation at the beginning of the calculation.In order to eliminate this phenomenon, the size of cylinder was expanded twice based on the principle of equivalent energy.When the radius  and height of the cylinder ℎ are given, the density  and initial pressure  of the expanded cylinder can be obtained by Three different doses of cylindrical TNT explosives were simulated, and all expanded to be twice size of the original.The parameters of the expanded cylindrical TNT explosives used in the numerical simulations are listed in Table 4.

Initial Condition and Boundary Condition.
The initial state before the explosion of cylindrical TNT explosives in the whole chamber is shown in Figure 9(a).The cylinder volume with high pressure is located at the center of the chamber.Four walls of the chamber are all set to be reflective boundary conditions, except for the venting hole on wall A with outflow boundary condition.In the present calculation, the grid number is set to be 180 × 80 × 80 to get balance between accuracy and consuming time.In order to observe the evolution of shock wave directly, five slices along the length direction (direction) of the chamber are selected as shown in Figure 9  ×1 = 0 mm, ×2 = 900 mm, ×3 = 1350 mm, ×4 = 1650 mm, and ×5 = 1800 mm.The slices ×1 and ×5 are located at the two end walls of the chamber to describe the phenomena that blast waves converge and reflect at the corner of the chamber.The slice ×2 located at the middle of the chamber is to show the expansion of the high pressure cylinder.The slice ×3 is used to describe the pressure change at the gauging points numbers 1, 3, and 4. The slice ×4 is used to describe the pressure change of the gauging point numbers 2 and 5.

Pressure Field.
In order to show the evolution of pressure field of the whole chamber, four typical simulation results of the Case 4 are presented in Figure 10.And the unit of pressure is Pa.From Figure 10(a), it can be seen that before the blast waves reach the walls of the chamber, the flow field is a three-dimensional cylindrically symmetric flow as the blast in the open air.When the blast waves reach the walls for the first time, reflection occurs as shown in Figure 10(b) and the pressure field on wall A is different from that on walls B, C, and D due to the influence of venting hole.In Figure 10(c), the blast waves reach the two end walls of the chamber along the length direction of the chamber for the first time and reflect from the walls.Figure 10(d) showed that the blast waves converge and reflect at the corner of the chamber for the first time.

Overpressure Histories.
Overpressure histories of gauging points of different simulation cases are presented in Figures 11,12,13, and 14, respectively.It is revealed that the simulated overpressure histories of numbers 1 and 4 are almost the same.The first overpressure peaks of numbers 1 and 4 are higher than that of number 3. The arrival times corresponding to the first peak of numbers 1 and 4 are earlier than that of number 3.These phenomena can be explained by analyzing the detonation process of the cylindrical TNT explosives.It is assumed that the cylindrical TNT explosive is instantaneously detonated at the very beginning of the simulation.The blast waves induced by the detonation gas are schematically shown in Figure 15.Due to the symmetric characteristics of blast wave, the pressure at blast wave front A is the same as that at the blast wave front C along the axial direction of explosives, and the pressure at blast wave front B is the same as that at blast wave front D along the radial direction of explosives [14].Gauging points numbers 1 and 4 are symmetrical about the explosives and presents almost the same overpressure histories.From Figure 10(a), it can be    seen that the shape of high pressure region transforms from cylinder into ellipsoid gradually, and the longer axis of the ellipsoid is located in the radial direction of the explosive, which indicated that the radial energy is higher than axial energy produced by the explosive.As the intensity of shock wave is proportional to the explosion energy, the intensity of shock wave in the radial direction is higher than that in the axial direction, and the propagation velocities of the blast waves along the radial direction are faster than those along the axial direction.Moreover, with the increase of the TNT charge mass, the differences of first peak value and arriving time between condition except for the different diameter of the venting hole on wall A of the chamber.It seems that the size of venting hole has little influences on the overpressure histories at early stage of the explosion.From Figures 16 and 17, it can be seen that the numerical results agree with the experimental results.From Figure 18, the numerical results are well consistent with the experimental data for gauging points numbers 1, 2, and 5, except for the value of the first overpressure peak of gauging points number 3 in Case 3.This great difference is because of the higher recorded first overpressure peak of gauging point number 3 in Case 3.  and comparison, it is concluded that the pressure transducer fails to capture the first overpressure peak of gauging point number 3 in Case 4.
The simulated overpressure peaks are slightly smaller than the experiment data and the arrival times of the shock waves are not in full accord with the experiments in all the four cases.Because the physical viscous effects and thermal conduction are not taken into consideration in the numerical code and the detonation gas is ideally assumed to follow perfect gas state.To some extent, these factors may reduce the intensity of the blast waves.Besides, there exist uncertain factors in the detonation process of explosives in the experimental tests.On the whole, the predicted shock waves match well with measured waves from the perspective of initial shock arrival time, peak pressure, and waveform and subsequent shock arrival times, peak pressures, and waveforms.

Conclusions
Blast waves generated by cylindrical TNT explosives in partially confined chamber have been studied numerically and experimentally.And the following conclusions can be drawn: (1) The in-house 3D WENO code developed in this paper was validated to be practicable in simulating the evolution of blat waves generated by cylindrical TNT explosives in partially confined space.

5. 1 .
Experimental Setup.The experimental device of a steel chamber with a venting hole is shown in Figure 6(a), and the size of the cuboid-shaped steel chamber is 1800 × 800 × Explosives Experimental chamber Venting hole

Figure 10 :
Figure 10: Pressure field of the whole chamber.

6. 4 .
Validation of the Numerical Method.Comparisons between numerical and experimental results of the four cases are shown in Figures 16, 17, 18, and 19.

Figure 18 :
Figure 18: Comparisons between numerical and experimental results of Case 3.
1 , and  = ( + )/.The diagonal matrix of  is Λ.The relation of , , , and Λ can be expressed as

Table 1 :
Parameters of the cylindrical TNT explosives.

Table 3 :
Initial parameters of detonation gas and air.

Table 4 :
Parameters of the expanded cylindrical TNT explosives.