The Finite Volume WENO with Lax–Wendroff Scheme for Nonlinear System of Euler Equations

We develop a Lax–Wendroff scheme on time discretization procedure for finite volume weighted essentially non-oscillatory schemes, which is used to simulate hyperbolic conservation law. We put more focus on the implementation of one-dimensional and two-dimensional nonlinear systems of Euler functions. The scheme can keep avoiding the local characteristic decompositions for higher derivative terms in Taylor expansion, even omit partly procedure of the nonlinear weights. Extensive simulations are performed, which show that the fifth order finite volume WENO (Weighted Essentially Non-oscillatory) schemes based on Lax–Wendroff-type time discretization provide a higher accuracy order, non-oscillatory properties and more cost efficiency than WENO scheme based on Runge–Kutta time discretization for certain problems. Those conclusions almost agree with that of finite difference WENO schemes based on Lax–Wendroff time discretization for Euler system, while finite volume scheme has more flexible mesh structure, especially for unstructured meshes.


Background
In this paper, a new scheme is proposed using Lax-Wendroff scheme to discretize the time derivative and finite volume WENO (Weighted Essentially Non-oscillatory) scheme to discretize the spatial derivative of PDEs (partial differential equations). The scheme is adopted to simulate one-dimensional and two-dimensional nonlinear Euler systems of compressible gas dynamics; in fact, the scheme and results can be generalized to other hyperbolic conservation laws easily.

Formulation of the Problem of Interest for This Investigation
It has been an important work to study nonlinear PDE [1] due to their rich mathematical structures and features [2][3][4][5] as well as important applications in fluid dynamics and plasma physics [6][7][8][9][10][11][12]. Although many theories and methods were proposed to discuss the PDE [13][14][15][16][17][18][19][20], however, most nonlinear PDE have no analytic solutions; numerical methods are necessary to study hydrodynamic characteristics of PDEs [21][22][23][24]. For the discretization of time derivative in time-dependent PDEs for numerical methods, the well-known ways is the ODE solver, such as multi-step method and Runge-Kutta method. These approaches are popular for their simplicity in concept and in coding, and better stability properties. Mostly, a total variation diminishing (TVD) time discretization is preferred.
Of course, there also exist disadvantages, the main disadvantage being that total variation diminishing (TVD) Runge-Kutta methods with positive coefficients cannot reach fifth order accuracy or a higher order. In many practical schemes, a third TVD Runge-Kutta method has also been used, even for a very high order in spatial. This combination gives satisfactory results as spatial accuracy is more crucial than temporal accuracy for many problems.
For the discretization of spatial derivative in PDEs for numerical methods, especially for hyperbolic conservation law, traditional finite difference methods and finite volume methods are not fit to simulate Riemann solutions. In recent years, even many numerical methods have been developed and applied in recent years, many studies have more interest on WENO scheme. WENO is one of the very important high accuracy numerical methods, as WENO provides a procedure of spatial discretization for PDE. WENO schemes are shown to be very robust, especially at those regions where solutions are discontinuities, or have sharp gradient and other complicated solution structures. The WENO schemes have been developed as a class of high order finite difference or finite volume methods for hyperbolic conservation laws [25,26] in recent years. Finite difference methods are easily programmed and are convenient for high dimensional cases, while weakness is only fit for structured mesh with regular computational domain. In fact, the areas of inshore coastal waters and rivers have irregularly shaped areas.

Literature Survey
WENO [27,28] is a numerical method [29] used to approximate the spatial derivative terms to solve PDEs. WENO schemes have been very important high accuracy numerical methods, since the first WENO scheme was constructed for the third order finite volume version based on ENO (Essentially Non-oscillatory) scheme [30]. The main idea of WENO schemes is a non-linear-weighted combination of several local reconstructions based on different stencils and the usage of it as a final WENO reconstruction. Now, many improved WENO scheme are developed, such as Hermite-WENO [31], hybrid-WENO [32], positivity-preserving WENO [33], ADER-WENO [34], CWENO [35] and so on. The finite difference and the finite volume WENO schemes are widely used in many areas [36][37][38][39][40][41]. Lu and Qiu [42] investigated the performance of finite volume WENO methods for shallow water equations on unstructured meshes and simulated the tidal bore on Qiantang River.
Another way to approximate the time derivative is Lax-Wendroff-type time discretization procedure [43]. The method is based on traditional Lax-Wendroff method [44], which is an alternative method for time discretization, referring to a Taylor expansion in time, also called Taylor type. Qiu and Shu [43] developed Lax-Wendroff time discretization procedure with finite difference WENO scheme to solve Euler system of compressible gas dynamics, which obtained an interesting conclusion of exploring a balance between reducing the cost and maintaining the non-oscillatory properties. Lu and Qiu [45] also developed Lax-Wendroff with finite difference WENO scheme for shallow water equations, and similar results were obtained.

Scope and Contribution of This Study
The Lax-Wendroff time discretization method is via the classical Lax-Wendroff procedure, which relies on the conversion that makes all the time derivatives into spatial derivatives in a temporal Taylor expansion, and by using the PDE, discrete spatial derivatives. When the two ways obtain the same high order accuracy, the Lax-Wendroff time discretization needs a smaller effective stencil than the ODE solver, and it makes full use of the original PDE. However, the formulation and coding of this method can be quite complex to understand. The Lax-Wendroff time discretization formula based on finite volume is easier than that of the finite difference, but the computation is little more complex than that of finite difference method, as cell averages have to be calculated for finite volume method, while finite difference method only needs to calculate the point value.
Thus, in this paper, we mainly discuss finite volume WENO scheme in space as it is more flexible than finite difference scheme. We follow the same ideas of the previous papers [46,47] about finite volume WENO schemes.

Organization of the Paper
We develop the fifth order WENO finite volume schemes based on the third order Lax-Wendroff-type time discretization in this paper, which is used to solve the one-dimensional and two-dimensional nonlinear systems of Euler equations to see if similar conclusions still hold. The main work of this paper is organized as follows, in Section 2, we describe details of the discretization for Euler system with finite volume WENO schemes and Lax-Wendroff-type time discretization. Numerical examples are given in Section 3. The results demonstrate the accuracy, resolution and efficiency of the constructed schemes. Concluding remarks are included in Section 4.

The Nonlinear Euler System
We introduce the nonlinear Euler system in two dimensions; it has the following form with U = (ρ, ρu, ρv, E) T ; where ρ stands for the density, E is the total energy, u, v are the velocity in x, y directions, p is the pressure, and p is related to E = p γ−1 + 1 2 ρ(u 2 + v 2 ), with γ related to specific heat ratio, for ideal gas γ = 1.4. Euler system is important and typical hyperbolic conservation laws. In the Euler system (Equation (1)), U is the conserved vector, F(U), G(U) are the flux terms. We rewrite the Euler system in a quasi-linear form as follows where A, B are Jacobian matrices, The eigenvalues of Jacobian matrices A and B are respectively, where c = γp ρ is the sound speed of the polytropic gas. The corresponding right and left eigenvectors of Jacobian matrices A and B as follows: Neglecting all the variables in y direction of the above formula, we can easily get the characteristic of one-dimensional Euler system (Equation (6)), thus the following procedure is done: with

Overview of the Fifth Order Finite Volume WENO Schemes
We use the finite volume WENO scheme to solve the spatial derivative of Euler equations. We take a brief look at WENO scheme [28] for the typical scalar one-dimensional hyperbolic conservation law Equation (8) u Denote ] the ith cell; the center is x i . We assume that the meshes are uniform and . For a finite volume scheme, we define the approximate average valueū i ); here, the simplest and monotone Lax-Friedrich flux is employed, which has the form  are computed through k neighboring mean valuesū j , based on stencil with r cells to the left on the order of accuracy k. Thus, there exist constants c kj andc kj , such that The constants c rj are obtained by the unique k − 1 degree polynomial p r (x), which is so the approximate point value is u L ).
For example, for the fifth order WENO scheme (k = 3), the r(r = 0, 1, 2) approximated point values of u L i+ 1 2 are given by The following part describes the computation of the smoothness indicators β s and the nonlinear weights ω m . The (2k − 1)th order WENO reconstruction is a convex combination based on all these k point values The nonlinear weights ω m have the following new definitions: where ε is a small positive constant to avoid zero denominator; usually, we take it as 10 −6 . d s are linear weights which yield (2k − 1)th order accuracy on area of smooth solution, β s are the so-called "smoothness indicators", which are usually used to measure the smoothness of the polynomials and they have the form: For the fifth order WENO scheme on uniform mesh, the linear weights d s have the form: and the smoothness indicators β s are given by On uniform meshes, we do not have to calculatec rj . The reconstruction of u R . For the Euler system, the reconstruction is performed in eigenspace, so local characteristic decompositions should be used. The reconstruction in eigenspace is more robust than a component by component method. We need to calculate the right eigenvectors and the left eigenvectors of the Jacobian matrices A, B and Roe's average in the local characteristic decomposition. For two-dimensional finite volume case, we cannot simply procede in a dimension-by-dimension fashion as in finite difference method. The spatial surface integration on rectangular cell can be converted to line integration by Gauss theorem. To obtain the fifth order, we can use three-point Gaussian quadrature for line integration of every rectangular cell as in Figure 1. For example, to calculate numerical flux in x direction, three Gaussian point values on every vertical lines should be reconstructed from cell average values in y direction first; then left and right Gaussian point values should be reconstructed by the WENO scheme; next is the numerical flux; and finally the cell average value can be approximate on next time level. When dissipation does not exist, we first suppose where c is a constant and denote the phase speed of the wave. Equation (8) implies that Based on the above transform, Equation (7) becomes − c ∂ ∂r where ). When we can easy infer that the relative vorticity ∇ 2 u is a function of stream function u only. Thus, we get

The Lax-Wendroff-Type Time Discretization Procedure for Euler Equations
The details of Lax-Wendroff-type time discretization are described based on finite volume WENO schemes for two-dimensional Euler system (Equation (1)) [43,45] in this section.
Define ∆t as the time step; then (n + 1)th time level t n+1 = t n + ∆t. We introduce U n i,j as the approximate average values of U on rectangular cell I i,j at the nth time level. For simplicity, we denote the rth order time derivative of U as U (r) rather than ∂ r U ∂t r , and U , U and U stand for the first, the second and the third time derivatives of U, respectively. The idea of Lax-Wendroff-type time discretization procedure is to replace time derivative by spatial derivatives using the original PDE (Equation (1)) [43]. By a temporal Taylor's expansion for conservative vector U(x, y, t + ∆t), we obtain U(x, y, t + ∆t) = U(x, y, t) + ∆tU (x, y, t) + (∆t) 2 2 U (x, y, t) + (∆t) 3 6 U (x, y, t) + (∆t) 4 24 U (4) (x, y, t) + ..., The sliding average is U(x, y, t + ∆t) = U(x, y, t) + ∆tU (x, y, t) + (∆t 2 ) 2 U (x, y, t) + (∆t) 3 6 U (x, y, t) + (∆t) 4 24 U (4) (x, y, t) + ..., where U (r) means mean value of U (r) in space. If we need to obtain kth degree of accuracy in time, we should approximate those terms until the kth time derivatives on the right hand side in Equation (17). Here, we calculate the third order in time and the fifth order in space, which means that we should approximate until the third order time derivatives on the right hand side in Equation (17), and for every time derivative the fifth order spatial approximation method should be used; actually, the fifth order WENO reconstruction is used.
Using Euler system (Equation (1)), by Lax-Wendroff type dime discretization, we replace the time derivatives in Equation (16) by the spatial derivatives as follows where brief notes are used, such as F (U) 2 = F (U)F (U), the details are given below. Thus, Equation (16) can be written as Then, the sliding average (Equation (17)) with third order in time on a rectangular cell I ij is approximated as whereF(U) andG(U) are approximated by numerical fluxF andĜ, such asF(U i+ 1 2 , * ) ≈F i+ 1 2 , * = F(U L  (21) leads to a one order lower flux approximation than that of the first time derivative U . In x-direction, the linear polynomial p r (x) and constant p r (x) are the approximations to U x and U xx . According to polynomial p r (x) (Equation (9)), on r(r = 0, 1, 2) small stencil, we have The same linear weights of U are used for U x , U xx . According to the definition of smooth, Equation (11) indicates constant β r and zero β r are used for U x and U xx , as Similar formulae can be written for U y and U yy in y-direction in WENO reconstruction. In Equation (21), the Jacobian matrices are F (U) = A(U), G (U) = B(U), which have been mentioned above. F (U) 2 U x , G (U) 2 U y are 4 × 1 vectors. If we denote where A 21 , A 22 , A 23 , A 24 are 4 × 4 matrices, which is also a 4 × 1 vector and have been calculated in a previous step, in which we known F (U) 2 (U x ) is a 4 × 1 vector, with where B 21 , B 22 , B 23 , B 24 are 4 × 4 matrices, and If a higher accuracy order in time is required, the procedure needs to be continued until the corresponding derivatives in a similar fashion. For convenience, in the following parts, WENO5-LW3 denotes a scheme coupling the fifth order finite volume WENO scheme in space and the third order Lax-Wendroff-type method in time and WENO5-RK3 denotes a scheme coupling the fifth order finite volume WENO schemes and the third order Runge-Kutta methods.

Numerical Results
In this section, extensive numerical examples are used for performance testing for the new finite volume WENO5-LW3 schemes. We compare the results of WENO5-LW3 for one-and two-dimensional Euler equations and that of the WENO5-RK3 schemes, focusing on the important issue of CPU time cost, non-oscillatory and efficiency. For comparison of CPU time cost, all computations are performed on a personal computer, Core (TM) i5-3470 CPU @ 3.20 GHz with 4 GB RAM. In these numerical simulations, we use CFL number 0.4 and 0.6 for WENO5-LW3 and WENO-RK3 schemes, respectively. The small positive constant ε = 10 −6 is used to avoid denominator zero in the nonlinear weight formula. All of these simulations are based on uniform meshes.

Example 1 (One-Dimensional Accuracy Test).
To check the order accuracy of the finite volume WENO5-LW3 schemes based on one-dimensional Euler system with a smooth solution (Equation (6)), for comparison, we use the same test as in [43]. In the computational domain [−1, 1], the initial condition is given below: Two-periodic boundary conditions are used; we simulate until the ending time t = 2, when the solutions are continuous and there are no shocks. The exact solution is ρ(x, t) = 1 + 0.2sin(π(x − t)), v(x, t) = 1.0, p(x, t) = 1.0.
The numerical errors and orders of accuracy based on L 1 and L ∞ norm metric for density ρ by WENO5-LW3 are shown in Table 1. N stands for the number of mesh grids; adjusted adaptive time step size should be used to get the fifth order as mesh refining. For comparison, the L 1 and L ∞ numerical errors and accuracy orders by WENO5-RK3 schemes are displayed in Table 1. In Table 1, we can find that the designed accuracy order are attained for both WENO5-LW3 and WENO5-RK3, and they generate similar numerical errors. In Figure 2, we give the numerical errors and CPU time cost for density ρ based on L 1 and L ∞ norm metric with cell numbers are 10 × 2 k−1 , k = 1, 2, ..., 7, and Log scale is used for numerical errors in y axis. In this figure, the line of WENO5-LW3 always underneath the line of WENO5-RK3, which implies that, with the same CPU time cost, WENO5-LW3 obtains a higher accuracy than WENO5-RK3. From the trend of the line, it seems that the phenomenon will be more serious when the test has more mesh points. Example 2 (One-Dimensional Test with Shocks). This problem is commonly applied to test the ability of shock capturing for schemes in one-dimensional Euler system [43]. In this case, we considertwo different initial conditions in the domain [−0.5, 0.5], with continuation boundary condition for two boundary; the first test is Lax Problem, (ρ, v, p) = (0.445, 0.698, 3.528) x ≤ 0 (0.5, 0, 0.571) x > 0 . (27) and the second test is Sod Problem, Both problems have Riemann initial condition. The simulation is performed up to final time t = 0.16 for the Lax Problem and ending time t = 0.1644 for the Sod Problem, when shocks have appeared in the solutions.
In this case, we put more concern on the resolution of numerical shocks. The computed density ρ based on the WENO5-LW3 and WENO5-RK3 schemes with N = 200 on uniform mesh are plotted in Figure 3 and compared with the exact solution. The numerical results and the exact solution fit very well, no spurious oscillations occur, and the results are comparable with that in [43]. In this test, we take the CPU time based on 2000 uniform mesh grids and tested for 10 times; the CPU time for WENO5-LW3 and WENO5-RK3 schemes are 23.41 s and 26.52 s, respectively. From the data, we find that the CPU time cost of WENO5-LW3 scheme is about 10% less than that of WENO5-RK3 scheme. Example 3 (One-Dimensional Test of Shocks Interaction with Entropy Waves). This problem is a typical one-dimensional test for higher order method to evaluate the non-oscillatory properties [43], as this case contains not only shocks but also complex smooth region structure; more specifically, it is a problem about the interaction with entropy waves of shocks. The given initial condition is: where = 0.2 is used. Actually, this is a problem about a moving Mach = 3 shock interacting with sine waves. The simulation is performed until the ending time t = 1.8 s with continuation boundary condition on the domain [−5, 5].
The computed density ρ based on the WENO5-LW3 schemes and WENO5-RK3 schemes with the mesh points N = 200 and N = 400 and zoomed in pictures are plotted in Figure 4 against a reference exact solution by 2000 mesh points, which show that the computed solution converges to the reference exact solution, and no spurious oscillations occur. The WENO5-LF3 scheme could get much better resolution than that of WENO5-RK3 scheme. The results are better confirmed with reference solution on N = 400 than that on N = 200, which also shows a converged solution computed by WENO5-LW3. All results agree with that in [43]. The time cost to run 10 times based on uniform mesh with N = 2000 is taken as the CPU time for every problem. The CPU time is 59.94 s for WENO5-LW3 and 66.24 s for WENO5-RK3, respectively, so the CPU time cost of WENO5-RK3 is around 8.5% more than that of WENO5-LW3.  Example 4 (Two-Dimensional Accuracy Order Test). To check the order accuracy of the schemes for two-dimensional Euler system with a smooth solution, for comparison, we use the test in [43]. In the computational domain [0, 2] × [0, 2], we take the following initial conditions: ρ(x, y, 0) = 1 + 0.2sin π(x + y) , u(x, y, 0) = 0.7, v(x, y, 0) = 0.3, p(x, y, 0) = 1.
Four-periodic boundary conditions are used; the final time t = 2 is when no shocks appear in the solution. The exact solution for the problem is u(x, y, t) = 0.7, v(x, y, t) = 0.3, p(x, y, t) = 1.
In Table 2, the numerical errors and accuracy orders for ρ by WENO5-LW3 based on L ∞ and L 1 norm metric are shown, while N is the number of mesh cells. The rectangle cell is used to avoid cancelations of some special error due to the symmetry axis being in the diagonals of cells. L ∞ and L 1 numerical errors and accuracy orders by WENO5-RK3 schemes are also displayed in Table 2 to compare with the results of WENO5-LW3. The two schemes have achieved the fifth order accuracy at the final time and generate similar numerical errors. On mesh 64 × 96, the corresponding CPU time costs for WENO5-LW3 and WENO5-RK3 are 986.73 s and 1103.93 s. We can see that WENO5-RK3 cost about 10% more CPU time than WENO5-LW3.
Example 5 (Two-Dimensional Double Mach Reflection). This problem is a two-dimensional double mach reflection test; complex boundary conditions are implemented. For the example, multiple shocks with complex self-similar structure are emergent, and the initial geometry is very simple [43].
In this test, the wedge of the oblique shock is set as a partially reflecting horizontal lower boundary. The computational domain is [0, 4] × [0, 1], the x axis is 4.0 units long. At the first 1/6 units, the outflow is free and the remaining units reflect the rest of the boundary. At the upper boundary, the flow values are used to describe the exact motion of a shock at Mach = 10. The left boundary is constant with the post oblique shock values, and there is no limit for the right boundary. In the whole calculation, the final time is t = 0.2.
At the ending time, we plot the computed density contours by the WENO5-LW3 scheme with 1920 × 480 meshes in Figure 5, which is almost the same as the results by WENO5-RK3, and comparable with previous results in [43]. When the CPU time is considered, we can find that WENO5-RK3 also costs a little more time than WENO5-LW3.

Concluding Remarks
In this paper, we take advantage of Lax-Wendroff-type time discretizations and high resolution finite volume WENO schemes to simulate one-dimensional and two-dimensional Euler equations. By comparing the addressing on CPU time, non-oscillatory, resolution and relevant efficiency of finite volume WENO schemes using time discetizations, the numerical results show that the WENO5-LW3 schemes takes smaller CPU costs than WENO5-RK3 schemes, while WENO5-LW3 schemes got the same order of accuracy with Runge-Kutta time discretization on the same mesh grids number. These conclusions agree with the theory of the finite difference WENO with Lax-Wendroff-type and Runge-Kutta time discretizations for Euler system of compressible gas dynamics.