Mathematical Problems in Engineering

Volume 2018, Article ID 2530239, 17 pages

https://doi.org/10.1155/2018/2530239

## Numerical Method for Predicting the Blast Wave in Partially Confined Chamber

^{1}Key Laboratory of High Performance Ship Technology (Wuhan University of Technology), Ministry of Education, Wuhan 430063, China^{2}Departments of Naval Architecture, Ocean and Structural Engineering, School of Transportation, Wuhan University of Technology, Wuhan 430063, China

Correspondence should be addressed to Xiang-shao Kong; nc.ude.tuhw@sxgnok

Received 27 August 2017; Accepted 22 January 2018; Published 27 February 2018

Academic Editor: Michele Betti

Copyright © 2018 Wei-zheng Xu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.