A new spectral difference method using hierarchical polynomial bases for hyperbolic conservation laws

Abstract

To solve hyperbolic conservation laws, a new method is developed based on the spectral difference (SD) algorithm. The new scheme adopts hierarchical polynomials to represent the solution in each cell instead of Lagrange interpolation polynomials used by the original one. The degrees of freedom (DOFs) of the present scheme are the coefficients of these polynomials, which do not represent the states at the solution points like the original method. Therefore, the solution points defined in the original SD scheme are discarded, while the flux points are preserved to construct a Lagrange interpolation polynomial to approximate flux function in each cell. To update the DOFs, differential operators are applied to the governing equation as well as the Lagrange interpolation polynomial of flux function to evaluate first and higher order derivatives of both solution and flux at the centroid of the cell. The stability property of the current scheme is proved to be the same as the original SD method when the same solution space is adopted. One dimensional methods are always stable by the use of zeros of Legendre polynomials as inner flux points. For two dimensional problems, the introduction of Raviart–Thomas spaces for the interpolation of flux function proves stable schemes for triangles. Accuracy studies are performed with one- and two-dimensional problems. p-Multigrid algorithm is implemented with orthogonal hierarchical bases. The results verify the high efficiency and low memory requirements of implementation of p-multigrid algorithm with the proposed scheme.

Introduction

In the recent years, some problems of great interest, such as computational aeroacoustics (CAA), direct numerical simulations (DNS) and large eddy simulation (LES) of turbulence, require numerical methods with low dispersion and dissipation. In response to these requests, some high-order methods have been proposed and developed, including essentially nonoscillatory (ENO) and weighted essentially nonoscillatory (WENO) schemes [1], [2], [3], [4], [5], [6], discontinuous Galerkin (DG) method [7], [8], [9], [10], [11], spectral volume (SV) method [16], [17], [18], [19], [20], [21], spectral difference (SD) method [24], [25], [26], [27]. Those works have been summarized comprehensively by Ekaterinaris [29] and Wang [30].

As the pioneer of high-order numerical methods for solution of partial difference equations, the DG method has become popular in CFD community in recent years because of some attractive features, such as local conservative, geometric flexibility, compactness which makes it ideally suited for parallel computing. The DG algorithm was first introduced by Reed and Hill in 1973 to solve neutron transport equation [7]. A major development of the DG scheme to solve hyperbolic conservation laws was carried out by Cockburn, Shu and their collaborators in a series of papers [8], [9], [10], [11] on the Runge–Kutta DG (RKDG) method, which built a framework that combine DG space discretization and Runge–Kutta time discretization to solve nonlinear time dependent hyperbolic conservation laws. However, the original RKDG method suffers a number of weaknesses including the huge memory requirements and high computational cost as order of accuracy is increasing [16]. To avoid these drawbacks, some variations with high efficiency, such as nodal DG method [12] and hybridizable DG method [13], [14], [15], have been proposed.

Wang and his collaborators proposed the SV method [16], [17], [18], [19], [20], [21] which can be viewed as an extension of the finite volume method to higher order by introducing universal local reconstruction concept inherent in the DG scheme. In the SV algorithm, each cell, referred to as spectral volume (SV), is partitioned into structured sub-cells called control volumes (CV). Therefore the DOFs are the CV averaged conservative variables. To update these DOFs, the integral form of governing equation is satisfied at each CV. However, as the order is increasing, the computational cost increases rapidly because of the growth of the number of interior facets and quadrature points for each facet.

To avoid quadrature, the SD fundamentals were introduced by Kopriva and Kolias [22], [23] and later adapted by Liu, Wang and other researchers [24], [25], [26], [27]. The SD method defines two sets of points, i.e., the solution points and flux points in each cell. The conservative variables defined at the solution points provide the DOFs. To evolve these DOFs, the differential form of the governing equation is applied at solution points and the divergence of flux is evaluated in terms of values at flux points. The SD scheme has recently emerged as a promising alternative because it seems easier to implement [32]. In fact, the SD algorithm can also be viewed as a quadrature-free DG method [28]. There are a lot of applications of the SD method which have proved its capabilities [31], [32], [33]. Kris Van den Abeele and his collaborators performed stability analysis of the SD algorithm [34]. They showed that the linear stability property of the SD scheme is independent of the position of the solution points and is strongly affected by the distribution of the flux points. However, whether this feature is also valid for nonlinear problems is not clear. The 1D SD schemes higher than second-order using the Chebyshev–Gauss–Lobatto (C–G–L) nodes as the flux points have a weak instability. For 2D SD schemes higher than second-order on triangular grids, a weak instability was also found. Recently, Jameson [35] proved the stability of the 1D SD method for all orders of accuracy in an energy norm of Sobolev type. He found that the use of zeros of the Legendre polynomial as inner flux points could always produce stable 1D SD scheme. May and Schöberl [36] constructed stable high-order SD algorithm on triangular element by introducing a new flux interpolation technique using Raviart–Thomas (RT) spaces.

All high-order methods mentioned above evolve a polynomial in each cell along time instead of cell-averaged variables as in the finite volume algorithm, so they are compact. Hierarchical polynomials are attractive for the convenience of handling p- or hp-multigrid strategy and p-adaptive algorithm [38]. To meet the requirement of these algorithms, Luo and some collaborators provided a DG framework using Taylor basis functions as bases of solution [38].

As an acceleration technique, the p-multigrid method takes advantages of the polynomials of the high-order schemes without generating a series of geometric meshes. This technique has been implemented in conjunction with the DG method [40], [41] and the SD method [42], [43], [44]. The motivation of the present paper is to provide a new scheme to implement the hierarchical bases in the framework of the SD method to facilitate the implementation of the p-multigrid algorithm.

In this article, a new conservative high-order numerical discretization technique based on the SD scheme is developed for conservation laws. This algorithm evolves hierarchical polynomials in each cell instead of Lagrange interpolation polynomials defined at solution points like the original SD method. Therefore the solution points are discarded. The flux points are preserved to construct a Lagrange polynomial for evaluation of the first and higher order derivatives of flux. Then, the governing equation and its derivatives are solved to evolve the DOFs.

The remainder of this article is organized as follows. Section 2 gives a brief review of the SD method. Section 3 describes details of implementation of hierarchical polynomials in the SD method. A discussion about linear stability of the new algorithm is also given in this section. To illustrate the advantage of the current scheme, p-multigrid technique is implemented. The details of the p-multigrid algorithm are given in Section 4. To reveal the order of accuracy of the new scheme and the efficiency of the p-multigrid method, some numerical cases are presented in Section 5. Concluding remarks and some possible future works are included in the last section.