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.
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 .
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
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:
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.
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 .
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.
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.
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.
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).