Dispersion Relation Preserving Combined Compact Difference Schemes for Flow Problems

Abstract

In this work, we have proposed two new combined compact difference (CCD) schemes for the solution of Navier–Stokes equations. These spatial discretization schemes have not only high spectral resolution for obtaining first and second derivative terms, but also have improved dispersion relation preserving properties when the fourth-order four-stage Runge–Kutta scheme is used for time integration. Out of the two proposed CCD schemes, the first scheme has upwind stencil, while the second scheme has a central stencil. Important numerical properties of these schemes have been analyzed and their effectiveness have been shown by solving the model wave equations, as well as Navier–Stokes equations. Results show that the upwind CCD scheme is suitable for high accuracy large eddy simulation of transitional and turbulent flowfields.

Appendix

In this section, we have provided the steps while developing an upwind combined compact difference scheme. The first and the second derivative terms ($\frac{\partial u}{\partial x}$ and $\frac{\partial ^2 u}{\partial x^2}$) are approximated as

$$\begin{aligned} a_1\frac{\partial u}{\partial x}|_{i-1} \!+\!\frac{\partial u}{\partial x}|_{i}\!+\!a_3\frac{\partial u}{\partial x}|_{i+1}&= \frac{1}{h}( c_{1} u_{i-2} \!+\! c_{2}u_{i-1}\!+\! c_{3}u_{i} )\nonumber \\&-\,h\left( b_1\frac{\partial ^2 u}{\partial x^2}|_{i-1}\!+\!b_2\frac{\partial ^2 u}{\partial x^2}|_{i} \!+\!b_3\frac{\partial ^2 u}{\partial x^2}|_{i+1} \right) , \end{aligned}$$

(30)

$$\begin{aligned} -\frac{1}{8}\frac{\partial ^2 u}{\partial x^2}|_{i-1} \!+\!\frac{\partial ^2 u}{\partial x^2}|_{i} \!-\!\frac{1}{8}\frac{\partial ^2 u}{\partial x^2}|_{i+1}&= \frac{3}{h^2}( u_{i-1} \!-\!2u_{i} \!+\! u_{i+1} ) \!-\!\frac{9}{8h}\left( -\frac{\partial u}{\partial x}|_{i-1} \!+\! \frac{\partial u}{\partial x}|_{i+1}\right) .\nonumber \\ \end{aligned}$$

(31)

Here, we have used a central stencil to approximate the second-order derivative term. Taylor series expansion has been used to obtain the coefficients shown in (31), which results in a sixth-order accurate stencil based on the leading truncation error term. Coefficients in Eq. (30) are obtained by applying the Taylor series expansion for the terms $u_{i-1},\,u_{i+1},\,\frac{\partial u}{\partial x}|_{i-1},\,\frac{\partial u}{\partial x}|_{i},\,\frac{\partial u}{\partial x}|_{i+1},\,\frac{\partial ^{2} u}{\partial x^{2}}|_{i-1},\,\frac{\partial ^{2} u}{\partial x^{2}}|_{i}$ and $\frac{\partial ^{2} u}{\partial x^{2}}|_{i+1}$ with respect to $u_i$.

In Eq. (30), there are in total eight unknown coefficients. Elimination of the leading truncation error terms gives the following set of seven algebraic equations.

$$\begin{aligned}&c_1 + c_2 + c_3 = 0, \end{aligned}$$

(32)

$$\begin{aligned}&-2c_1-c_2-a_1-a_3=1, \end{aligned}$$

(33)

$$\begin{aligned}&4c_1+c_2+2a_1-2a_3-2b_1-2b_2-2b_3 = 0, \end{aligned}$$

(34)

$$\begin{aligned}&8c_1+c_2+3a_1+3a_3-6b_1+6b_3 = 0, \end{aligned}$$

(35)

$$\begin{aligned}&16c_1+c_2+4a_1-4a_3-12b_1-12b_3 = 0, \end{aligned}$$

(36)

$$\begin{aligned}&32c_1+c_2+5a_1+5a_3-20b_1+20b_3 = 0, \end{aligned}$$

(37)

$$\begin{aligned}&64c_1+c_2+6a_1-6a_3-30b_1-30b_3 = 0. \end{aligned}$$

(38)

Derivation of one more algebraic equation is needed to determine all the eight coefficients in Eq. (30) uniquely for $\frac{\partial u}{\partial x}|_{i}$. The first-order derivative term in Eq. (30) can be better approximated if the dispersive nature of the term $\frac{\partial u}{\partial x}$ is well retained [1].

The expressions of the actual wavenumber for Eqs. (30) and (31) can be derived as

$$\begin{aligned} \mathbf {i}\alpha h~( a_1 e^{-\mathbf {i}\alpha h} + 1 + a_3 e^{\mathbf {i}\alpha h })&\simeq \left( c_1e^{-2\mathbf {i}\alpha h}+c_2e^{-\mathbf {i}\alpha h}+c_3\right) \nonumber \\&- \left( \mathbf {i}\alpha h\right) ^2\left( b_1e^{-\mathbf {i}\alpha h} + b_2 + b_3e^{\mathbf {i}\alpha h}\right) , \end{aligned}$$

(39)

$$\begin{aligned} (\mathbf {i}\alpha h)^2\left( -\frac{1}{8}e^{-\mathbf {i}\alpha h} +1-\frac{1}{8} e^{\mathbf {i}\alpha h}\right) \simeq ~ (3e^{-\mathbf {i}\alpha h}-6+3 e^{\mathbf {i}\alpha h}) \nonumber \\ -\,\mathbf {i}\alpha h~\left( -\frac{9}{8}e^{-\mathbf {i}\alpha h} +\frac{9}{8}e^{\mathbf {i}\alpha h}\right) . \end{aligned}$$

(40)

While approximating the first-order derivative term, dispersion error is minimized if the exact and the numerical wavenumbers are matched excellently over a complete wavenumber range. This amounts to equating the effective wavenumbers $\alpha ^{'}$ and $\alpha ^{''}$ to those shown in the right-hand sides of Eqs. (41) and (42) [1]. Thus, we obtained following two equations

$$\begin{aligned}&\mathbf {i}\alpha ^{'} h~\left( a_1 e^{-\mathbf {i}\alpha h} + 1 + a_3 e^{\mathbf {i}\alpha h }\right) =\left( c_1e^{-2\mathbf {i}\alpha h}+c_2e^{-\mathbf {i}\alpha h}+c_3\right) \nonumber \\&\quad -\,\left( \mathbf {i}\alpha ^{''} h\right) ^2\left( b_1e^{-\mathbf {i}\alpha h} + b_2 + b_3e^{\mathbf {i}\alpha h}\right) ,\end{aligned}$$

(41)

$$\begin{aligned}&(\mathbf {i}\alpha ^{''} h)^2\left( -\frac{1}{8}e^{-\mathbf {i}\alpha h} +1-\frac{1}{8} e^{\mathbf {i}\alpha h}\right) =(3e^{-\mathbf {i}\alpha h}-6+3 e^{\mathbf {i}\alpha h})\nonumber \\&\quad -\,\mathbf {i}\alpha ^{'} h~\left( -\frac{9}{8}e^{-\mathbf {i}\alpha h} +\frac{9}{8}e^{\mathbf {i}\alpha h}\right) . \end{aligned}$$

(42)

Equations (41) and (42) are solved to get the expression for $\alpha ' h$ which has been used subsequently to minimize the dispersion error. This expression for $\alpha ' h$, which is in general complex with the real and imaginary parts of the numerical modified wavenumber $\alpha 'h$, provides information regarding the dispersion error (phase error) and the dissipation error (amplitude error), respectively. For getting better dispersive accuracy of $\alpha '$, the value of $\alpha h$ should be closer to $\mathfrak {R}[\alpha ' h]$, where $\mathfrak {R}[\alpha ' h]$ denotes the real part of $\alpha ' h$. Thus the error function $E(\alpha )$ as defined below should be very small. It has been evaluated over the integration range given below as

$$\begin{aligned} E(\alpha ) = \int \limits _{0}^{\frac{7\pi }{8}} \left[ \left( \alpha \,h-\mathfrak {R}[\alpha '\,h]\right) \right] ^2 d(\alpha h). \end{aligned}$$

(43)

To make the error function defined in $0\le \alpha h\le \frac{7\pi }{8}$ to be positive and minimal, the extreme condition $\frac{\partial E}{\partial c_{3}} = 0$ is enforced to minimize the numerical wavenumber error. This constraint equation has been used together with previously derived seven algebraic equations to obtain all the eight unknowns. The resulting eight introduced unknown coefficients can be uniquely determined as

$$\begin{aligned} a_1&= 0.888251792581,\,a_3 = 0.049229651564, b_1 = 0.150072398996,\\ b_2&= -0.250712794122,\, b_3 = -0.012416467490,\,c_1 = 0.016661718438,\\ c_2&= -1.970804881023 \,and\, c_3 = 1.954143162584. \end{aligned}$$

For $u<0$, the proposed three-point stencil non-centered combined compact difference scheme can be similarly derived below for the approximation of the derivative term $\frac{\partial u}{\partial x}$

$$\begin{aligned}&0.049229651564\frac{\partial \upphi }{\partial x}|_{i-1} +\frac{\partial \upphi }{\partial x}|_{i} +0.888251792581\frac{\partial \upphi }{\partial x}|_{i+1}\nonumber \\&\quad +\,h \Big (0.012416467490\frac{\partial ^2 \upphi }{\partial x^2}|_{i-1} +0.250712794122\frac{\partial ^2 \upphi }{\partial x^2}|_{i} -0.150072398996\frac{\partial ^2 \upphi }{\partial x^2}|_{i+1}\Big ) \nonumber \\&\quad =\frac{1}{h}\left( -1.954143162584{\upphi }_{i} + 1.970804881023{\upphi }_{i+1} - 0.016661718438{\upphi }_{i+2}\right) . \end{aligned}$$

(44)

One should pay a careful attention to the following important points. Significant improvements in spectral resolution as shown in Fig. 1e and the DRP region shown in Fig. 5 are possible using a small stencil because we have derived the scheme for first and second derivatives in a coupled fashion. The expression for $\alpha ^{'} h$ is derived by considering Eqs. (41) and (42) together. Figure 1f shows that efficiency in the evaluation of the second derivative is very close to the exact value. Thus, while designing CCD schemes, use of one optimization equation given by $\frac{\partial E}{\partial c_{3}} = 0$ is justified.