Combined prefactored compact schemes for first- and second-order derivatives: conceptual derivation

TL;DR

This paper presents a conceptual derivation of combined prefactored compact schemes that enhance the accuracy of first and second derivatives from sixth to tenth order without increasing stencil size, using Fourier analysis.

Contribution

It introduces a novel derivation method for higher-order compact schemes that maintains computational efficiency while improving accuracy.

Findings

Order of accuracy increased from sixth to tenth

Same stencil size used for higher accuracy

Closed set of equations for scheme weights derived

Abstract

The derivation of combined prefactored compact schemes for first and second order derivatives is described here, relying on the Fourier analysis of the original prefactored compact schemes. By this approach, the order of accuracy of the original schemes can be increased from sixth to eight, or from eight to tenth (depending on the order of the original scheme), while the number of grid points in the stencil is kept the same. Here, we only frame the conceptual derivation of the schemes, leading to a closed set of equations for the weights.

Equations

$$u_i^{\prime }+\frac{7}{16}\left(u_{i+1}^{\prime }+u_{i-1}^{\prime }\right)-\frac{h}{16}\left(u_{i+1}^{\prime \prime }-u_{i-1}^{\prime \prime }\right)=\frac{15}{16h}\left(u_{i+1}-u_{i-1}\right)$$

$$u_i^{\prime \prime }+\frac{9}{8h}\left(u_{i+1}^{\prime }-u_{i-1}^{\prime }\right)-\frac{1}{8}\left(u_{i+1}^{\prime \prime }+u_{i-1}^{\prime \prime }\right)=\frac{3}{h^2}\left(u_{i+1}+u_{i-1}\right)-\frac{6}{h^2}{u}_i$$

$$D_i+\frac{7}{16}\left(D_{i+1}+D_{i-1}\right)-\frac{h}{16}\left(D_{i+1}^2-D_{i-1}^2\right)=\frac{15}{16h}\left(u_{i+1}-u_{i-1}\right)$$

$$D_i^2+\frac{9}{8h}\left(D_{i+1}-D_{i-1}\right)-\frac{1}{8}\left(D_{i+1}^2+D_{i-1}^2\right)=\frac{3}{h^2}\left(u_{i+1}+u_{i-1}\right)-\frac{6}{h^2}{u}_i$$

$$D_i+\frac{17}{36}\left(D_{i+1}+D_{i-1}\right)-\frac{h}{12}\left(D_{i+1}^2-D_{i-1}^2\right)=\frac{107}{108h}\left(u_{i+1}-u_{i-1}\right)-\frac{1}{108h}\left(u_{i+2}-u_{i-2}\right)$$

$$D_i^2+\frac{23}{18h}\left(D_{i+1}-D_{i-1}\right)-\frac{1}{6}\left(D_{i+1}^2+D_{i-1}^2\right)=\frac{88}{27h^2}\left(u_{i+1}+u_{i-1}\right)-\frac{1}{108h^2}\left(u_{i+2}-u_{i-2}\right)-\frac{13}{2h^2}{u}_i$$

$$[B]\{D\}=[C]\{u\}$$

$$D_i=\frac{1}{2}(D_i^F+D_i^B)$$

$$D_i^2=\frac{1}{2}(D_i^{2F}+D_i^{2B})$$

$$D_i^F+{\beta }_F^I{D}_{i+1}^F+{\theta }_F^{I}h{D}_{i+1}^{2F}=\frac{1}{h}\left(a_F^I{u}_{i-1}+b_F^I{u}_i+c_F^I{u}_{i+1}\right)$$

$$D_i^B+{\beta }_B^I{D}_{i-1}^B+{\theta }_B^{I}h{D}_{i-1}^{2B}=\frac{1}{h}\left(a_B^I{u}_{i-1}+b_B^I{u}_i+c_B^I{u}_{i+1}\right)$$

$$D_i^{2F}+{\beta }_F^{II}\frac{1}{h}{D}_{i+1}^F+{\theta }_F^{II}{D}_{i+1}^{2F}=\frac{1}{h^2}\left(a_F^{II}{u}_{i-1}+b_F^{II}{u}_i+c_F^{II}{u}_{i+1}\right)$$

$$D_i^{2B}+{\beta }_B^{II}\frac{1}{h}{D}_{i-1}^B+{\theta }_B^{II}{D}_{i-1}^{2B}=\frac{1}{h^2}\left(a_B^{II}{u}_{i-1}+b_B^{II}{u}_i+c_B^{II}{u}_{i+1}\right)$$

$$u(x)=\sum _{k=-N/2}^{k=N/2}{\hat{u}}_k{e}^{\frac{2\pi ikx}{L}}$$

$${\hat{u}}_k^{\prime }=\frac{iw}{h}{\hat{u}}_k,{\hat{u}}_k^{\prime \prime }=-{\left(\frac{w}{h}\right)}^2{\hat{u}}_k$$

$$({\hat{u}}_k^{\prime }{)}_{CCD}=\frac{i{w}^{\prime }}{h}{\hat{u}}_k,({\hat{u}}_k^{\prime \prime }{)}_{CCD}=-{\left(\frac{w^{\prime \prime }}{h}\right)}^2{\hat{u}}_k$$

$$w^{\prime }(w)=\frac{9\sin w[4+\cos w]}{24+20\cos w+\cos 2w}$$

$$w^{\prime \prime 2}(w)=\frac{81-48\cos w-33\cos 2w}{48+40\cos w+2\cos 2w}$$

$$w^{\prime }(w)=\frac{\sin w[293+126\cos w+\cos 2w]}{6(34+33\cos w+3\cos 2w)}$$

$$w^{\prime \prime 2}(w)=\frac{1730-675\cos w-10870\cos 2w-29\cos 3w}{36(34+33\cos w+3\cos 2w)}$$

$$i{w}_F^{\prime }+i{\beta }_F^I{w}_F^{\prime }{e}^{iw}-{\theta }_F^I(w_F^{\prime \prime }{)}^2{e}^{iw}=a_F^I{e}^{-iw}+b_F^I+c_F^I{e}^{iw}$$

$$-(w_F^{\prime \prime }{)}^2+i{\beta }_F^{II}{w}_F^{\prime }{e}^{iw}-{\theta }_F^{II}(w_F^{\prime \prime }{)}^2{e}^{iw}=a_F^{II}{e}^{-iw}+b_F^{II}+c_F^{II}{e}^{iw}$$

$$i{w}_B^{\prime }+i{\beta }_B^I{w}_B^{\prime }{e}^{-iw}-{\theta }_B^I(w_B^{\prime \prime }{)}^2{e}^{-iw}=a_B^I{e}^{-iw}+b_B^I+c_B^I{e}^{iw}$$

$$-(w_B^{\prime \prime }{)}^2+i{\beta }_B^{II}{w}_B^{\prime }{e}^{-iw}-{\theta }_B^{II}(w_B^{\prime \prime }{)}^2{e}^{-iw}=a_B^{II}{e}^{-iw}+b_B^{II}+c_B^{II}{e}^{iw}$$

$$[A]\{X\}=\{R\}$$

$$\displaystyle[A]=\left[\begin{array}[]{ccccc}-\beta_{F}^{I}\sin{w}&-1-\beta_{F}^{I}\cos{w}&-\theta_{F}^{I}\cos{w}&\theta_{F}^{I}\sin{w}\\ 1+\beta_{F}^{I}\cos{w}&-\beta_{F}^{I}\sin{w}&-\theta_{F}^{I}\sin{w}&-\theta_{F}^{I}\cos{w}\\ -\beta_{F}^{II}\sin{w}&-\beta_{F}^{II}\cos{w}&-1-\theta_{F}^{II}\cos{w}&\theta_{F}^{II}\sin{w}\\ \beta_{F}^{II}\cos{w}&-\beta_{F}^{II}\sin{w}&-\theta_{F}^{II}\sin{w}&-1-\theta_{F}^{II}\cos{w}\end{array}\right],$$

$$\displaystyle\{X\}=\{\begin{array}[]{c}\Re{(w_{F}^{\prime})}\hskip 11.38109pt\Im{(w_{F}^{\prime})}\hskip 11.38109pt\Re{((w_{F}^{\prime\prime})^{2})}\hskip 11.38109pt\Im{((w_{F}^{\prime\prime})^{2})}\end{array}\}^{T},$$

$$\displaystyle\{R\}=\{\begin{array}[]{c}(c_{F}^{I}+a_{F}^{I})\cos{w}+b_{F}^{I}\hskip 11.38109pt(c_{F}^{I}-a_{F}^{I})\sin{w}\hskip 11.38109pt(c_{F}^{II}+a_{F}^{II})\cos{w}+b_{F}^{II}\hskip 11.38109pt(c_{F}^{II}-a_{F}^{II})\sin{w}\end{array}\}^{T}$$

$$\Re (w_F^{\prime })=\frac{\sin w[f_1^I+f_2^I\cos w+f_3^I\cos 2w]}{g_1+g_2\cos w+g_3\cos 2w}$$

$$\Re ((w_F^{\prime \prime }{)}^2)=\frac{f_1^{II}+f_2^{II}\cos w+f_3^{II}\cos 2w+f_4^{II}\cos 3w}{g_1+g_2\cos w+g_3\cos 2w}$$

Full text

1 Introduction

Compact difference schemes possess the advantage of attaining higher-order of accuracy with fewer grid points per stencil. They are preferred in applications where high accurate results are desired, such as direct numerical simulations, large eddy simulations, computational aeroacoustics or electromagnetism, to enumerate few. One of the disadvantages of compact schemes is that an implicit approach is required to determine the grid functions, wherein a large matrix has to be inverted. A comprehensive study of high-order compact schemes approximating both first and second derivatives on a uniform grid was performed by Lele [11]. A wavenumber based optimization was introduced wherein the dispersion error was reduced significantly, achieving spectral-like resolution. Over the next years, compact schemes have been studied by many research groups, and applied to various engineering problems (see for example, Li et al. [12], Adams and Shariff [1], Liu [13], Deng and Maekawa [4], Fu and Ma [6, 7], Meitz and Fasel [16], Shen et al. [19], Shah et al. [18]). Other worth-mentioning examples include Kim and Lee [10] who performed an analytic optimization of compact finite difference schemes, Mahesh [15] who derived a family of compact finite difference schemes for the spatial derivatives in the Navier-Stokes equations based on Hermite interpolations (see also, Chu and Fan [3] for a similar prior analysis), or Deng and Zhang [5] who developed compact high-order nonlinear schemes which are equivalent to fifth-order upwind biased explicit schemes in smooth regions.

Hixon [8, 9] derived prefactored high-order compact schemes that use three-point stencils and returns up to eighth-order of accuracy. These schemes combine the tridiagonal compact formulation with the optimized split derivative operators of an explicit MacCormack type scheme. The optimization of Hixon’s [8, 9] schemes in terms of reducing the dispersion error was performed by Ashcroft and Zhang [2] who used Fourier analysis to select the coefficients of the biased operators such that the dispersion characteristics match those of the original centered compact scheme and their numerical wavenumbers have equal and opposite imaginary components. Today, compact schemes are widely used in numerical simulations of turbulent flows (e.g., direct numerical simulations), computational aeroacoustics, or computational electromagnetics. In order to increase the speed of such numerical simulations it is desirable to derive more computational efficient compact schemes without affective the order of accuracy and the wavenumber characteristics.

In this work, we frame the conceptual derivation of combined prefactored compact schemes for first and second order derivatives, aimed at increasing the resolution accuracy with fewer points per stencil. They are based on the type of prefactorization introduced previously by Hixon [8, 9]. One of the advantages of these schemes is that the derivatives are explicitly determined by sweeping from one boundary to the other, thus avoiding the inversion of matrices which can increase the computational time significantly.

2 Derivation of the combined prefactored schemes

Chu and Fan [3] derived a combined compact difference scheme (over a three-point stencil), which is sixth-order accurate (see also Mahesh [15] for a generalization). The scheme consists of two coupled equations for the first and second derivatives and can be written as:

[TABLE]

[TABLE]

or, following the notation used in [8],

[TABLE]

[TABLE]

where $h$ is the grid step, $u$ is the grid function, and $D$ and $D^{2}$ stand for first and second order derivatives. Mahesh [15] extended the idea to higher-order, such as the eight-order coupled stencils for the first and second derivatives in the form:

[TABLE]

[TABLE]

Equations (3) and (4) or (5) and (6) can be jointly written in matrix form:

[TABLE]

where $[B]$ and $[C]$ are $2N\times 2N$ matrices, $\{D\}$ is a vector of $2N$ elements containing the first and the second derivatives for $i=1,2,...,N$ and $\{u\}$ is a vector of $2N$ elements containing the dependent variables twice ($N$ is the number of grid points along a grid line).

To derive the prefactored compact schemes, the forward and backward operators $D_{i}^{F}$ and $D_{i}^{B}$ for the first derivative and $D_{i}^{2F}$ and $D_{i}^{2B}$ for the second derivative are defined as

[TABLE]

and

[TABLE]

The spatial stencils for the forward and backward derivative operators are defined as a combination of both the first and second derivatives:

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

where the coefficients $a$, $b$, and $c$ must be chosen such that when the two biased stencils are added, the original stencils of Chu and Fan [3] or Mahesh [15] are recovered. Following the analysis of Hixon and Turkel [8], the real components of the numerical wavenumbers of the forward and backward stencils are required to be equal and identical to the numerical wavenumber of the original scheme, and the imaginary components of the numerical wavenumbers are required to be equal and opposite.

To apply the Fourier analysis, the dependent variable $u$ is assumed to be periodic in the domain $[0,L]$. The Fourier decomposition of $u$ is:

[TABLE]

where $i=\sqrt{-1}$. To simplify the analysis, a scaled wavenumber $w=2\pi kh/L$ and a scaled coordinate $s=x/h$ are introduced, such that the Fourier modes are simply $exp(iws)$. The first and the second derivatives of the exact Fourier coefficients are given by:

[TABLE]

while the Fourier coefficients of the derivatives obtained from the differencing schemes are:

[TABLE]

where $w^{\prime}=w^{\prime}(w)$ and $w^{\prime\prime}=w^{\prime\prime}(w)$ are the modified wavenumbers for the first and second order derivatives, respectively.

The numerical wavenumbers of the original combined compact scheme are given by [3]

[TABLE]

and

[TABLE]

for stencils (3) and (4) (Chu and Fan [3]) and

[TABLE]

and

[TABLE]

for stencils (5) and (6) (Mahesh [15]).

For the prefactored CCD schemes the modified wavenumbers $w_{F}^{\prime}(w)$, $w_{B}^{\prime}(w)$, $w_{F}^{\prime\prime}(w)$ and $w_{B}^{\prime\prime}(w)$ can be determined from the equations obtained by applying the Fourier transform to (10)-(13):

[TABLE]

[TABLE]

for the forward operators and

[TABLE]

[TABLE]

for the backward operators. Equations (21) and (22) can be written in matriceal form as

[TABLE]

where $[A]$ is a $4\times 4$ matrix given by

[TABLE]

and $\{X\}$ is the unknown vector of 4 elements given by

[TABLE]

where $\Re()$ and $\Im()$ stands for real and imaginary parts, respectively. $\{R\}$ in (25) is the right hand side vector given by

[TABLE]

The solution to the set of equations (25) gives the real parts of the numerical wavenumbers in the form:

[TABLE]

and

[TABLE]

where $f_{1}^{I}$, $f_{2}^{I}$, $f_{3}^{I}$, $g_{1}$, $g_{2}$, $g_{3}$, $f_{1}^{II}$, $f_{2}^{II}$, $f_{3}^{II}$, $f_{4}^{II}$ are functions of the weights in equations (21) and (22). Comparing the equations (35) and (36) to the equations (19) and (20), the weights in equations (21) and (22) can be determined by solving the next polynomial equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The next open problem is to solve the above closed set of polynomial equations by an iterative method to obtain the weights of the scheme in (21) and (22).