JAMP  Vol.5 No.2 , February 2017
Postshock Oscillations on Non-Uniform Mesh
Abstract: An investigation of postshock oscillations on non-uniform grids is performed in this paper. These oscillations are generated as shock passes through the grid interfaces. The LLF scheme is checked for 1D and 2D problems on the discontinuous grids. Oscillations are observed only for nonlinear systems and the solutions of the scalar conservation laws and linear systems behave logically. The integral curves suggest underlying properties of these oscillations. The results of the paper reveal a flaw that adaptive methods for conservation laws have to refine grids at each time step.

1. Introduction

This paper considers the numerical solutions of hyperbolic systems of conservation laws on discontinuous grids. The discontinuity in the grid is often due to the overlapping of different mesh systems, or to the adaptive mesh refinement. In either case, if the discontinuities in solutions are genuine nonlinear, some oscillations will be conserved around the grid interface for nonlinear systems.

Computations on discontinue grids are becoming more common for two reasons. First, many multidimensional problems of practice interest involve complex geometry, and in general it is not sufficient to be able solve hyperbolic equations on a uniform Cartesian grid in a rectangular domain. As the configurations that can be modeled become more complex, so does the grid generation problem. Generally, it is very difficult to generate one smooth body-fitted grid to cover the whole complex domain. To simplify this procedure, it is accustomed to use multigrid to fit the complex domain. Each part of the domain will still have a smooth grid, but now the component grids will in general overlap rather than in an irregular fashion.

The second reason of these discontinue grids comes from the use of adaptive methods [1] [2] [3] [4] , where a grid is abruptly refined in order to gain high- resolution solution around areas of large solution variation. In this approach, many fine grids are distributed on those regions where high resolution is required. Efficient implementation of the adaptive strategy can increase the accuracy of the numerical approximations and meanwhile decrease the computational cost.

The main objective of this paper is to investigate postshock oscillations around the grid interface. Some linear and nonlinear equations are tested with LLF scheme, and this kind of oscillations only appears in nonlinear equation system. In addition, Godunov scheme behaves similarly to the LLF scheme, and so we don’t give the numerical results of Godunov scheme for saving space. We incline to believe that the postshock oscillations are an inevitable feature of shocks captured by currently frequently-used methods. The postshock oscillations due to discontinue grids are different from common nonphysical oscillations caused by high-order interpolation, but it is very similar to that appear behind the slow shocks [5] [6] . The amplitude and concave-convex shape of oscillations are determined by the ratio between coarse- and fine-grid length. From the point of Riemann invariant, these oscillations have similar feature with rarefaction. That is Riemann invariant has the same value at all points on the corresponding oscillation.

The rest of the paper is organized as follows. In Section 2, we briefly review some model problems and give the mesh distribution. In Section 3, the LLF scheme is described. The numerical results and the observed behavior are demonstrated in Section 4. However, we have yet found the reason which leads to the postshock oscillations observed in the paper.

2. Equations and the Grid Distribution

The postshock oscillations are investigated under the following five sets of equations: the advection equation,

u t + u x = 0 ; (1)

the inviscid Burgers equation,

u t + ( u 2 2 ) x = 0 ; (2)

the shallow water equation,

[ h h u ] t + [ h u h u 2 + 1 2 g h 2 ] x = 0 , (3)

where h is the depth of the fluid, u is the horizontal velocity and g is acceleration of gravity; and full Euler equations for an ideal gas with constant specific heats,

[ ρ ρ u E ] t + [ ρ u ρ u 2 + p u ( E + p ) ] x = 0 , (4)

where ρ , u , E , p are the density, velocity, total energy and pressure, respectively. The system of Euler equations is closed by the equation of state for an ideal polytropic gas

E = p γ 1 + 1 2 ρ u 2 ,

here γ = 1.4 is the ratio of specific heat.

The last example is 2D Euler equations,

[ ρ ρ u ρ v E ] t + [ ρ u ρ u 2 + p ρ u v u ( E + p ) ] x + [ ρ v ρ u v ρ v 2 + p v ( E + p ) ] y = 0 , (5)

where the above system is closed by the equation of state,

E = p γ 1 + 1 2 ρ ( u 2 + v 2 ) .

For ease of presentation, we consider the case of mesh refinement by a factor θ . Figure 1 shows a grid consists of two uniform subgrids and the factor θ = 2 . Suppose the problem is on the physical domain [ a , b ] , which is subdivided into cells C 1 , , C N , C N + 1 , , C N + M with x 1 2 = a and x N + M + 1 2 = b , so that

Δ x 1 = = Δ x N = b a N + M θ and Δ x N + 1 = = Δ x N + M = b a N + M θ θ . Apparently, the interface between the coarse and fine grids is located at x N + 1 2 .

3. Numerical Methods

The equations presented in the previous section will be solved using first order accurate LLF [7] method. The scheme for solving Equations (1)-(5) can be written in the form

U i n + 1 = U i n Δ t Δ x ( F i + 1 / 2 n F i 1 / 2 n ) , (6)

where F j + 1 / 2 n is the numerical flux function at the interface between cells i and i + 1 at time level n . Depending on the choice of the numerical flux formula, Equation (6) is referred to as LLF scheme.

LLF scheme is the improved version of the classical Lax-Friedrichs scheme by replacing the value a = Δ t Δ x by a locally determined value,

F i + 1 / 2 = 1 2 [ F ( U i ) + F ( U i + 1 ) a i + 1 / 2 ( U i + 1 U i ) ] (7)

Figure 1. Interface between the coarse and fine grids.

where a i + 1 / 2 = max ( | λ ( U i ) | , | λ ( U i + 1 ) | ) .

4. Numerical Result

4.1. Advection Equation

Firstly, advection Equation (1) is considered on a non-uniform grid which is composed of two uniform subgrids with initial condition

u ( x , 0 ) = { 1, if 2 x 1, 0, if 1 x 1.

The ratio between the coarse and fine grids is θ = 2 , and the interface in the physics domain is located at x = 0 . The boundary of domain a = 2 , b = 1 , and the number of cell N = 200 , M = 200 . Obviously, the discontinuity will arrive the interface when t = 1 and wholly pass through it when t = 1.5 .

For the linear advection equation, LLF scheme reduces to upwind scheme. The solutions to advection equation can be seen in Figure 2. From this figure, we find that the solution has monotone profile as expected and no oscillation can be observed around x = 0 .

4.2. Inviscid Burgers Equation

The second example is the inviscid Burgers Equation (2) subject to the initial data

u ( x , 0 ) = { 2, if 2 x 1, 0, if 1 x 1.

The solution propagates to the right with shock speed s = 1 . Figure 3 shows the solution at t = 1.5 , and the shocks pass through the interface at t = 1 . The zoomed version around x = 0 is shown in Figure 3. As in the previous problem, LLF scheme shows monotone solutions which are similar to the solutions on the uniform grids.

4.3. Shallow Water Equations

Consider the shallow water equations (3) with the piecewise-constant initial conditions

h ( x , 0 ) = { 3 , if 5 x 0 , 1 , if 0 x 5 , u ( x , 0 ) = 0. (8)

This is a special case of the Riemann problem in which u l = u r = 0 , and is called the dam-break problem because it models what happens if a dam separating two levels of water bursts at time t = 0 . This is the shallow water equivalent of the shock tube problem of gas dynamics.

Note that the Jacobian matrix F ( u ) of the shallow water equations is

F ( u ) = [ 0 1 u 2 + g h 2 u ] .

The eigenvalues of F ( u ) are

Figure 2. Left: numerical solution at t = 1.5 with LLF scheme. Right: zoomed region around the interface between the coarse and fine grids.

Figure 3. Left: numerical solution of Burgers equation at t = 1.5 with LLF scheme. Right: zoomed region around the interface x = 0 between the coarse and fine grids.

λ 1 = u g h , λ 2 = u + g h ,

with the corresponding eigenvectors

r 1 = [ 1 u g h ] , r 2 = [ 1 u + g h ] .

Firstly, a uniform grid is considered. Figure 4 shows the structure of solution of this Riemann problem. Note that the shallow water equations are a system of two equations, and so the Riemann solutions generally contain two waves. For this initial-value problem (3) (8), these always consist of one 1-rarefaction wave on the left and one 2-shock on the right with shock speed s 1.6226 . In this case, we find that the solutions behave well and no postshock oscillation is introduced.

For the shallow water equations, Leveque [8] gives detailed formula of Hugoniot

Figure 4. Left: solution of the dam-break Riemann problem for the shallow water equations with u l = u r = 0 on uniform grids; Right: zoomed version around x = 3 .

loci and integral curves. The Hugoniot curves correspond to the 1-shock and 2-shock are

h u = u l h g 2 ( h l h h h l ) ( h l h )


h u = u r + h g 2 ( h r h h h r ) ( h r h )


The integral curves correspond to 1-rarefaction and 2-rarefaction waves are

h u = h u l + 2 h ( g h l g h )


h u = h u r 2 h ( g h r g h )


As shown in previous paragraph, the solution of this initial-value problem contains both shock and rarefaction. So, Hugoniot and integral curves originate from the left state U r and right state U l , and the intersection of them is intermediate state U m . Figure 5 shows that numerical solution by LLF scheme almost lies on the Hugoniot and integral curves.

In Figure 6, physics domain [ 5,5 ] is divided into two subdomain [ 5,3 ] and [ 3,5 ] , and each subdomain is distributed some sort of uniform grids. Denote the length of the cells on interval [ 5,3 ] and [ 3,5 ] are Δ x l and Δ x r respectively. Figure 6 shows the solution on nonuniform grids with Δ x l / Δ x r = 2 at t = 2.5 . Amplifying the figure around interface at x = 3 , an obvious oscillation can be observed behind the shock (on the left of the interface). Also, we find that the shock is in the correct location and the shock speed remain unchanged

Figure 5. Left: Hugoniot loci and integral curve of the eigenvector r 1 for the shallow water equations in the state space ( h , h u ) ; Right: zoomed version.

Figure 6. Solution of the dam-break Riemann problem for the shallow water equations with nonuniform grids which consist of two uniform grids with Δ x l / Δ x r = 2 .

after it pass through the interface. The oscillation is the wave that arose from the initial discontinuity at x = 3 when the shocks pass the interface.

A interesting phenomenon can be observed in Figure 7. The hump (oscillation) in Figure 6 changes into dip if we shift Δ x l / Δ x r from 2 into 0.5. Note that the hump and dip both belong to acoustic wave moving at speed u m c m 0.6148 , where c m = g h m .

Note that, for the shallow water equations, the function of 1-Riemann invariant and 2-Riemann invariant are

w 1 = u + 2 g h


w 2 = u 2 g h


Figure 8 shows the Riemann invariants of the equations. Generally, the function of 1-Riemann invariant takes the same value at all points on the integral curve of r 1 . Similarly, the value of 2-Riemann invariant is constant along any integral curve of r 2 . We find that the acoustic wave (oscillation behind the shock) takes the same value as 1-rarefaction wave for 1-Riemann invariant.

Figure 9 shows Hugoniot, integral curves and numerical orbit in state space on nonuniform grids. Note that the numerical orbit deviating the integral curve in Figure 9 is due to the departure between numerical solution and exact solution at the left of the shock. In addition, a redundant tip appears in this figure which is different from the distribution of numerical orbit in Figure 5.

4.4. 1D Euler Equations

Euler equations of gas dynamics is discussed in this section. We test the postshock oscillations with initial data

( ρ , v , p ) = { ( 1 , 0 , 1 ) , if x 0 , ( 0.125 , 0 , 0.1 ) , if x > 0. (9)

This is a well known test problem proposed by Sod [9] . The physical domain

Figure 7. Solution of the dam-break Riemann problem for the shallow water equations with nonuniform grids which consist of two uniform grids with Δ x l / Δ x r = 0.5 .

Figure 8. The numerical orbit of the Riemann invariants for the shallow water equations. Left: 1-Riemann invariant; Right: 2-Riemann invariant.

Figure 9. Left: Hugoniot loci and integral curve of the eigenvector r 1 for the shallow water equations in the state space ( h , h u ) ; Right: zoomed version.

Figure 10. Numerical solution of Sod problem on nonuniform grids at t = 2.0 .

[ 5,5 ] is also divided into two uniform subdomain [ 5,3 ] and [ 3,5 ] .

Except for the first and third characteristic fields are genuinely nonlinear which have similar behavior to the two characteristic fields in the shallow water equations, the second field is linearly degeneration. So, the solution to Riemann problem for Euler problem typically has two nonlinear waves and a contact discontinuity.

Figure 10 shows the numerical solution on nonuniform grids with Δ x l / Δ x r = 4 . The shock has already passed through the interface at t = 2.0 , and two humps arise near x = 3 . This is different from what we observed in shallow water equations. In that system, there is just one wave arises when the shock passes through the interface between the coarse and fine grids, but here two waves arise for Euler equations. As mentioned previously, the Riemann problem for shallow water and Euler equations generally consist of two and three waves respectively. So, the phenomena that postshock oscillations arose when the shock through the interface is similar to the solution of Riemann problem.

From Figure 10, we find that two waves propagate with different direction. Note that these two waves are different from some nonphysical oscillations. The left wave propagate to the left with the speed u * c * = 0.3367 , and the left wave propagate to the other direction with the speed u * = 0.9275 , where asterisk denotes intermediate state and c * = γ p * / ρ r * . In addition, the amplitude of vibration is decided by the ratio Δ x l / Δ x r . The bigger ratio between the coarse and fine grids, the bigger amplitude the postshock will have.

For this initial value problem, all of the Riemann invariants for a polytropic ideal gas are summarized below:

1 -Riemann invariants , s , u + 2 c γ 1 , 2 -Riemann invariants , u , p , 3 -Riemann invariant , u 2 c γ 1 .

Figure 11 shows the 1-Riemann invariants for this initial value problem. We

Figure 11. 1-Riemann invariants for Euler equations with initial data (9).

Figure 12. 2-Riemann invariants.

find that the 1-Riemann invariants are constant across the rarefaction and the left wave which generate from x = 3 , and only one oscillation was left. Similarly, Figure 12 shows the 2-Riemann invariants which are constant across the second oscillation at the right of x = 3 . And 3-Riemann invariant is given in Figure 13. From these three figures, the oscillations behind the shock behave similarly to the rarefaction for Riemann invariants.

4.5. 2D Euler Equations

The last example is 2D Euler equations (5), with initial data

( ρ , u , v , p ) = { ( 1.1 , 0.0 , 0.0 , 1.1 ) , if x > 0.5 , y > 0.5 , ( 0.5065 , 08939 , 0.0 , 0.35 ) , if x < 0.5 , y > 0.5 , ( 1.1 , 0.8939 , 0.8939 , 1.1 ) , if x < 0.5 , y < 0.5 , ( 0.5065 , 0.0 , 0.8939 , 0.35 ) , if x > 0.5 , y < 0.5 , (10)

which corresponds to the case of left forward shock, right backward shock, upper backward shock, and lower forward shock. We refer the readers to [10] [11] [12] for details. This example can be used to test the postshock oscillations as shock across the interface.

Figure 14 shows the cross-sectional solution at x = 0.1 . One solution on uniform grids and another solution on nonuniform grids are shown in this figure. The obvious difference between these two solutions is that two additional humps appear in the solution under nonuniform grids. This phenomenon is similar to what appear in the previous nonlinear system.

5. Conclusion

The postshock oscillations appeared in the nonlinear system has been elucidated numerically. We found that these oscillations would not appear in the nonlinear scalar equation and linear equations system. Although all the numerical results were computed with LLF method, these oscillations also appeared with other schemes, such as Godunov, Roe schemes and so on. So the postshock oscillations

Figure 13. 3-Riemann invariant.

Figure 14. Cross section at x = 0.1 . Left: the amplified figure.

have nothing to do with the specific methods which are used to compute solution on nonuniform grids. In addition, this kind of oscillations is different from those nonphysical oscillations arose near the shock, and they are similar to those emerge behind the slow-moving shocks which can be found in [5] [6] .


This work was supported by the National Natural Science Foundation of China (No. 11501238), Natural Science Foundation of Guangdong Province (No. 2012A030313119) and Supported by the Major Project Foundation of Guangdong Province Education Department (No. 2014KZDXM070).

Cite this paper: Hu, F. (2017) Postshock Oscillations on Non-Uniform Mesh. Journal of Applied Mathematics and Physics, 5, 481-493. doi: 10.4236/jamp.2017.52043.

[1]   Berger, M.J. and Colella, P. (1989) Local Adaptive Mesh Refinement for Shock Hydrodynamics. Journal of Computational Physics, 82, 64-84.

[2]   Berger, M.J. and Oliger, J. (1984) Adaptive Mesh Refinement for Hyperbolic Partial Differential Equations. Journal of Computational Physics, 53, 484-512.

[3]   Huang, W., Ren, Y. and Russell, R.D. (1994) Moving Mesh Partial Differential Equations (MMPDES) Based on the Equidistribution Principle. SIAM Journal on Numerical Analysis, 31, 709-730.

[4]   Stockie, J.M., Mackenzie, J.A. and Russell, R.D. (2001) A Moving Mesh Method for One-Dimensional Hyperbolic Conservation Laws. SIAM Journal on Scientific Computing, 22, 1791-1813.

[5]   Roberts, T.W. (1990) The Behavior of Flux Difference Splitting Schemes near Slowly Moving Shock Waves. Journal of Computational Physics, 90, 141-160.

[6]   Arora, M. and Roe, P.L. (1997) On Postshock Oscillations Due to Shock Capturing Schemes in Unsteady Flows. Journal of Computational Physics, 130, 1-24.

[7]   Rusanov, V.V. (1962) The Calculation of the Interaction of Non-Stationary Shock Waves and Obstacles. USSR Computational Mathematics and Mathematical Physics, 1, 304-320.

[8]   Leveque, R.J. (2002) Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, Cambridge.

[9]   Sod, G.A. (1978) A Survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws. Journal of Computational Physics, 107, 1-31.

[10]   Lax, P.D. and Liu, X.D. (1998) Solutions of Two-Dimensional Riemann Problems of Gas Dynamics by Positive Schemes. SIAM Journal on Scientific Computing, 19, 319-340.

[11]   Schulz-Rinne, C.W., Collins, J.P. and Glaz, H.M. (1993) Numerical Solution of the Riemann Problems for Two-Dimensional Gas Dynamics. SIAM Journal on Scientific Computing, 14, 1394-1414.

[12]   Tang, H.Z. and Tang, T. (2003) Adaptive Mesh Methods for One-and Two-Dimensional Hyperbolic Conservation Laws. SIAM Journal on Numerical Analysis, 41, 487-515.