Two Dimensional Analysis of Hybrid Spectral/Finite Difference Schemes for Linearized Compressible Navier–Stokes Equations

Abstract

This study seeks to compare different combinations of spatial dicretization methods under a coupled spatial temporal framework in two dimensional wavenumber space. The aim is to understand the effect of dispersion and dissipation on both the convection and diffusion terms found in the two dimensional linearized compressible Navier–Stokes Equations (LCNSE) when a hybrid finite difference/Fourier spectral scheme is used in the x and y directions. In two dimensional wavespace, the spectral resolution becomes a function of both the wavenumber and the wave propagation angle, the orientation of the wave front with respect to the grid. At sufficiently low CFL number where temporal discretization effects can be neglected, we show that a hybrid finite difference/Fourier spectral schemes is more accurate than a full finite difference method for the two dimensional advection equation, but that this is not so in the case of the LCNSE. Group velocities, phase velocities as well as numerical amplification factor were used to quantify the numerical anisotropy of the dispersion and dissipation properties. Unlike the advection equation, the dispersion relation representing the acoustic modes of the LCNSE contains an acoustic terms in addition to its advection and viscous terms. This makes the group velocity in each spatial direction a function of the wavenumber in both spatial directions. This can lead to conditions for which a hybrid Fourier spectral/finite difference method can become less or more accurate than a full finite difference method. To better understand the comparison of the dispersion properties between a hybrid and full FD scheme, the integrated sum of the error between the numerical group velocity $V^{*}_{grp,full}$ and the exact solution across all wavenumbers for a range of wave propagation angle is examined. In the comparison between a hybrid and full FD discretization schemes, the fourth order central (CDS4), fourth order dispersion relation preserving (DRP4) and sixth order central compact (CCOM6) schemes share the same characteristics. At low wave propagation angle, the integrated errors of the full FD and hybrid discretization schemes remain the same. At intermediate wave propagation angle, the integrated error of the full FD schemes become smaller than that of the hybrid scheme. At large wave propagation angle, the integrated error of the full FD schemes diverges while the integrated error of the hybrid discretization schemes converge to zero. At high reduced wavenumber and sufficiently low CFL number where temporal discretization error can be neglected, it was found that the numerical dissipation of the viscous term based on the CDS4, DRP4, CCOM6 and isotropy optimized CDS4 schemes ($\hbox {CDS4}_{{opt}}$) schemes was lower than the actual physical dissipation, which is only a function of the cell Reynolds number. The wave propagation angle at which the numerical dissipation of the viscous term approaches its maximum occurs at $\pi /4$ for the CDS4, DRP4, CCOM6 and $\hbox {CDS4}_{{opt}}$ schemes.

1 Introduction

In direct numerical simulation of turbulent flow, the ability to resolve all spatial and temporal scales of turbulence is directly associated with the choice of grid size, time step size as well as the resolution characteristic of the numerical schemes. Furthermore, when one solves a non linear equation such as the Navier–Stokes equation, the direct discretization of the non linear terms can cause aliasing errors, leading to non linear instability issues. Although de-aliasing techniques such as a low pass filter can be used to damp/remove the high wavenumber content which are the leading sources of aliasing error, this may inadvertently filter out the physically correct high wavenumber content. In addition, the accurate simulation of aero-acoustics problem requires the solver to capture the compressibility effects by correctly estimating the pressure fluctuation as well as having sufficient accuracy to propagate information from the source to the far field correctly [14]. Since acoustic fluctuations are very weak as compared to the aerodynamic ones, this necessitates the use of high order accurate schemes with low dispersive and dissipative properties. For these purposes, it is important that high resolution schemes are used to obtain accurate results. This motivates the development of tools which can be used to optimise or analyse numerical schemes for their appropriate usage.

Assessment of numerical schemes based on the use of Fourier analysis is well established and its earliest ideas date back to [30, 32]. The use of the 1D advection equation provides a suitable model for the calibration of schemes initially. However, assessment of schemes based on the 1D advection equation does not provide the full picture in terms of their resolution characteristics, especially when one considers solving Navier–Stokes Equations in two-dimensional space. Moreover, numerical anisotropy is introduced when finite difference (FD) schemes are applied in multi-dimensional problems. It represents a phenomenon for which the resolution characteristic of the numerical scheme differs in different spatial directions. The resolution characteristic of the numerical scheme becomes a function of both the reduced wavenumber ($k \varDelta$) as well as the wave propagation angle, $\theta _{w}$. The wave propagation angle represents the orientation of the wave front with respect to the grid and it does not necessary coincide with the convection angle ($\theta _{a}$) of the velocity components (see Fig. 1).

From an error dynamics point of view, the formal accuracy of a numerical scheme determined through its leading order truncation error is often limited because no information is known in terms of the phase or amplitude error as a function of the wavenumber. As such, Fourier analysis of FD schemes was introduced in [32]. Fourier analysis of FD schemes in one dimensional wavenumber space has been extensively studied in majority of the literature such as [5, 6, 10, 12, 13, 15, 16, 20, 22, 24, 27, 28, 33,34,35]. The analysis of numerical schemes under full discretization can be found in [4, 20, 21, 24], from which group velocities, phase velocities and numerical amplification factor of the one dimensional advection equation was analysed under different spatial and temporal discretization schemes. Furthermore, the analysis of the one dimensional advection diffusion equation can be found in [27]. Studies concerning two dimensional wavenumber analysis can be found in [13] and [26] but they were limited to spatial discretization considerations. A further study [23] by the same author has looked into numerical anisotropy under full discretisation in two dimensional wavenumber space. The study assesses the effect of changing grid aspect ratio and propagation angle on the solution of the linearized rotating shallow water equations and the two dimensional advection equation. It was conclusively shown that wave-aligned grids are superior over grids with equal spacing in both directions. However, the study did not consider the effect of hybrid Fourier spectral/finite difference spatial discretization and primarily considered the convection term. As for the anisotropy error of the diffusion term, the literature which discuses particularly about the error of the diffusion term can be found in [18] but the study was limited to the two dimensional diffusion equation and the use of finite difference schemes in both spatial directions. Further analysis of the diffusion term can be found in [31] where the authors have shown that the focusing phenomenon, a condition where numerical methods producing acceptable results for a long time abruptly blow up, was due to the anti-diffusion caused by the FD discretization of the diffusion term. The importance of understanding the error of the diffusion term stems from the presence of the diffusion term in parabolic PDEs which is used to model many important physical phenomenon such as heat diffusion and boundary layer flow over a flat plate. An additional consideration in this analysis is related to the use of non uniform grid in the simulation of practical problems. Traditionally, this can be accounted for through the grid transformation metric in the analysis equation [3]. However the finite difference discretization is usually applied on a uniform computational grid as the Navier–Stokes Equations are integrated in time. Therefore, the accumulation of discretization error is mainly considered in the case of a uniform grid. It should be noted that the discretization error coming from the grid transformation metric terms accumulates only in the pre and post processing step for a non moving grid simulation problem. As such, the analysis in this study will only consider a uniform grid.

To address the research gaps mentioned above, the dispersion relation of the two dimensional LCNSE was used to analyse the dispersion and dissipation property of the numerical schemes in two dimensional wavenumber space under full discretization. The use of the dispersion relation of the LCNSE allows for the coupled analysis of the advective, diffusive and acoustic effects. The aim was to better understand its implication on numerical anisotropy when different spatial discretization methods are used in both spatial directions for the convection and diffusion terms of the Navier–Stokes Equations. In particular, the case comparing a hybrid FD/FS method to a full FD method in two spatial directions is scrutinized.

Firstly, the modified wavenumber relation of various FD schemes are introduced. This is followed by a Fourier analysis of the 1D advection equation under full discretization from which the concept of group velocity, phase velocity and numerical amplification factor are introduced. The Fourier analysis of the 2D advection equation under full discretization was also analysed in order to showcase the resolution characteristic of the numerical schemes in two-dimensional space. Then the dispersion relations of the LCNSE were presented. The dispersion and dissipation property of the coupled spatial temporal system for different combinations of spatial discretization schemes for both the convective and viscous terms are discussed. For verification purposes, physical and numerical group velocities calculated based on the dispersion relations of the LCNSE are then compared with numerical solutions of the acoustic waves under different initial conditions.

Fig. 1

Illustration of wave propagation angle ($\theta _{w}$) and convection angle ($\theta _{a} = \tan ^{-1} (\frac{v}{u})$)

2 Modified Wavenumber of the Finite Difference and Fourier Spectral Schemes

In this paper, finite difference (FD) and Fourier spectral (FS) schemes are used in combination. For the FD schemes considered, standard central difference, dispersion relation preserving (DRP) schemes, central compact schemes [2] and isotropic optimized central difference schemes [25] are used. These schemes are introduced and their key characteristics are discussed. An example of the standard central difference scheme is the fourth order central differencing scheme (CDS4):

$$\begin{aligned} \varPhi ^{'}_{i} = \frac{\partial \varPhi _{i}}{\partial x} \approx \frac{8(\varPhi _{i+1} - \varPhi _{i-1}) -1(\varPhi _{i+2} - \varPhi _{i-2}) }{12 \ \varDelta x}. \end{aligned}$$

(1)

The modified wavenumber relation of central type schemes does not have an imaginary component since the use of the same number of nodal points with equal but opposite signs coefficients on each side of its stencils results in the cancellation of its imaginary component. This implies that the spatial discretization scheme does not provide any numerical dissipation. Numerical schemes with better dispersion properties can be achieved by minimizing the dispersion error through the optimization of the coefficients based on the dispersion relation of the one dimensional advection equation. Such schemes are known as the dispersion relation preserving schemes in literature. The fourth order DRP scheme (DRP4) [28] is given by:

$$\begin{aligned} \begin{aligned} \varPhi ^{'}_{i} = \frac{\partial \varPhi _{i}}{\partial x} \approx&\frac{1}{\ \varDelta x} \Big ( 0.79926643(\varPhi _{i+1} - \varPhi _{i-1}) -0.18941314(\varPhi _{i+2} - \varPhi _{i-2}) \\&+ 0.02651995(\varPhi _{i+3} - \varPhi _{i-3}) \Big ). \end{aligned} \end{aligned}$$

(2)

Such optimization considers only the dispersion error of the spatial discretization scheme and are therefore applied on central type schemes. The only drawback of such an approach is the loss in the formal order of accuracy. Another approach to improve the dispersion property is to increase the order of the truncation error through the use of more neighbouring nodes. However, the amount of improvement diminishes as the order of accuracy is increased. To alleviate this issue, compact schemes were introduced [2]. The gain in “resolution accuracy” of the compact schemes is the consequence of having an effective stencil that virtually ranges over the considered computational domain, since a linear system of equations has to be solved to get the derivative vector. Hence, its characteristics tend to the behaviour of a discretization with more global than local features [11]. Compact schemes are also termed implicit spatial discretization method since the derivatives are calculated implicitly through a system of equations involving penta/tri diagonal (banded) matrices. Their computational cost is relatively higher as compared to explicit FD schemes involving the same stencil width. However, their higher computational cost is offset by the provision of better spectral accuracy. It should be noted that compact FD schemes can be far less efficient than explicit FD schemes without the use of appropriate parallelization [7,8,9]. In the context of parallelization, compact FD schemes can be further classified into coupled and decoupled approaches depending on whether the operation of numerical differentiation inside a subdomain or processor is independent of other parts of the domain at a particular computational step. The first group (coupled approach) focuses on the method of solving a linear system of the full domain in parallel while the second group (decoupled approach) essentially decouples compact schemes to enable them to be solved independently on each processor. The decoupling methods inevitably introduce interdomain interfaces, on which a boundary method is required to close the compact scheme inside of each subdomain [7]. An example of the standard compact scheme is the sixth order central compact difference scheme (CCOM6) [12]:

$$\begin{aligned} \frac{1}{3}\varPhi ^{'}_{i-1} + \varPhi ^{'}_{i} + \frac{1}{3}\varPhi ^{'}_{i+1} \approx \frac{1}{9} \frac{\varPhi _{i+2} - \varPhi _{i-2}}{4 \varDelta x} + \frac{14}{9} \frac{\varPhi _{i+1} - \varPhi _{i-1}}{2 \varDelta x}. \end{aligned}$$

(3)

Further variants of the split type compact schemes can be found in [11] where the scheme is derived for a mixed Fourier Spectral/FD method that is designed for the direct numerical simulation (DNS) of boundary layer transition and turbulence. Split type compact schemes are essentially a combination of positively and negatively biased stencils, that add-up to a central FD. Their main purpose is to use the inherent numerical dissipation of the biased FD sub schemes to remove the aliasing error (at higher wavenumber) caused by the discretization of a non linear term. Finally, the last type of FD scheme considered is the fourth order isotropic optimized central difference scheme [25]. This type of scheme is derived from the standard fourth order central difference (CDS4) scheme and it becomes identical to the standard CDS4 scheme for 1D problems. The key difference as compared to the standard CDS4 scheme lies in the use of nodal points from more than one direction. Such a scheme is derived based on the dependency of the derivatives on both x and y derivatives and the condition that the phase or group velocities are the same on some specific spatial direction [25]. The fourth order isotropic optimized central difference scheme ($CDS4_{opt}$) is expressed as:

$$\begin{aligned}&\begin{aligned} \left( \frac{\partial \varPhi }{\partial x} \right) _{i,j} \approx&\frac{1}{12 \varDelta x (1 + \beta )} \Bigl [ \varPhi _{i-2,j} - 8\varPhi _{i-1,j} + 8\varPhi _{i+1,j} - \varPhi _{i+2,j} \\&+ \frac{\beta }{2}( \varPhi _{i-2,j-2} - 8\varPhi _{i-1,j-1} + 8\varPhi _{i+1,j+1} - \varPhi _{i+2,j+2} - \varPhi _{i+2,j-2} \\&+ 8\varPhi _{i+1,j-1} - 8\varPhi _{i-1,j+1} + \varPhi _{i-2,j+2}) \Bigr ], \end{aligned} \end{aligned}$$

(4)

$$\begin{aligned}&\begin{aligned} \left( \frac{\partial \varPhi }{\partial y} \right) _{i,j} \approx&\frac{1}{12 \varDelta y (1 + \beta )} \Bigl [ \varPhi _{i,j-2} - 8\varPhi _{i,j-1} + 8\varPhi _{i,j+1} - \varPhi _{i,j+2} \\&+ \frac{\beta }{2}( \varPhi _{i-2,j-2} - 8\varPhi _{i-1,j-1} + 8\varPhi _{i+1,j+1} - \varPhi _{i+2,j+2} + \varPhi _{i+2,j-2} \\&- 8\varPhi _{i+1,j-1} + 8\varPhi _{i-1,j+1} - \varPhi _{i-2,j+2}) \Bigr ], \end{aligned} \end{aligned}$$

(5)

where $\beta$ is the isotropy correction factor. It is found to be 0.26 for the fourth order central difference scheme [25]. The calculation of the modified wavenumber ($k^{*}$) of the FD schemes described above are presented in the appendix section. A comparison of the dispersion property of the different spatial differencing scheme is provided in Fig. 2. As seen in Fig. 2, the phase error of the CDS4 scheme is represented by its difference with respect to the black solid line (Fourier spectral scheme) along the entire spectrum of wavenumbers. The Fourier spectral scheme represents the ideal condition for which the modified wavenumbers matches the actual wavenumbers exactly ($\mathfrak {R}(k^{*}_{x} \varDelta x)$ = $\mathfrak {R}(k_{x} \varDelta x)$). It can be seen that compact difference scheme provides the best dispersion property among the different FD schemes considered. For the $CDS4_{opt}$ scheme, the modified wavenumber relation is a function of the reduced wavenumber in both spatial directions. In the case where $k_{y} \varDelta y$ is 0, the modified wavenumber relation of the standard CDS4 scheme is recovered. In the case where $k_{y} \varDelta y > 0$, the dispersion property of the $CDS4_{opt}$ scheme becomes poorer than that of the standard CDS4 scheme.

Fig. 2

Comparison of dispersion, $\mathfrak {R}(k_{x}^{*} \varDelta x)$ property for different FD schemes

3 Fourier Analysis of Advection Equation Under Full Discretization

The modified wavenumber relation of the spatial differencing schemes was established in the previous section. These relations provide an idea of the dispersion and dissipation property of the spatial discretization scheme. In order to assess the accuracy of the numerical schemes under full discretization, the numerical amplification factor of the temporal discretization scheme is introduced.

3.1 1D Advection

To illustrate this, the linear model equation for inviscid flow (1D advection) will be used. Equation (6) is a non-dispersive equation that convects the initial solution to the right at a group velocity, ${\overline{u}}$ equal to the phase speed at all times. In solving Eq. (6) numerically, the spatial and temporal derivatives are replaced by numerical approximations. The field variable, $\varPhi$, is written in terms of its inverse Fourier transform:

$$\begin{aligned}&\frac{\partial \varPhi }{\partial t} + {\overline{u}} \frac{\partial \varPhi }{\partial x} = 0. \end{aligned}$$

(6)

$$\begin{aligned}&\varPhi (x,t) = \frac{1}{2 \pi } \int _{-\pi /\varDelta x}^{\pi /\varDelta x} {\widehat{\varPhi }} (k_{x},t)e^{ik_{x}x} dk_{x}, \end{aligned}$$

(7)

Substituting Eq. (7) into (6) and noting that the introduction of the inverse Fourier transform of the field variable into the numerical approximation of the spatial derivative term introduces a modified wavenumber $k_{x}^{*}$ [32], the resulting equation is:

$$\begin{aligned} \frac{\partial {\widehat{\varPhi }}(k_{x},t) }{\partial t} = -i {\overline{u}} k_{x}^{*} {\widehat{\varPhi }}(k_{x},t). \end{aligned}$$

(8)

In this work, the 4th order Runge Kutta (RK4) temporal scheme is considered [17]. Using the semi-discretized 1D advection equation expressed in Fourier space (Eq. 8), the time derivative term can be expressed as follows:

$$\begin{aligned} \frac{\partial {\widehat{\varPhi }} }{\partial t} = -i {\overline{u}} k_{x}^{*} {\widehat{\varPhi }} = -i \omega ^{*}_{semi} {\widehat{\varPhi }}, \end{aligned}$$

(9)

where $\omega ^{*}_{semi} = {\overline{u}} k_{x}^{*}$ is the semi-discretized form of the numerical dispersion relation. It does not take into account the effect of the temporal discretization. In order to account for the temporal discretization effect, Eq. (9) is substituted back into the numerical amplification factor of the RK4 scheme [17]. The resulting numerical amplification factor, Z of the RK4 scheme can thus be expressed in terms of $\omega ^{*}_{semi}$ as follows:

$$\begin{aligned} \frac{{\widehat{\varPhi }}(k_{x}, t + \varDelta t)}{{\widehat{\varPhi }} (k_{x},t)} = \underbrace{(1 - i \varDelta t \omega ^{*}_{semi} + \frac{1}{2} (i \varDelta t \omega ^{*}_{semi})^{2} - \frac{1}{6} (i \varDelta t \omega ^{*}_{semi})^{3} + \frac{1}{24} (i \varDelta t \omega ^{*}_{semi})^{4} )}_{Z}, \nonumber \\ \end{aligned}$$

(10)

where

$$\begin{aligned} \omega ^{*}_{semi} \varDelta t = {\overline{u}}k^{*}_{x} \varDelta t = Rk^{*}_{x} \varDelta x, \quad R = \frac{{\overline{u}} \varDelta t}{\varDelta x}. \end{aligned}$$

(11)

$\textit{R}$ represents the CFL number. The solution, ${\widehat{\varPhi }}$ in the next instance is represented by the product of a complex number, $\textit{Z}$ and the solution of the previous timestep. Using this information and Eq. (7), the general solution of $\varPhi$ considered after $\textit{n}$ time steps can be written as follows:

$$\begin{aligned} \varPhi (x,t_{0} + n\varDelta t) = \frac{1}{2 \pi } \int _{-\pi /\varDelta x}^{\pi /\varDelta x} Z^{n} \ {\widehat{\varPhi }} (k_{x},0)e^{ik_{x} x} dk_{x}, \end{aligned}$$

(12)

where $\textit{Z}$ is a complex number which can be expressed in terms of its modulus and argument [32], $Z = |Z|e^{-i\psi }$. The resulting equation is:

$$\begin{aligned} \varPhi (x,t_{0} + n\varDelta t) = \frac{1}{2 \pi } \int _{-\pi /\varDelta x}^{\pi /\varDelta x} |Z|^{n} \ {\widehat{\varPhi }} (k_{x},0)e^{i(k_{x}x - n \psi )} dk_{x}. \end{aligned}$$

(13)

The general numerical solution to Eq. (6) [21] is given as follows:

$$\begin{aligned} \varPhi (x,t) = \frac{1}{2 \pi } \int _{-\pi /\varDelta x}^{\pi /\varDelta x} |Z|^{t/\varDelta t}{\widehat{\varPhi }}(k_{x},0) e^{i(k_{x}x - U^{*}_{phase, full} k_{x}t)} dk_{x} \end{aligned}$$

(14)

Instead of having a constant physical phase speed, ${\overline{u}}$ as in the case of an exact solution, the numerical solution is represented by a wavenumber dependent numerical phase speed, $U^{*}_{phase, full}$. This is because the numerical solution is dispersive [21]. By comparing the numerical solution to the exact solution of the 1D advection equation:

$$\begin{aligned} \begin{aligned} \varPhi (x,t)&= \frac{1}{2 \pi } \int _{-\pi /\varDelta x}^{\pi /\varDelta x} |Z|^{t/\varDelta t} {\widehat{\varPhi }}(k_{x},0) e^{i(k_{x}x - U^{*}_{phase, full} k_{x}t)} dk_{x} \\&= \frac{1}{2 \pi } \int _{-\pi /\varDelta x}^{\pi /\varDelta x} |Z|^{t/\varDelta t} {\widehat{\varPhi }}(k_{x},0) e^{i(k_{x}x - \omega t)} dk_{x}, \end{aligned} \end{aligned}$$

(15)

it can be seen that the numerical dispersion relation is given as $\omega = U^{*}_{phase, full} k_{x}$. The numerical dispersion relation in this form considers both space and time discretization and will be denoted as $\omega ^{*}_{full}$ thereafter. The relation between $\omega ^{*}_{full}$ and the numerical amplification factor, Z is established next. By comparing Eq. (15) to Eq. (13), it can be seen that $\omega ^{*}_{full}$ of the 1D advection equation is related to the phase shift of the amplification factor via:

$$\begin{aligned} \begin{aligned} \omega ^{*}_{full} t&= U^{*}_{phase, full} k_{x} n \varDelta t = n \psi , \\ \omega ^{*}_{full}&= U^{*}_{phase, full} k_{x} = \frac{\psi }{\varDelta t}, \\ \frac{\omega ^{*}_{full}}{k_{x}}&= U^{*}_{phase, full} = \frac{\psi }{k_{x} \varDelta t}. \end{aligned} \end{aligned}$$

(16)

The physical dispersion relation of the 1D advection is given by $\omega = {\overline{u}} k_{x}$ and the corresponding physical phase speed is calculated as $U_{phase, exact} = \frac{\omega }{k_{x}} = {\overline{u}}$. The normalized numerical phase speed can be expressed as follows:

$$\begin{aligned} \frac{U^{*}_{phase, full}}{U_{phase, exact}} = \frac{\psi }{{\overline{u}} k_{x} \varDelta t} = \frac{\psi }{R k_{x} \varDelta x}. \end{aligned}$$

(17)

The normalized numerical group velocity based on full discretization consideration is obtained using the fully discretized numerical dispersion relation as follows:

$$\begin{aligned} \frac{U^{*}_{grp, full}}{U_{grp, exact}} = \frac{1}{{\overline{u}}} \frac{\partial \omega ^{*}_{full}}{\partial k_{x}} = \frac{1}{{\overline{u}} \varDelta t} \frac{\partial \psi }{\partial k_{x}} = \frac{1}{R \varDelta x} \frac{\partial \psi }{\partial k_{x}}, \end{aligned}$$

(18)

where $\psi = - \tan ^{-1} \bigg ( \frac{Im(Z(\varDelta t \omega ^{*}_{semi}))}{Real(Z(\varDelta t \omega ^{*}_{semi}))} \bigg )$. It should be noted that $U_{grp, exact} = \frac{\partial \omega }{\partial k_{x}} ={\overline{u}}$ for the 1D advection equation. For group velocity calculated under full discretization, the fully discretized numerical dispersion relation is considered. In order to illustrate this, plots representing the modulus of the numerical amplification factor (|Z|), normalized group velocities ($\frac{U^{*}_{grp, full}}{{\overline{u}}}$) as well as the normalised numerical phase speed ($\frac{U^{*}_{phase, full}}{{\overline{u}}}$) for the CDS4, DRP4, CCOM6 and FS schemes along with the RK4 temporal discretization schemes are presented in Figs. 3, 4, 5 and 6. These figures are plotted as a function of the CFL number and reduced wavenumber. ${\overline{u}}$ is defined as 1. The normalized group velocity plots are all scaled similarly from 0 to 5. The normalized phase velocity plots are all scaled similarly from 0 to 2 for comparison sake. For the normalized group velocity plots, hatched regions are denoted as the q-wave region where $\frac{U^{*}_{grp, full}}{{\overline{u}}} < 0$ [19]. q-waves are essentially non-physical waves which arises due to numerical methods. They represent numerical waves which travel in opposite direction to the physical waves.

Fig. 3

Comparison of a numerical amplification factor (|Z|), b normalized numerical group velocity ($\frac{U^{*}_{grp, full}}{{\overline{u}}}$) and c normalized numerical phase velocity ($\frac{U^{*}_{phase, full}}{{\overline{u}}}$) for CDS4 spatial and RK4 temporal scheme. Black contour lines are ranged from 0.99 to 1 for figure (a). Black contour lines are ranged from −1 to 1 for figure (b). Black contour lines are ranged from 0 to 2 for figure (c). Hatched region in figure (b) represents the q-wave region

Fig. 4

Comparison of a numerical amplification factor (|Z|), b normalized numerical group velocity ($\frac{U^{*}_{grp, full}}{{\overline{u}}}$) and c normalized numerical phase velocity ($\frac{U^{*}_{phase, full}}{{\overline{u}}}$) for DRP4 spatial and RK4 temporal scheme. Black contour lines are ranged from 0.99 to 1 for figure (a). Black contour lines are ranged from −1 to 1 for figure (b). Black contour lines are ranged from 0 to 2 for figure (c). Hatched region in figure (b) represents the q-wave region

Fig. 5

Comparison of a numerical amplification factor (|Z|), b normalized numerical group velocity ($\frac{U^{*}_{grp, full}}{{\overline{u}}}$) and c normalized numerical phase velocity ($\frac{U^{*}_{phase, full}}{{\overline{u}}}$) for CCOM6 spatial and RK4 temporal scheme. Black contour lines are ranged from 0.99 to 1 for figure (a). Black contour lines are ranged from −1 to 1 for figure (b). Black contour lines are ranged from 0 to 2 for figure (c)

Fig. 6

Comparison of a numerical amplification factor (|Z|), b normalized numerical group velocity ($\frac{U^{*}_{grp, full}}{{\overline{u}}}$) and c normalized numerical phase velocity ($\frac{U^{*}_{phase, full}}{{\overline{u}}}$) for FS spatial and RK4 temporal scheme. Black contour lines are ranged from 0.99 to 1 for figure (a). Black contour lines are ranged from −1 to 1 for figure (b). Black contour lines are ranged from 0 to 2 for figure (c)

In Figs. 3b, 4b and 5b, the hatched region represents values for which the normalized numerical group velocity is negative. They are denoted as the q-waves region. The dark red region represents positive normalized numerical group velocity for which its value is much greater than 1. The light red region represents positive normalized numerical group velocity for which its value is close to 1. For the Fourier spectral scheme (Fig. 6), it can be seen that the numerical amplification factor, normalized numerical group and phase velocity remains $> rsim 0.99$ for the entire spectrum of reduced wavenumbers up to a CFL number of approximately 0.4. The non red contour region in Fig. 6a represents the stable region where $|Z| < 1$. Within the stable region, there exists a critical CFL number, R for which the numerical amplification becomes significantly less than 1 beyond the critical value, leading to dissipation error. In this context, the critical CFL number, $R_{crit}$ is defined as the value for which the amplification factor remains close to 1.0 (or $|Z| > rsim 0.99$) for the full range of reduced wavenumber considered. For example, $R_{crit}$ is found to be approximately 0.4 for the case representing the FS spatial and RK4 temporal discretization scheme in Fig. 6. At $R \lesssim$ 0.4, the numerical solution of the 1D advection equation exhibit minimal numerical dissipation since $|Z| > rsim 0.99$ for all wavenumbers considered. For $0.4 \lesssim R \lesssim 0.8$, the numerical solution will exhibit numerical dissipation despite being stable. It should be noted that the dissipative error originates from the temporal discretization scheme since the Fourier spectral scheme is non dissipative. A comparison of $R_{crit}$ for the different spatial discretization schemes tested under full discretization (Figs. 3a, 4a, 5a and 6a) consideration shows that $R_{crit}$ becomes lower when a higher accuracy spatial discretization scheme is used. This suggests that the use of higher accuracy spatial discretization schemes comes at a cost of longer time integration period.

3.2 2D Advection

For the purpose of demonstrating the 1D full discretization analysis to two-dimensional space, the 2D advection equation is considered. In order to compare the case between a hybrid FD/FS method and a full FD method in two spatial directions, the CDS4/FS and CDS4/CDS4 cases are used in this subsection. The fully discretized numerical dispersion relation of the 2D advection equation is given as $\omega ^{*}_{full} = U^{*}_{phase, full}k_{x} + V^{*}_{phase, full}k_{y}$ [23]. In principle, group and phase velocities computed under full discretization consideration are determined from the fully discretized numerical dispersion relation, $\omega ^{*}_{full}$. Since $\omega ^{*}_{full}$ is related to $\psi$ and consequently Z, the calculation procedure of group and phase velocities under full discretization starts first from determining the semi-discretized numerical dispersion relation of the governing equation. The semi-discretized numerical dispersion relation of the 2D advection equation can be shown to be $\omega ^{*}_{semi} = {\overline{u}}k^{*}_{x} + {\overline{v}}k^{*}_{y}$ using similar steps as outlined for the 1D advection equation. By substituting the semi-discretized numerical dispersion relation into the numerical amplification factor in Eq. (10), one can obtain the group and phase velocities for full discretization with the phase shift, $\psi$. In this instance, the CFL number used in the x and y directions are the same ($R_{x}$ = $R_{y} = R$).

$$\begin{aligned}&\begin{aligned} \frac{U_{grp, full}^{*}}{U_{grp, exact}}&= \frac{1}{{\overline{u}}} \frac{\partial \omega ^{*}_{full} }{\partial k_{x}} = \frac{1}{{\overline{u}} \varDelta t} \frac{\partial \psi }{\partial k_{x}} = \frac{1}{R \varDelta x} \frac{\partial \psi }{\partial k_{x}},\\ \frac{V_{grp, full}^{*}}{V_{grp, exact}}&= \frac{1}{{\overline{v}}} \frac{\partial \omega ^{*}_{full} }{\partial k_{y}} = \frac{1}{{\overline{v}} \varDelta t} \frac{\partial \psi }{\partial k_{x}} = \frac{1}{R \varDelta y} \frac{\partial \psi }{\partial k_{y}}, \end{aligned} \end{aligned}$$

(19)

$$\begin{aligned}&\begin{aligned} \frac{U_{phase, full}^{*}}{U_{phase, exact}}&= \frac{1}{U_{phase,exact}} \frac{\psi }{k_{x} \varDelta t},\\ \frac{V_{phase, full}^{*}}{V_{phase, exact}}&= \frac{1}{V_{phase,exact}} \frac{\psi }{k_{y} \varDelta t}, \end{aligned} \end{aligned}$$

(20)

where $U_{phase, exact} = \frac{{\overline{u}}k_{x} + {\overline{v}}k_{y}}{k_{x}}$ and $V_{phase, exact} = \frac{{\overline{u}}k_{x} + {\overline{v}}k_{y}}{k_{y}}$. In order to neglect the effect of the temporal discretization error, R is kept at a sufficiently low value of 0.1. For illustration sake, this can be seen in Fig. 3b, where the normalised group velocity contour plot of the RK4-CDS4 case based on the 1D advection remain approximately the same even when the CFL number is varied from 0 up to approximately 0.4. For sufficiently low CFL number, the magnitude of the numerical amplification factor, |Z| would remain very close to 1.0 for all combinations of reduced wavenumber. Also, the dispersion error present in the group velocity plots would come solely from the spatial discretization scheme under such consideration.

Fig. 7

Normalized group velocities a $\frac{U^{*}_{grp, full}}{U_{grp, exact}}$ and b $\frac{V^{*}_{grp, full}}{V_{grp, exact}}$ plots of 2D advection equation using RK4 temporal scheme and CDS4 spatial discretization scheme applied in both directions for $R_{x} = R_{y} = 0.1$. Hatched region represents the q-wave region ($\frac{U^{*}_{grp, full}}{U_{grp, exact}}, \frac{V^{*}_{grp, full}}{V_{grp, exact}} < 0$)

Fig. 8

Normalized phase velocities a $\frac{U^{*}_{phase, full}}{U_{phase, exact}}$ and b $\frac{V^{*}_{phase, full}}{V_{phase, exact}}$ plots of 2D advection equation using RK4 temporal scheme and CDS4 spatial discretization scheme applied in both directions for $R_{x} = R_{y} = 0.1$

Fig. 9

Normalized group velocities a $\frac{U^{*}_{grp, full}}{U_{grp, exact}}$ and b $\frac{V^{*}_{grp, full}}{V_{grp, exact}}$ plots of 2D advection equation using RK4 temporal scheme and CDS4/FS spatial discretization scheme (CDS4 in $\textit{x}$ and FS in $\textit{y}$) for $R_{x} = R_{y} = 0.1$. Hatched region represents the q-wave region ($\frac{U^{*}_{grp, full}}{U_{grp, exact}}, \frac{V^{*}_{grp, full}}{V_{grp, exact}} < 0$)

Fig. 10

Normalized phase velocities a $\frac{U^{*}_{phase, full}}{U_{phase, exact}}$ and b $\frac{V^{*}_{phase, full}}{V_{phase, exact}}$ plots of 2D advection equation using RK4 temporal scheme and CDS4/FS spatial discretization scheme (CDS4 in $\textit{x}$ and FS in $\textit{y}$) for $R_{x} = R_{y} = 0.1$

In Fig. 7, it can be seen that the normalised group velocities remain close to 1.0 ($> rsim 0.99$) in either spatial direction for $k_{x} \varDelta x \lesssim \pi /4$ or $k_{y} \varDelta y \lesssim \pi /4$. Similar to the 1D case, the group velocity decreases and becomes negative with increasing reduced wavenumber in the x or y direction respectively. It can be seen that the group velocities, $U^{*}_{grp}$ and $V^{*}_{grp}$ display variations only in the $k_{x}$ or $k_{y}$ direction separately. The implication of this is that the resolution characteristic of the numerical scheme in x or y spatial direction is only dependent on the type of scheme used in the direction considered. As such, the numerical schemes used for the 2D advection equation can be analysed in 1D wavenumber space with no loss of information. In order to understand this phenomenon, the results are compared to group velocities calculated based on semi-discretized numerical dispersion relation since the group velocities calculated based on fully discretized numerical dispersion relation would approach those calculated based on the semi-discretized numerical dispersion relation as the time step approaches zero ($U_{grp, full}^{*}$ $\rightarrow$ $U_{grp, semi}^{*}$ as $\varDelta t$ $\rightarrow$ 0). The normalized group velocities calculated based on the semi-discretized numerical dispersion relation are given as follows:

$$\begin{aligned} \begin{aligned} \frac{U_{grp, semi}^{*}}{U_{grp, exact}}&= \frac{1}{{\overline{u}}} \frac{\partial \omega ^{*}_{semi} }{\partial k_{x}} = \frac{\partial k^{*}_{x}}{\partial k_{x}},\\ \frac{V_{grp, semi}^{*}}{V_{grp, exact}}&= \frac{1}{{\overline{v}}} \frac{\partial \omega ^{*}_{semi} }{\partial k_{y}} = \frac{\partial k^{*}_{y}}{\partial k_{y}}. \end{aligned} \end{aligned}$$

(21)

In Eq. (21), it can be seen that the normalized group velocity in one spatial direction is only dependent on the gradient of the modified wavenumber in one specific spatial direction. Since the modified wavenumber of the CDS4 scheme is a function of wavenumber in one spatial direction, the normalized group velocity will exhibit variation in only one spatial direction. In Fig. 9, the case representing the RK4 temporal and CDS4/FS spatial discretization scheme is plotted. It can be seen that the normalised group velocity in the x direction remains the same as the CDS4/CDS4 case while the normalised group velocity in the y direction become 1.0 for all reduced wavenumber combinations (since $\frac{\partial k^{*}_{y}}{\partial k_{y}} = 1$). The results of the hybrid discretization case illustrated here shows that the hybrid spectral/finite difference case will always produce results that are more accurate than the full finite difference case for the two dimensional advection equation since there is no dispersive error for $\frac{V^{*}_{grp, full}}{V_{grp, exact}}$ for all wavenumber combination. The dispersive error from the spatial scheme used in the x direction does not influence $\frac{V^{*}_{grp, full}}{V_{grp, exact}}$. The numerical results of the 2D advection equation that correspond to this theoretical assessment can be found in [29]. The normalized numerical phase speed plots are given in Figs. 8 and 10. For the CDS4/CDS4 case, the normalized phase speed ($\frac{U^{*}_{phase, full}}{U_{phase, exact}}$, $\frac{V^{*}_{phase, full}}{V_{phase, exact}}$) plots showcase variation in both spatial directions. For the CDS4/FS case, the normalized phase speed ($\frac{U^{*}_{phase, full}}{U_{phase, exact}}$, $\frac{V^{*}_{phase, full}}{V_{phase, exact}}$) plots showcase variation mainly in one spatial directions. The variation of group velocities in only one spatial direction is also valid for the DRP4 and CCOM6 since the modified wavenumber of the DRP4 and CCOM6 schemes are only a function of the wavenumber in one spatial direction. However, this is not necessarily the case for the $\hbox {CDS4}_{{opt}}$ scheme. The main reason being that its modified wavenumber is a function of the wavenumber in both spatial direction. Hence its group velocities ($V^{*}_{grp, full}$, $U^{*}_{grp, full}$) will showcase variation in both spatial directions for both the $\hbox {CDS4}_{{opt}}$/$\hbox {CDS4}_{{opt}}$ and $\hbox {CDS4}_{{opt}}$/FS cases at low R when temporal discretization effects can be neglected. Although it should be noted that the results of the hybrid $\hbox {CDS4}_{{opt}}$/FS case is still more accurate than the $\hbox {CDS4}_{{opt}}$/$\hbox {CDS4}_{{opt}}$ case since there is no dispersive error from the spatial scheme in the y direction. The plot of group and phase velocities of both the $\hbox {CDS4}_{{opt}}$/$\hbox {CDS4}_{{opt}}$ and $\hbox {CDS4}_{{opt}}$/FS cases are shown in Figures 11, 1213 and 14.

Fig. 11

Normalized group velocities a $\frac{U^{*}_{grp, full}}{U_{grp, exact}}$ and b $\frac{V^{*}_{grp, full}}{V_{grp, exact}}$ plots of 2D advection equation using RK4 temporal scheme and $\hbox {CDS4}_{{opt}}$/$\hbox {CDS4}_{{opt}}$ spatial discretization scheme for $R_{x} = R_{y} = 0.1$. Hatched region represents the q-wave region ($\frac{U^{*}_{grp, full}}{U_{grp, exact}}, \frac{V^{*}_{grp, full}}{V_{grp, exact}} < 0$)

Fig. 12

Normalized phase velocities a $\frac{U^{*}_{phase, full}}{U_{phase, exact}}$ and b $\frac{V^{*}_{phase, full}}{V_{phase, exact}}$ plots of 2D advection equation using RK4 temporal scheme and $\hbox {CDS4}_{{opt}}$/$\hbox {CDS4}_{{opt}}$ spatial discretization scheme for $R_{x} = R_{y} = 0.1$

Fig. 13

Normalized group velocities a $\frac{U^{*}_{grp, full}}{U_{grp, exact}}$ and b $\frac{V^{*}_{grp, full}}{V_{grp, exact}}$ plots of 2D advection equation using RK4 temporal scheme and $\hbox {CDS4}_{{opt}}$/FS spatial discretization scheme for $R_{x} = R_{y} = 0.1$

Fig. 14

Normalized phase velocities a $\frac{U^{*}_{phase, full}}{U_{phase, exact}}$ and b $\frac{V^{*}_{phase, full}}{V_{phase, exact}}$ plots of 2D advection equation using RK4 temporal scheme and $\hbox {CDS4}_{{opt}}$/FS spatial discretization scheme for $R_{x} = R_{y} = 0.1$

4 Fourier Analysis of the Linearized Compressible Navier–Stokes Under Full Discretization

Following a full discretization analysis of the numerical schemes in two dimensional wavenumber space based on the advection equation, the full discretization analysis framework is now extended to the linearized compressible Navier–Stokes Equations (LCNSE). Similar procedures as established earlier will be first used to derive the physical dispersion relation of the linearized compressible Navier–Stokes. Following that, one can obtain the semi-discretized numerical dispersion relation which is used within the numerical amplification factor. The numerical group and phase velocities under full discretization can then be obtained using the fully discretized numerical dispersion relation which is a function of the phase shift of the numerical amplification factor.

4.1 Dispersion Relation of the LCNSE

In order to derive the semi-discretized numerical dispersion relation ($\omega ^{*}_{semi}$) of the LCNSE, the derivation of the physical dispersion relation of the LCNSE is first introduced. The primitive form of the two dimensional compressible Navier–Stokes Equations is given by [1]:

$$\begin{aligned} \frac{\partial \mathbf{U} }{\partial t} + A_{1}\frac{\partial \mathbf{U} }{\partial x} + A_{2}\frac{\partial \mathbf{U} }{\partial y} - B_{1} \frac{\partial ^{2} \mathbf{U} }{\partial x^{2}} -B_{2} \frac{\partial ^{2} \mathbf{U} }{\partial y^{2}} - B_{3} \frac{\partial ^{2} \mathbf{U} }{\partial x \partial y} = \mathbf{0} , \end{aligned}$$

(22)

where

$$\begin{aligned} \begin{aligned}&\mathbf{U} = \begin{bmatrix} \rho \\ u \\ v \\ p \end{bmatrix} \quad \mathbf{A} _{1} = \begin{bmatrix} u &{} \rho &{} 0 &{} 0 \\ 0 &{} u &{} 0 &{} \frac{1}{\rho } \\ 0 &{} 0 &{} u &{} 0 \\ 0 &{} \rho a^{2} &{} 0 &{} u \end{bmatrix} \quad \mathbf{A} _{2} = \begin{bmatrix} v &{} 0 &{} \rho &{} 0 \\ 0 &{} v &{} 0 &{} 0 \\ 0 &{} 0 &{} v &{} \frac{1}{\rho }\\ 0 &{} 0 &{} \rho a^{2} &{} v \end{bmatrix}\\\\&\mathbf{B} _{1} = \begin{bmatrix} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} \frac{4}{3} \frac{\mu }{\rho } &{} 0 &{} 0 \\ 0 &{} 0 &{} \frac{\mu }{\rho } &{} 0 \\ -\frac{\gamma \mu p}{P_{r} \rho ^{2}} &{} 0 &{} 0 &{} \frac{\gamma \mu }{P_{r} \rho } \end{bmatrix} \quad \mathbf{B} _{2} = \begin{bmatrix} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} \frac{\mu }{\rho } &{} 0 &{} 0 \\ 0 &{} 0 &{} \frac{4}{3} \frac{\mu }{\rho } &{} 0 \\ -\frac{\gamma \mu p}{P_{r} \rho ^{2}} &{} 0 &{} 0 &{} \frac{\gamma \mu }{P_{r} \rho } \end{bmatrix} \quad \mathbf{B} _{3} = \begin{bmatrix} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{}\frac{1}{3} \frac{\mu }{\rho } &{} 0 \\ 0 &{} \frac{1}{3} \frac{\mu }{\rho } &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{bmatrix}. \end{aligned} \end{aligned}$$

(23)

$\rho$, p, a, u and v represent the density, pressure, speed of sound as well as the velocity components in the x and y directions, respectively. Firstly, linearization of the two dimensional Compressible Navier Stokes Equations is carried out by applying a perturbation to the mean flow, ${\mathbf {U}} = \mathbf {U}_{\mathbf {0}} + \mathbf {U}^{\mathbf {'}}$ and assuming that the perturbation, $\mathbf {U}^{\mathbf {'}}$ is sufficiently small. This assumption allows one to neglect all non linear terms associated with either the squared of perturbed quantities or products of the perturbed quantities and perturbed gradients from the resulting expanded equations. Secondly, it is assumed that the mean flow is homogeneous, which allows one to neglect all mean spatial gradients. Thirdly, the mean flow is considered to be weakly compressible, which means that an isothermal assumption (no temperature gradients) is used on the energy equation. Finally, the fluid medium is considered to be mono-atomic and hence bulk viscosity is taken to be zero. After subtracting the mean flow terms from the expanded equations, the two dimensional linearized compressible Navier Stokes Equations can be written as follows [14]:

$$\begin{aligned} \frac{\partial \mathbf {U}^{\mathbf {'}}}{\partial t} + \overline{A_{1}}\frac{\partial \mathbf {U}^{\mathbf {'}}}{\partial x} + \overline{A_{2}}\frac{\partial \mathbf {U}^{\mathbf {'}}}{\partial y} - \overline{B_{1}} \frac{\partial ^{2} \mathbf {U}^{\mathbf {'}}}{\partial x^{2}} - \overline{B_{2}} \frac{\partial ^{2} \mathbf {U}^{\mathbf {'}}}{\partial y^{2}} - \overline{B_{3}} \frac{\partial ^{2} \mathbf {U}^{\mathbf {'}}}{\partial x \partial y} = \mathbf{0} , \end{aligned}$$

(24)

where

$$\begin{aligned} \begin{aligned}&\mathbf {U}^{\mathbf {'}} = \begin{bmatrix} \rho ^{'} \\ u^{'} \\ v^{'} \\ p^{'} \end{bmatrix} \quad {\overline{\mathbf {A}}}_{\mathbf {1}} = \begin{bmatrix} {\overline{u}} &{} {\overline{\rho }} &{} 0 &{} 0 \\ 0 &{} {\overline{u}} &{} 0 &{} \frac{1}{{\overline{\rho }}} \\ 0 &{} 0 &{} {\overline{u}} &{} 0 \\ 0 &{} {\overline{\rho }} a^{2} &{} 0 &{} {\overline{u}} \end{bmatrix} \quad \overline{\mathbf {A}_{\mathbf {2}}} = \begin{bmatrix} {\overline{v}} &{} 0 &{} {\overline{\rho }} &{} 0 \\ 0 &{} {\overline{v}} &{} 0 &{} 0 \\ 0 &{} 0 &{} {\overline{v}} &{} \frac{1}{\rho }\\ 0 &{} 0 &{} {\overline{\rho }} a^{2} &{} {\overline{v}} \end{bmatrix}\\\\&\overline{\mathbf {B}_{\mathbf {1}}} = {\overline{\nu }} \begin{bmatrix} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} \frac{4}{3} &{} 0 &{} 0 \\ 0 &{} 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{bmatrix} \quad \overline{\mathbf {B}_{\mathbf {2}}} = {\overline{\nu }} \begin{bmatrix} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} \frac{4}{3} &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{bmatrix} \quad \overline{\mathbf {B}_{\mathbf {3}}} = {\overline{\nu }} \begin{bmatrix} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{}\frac{1}{3} &{} 0 \\ 0 &{} \frac{1}{3} &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{bmatrix}. \end{aligned} \end{aligned}$$

(25)

In Eq. (25), ${\overline{\nu }}$ = $\frac{ {\overline{\mu }} }{ {\overline{\rho }} }$. Field variables denoted with $(..)^{'}$ represent the perturbed field while those with $\overline{(..)}$ represent the time averaged mean field. A Fourier-Laplace transform of $\mathbf{U} ^{\prime }$ is given by:

$$\begin{aligned} \widehat{\mathbf{U }}^{\prime }(k_{x},k_{y}, \omega ) = \frac{1}{(2 \pi )^{3}} \int _{0}^{\infty } \int _{-\infty }^{\infty } \int _{-\infty }^{\infty } \mathbf{U} ^{\prime } e^{i(k_{x}x + k_{y}y -\omega t)} dx dy dt, \end{aligned}$$

(26)

where the components of the wavenumber and frequency are denoted as $k_{x}$, $k_{y}$ and $\omega$ respectively. By substituting Eq. (26) into Eq. (24):

$$\begin{aligned} \begin{aligned} \frac{\partial \widehat{\mathbf{U }}^{\prime }}{\partial t}&= -i \omega \widehat{\mathbf{U }}^{\prime }, \quad \frac{\partial \widehat{\mathbf{U }}^{\prime }}{\partial x} = i k_{x} \widehat{\mathbf{U }}^{\prime }, \quad \frac{\partial \widehat{\mathbf{U }}^{\prime }}{\partial y} = i k_{y} \widehat{\mathbf{U }}^{\prime }, \\ \frac{\partial ^{2} \widehat{\mathbf{U }}^{\prime }}{\partial x^{2}}&= -(k_{x})^{2} \widehat{\mathbf{U }}^{\prime }, \quad \frac{\partial ^{2} \widehat{\mathbf{U }}^{\prime }}{\partial y^{2}} = -(k_{y})^{2} \widehat{\mathbf{U }}^{\prime }, \quad \frac{\partial ^{2} \widehat{\mathbf{U }}^{\prime }}{\partial x \partial y} = - k_{x}k_{y} \widehat{\mathbf{U }}^{\prime }. \end{aligned} \end{aligned}$$

(27)

Equation (24) can be re-written as:

$$\begin{aligned} \mathbf{I} \omega \hat{\mathbf{U }} = \underbrace{(k_{x} \overline{A_{1}} + k_{y} \overline{A_{2}} - i(k_{x})^{2} \overline{B_{1}} - i(k_{y})^{2} \overline{B_{2}} - i(k_{x}k_{y}) \overline{B_{3}} )}_\mathbf{C } \hat{\mathbf{U }}, \end{aligned}$$

(28)

where $\mathbf{I}$ represents the identity matrix and

$$\begin{aligned} \mathbf{C} = \left[ \begin{matrix} k_{x}{\overline{u}}+k_{y}{\overline{v}} &{} k_{x}{\overline{\rho }} &{} k_{y}{\overline{\rho }} &{} 0 \\ 0 &{} k_{x}{\overline{u}}+k_{y}{\overline{v}}-\frac{4}{3}i {\overline{\nu }} (k_{x})^{2} - i {\overline{\nu }} (k_{y})^{2} &{} -\frac{1}{3} i {\overline{\nu }} k_{x} k_{y} &{} \frac{k_{x}}{{\overline{\rho }}}\\ 0 &{} -\frac{1}{3} i {\overline{\nu }} k_{x} k_{y} &{} k_{1}{\overline{u}}+k_{y}{\overline{v}} - i{\overline{\nu }} (k_{x})^{2} - \frac{4}{3} i{\overline{\nu }} (k_{y})^{2} &{} \frac{k_{y}}{{\overline{\rho }}} \\ 0 &{} {\overline{\rho }} {\overline{a}}^{2}k_{x} &{} {\overline{\rho }} {\overline{a}}^{2}k_{y} &{} k_{x}{\overline{u}}+k_{y}{\overline{v}} \end{matrix} \right] . \end{aligned}$$

(29)

Solving for the four eigenvalues of $\mathbf{C}$ leads to the following solution:

$$\begin{aligned} \begin{aligned} \omega _{1}&= k_{x}{\overline{u}} + k_{y}{\overline{v}} - i \nu |\mathbf{k} |^{2} + |\mathbf{k} | a, \\ \omega _{2}&= k_{x}{\overline{u}} + k_{y}{\overline{v}} - i \frac{4}{3} \nu |\mathbf{k} |^{2} - |\mathbf{k} | a, \\ \omega _{3}&= \omega _{4} = k_{x}{\overline{u}} + k_{y}{\overline{v}}, \end{aligned} \end{aligned}$$

(30)

where $|{\mathbf {k}}| = \sqrt{ (k_{x})^{2} + (k_{y})^{2} }$. The first two eigenvalues represent the dispersion relation for the acoustic waves while the last two eigenvalues represent the dispersion relation for the entropy and vortical waves. The velocity components can be expressed in terms of the product of the speed of sound and Mach number of respective spatial directions. The resulting dispersion relation of the different modes can be re-written as:

$$\begin{aligned} \begin{aligned} \omega _{1}&= k_{x}M_{x}a + k_{y}M_{y}a - i \nu |\mathbf{k} |^{2} + |\mathbf{k} | a, \\ \omega _{2}&= k_{x}M_{x}a + k_{y}M_{y}a - i \frac{4}{3} \nu |\mathbf{k} |^{2} - |\mathbf{k} | a, \\ \omega _{3}&= \omega _{4} = k_{x}M_{x}a + k_{y}M_{y}a. \end{aligned} \end{aligned}$$

(31)

For the dispersion relation belonging to the acoustic modes ($\omega _{1,2}$), the first two term on the right hand side represent the advection contribution. The third and fourth term represent the diffusive and acoustic contribution respectively. For the dispersion relation of the entropy and vortical modes ($\omega _{3,4}$), it only comprises of the advection contribution.

4.2 Physical Group and Phase Velocities of the LCNSE

For the purpose of having a reference solution for comparison and normalization sake, the physical group and phase velocity are computed. The calculation of the physical group and phase velocities for one of the acoustic modes ($\omega _{1}$) are demonstrated below:

$$\begin{aligned}&\begin{aligned} U_{grp, exact}&= \frac{\partial \omega _{1}}{\partial k_{x}} = M_{x}a - i2 \nu k_{x} + \frac{ak_{x}}{\sqrt{k^{2}_{x} + k^{2}_{y}}},\\ V_{grp, exact}&= \frac{\partial \omega _{1}}{\partial k_{y}} = M_{y}a - i2 \nu k_{y} + \frac{ak_{y}}{\sqrt{k^{2}_{x} + k^{2}_{y}}}, \end{aligned} \end{aligned}$$

(32)

$$\begin{aligned}&\begin{aligned} U_{phase, exact}&= \frac{\omega _{1}}{k_{x}} = M_{x}a + \frac{k_{y}M_{y}a}{k_{x}} - \frac{i\nu |\mathbf{k} |^{2}}{k_{x}} + \frac{|\mathbf{k} |a}{k_{x}},\\ V_{phase, exact}&= \frac{\omega _{1}}{k_{y}} = \frac{k_{x}M_{x}a}{k_{y}} + M_{y}a - \frac{i\nu |\mathbf{k} |^{2}}{k_{y}} + \frac{|\mathbf{k} |a}{k_{y}}. \end{aligned} \end{aligned}$$

(33)

4.3 Numerical Group and Phase Velocities of the LCNSE

In this section, numerical group and phase velocities computed under full discretization consideration will be presented. As before, the results presented here are demonstrated for only one of the acoustic modes ($\omega ^{*}_{semi,1}$). The procedure is applicable to the other modes in Eq. (31). In order to obtain the group velocities under full discretization, the numerical amplification factor is introduced. The semi-discretized numerical dispersion relation ($\omega ^{*}_{semi,1}$) of the LCNSE is substituted into the numerical amplification factor of the RK4 scheme, Eq. (10). The semi-discretized dispersion relation is obtained by replacing the actual wavenumber with the modified wavenumber, $k^{*}_{x,y}$ = $k_{x,y}$) to account for the spatial discretization effects.

$$\begin{aligned} \varDelta t \omega _{semi,1}^{*} = \varDelta t (k^{*}_{x}M_{x}a + k^{*}_{y}M_{y}a - i \nu |\mathbf{k} ^{*}|^{2} + |\mathbf{k} ^{*}| a) \end{aligned}$$

(34)

and defining the grid spacing, CFL number and cell Reynolds number in terms of the speed of sound:

$$\begin{aligned}&R_{a} = \frac{a \varDelta t}{\varDelta x}, \quad Re_{a} = \frac{a \varDelta }{\nu }, \quad \varDelta = \varDelta x = \varDelta y, \end{aligned}$$

(35)

$$\begin{aligned}&\begin{aligned} \varDelta t \omega _{semi,1}^{*}&= R_{a} \bigg ( M_{x} k^{*}_{x} \varDelta + M_{y} k^{*}_{y} \varDelta - i \frac{\nu }{a \varDelta } |\mathbf{k} ^{*} \varDelta |^{2} + |\mathbf{k} ^{*} \varDelta | \bigg )\\&= R_{a} \bigg ( M_{x} k^{*}_{x} \varDelta + M_{y} k^{*}_{y} \varDelta - i \frac{1}{Re_{a}} |\mathbf{k} ^{*} \varDelta |^{2} + |\mathbf{k} ^{*} \varDelta | \bigg ). \end{aligned} \end{aligned}$$

(36)

Numerical group and phase velocities under full discretization are computed from the fully discretized numerical dispersion relation ($\omega ^{*}_{full}$) as follows:

$$\begin{aligned}&U_{grp, full}^{*} = \frac{\partial \omega ^{*}_{full}}{\partial k_{x}} = \frac{1}{\varDelta t} \frac{\partial \psi }{\partial k_{x}}, \quad V_{grp, full}^{*} = \frac{\partial \omega ^{*}_{full}}{\partial k_{y}} = \frac{1}{\varDelta t} \frac{\partial \psi }{\partial k_{y}}, \end{aligned}$$

(37)

$$\begin{aligned}&\begin{aligned}&U_{phase, full}^{*} = \frac{\omega ^{*}_{full}}{k_{x}} = \frac{\psi }{k_{x} \varDelta t}, \quad V_{phase, full}^{*} = \frac{\omega ^{*}_{full}}{k_{y}} = \frac{\psi }{k_{y} \varDelta t},\\&\psi = -tan^{-1} \bigg ( \frac{Im(Z(\varDelta t \omega _{semi,1}^{*}) )}{Real( Z(\varDelta t \omega _{semi,1}^{*}) )} \bigg ). \end{aligned} \end{aligned}$$

(38)

An inspection of Eq. (36) shows that the numerical amplification factor, group and phase velocities under full discretization are dependent on the following factors:

• CFL number ($R_{a}$); as $R_{a}$ increases, the magnitude of the amplification factor as well as the region of stability in the two dimensional wavenumber plane changes. Beyond a certain value of $R_{a}$, dispersion error may be introduced as well, reflecting changes in the group and phase velocities.

• Mach number; a change in the Mach number regime not only affects the amplitude of the amplification factor but also the characteristics of the anisotropy map at different Mach number. At low Mach number, the acoustic and viscous terms dominate over the convective term. On the other hand, the convective term dominates over the viscous and acoustic terms at high Mach number.

• Viscous effects are related to the cell Reynolds number, $Re_{a}$. As $Re_{a}$ reduces, viscous effects become more dominant.

5 Analysis of Linearized Compressible Navier–Stokes in Two Dimensional Wavenumber Plane

In this section, analytical solutions of numerical amplification factor, group and phase velocities for full discretization are plotted and compared for different combinations of spatial discretization schemes in Table 1. The RK4 temporal scheme is used and a sufficiently low CFL number, $R_{a}$ of 0.1 was used so as to neglect any error from the temporal discretization scheme. Numerical group and phase velocities under full discretization are computed based on Eqs. (37) and (38). The parameters and type of schemes used for each case are listed in Table 1. The speed of sound, a is taken to be 1.0 for all cases. In order to demonstrate the coupling effect in the group velocities, $M_{y}$ is kept zero to illustrate the contribution of the acoustic term solely for $V^{*}_{grp}$. This is used especially in the context of a comparison between the hybrid and full FD discretization case. $M_{x}$ is set to 0.5 for all cases illustrated in this paper. It should be pointed out that the coupling effect of the numerical schemes in two spatial directions is independent of the change in Mach number. A change in Mach number only influences the magnitude of the group velocities. This is evident from the dispersion relation given in Eq. 36.

Table 1 Table representing parameter values and numerical schemes for various cases

5.1 Effect of Varying Spatial Discretization Schemes for the Convective Term

In order to assess the effect of the spatial discretization on numerical anisotropy of the convection term, a comparison of the normalized group velocities using different combinations of spatial discretisation schemes was carried out. The parameters used for the plots in Fig. 15a–h correspond to cases 1a–1h of Table 1. As before, the CFL number is kept sufficiently small to minimise the effect of the temporal discretization error. Normalized numerical group velocity plots belonging to the acoustic mode ($\omega ^{*}_{1}$) for different combinations of spatial discretization methods are presented in Fig. 15. Figures 15a, c, e and g represent the hybrid discretization case. Figures 15b, d, f, h represent the full finite difference cases. As the normalized numerical group velocity, $\frac{V^{*}_{grp, full}}{V_{grp, exact}}$ encompasses solely the effect of the acoustic term, $\frac{V^{*}_{grp, full}}{V_{grp, exact}}$ of the hybrid discretization cases show how the dispersive error originating from the spatial discretization scheme in the x direction directly influences the group velocity in the y direction, especially in the region where $k_{x} \varDelta x$ is greater than 0.6$\pi$ across the full range of $k_{y} \varDelta y$. For regions where $k_{x} \varDelta x$ is much less than 0.5$\pi$, $\frac{V^{*}_{grp, full}}{V_{grp, exact}}$ of the hybrid discretization cases have values close to 1.0. This implies that the dispersive error from the spatial scheme in the x direction has little effect on $\frac{V^{*}_{grp, full}}{V_{grp, exact}}$ for low values of $k_{x} \varDelta x$. As for $\frac{U^{*}_{grp, full}}{U_{grp, exact}}$, the dispersive error from the FD scheme used in the x direction in clearly reflected in all cases. For the CDS4/CDS4 and CDS4/FS cases, it can be clearly seen that $\frac{U^{*}_{grp, full}}{U_{grp, exact}}$ becomes negative when $k_{x} \varDelta x > rsim 0.625 \pi$. This represents the q-wave region. Differences in $\frac{U^{*}_{grp, full}}{U_{grp, exact}}$ between the hybrid and full FD discretization cases can be found in regions of high $k_{y} \varDelta y$ values ($k_{y} \varDelta y > rsim 0.8 \pi$). This is evidence of dispersive errors propagating from the spatial scheme used in the y direction to the x direction. For $k_{y} \varDelta y \lesssim 0.4 \pi$, $\frac{U^{*}_{grp, full}}{U_{grp, exact}}$ remains almost the same for both the CDS4/FS and CDS4/CDS4 discretization cases. The main difference between the different set of spatial schemes lies in the critical wavenumber at which the normalized group velocities become negative. Normalized numerical phase velocity plots belonging to the acoustic mode ($\omega ^{*}_{1}$) for different combinations of spatial discretization methods are presented in Fig. 16. The key difference between the hybrid and full finite difference discretization cases lies in the variation of the phase velocity. The hybrid discretization cases display variation in phase velocity mainly in one spatial direction. The full finite difference cases display variation in phase velocity in two spatial directions.

In general, depending on the wave propagation angle ($\theta _{w}$) and the reduced wavenumber ($k \varDelta$), the normalized numerical group velocity $\frac{V^{*}_{grp, full}}{V_{grp, exact}}$ of the hybrid discretization case may be more or less accurate than the full FD discretization case as compared to the exact solution. To illustrate this more clearly, Fig. 17 compares the solution of the numerical group velocity $V^{*}_{grp, full}$ as a function of $\theta _{w}$ for different spatial discretization methods at a fixed $k \varDelta$. The group velocities are extracted from the two dimensional wavenumber plot based on the schematic shown in Fig. 22. The numerical results are compared against the exact group velocity solution, $V_{grp, exact}$. The plot in Fig. 17a can be demarcated into three regions for the hybrid schemes. Below a certain value of $\theta _{w}$, $V^{*}_{grp, full}$ of the hybrid and full finite difference cases give the same results. At intermediate values of $\theta _{w}$, $V^{*}_{grp, full}$ of the full finite difference case come closer to the exact solution. For large values of $\theta _{w}$, $V^{*}_{grp, full}$ of the hybrid discretization case come closer to the exact solution. The exact values of $\theta _{w}$ separating these three region is scheme dependent and $k \varDelta$ dependent. A comparison between the hybrid cases in Fig. 17a shows that $V^{*}_{grp, full}$ of the CCOM6/FS case has a larger range of $\theta _{w}$ for which its result is the same as its respective full FD case. Also, $V^{*}_{grp, full}$ of the CCOM6/FS case is closer to the exact solution than the other hybrid discretization cases for the range of wave propagation angle where $V^{*}_{grp, full}$ is not the same as the exact solution. At lower values of $k \varDelta$ (Fig. 17b, c), these three regions become less distinct. As $k \varDelta \rightarrow 0$, all the different cases would give the same results as the exact solution. To better understand the comparison of the dispersion property between a hybrid and full FD scheme, the integrated sum of the error between the numerical group velocity $V^{*}_{grp,full}$ and the exact solution across all wavenumbers for a range of wave propagation angle is examined. This is illustrated in Fig. 18. In the comparison between a hybrid and full FD discretization schemes, the CDS4, DRP4 and CCOM6 schemes share the same characteristics. At low wave propagation angle, the integrated error of the full FD and hybrid discretization schemes remain the same. At intermediate wave propagation angle, the integrated error of the full FD schemes become smaller than the hybrid scheme. At large wave propagation angle, the integrated error of the full FD schemes diverges while the integrated error of the hybrid discretization schemes converges to zero. This is expected since the integrated error of the hybrid discretization schemes should reduce as we approach the spatial direction where the Fourier spectral method is employed. A slightly different trend is observed for the $\hbox {CDS4}_{{opt}}$ scheme. At low wave propagation angle, the integrated error of the full FD scheme is lower than the hybrid discretization scheme. However, at intermediate and large wave propagation angle, the integrated error of the hybrid scheme remain lower than the full FD scheme and converges toward zero with increasing wave propagation angle.

Fig. 15

Comparison of normalized numerical group velocities ($\frac{U^{*}_{grp, full}}{U_{grp, exact}}$, $\frac{V^{*}_{grp, full}}{V_{grp, exact}}$) under different spatial discretization method for convective terms; a CDS4/FS, b CDS4/CDS4, c DRP4/FS, d DRP4/DRP4, e CCOM6/FS and f CCOM6/CCOM6, g $\hbox {CDS4}_{{opt}}$/FS, h $\hbox {CDS4}_{{opt}}$/$\hbox {CDS4}_{{opt}}$. They represent cases 1a to 1h of Table 1. All colour contour plots are ranged from 0.0 to 8. Hatched region represents the q-wave region where $\frac{U^{*}_{grp, full}}{U_{grp, exact}}, \frac{V^{*}_{grp, full}}{V_{grp, exact}} < 0$

Fig. 16

Comparison of normalized numerical phase velocities ($\frac{U^{*}_{phase, full}}{U_{phase, exact}}$, $\frac{V^{*}_{phase, full}}{V_{phase, exact}}$) under different spatial discretization method for convective terms; a CDS4/FS, b CDS4/CDS4, c DRP4/FS, d DRP4/DRP4, e CCOM6/FS and f CCOM6/CCOM6, g $\hbox {CDS4}_{{opt}}$/FS, h $\hbox {CDS4}_{{opt}}$/$\hbox {CDS4}_{{opt}}$. They represent cases 1a–1h of Table 1. All contour plots are ranged from 0.0 to 1.0

Fig. 17

Comparison of numerical group velocities $V^{*}_{grp, full}$ as a function of $\theta _{w}$ for different spatial discretization method; a $k \varDelta = 0.9\pi$, b $k \varDelta = 0.75\pi$ and c $k \varDelta = 0.6\pi$. They represent cases 1a–1h of Table 1

Fig. 18

Comparison of integrated error of group velocities across all wavenumber between full FD and hybrid discretization for convective terms. They represent cases 1a–1h of Table 1

5.2 Effect of Varying Spatial Discretization Schemes for the Diffusive Term

In order to assess the effect of the spatial discretization on numerical anisotropy of the diffusion term, a comparison of the numerical amplification factor using different combinations of spatial discretization schemes was carried out. Under full discretization considerations, the imaginary viscous term will contribute to both the amplitude of the numerical amplification factor as well as its phase shift, reflecting changes to the dispersion and dissipation property. However, this is not the case in the limit of spatial discretization assumptions. The semi-discretized group velocities contains an imaginary viscous term, from which dispersive error is only introduced when a one-sided biased scheme is used. Since only the real part of the dispersion relation contributes to the group or phase velocities, one sided biased scheme possesses an imaginary component in its modified wavenumber which makes the viscous term real. As such, the use of central type scheme for the viscous terms only affect the amplitude of the numerical amplification factor under semi-discretization assumptions. In Fig. 19a, it can be seen that the numerical amplification factor remains very close to 1.0 for different combinations of reduced wavenumbers when $\nu = 0$ ($Re_{a} = \infty$). With increasing physical viscosity ($\nu$ $\uparrow$), the amplification factor becomes smaller than 1 (Fig. 19b–d). Between values of $Re_{a} = \infty$ and $Re_{a} = 500$, the numerical amplification displays almost no change in magnitude. When $Re_{a}$ becomes smaller than 500, changes in the numerical amplification factor become more significant. It is noted that the physical diffusion reflected in the numerical amplification plot varies according to the spatial discretization scheme used on the viscous term (see Fig. 20). There exist conditions for which the numerical dissipation is either less or more than the physical dissipation. This phenomenon is known as sub/super-dissipation [18]. The case with the FS/FS spatial discretization scheme (Fig. 20a) represents the exact physical viscous dissipation and is only a function of $Re_{a}$. For the CDS4/CDS4, DRP4/DRP4, CCOM6/CCOM6 and $\hbox {CDS4}_{{opt}}$/$\hbox {CDS4}_{{opt}}$ schemes, the maximum dissipation occurs at $\theta _{w} \approx \pi /4$. This is reflected in Fig. 20b–e and Fig. 21a–c, where the minimum |Z| is found to be at $\theta _{w} \approx \pi /4$ for various $k \varDelta$ considered. The region of lowest numerical amplification represent the highest numerical dissipation. Similarly, this result can be interpreted in another manner. Considering a flow with one wave number component with $k_{x} = k_{y}$, the wave propagation angle would instead represent the grid aspect ratio, where $\frac{\varDelta y}{\varDelta x} = \tan (\frac{\pi }{4}) = 1$. This would imply that the maximum numerical dissipation of the viscous term occurs at a grid aspect ratio of 1.0. Also, it can be seen that |Z| of the hybrid discretization cases approaches that of the FS/FS case as the wave propagation angle increases. The exact value of $\theta _{w}$ at which |Z| of the hybrid discretization case become the same as the FS/FS case is scheme dependent, although it is noted that this value is smaller for scheme with better dispersion property.

Fig. 19

Comparison of |Z| with increasing viscosity (or decreasing $Re_{a}$) for CDS4/CDS4 spatial discretization on the viscous term; a $Re_{a}$ = $\infty$, b $Re_{a}$ = 500, c $Re_{a}$ = 200, d $Re_{a}$ = 50. They represent cases 2a–2d of Table 1. All colour contour plots are ranged from 0.989 to 1.0

Fig. 20

Comparison of |Z| using different spatial discretization schemes for the diffusive term at $Re_{a}$ = 200. a FS/FS, b CDS4/CDS4, c DRP4/DRP4, d CCOM6/CCOM6, e $\hbox {CDS4}_{{opt}}$/$\hbox {CDS4}_{{opt}}$, f CDS4/FS, g DRP4/FS, h CCOM6/FS, i $\hbox {CDS4}_{{opt}}$/FS. They represent cases 3a–3i of Table 1. All colour contour plots are ranged from 0.989 to 1.0

Fig. 21

Comparison of |Z| as a function of $\theta _{w}$ for different spatial discretization method; a $k \varDelta = 0.9\pi$, b $k \varDelta = 0.75\pi$ and c $k \varDelta = 0.6\pi$. They represent cases 3a–3i of Table 1

6 Numerical Experiments

The analytical (exact and numerical) group velocity solutions determined through Fourier analysis will be verified in this section. Results from numerical experiments will be used to verify the results obtained from analytical prediction. In these numerical experiments, the boundary conditions were kept periodic. A single frequency modulated Gaussian is used as the initial condition in order to localise the wave packets for easy comparison. A 2D uniform Cartesian grid (300 x 300) is used for most of the cases for grid converged solution. In the case of numerical solutions where the Fourier spectral spatial discretization method is used, a more refined grid of (600 x 600) was used for solutions with large reduced wavenumber initial conditions. The computational domain is [−2$\pi$, 2$\pi$] $\times$ [−2$\pi$, 2$\pi$]. The linearized compressible Navier–Stokes Equations (Eq. 24) are solved numerically and compared with both physical and numerical group velocity solutions. For all numerical experiments, the single frequency modulated Gaussian initial condition was used. The single frequency modulated Gaussian initial condition is mathematically expressed as:

$$\begin{aligned} \phi (x,y,t = 0) = e^{-0.5(x^{2} + y^{2})} \sin (k_{x}x + k_{y}y), \end{aligned}$$

(39)

where

$$\begin{aligned} k_{x} = k \cos (\theta _{w}) \quad k_{y} =k \sin (\theta _{w}). \end{aligned}$$

(40)

$\theta _{w}$ represents the wave propagation angle (or alignment of the wave front with respect to the grid) which defines the wavenumber in both spatial directions. The variable $\phi$ represents either the perturbed pressure or density fields depending on the type of solution sought. For a pressure initial condition, the solution of the acoustic waves is obtained. x and y represent the grid coordinates in Eq. (39). For a given grid spacing, the specified wavenumber is adjusted to match the reduced wavenumber of the initial condition. Analytical solutions of $U^{*}_{grp, full}$, $V^{*}_{grp, full}$, $U_{grp, exact}$, $V_{grp, exact}$ and |Z| can be extracted from the two dimensional wavenumber plot as shown in the previous section for each initial condition. This is illustrated in Fig. 22.

Fig. 22

Extraction of group velocity/numerical amplification factor from 2D wavenumber plot for a given initial condition. Note $k = \sqrt{k_{x}^{2} + k_{y}^{2}}$ and $\varDelta = \varDelta x = \varDelta y$

6.1 Numerical Verification

For verification purpose, the initial conditions used in this sub section are listed in Table 2. Numerical solutions were obtained using the CDS4/CDS4 scheme. The effects of varying $k \varDelta$ and $\theta _{w}$ were compared. In Fig. 23, the shaded patches represent the numerical solution of the CDS4/CDS4 scheme. The black arrows represent the correct physical solution obtained using Fourier spectral spatial discretization and the green arrows represent the analytical group velocity solution corresponding to the spatial discretization schemes used for the numerical solution. The solution based on the Fourier spectral spatial discretization (black arrows) are calculated from the group velocities solution obtained from Eq. (32) corresponding to the initial conditions given in Table 2. The magnitude of the numerical results are normalized such that the values are scaled from 0 to 1. It can be seen that the numerical solution for the acoustic waves (shaded patches) matches the analytical group velocity solution (green arrows). The analytical group velocity solutions are obtained for specific values of $k_{x}$ and $k_{y}$ corresponding to the initial conditions. The product of the group velocities with the final time step gives the position of the wave packets. As seen in Fig. 23a, the wave packets at $t_{max}=3$ are located almost exactly at the same position as the reference solution predicted by the Fourier spectral spatial discretization (black arrows).

Fig. 23

Comparison of physical group velocities (FS/FS) and numerical group velocities (CDS4/CDS4) of wave packets at $t_{max}=3$ for different initial conditions. Numerical solutions are computed using RK4 temporal and CDS4/CDS4 spatial discretization in x and y, $M_{x}$ = 0.5, $M_{y}$ = 0.0, $R_{a}$ = 0.1, $Re_{a}$ = $\infty$. a Case A1, b Case A2, c Case A3, d Case A4. Wave packets are denoted as either $\omega ^{*}_{1}$ or $\omega ^{*}_{2}$ depending on the acoustic mode it represents. In figure d, the wave packets of the different modes of the CDS4/CDS4 case are located in opposite y direction as compared to the FS/FS solution

Table 2 Initial conditions of reduced wavenumbers and wave propagation angles

With increasing $k \varDelta$, the wave packets at $t_{max}=3$ are positioned further away from their exact solution (black arrows) as a result of the dispersive error. This is reflected in Fig. 23b–d. A change in $\theta _{w}$ also alters the values of $k_{x}$ and $k_{y}$, which consequently change the final location of the wave packets. This is illustrated in a comparison between Fig. 23c and d where $\theta _{w}$ was changed from $\pi /6$ to $\pi /4$. Significant changes in the x and y position of the wave packets as predicted by the analytical group velocities (green arrow) were found. The error in group velocities leads to directional error in two-dimensional space.

6.2 Numerical Results—Effect of Varying Spatial Discretization Schemes for the Convective Term

In this subsection, the case which uses the hybrid spectral/finite difference discretization (CDS4/FS) will be compared to the cases using full Fourier spectral and full finite difference spatial discretization (CDS4/CDS4). The analytical group velocity results shown here correspond to the results presented in Fig. 15a–b. The aim is to illustrate the conditions for which the hybrid FD/FS case is more or less accurate than the full FD case when compared to solutions based on the FS case numerically. In two-dimensional space, the resolution characteristic of the numerical schemes is a function of both the reduced wavenumber as well as the wave propagation angle. In the case of the hybrid discretization, the propagation of dispersive errors from one spatial direction to the spatial direction based on the Fourier spectral discretization method is dependent on the wave propagation angle. The initial conditions used for these numerical experiments are listed in Table 3. A large value of reduced wavenumber, $k \varDelta = 0.9$ was chosen for these initial conditions in order to better reflect the dispersive error from the spatial discretization scheme. The numerical results of the CDS4/FS case alongside the analytical group velocity results of the FS/FS, CDS4/FS and CDS4/CDS4 cases are illustrated in Fig. 24.

Fig. 24

Comparison of physical group velocities (FS/FS) and numerical group velocities (CDS4/FS , CDS4/CDS4) of wave packets at $t_{max}=3$ for different initial conditions. Numerical solutions are computed using RK4 temporal and CDS4/FS spatial discretization in x and y, $M_{x}$ = 0.5, $M_{y}$ = 0.0, $R_{a}$ = 0.1, $Re_{a}$ = $\infty$. a Case B1, b Case B2, c Case B3, d Case B4 e Case B5, f Case B6. Wave packets are denoted as either $\omega ^{*}_{1}$ or $\omega ^{*}_{2}$ depending on the acoustic mode it represents

Table 3 Initial conditions of reduced wavenumbers and wave propagation angle

For comparison sake, the group velocity results of the full FD cases for $\omega ^{*}_{1}$ (blue arrows) are plotted beside the group velocity results of the hybrid discretization cases (green arrows) in Fig. 24. For $\theta _{w} = \pi /24$ and $\theta _{w} = \pi /12$, it was found that both the full FD and hybrid discretization case have the same group velocity in the y direction. This is illustrated in Fig. 24a, b where the blue arrow belonging to the CDS4/CDS4 case overlaps with the green arrows of the CDS4/FS case. This is because at such angles, the dispersive error from the CDS4 scheme in the y direction for the full FD case is almost negligible. The reason for the difference in y position of the wave packet for the hybrid discretization and full FS discretization case lies in the propagation of dispersion error from the FD scheme in the x direction to the y direction. For $\theta _{w}$ equal to $\pi /8$ and $\pi /6$, it was found that the y position of the wave packet for the full FD discretization case is closer to the exact (FS/FS) solution as compared to the hybrid discretization case. This is illustrated in Fig. 24c, d where the y position of the blue arrows are located closer to the black arrows. As before, the reason for the difference in the y position of the wave packet for the hybrid discretization and full FS discretization case lies in the propagation of the dispersion error from the FD scheme in the x spatial direction to the y spatial direction. For wave propagation angle of $5 \pi /24$ and $\pi /4$, it was found that the y position of the wave packet for the hybrid discretization case is closer to the exact (FS/FS) solution as compared to the full FD discretization case. This is illustrated in Fig. 24e and f where the y position of the green arrows are located closer to the black arrows. Although the dispersion error propagates from the FD scheme in the x direction to the y direction for the hybrid discretization case, the dispersion error originating from the FD scheme in the y direction of the full FD discretization case is larger and hence the results of the full FD discretization case is less accurate. A similar trend is observed for the comparison between the $V^{*}_{grp,full}$ of the DRP4/DRP4 and DRP4/FS cases in Fig. 25. For $\theta _{w}$ = $\pi /24$, $\pi /12$ and $\pi /8$, $V^{*}_{grp,full}$ of the DRP4/DRP4 and DRP4/FS cases have the same results. For $\theta _{w} = \pi /6$, $V^{*}_{grp,full}$ of the DRP4/DRP4 is closer to the exact result as compared to the DRP4/FS case. At large wave propagation angle ($\theta _{w} = 5\pi /24$ and $\pi /4$), $V^{*}_{grp,full}$ of the DRP4/FS is closer to the exact result as compared to the DRP4/DRP4 case. It should be pointed out that the observation made here is only valid under the consideration of large reduced wavenumber.

From a physical perspective, if the grid is sufficiently refined in both spatial direction, then the wave would propagate at the correct speed and direction irrespective of the type of spatial discretization methods employed and wave propagation angle considered. However, differences in the solutions of the different discretization method become obvious once the grid become sufficiently coarse. For a sufficiently coarse grid where $k \varDelta > rsim 0.5 \pi$, differences in $y_{pos}$ between the exact and full FD case can become significant for $\theta _{w}$ close to and greater than $\pi /4$. This implies that the waves travel in the wrong direction with increasing misalignment from the grid.

Fig. 25

Comparison of physical group velocities (FS/FS) and numerical group velocities (DRP4/FS , DRP4/DRP4) of wave packets at $t_{max}=3$ for different initial conditions. Numerical solutions are computed using RK4 temporal and DRP4/FS spatial discretization in x and y, $M_{x}$ = 0.5, $M_{y}$ = 0.0, $R_{a}$ = 0.1, $Re_{a}$ = $\infty$. a Case B1, b Case B2, c Case B3, d Case B4 e Case B5, f Case B6. Wave packets are denoted as either $\omega ^{*}_{1}$ or $\omega ^{*}_{2}$ depending on the acoustic mode it represents

6.3 Numerical Results—Effect of Varying Spatial Discretization Schemes for the Diffusive Term

In this sub section, the effect of varying spatial discretization schemes for the diffusive term on the numerical amplification is investigated numerically. For sufficiently low CFL number, the use of central type schemes for the viscous term only causes numerical dissipation. In order to understand the effect of the viscous term on the amplitude of the numerical solution, different spatial discretization schemes were used for the viscous terms of the LCNSE. The evolution of the amplitude of the numerical solution was analysed. The numerical parameters and schemes used for this analysis is listed in Table 4.

Table 4 Numerical parameters and schemes of different cases

Fig. 26

Comparison of $P^{\prime }_{max}$ as a function of time for different spatial discretization schemes used on the viscous term, a $k \varDelta = 0.9 \pi$ and $\theta _{w} = \pi /8$, b $k \varDelta = 0.9 \pi$ and $\theta _{w} = \pi /4$. FS/FS scheme is used for the convective terms. $Re_{a} = 200$, $R_{a} = 0.1$ and RK4 temporal scheme

In order to verify the analytical results shown in Fig. 21, the amplitude of the numerical solution, $P^{\prime }_{max} (t)$ is compared for different spatial discretization scheme used for the viscous term. Two different initial conditions are considered; $k \varDelta = 0.9 \pi , \theta _{w} = \pi /8$ and $k \varDelta = 0.9 \pi , \theta _{w} = \pi /4$. The aim is to compare the amplitude of the numerical solution for different wave propagation angle at high reduced wavenumber. The numerical results shown here correspond to the analytical results shown in Fig. 21a. For the case of $\theta _{w} = \pi /8$, the analytical results show that the case of the CDS4/CDS4 and CDS4/FS scheme have the same numerical dissipation while the case of the DRP4/DRP4 and DRP4/FS scheme have the same numerical dissipation as well. However, both the DRP4/DRP4 and DRP4/FS cases have a larger numerical dissipation than the case of the CDS4/CDS4 and CDS4/FS scheme. As before, the case with the FS/FS scheme have the most numerical dissipation. This same trend can be seen in Fig. 26a, where the maximum pressure fluctuation amplitude is plotted as a function of time for the different cases. For the case of $\theta _{w} = \pi /4$, the analytical results show that the case with the CDS4/CDS4 scheme have the least numerical dissipation. This is followed by the CDS4/FS, DRP4/DRP4, DRP4/FS and finally the FS/FS scheme. The case with the FS/FS scheme have the most numerical dissipation. This same trend can be seen in Fig. 26b, where the maximum pressure fluctuation amplitude is plotted as a function of time for the different cases. The analytical and numerical results illustrated here confirms the finding that the numerical dissipation is a function of both the reduced wavenumber and wave propagation angle.

7 Conclusions

In this paper, it was shown that, in the limit of spatial discretization assumption, a hybrid spectral/FD discretization method always gives a more accurate solution than a full FD discretization method for the 2D advection equation. But this is not necessarily true for the LCNSE. This is primarily because of the acoustic term within the dispersion relation of the LCNSE which leads to a coupling of the resolution of the discretization schemes in both spatial directions. To better understand the comparison of the dispersion properties between a hybrid and full FD scheme, the integrated sum of the error between the numerical group velocity $V^{*}_{grp,full}$ and the exact solution across all wavenumbers for a range of wave propagation angle is examined. In the comparison between a hybrid and full FD discretization schemes, the CDS4, DRP4 and CCOM6 schemes share the same characteristics. At low wave propagation angle, the integrated error of the full FD and hybrid discretization schemes remain the same. At intermediate wave propagation angle, the integrated error of the full FD schemes become smaller than the hybrid scheme. At large wave propagation angle, the integrated error of the full FD schemes diverges while the integrated error of the hybrid discretization schemes converges to zero.

At sufficiently low CFL number, the use of a central difference discretization on the viscous term only leads to error in the numerical dissipation. It was found that the numerical dissipation of the viscous term which is represented by its numerical amplification factor is lesser than the actual physical dissipation at high reduced wavenumber ($k \varDelta > rsim 0.6 \pi$ ) for the different central spatial discretization schemes. The numerical dissipation in this context represents the total diffusion obtained by the FD discretization of the diffusion term. The actual physical dissipation which is represented by the FS/FS discretization case is only a function of the cell Reynolds number. For the CDS4/CDS4, DRP4/DRP4, CCOM6/CCOM6 and $\hbox {CDS4}_{{opt}}$/$\hbox {CDS4}_{{opt}}$ cases, the wave propagation angle at which the numerical diffusion of the viscous term approaches the true physical diffusion (FS/FS case) occurs at $\pi /4$. As for the hybrid discretization case, the numerical diffusion of the viscous term matches the exact physical diffusion after a specific wave propagation angle for every $k \varDelta$ considered.

Appendix

The modified wavenumber relation of the CDS4 scheme is given as

$$\begin{aligned} (k_{x}^{*} \varDelta x)_{CDS4} = \frac{4}{3}\sin (k_{x} \varDelta x) - \frac{1}{6} \sin (2 k_{x} \varDelta x). \end{aligned}$$

(41)

The modified wavenumber relation of the DRP4 scheme is given as

$$\begin{aligned} \begin{aligned}&(k_{x}^{*} \varDelta x)_{DRP4} = 2a \sin (k_{x} \varDelta x) - 2b \sin (2k_{x} \varDelta x) + 2c \sin (3k_{x} \varDelta x),\\&a = 0.79926643 \quad b = 0.18941314 \quad c = 0.02651995. \end{aligned} \end{aligned}$$

(42)

The corresponding modified wavenumber in the y direction is obtained by taking $\varDelta y$ to be $\varDelta x$ and replacing $k_{x}$ with $k_{y}$. The modified wavenumber relations of the CDS4 and DRP4 schemes contains no imaginary component which makes the scheme inherently non-dissipative. For compact schemes, the numerical approximation of the spatial derivative takes on a more general form:

$$\begin{aligned} \sum _{m=-M}^{+M} b_{m} \ \bigg ( \frac{\partial \varPhi }{\partial x} \bigg )_{i+m} = \frac{1}{\varDelta x} \sum _{j=-L}^{+L} a_{j} \varPhi _{i+j}. \end{aligned}$$

(43)

By applying a Fourier transform of the variables on both sides of Eq. (43), one can obtain a modified wavenumber relation of the following form for the CCOM6 scheme [12]:

$$\begin{aligned} (k_{x}^{*} \varDelta x)_{CCOM6} =\frac{\frac{14}{9}\sin (k_{x} \varDelta x) + \frac{1}{18} \sin (2 k_{x} \varDelta x)}{1 + \frac{2}{3}\cos (k_{x} \varDelta x)}. \end{aligned}$$

(44)

The corresponding modified wavenumber in the y direction is obtained by taking $\varDelta y$ to be $\varDelta x$ and replacing $k_{x}$ with $k_{y}$. The modified wavenumber relation of the $CDS4_{opt}$ scheme is given as [25]:

$$\begin{aligned}&\begin{aligned} (k_{x}^{*} \varDelta x)_{CDS4_{opt}} =&\frac{1}{6 (1 + \beta )} \Bigr \{ 8 \sin (k_{x} \varDelta x) - \sin (2 k_{x} \varDelta x) + \frac{\beta }{2} \Bigr ( 8 \sin (k_{x} \varDelta x + k_{y} \varDelta y) \\&-\sin (2k_{x} \varDelta x + 2k_{y} \varDelta y) + 8 \sin (k_{x} \varDelta x - k_{y} \varDelta y) \\&-\sin (2k_{x} \varDelta x - 2k_{y} \varDelta y) \Bigl ) \Bigl \}, \end{aligned} \end{aligned}$$

(45)

$$\begin{aligned}&\begin{aligned} (k_{y}^{*} \varDelta y)_{CDS4_{opt}} =&\frac{1}{6 (1 + \beta )} \Bigr \{ 8 \sin (k_{y} \varDelta y) - \sin (2 k_{y} \varDelta y) + \frac{\beta }{2} \Bigr ( 8 \sin (k_{x} \varDelta x + k_{y} \varDelta y) \\&-\sin (2k_{x} \varDelta x + 2k_{y} \varDelta y) - 8 \sin (k_{x} \varDelta x - k_{y} \varDelta y) \\&+ \sin (2k_{x} \varDelta x - 2k_{y} \varDelta y) \Bigl ) \Bigl \}. \end{aligned} \end{aligned}$$

(46)