Comparative study of numerical schemes for strong shock simulation using the Euler equations

A numerical study of extremely strong shocks was presented. Various types of numerical schemes with first-order accuracy and higherorder accuracy with adaptive stencils were implemented to solve the one and twodimensional Euler equations based on the explicit finite difference method, including Roe’s first-order upwind, Steger-Warming Flux Vector splitting (FVS), Sweby’s flux-limited and Essentially Non-oscillatory (ENO) scheme. The result comparisons were carried out to discuss which scheme is the most suitable for strong shock problem. The dissipative nature of the firstorder scheme can be easily seen from the numerical solutions. High order ENO scheme had the best resolution for the case having weak discontinuity, but it overpredicted the shock wave location for the case of strong discontinuity.


INTRODUCTION
In physics, shock waves are small transition layers of abrupt change of physical states such as density, pressure or temperature. For engineering problems, shock waves are often observed in dam break problem or from an explosion. In order to predict and evaluate the effects of strong shocks, computational approach is a robust and low cost method to investigate the nature of the shock waves. This study concerned strong shock, e.g. shock from an explosion, which has large discontinuities. Many robust, stable and accurate numerical methods have been developed for shock-capturing problem including some basic methods such as Lax-Friedrichs's method, Lax-Wendroff's method, MacCormark's method Godunov's method and some modern methods such as Flux-limited method, Flux-corrected method. While the basic methods have the linear numerical dissipation and distributed evenly for all grid points, which makes it not possible to capture strong shock, whereas the modern methods can adaptively distribute the non-linear numerical dissipation to each grid point, thus the shock waves can be moderately captured. In the past, numerical simulation involving strong shock has been documented by many researchers.
Collela and Glaz [16] have performed solution of complicated flow field by using the Riemann problem for gasdynamics. D. Kh. Ofengeim, D. Drikakis [17] have studied the propagation of plane blast wave over a cylinder by solving the Euler equations and the Navier-Stokes equations with an adaptive-grid method and a second-order Godunov scheme. Emre Alpman [6] has simulated spherical blast wave by using the Euler equations and adaptive grid. In [8], the blast wave is also simulated by the Runge-Kutta Discontinuous Galerkin Method. This study presented the implementation of some typical finite difference schemes solving the compressible Euler equations. As categorized in [1], the numerical methods will be studied are the Flux Difference scheme, whose Roe's Approximate Riemann Solver [9] is the representative; the Flux Vector Splitting scheme, whose Steger-Warming method [10] is the representative; the Flux Limited Method, whose TVD Sweby's method [11] is the representative and the Flux-Corrected Method whose Essentially Non-Oscillatory Scheme [12], [13] is the representative.

THEORY BACKGROUND
The hyperbolic partial differential equations describe many physical phenomena such as wave, heat, fluid flows, elasticity, etc. In fluid dynamics, the compressible flow is governed by the Euler equations. These equations are simplified from the Navier-Stokes equations by neglecting the effects of viscosity. The Euler equations are generally presented as a system of hyperbolic conservation laws. Due to the mathematical properties of the system of hyperbolic conservation laws, the solutions consist of waves traveling with the characteristic speed. By casting in the conservation form, the Euler equations allow shock waves or discontinuities be part of the solutions.

Governing Equations
The two-dimensional Euler equations can be written in conservation form as follows [2] , t x y where u are the conserved variables and f(u) and g(u) are conservative fluxes in x and y direction, respectively given by Here ρ is the density, p is the pressure, u and v are the respective velocity components in x and y direction and E is the total energy per unit volume and the total enthalpy is defined as

Solution Method
In order to solve the multidimensional Euler equations, we apply the dimensional splitting technique [2] by solving two consecutive extended one-dimensional problems with extra velocity components.
Based on the dimensional splitting method, all numerical methods will be presented in onedimensional form for the x-split Euler equations. For the calculation of y-split, the role of the two velocity components are interchanged. Consider the integral form of the x-split Euler equations of 1 , The equation (7) can be written in the form with approximate solution operator X as

Roe's Approximate Riemann Solver
The idea was first introduced by Godunov [4], the numerical fluxes at cell interface are evaluated by using the solution of the Riemann problem. The Roe's scheme locally linearize the nonlinear flux function. In the control volume, the Jacobian matrix of the Euler equations is replaced by the Roe-averaged matrix, which is based on the left and the right states. The inner-cell fluxes for the xsplit two-dimensional Euler equations are evaluated as [1,9]       with the Roe-averaged quantities are calculated as follows The wavelength The eigenvalues i λ of the Roe-averaged matrix are given as Due to the linear approximation of the flux, the original Roe scheme is not accurate at sonic points. In order to overcome this problem, Harten's entropy correction is implemented [1,19]

First-order Upwind Flux Vector Splitting
Another approach of the upwind method is flux vector splitting (FVS), which is based on the special property of the Euler equations, namely the homogeneity property. The numerical flux is decomposed into left-running and right-running wave The homogeneity property read With the assumption of hyperbolicity, the Jacobian matrix may be expressed as 1 , (19) where Λ is the diagonal matrix formed by the eigenvalues of the Jacobian matrix A. The column of matrix K is the right eigenvectors of matrix A.
For the x-split Euler equations, the eigenvalues and the eigenvectors read Therefore, the flux vector splitting can be related to the wave splitting by According to Steger and Warming's wave splitting scheme [10], the eigenvalues are defined as follows And the numerical flux is given by [2] 1 2 4

Flux Limited Method-Sweby's TVD Scheme
A high resolution scheme with the nonlinear stability condition of total variation diminishing (TVD) is considered. The idea is using two complementary numerical methods, one near shock and a different one in the smooth region. For a scalar conservation law, the numerical flux at the cell interface was evaluated by (1) (2) (1) The above equation can be interpreted as an overly dissipative flux (1) 1 2 i f  plus a limited amount of an antidissipative flux (2) (1) is the numerical flux calculated by Roe scheme as in (10). The right eigenvector the and the wave strengths ∆v were evaluated as in (11), (13) and (14). The ratio flux difference is defined as In this study, the superbee limiter function [1,11] is employed. The role of limiter function is to limit the gradients and distribute a certain amount of dissipation to the scheme corresponding to the ratio flux difference r in

Essentially Non-oscillatory Scheme
Consider the semi-discrete equation The numerical flux f that is k−th order accuracy approximate the physical flux f has to be satisfied the relation The approximated flux of a cell is reconstructed as a polynomial of a degree at most k − 1 via the primitive function F by using the appropriate chosen stencils (e.g. where r and s are the number of the chosen stencils on the right and the left of cell i ). For the uniform grid, the numerical flux at cell interface is given by [13] where the constants rj c and rj c for uniform grid were evaluated depends on the order of accuracy k and the number of 'smoother' stencil r as in [13]. The procedure for choosing appropriate stencil requires the divided differences. The large divided difference indicates large or discontinuous derivatives, therefore the 'smoother' stencil is the stencil having less divided difference in absolute value. The divided difference is defined as follow The first-order accurate ENO scheme is identical with the first-order upwind scheme The numerical flux of the second-order accurate ENO scheme is given by [13]     The ENO-Roe approach is used to ensure the upwinding by using the eigenvalues and eigenvectors of the Roe-averaged matrix to locally form the characteristic form [13]. Besides, the Lax-Friedrichs splitting is also considered for the sonic points [13]   over the relevant range of u.

NUMERICAL EXPERIMENTS
This section contains the results obtained by the developed code written in Fortran based on [3].

Planar Shock Tube Problem
The first numerical experiment was to validate the code by solving the classical shock tube problem, in which the driver section with high pressure is separated from the driven section with low pressure by a barrier. When the barrier is removed the shock wave and the contact discontinuity propagate to the lower pressure section and the expansion wave propagates to the higher pressure section. The initial conditions are given as in [2]  The initial conditions are modified from the original initial conditions studied by Sod [5] which is a very common test case for numerical methods in gas dynamics . The problem was solved in Cartesian coordinate with the number of cells is 100, makes the grid width to be ∆x = 0.01.
The CFL coefficient for all four numerical methods was chosen as λ = 0.9 to satisfy the CFL condition. The fluid is assumed to be a calorically perfect gas with ratio of specific heat γ = 1.4. The exact solution obtained from [2], which consists of left rarefaction, contact discontinuity and right shock wave, is also presented for comparison. The density and the pressure distribution at time t = 0.2sec are shown in Figure 1 and Figure 2. All four numerical schemes yielded similar results in the smooth regions. High order schemes, i.e. Sweby's and ENO scheme, which satisfy the TVD condition, exhibit oscillation-free in the vicinity of flow discontinuities. Whereas the first-order schemes expressed the dissipative nature over the flow discontinuities. Figure 3 shows the density distribution in the vicinity of left rarefaction wave, the solution of FVS method was smeared the most and even had unphysical expansion shock at the expansive sonic point. By applying Harten's entropy correction, the Roe scheme has overcome this error. Density distribution in the vicinity of contact discontinuity are presented in Figure 4 , the numerical solution of FVS scheme was dissipated the most, smeared out over 20 numerical cells, while the Third-order accuracy ENO scheme had the least dissipative solution, was spread over 14 numerical cells. The numerical results in the vicinity of shock wave is presented in Figure 5, the shock wave is more sharply resolved by the FVS and Sweby's scheme than Roe's and ENO scheme. Next, another Riemann problem with extremely strong discontinuity is investigated. This problem is chosen to examined due to is resemble with the blast wave problem that would be investigated in the next section. The initial conditions is taken as in [8]  This problem was also solved in Cartesian coordinate with the same grid width as well as other assumptions in the above problem ∆x = 0.01. In this problem there are two computing cells between the tail of the left rarefaction wave and the con-tact discontinuity and between the contact discontinuity and the shock wave there are three computing cells, therefore it was difficult for the numerical schemes to capture correctly the location of contact discontinuity and shock wave. In order to see the difficulties of numerical schemes when resolving extremely strong shock, the density and pressure distribution in the vicinity of the contact discontinuity and the shock wave is presented in Figure 6 and Figure 7. Although the predictions from all four numerical methods preserved monotonicity near the shock, the numerical solutions were strongly dissipated and overestimated the location of contact discontinuity and the shock wave. As expected, the ENO scheme was the least dissipative scheme, however it over-estimated shock wave location more than the other methods. The Sweby's scheme dissipated as much as the two first-order scheme but had the best pre-diction of shock wave location of the four schemes.

Blast Wave Simulation
In this section, the blast wave from a spherical explosion of 1kg charge TNT as in [8] was simulated by solving the Euler equations in the spherical co-ordinate system. This problem could be considered as a Riemann problem by assuming the charge had spherical shape with the radius TNT r calculated from the specified charge weight and the charge density, taken as 3 1600 / kg m .
Inside this radius was the high pressure region obtained by assuming the explosive material TNT transformed instantly into gas phase when detonated, the pressure was taken as 8.381GPa, obtained using the blast energy of TNT. The temperature effect was not taken into account, the ambient air outside the sphere was treated as a calorically perfect gas with density and pressure of omitted. The calculation domain was r = 10m with grid width ∆r = 5mm, which was sufficient to model the explosive charge by approximately 10 computing cells. Figure 8 presents the variation of over-pressure, the pressure above (or below) ambient atmospheric pressure, measured from the center of the explosive. The numerical results are compared with the incident pressure of spherical free-air burst obtained from the Conventional Weapons Effects pro-gram ConWep [18] and the numerical solutions of the Runge-Kutta Discontinuous Galerkin Method (RKDG) obtained from [8]. All numerical results agree well with numerical results in [8], however the overpressure was overestimated compared to data from ConWep. Figure 9 presents the pressure distribution at different times obtained by different numerical methods. It was well agreed with [8] that the formation of secondary shock wave following the primary shock wave then it propagated inward and being swept out by the rarefaction wave.

Two Dimensional Blast Wave
The performance of different numerical schemes in two dimensional using the dimensional splitting method is assessed for the cylindrical blast wave.
The two dimensional Euler equations were solved on a Cartesian coordinate system. The computational domain was a 2 × 2 square domain consists of two regions separated by a circle at center with radius r = 0.5. The initial flow variables were constant in each of the two regions and had a circular discontinuity at time t = 0.
These initial values of computing cells cutting the circular discontinuity were modified areaweighted in order to avoid the staircase configuration of the data. The initial conditions are depicted as in Figure Figure 11 presents the density distribution results of numerical methods at the time t = 0.25. Similar as in the one dimensional problem, the circular shock wave and the circular contact discontinuity propagate outwards from the center and the rarefaction wave travels inwards to the center. Figure 12 shows the corresponding pressure distribution. Since the pressure is continuous across the contact discontinuity, the solutions exhibit a circular shock and circular rarefaction wave. The density distribution and the pressure distribution along the radial line are also compared with the solution obtained from [2] in Figure 13 and Figure

CONCLUSION
The moving shock wave problem was simulated with finite difference method. The fluxes at cell interfaces were calculated by four different schemes. The linearization of the Jacobian matrix, i.e. utilizing the Roe-averaged matrix, plays an important role in the Roe's Firstorder Upwind scheme, the Sweby's TVD scheme and the ENO scheme. While the presence shocks in the two first-order upwind schemes and ENO scheme is determined by the sign of wave speed, in the Sweby's TVD scheme it is the ratio of solution differences that indicates the shocks.
The adapted code was first tested in the modified planar shock tube problem with weak discontinuities and then solved for another onedimensional Riemann problem with extremely strong discontinuities which resembles the explosion problem. The obtained solution showed the disadvantages of two first-order schemes. The Sweby's TVD scheme although had less computational cost than the ENO scheme but had better resolved shock solution for the extremely strong shock problem. The simulation of the explosion from 1kg TNT was well agreed with other author but still overpredicted compared to empirical results.