Open access

Blast Wave Simulations with a Runge-Kutta Discontinuous Galerkin Method

Written By

Emre Alpman

Submitted: 08 December 2011 Published: 07 November 2012

DOI: 10.5772/48613

From the Edited Volume

Advances in Modeling of Fluid Dynamics

Edited by Chaoqun Liu

Chapter metrics overview

3,259 Chapter Downloads

View Full Metrics

1. Introduction

An explosion generates high pressure and temperature gases which expand into the surrounding medium to generate a spherical shock wave called a blast wave [1]. Experimental methods to simulate blast waves require handling real explosives. Hence, they may be dangerous and very expensive. On the other hand, computational fluid dynamics (CFD) could be used to generalize and support experimental results in simulating blast waves. Sharma and Long [2] used Direct Simulation Monte Carlo (DSMC) method to simulate blast waves. DSMC is a very powerful method for detonations and flows with chemical reactions [3], [4], [5]. However, being originally designed for rarefied gas dynamics, it can be very expensive for continuum flows. Also, if one makes a local thermodynamic equilibrium assumption in DSMC, the method becomes equivalent to solving Euler equations [1]. In references [6] - [10] blast wave simulations were performed by solving Euler equations using a finite volume method where a second order accuracy was achieved by extrapolating density, velocity and pressure at the cell interfaces. The fluxes at cell faces were calculated using the AUSM+ method [11] due to its algorithmic simplicity and suitability for flows with discontinuities. Resulting semi-discrete equations were solved using a 2 stage modified Runge-Kutta method [12].

Discontinuous Galerkin (DG) methods, which have the features of both finite volume and finite element methods, have become popular for the solution of hyperbolic conservation laws like Euler equations [13]. DG method represents solution on elements by a collection of piecewise discontinuous functions hence sometimes considered as a high order accurate extension of finite volume method [13]. One class of DG methods is the Runge-Kutta DG method (RKDG) where spatial discretization is performed using polynomials which are discontinuous across element faces. Then the resulting system of ordinary differential equations is solved using a total variation diminishing (TVD) Runge-Kutta scheme [14]. The details of the RKDG scheme can be found in references [15]-[20]. An alternative to the RKDG method is the Space-Time Discontinuous Galerkin Method (STDG) in which the space and time discretizations are not separated [21], [22]. This method is particularly suitable for numerical solutions on adaptive meshes.

In this study RKDG method was selected for numerical solutions because of its relative programming simplicity and its similarity to the semi-discrete approach adopted for the finite volume solutions described above. Simulations were performed for one-dimensional moving rectangular and spherical shock waves (blast waves). Planar shock waves were simulated by solving planar shock tube problem with moderate and extremely strong discontinuities [23]. These two cases were compared in order to display the difficulties numerical methods may experience when the discontinuities in the flowfield are strong, as in the case of blast waves due to explosions. The predictions were also compared to the exact solutions by Sod [24]. Spherical shock waves were simulated by solving a spherical shock tube problem and blast waves generated from explosion of 1 kg of TNT. Here explosive material was assumed to transform gas phase instantly after the detonation, thus explosive was modeled as a high pressure sphere. This way, blast wave problem also became a spherical shock tube problem where an isobaric sphere of a high pressure gas expanded toward surrounding low pressure air.

Section 2 contains the theory of RKDG method and Section 3 describes the methodology followed for numerical solutions. Numerical solutions are presented and discussed in Section 4 and finally the conclusions are drawn in Section 5.

Advertisement

2. Theory

One-dimensional Euler equations are written as:

qt+1rn(rnf)r+h=0E1
with
q=[ρρuρE]TE2
f=[ρρu2+pρuH]TE3
h=nr[0p0]TE4

In the above equations ρ is fluid density, u is fluid velocity, p is static pressure, E and H are total energy and total enthalpy per unit mass.

E=pρ(γ1)+VV2E5
H=E+pρE6

In equation (4) n = 0 for Cartesian coordinates, n = 1 for cylindrical coordinates and n = 2 for spherical coordinates.

The DG method is derived by multiplying equation (1) by a test function (r) and integrating over a control volume (or element) (K)[20]

Kqtϕ(r)dK+K1rn(rnf)rϕ(r)dK+Khϕ(r)dK=0E7

where dK=rndr for one-dimensional flow. Using this relation fordK, the second integral in equation (7) can be rewritten as for one-dimensional flow as:

K1rn(rnf)rϕ(r)dK=ri1/2ri+1/2(rnf)rϕ(r)dr=ri1/2ri+1/2(rnfϕ(r))rdrri1/2ri+1/2rnfdϕ(r)drdrE8

Therefore, for one space dimension equation (7) takes the following form for an element i [20]:

ri1/2ri+1/2qtϕ(r)rndr+[frnϕ(r)]i+1/2[frnϕ(r)]i1/2ri1/2ri+1/2fdϕ(r)drrndr+ri1/2ri+1/2hϕ(r)rndr=0E9

In the DG method the numerical solution q is represented using a collection of polynomials on an element. Since polynomial continuity across element faces is not enforced the solution may be discontinuous at the element faces and which results in a multi-valued the flux vector f at these locations. This problem is overcome by using numerical flux f^ at the faces, which must be a conservative and upwind flux [19], [20]. Note that equation (9) reduces to finite volume formulation for (r) = 1 [19].

Normally in a DG method the solution q and the test function (r) belongs to the space of polynomials of degree smaller or equal to some k. This way a (k+1)st order accurate approximation to the solution can be obtained [19],[20] The approximate solution can then be written as:

q(r,t)=j=1kqj(t)ϕj(r)E10

where qj(t) are called degrees of freedom of q and j(r) are the bases of the polynomial solution space. Using the approximation given in equation (10) the first integral in equation (9) is written as:

ri1/2ri+1/2qtϕ(r)rndr=j=1kd(qj(t))dtri1/2ri+1/2ϕj(r)ϕl(r)rndrE11

Substituting the right hand side of equation (11) into equation (9) one can get a system of ordinary differential equations for the degrees of freedom of q. The integral on the right hand side of equation (11) is an inner product of the bases of the polynomial solution space with respect to the weight function r n and it yields a matrix called the mass matrix [20]. When orthogonal polynomials are chosen as bases then this matrix becomes diagonal and can easily be inverted [13]. For solution in Cartesian coordinates (n = 0), the weighting function becomes unity therefore Legendre polynomials are typically used as bases [19] - [22]. For this case (n = 0) the mass matrix, M, for an element i becomes

Mjl=ri1/2ri+1/2ϕj(r)ϕl(r)dr=Δri2j+1δjlE12

where δjl is the Kronecker’s delta and Δri=ri+1/2ri1/2 is the grid size. Here the bases are defined as [20]

ϕj(r)=Pj(2(rri)Δri)E13

where Pj(r) is the Legendre polynomial of degree j. Therefore, a set of ordinary differential equations for the degrees of freedom can be obtained for Cartesian coordinates as:

dqjdt=2j+1Δri[[f^ϕj(r)]i+1/2[f^ϕj(r)]i1/2ri1/2ri+1/2fdϕj(r)drdr+ri1/2ri+1/2hϕj(r)dr]E14

Using Legendre polynomials as bases not only yields a diagonal mass matrix, but also the diagonal elements have a simple form as displayed in equation (13) Unfortunately, this simplicity Legendre polynomials provide for Cartesian coordinates does not extend to cylindrical or spherical coordinates because Legendre polynomials are not orthogonal with respect to the weight function r n when n is not equal to zero. The mass matrix corresponding to the weight r n and k = 1 is given in equation (15)

M=[(rn+1n+1)|ri1/2ri+1/22Δri(rn+2n+2rirn+1n+1)|ri1/2ri+1/22Δri(rn+2n+2rirn+1n+1)|ri1/2ri+1/24Δri2(rn+3n+32rirn+2n+2+ri2rn+1n+1)|ri1/2ri+1/2]E15

For k = 2 the corresponding matrix will be 3x3 and the additional elements will be

M13=M31=6Δri2(rn+3n+32rirn+2n+2+(ri2Δri212)rn+1n+1)|ri1/2ri+1/2E16
M23=M32=12Δri3(rn+4n+43rirn+3n+3+(3ri2Δri212)rn+2n+2(12ri3riΔri212)rn+1n+1)|ri1/2ri+1/2E17
M33=36Δri4(rn+5n+54rirn+4n+4+(ri26Δri26)rn+3n+3(ri39riΔri23)rn+2n+2++(ri4ri2Δri26+Δri4144)rn+1n+1)|ri1/2ri+1/2E18

Even for k = 1 the mass matrix is complicated for both cylindrical and spherical cases. Furthermore, the elements of the matrix are functions of element locations. Therefore, the elements of the inverse matrix must be computed for each element.

As an alternative approach one could search for a set of orthogonal polynomials with respect to the weight r n. Reference [25] contains a method to derive such a set of polynomials starting from known orthogonal polynomials. Starting from Legendre polynomials and following the procedure described in [25], polynomials orthogonal with respect to the weight r n can be derived as

ϕ0(r)=1E19
ϕ1(r)=2Δri(rn+1n+2ri+1/2n+2ri1/2n+2ri+1/2n+1ri1/2n+1)E20
ϕ2(r)=C1+C22(rri)Δri+12(12(rri)2Δri21)E21

where [25]

C1=I0,1I1,2I0,2I1,1I0,0I1,1I0,1I1,0E22
C2=I0,2I1,0I0,0I1,2I0,0I1,1I0,1I1,0E23

and [25]

Ir,s=ri1/2ri+1/2Ps(2(rri)Δri)rn+rdrE24

In equation (24) Ps(r) is the Legendre polynomial of degree s. One can clearly see from equations (19), (20) and (21) that these set of polynomials, especially the second order one, are more complicated compared to Legendre polynomials. Their advantage due to their orthogonality is offset by the complexity they bring to the terms on right hand side of the equation (9) for spherical or cylindrical coordinates. Therefore, this approach does not produce a good alternative.

Another approach for a simple formulization of equation (9) can be considered by changing the definition of the approximate solution originally given in equation (10). The first integral in equation (9) can be rewritten as:

ri1/2ri+1/2qtϕ(r)rndr=ri1/2ri+1/2(rnq)tϕ(r)drE25

Furthermore, if the definition of approximate solution is changed to be:

rnq(r,t)=j=1kqj(t)ϕj(r)E26

Then the following set of ordinary differential equations for the degrees of freedom can be obtained

dqjdt=2j+1Δri[[f^ϕjrn]i+1/2[f^ϕjrn]i1/2ri1/2ri+1/2fdϕj(r)drrndr+ri1/2ri+1/2hϕj(r)rndr]E27

Equation (27) provides a simple and similar formulation for all coordinate systems. Therefore, it looks like a very good approach. However, during the solutions, this approach leaded to difficulties in applying limiters (see Section 3 about limiters) during the calculation of numerical fluxes at the element faces.

Considering the approaches described above an approximation in the formulations is proposed which will use equation (10) and Legendre polynomials to define the solution in the elements however, will yield a set of ordinary differential equations similar to equation (27). This approximation is done in the first integral of equation (9). Here, the term r n inside the integral is replaced by its average inside the element, which is nothing but the location of the center of that element ri. This integral is approximated as:

ri1/2ri+1/2qtϕ(r)rndrrinri1/2ri+1/2qtϕ(r)drE28

Then using equation (10) to describe the solution in an element, the following equations are obtained for the degrees of freedom:

dqjdt=2j+1rinΔri[[f^ϕjrn]i+1/2[f^ϕjrn]i1/2ri1/2ri+1/2fdϕj(r)drrndr+ri1/2ri+1/2hϕj(r)rndr]E29

The approximation in equation (29) not only provides a simple formulation for cylindrical and spherical coordinates, but also turns in to the exact formulation (equation (14)) for Cartesian coordinates.

Advertisement

3. Methodology

Numerical solutions were performed using an in-house computer code written in C language. The code solves equation (29) in Cartesian, cylindrical, and spherical coordinates and then constructs the numerical solution using equation (10). Numerical solution of equation (29) requires a numerical upwind flux function. The selected flux function affects the quality of the solution. Performance of RKDG method with different numerical fluxes can be found in [26]. Numerous different methods were used in literature for flow problems involving shock waves. Among these methods AUSM (Advection Upstream Splitting Method) -family of methods provide algorithmic simplicity by using a scalar diffusion term [27]. Hence, unlike many other upwind type methods, AUSM-family methods do not require the knowledge of the eigenstructure of the flow problem and this helps generate very efficient implementations. In this study the AUSM+ method [11] was used to calculate the numerical fluxes mainly because this method has previously been used by the author for moving shock wave simulations with satisfactory results [6] - [10]. Performance of this method was also tested against performance of van Leer [28], HLLE [29] and Roe’s [30] upwind methods. Numerical solutions were obtained for k =1 and 2 which correspond to second and third order spatial accuracies, respectively.

The integrals in equation (29) were obtained using a quadrature method. According to [19], for a given value of k, the applied quadrature rule must be exact for polynomials of degree 2k for the interior of the elements and three-point Gauss-Legendre rule [31] was used in that reference for k = 2. This method however, required calculation of flow variables at the interior locations of an element. The three-point Simpson’s rule [31], on the other hand, used the flow variables at element faces, which were already computed for flux calculations, and element center. Thus, it provided programming simplicity. Therefore, this rule was also used in the solutions although it was exact up to the cubic polynomials only. Numerical solutions obtained using these two integration techniques were compared with each other.

A numerical method used for a flow with discontinuities should yield entropy satisfying weak solutions [32]. The so-called monotone methods yield such solutions, but they are only first order accurate, they perform poorly near discontinuities due to high amounts of artificial viscosity present in them [32]. High order accurate methods can be used to decrease artificial viscosity; however they yield dispersion errors which cause spurious oscillations in the vicinity of flow discontinuities. Therefore, high resolution methods, which provide high accuracy in smooth regions of the flow and apply some kind of stabilization near discontinuities, are required to get an accurate description of the flow [32]. Typically, DG method is monotone for k = 0 only and solutions for k > 0 generate spurious oscillations in the vicinity of flow discontinuities [33]. Therefore, some kind of stabilization strategy is needed. Cockburn et.al. [13], [15], [20] adopted a generalized slope limiter with a TVB corrected minmod limiter [16], Calle et. al. used an artificial diffusion term in their stabilized discontinuous Galerkin method [34], van der Vegt and van der Ven [21], [22] used a stabilization operator, and Qui and Shu [35] used WENO type limiters to avoid spurious oscillations in their studies. Among these limiting approaches, slope limiting approach, which is also used in finite volume schemes, is computationally the simplest. Therefore, a slope limiting procedure was used in this study. In the slope limiting approach, gradients in each element (or cell) were limited by using a limiter function. There are many limiter functions available in the literature [36]. Details of limiting procedure and some comparisons can be found in [37].

It is known that minmod limiter is the most dissipative limiter [30], [37] and the TVB correction mentioned in [16], [19], [20] is problem dependent. On the other hand, Superbee limiter [30] is the least dissipative and is suitable for flows with very strong shocks like a blast wave [23]. In this study minmod and Superbee limiters along with the Sweby limiter [37], which is between Superbee and minmod limiters, were used for numerical solutions. During the limiting procedure, the ratios the successive gradients of (j-1)st degree of freedom of an element to the j th degree of freedom of the same element were calculated, and minimum of these ratios were supplied as input to limiting function. Then the j th degree of freedom was multiplied by the output of the limiter. This approach is similar to the one followed in [38].

After the above steps were followed, the resulting system of ordinary differential equations for the degrees of freedom (equation(29)) could be solved using a Runge-Kutta technique. It was shown in Ref. [14] if the Runge-Kutta method used does not satisfy the TVD property, spurious oscillations may appear in the solution of nonlinear problems, even when a TVD space discretization is used. Therefore, equation (29) was solved using a 2nd order accurate TVD Runge-Kutta scheme described in [14].

Advertisement

4. Results and discussion

This section includes numerical simulations of planar and spherical shock waves using the RKDG method described above. Since the planar shock tube problem had exact solution; effects of limiters, flux functions and quadrature rule on solutions were analyzed and presented for that case.

4.1. Planar shock tube problem

Since blast wave simulation is basically a moving shock wave problem, numerical solution procedure was first tested for Sod’s shock tube problem [24]. In this problem a stationary high pressure fluid was separated from a stationary low pressure fluid by a barrier. At t = 0 the barrier was removed. This leaded to a shock wave and a contact discontinuity move towards the low pressure region and an expansion fan move towards the high pressure region. In this study two shock tube problems with different initial conditions were solved. Predictions were compared with the exact solutions [24].

4.1.1. Shock tube with moderate discontinuities

In this case the initial conditions were given as follows:

ρ = 1 kg/m3, u = 0, p = 1x105 Pa, r ≤ m

ρ = 0.125 kg/m3, u = 0, p = 1x104 Pa, r > m

These initial conditions were the original initial conditions studied by Sod [24] and this problem is a very popular test case for numerical methods for gas dynamics. The problem was solved in Cartesian coordinates and size of the solution domain was 5m. Fluid in the problem was calorically perfect air. Figure 1 shows density distribution at t = 0.002 seconds obtained with 1001 grid points. This corresponded to a mesh spacing of 0.5 cm. Predictions included RKDG method, finite volume (FV) [6], [7] solutions with minmod and superbee limiters, and first order accurate solutions. All the second order accurate methods yielded similar results except in the vicinity of flow discontinuities and the head expansion wave. On the other hand, dissipative nature of first order accurate solution could be easily seen in the vicinity of the contact surface. In order to better see the differences between the predictions of different methods, density distribution in the vicinity of the head expansion wave, contact surface and shock wave were displayed in Figure 2, Figure 3 and Figure 4. A close examination of these figures clearly indicates that RKDG method is less dissipative than the finite volume method using the same numerical flux function. Among the solutions RKDG method with superbee limiter is the least dissipative one. Spurious oscillations are evident for this case in the vicinity of the head expansion wave (Figure 2). But due its least dissipative nature, this method was able to resolve the contact discontinuity (Figure 3) and the shock wave (Figure 4) much better than the other methods did.

Figure 1.

Density distribution at t = 0.002s. (ρ0 = 0.125 kg/m3)

Figure 2.

Density distribution at t = 0.002s in the vicinity of the head expansion wave (ρ0 = 0.125 kg/m3)

Figure 3.

Density distribution at t = 0.002s in the vicinity of the contact surface (ρ0 = 0.125 kg/m3)

Figure 4.

Density distribution at t = 0.002s in the vicinity of the shock wave (ρ0 = 0.125 kg/m3)

4.1.2. Shock tube with extremely strong discontinuities

In this case the initial conditions were given as follows:

ρ = 1 kg/m3, u = 0, p = 1x105 Pa, r ≤ 2m

ρ = 0.0001 kg/m3, u = 0, p = 10 Pa, r > 2m

These were the same conditions studied in [23] as an extremely strong discontinuity case. This case was selected because it more closely resembled a blast wave simulation problem and it constituted an extremely difficult case for numerical methods because of large gradients in the flowfield. The problem was solved in Cartesian coordinates, size of the solution domain was 7m and fluid was calorically perfect air. Numerical solutions were obtained using a mesh with N = 701 grid points. This corresponds to a mesh spacing of 1 cm. Figure 5 shows density distribution at t = 1.2 ms in the vicinity of the contact surface and shock wave. Predictions were obtained using FV and RKDG method with k = 1 using Superbee (sb), minmod (mm) and Sweby (sw) limiters. This way both methods would have the same spatial accuracy. Exact solution was also displayed for comparison. This figure clearly shows the difficulty the methods had experienced while resolving the strong discontinuities. Note the logarithmic scale used for the vertical axis. Predictions obtained using minmod limiter was nearly the same for FV and RKDG methods, in which, results showed a much smeared contact surface and over-predicted the shock location. FV method with Superbee limiter performed relatively well although the numerical solutions are on the overall poor. In order to improve the results solutions were performed on a finer mesh which was obtained by halving the mesh spacing using N = 1401 grid points. Density distributions for this case were displayed in Figure 6. Greatest improvement was observed for FV method with the Superbee limiter. Other solutions were also improved slightly but not as much. Overall, the RKDG method predictions were poor compared to FV method predictions obtained using the same limiter despite the fact that RKDG method with Superbee limiter had yielded the best predictions for the case with moderate discontinuities. It was also interesting to note that RKDG solutions took nearly 1.5 times more CPU time than FV solutions although this was not a big issue for one-dimensional flow simulations.

Figure 5.

Density distribution at t = 1.2 ms in the vicinity of the contact surface and shock wave. (ρ0 = 0.0001 kg/m3, N = 701)

Figure 6.

Density distribution at t = 1.2 ms in the vicinity of the contact surface and shock wave. (ρ0 = 0.0001 kg/m3, N = 1401)

Pressure distributions at t = 1.2 ms obtained using N = 1401 grid points were shown in Figure 7. In this figure all numerical solutions experienced a spurious oscillation right after the expansion region. This location and the location of the shock wave were the places where discrepancies between numerical solutions were most severe. They were pretty much in agreement everywhere else in the solution domain. However, in Figure 6 there were also considerable differences between numerical solutions in the vicinity of the contact surface location which was not observed in Figure 7. Therefore, it was concluded that prediction of shock location could be improved if the contact surface could be resolved better.

Figure 7.

Pressure distribution at t = 1.2 ms. (p0 = 10 N/m2)

Comparing Figure 6 and Figure 7 one can see that spurious oscillations were observed in the vicinity of the contact surface rather than the shock wave. These oscillations could be suppressed by using a more dissipative limiter however, this was shown to negatively affect the overall solution quality; solutions obtained using minmod limiter were the worst compared to others. Knowing that only density experiences a discontinuity across the contact surface while pressure and velocity remain continuous, an alternative limiting strategy was tested to see if it could improve the predictions. In this study limiting was performed on the primitive flow variables; density, velocity and pressure. Therefore, in this alternative strategy minmod limiter was applied for density only, while less dissipative Superbee and Sweby limiters were applied for pressure and velocity. Density distributions obtained using this limiting strategy was displayed in Figure 8. In the figure, legend first limiter was for density and second limiter was for pressure and velocity. The alternative limiting strategy, especially using minmod and Sweby limiters, improved RKDG predictions greatly. This strategy also suppressed oscillations observed in the previous FV solutions however; it also negatively affected the prediction of shock location. So far, the best predictions were obtained with the Superbee limiter for the FV method and with the minmod-Sweby combination for the RKDG method. These predictions were compared with exact solution in Figure 9. It is clear that RKDG method with this alternative limiting successfully suppressed oscillations after the expansion region and in the vicinity of the shock wave. Also it made a slightly better shock prediction compared to FV method. Considering this improvement, RKDG solutions will employ this alternative limiting with minmod and Sweby limiters from this point forward.

Figure 8.

Density distribution at t = 1.2 ms in the vicinity of the contact surface and shock wave. (ρ0 = 0.0001 kg/m3, N = 1401)

In order to see the effect of mesh spacing, solutions displayed in Figure 9 were repeated on the previously used coarse mesh and results were displayed in Figure 10. It is clear that RKDG solution was less sensitive to mesh spacing compared to FV method.

Spatial accuracy of RKDG methods can be easily increased by using a high order polynomial in the element. This avoids using a larger stencil as in the case of FV or finite difference method [16]. For the problem considered density distribution obtained using k = 1 and k = 2 were displayed in Figure 11. Results obtained with these two polynomials are nearly the same except k = 2 case placed the shock slightly ahead of the k = 1 case. The reason using high order polynomial did not bring much benefit was that the limiters applied degraded the high order accuracy in the vicinity of the discontinuities. This might have been overcome by using a limiter which preserves the accuracy like the TVB corrected minmod limiter of reference [19] however that correction was problem dependent.

Figure 9.

Density distribution at t = 1.2 ms in the vicinity of the contact surface and shock wave. (ρ0 = 0.0001 kg/m3, N = 1401)

Figure 10.

Density distribution at t = 1.2 ms in the vicinity of the contact surface and shock wave. (ρ0 = 0.0001 kg/m3, N = 701)

Figure 11.

Density distribution at t = 1.2 ms obtained using k= 1 and 2. (ρ0 = 0.0001 kg/m3, N = 701)

4.1.3. Effect of the quadrature rule

Equation (29) contained integrals which were evaluated using a quadrature rule. The predictions displayed previously were obtained using three-point Simpson’s rule of integration which is third-order accurate [31]. In order to see the effect of the quadrature rule on predictions, RKDG solution with k = 2 was also performed using three-point Gauss-Legendre rule (5th order accurate [31]) and compared to the one with Simpson rule in Figure 12. According to this figure, numerical solution was nearly insensitive to the quadrature rule employed. Since three-point Simpson rule used flow variables at the faces (already calculated for flux calculations) and the center of the element it did not require calculation of flow variables at the interior locations of as did three-point Gauss-Legendre rule. Therefore, it was mainly preferred and employed in the solutions.

Figure 12.

Density distribution at t = 1.2 ms obtained using different quadrature rules.(ρ0 = 0.0001 kg/m3, N = 701)

4.1.4. Effect of numerical flux function

It is also known that the numerical flux function used in the computations affects the accuracy of the DG solutions. In order to see this effect, performance of the AUSM+ method [11] was compared to those of van Leer (VL) [28], HLLE [29] and Roe’s [30] methods. Density predictions obtained using these flux functions and their comparisons with the exact solution were displayed in Figure 13. According to this figure, performances of these methods were close to each other except the HLLE method which clearly underpredicted the shock location. Among the flux functions tested, AUSM+ showed the best performance as expected.

4.2. Blast wave simulations

In this section blast waves generated by an explosion were analyzed. Here explosive material was assumed to transform gas phase instantly after the detonation, thus explosive was modeled as a high pressure sphere. This way, blast wave problem became a spherical shock tube problem where an isobaric sphere of a high pressure gas expanded toward surrounding low pressure air. Numerical solutions were obtained using RKDG method with k = 1, AUSM+ flux function and the alternative limiting technique which gave the best results for the planar shock tube problem analyzed in the previous section. A grid was used with a mesh spacing of 4 mm in a 10 m solution domain.

Figure 13.

Density distribution at t = 1.2 ms obtained using different numerical flux functions.(ρ0 = 0.0001 kg/m3, N = 701)

Numerical solutions were first compared with the approximate analytic solutions of blast waves by Friedman [39]. The problem considered contained a sphere of air compressed to 22 times the ambient air pressure. This sphere had a non-dimensional radius of one [40]. Temperature was taken to be uniform in the entire solution domain and air was assumed to be calorically perfect. Figure 14 contained loci of resulting primary shock, contact surface and secondary shock. Predictions were also compared with numerical solutions by Brode [41]. In this figure horizontal axis represented distance from the center normalized by the initial sphere radius, r0 and vertical axis represented time normalized by r0 and ambient speed of sound, a0.

Figure 14.

Loci of primary shock, contact surface and secondary shock.

According to this figure RKDG solutions showed only qualitative agreement with Friedman’s approximate analytic result; however agreement with Brode’s results were much better especially during the implosion of the secondary shock and its reflection from the center. Outward motion of the secondary shock after its reflection and its interaction with the contact surface were also observed in Figure 14.

Next, blast waves generated by explosion of 1 kg of TNT were analyzed and simulated. The explosive was modeled as an isobaric high pressure sphere in which the density was 1600 kg/m3 and pressure was 8.8447 GPa which was obtained using the blast energy of TNT [42] and the JWL (Jones-Wilkins-Lee) equation of state [43]. Outside this sphere ambient air density and pressure were 1.225 kg/m3 and 101320 Pa, respectively. Simulations involved two different fluids; detonation products and ambient air for which different state equations were used. For detonation products JWL equation of state was used mainly due to its popularity [44]. For surrounding air two different cases were followed. First case assumed ambient air to be calorically perfect. However, after an explosion like this one, ambient air temperature may easily become very high so that calorically perfect gas assumption might cease to be valid due to gas dissociations [45]. Therefore, as the second case, high temperature effects were included by assuming local chemical equilibrium meaning that the chemical reactions occur instantaneously [45]. Equilibrium relations given in reference [46] were used for air to calculate ratio of specific heat capacities and temperature in terms of pressure and density.

One of the important variables in blast wave simulation is the over-pressure, which is the rise of pressure above ambient pressure downstream of the primary shock wave. Figure 15 showed comparison of over-pressure predictions obtained for calorically perfect and equilibrium air cases with the data obtained from [42] which states that these results are curve fits to the data used in the weapons effect calculation program CONWEP [47]. Both cases over-predicted over-pressure with calorically perfect case being slightly better. The main discrepancy between both cases occurred between r = 1 and 3 m. Except this region both curves were nearly parallel.

Another important variable in blast wave simulation is the speed of the primary shock, which can be used to calculate shock arrival times. Primary shock speeds computed using RKDG method were compared with data from reference [42] in Figure 16. Here numerical solutions gave higher shock speeds close to the explosive but agreement with data from reference [42] became much better after r = 2 m, especially for calorically perfect air solution.

Overall solutions performed by treating air as a calorically perfect gas were in better agreement with data from reference [42]. But it is known that high temperatures encountered during detonation of high explosives like TNT makes calorically perfect assumption invalid. Since this assumption ignores any chemical reaction that might occur, it may yield unrealistically high temperatures after a shock wave. Such problems are usually encountered for hypersonic flows over reentry vehicles (see reference [48], chapter 16). In order to check this for our problem, temperature behind the primary shock wave predicted by the two cases of air were displayed in Figure 17. It was clear from this figure that temperatures predicted using the two assumptions were very close except in the region between r = 1 and 3 m, where major discrepancies between these numerical solutions took place. Since temperature predicted using calorically perfect and equilibrium air assumptions were very close, the potential negative consequences of using the former, was not experienced in this problem. Nevertheless, the predicted temperatures at lower r values were still high enough to make non-constant. In order to see the effect of temperature on , its values calculated using the equilibrium air assumption was compared with the constant value of 1.4 for calorically perfect gas assumption. This comparison was displayed in Figure 18 where a considerable drop in can be easily observed when r is less than 3m. This explains the relatively less pressure and temperature drop yielded by the equilibrium air assumption for the same amount of outward expansion. (See Figure 15 and Figure 17).

Figure 15.

Variation of over-pressure with distance measured from the center.

Figure 16.

Variation of primary shock speed with distance from the center.

RKDG predictions for blast waves represented above were obtained using the alternative limiting approach because it yielded better results for planar shock wave problem where the strength of the discontinuities remain constant. However, in a spherical shock problem the strength of the discontinuities decrease as high pressure gases planar shock tube problem with moderate discontinuity strengths showed that the performances of different limiting strategies became closer to each other as the discontinuity strengths decreased. In order to see the effect of limiting procedure on blast wave simulations density distributions at t=0.1 and 1.6 ms were plotted in Figure 19 and Figure 20 for the blast wave generated by the explosion of 1 kg of TNT. Density was plotted in order to see the contact surface. Here the ambient air was assumed to be calorically perfect air. Predictions include RKDG solutions with minmod (mm), Sweby (sw), superbee (sb) and alternative (mm for density, sw for pressure and velocity) limiters. Finite volume solution with superbee limiter was also included in the figures. According to Figure 19 RKDG solutions were very close to each other especially for secondary shock wave. The major differences were in the vicinity of the contact surface where solutions with mm and mm-sw limiters yielded more smeared contact surfaces as expected. At the same time finite volume (FV) solution showed a considerable discrepancy with RKDG solutions. Nevertheless, the strength of the discontinuity at the contact surface was small compared to that of the planar shock tube problem studied above. Therefore, the discrepancies at the contact surface location did not affect the shock predictions considerably. This was also supported in Figure 20 which showed density distributions well after the secondary shock wave reflected from the origin and crossed the contact surface. Again RKDG predictions are very close to each other and disagreement with the FV solution is evident.

Figure 17.

Variation of temperature behind primary shock with distance from the center

Figure 18.

Variation of ratio of specific heat capacities behind primary shock with distance from the center.

Figure 19.

Density Distribution at t=0.1 ms

Figure 20.

Density Distribution at t=1.6 ms

The blast wave solutions presented above were obtained by solving equation (29) which was derived by making the approximation described in equation (28). In order to see the effect of this approximation on the predictions, solutions of this approximate formulation were compared to the solutions obtained using exact formulation, i.e., by inverting the mass matrix given in equation (15). The loci of primary shock, contact surface and secondary shock wave predicted using these two approaches were displayed in Figure 21. Here solid lines corresponded to solutions of equation (29) and dashed lines corresponded to solutions of exact formulation. Examination of Figure 21 revealed that the approximation made in equation (29) had very little effect on the primary shock wave. Here shock speeds were nearly identical and this meant that the shock strengths were also predicted nearly equal. There were, however, visible discrepancies in the secondary shock predictions. These were evident especially after this shock reflected from the origin. Nevertheless, after the predicted secondary shocks interacted with the contact surface, their trajectories became nearly parallel. This also meant very nearly identical shock speeds and strengths. The most serious difference between predictions was observed in contact surface predictions. However, this had very little or no effect on the primary shock strength and speed which are the two most important parameters in blast wave simulations.

Figure 21.

Loci of primary shock, contact surface and secondary shock wave. (Solid lines; approximate formulation (equation (29), dashed lines; exact formulation)

Advertisement

5. Conclusions

An implementation of RKDG method was performed to simulate moving planar and spherical shock waves (blast waves). The Legendre polynomials which were mainly used as basis polynomials for RKDG solutions in Cartesian coordinates did not yield a diagonal mass matrix when the formulation was extended to cylindrical and spherical coordinates. Attempts on inverting this non-diagonal mass matrix or finding a suitable set of orthogonal polynomials for cylindrical and spherical coordinates failed to yield a simplistic formulation which not only had the simplicity of the formulation in Cartesian coordinates, but also could be easily generalized for all coordinate systems. Therefore, an approximation to the integral containing the unsteady term was introduced which provided a formulation which possessed these properties. The resulting formulation was approximate for cylindrical and spherical coordinates; nevertheless, it was exact for Cartesian coordinates. In order to solve the resulting system of equations, an in-house computational fluid dynamics code which was able to solve Euler equations using a finite volume technique was further modified to include an RKDG method. The code was first tested with a Sod’s shock tube problem with the so-called moderate extremely strong discontinuities [23]. Predictions were compared with the predictions obtained using a finite volume technique and to the exact solutions [24]. It was observed that the limiter used in the solutions clearly affected the overall quality of the predictions. The RKDG method was not able to produce satisfactory results with neither of the limiter functions used and was out-performed by FV method using the same limiter. After observing that spurious oscillations occur around the contact surface rather than the shock wave and alternative limiting strategy which used more dissipative minmod limiter [37] for density and less dissipative Superbee or Sweby limiter [37] for pressure and velocity was tested. This strategy clearly improved the solutions of RKDG method and RKDG with minmod-Sweby combination gave the best results. It was interesting to see that this alternative limiting strategy did not produce a similar improvement for the FV method. Overall RKDG method with this so-called alternative limiting outperformed FV method for the analyzed shock tube problem with very strong discontinuities and it was shown to be less sensitive to grid coarsening compared to FV method.

The formulation derived contained integrals which were evaluated using a quadrature method. Here the performance of three-point Simpson’s rule was compared to that of three-point Gauss-Legendre rule and solutions were found to be in very close agreement with each other. Therefore, the simpler Simpson’s rule was mainly employed. The resulting method was tested for different flux functions and comparison of the results showed that the AUSM+ [11] method performed relatively better compared to the HLLE [29], van Leer [28] and Roe’s [30] methods.

Blast waves were simulated by modeling the explosion problem as a spherical shock tube problem. For blast waved generated by explosion of 1kg of TNT, the numerical solutions were performed with and without including the high temperature effects for the surrounding air. Although high temperatures were encountered in the solutions, using calorically perfect gas assumption for air did not produce negative effects, and even gave better agreement with blast wave data taken from Ref. [42]. This may be due to the fact that equilibrium relations given in [46] were obtained for hypersonic flow at the upper levels of the atmosphere where pressure and density is low.

Finally, the effect of the approximation made in the formulation on the predictions was analyzed by comparing the blast wave solutions obtained using exact and approximate formulations. There was very little influence on the primary shock wave, although some discrepancies were observed for the secondary shock and the contact surface. Nevertheless, the approximate formulation provided a simple formulation for all coordinates and yielded satisfactory results.

References

  1. 1. DeweyJ. M.2001Expanding Spherical Shocks (Blast Waves). Handbook of Shock Waves. 2441481
  2. 2. SharmaA.LongL. N.2004Numerical Simulation of the Blast Impact Problem using the Direct Simulation Monte Carlo (DSMC) Method. Journal of Computational Physics. 200211237
  3. 3. LongL. N.AndersonJ. B.2000The simulation of Detonations using a Monte Carlo Method. Rarefied Gas Dynamics Conference, Sydney, Australia.
  4. 4. AndersonJ. B.LongL. N.2003Direct Monte Carlo Simulation of Chemical Reaction Systems: Prediction of Ultrafast Detonations. Journal of Chemical Physics. 11831023110
  5. 5. AndersonJ. B.LongL. N.2002Direct Simulation of Pathological Detonations. 18th International Symposium on Rarefied Gas Dynamics, Vancouver, Canada.
  6. 6. AlpmanE.2009Blast Wave Simulations using Euler Equations and Adaptive Grids. 5th Ankara International Aerospace Conference, AIAC-2009093Ankara, Turkey.
  7. 7. AlpmanE.2009Simulation of Blast Waves using Euler Equations. 17th National Thermal Science and Technology Conference, Sivas, Turkey. (in Turkish )
  8. 8. ChenC. C.AlpmanE.LinzellD. G.LongL. N.2008Effectiveness of Advanced Coating Systems for Mitigating Blast Effects on Steel Components. Structures Under Shock and Impact X, Book Series: WIT Transactions on the Build Environment. 988594
  9. 9. ChenC. C.LinzellD. G.LongL. N.AlpmanE.2007Computational studies of polyurea coated steel plate under blast loads. 9th US National Congress on Computational Mechanics, San Francisco, CA.
  10. 10. AlpmanE.LongL. N.ChenC.C.LinzellD. G.2007Prediction of Blast Loads on a Deformable Steel Plate Using Euler Equations. 18th AIAA Computational Fluid Dynamics Conference, Miami, FL.
  11. 11. LiouM.S.1996A sequel to AUSM: AUSM+. Journal of Computational Physics. 129364382
  12. 12. HoffmannK. A.ChiangS. T.2004Computational Fluid Dynamics for Engineers, Volume I. Engineering Education SystemTM.
  13. 13. CockburnB.2001Devising Discontinuous Galerkin Methods for Non-Linear Hyperbolic Conservation Laws. Journal of Computational and Applied Mathematics. 128187204
  14. 14. GottliebS.ShuC. W.1998Total Variation Diminishing Runge-Kutta Schemes. Mathematics of Computation. 677385
  15. 15. CockburnB.ShuC. W.1991The Runge-Kutta Local Projection 1Galerkin Method for Scalar Conservation Laws. M2 AN. 23:337.
  16. 16. CockburnB.ShuC. W.1989TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Scalar Conservation Laws II: General Framework. Mathematics of Computation. 52:411.
  17. 17. CockburnB.LinS. Y.ShuC. W.1989TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Scalar Conservation Laws III: One-dimensional Systems. Journal of Computational Physics, 84:90.
  18. 18. CockburnB.HouS.ShuC. W.1990TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Scalar Conservation Laws IV; The Multidimensional Case. Mathematics of Computation. 54:545.
  19. 19. CockburnB.ShuC. W.1998The Runge-Kutta Discontinuous Galerkin Methods for Conservation Laws V: Multidimensional Systems. Journal of Computational Physics. 141199224
  20. 20. CockburnB.ShuC. W.2001The Runge-Kutta Discontinuous Galerkin Methods for Convection-Dominated Problems. Review Article, Journal of Scientific Computing. 16173261
  21. 21. van der VegtJ. J. W.van der VenH.2002Space-Time Discontinuous Galerkin Finite Element Method with Dynamic Grid Motion for Inviscid Compressible Flows: I. General Formulation. Journal of Computational Physics. 182546585
  22. 22. van der VegtJ. J. W.van der VenH.2002Space-Time Discontinuos Galerkin Finite Element Method with Dynamic Grid Motion for Inviscid Compressible Flows: II. Efficient Flux Quadrature. Computer Methods in Applied Mechanics and Engineering. 19147474780
  23. 23. TaiC. H.ChiangD. C.SuY. P.1997Explicit Time Marching Methods for The Time-Dependent Euler Computations. Journal of Computational Physics. 130191202
  24. 24. SodG. A.1978A survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws. Journal of Computational Physics. 27131
  25. 25. PriceT. E.1979Orthogonal Polynomials for Nonclassical Weight Functions. SIAM Journal on Numerical Analysis. 169991006
  26. 26. QiuJ.KhooB. C.ShuC. W.2006A Numerical Study for the Performance of the Runge-Kutta Discontinuous Galerkin Method based on Different Numerical Fluxes. Journal of Computational Physics. 212540565
  27. 27. LiouM.S.2000Mass Flux Schemes and Connection to Shock Instability. Journal of Computational Physics. 160623648
  28. 28. van LeerB.1982Flux-Vector Splitting for the Euler Equations. Lecture Notes in Physics, Springer- Verlag, 170507512
  29. 29. EinfeldtB.1988On Godunov Type Methods for Gas-Dynamics. SIAM Journal of Numerical Analysis. 25, 294 EOF
  30. 30. RoeP. L.1986Characteristic-based Schemes for Euler Equations. Annual Review of Fluid Mechanics. 18337365
  31. 31. ChapraS. C.CanaleR. P.2006Numerical Methods for Engineers. 5th Edition, Mc Graw-Hill, 626
  32. 32. Le VequeR. J.1992Numerical Methods for Conservation Laws. 2nd Edition, BirkHauser.
  33. 33. FlahertyJ. E.KrivodonovaL.RemacleJ. F.ShephardM. S.2002Aspects of Discontinuous Galerkin Methods for Hyperbolic Conservation Laws. Finite Elements in Analysis and Design. 38889908
  34. 34. CalleJ. L. D.DevlooP. R. B.GomesS. M.2005Stabilized Discontinuous Galerkin Methods for Hyperbolic Equations. Computer Methods in Applied Mechanics and Engineering. 19418611874
  35. 35. QiuJ.ShuC. W.2005Runge-Kutta Discontinuous Galerkin Method using WENO limiters. SIAM Journal of Scientific Computing. 26907929
  36. 36. http://en.wikipedia.org/wiki/Flux_limiter
  37. 37. SwebyP. K.1984High Resolution Schemes using Flux Limiters for Hyperbolic Conservation Laws. SIAM Journal on Numerical Analysis. 219951011
  38. 38. BiswasR.DevineK. D.FlahertyJ. E.1994Parallel, Adaptive, Finite Element Methods for Conservation Laws. Appl. Numer. Math. 14: 255, 1994.
  39. 39. FriedmanM. P.1 EOF15 EOF1961A simplified Analysis of Spherical and Cylindrical Blast Waves. Journal of Fluid Mechanics. 11: 1.
  40. 40. SachdevP. L.2004Shock Waves and Explosions. Monographs and Surveys in Pure and Applied Mathematics. Chapman & Hall/CRC, 208224
  41. 41. BrodeH. L.1957Theoretical Solutions of Spherical Shock Tube Blasts. Rand Corporation Report. RM- 1974.
  42. 42. SmithP. D.HetheringtonJ. D.1994Blast and Ballistic Loading of Structures, Butterworth-Heinemann. 2462
  43. 43. DobratzB. M.CrawfordP. C.1985LLNL Explosives Handbook- Properties of Chemical Explosives and Explosive Simulants. Lawrence Livermore National Laboratory.
  44. 44. KubotaS.NagayamaK.SaburiT.OgataY.2007State Relations for a Two-Phase Mixture of Reacting explosives and Applications. Combustion and Flame. 7484
  45. 45. HoffmanK.A.ChiangS. T.2000Computational Fluid Dynamics, Volume II, Engineering Education SystemTM.
  46. 46. TannehillJ. C.MuggeP. H.1974Improved Curve Fits for the Thermodynamic Properties Equilibrium Air Suitable for Numerical Computation using Time Dependent and Shock-Capturing Methods. NASA CR-2470, 1974.
  47. 47. HydeD. W.1991CONWEP: Conventional Weapons Effects Program. US Army Waterways Experimental Station.
  48. 48. AndersonJ. D.2004Modern Compressible Flow with Historical Perspective, 3rd Edition, McGraw Hill, 587

Written By

Emre Alpman

Submitted: 08 December 2011 Published: 07 November 2012