Journal of Applied Mathematics

Volume 2014 (2014), Article ID 275425, 13 pages

http://dx.doi.org/10.1155/2014/275425

Research Article

Simple and High-Accurate Schemes for Hyperbolic Conservation Laws

Renzhong Feng and Zheng Wang

LMIB and School of Mathematics and Systems Science, Beijing University of Aeronautics and Astronautics, Beijing 100191, China

Received 1 May 2013; Revised 4 January 2014; Accepted 6 January 2014; Published 2 March 2014

Academic Editor: Vit Dolejsi

Copyright © 2014 Renzhong Feng and Zheng Wang. 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.

The paper constructs a class of simple high-accurate schemes (SHA schemes) with third order approximation accuracy in both space and time to solve linear hyperbolic equations, using linear data reconstruction and Lax-Wendroff scheme. The schemes can be made even fourth order accurate with special choice of parameter. In order to avoid spurious oscillations in the vicinity of strong gradients, we make the SHA schemes total variation diminishing ones (TVD schemes for short) by setting flux limiter in their numerical fluxes and then extend these schemes to solve nonlinear Burgers’ equation and Euler equations. The numerical examples show that these schemes give high order of accuracy and high resolution results. The advantages of these schemes are their simplicity and high order of accuracy.

1. Introduction

In designing numerical schemes of very high order of accuracy for solving hyperbolic conservation laws one faces at least three major difficulties. One of them concerns the preservation of high accuracy in both space and time for multidimensional problems. Another one concerns conservation; this is mandatory in the presence of shock waves. The other very important issue relates to the generation of spurious oscillations in the vicinity of strong gradients. According to Godunov’s theorem [1] these are unavoidable by linear schemes of accuracy greater than one. Each one of these difficulties is in itself not easy to resolve; the simultaneous resolution of all three difficulties above represents a formidable task in the numerical analysis of hyperbolic conservation laws. At present, there are various approaches for constructing numerical schemes that attempt to overcome the above difficulties. State-of-the-art very high order methods (at least third order) for hyperbolic conservation laws include the class of ENO/WENO schemes [2–4], Spectral Method [5], the class of Compact Difference Methods [6], Discontinuous Galerkin Finite Element Methods [7], and ADER Methods [8]. These schemes meet the requirement of very high order spatial accuracy; matching time accuracy to space accuracy, however, remains an issue in all of the above approaches. As regards the ENO/WENO/MPWENO approach, the most accurate scheme reported so far uses 9th order spatial discretisation with Runge-Kutta methods for time evolution. To avoid spurious oscillations these Runge-Kutta methods must be TVD. This leads to accuracy barriers: the accuracy of such methods cannot be higher than fifth. Moreover, fourth and fifth order methods are quite complicated and have reduced stability range. In most practical implementation reported, when the solution is not smooth, a third order TVD Runge-Kutta method has been used, for example, [2]. Although increased order of spatial discretization improves accuracy for some problems, such schemes converge with third order only when the mesh is refined. Therefore when numerical simulation of compressible turbulence and wave propagation problems involves long-time evolution, it would be beneficial to use schemes which are simple and efficient and converge with higher order both in time and space.

How to construct simple and high-accurate schemes in both space and time by simple approach is an interesting work. There are various ways to construct high order numerical schemes for conservation laws. One of them is the MUSCL (Monotone Upstream-Centred Scheme for Conservation Laws) [9–11] approach which is a very popular and simple one to construct explicit, fully discrete schemes that lies on the solution of Riemann problem and could be extended to nonlinear systems. The MUSCL-type methods based on piecewise linear, local reconstruction of data are 2nd order of accuracy and their available highest order is 3rd for linear advection equation. In this paper, we generalize the MUSCL approach by modifying the piecewise linear reconstruction with a gradual slope, solving a Riemann problem at each interface and using Lax-Wendroff approach to compute the numerical flux. The new schemes constructed can reach 3rd order of accuracy in time and space, and its available highest order of accuracy is 4th. The flux limiter approach is used to obtain TVD version of the new schemes to avoid unphysical oscillations in the vicinity of large gradients.

The rest of the paper is organized as follows. In Section 2 we describe the construction of the class of schemes SHA for linear advection equation and prove its convergent order; next we set the limiter to the fluxes of these schemes to make them TVD ones in Section 3 and carry out some numerical examples in Section 4 and then extend them to solve nonlinear Burgers' equation in Section 5 and carry out some numerical examples in Section 6 and finally apply the schemes to solve one-dimensional Euler equations in Section 7. The conclusions are drawn in Section 8.

2. The Construction of Simple High-Accurate Schemes (SHA Schemes) for Advection Equations

In this section we consider the Initial Boundary Value Problem (IBVP) for the linear advection equation in the domain on the plane. This consists of a PDE together with initial condition (IC) and boundary conditions (BCs); namely,

with suitable boundary conditions.

A mesh for discretizing the domain is a regular grid of dimension (spacing of grid points in the -direction) by (spacing in the -direction). We suppose that is discretized by equally spaced grid points. Then

The mesh points of the mesh are positioned at on the plane, with and . We will use the notation , . The discrete values of the function at will be denoted by

We will also use the symbol to denote an approximation to the exact mesh value . The distinction will be made at appropriate places.

We could construct a class of simple high-accurate schemes (SHA schemes) for advection equation as follows.

Step 1. Data reconstruction with boundary extrapolated values

We choose , , where

Step 2. Consider evolution of , by a time ( according to

Step 3. Solution of piecewise constant data Riemann problem

This step solves the conventional Riemann problem at the interface with data

the solution of which for linear advection equation (1) is

where . Then applying the solution (9) to the Lax-Wendroff flux for linear advection equation

we get a new flux

After some algebra operations, becomes

Also, we can substitute (9), the solution of (7), into the following type of flux;

or two-step Richtmyer type flux of method

to get the flux as (12), where , .

Conservative Scheme. Substitution of as numerical flux into conservative method gives a new scheme:

where

As for linear advection equation , we have

The following introduces a theorem on how to verify the accuracy of a linear scheme:

Theorem 1 (Roe [12]). A scheme of the form (17) is order accurate in space and time if and only if

According to Roe’s theorem, it is easy to verify that scheme (16) is third order accurate in space and time for any value and the scheme is fourth order accurate when .

3. TVD Version of SHA Scheme

According to Godunov’s Theorem [1] there are no monotone, linear schemes (17) of second or higher order of accuracy. This means that any linear scheme of second or higher is nonmonotone and there will be spurious oscillation in the vicinity of high gradients. The scheme must be modified in order to avoid spurious oscillation near the high gradients. Harten [13] proposed the TVD property of conservation schemes. The numerical solutions will be free from spurious oscillations if the scheme has the TVD property. Now we will modify the SHA method above to a scheme with the TVD property according to Harten’s Theorem.

Considering the class of nonlinear schemes

where and the coefficients , are in general assumed to be functions of the data , Harten had proved the following important result.

Theorem 2 (Harten [13]). For any scheme of the form (19) to solve (1), a sufficient condition for the scheme to be TVD is that the coefficients satisfy

The conservative scheme SHA can be written as follows when applied to linear advection equation:

where

Thus we obtain

or

Let

or

Then we could modify (23) as

or

where is a function

Obviously ; schemes (27a) and (27b) modified above are TVD according to Harten’s result. In addition formula (27a) can be rewritten into

Formula (27b) also can be done similarly. So it can be seen that the conservativity of the TVD schemes (27a) and (27b) is not destroyed.

4. Numerical Experiments for Linear Advection Equation

4.1. Accuracy Test for SHA Scheme

In the subsection we first test the order of accuracy of scheme SHA (16). Let us consider the linear equation subject to periodic initial data:

This problem admits the global smooth solution that was computed at time with the varying number of grid points . In Table 1 the accuracy of the scheme is shown. The table lists the errors at unit using the fourth order scheme SHA (16) (); notice that we have the full order of accuracy, 4, in both and norms.

Table 1: Initial value problem (30) and and norms of the errors.

-error
Order
-error
Order

10

—

—

20

3.957

3.934

40

3.985

3.984

80

3.996

3.996

160

3.999

3.999

4.2. Numerical Experiment I

Next we will use the fourth order scheme SHA (16) () and TVD schemes (27a) and (27b) to compute the continuous initial value test problem (31) and discontinuous initial value test problem (32) for the linear advection (1). Computed results are shown, respectively, in Figures 1, 3, and 5 and Figures 2, 4, and 6. The first test problem to be solved is

where .

Figure 1: SHA method applied to linear advection equation with a continuous initial value; numerical (symbols) and exact (line) solutions are compared at units.

Figure 2: SHA method applied to linear advection equation for a square wave; numerical (symbols) and exact (line) solutions are compared at units.

Figure 3: SHA method with TVD (27a) applied to linear advection equation with a continuous initial value; numerical (symbols) and exact (line) solutions are compared at units.

Figure 4: SHA method with TVD (27a) applied to linear advection equation for a square wave; numerical (symbols) and exact (line) solutions are compared at units.

Figure 5: SHA method with TVD (27b) applied to linear advection equation with a continuous initial value; numerical (symbols) and exact (line) solutions are compared at units.

Figure 6: SHA method with TVD (27b) applied to linear advection equation for a square wave; numerical (symbols) and exact (line) solutions are compared at units.

In this experiment the Courant number is . Figure 1 shows the excellent behavior of scheme (16) after long time. Figures 3 and 5 show the results computed by scheme (27a) and scheme (27b), respectively, at units and reveal that schemes ((27a)-(27b)) also have high approximate accuracy to smooth initial problem.

4.3. Numerical Experiment II

We solve the same problem as in experiment I, except that we replace the smooth initial condition with a discontinuous initial condition:

where .

In experiment II the is still 0.9 and . Figure 2 shows the result computed by the 4th order scheme SHA (16) () at units and reveals the 4th order scheme SHA without TVD products spurious oscillation in the vicinity of discontinuity. Figure 4 shows the result computed by the scheme SHA with TVD (27a) at units and Figure 6 shows the result computed by scheme (27b) at units. From the results shown by Figures 4 and 6 we see that schemes (27a) and (27b) both control the spurious oscillation, but the comparison of Figures 4 and 6 reveals that scheme (27b) has better resolution to discontinuous initial problem. The period boundary conditions are used in both experiments I and II.

5. The Application of SHA Scheme to Burgers’ Equation

We can see the SHA schemes perform well in solving linear advection equation from the numerical results above. It can reach fourth order accuracy in space and time where the solution is smooth enough and the implementation of the scheme is so easy that the method can be extended to nonlinear conservative laws conveniently.

In the following we will extend the schemes to inviscid Burgers’ equation

with suitable boundary conditions.

5.1. Number and Numerical Boundary Conditions

Assume that the space is divided into cells ; the spacing of grid points is , which is . The full discrete conservative formula for the above equation is

In the application of SHA scheme to the inviscid Burgers’ equation we use the following transmissive boundary conditions:

to provide numerical fluxes , . Denoting by the maximum wave speed throughout the domain at time level we define the maximum Courant number :

where is such that

The time step follows as

where

The calculation of the max wave speed value in (39) needs to use all of the speed values in the cells, including the boundary cells.

The Riemann problem for Burgers’ equation is

Its exact solution is

5.2. SHA Schemes for Solving Nonlinear Scalar Equation

The high order scheme is constructed as follows.

Step 1. Data reconstruction

is the slope on cell for , and

Then we get the values on the interface of cells:

Step 2. Consider evolution of , by a time ( according to

Step 3. Solution of piecewise constant data Riemann problem

Step 4. Compute the integral

Then the numerical flux is computed by .

The is required for computing the numerical flux ; that means the integral of the Riemann solution needs to be calculated. The following gives the exact integral with exact Riemann solution of inviscid Burgers’ equation. There are two cases.

Case 1. See Figure 7, ; the shock wave is formed; the exact solution for Riemann problem is

where . Substitute (48) into (47); then

Figure 7: shock wave.

Case 2. See Figure 8, , and the solution is a centred rarefaction wave; the exact solution for Riemann problem is

Substitute the above solution (50) into (47), then

We can get the numerical flux at th cell interface. As we can see, the flux is the same as two-step Richtmyer flux versions (13a) and (13b). Then we can give the conservative formula for solving Burgers’ equation:

Figure 8: rarefaction wave.

5.3. TVD Version of SHA for Burgers’ Equation

In this section we will construct the TVD version of SHA (52) by the combination of a low order monotone flux and a high order flux as

where is a flux limiter. Here a specific Flux Limiter Centred scheme is constructed by taking

where is a First Order Centred scheme and is given as follows for nonlinear systems.

A classical scheme of first order accuracy to solve hyperbolic conservation laws is that of Lax-Friedrichs, whose numerical flux at the interface of two states , is

A second order accurate scheme Richtmyer scheme is given by computing a numerical flux

where the intermediate state

Hence is given as

As we know in [14], the relationship between conventional upwind flux limiters and centred flux limiters is

with

where is the maximum Courant number. By choosing the upwind flux limiters SUPERBEE [15], VANLEER [16], VANALBADA [17], and MINBEE [15] we can obtain corresponding centred flux limiters. Here we recommend the use of the SUPERBEE flux limiter to obtain a centred flux limiter:

Set

A single flux limiter is computed as

Here is the limiter function (61), and we will apply to the flux in (53). As we see, the method can be easily extended to a system of conserved variables; we just need to apply (62)-(63) to every component and obtain a corresponding limiter , . Then, one can select the final limiter as

In next section, we present some numerical results for above methods.

6. Numerical Experiment for Inviscid Burgers’ Equation

Apply the SHA scheme (52) to the inviscid Burgers’ equation as follows:

the exact solution of which is .

Firstly, we test the accuracy of the scheme when the solution is still continuous. and errors are listed in Table 2 using SHA scheme (52) at units. We have the order of accuracy, 2, in both and norms.

Table 2: Initial value problem (65) and and norms of the errors.

-error
Order
-error
Order

10

—

—

20

2.076

1.885

40

1.909

1.836

80

2.025

2.021

160

2.007

2.001

Next we compute the shock by SHA scheme (52) with and . The comparison of exact solution and numerical solution for initial value problem (65) is given in Figure 9. Then we use the SHA scheme (53) with the limiter in (63) to solve the initial value problem (65). The comparison is shown in Figure 10. From Figures 9 and 10 we can see that SHA scheme gives satisfactory solution when shock develops at .

Figure 9: SHA method applied to inviscid Burgers’ equation for continuous initial value; numerical (symbols) and exact (line) solutions are compared at units.

Figure 10: SHA method with TVD applied to inviscid Burgers’ equation for continuous initial value; numerical (symbols) and exact (line) solutions are compared at units.

Considering the discontinuous initial value problem for Burgers’ equation

The comparison of exact solution and numerical solution for initial value problem (66) computed by SHA scheme (52) is given in Figure 11; then we use the SHA scheme (53) with TVD to solve the initial value problem (66); the comparison is shown in Figure 12. From Figures 11 and 12 we see that SHA scheme (52) cannot control the spurious oscillation, but SHA scheme (53) not only can control the spurious oscillation but also has better resolution to discontinuity.

Figure 11: SHA method applied to inviscid Burgers’ equation for discontinuous initial value; numerical (symbols) and exact (line) solutions are compared at units.

Figure 12: SHA method with TVD applied to inviscid Burgers’ equation for continuous initial value; numerical (symbols) and exact (line) solutions are compared at units.

7. The Application of SHA Scheme to Euler Equations

We have presented every procedure of SHA scheme applying to Burgers’ equation; now we will turn our focus on the application of SHA scheme to the following Euler equations:

with

where

and is conservative variable, is vector of convective flux, is specific internal energy, is total energy, and is ratio of specific heats, the value of which is 1.4 for ideal gas.

The process of applying SHA to Euler equations is similar to that of Burgers’ equation. The time step is determined by

where is the maximum wave speed at the time . A practical choice of is

where is the sound speed, for ideal gas .

The data reconstruction and evolution part is the same as Burgers’ equation. When we obtain the and at each interface, there are two ways to construct the flux for each interface.

The Riemann solver free way is to use two-step Richtmyer version of Lax-Wendroff method to obtain the SHA numerical flux:

where , . Then (53) can be used to construct TVD scheme for the Euler equations which is similar to the procedure for Burgers’ equation.

The other way to obtain the flux at the cell is according to (13a); for the Euler equations, the flux is given as

where is the solution of the Riemann problem with piecewise constant data , . And approximate Riemann solvers may be applied to get the solution of the Riemann problem. Assume the four constant states of the solution are , , , in order; in that way, the integral in (73) becomes [14]

where is the value of in region and is the Courant number for wave of speed . And then the flux is obtained from (73). The TVD version of the flux of the scheme is constructed as follows:

where , and , for ; , for . The quantity is known to change across each wave family in the solution of the Riemann problem. For the Euler equations the choice (density) gives very satisfactory results. And a recommended limiter can be found in [14]:

Finally, the TVD flux is .

In order to test the effectiveness of the scheme, we apply the TVD schemes (72) and (53) (denoted as SHA_1) and the TVD schemes (73), (75), and (76) (denoted as SHA_2) for Euler equations to Sod's tube problem and Lax problem [18, 19]. HLLC approximate Riemann solver [20] is used to obtain the solution of the Riemann problem at the interface. At the same time, the MUSCL-Hancock scheme [21], WEN (3rd order) [22, 23], and RKDG (3rd order) [24] are also applied for comparison.

Sod’s Tube Problem. The initial value consists of two constant states:

This is a very common test for numerical scheme and its solution consists of a left rarefaction, a contact, and a right shock. The mesh used is and the Courant number is 0.9. The comparisons between the numerical solution obtained by the five schemes and the exact solution are shown in Figures 13, 14, and 15. We can see that both kinds of SHA schemes with TVD give satisfactory results. The numerical results from the SHA_1 without Riemann solver are similar to that of MUSCL-Hancock methods, and both methods give more dissipation than the other three methods, and SHA_2 with Riemann Solver scheme gives least dissipation at the contact wave and shock wave.

Figure 13: Density comparison of schemes at time 0.2 units for Sod’s tube problem.

Figure 14: Velocity comparison of schemes at time 0.2 units for Sod’s tube problem.

Figure 15: Pressure comparison of schemes at time 0.2 units for Sod’s tube problem.

Lax Problem. The initial states are given by

This is a tough test for numerical schemes. The Courant number in this test is still 0.9 and . The numerical results computed by the five schemes are shown in Figures 16, 17, and 18 comparing with the exact solution. We can see that the numerical results from both SHA schemes are good, the results of SHA_2 are comparable with results obtained by WENO or DG, but SHA_2 is more simple and requires less computational cost. These show the effectiveness of our method.

Figure 16: Density comparison of schemes at time 0.16 units for Lax problem.

Figure 17: Velocity comparison of schemes at time 0.16 units for Lax problem.

Figure 18: Pressure comparison of schemes at time 0.16 units for Lax problem.

8. Conclusion

This paper introduces a class of SHA schemes with 3rd order approximation accuracy only for linear advection equation by linear data reconstruction, solving Riemann problem and using Lax-Wendroff flux, and its available highest order is 4th. The general flux limiter approach is used to make the SHA schemes TVD in order to avoid spurious oscillations in the vicinity of strong gradients for nonlinear equations and systems. The numerical results show that these schemes can solve hyperbolic conservation laws with high order of accuracy, high resolution, and less computational cost. The next work is to extend SHA schemes to high-dimensional hyperbolic equations and apply them to compute some practical flows, for example, around an airfoil and wing.

Conflict of Interests

The authors declare that there is no conflict of interests concerning the publication of this paper.

Acknowledgments

The authors thank the reviewers very much for giving many valuable comments on their paper. The work was supported by the National Natural Science Foundation of China (no. 11271041), National Natural Science Key Foundation of China (no. 10931004), ISTCP of China Grant (no. 2010DFR00700), and CASP of China Grant (no. MJ-F-2012-04).

References

S. K. Godunov, “A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics,” Matematicheskii Sbornik, vol. 47–89, no. 3, pp. 271–306, 1959.

D. S. Balsara and C. W. Shu, “Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy,” Journal of Computational Physics, vol. 60, no. 2, pp. 405–452, 2000.

G. S. Jiang and C. W. Shu, “Efficient implementation of weighted ENO schemes,” Journal of Computational Physics, vol. 126, no. 1, pp. 202–228, 1996.

A. Harten, B. Engquist, S. Osher, and S. R. Chakravarthy, “Uniformly high order accurate essentially non-oscillatory schemes, III,” Journal of Computational Physics, vol. 71, no. 2, pp. 231–303, 1987.

C. Canuto and A. Quarteroni, Eds., Spectral and Higher Order Methods for Partial Differential Equations, North-Holland, Amsterdam, The Netherlands, 1990.

A. I. Tolstykh, High Accuracy Non-Centered Compact Difference Schemes for Fluid Dynamics Applications, World Scientific Publishing, Hackensack, NJ, USA, 1994.

B. Cockburn and C. W. Shu, “The runge-kutta discontinuous Galerkin method for conservation laws V: multidimensional systems,” Journal of Computational Physics, vol. 141, no. 2, pp. 199–224, 1998.

E. F. Toro and V. A. Titarev, “TVD fluxes for the high-order ADER schemes,” Journal of Scientific Computing, vol. 24, no. 3, pp. 285–309, 2005.

B. van Leer, “Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow,” Journal of Computational Physics, vol. 23, no. 3, pp. 263–275, 1977.

B. van Leer, “Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection,” Journal of Computational Physics, vol. 23, no. 3, pp. 276–299, 1977.

B. van Leer, “Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method,” Journal of Computational Physics, vol. 32, no. 1, pp. 101–136, 1979.

P. L. Roe, “Numerical algorithms for the linear wave equation,” Tech. Rep. 81047, Royal Aircraft Establishment, Bedford, UK, 1981.

A. Harten, “High resolution schemes for hyperbolic conservation laws,” Journal of Computational Physics, vol. 49, no. 3, pp. 357–393, 1983.

E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, New York, NY, USA, 1997.

P. L. Roe, “Some contributions to the modelling of discontinuous flows,” in Proceedings of the SIAM/AMS Seminar, San Diego, Calif, USA, 1983.

B. van Leer, “Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme,” Journal of Computational Physics, vol. 14, no. 4, pp. 361–370, 1974.

G. D. van Albada, B. van Leer, and W. W. Roberts, “A comparative study of computational methods in cosmic gas dynamics, astron,” Astrophysics, vol. 108, pp. 76–84, 1982.

G. A. Sod, “A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws,” Journal of Computational Physics, vol. 27, no. 1, pp. 1–31, 1978.

B. Thornber, D. Drikakis, R. J. R. Williams, and D. Youngs, “On entropy generation and dissipation of kinetic energy in high-resolution shock-capturing schemes,” Journal of Computational Physics, vol. 227, no. 10, pp. 4853–4872, 2008.

E. F. Toro, M. Spruce, and W. Speares, “Restoration of the contact surface in the HLL-Riemann solver,” Shock Waves, vol. 4, no. 1, pp. 25–34, 1994.

B. van Leer, “On the relation between the upwind-differencing schems of Godunov, Enguist-Osher and Roe,” SIAM Journal on Scientific Computing, vol. 5, no. 1, pp. 1–20, 1985.

X. D. Liu, “Weighted essentially non-oscillatory schemes,” Journal of Computational Physics, vol. 115, no. 1, pp. 200–212, 1994.

G.-S. Jiang and C. W. Shu, “Efficient implementation of weighted ENO schemes,” Journal of Computational Physics, vol. 126, no. 1, pp. 202–228, 1996.

B. Cockburn and C. W. Shu, “TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general frame-work,” Mathematics of Compution, vol. 52, pp. 411–435, 1989.