Abstract
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.
1. 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–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 second-order 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 one-dimensional 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–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.
2. 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 aswhere 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 , , 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 three-dimensional Euler equations:, , , and are defined bywhere is the fluid density and , , 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.
2.2. Dimension Split Method
Generally, the dimension split method is used to solve the multidimensional Euler equations. In the present study, Strang’s operator splitting method is employed, and (4) is split, respectively, into , , and directions as follows:
In each direction, (7) or (8) or (9) can be considered as one-dimensional hyperbolic equations, respectively, as follows:
3. Numerical Methods
3.1. Fifth-Order WENO Schemes
Consider an uniform grid defined by the points , , which are also called cell centers, with cell boundaries given by , 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 assuch that the spatial derivative in (11) is exactly approximated by a conservative finite difference formula at the cell boundaries,where .
High-order polynomial interpolations to are computed by using known grid values of. The classical fifth-order WENO scheme uses a 5-point stencil , which is subdivided into 3-point stencils : , , and . The fifth-order polynomial approximation is built through the convex combination of the interpolated values , in which 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 , , and are the ideal weights since they generate the central upstream fifth-order scheme for the 5-point stencil . is referred to as the nonlinear weights. The parameter is used to avoid the division by zero in the denominator and 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
3.2. 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 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],withFrom (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, one-dimensional 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 , , , , and . The diagonal matrix of is . The relation of , , , and can be expressed as follows: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.
3.4. 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 third-order total variation diminishing Runge-Kutta scheme (RK-TVD):
4. 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.
4.1. One-Dimensional Simulation Examples
4.1.1. Riemann Problem of Sod
Sod shock tube problem [17] is the typical case to verify numerical methods. It contains rarefaction wave, shock wave, and contact wave. The length of computational domain is set to be 1 and the initial condition is defined as follows:
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.
(a) Density curve at
(b) Pressure curve at
4.1.2. 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:
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.
(a) Density curve at
(b) Pressure curve at
4.1.3. 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:
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.
(a) Density curve at
(b) Pressure curve at
4.2. Two-Dimensional Simulation Example
4.2.1. Double Mach Reflection
This problem was firstly proposed and studied in detail by Woodward and Colella [18]. The computational domain is as shown in Figure 4, and the reflecting wall lies at the bottom of the computational domain for . A right-moving shock with Mach number of 10 is initially positioned at , 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 to , 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 gas has the density of 1.4 and the pressure of 1. The problem was run till and the grid number is set to be . The ratio of specific heats . The results in the region of are displayed.
The density contours with 40 equally spaced contour lines from to are drawn in Figure 5(a). A “blown-up” portion around the double Mach is shown in Figure 5(b) to exhibit the significant resolution structures.
(a) The whole region
(b) Blown-up region
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.
5. Experimental Designs
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 mm. What is more, the radius of the venting hole is within the range of 50 mm to 100 mm. The walls of chamber are marked in Figure 6(b). Two ends of the chamber are sealed with thick plates. In the experimental tests, the charge was hung at the center of the chamber, and the cylindrical TNT explosives were detonated at one end close to the wall D of the chamber, as shown in Figure 7. Two PCB pressure transducers 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 /s.
(a) Experiment chamber
(b) Wall numbers of experiment equipment
(a) On wall A
(b) On wall B
(c) On wall C
5.2. Test Cases
Three different doses of cylindrical TNT explosives with the density of 1630 kg/m3 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.
6. 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.
6.1. 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.
6.2. 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 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(b). The locations of the slices are ×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.
(a) The chamber and initial cylinder
(b) Initial slices of the chamber
6.3. Simulation Results
6.3.1. 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.
(a) ms
(b) ms
(c) ms
(d) ms
6.3.2. 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.
(a) Numbers 1, 3, and 4
(b) Numbers 2 and 5
(a) Numbers 1, 3, and 4
(b) Numbers 2 and 5
(a) Numbers 1, 3, and 4
(b) Numbers 2 and 5
(a) Numbers 1, 3, and 4
(b) Numbers 2 and 5
(a) In Cases 2, 3, and 4
(b) In Case 1
Moreover, with the increase of the TNT charge mass, the differences of first peak value and arriving time between gauging points numbers 1, 4, and 3 become more apparent. This is mainly due to the larger ratios between the height and diameter of cylindrical explosive, of which the ratios are listed as 1.06, 1.35, and 1.26 corresponding to the charge mass of 55 g, 110 g, and 200 g. The blast waves induced by the detonation gas are shown in Figure 15(b). Blast wave fronts A and C along the axial direction of the cylindrical TNT explosives are approximately the same as the blast wave fronts B and D along the radial direction. Thus, the calculated pressure histories of gauging points 1, 3, and 4 are almost the same as shown in Figure 11(a).
From Figures 11(b), 12(b), 13(b), and 14(b), it can be found that the overpressure histories of gauging point numbers 2 and 5 are almost the same. The reason is that these two points are symmetrical to the position of explosives.
From Figures 12 and 13, it can be found that the overpressure histories of all gauging points in Case 2 are almost the same as those in Case 3, which share the same initial 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.
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.
(a) Gauging point number 2
(b) Gauging point number 3
(c) Gauging point number 4
(d) Gauging point number 5
(a) Gauging point number 1
(b) Gauging point number 2
(c) Gauging point number 3
(d) Gauging point number 5
(a) Gauging point number 1
(b) Gauging point number 2
(c) Gauging point number 3
(d) Gauging point number 5
(a) Gauging point number 2
(b) Gauging point number 3
(c) Gauging point number 4
(d) Gauging point number 5
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.
From Figures 17(a) and 17(c), it can be found that the recorded first overpressure peak of gauging point number 1 is larger than that of gauging point number 3 in Case 2, which is also close to numerical result presented in Section 6.3.2. The recorded first overpressure peak of gauging point number 1 is close to that of gauging point number 3 in Case 3 as shown in Figures 18(a) and 18(c). However, it is inconsistent with the numerical results. The reason is mainly due to the influence of uncertainty in the detonation process of the cylindrical TNT explosives in the experiments. In the experiment, the cylindrical TNT explosives were detonated at one end of the cylinder close to the wall D of the chamber (Figure 7). The blast waves induced by the detonation gases in the experiment are schematically shown in Figure 20. During the detonation process, nonuniform detonation wave would yield in the explosive and thus has consequentially effect on the subsequent pressure field, causing almost the same intensity of the blast wave fronts A, B, and D in Figure 20. When the explosives are detonated ideally, the intensity of the blast wave fronts B and D is higher than that of blast wave front A. As shown in Figure 20(a), the recorded first overpressure peak of gauging point number 1 is larger than that of number 3 in Case 2.
(a)
(b)
From Figure 19, it can be seen that the simulation results agree well with the experiment results in Case 4 except for the gauging point number 3. The recorded first overpressure peak is much smaller than the expected value, which matches the conditions in other experimental cases. Based on analyzing 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.
7. 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.
(2) The assumption of the detonation gases obeying the perfect gas state equation can be used to simulate the evolution of blast waves generated by the TNT explosives in partially confined chamber. The ratio of height to diameter of the cylindrical TNT explosives has great influences on the first overpressure histories near the explosion source and has little influences on the overpressure histories far from the explosion source.
(3) From the experimental and numerical results, it is found that the area of the venting hole has little influences on the overpressure histories at its early stage. Overall, the studies in this paper provide a base to conduct future researches on the blast waves in partially confined chambers.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge the support of National Defense Fundamental Research Project (B1420133057), Equipped with the Joint Fund of the Ministry of Education (6141A020331), National Natural Science Foundation of China (51409202 and 11502180), and the Fundamental Research Funds for the Central Universities (2016-YB-016 and 2015-YB-005). The authors wish to express gratitude to Dr. Hou Hai-liang for his valuable help during experimental tests.