Weighted essentially non-oscillatory scheme on unstructured quadrilateral and triangular meshes for hyperbolic conservation laws

Weighted essentially non-oscillatory scheme on unstructured quadrilateral and triangular meshes for hyperbolic conservation laws

Fengxiang Zhao kobezhao@126.com Liang Pan panliangjlu@sina.com Shuanghu Wang wang_shuanghu@iapcm.ac.cn The Graduate School of China Academy of Engineering Physics, Beijing, China Institute of Applied Physics and Computational Mathematics, Beijing, China
Abstract

In this paper, a third-order weighted essentially non-oscillatory (WENO) scheme is developed for hyperbolic conservation laws on unstructured quadrilateral and triangular meshes. As a starting point, a general stencil is selected for the cell with any local topology, and a unified linear scheme can be constructed. However, in the traditional WENO scheme on unstructured meshes, the very large and negative weights may appear for the mesh with lower quality, and the very large weights make the WENO scheme unstable even for the smooth tests. In the current scheme, an optimization approach is given to deal with the very large linear weights, and the splitting technique is considered to deal with the negative weights obtained by the optimization approach. The non-linear weight with a new smooth indicator is proposed as well, in which the local mesh quality and discontinuities of solutions are taken into account simultaneously. Numerical tests are presented to validate the current scheme. The expected convergence rate of accuracy is obtained, and the absolute value of error is not affected by mesh quality. The numerical tests with strong discontinuities validate the robustness of current WENO scheme.

keywords:
Third-order WENO scheme, Unstructured meshes, Hyperbolic conservation laws, Finite volume method.

1 Introduction

In recent decades, there have been continuous interests and efforts on the development of high-order schemes. For engineering applications, the construction of high-order numerical schemes on unstructured meshes becomes extremely demanding due to the complex computational domains. There are a gigantic amount of publications about the introduction and survey of high-order schemes, including discontinuous Galerkin (DG), spectral volume (SV), spectral difference (SD), correction procedure using reconstruction (CPR), essential non-oscillatory (ENO), weighted essential non-oscillatory (WENO), Hermite WENO (HWENO) and -exact method, etc.

The DG scheme was first proposed DG1 () to solve the neutron transport equation. The major development of DG method was carried out for the compressible Euler equations DG2 (); DG3 (); DG5 (). High-order accuracy is achieved by means of high-order polynomial approximation within each element rather than by means of wide stencils. Because only neighboring elements interaction is included, it becomes efficient in the application of complex geometry. Recently, a correction procedure via reconstruction framework (CPR) was developed CPR1 (); CPR2 (); CPR3 (); CPR4 (). The CPR formulation is based on a nodal differential form, with an element-wise continuous polynomial solution space. By choosing certain correction functions, the CPR framework can unify several well known methods, such as the DG, SV SV (), and SD SD () methods, and lead to simplified versions of these methods, at least for linear equations. The ENO scheme was proposed in Harten (); ENO-Shu () and successfully applied to solve the hyperbolic conservation laws and other convection dominated problems. The “smoothest” stencil is selected among several candidates in reconstruction to achieve high order accuracy and keep essentially non-oscillatory near discontinuities. Based on unstructured triangular meshes, the ENO scheme was developed as well Harten-unstructured (); Abgrall (); Sonar (). Following the ENO scheme, the WENO scheme WENO-Liu (); WENO-JS (); WENO-M (); WENO-Z (); Zhao () was developed. With the nonlinear convex combination of candidate polynomials, WENO scheme achieves higher order of accuracy and keeps essentially non-oscillatory property. WENO scheme improves upon ENO scheme in robustness, smoothness of fluxes, steady-state convergence, provable convergence properties, and more efficiency. The classical development of WENO scheme includes WENO-JS scheme WENO-JS (), WENO-M scheme WENO-M (), WENO-Z scheme WENO-Z (), et al. Based on the Hermite interpolation, HWENO schemes were developed HWENO1 (); HWENO2 (); HWENO3 () for the hyperbolic conservation laws. On the unstructured meshes, the WENO scheme with the same order of accuracy as the fundamental ENO scheme was constructed by Friedrichs Friedrich (). Compared with the fundamental ENO scheme, the WENO scheme is very smooth and stable in smooth regions, while the optimal order of accuracy is not obtained WENO-JS (). Hu and Shu Hu-Shu () presented the third-order and fourth-order WENO schemes on unstructured triangular meshes. The high-order of accuracy is obtained by the combination of lower order polynomials, which is similar with the one-dimensional WENO scheme. However, its successful application was limited by the appearance of the negative linear weights, and the problem appears commonly on the unstructured meshes. For a mesh that is close to the regular meshes, such WENO scheme works well by a regrouping approach to avoid negative weights. However, for a mesh with lower quality, the large linear weights appears and WENO schemes become unstable even for the smooth flows. A splitting technique dealing with the negative linear weights was proposed in splitting-weights (). This technique makes the WENO schemes Hu-Shu () generally work well on better quality meshes. However, the linear weights can be very large with the lower quality mesh, and the WENO schemes become unstable even if the splitting technique is adopted. In order to avoid the negative linear weights and very large linear weights, a class of WENO schemes were proposed in Dumbser (). Due to the very large stencils, the boundary treatment becomes complex and consuming computer memory, and the computational efficiency is reduced because of the inefficient combination in reconstruction. Based on the two classes of WENO schemes in Hu-Shu (); Dumbser (), a hybrid scheme is presented in hybrid-WENO () to deal with the problem with very large linear weights, while the scheme contains artificial parameter and the reconstruction is more complex. Another type of higher-order finite volume scheme is based on the so-called -exact reconstruction k-exact1 (); k-exact2 (); k-exact3 (). To guarantee the non-singularity of the reconstruction, the -exact reconstruction relies on the least-square approach.

In this paper, a third-order WENO scheme is developed on unstructured quadrilateral and triangular meshes for the hyperbolic conservation laws. As a starting point of WENO reconstruction, a general stencil is selected for a cell with any local topology. With the selected stencil, a unified linear scheme was constructed. With the linear candidate polynomials from the candidate sub-stencils, the linear wights can be obtained. However, the appearance of very large linear wights on unstructured meshes causes the WENO schemes to be unstable Hu-Shu (); hybrid-WENO (). For both quadrilateral and triangular meshes, an optimization approach is introduced to deal with the very large linear weights. The splitting technique is considered to deal with the negative weights obtained by the optimization approach. The non-linear weight with a new smooth indicator is defined for the solutions with discontinuities. With the optimization approach for very large weights and the splitting technique for negative weights, the scheme becomes more robust compared with the previous WENO scheme Hu-Shu () with the low quality meshes. A variety of numerical tests from the accuracy test to the solutions with strong discontinuities are presented to validate the accuracy and robustness of the current scheme.

This paper is organized as follows. In Section 2, the finite volume scheme on unstructured meshes are introduced. The third-order WENO schemes on unstructured quadrilateral and triangular meshes are presented in section 3. Section 4 includes numerical tests to validate the current algorithm. The section 5 is the conclusion.

2 Finite volume method

In this paper, the two-dimensional Euler equations are considered

(1)

where , , , and is the specific heat ratio. For a polygon cell , the boundary can be expressed as

where is the number of cell interface for cell . Integrating Eq.(1) over the cell , the semi-discretized form of finite volume scheme can be written as

(2)

where is the cell averaged value over cell , is the area of , and is the outer normal direction of .

In this paper, a third-order spatial reconstruction will be introduced, and the line integral over is discretized according to Gaussian quadrature as follows

(3)

where is the length of the cell interface , are the Gaussian quadrature weights, and the Gaussian quadrature points for are defined as

where are the endpoints of .

For the two-dimensional Euler equations, the rotational invariance property Riemann-Solver () can be expressed as

(4)

where is the rotation matrix

The numerical fluxes in Eq.(3) at the Gaussian quadrature points in the global coordinate can be given with the fluxes in the local coordinate according to Eq.(4), and the flux in the local coordinate can be approximated by Riemann solvers

where and are the reconstructed variables at both sides of Gaussian quadrature points of cell interface. In this paper, the HLLC Riemann solver Riemann-Solver () is adopted for numerical tests.

3 Third-order WENO scheme on unstructured meshes

The reconstruction procedure is the main theme of this paper, and it will be presented on both unstructured quadrilateral and triangular meshes in this section.

3.1 Linear reconstruction on quadrilateral meshes

As a starting point of WENO reconstruction, a linear reconstruction will be presented. For a piecewise smooth function over cell , a polynomial with degree can be constructed to approximate as follows

In order to achieve the third-order accuracy and satisfy conservative property, the following quadratic polynomial on the cell is constructed

(5)

where is the cell average value of over cell and are basis functions, which are given as follows

(6)
Figure 1: Stencil for third-order scheme on quadrilateral meshes.

In order to construct the polynomial uniquely, the stencil for reconstruction is given in Figure.1. For arbitrary unstructured quadrilateral meshes, the cells must exist, and the cells may not exist since it is common that only three or four cells share a common node. Thus, preconditioning is needed for a general quadrilateral mesh, and the procedure of selecting stencil for cell is given as follows

  1. initialize the parameter , which is corresponding to the cell ;

  2. read the neighboring cell of cell from the mesh file;

  3. read the neighboring cell of cell , if each of them is not exist, record and is replaced by ,
    read the neighboring cell of cell , if each of them is not exist, record and is replaced by ;

  4. read the neighboring cell of cell , if each of them is not exist, record and is replaced by ,
    read the neighboring cell of cell , if each of them is not exist, record and is replaced by .

To determine the polynomial , the following constrains need to be satisfied for the cells

(7)

where is the cell averaged value over . According to the definition of , Eq.(7) can be written into the following linear system

where , , and is given as

Traditionally, the following least square method is used to obtain the coefficient

(8)

However, the matrix may not be invertible. In order to obtain the coefficients in the quadratic polynomial uniquely, a more robust least square method is adopted

(9)

where I is the unit matrix, is a small number, and takes in the computation. The matrix is invertible, and the quadratic polynomial can be determined uniquely. Meanwhile, if the matrix is invertible, the difference between the solution of Eq.(9) and Eq.(8) is . With the coefficient , the point value of at the Gaussian quadrature point can be written as a linear combination of the cell averaged values

(10)

where if .

Remark 1

Denote , and the inverse of can be expressed as

(11)

Multiplying on both side of Eq.(11), we have

According to the Hamilton-Caylay theorem, can be expressed as

(12)

where are given by

Substituting Eq.(12) into Eq.(11), the linear system for can be written as

If the matrix is invertible, and the solution of the linear system above can be given as

where corresponds to the solution of the system with . Substituting c into Eq.(11), and supposing the solution of Eq.(8) is , the solution of Eq.(9) can be given as

Thus, the difference between solutions of Eq.(9) and Eq.(8) is the order of .

Similar with the standard WENO reconstruction WENO-Liu (); Hu-Shu (), twelve sub-stencils are selected from the large stencil given in Figure.1

Twelve candidate linear polynomials corresponding to the candidate sub-stencils can be constructed as follows

where is the -th cell in the sub-stencils , and is the cell averaged value over cell . With the robust least square method introduced above, twelve linear polynomials can be fully determined. Similarly, the point value of at the Gaussian quadrature point can be written as a linear combination of the cell averaged values . However, in order to calculate the linear weights, it must be taken into account that the cell may not exist for arbitrary unstructured quadrilateral meshes. According to the local topology, the candidate polynomials is rewritten as follows

  1. For , the linear polynomial can be expressed as

  2. For , the linear polynomial can be expressed as

For the linear scheme, the linear combination of

satisfies

(13)

where is the Gaussian quadrature point, and is the linear weights only depending on the local geometry of mesh. Substituting into at , we have

(14)

Comparing the coefficient of the cell average value in Eq.(10) and Eq.(14), the following thirteen linear equations can be obtained

and its matrix form is

(15)

where , , and . Similar with the analysis in Hu-Shu (), the linear system for Eq.(15) contains thirteen linear equations, and generally the system is under-determined with . and reproduce the linear function exactly and Eq.(13) is valid for , and under the following identical constraint

Two equations can be eliminated from the linear system with thirteen equations. The traditional least square becomes under-determined, and the following least square method are used to obtain the linear weights

3.2 Linear reconstruction on triangular meshes

For the reconstruction on unstructured triangular meshes, the procedure is similar with the reconstruction on quadrilateral meshes, except the choice for stencil and sub-stencils. The stencil is given in Figure.2, and the quadratic polynomial Eq.(5) can be constructed on cell as well. For the cell , the neighboring cell must exist and their neighboring cells may not exist, and the similar preconditioning is needed as well, and the procedure is given as follows

  1. initialize the parameter , which is corresponding to the cell ;

  2. read the neighboring cell of cell from the mesh file;

  3. read the neighboring cell of cell , if each of them is not exist, record and is replaced by ;

  4. read the neighboring cell of cell , if each of them is not exist, record and is replaced by .

Figure 2: Stencil for third-order scheme on triangular meshes.

In order to obtain the quadratic polynomial, the following constrains need to be satisfied

and the matrix form can be written as

where , , and is given as

Similar with the quadrilateral meshes, the polynomial can be determined. For the triangular meshes, nine candidate sub-stencils are selected

and nine linear polynomials can be fully determined. Similarly, the point value of at the Gaussian quadrature point can be written as a linear combination of the cell averaged values . However, the cell may not exist for arbitrary unstructured meshes. In order to calculate linear weights later, the candidate polynomials is given according to the local topology as follows

  1. For , the linear polynomial can be expressed as

  2. For , the linear polynomial can be expressed as

For the linear scheme, the linear combination

satisfies

(16)

where is the Gaussian quadrature point, and is the linear weights. The procedure of determining the linear weights for quadrilateral meshes can be extended to triangular meshes directly.

3.3 Optimization approach for linear weights

With the procedure introduced above, the linear weights for WENO reconstruction can be obtained. To deal with the flow with discontinuities, the non-linear weights need to be constructed, which is one of the most important step in the WENO schemes. However, for the unstructured meshes, the very large linear weights and negative linear weights appear commonly, which makes the WENO schemes unstable even for the smooth flow problems. In order to construct a robust WENO schemes on unstructured meshes, the very large linear weights and negative linear weights need to be treated carefully. To illustrate the very large weights, the quadrilateral and triangular meshes in Figure.4 and Figure.4 are considered, where the meshes with cell size are given. The meshes in the left column are closer to the uniform quadrilateral and triangular meshes, and they denoted as the regular meshes. Meanwhile, the meshes in the right column are denoted as irregular meshes. The maximum of the liner weights for these meshes are given in Table.2, where the maximum becomes larger with the mesh refinement. Such very large weights will trigger the instability of WENO scheme.

Figure 3: The quadrilateral meshes with cell size : regular mesh (left) and irregular mesh (right).
Figure 4: The triangular meshes for with cell size : regular mesh (left) and irregular mesh (right).
Figure 3: The quadrilateral meshes with cell size : regular mesh (left) and irregular mesh (right).

  mesh size regular quad irregular quad regular tri irregular tri 4.3589E01 4.9770E02 2.1023E02 2.6012E02 2.9415E02 8.4584E02 1.0262E02 2.9049E05 1.9161E02 2.7036E03 9.6171E02 2.6894E03 6.6936E02 2.7012E04 1.9599E03 1.1906E04  

Table 1: Maximum of the liner weights for the quadrilateral and triangular meshes.

  mesh size regular quad irregular quad regular tri irregular tri 1.3897E00 1.4613E00 1.6410E00 2.0553E00 1.3719E00 1.5690E00 1.7563E00 2.0730E00 1.4029E00 1.5359E00 1.6214E00 8.5631E00 1.3910E00 1.9566E00 1.6269E00 6.4324E00  

Table 2: Maximum of the optimal liner weights for the quadrilateral and triangular meshes.

The meshes with lower quality will result in the very large linear weights, and an optimization approach is introduced to deal with them. In the current optimization approach, the weighting parameters is introduced for each cell, such that the ill cell contributes little to the quadratic polynomial in Eq.(10). For quadrilateral meshes, the weighting parameters are defined

For triangular meshes, the weighting parameters are defined

The optimized coefficients of the quadratic polynomial in Eq.(10) are given by the following weighted linear system

where for the quadrilateral meshes, and for the triangular meshes. The maximum of optimized liner weights for the meshes in Figure.4 and Figure.4 are given in Table.2. With the mesh refinement, the maximum of liner weights is reduced greatly. This linear system is solved with the least square method proposed above. With the procedure above, the maximum becomes the order .

3.4 Nonlinear weights

Before implementing the WENO reconstruction, the splitting technique splitting-weights () is considered for the negative linear weights which is obtained by the optimized approach as follows

where is taken in numerical tests, and the scaled non-negative weights are given by

With the non-negative weights, the splitting polynomial is written as

which satisfies

Based on the non-negative weights obtained by splitting approach, the nonlinear weights are obtained as

where is a small positive number, is a new smooth indicator defined on unstructured meshes

and is defined as

where is a multi-index and is the derivative operator. In the new smooth indicator, the accuracy keeps the original order in smooth regions with . In addition, the linear weights are considered in constructing the new smooth indicator. Suppose that is a very large linear weight and the corresponding contains the discontinuity. The ratio of nonlinear weights of and becomes

The contribution of the candidate polynomial becomes very little, and the essentially non-oscillatory property retains for the very large linear weight. The final reconstructed value at the Gaussian quadrature points can be written as

4 Numerical tests

In this section, the numerical tests will be presented to validate the current WENO scheme. For the temporal discretization, the classical third-order TVD Runge-Kutta method TVD-RK () is used for Eq.(2), which can be written as follows

where is the operator corresponding to the spatial discretization.

In order to eliminate the spurious oscillation and improve the stability, the reconstruction for Euler system can be performed on the characteristic variables. The Jacobian matrix at one quadrature point is , and the Roe’s mean matrix is used, which can be diagnoalized by the right eigenmatrix . The variables for reconstruction are defined as . With the reconstructed values, the conservative variables can be obtained by the inverse projection.

  scheme mesh error Order error Order WENO-3 2.533E-01 1.058E-01 regular 3.889E-02 2.70 2.348E-02 2.17 meshes 3.730E-03 3.38 2.399E-03 3.29 2.127E-04 4.13 1.480E-04 4.02 WENO-3 2.220E-01 9.791E-02 irregular 3.262E-02 2.77 2.122E-02 2.21 meshes 2.602E-03 3.65 1.877E-03 3.50 1.544E-04 4.07 1.122E-04 4.06  

Table 3: Accuracy tests: linear wave propagation on quadrilateral with WENO scheme.

  Schemes mesh error Order error Order WENO-3 6.140E-02 3.372E-02 regular 6.901E-03 3.15 4.591E-03 2.88 meshes 4.822E-04 3.84 3.410E-04 3.75 5.633E-05 3.10 2.770E-05 3.62 WENO-3 5.419E-02 3.211E-02 irregular 6.226E-03 3.12 4.505E-03 2.83 meshes 4.227E-04 3.88 3.288E-04 3.78 4.845E-05 3.13 2.550E-05 3.69  

Table 4: Accuracy tests: linear wave propagation on triangular with WENO scheme.
Figure 5: Advection of density perturbation: density distribution on the quadrilateral meshes with cell size , regular mesh (left) and irregular mesh (right).
Figure 6: Advection of density perturbation: density distribution on the triangular meshes with cell size , regular mesh (left) and irregular mesh (right).
Figure 5: Advection of density perturbation: density distribution on the quadrilateral meshes with cell size , regular mesh (left) and irregular mesh (right).

  scheme mesh error Order error Order WENO-3 2.291E-01 4.343E-02 regular 3.013E-02 2.93 6.597E-03 2.72 meshes 2.984E-03 3.34 6.715E-04 3.30 3.431E-04 3.12 6.543E-05 3.36 WENO-3 1.637E-01 3.605E-02 irregular 2.094E-02 2.97 5.355E-03 2.75 meshes 2.346E-03 3.16 6.094E-04 3.14 2.794E-04 3.07 6.976E-05 3.13  

Table 5: Accuracy tests: isentropic vortex propagation on quadrilateral with WENO scheme.

  Schemes mesh error Order error Order WENO-3 5.515E-02 1.236E-02 regular 6.256E-03 3.14 1.445E-03 3.10 meshes 7.219E-04 3.12 1.368E-04 3.40 1.328E-04 2.44 1.483E-05 3.21 WENO-3 5.248E-02 2.870E-02 irregular 5.200E-03 3.34 1.370E-03 4.39 meshes 6.758E-04 2.94 1.422E-04 3.27 1.270E-04 2.41 1.997E-05 2.83  

Table 6: Accuracy tests: isentropic vortex propagation on triangular with WENO scheme.

4.1 Accuracy tests

The first case is the advection of density perturbation, and the initial condition is given as follows

The computational domain is , and the periodic boundary conditions are applied in both directions. To keep the third-order accuracy, is used for the third-order scheme where is the scale of meshes. In order to validate the accuracy of the current scheme, both quadrilateral and triangular meshes are tested. The meshes with mesh size in Figure.4 and Figure.4 are shown as example, and the density distribution with cell size are given in Figure.6 and Figure.6. The and errors and convergence orders are presented in Table.4 and Table.4 for quadrilateral and triangular meshes at . The expected order of accuracy are obtained, and the order of error is not affected by the quality of meshes.

The second accuracy test is the isentropic vortex propagation problem. The mean flow is . An isentropic vortex is added to the mean flow with the following perturbations

where , , and the vortex strength . The computational domain is , and periodic boundary condition is applied to all boundaries. The regular and irregular quadrilateral meshes, regular and irregular triangular meshes are tested respectively. The and errors and convergence orders for current scheme are presented in Table.6 and Table.6 for quadrilateral and triangular meshes at . The expected order of accuracy are obtained, and the order of error is not affected by the quality of meshes as well. With the fixed mesh size , the long-time evolution of the vortex is computed. The density distribution at center line with and time periods on quadrilateral and triangular irregular meshes are given in Figure.9 with mesh size . The dissipation is well controlled by the current scheme even after long time evolution. Due to more number of cell with the same mesh scale, the scheme performs better on triangular meshes than quadrilateral meshes.

Figure 7: Isentropic vortex propagation: density distribution at center line with and time periods on quadrilateral and triangular irregular meshes with mesh size .
Figure 8: One dimensional Riemann problem: quadrilateral regular (left) and quadrilateral irregular (right) meshes.
Figure 9: One dimensional Riemann problem: triangular regular (left) and triangular irregular (right) meshes with.
Figure 8: One dimensional Riemann problem: quadrilateral regular (left) and quadrilateral irregular (right) meshes.
Figure 10: One dimensional Riemann problem: the density, velocity and pressure distribution for Sod problem at (left) and Lax problem at (right) on quadrilateral and triangular meshes.

4.2 One dimensional Riemann problem

In this case, four one-dimensional Riemann problems are tested by the third-order WENO scheme on the quadrilateral and triangular meshes. The first one is the Sod problem, and the initial condition is given as follows

The second one is the Lax problem, and the initial condition is given as follows

For these two cases, the computati