V.K. SUMAN et al

*V.K. Suman, CTFD Division, CSIR-NAL.

CTFD Division, CSIR-NAL, India.

# Determining Parameter Ranges for High Accuracy Large Eddy Simulation by Lax-
Wendroff Method††thanks: Parameters for High Accuracy LES for Lax-Wendroff
Method

V.K. Suman  Soumyo Sengupta  P. Sundaram  Aditi Sengupta  Tapan K. Sengupta
Computational & Theoretical Fluid Dynamics Division, CSIR-NAL, Karnataka,
India  High Performance Computing Laboratory, Aerospace Engineering, Indian
Institute of Technology Kanpur, Uttar Pradesh, India  CERFACS, 42 Avenue G.
Coriolis, 31057, Toulouse Cedex 1, France  Dept. of Mechanical Engineering,
Indian Institute of Technology (ISM) Dhanbad, Jharkand, India
<EMAIL_ADDRESS>
(9 August 2022; <Date> <Month> <Year>; <Date> <Month> <Year>)

###### Abstract

[Abstract]The analysis of Lax-Wendroff (LW) method is performed by the generic
modified differential equation (MDE) approach in the spectral plane using
Fourier transform. In this approach, the concept of dispersion relation plays
a major role relating spatial and temporal dependence of the governing
differential equation, including initial and boundary conditions in developing
high accuracy schemes. Such dispersion relation preserving schemes are
calibrated in the spectral plane using the global spectral analysis for the
numerical method in the full domain. In this framework, the numerical methods
are calibrated by studying convection and diffusion as the underlying physical
processes for this canonical model problem. In the LW method spatial and
temporal discretizations are considered together, with time derivatives
replaced by corresponding spatial derivatives using the governing equation.
Here the LW method is studied for the convection-diffusion equation (CDE) to
establish limits for numerical parameters for an explicit central difference
scheme that invokes third and fourth spatial derivatives in the MDE, in its
general form. Thus, for the LW method, two different MDEs are obtained,
depending on whether the LW method is applied only on the convection operator,
or both on the convection and diffusion operators. Motivated by a one-to-one
correspondence of the Navier-Stokes equation with the linear CDE established
in “Effects of numerical anti-diffusion in closed unsteady flows governed by
two-dimensional Navier-Stokes equation- Suman et al. Comput. Fluids, 201,
104479 (2020)”, an assessment is made here to solve flow problems by these two
variants of the LW method. Apart from mapping the numerical properties for
performing large eddy simulation for the LW methods, simulations of the
canonical lid-driven cavity problem are performed for a super-critical
Reynolds number for a uniform grid.

###### keywords:

Global spectral analysis; Error dynamics; Lax-Wendroff method; Large eddy
simulation

††articletype: Article Type00footnotetext: Abbreviations: CDE, convection-
diffusion equation; GSA, global spectral analysis; NSE, Navier-Stokes
equation; LW, Lax-Wendroff; MDE, modified differential equation

##  1 Introduction

Progress in scientific computing has been made possible due to the
developments in hardware and analysis of numerical methods. Classical analysis
use the order of the numerical methods based on Taylor’s series expansion for
the spatial and temporal discretization adopted in solving the governing
partial differential equation. Different governing equations require a priori
analysis for suitability of a numerical method. For example, using uniformly
second order central differencing scheme for space and time derivatives led to
a solution for the canonical one-dimensional (1D) convection equation,

$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0$ (1)

with the solution remaining bounded strictly 1 for this mid-point leap frog
method 2. But, for the canonical heat equation,

$\frac{\partial u}{\partial t}=\alpha\frac{\partial^{2}u}{\partial x^{2}}$ (2)

the numerical method becomes unconditionally unstable for the same second
order discretization 2, 3 for this parabolic partial differential equation,
with the method attributed to Richardson 4. The hyperbolic partial
differential equation (1) with the midpoint leap frog method played a major
role in earlier weather prediction.

The Fourier series analysis of heat equation by von Neumann and Richtmyer 5,
6, 7 explained why the second order discretization of heat equation is
unconditionally unstable, which is due to the creation of a spurious unstable
numerical mode. Despite the early success of this numerical stability study,
it is noted 3 that Fourier series analysis has some unresolved issues that
prompted Zingg 8 to state that through Fourier analysis, one can evaluate the
phase and amplitude error of a given method as a function of wavenumber.
However, this information can be difficult to interpret. The author is
referring to the inability of Fourier series analysis method to explain the
correct dispersion error of any numerical method for space-time dependent
governing equation.

The main limitation of Fourier analysis is due to its strict validity for
spatially periodic problem only, for the governing differential equations with
constant coefficient. Also one performs a normal mode analysis, where modes do
not interact with each other, even for a linear governing equation. Thus,
despite the Fourier analysis method being easy to apply, it is not considered
further due to its inability to solve general non-periodic problems while
including effects of boundary conditions.

Analysis of the heat equation (Eq. (2)) is mainly used for numerical stability
at the steady state, without considering about the transient state. This is
not like the convection equation (Eq. (1)), for which the numerical solution
must propagate at the correct speed, without dissipation and dispersion. Thus,
this is a more relevant canonical problem for numerical analysis, as noted in
1. The main application for this equation was for numerical weather prediction
9, where tracking the propagating disturbances (as resolved into acoustic,
vortical and entropic waves) are vitally important 3, 9, 10, 11. Relevance of
this equation is explained briefly in terms of the concept of dispersion
relation preservation (DRP) property in the following section.

The numerical analysis of Eq. (1) involves the dependence of the unknown on
space and time simultaneously, and one writes the general representation of
the unknown in terms of Fourier-Laplace transform as,

$u(x,t)=\int\int\hat{U}(k,\omega_{0})e^{i(kx-\omega_{0}t)}dk\;d\omega_{0}$ (3)

Substituting Eq. (3) in (1), one gets the physical dispersion relation as,

$\omega_{0}=kc$ (4)

This relates the spatial and temporal scales in the spectral plane. For the
computation of Eq. (1), one notes the spatial and temporal discretizations to
be such that Eq. (4) is satisfied for every wavenumber, with a corresponding
circular frequency for the DRP schemes. For such schemes, one can show that
the energy of the system travels at the group velocity defined 3 by,

$v_{g}=\frac{d\omega_{0}}{dk}$

Comparing Eq. (1) with Eq. (4), it is noted that $v_{g}=c$. One can also show
that if the initial condition to solve Eq. (1) is given by,
$u_{0}(x)=\int\hat{U}_{0}(k)e^{ikx}dk$, then the general solution at a later
time is given by, $u(x,t)=u_{0}(x-ct)$. This implies that the initial
amplitude remains unchanged, as the solution propagates downstream with the
speed $c$, for all the spatial scales associated with the initial condition.
Thus, the solution of Eq. (1) is non-dissipative and non-dispersive. One notes
that $u(x,t)$ in Eq. (3) via the integration of the variable (along its strip
of convergence) ensures all modal and nonmodal interactions. This is the basis
of the global spectral analysis (GSA), which is distinct from the von Neumann
analysis. GSA has been advocated earlier in 11, 12 and comprehensive details
can be noted in 3.

In 3, 9, the major achievement was in identifying the correct numerical
dispersion relation. The authors underscored the importance of spatial and
temporal discretizations to obtain a numerical phase speed which is different
from the physical phase speed, such that the numerical dispersion relation is
given by,

$\omega_{num}=kc_{num}$ (5)

Note that in GSA, $x$ or $k$ ($x$ is the space variable and $k$ is wavenumber)
is the independent variable, while the temporal scale is related to $k$ via
the numerical dispersion relation. A DRP method is one, for which the
departure of the numerical dispersion relation from the physical dispersion
relation is negligible for accessible ranges of numerical parameters. The
quantitative analysis used in GSA 3, 13, 14 is based on the modified
differential equation (MDE) approach that is described next.

Notable examples of the MDE are for (i) more accurate numerical methods for
parabolic equations in Milne’s method 15, and (ii) for ‘improved’ numerical
stability in Dufort-Frankel’s method 16. It is known that increased stability
for Dufort-Frankel’s method 16 also brings in the problem of inconsistency 2,
3 for some numerical parameter combinations. In the MDE, one reconverts the
discretized equation to an equivalent differential equation, as shown in 17.
Linear problems with variable coefficients have been studied in 19 and
nonlinear problems are shown in 20 by the MDE. However, Li and Yang 18 made a
generalization that the MDE is very heuristic, unfortunately just valid for
solutions in smooth regions or at low frequency modes 17: Therefore the
connection with the von Neumann analysis is only restricted there, and cited
21 that there is a lack of theoretical foundation for the MDE.

The MDE continues to be used in Lax-Wendroff (LW) method for hyperbolic
partial differential equations 22, 23. Apart from these, one notes that the
Russian school 24, 25 explained different aspects of the MDE. Acoording to
their classification, if one converts the difference equation into an
equivalent differential equation form by retaining both space and time
derivatives, then it is called the $\Gamma$\- form analysis. On the other
hand, if the discrete equation is converted back in the differential form with
all the truncation terms converted in terms of the spatial derivatives, then
that is termed as the $\Pi$-form analysis. Variants of the $\Gamma$\- and
$\Pi$-forms of the MDE along with Fourier-Laplace representation of the
unknowns are used in development of the GSA and explained in the following.

The LW method is one of the popular methods for higher order solution of PDEs
and continues to be used in the CFD community 39, 40, 41, 42, 52, 53, 54.
Despite its popularity, an accurate and detailed analysis of the method is
lacking in the literature. Recently, Soumyo et al. have analyzed the LW method
in 1D and 2D for the linear convection-diffusion problem using GSA and they
have identified the stability limits and acceptable range of simulation
parameters for the method 52. They have also validated their analyses by
solving 2D Navier-Stokes equation for the Taylor-Green vortex problem.
Although the authors have performed detailed analyses with respect to the
characterization of the method for explicit CD2 scheme, important questions
remain pertaining to its variants and general characteristics for
multidimensional problems.

It is known that for governing equations of the convection-diffusion type, two
variants of the LW scheme can be derived based on the treatment of convection
or convection-diffusion terms with the former variant analyzed in 52. However,
these have not been analyzed in the literature with a view to compare/suggest
the differences between the two. This is addressed in detail in the present
manuscript using GSA and comparing/contrasting the relevant metrics. Also, it
is important to note that the method has not been correctly analyzed with
respect to error dynamics in the multidimensional space. This is also
addressed in the present research where it is shown that the LW method
includes an additional cross-derivative term which can impact its
stability/accuracy properties. Finally, optimal parameters for the method are
determined for accurate solution of the 2D Navier-Stokes equations and they
are corroborated by computing the flow inside a square lid driven cavity at
supercritical Reynolds numbers of 10,000. The computations with the present
scheme are compared with available benchmark solutions 48, 49 which to the
author’s knowledge has never been done before.

The paper is organized as follows. In section 2, the basic principles and
rationale for GSA is presented. This is followed by section 3 describing the
Lax-Wendroff method for the 1D linear convection equation. In section 4, the
Lax-Wendroff method is derived for the 1D linear convection-diffusion equation
and two of its variants are analyzed using GSA. Multidimensional analysis of
the method is presented in section 5 for the 2D convection-diffusion equation
and the optimal parameters are identified. The GSA results are corroborated in
section 6 by solving the 2D Navier-Stokes equations for the flow inside a
square lid driven cavity and comparing the results with benchmark data. The
paper ends with summary and conclusions in section 7.

##  2 Quantifying Methods by Global Spectral Analysis (GSA)

In GSA, the unknown is represented by a hybrid-spectral form as,

$u(x,t)=\int\hat{U}(k,t)e^{ikx}dk$ (6)

Here $\hat{U}$ is the time-dependent Fourier amplitude and $k$ is the
independent variable. The exact spatial derivative is obtained from Eq. (6)
as, ${\frac{\partial u}{\partial x}\biggl{|}}_{{}_{{}_{exact}}}=\int
ik\hat{U}e^{ikx}dk$. One can write an equivalent numerical spatial derivative
as, ${\frac{\partial u}{\partial x}\biggl{|}}_{{}_{{}_{num}}}=\int
ik_{eq}\hat{U}e^{ikx}dk$. The Fourier-Laplace amplitude is multiplied by
$ik_{eq}$ for the numerical derivative (instead of $ik$ for the exact
derivative). This notation for the numerical derivative was originally made
popular by Vichnevetsky and Bowles 26 and subsequently in 27, 28, 29, among
many other references.

There are two main reasons for using $k_{eq}$ in the spectral representation:
first, to provide the yardstick in comparing different discretization methods
used in any computation by evaluating an equivalent $k_{eq}$ for the method.
Expressions for this quantity have been presented for finite volume and
Galerkin finite element methods for the solution of Eq. (1) in 3. Ideally the
ratio $k_{eq}/k$, should be equal to one and it is a measure of the resolution
of the discretization scheme. One notes that this ratio is presented in the
literature as a function of $kh$, with $h$ as the uniform grid spacing. The
second use of $k_{eq}$ is due to its erroneous use for numerical dispersion
relation. If one ignores error due to temporal discretization, then the form
of the spatial derivative given in terms of $k_{eq}$ for Eq. (1) has been used
incorrectly by few researchers to write the numerical dispersion relation as,

$\omega_{num}=k_{eq}c$ (7)

This has also been noted as correct for the semi-discrete numerical stability
analysis 30, 31, 32. It may appear as correct to treat $c$ as a constant in
Eq. (7). However, a correct quantitative analysis is impossible by the semi-
discrete approach using the numerical dispersion relation given by Eq. (7). It
is noted that one can compute numerical group velocity for this approach
following the basic definition as 3, 33, 34, 35,

$v_{g,num}=\frac{d\omega_{num}}{dk}=c\frac{dk_{eq}}{dk}$ (8)

However, the obtained group velocity is incorrect, as it is independent of
time discretization scheme for Eq. (1). Propagation of wave-packet studied in
3 show the group velocity is a strong function of both space and time
discretizations. The authors in 36 have used this numerical dispersion
relation in Eq. (8) to derive a DRP scheme by considering spatial
discretization alone. A four time-level method 36 was used to propose their
DRP scheme, failing to note that such methods will invoke two spurious
numerical modes, a topic explained in 3, 37, 38.

Instead of using Eq. (7), authors in 9, 13 have shown that when one solves Eq.
(1), the phase speed can no longer be a constant. This may appear paradoxical,
but it is clearly explained in 13, 14, 37 that the space-time discretization
method fixes the phase shift per time step, which in turn determines the
numerical phase speed ($c_{num}$), that is different from the physical phase
speed, $c$. Then $c_{num}$ becomes function of wavenumber/circular frequency,
giving rise to phase and dispersion error. This simple, yet subtle cause for
$c_{num}\neq c$, is one of the central results of GSA in the expression for
the numerical dispersion relation given in Eq. (5). The correct group velocity
accounting for the space-time discretization is therefore given by,

$v_{g,num}=c_{num}+k\frac{dc_{num}}{dk}$ (9)

We note that $k$ is the independent variable, and the spatial and temporal
discretization schemes fix the numerical dispersion relation, which in turn
fixes other numerical parameters and error dynamics 3, 43. The fact that $c$
changes to $c_{num}$, applies equally to other coefficients of many other
transport and diffusion equations 43, 44, 47.

###  2.1 Rationale for the GSA

To justify the use of GSA, an explanation is provided with the help of an
example to demonstrate the utility of Eq. (5) for the 1D convection equation
solved by midpoint leap-frog scheme 4 on a uniform grid of spacing $h$.
Discretizing Eq. (1), one obtains the difference equation for a node at
($x_{j},t^{n}$) as,

$\frac{u^{n+1}_{j}-u^{n-1}_{j}}{2\Delta
t}+c\left(\frac{u^{n}_{j+1}-u^{n}_{j-1}}{2h}\right)=0$ (10)

where $\Delta t$ is the time-step used. We note that the numerical phase shift
derived using Eq. (6) is different from that given by the $\Pi$-form analysis
24, due to the fact that the numerical dispersion relations of these two
approaches are different.

For the discrete Eq. (10), the unknown is next represented by Eq. (6). In the
$\Pi$-form analysis, the numerical phase speed is computed from the dispersion
relation as $c_{{num}_{\Pi}}=\frac{k_{eq}}{k}c$ where
$\omega_{\Pi}=kc_{{num}_{\Pi}}$. For the second order central difference
scheme, $k_{eq}=\frac{\sin(kh)}{h}$. Thus, the numerical phase speed from
$\Pi$-form for the leap-frog scheme is given by

$\frac{c_{{num}_{\Pi}}}{c}=\frac{\sin(kh)}{kh}$ (11)

The $\Pi$-form analysis in 26 reports the above expression (cf. Eq. (2.13)).
There are two distinctive features of this result. First, the time-integration
is by a three-time level method, and one gets two distinct numerical phase
speeds (not given here). Secondly, and most importantly, the assumption that
the numerical phase speed is independent of time discretization is used in
this $\Pi$-form analysis. A correct approach based on the $\Gamma$-form
approach is given next.

The circular frequency and phase speed is obtained as per the governing
differential equation and its discretization in writing the numerical
dispersion relation by Eq. (5). Based on the difference equation, one obtains
the numerical amplification factor following the representation in Eq. (6).
Such a numerical amplification factor fixes phase shifts per time step, and
provides the numerical phase speeds. For the $\Gamma$-form analysis, the
unknown is represented by Eq. (6). One can represent the initial condition for
Eq. (1) as,

$u(x_{j},t=0)=u_{j}^{0}=\int\hat{U}_{0}(k)\;e^{ikx_{j}}dk$ (12)

with the subscript and superscript denoting the spatial and temporal indices,
respectively. The solution at any time, $t=n\Delta t$, is written using the
numerical amplification factor as

$u_{j}^{n}=\int\hat{U}_{0}(k)\;|G_{num,j}|^{n}\;e^{i(kx_{j}-n\phi_{j})}dk$
(13)

where $G_{num,j}=\left(\frac{\hat{U}(k,t^{n}+\Delta
t)}{\hat{U}(k,t^{n})}\right)$ is the complex amplification factor, i.e.
$G_{num,j}=G_{num,rj}+iG_{num,ij}$, such that
$|G_{num,j}|=(G_{num,rj}^{2}+G_{num,ij}^{2})^{1/2}$. The phase shift per time
step is calculated as $\tan\phi_{j}=-{G_{num,ij}}/{G_{num,rj}}$. From this the
numerical phase speed ($c_{{num}_{\Gamma}}$) is obtained as,

$c_{{num}_{\Gamma}}=\frac{\phi_{j}}{k\Delta t}$ (14)

Here, $c_{{num}_{\Gamma}}$ depends on $k$, i.e. the numerical solution is
dispersive, as opposed to the non-dispersive nature of physical solution.
Thus, both the $\Pi$\- and $\Gamma$-form analyses show the numerical phase
speed to depend on $k$, with the main difference that the latter uses the
temporal discretization information. In contrast, the $\Pi$-form analysis
(like the semi-discrete analysis) ignores any information originating from the
temporal discretization.

Substituting the variable $u$ given by Eq. (6) in Eq. (10), one obtains the
following,

$\hat{U}^{n+1}_{j}-\hat{U}^{n-1}_{j}+\frac{c\Delta
t}{h}\left(e^{ikh}-e^{-ikh}\right)\hat{U}^{n}_{j}=0$ (15)

where the variables with the hat denote spectral amplitudes. With the
definition of numerical amplification factor
$G_{num,j}=\left(\frac{\hat{U}^{n+1}_{j}}{\hat{U}^{n}_{j}}\right)=\left(\frac{\hat{U}^{n}_{j}}{\hat{U}^{n-1}_{j}}\right)$,
a quadratic for $G_{num,j}$ is obtained with the two roots obtained as,

${G_{num,j}}_{1,2}=-iN_{c}\sin(kh)\pm\sqrt{1-{N_{c}}^{2}\sin^{2}(kh)}$ (16)

where $N_{c}=\frac{c\Delta t}{h}$ is the Courant-Friedrich-Lewy (CFL) number.
From the above, $c_{{num}_{\Gamma}}$ is calculated from Eq. (14) as

$\frac{c_{{num}_{\Gamma_{1,2}}}}{c}=\frac{\phi_{j}}{(kh)N_{c}}=\left(\frac{1}{(kh)N_{c}}\right)\tan^{-1}\left({\frac{N_{c}\sin(kh)}{\pm\sqrt{1-{N_{c}}^{2}\sin^{2}(kh)}}}\right)$
(17)

Vichnevetsky and Bowles 26 also used the $\Gamma$-form analysis to determine
the numerical phase speed. The unknown is represented by the Fourier-Laplace
transform and substituting in Eq. (10), one obtains the spectral plane
representation for the leap-frog scheme as,

$e^{-i\omega_{0}\Delta t}-e^{i\omega_{0}\Delta
t}+N_{c}\left(e^{ikh}-e^{-ikh}\right)=0$

which upon simplification yields, $\sin(\omega_{0}\Delta t)=N_{c}\sin(kh)$.

This provides the circular frequency as,

$\omega_{0}=\frac{1}{\Delta t}\sin^{-1}(N_{c}\sin(kh))$ (18)

With the expression for $\omega_{0}$ in terms of $k$, numerical phase speed is
computed as, $c_{{num}_{VB}}=\frac{\omega_{0}}{k}$, which upon simplification
yields

$\frac{c_{{num}_{VB}}}{c}=\frac{1}{(kh)N_{c}}\sin^{-1}(N_{c}\sin(kh))$ (19)

We note that Vichnevetsky and Bowles obtained the amplification factors given
by Eq. (4.8c) in 26 and show the above expression for numerical phase speed in
Table 4.3 26. At a first glance the numerical phase speed computed from the
$\Gamma$-form analysis (Eq. (17)) and the expression given by Vichnevetsky and
Bowles 26 (Eq. (19)) may appear as different. But, these are equivalent
following the trigonometric identity:
$\sin^{-1}(x)=\tan^{-1}\left(\frac{x}{\sqrt{1-x^{2}}}\right)$.

Having described the physical and different form of numerical dispersion
relations, and appearance of additional spurious numerical modes for three/
higher time level methods for solving Eq. (1), we discuss about the LW method,
which is unambiguous and helps explain some of the fundamental principles of
the MDE.

##  3 The Lax-Wendroff Method for the Convection Equation

In this method, the proponents retained second order accuracy for time
integration in discretizing the governing equation. For the MDE, we have an
example for which spurious modes are not invoked, while the time integration
is by a higher order method, with the DRP property improved. The difference
equation is obtained from the Taylor series expansion,

$u_{j}^{n+1}=u_{j}^{n}+\Delta t\frac{\partial u}{\partial t}+\frac{(\Delta
t)^{2}}{2}\frac{\partial^{2}u}{\partial t^{2}}$ (20)

In the LW method, the second time derivative is obtained from Eq. (1) as,
$\frac{\partial^{2}u}{\partial t^{2}}=c^{2}\frac{\partial^{2}u}{\partial
x^{2}}$. When the first and second time derivatives are substituted in Eq.
(20), one obtains the MDE given by,

$u(t+\Delta t)=u(t)-c\Delta t\frac{\partial u}{\partial x}+\frac{c^{2}(\Delta
t)^{2}}{2}\frac{\partial^{2}u}{\partial x^{2}}$ (21)

One immediately notices the presence of the last term on the right hand side
as a strictly diffusive term that was absent in the original governing
equation. This shows that a pure convection equation (a hyperbolic PDE) is
converted into CDE (a parabolic PDE). This added diffusive term may provide
numerical stabilization, as often used in CFD. However by Lax’s theorem, this
is an issue of inconsistency, as a hyperbolic PDE is converted to a parabolic
PDE 2, 3. This has been noted for Dufort-Frankel method 16, 3 also, where the
heat equation (parabolic PDE) is converted into a wave propagation problem –
another example of inconsistency, if the numerical parameters are not chosen
with care (see Fig. 6.4 in 3 for details). In defense of the LW method, one
should note that with refined time steps, the inconsistency due to added
diffusion can be made progressively sub-dominant.

##  4 The Lax-Wendroff Method for the 1D Convection-Diffusion Equation

In the previous section, the LW method is presented for the 1D convection
equation and a key issue of inconsistency of the MDE is demonstrated. Here, we
evaluate using GSA, the same method applied to a 1D convection-diffusion
equation (CDE). The objective in doing this study is two-fold. First objective
is to check if the MDE is consistent with the original PDE. The second
objective is to briefly analyze and contrast two forms of application of the
LW method- application to only convection term and both convection-diffusion
terms, respectively. The former strategy is employed by a popular and well
established research code- ABVP, developed at CERFACS, France, for problems of
combustion 39, 40, 41, 42.

However, no analysis or quantification is provided in the literature to
compare the above two strategies. This is established here by using GSA and
comparing the numerical properties of the two strategies for representative
simulation parameters. This is one of the contributions of the present work
and has never been done before.

We note that unlike for the convection equation, if the LW method is used in
the 1D CDE, the problem of inconsistency is not present as shown next.
However, the effects of additional diffusion terms will determine the accuracy
of space-time discretization simultaneously. For the CDE, the governing
canonical equation is given by,

$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial
x}=\alpha\frac{\partial^{2}u}{\partial x^{2}}$ (22)

For the LW method substitution of Eq. (22) in Eq. (20) yields the following,

$u(t+\Delta t)=u(t)-c\Delta t\frac{\partial u}{\partial x}+\alpha\Delta
t\frac{\partial^{2}u}{\partial x^{2}}+\frac{(\Delta
t)^{2}}{2}\frac{\partial^{2}u}{\partial t^{2}}$ (23)

Furthermore, deriving $\frac{\partial^{2}u}{\partial t^{2}}$ from Eq. (22) one
obtains

$\frac{\partial^{2}u}{\partial t^{2}}=c^{2}\frac{\partial^{2}u}{\partial
x^{2}}-2\alpha c\frac{\partial^{3}u}{\partial
x^{3}}+\alpha^{2}\frac{\partial^{4}u}{\partial x^{4}}$ (24)

Substitution of Eq. (24) in Eq. (23), one notices added diffusive terms
arising from the pure convection and diffusion terms. Additionally there is a
dispersive term due to the interaction between convection and diffusion terms.
A higher order diffusion term also appears from the basic diffusion term,
which can act like the hyper-viscosity term added in pseudo-spectral method to
provide numerical stabilization 51, when used with lower order Runge-Kutta
schemes.

###  4.1 The Lax-Wendroff method for explicit central difference scheme

One can use central difference schemes for various derivatives in Eqs. (22)
and (24), and substituting these in Eq. (23), one gets the corresponding
difference equation. Using the following notations:
$\frac{\partial^{2}u}{\partial x^{2}}=h^{2}D^{2}u$,
$\frac{\partial^{3}u}{\partial x^{3}}=h^{3}D^{3}u$ and
$\frac{\partial^{4}u}{\partial x^{4}}=h^{4}D^{4}u$, and the non-dimensional
parameters: $N_{c}=\frac{c\Delta t}{h}$ and $Pe=\frac{\alpha\Delta t}{h^{2}}$,
as the CFL and the Peclet numbers, the difference equation becomes,

$u_{j}^{n+1}=u_{j}^{n}-\frac{N_{c}}{2}(u^{n}_{j+1}-u^{n}_{j-1})+(Pe+\frac{N_{c}^{2}}{2})(u^{n}_{j+1}-2u_{j}^{n}+u^{n}_{j-1})-PeN_{c}D^{3}u_{j}^{n}+\frac{Pe^{2}}{2}D^{4}u_{j}^{n}$
(25)

Various differences appearing in Eq. (25) can be evaluated using full and
half-node locations of grid points. For example, half-node locations provide
the following, $D^{2}u_{j}=u_{j-1}-2u_{j}+u_{j+1}$ and
$D^{4}u_{j}=16(u_{j+1}-4u_{j+1/2}+6u_{j}-4u_{j-1/2}+u_{j-1})$. For
$D^{3}u_{j}$, a combination of full and half-node representations provide,
$D^{3}u_{j}=(u_{j+2}-2u_{j+1}+2u_{j-1}-u_{j-2})/2$. A formal Taylor series
expansion would reveal the consistency of these difference expressions. Note
that in some applications, the LW method is used for convection terms only 39,
40, 41, 42. In such cases, one would switch off the $D^{3}u_{j}$ and
$D^{4}u_{j}$ terms. In the following, the numerical properties of methods to
solve the CDE, with and without $D^{3}u_{j}$ and $D^{4}u_{j}$ terms, are
compared following the analysis given in 43.

For the 1D CDE, $c$ and $\alpha$ are constant real numbers and for a generic
non-periodic problem one can represent $u$ by Eq. (6) so that the CDE becomes,

$\frac{d\hat{U}}{dt}+ick\hat{U}=-\alpha k^{2}\hat{U}.$ (26)

which can be solved analytically for a general initial condition in Eq. (12)
to yield,

$\hat{U}(k,t)=\hat{U}_{0}(k)\>e^{-\alpha k^{2}t}e^{-ikct}.$ (27)

The physical dispersion relation of the CDE is obtained using Eq. (3) as,

$\omega_{0}=c\>k-i\>\alpha\>k^{2}.$ (28)

The physical phase speeds is obtained from,

$c_{phys}=\frac{\omega_{0}}{k}=c-i\>\alpha k,$ (29)

The phase speed is taken as the real part of the above, as reported in the
literature. The physical group velocity is the energy propagation velocity of
a wave-packet and is given following the definition in 3, 33 as,

$v_{g,phys}=\frac{d\omega_{0}}{dk}=c-2\>i\>\alpha k.$ (30)

The implication of this complex group velocity is given in 43 and for real
coefficient of diffusion, one obtains the real part of the group velocity
given by, $c$, and the same is used here for the CDE.

The analytical solution of the CDE can be interpreted by the physical
amplification factor, i.e. comparing the solution amplitude at two distinct
instants separated by $\Delta t$ so that,

$G_{phys}=\frac{\hat{U}(k,t+\Delta t)}{\hat{U}(k,t)}=e^{-\alpha\>k^{2}\Delta
t}e^{-i\>k\>c\>\Delta t}=e^{-i\>\omega_{0}\>\Delta t}$ (31)

which can be expressed in terms of the numerical parameters, $N_{c}$ and $Pe$
as,

$G_{phys}=e^{-Pe\>(kh)^{2}}\>e^{-i\>N_{c}\>(kh)}.$ (32)

It has been shown in 43, 44, 45 that contrary to the popular assumption of
von-Neumann analysis, the numerical convection speed, $c_{num}$, and the
numerical diffusion, $\alpha_{num}$, are dependent on $k$, $h$, $N_{c}$ and
$Pe$. All numerical methods have their numerical amplification factor,
$G_{num}$, and numerical dispersion relation which dictate the evolution of
the solution in time. It is essential that for any numerical scheme, $G_{num}$
should be as close as possible to $G_{phys}$.

The numerical discretization of the continuous problem would result in a
numerical dispersion relation equivalent to Eq. (28) as derived in 43, and
expressed as,

$\omega_{num}=c_{num}\>k-i\>\alpha_{num}\>k^{2}.$ (33)

where $\omega_{num}$ is complex and differs from the physical dispersion
expression, since $c_{num}$ and $\alpha_{num}$ vary with $kh$, $N_{c}$ and
$Pe$. Exploiting the similarity to the exact solution, the numerical
amplification factor, $G_{num}$ is obtained as a function of $\alpha_{num}$
and $c_{num}$ as,

$G_{num}=e^{-\alpha_{num}\>k^{2}\Delta t}\>e^{-i\>k\>c_{num}\Delta
t}=e^{-i\>\omega_{num}\Delta t}.$ (34)

From Eq. (25), for the LW method based on explicit central differences for the
1D CDE, one obtains the complex $G_{num}$ as

$\begin{split}G_{num}=\,&1-\frac{N_{c}}{2}\left(e^{ikh}-e^{-ikh}\right)+\left(\frac{N_{c}^{2}}{2}+Pe\right)\left(e^{ikh}-2+e^{-ikh}\right)\\\\[6.45831pt]
&-\frac{Pe\,N_{c}}{2}\left(e^{2ikh}-2e^{ikh}+2e^{-ikh}-e^{-2ikh}\right)\\\\[6.45831pt]
&+8Pe^{2}\left(e^{2ikh}-4e^{ikh/2}+6-4e^{-ikh/2}+e^{-2ikh}\right)\end{split}$
(35)

Using similar previous relations, one notes that:

$c_{num}=\Re\left(\frac{\omega_{num}}{k}\right),$ (36)

and,

$v_{g,num}=\Re\left(\frac{d\omega_{num}}{dk}\right).$ (37)

One can evaluate $c_{num}$ from the numerical phase shift per time step, as
the ratio between the imaginary and real part of $G_{num}$ given by,

$tan(\beta)=-\frac{\Im(G_{num})}{\Re(G_{num})}=tan(c_{num}k\Delta t).$ (38)

where the real and imaginary parts of $G_{num}$ are represented in the above
as $\Re(G_{num})$ and $\Im(G_{num})$, respectively. Also, note that for the
discussion of convection equation we have termed this phase shift per time
step as $\beta$.

This allows one to express the non-dimensional effective numerical phase speed
as,

$\frac{c_{num}}{c_{phys}}=\frac{\beta}{kc\Delta
t}=-\frac{1}{(kh)N_{c}}\>tan^{-1}\left[\frac{\Im(G_{num})}{\Re(G_{num})}\right].$
(39)

Similarly, the numerical group velocity is obtained as,

$v_{g,num}=\Re\left(\frac{\partial\omega_{num}}{\partial
k}\right)=\frac{1}{\Delta t}{\partial\beta\over\partial k},$ (40)

which can be written in normalized form as,

$\frac{v_{g,num}}{v_{g,phys}}=\frac{1}{N_{c}}\frac{d\beta}{d(kh)}.$ (41)

Following Eq. (32), one directly notes that the modulus of numerical
amplification factor is solely dependent on $\alpha_{num}$ by analogy as,

$|G_{num}|=e^{-\alpha_{num}k^{2}\Delta t},$ (42)

Alternately in terms of $Pe$, one can write this modulus also as,

$ln|G_{num}|=-\frac{\alpha_{num}}{\alpha}(kh)^{2}\>Pe.$ (43)

Equation (43) can also be used to evaluate the numerical diffusion coefficient
in non-dimensional form as,

$\frac{\alpha_{num}}{\alpha}=-\frac{ln|G_{num}|}{(kh)^{2}Pe}.$ (44)

This clearly establishes that a numerical scheme for the CDE changes the
diffusion process differently for different length scales ($kh$), depending on
the numerical scheme. For the LW method this corresponds to choosing the
spatial discretization only, and that is a special feature of the LW method.
From Eq. (25) one notes that the retention of higher order term for the time
integration will also involve the terms $D^{3}u_{j}^{n}$ and $D^{4}u_{j}^{n}$,
when the diffusion term is expressed with the second order time
discretization.

Using Eqs. (35), (44), (39) and (41), the properties
$\frac{|G_{num}|}{|G_{phys}|}$, $\frac{\alpha_{num}}{\alpha}$,
$\frac{c_{num}}{c_{phys}}$ and $\frac{v_{g,{num}}}{v_{g,phys}}$ are obtained
for any ($N_{c}$, $Pe$) combination. These property charts are then used to
gauge the accuracy of the developed scheme.

Figure 1: The ratio of numerical and physical amplification factors
$\left(\frac{|G_{num}|}{|G_{phys}|}\right)$ and the ratio of numerical and the
physical diffusion coefficients $\left(\frac{\alpha_{num}}{\alpha}\right)$ for
the Lax-Wendroff method based on explicit CD scheme applied on the convection
term only of the CDE plotted in the $(N_{c},kh)$-plane for the representative
Peclet numbers $Pe=0.01$ and $0.02$. Regions of numerical instability are as
marked in the panels. Figure 2: The ratio of numerical and physical phase
speeds $\frac{c_{num}}{c_{phys}}$ and the ratio of numerical and the physical
group velocity $\frac{v_{g,num}}{v_{g,phys}}$ for the Lax-Wendroff method
based on explicit CD scheme applied on the convection term only of the CDE
plotted in the $(N_{c},kh)$-plane for the representative Peclet numbers
$Pe=0.01$ and $0.02$.

In Figs. 1 and 2, property charts are shown in the $(N_{c},kh)$-plane for the
Peclet numbers, $Pe=0.01$ and $0.02$, for the LW method based on explicit
central differences as applied on the convection term of the governing CDE,
i.e. $D^{3}$ and $D^{4}$ terms in Eq. (25) are not included. From Fig. 1, one
notes that at low $N_{c}$ values, the ratio $\frac{|G_{num}|}{|G_{phys}|}$ is
greater than unity for all $k$, implying that the numerical solution has lower
numerical diffusion. This is corroborated from the plot of
$\frac{\alpha_{num}}{\alpha}$ contours where the corresponding values are
lower than one. Numerical diffusion is also noted to decrease progressively
with increasing $k$. However adjacent to this region, another region exists
with higher numerical diffusion compared to the physical value at all
wavenumbers with the diffusion strength increasing with $k$. Beyond $N_{c}=1$,
$\alpha_{num}$ is negative implying presence of anti-diffusion for some range
of resolved $k$ that will lead to catastrophic numerical instability as shown
in the marked regions in red in the bottom frames of Fig. 1.

The numerical phase speed ($c_{num}$) and group velocity ($v_{g,num}$)
contours plotted in Fig. 2, show poorer performance of the LW scheme in
preserving a signal’s physical dispersion and propagation characteristics as
compared to many compact schemes used with higher order Runge-Kutta methods 3.
This performance degradation is attributed to the poor resolution of explicit
CD scheme in calculating the spatial derivatives. However, some specific
ranges of $kh$ and $N_{c}$ can be found, where this method can be used for
very low error tolerance. This is one of the reasons for the reported analysis
here, which helps in locating specific $N_{c}$ and $kh$ ranges for a fixed
$Pe$.

Figure 3: The ratio of numerical and physical amplification factors
$\left(\frac{|G_{num}|}{|G_{phys}|}\right)$ and the ratio of numerical and the
physical diffusion coefficients $\left(\frac{\alpha_{num}}{\alpha}\right)$ for
the Lax-Wendroff method based on explicit CD stencils for D3 and D4 terms
included in CDE, plotted in the $(N_{c},kh)$-plane for the representative
Peclet numbers $Pe=0.01$ and $0.02$. Regions of numerical instability are as
marked in the panels. Figure 4: The ratio of numerical and physical phase
speeds $\frac{c_{num}}{c_{phys}}$ and the ratio of numerical and the physical
group velocity $\frac{v_{g,num}}{v_{g,phys}}$ for the Lax-Wendroff method
based on explicit CD stencils for D3 and D4 terms included in CDE plotted in
the $(N_{c},kh)$-plane for the representative Peclet numbers $Pe=0.01$ and
$0.02$.

The relative advantages of using the LW method including the $D^{3}u$ and
$D^{4}u$ terms, based on explicit central differences in the CDE are noted in
Figs. 3 and 4, with property charts plotted in the $(N_{c},kh)$-plane for the
same two $Pe$ values used in Figs. 1 and 2. As noted earlier for the method
without $D^{3}u$ and $D^{4}u$ terms, one notes similar behavior for
$\alpha_{num}$, i.e. lower $\alpha_{num}$ occurring at lower $N_{c}$ values,
whereas higher $\alpha_{num}$ occurs for an adjacent range of $N_{c}$, as
compared to the physical value of $\alpha$. Comparing Figs. 1 and 3, one notes
the LW method with $D^{3}u$ and $D^{4}u$ terms to have an increased stable
region, as compared to the LW method applied only for the convection term. The
full LW scheme also has better numerical diffusion at all $kh$. This is a
consequence of retaining the $D^{4}u$ term that moderates the added/reduced
numerical diffusion across the $N_{c}$ range. The role of positive fourth
diffusion term in Eq. (25) can be understood by noting that
$\partial^{4}u\over\partial x^{4}$ term varies as $k^{4}$ in the spectral
plane, thereby reducing the added numerical diffusion by the introduced second
derivative term in the Taylor series expansion of the convection term.

Comparing the numerical phase speed and group velocity contours between the
convection alone and full LW method based on CD schemes from Figs. 2 and 4,
one notes the latter to have properties extended slightly in the $N_{c}$
direction due to increased stability. Apart from this, only marginal
differences are noted in terms of accuracy between the two strategies for
dispersion errors.

We note that both strategies with CD scheme are less suitable for high
accuracy solutions of the CDE due to present errors at all $kh$ noted for
$\alpha_{num}$, $c_{num}$ and $v_{g,num}$. This implies that these strategies
are strictly less suitable for DNS, due to noted one-to-one correspondence
between the 2D linear CDE and 2D Navier-Stokes equations demonstrated in 44.
However, it is feasible to quantify the parameters for the LW method on the
CDE and look for some combinations of $N_{c}$ and $kh$ ranges, which is shown
next for the various differences appearing in the LW method used for the 2D
CDE.

##  5 Quantifying LW method using 2D CDE for LES

Having established that the convection-only LW discretization of the 1D CDE is
comparable in accuracy to the LW method retaining $D^{3}$, $D^{4}$ terms, in
this section, the properties of the same scheme are further analyzed for the
2D CDE. Previous research work 44 has demonstrated a one-to-one correspondence
of the linear 2D CDE with the 2D Navier-Stokes equation. Hence, its role to
analyze schemes for the Navier-Stokes equation cannot be over-emphasized, and
is adopted here for the same purpose.

The essentials of the GSA of 2D CDE required to analyze the 2D LW scheme are
given below. Interested readers can consult 44 for details and discussion
regarding the expressions obtained from GSA.

The 2D CDE is given as,

$\frac{\partial u}{\partial t}+c_{x}\frac{\partial u}{\partial
x}+c_{y}\frac{\partial u}{\partial
y}=\alpha\left(\frac{\partial^{2}u}{\partial
x^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\right)$ (45)

For the 2D CDE, the LW scheme applied to the convection terms can be obtained
as

$u(t+\Delta t)=u(t)-\Delta t\left(c_{x}\frac{\partial u}{\partial
x}+c_{y}\frac{\partial u}{\partial y}\right)+\alpha\Delta
t\left(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial
y^{2}}\right)+\frac{(\Delta
t)^{2}}{2}\left(c_{x}^{2}\frac{\partial^{2}u}{\partial
x^{2}}+2c_{x}c_{y}\frac{\partial^{2}u}{\partial x\partial
y}+c_{y}^{2}\frac{\partial^{2}u}{\partial y^{2}}\right)$ (46)

To perform GSA of Eq. (45), the unknown is expressed in the hybrid-spectral
plane as,

$u(x,y,t)=\iint\hat{U}(k_{x},k_{y},t)e^{i(k_{x}x+k_{y}y)}dk_{x}dk_{y}$ (47)

where $\hat{U}$ is the Fourier-Laplace amplitude and $k_{x}$, $k_{y}$ are the
wavenumber components in the $x$\- and $y$-directions, respectively.
Substituting the expression for $u$ in the governing equation one obtains,

$\hat{U}(k_{x},k_{y},t)=\hat{U}_{0}(k_{x},k_{y})\;e^{-\alpha(k_{x}^{2}+k_{y}^{2})t}e^{-i(k_{x}c_{x}+k_{y}c_{y})t}$
(48)

where $u(x,y,0)=\iint\hat{U}_{0}(k_{x},k_{y})e^{i(k_{x}x+k_{y}y)}dk_{x}dk_{y}$
is the initial solution.

The physical dispersion relation, an important property that must be obeyed by
numerical schemes to minimize phase and dispersion errors, is obtained by
expressing $u$ in the full spectral space as

$\omega_{0}=c_{x}k_{x}+c_{y}k_{y}-i\alpha(k_{x}^{2}+k_{y}^{2})$ (49)

From this relation the physical phase speed is obtained as,

$c_{phys}=\frac{\omega_{0}}{\sqrt{k_{x}^{2}+k_{y}^{2}}}=\frac{c_{x}k_{x}+c_{y}k_{y}}{\sqrt{k_{x}^{2}+k_{y}^{2}}}$
(50)

and the physical group velocity components are obtained as,

$v_{gx,phys}=\frac{\partial\omega_{0}}{\partial k_{x}}=c_{x}$ (51)
$v_{gy,phys}=\frac{\partial\omega_{0}}{\partial k_{y}}=c_{y}$ (52)

The physical amplification factor for the CDE is

$\displaystyle
G_{phys}=e^{-[Pe_{x}(k_{x}h_{x})^{2}+Pe_{y}(k_{y}h_{y})^{2}]}e^{-i[N_{cx}k_{x}h_{x}+N_{cy}k_{y}h_{y}]}$
(53)

where $\Delta t$ is the discrete time-step and $h_{x}$ and $h_{y}$ are the
grid spacings in $x$\- and $y$-directions, respectively. The CFL and Peclet
numbers in 2D are: ${N}_{cx}=\frac{c_{x}\Delta t}{h_{x}}$;
${N}_{cy}=\frac{c_{y}\Delta t}{h_{y}}$; $Pe_{x}=\frac{\alpha\Delta
t}{{h_{x}}^{2}}$; $Pe_{y}=\frac{\alpha\Delta t}{{h_{y}}^{2}}$.

For a numerical scheme the corresponding numerical dispersion relation is
obtained by analogy as,

$\omega_{num}=\left(\sqrt{{k_{x}}^{2}+{k_{y}}^{2}}\right)c_{num}-i\alpha_{num}(k_{x}^{2}+k_{y}^{2})$
(54)

where the subscript $num$ denotes that the corresponding property is for the
numerical scheme.

From this numerical dispersion relation, one obtains numerical amplification
factor, $G_{num}$, as

$\displaystyle G_{num}=e^{-i\omega_{num}\Delta
t}=e^{-\alpha_{num}({k_{x}}^{2}+{k_{y}}^{2})\Delta
t}\,\,e^{-i\left(\sqrt{{k_{x}}^{2}+{k_{y}}^{2}}\right)c_{num}\Delta t}$ (55)

Employing second order, explicit central difference schemes for all the
derivatives in Eq. (46), the numerical amplification factor for the 2D LW
scheme is

$\begin{split}G_{num}=&1-iN_{cx}\sin(k_{x}h_{x})-iN_{cy}\sin(k_{y}h_{y})+(N_{cx}^{2}+2Pe_{x})[\cos(k_{x}h_{x})-1]\\\\[6.45831pt]
&+(N_{cy}^{2}+2Pe_{y})[\cos(k_{y}h_{y})-1]-N_{cx}N_{cy}\sin(k_{x}h_{x})\sin(k_{y}h_{y})\end{split}$
(56)

For accuracy, we require $G_{num}/G_{phys}\approx 1$, with $|G_{num}|<1$. From
$G_{num}$, numerical phase speed is obtained as

$\frac{c_{num}}{c_{phys}}=-\left[\frac{1}{N_{cx}(k_{x}h_{x})+N_{cy}(k_{y}h_{y})}\right]{\tan}^{-1}\left[\frac{(G_{num})_{Img}}{(G_{num})_{Real}}\right]$
(57)

where $c_{phys}$ is the physical phase speed.

The numerical group velocity components are then obtained as

$\frac{v_{gx,num}}{c_{x}}=\frac{1}{N_{cx}}\frac{\partial\beta_{num}}{\partial(k_{x}h_{x})}$
(58)

$\frac{v_{gy,num}}{c_{y}}=\frac{1}{N_{cy}}\frac{\partial\beta_{num}}{\partial(k_{y}h_{y})}$
(59)

with
$\tan(\beta_{num})=-\left[\frac{(G_{num})_{Img}}{(G_{num})_{Real}}\right]$.

The numerical diffusion coefficient $\alpha_{num}$ is given as

$\frac{\alpha_{num}}{\alpha}=-\frac{\ln{|G_{num}|}}{[{Pe}_{x}(k_{x}h_{x})^{2}+{Pe}_{y}(k_{y}h_{y})^{2}]}$
(60)

The ratio $\frac{\alpha_{num}}{\alpha}$ determines the numerical diffusion
offered by the scheme. If it is unity, then the scheme models the physical
diffusion exactly. If the ratio is greater or lesser than unity, then the
numerical diffusion is higher or lower than the physical diffusion. Negative
value(s) denote anti-diffusion which leads to numerical instability. For
accuracy of solution, $\frac{\alpha_{num}}{\alpha}$, $\frac{c_{num}}{c}$,
$\frac{v_{gx,num}}{c_{x}}$ and $\frac{v_{gy,num}}{c_{y}}$ should be equal to
unity.

An interesting observation can be made by comparing the 2D LW scheme with its
1D counterpart. One notes the presence of an additional cross-derivative term
$u_{xy}$ for the 2D case which is absent for the 1D case. From spectral
analysis, this cross-derivative term introduces additional numerical
dissipation when both $k_{x}$ and $k_{y}$ are of the same sign and anti-
diffusion when the wavenumbers are of the opposite sign. This behavior of the
2D LW scheme has not been reported before. Hence, it becomes necessary to
evaluate the properties of the scheme in order to assess it for simulating
fluid flows.

In Fig. 5, the ratio of the numerical amplification factor to the physical
amplification factor $|G_{num}/G_{phys}|$ is plotted for Peclet numbers
$Pe_{x}=Pe_{y}=0.01$ and CFL numbers $N_{cx}=N_{cy}=0.08$ to $0.11$. The
property charts correspond to the simulation parameters for a uniform grid
with equal spacing in $x$ and $y$\- directions ($h_{x}=h_{y}$) and the
solution propagating at an angle of $45^{o}$ ($\tan(c_{y}/c_{x})=1$). For
these simulation parameters, anti-diffusion is absent (not shown in the
figure) and hence, the numerical simulations will be stable, i.e. $G_{num}\leq
1$. Only two contour values are plotted which correspond to errors of
$10^{-4}$ and $10^{-6}$ from the ideal value of $|G_{num}/G_{phys}|=1$. Thus,
values of $0.9999$ and $1.0001$ correspond to the tolerance limit of error of
$10^{-4}$ while the values $0.999999$ and $1.000001$ denote a lower error of
$10^{-6}$, respectively. Comparing the contours for different CFL values, one
notes that a maximum extent of region with low error is obtained for
$N_{cx}=N_{cy}=0.09$ and hence, this is used as optimal values for the
considered $Pe$ values.

Figure 5: Ratio of numerical to physical amplification factor
($|G_{num}/G_{phys}|$) of the Lax-Wendroff method based on explicit CD
stencils for the convection terms in 2D CDE, plotted in the
$(k_{x}h_{x},k_{y}h_{y})$-plane for the representative Peclet and CFL numbers.
Properties are shown for a uniform grid with equal spacing in both directions
and the wave propagating at $45^{o}$.

The importance of determining the optimal values for any numerical scheme
cannot be underscored as they will help the research community to utilize the
scheme in an efficient manner. By identifying optimal $(N_{cx},N_{cy})$ for
given Peclet numbers $(Pe_{x},Pe_{y})$, timestep $\Delta t$ or grid spacings
$h_{x},h_{y}$ can be fixed. One can fully appreciate the vital and critical
role played by GSA due to its accurate characterization of numerical schemes
for solving specific governing equations. For the present case, it should be
noted that only one free variable exists $(N_{cx})$ for determining the
optimal conditions as the numerical setup involves an equispaced, uniform grid
with signal propagating at $45^{o}$.

The optimal simulation parameters, evaluated from the analyses based on
acceptable errors in $|G_{num}/G_{phys}|$ contours, are reinforced by plotting
the ratio of numerical to physical diffusion $\alpha_{num}/\alpha$ for the LW
scheme in Fig. 6 for the same parameters as before. One notes that a maximum
resolution for performing fine simulations is obtained for
$Pe_{x}=Pe_{y}=0.01$ when $N_{cx}=N_{cy}=0.09$. An important observation can
be immediately drawn from the figure that the numerical diffusion increases in
strength as the CFL value is increased for fixed Peclet values. This is
expected as the LW scheme introduces additional numerical dissipation as
established in Eq. (46). Further, with increase in CFL values, the
effectiveness of numerical diffusion increases for higher wavenumbers which
can be an advantage in controlling high wavenumber numerical instabilities
such as aliasing, arising due to a lack of resolution.

Figure 6: Ratio of numerical to physical diffusion coefficient
($\alpha_{num}/\alpha$) for the Lax-Wendroff method based on explicit CD
stencils for the convection terms in 2D CDE, plotted in the
$(k_{x}h_{x},k_{y}h_{y})$-plane for the representative Peclet and CFL numbers.
Properties are shown for a uniform grid with equal spacing in both directions
and the wave propagating at $45^{o}$.

In Fig. 7, the regions representing the two tolerance levels are shown for the
optimal CFL values for Peclet numbers of $0.01$ and $0.02$, respectively. The
boundaries of the regions, denoted by OABC, are chosen so as to maximize the
resolutions $k_{x}$, $k_{y}$ satisfying the error constraint. The top frames
correspond to the higher error tolerances of $10^{-4}$ with the regions marked
by solid blue lines in the form of a rectangle and the bottom frames are for
the lower error tolerance of $10^{-6}$ with the regions marked by dashed blue
rectangles, respectively. In the present work, the former condition is
representative of coarse simulations while the latter case denotes fine
simulations such as LES/DNS. This is the rationale behind the choice of the
two tolerance values. As the tolerance level is reduced, the region satisfying
the error criterion shrinks i.e one obtains lower errors by moving towards the
continuum limit. It is also interesting to note that the regions for coarse
and fine simulations decrease with increasing $Pe$ values. In determining the
optimal limits of CFL numbers and the accuracy of the scheme, we have only
considered the error in representing the physical amplification factor.
However, error is also contributed due to dispersion errors as numerical phase
speed $c_{num}$ need not be equal to the physical phase speed $c_{phys}$.

Figure 7: Zones of coarse (top panels) and fine simulations (bottom panels)
for the Lax-Wendroff method for explicit CD schemes used for the convection
term in 2D CDE, for the indicated Peclet numbers and optimal CFL values of
$0.09$. The solid and dashed blue rectangles denote regions with errors
$|1-|G_{num}/G_{phys}||\leq 10^{-4}$ and $|1-|G_{num}/G_{phys}||\leq 10^{-6}$,
respectively.

In the analyses described earlier, we have evaluated the performance of the LW
scheme for a small operating range of Pe and $N_{c}$ values. Therefore, it
would only be natural and also essential to quantify the scheme for a wider
range of simulation parameters, i.e. determine the accuracy limits of the
scheme. It is important to note that in the process of determining these
limits, the information on accuracy should also be simultaneously available to
the researcher. A direct presentation of this information is a complicated
task due to the multidimensional ($\geq 3D$) nature of the data. In this
regard, we present a simple and a tractable approach that can be followed to
achieve the desired goal. First, $\alpha_{num}/\alpha$ property is extracted
along a diagonal line $(k_{x}h_{x}=k_{y}h_{y})$ for each $N_{c}$ value and for
a fixed Pe number. The resulting data can be then plotted as contours in the
$[N_{cx}(=N_{cy}),k_{x}h_{x}(=k_{y}h_{y})]$-plane for a fixed Pe value.

This approach is demonstrated in Fig. 8 for the LW scheme for Pe values $0.01$
and $0.02$, respectively. One notes a striking similarity of this figure with
the 1D results (Figs. 1, 3) in the previous section. From the results, the
stability limits are determined by noting the CFL values ($N_{cx,crit}$) for
which anti-diffusion appears. This is noted to be $0.4515$ and $0.4746$ for
$Pe_{x}=0.01$ and $0.02$, respectively. Thus, increasing Pe value has the
effect of increasing $N_{cx,crit}$. It is also interesting to note that the
results can also be used to determine the optimal $N_{cx}$ values by noting
the location(s) at which $\alpha_{num}/\alpha=1$ contour is vertical. This is
determined from the figure as $N_{cx}=0.089$ and $0.14$ for $Pe_{x}=0.01$ and
$0.02$, respectively.

Figure 8: Ratio of numerical to physical diffusion coefficient
($\alpha_{num}/\alpha$) for the Lax-Wendroff method based on explicit CD
stencils for the convection terms in 2D CDE, plotted in the
$(N_{cx},k_{x}h_{x})$-plane for the representative Peclet numbers. Properties
are shown for a uniform grid with equal spacing in both directions and the
wave propagating at $45^{o}$.

Apart from the error of representing the numerical amplification factor
accurately, simulations suffer from additional source(s) of error. In 14, the
correct error dynamics equation is derived for the 1D convection equation for
the first time using GSA which reveals all the sources of error in computing.
Another important contribution of 14 is to dispel the notion that the
evolution equation for the error has the same form as the governing equation.
In 55, the governing equation for error evolution is derived for the 1D
diffusion equation using GSA showing the contribution to error due to
numerical diffusion being different from physical diffusion. Following the
same approach, the error dynamics equation for the 2D CDE is presented for the
first time as

$\begin{split}e_{t}+c_{x}e_{x}+c_{y}e_{y}&=\alpha(e_{xx}+e_{yy})\\\
&+\iint\limits_{-k_{max}}^{k_{max}}(\alpha_{num}-\alpha)\;(k_{x}^{2}+k_{y}^{2})\;|G_{num}|^{t^{n}/\Delta
t}\;\hat{U}_{0}(k_{x},k_{y})\;e^{i[k_{x}x+k_{y}y-(\sqrt{k_{x}^{2}+k_{y}^{2}})c_{num}t^{n}]}\,dk_{x}dk_{y}\\\
&\iint\limits_{-k_{max}}^{k_{max}}i(\sqrt{k_{x}^{2}+k_{y}^{2}})(c_{num}-c_{phys})\;|G_{num}|^{t^{n}/\Delta
t}\;\hat{U}_{0}(k_{x},k_{y})\;e^{i[k_{x}x+k_{y}y-(\sqrt{k_{x}^{2}+k_{y}^{2}})c_{num}t^{n}]}\,dk_{x}dk_{y}\end{split}$
(61)

Here, $e$ denotes the error and is defined by the difference between exact
($u_{phys}$) and numerical solution ($u_{num}$). From the above equation one
notes contribution to error due to incorrect numerical diffusion and phase
speed, respectively. The former contribution is also denoted as error due to
spurious numerical diffusion/amplification and the latter is the dispersive
error. Comparing the error dynamics equation with the MDE for 2D LW method for
CDE given by Eq. (46), the cross-derivative term is noted to be missing which
leads one to believe that the error dynamics equation is incomplete. However,
this assertion is incorrect as the missing term is already accounted by
$G_{num}$ and $c_{num}$ terms which incorporate the effects of all the terms
in MDE including the cross-derivative term.

The error properties are quantified for the Lax-Wendroff scheme applied to the
convection term in the 2D CDE in Tables 1 and 2 for two Peclet numbers $0.01$
and $0.02$, respectively. These properties are obtained from the rectangular
regions OABC identified in Fig. 7. The tables present the maximum resolution,
dispersion and dissipation errors exhibited by the numerical scheme. From
these results, it is apparent that the maximum resolution satisfying the error
constraint is higher for the coarse simulation case when compared to the fine
simulation case. We note that the dispersion errors are also maximum for the
coarse simulation case. This can be attributed to the poor spectral resolution
of the CD2 stencils employed for the Lax-Wendroff scheme. A comparison of the
error properties for the two Peclet number cases reveals significant
improvement in dispersion and group velocity error properties for the higher
Pe value. However, a substantial difference is noted in the maximum resolution
for the coarse simulation cases with the higher Pe case showing a reduced
resolution.

Lax-Wendroff | Simulation | Maximum | Maximum |  $\alpha_{num}$/ $\alpha$ |  $c_{num}$/$c_{phys}$ |  $v_{gx,num}$/$c_{x}$ |  $v_{gy,num}$/$c_{y}$
---|---|---|---|---|---|---|---
Method for | type | resolution | resolution | range | range | range | range
|  | $(k_{x}h_{x})_{max}$ | $(k_{y}h_{y})_{max}$ |  |  |  |
Convection term | Coarse | 1.4 | 1.4 | 0.994641-1.00691 | 0.72529-1 | 0.205689-1.02329 | 0.211452-1.02393
Convection term | Fine | 0.34375 | 0.34375 | 1-1.00125 | 0.95308-1 | 0.95507-1.00081 | 0.953638-1.00082

Table 1: Quantification of DRP properties for coarse
($|1-|G_{num}/G_{phys}||\leq 10^{-4}$) and fine simulations
($|1-|G_{num}/G_{phys}||\leq 10^{-6}$) for 2D CDE using the Lax-Wendroff
method based on CD schemes for $Pe_{x}=Pe_{y}=0.01$ and $N_{cx}=N_{cy}=0.09$.

Lax-Wendroff | Simulation | Maximum | Maximum |  $\alpha_{num}$/ $\alpha$ |  $c_{num}$/$c$ |  $v_{gx,num}$/$c_{x}$ |  $v_{gy,num}$/$c_{y}$
---|---|---|---|---|---|---|---
Method for | type | resolution | resolution | range | range | range | range
|  | $(k_{x}h_{x})_{max}$ | $(k_{y}h_{y})_{max}$ |  |  |  |
Convection term | Coarse | 0.64264 | 0.64264 | 0.98811-1 | 0.940531-1 | 0.82691-1.00961 | 0.829639-1.00994
Convection term | Fine | 0.20498 | 0.20498 | 0.998821-1 | 0.994327-1 | 0.983572-1.00076 | 0.98273-1.0008

Table 2: Quantification of DRP properties for coarse
($|1-|G_{num}/G_{phys}||\leq 10^{-4}$) and fine simulations
($|1-|G_{num}/G_{phys}||\leq 10^{-6}$) for 2D CDE using the Lax-Wendroff
method based on CD schemes for $Pe_{x}=Pe_{y}=0.02$ and $N_{cx}=N_{cy}=0.09$.

From the presented property charts and error quantification for the Lax-
Wendroff scheme for 2D CDE it can be inferred that accurate computation of
Navier-Stokes equation is possible provided one employs timestep as dictated
by the optimal CFL values and a finely resolved grid such that all the
relevant flow scales are within the maximum resolution dictated by the fine
simulation error tolerances. This is demonstrated in the next section by
solving the steady and unsteady flow inside a square lid driven cavity (LDC).

##  6 LES of 2D NSE by LW Method

In this section, we demonstrate the utility of the Lax-Wendroff scheme in
accurately solving the 2D incompressible Navier-Stokes equation for the flow
inside a square LDC for post-critical Reynolds number. This problem is chosen
specifically due to its simple, unambiguous boundary conditions and the
availability of benchmark solutions for unsteady case 48, 49. The post-
critical Reynolds number considered in the present study is 10,000. A concise
summary of the governing equations and the methodology adopted for the
numerical solution is discussed next followed by numerical simulations using
the 2D LW scheme.

###  6.1 Governing Equations and Solution Methodology

The governing 2D incompressible NSE is solved in the streamfunction-vorticity
($\psi$-$\omega$) formulation. This formulation is adopted due to its
advantages in satisfying mass conservation automatically in the computational
domain. In this approach, one solves two equations- a Poisson equation for
$\psi$ and a transport equation for $\omega$, which are given below.

$\nabla^{2}\psi=-\omega$ (62) $\frac{\partial\omega}{\partial
t}+(\vec{V}\cdot\vec{\nabla})\omega=\frac{1}{Re}\nabla^{2}\omega$ (63)

These equations are given in their non-dimensional form with $Re$ denoting the
reference Reynolds number based on the side of the cavity and the speed of the
upper lid which moves from left to right. Velocity vector
$\vec{V}=u\hat{i}+v\hat{j}$ is computed from the stream-function by
$\vec{V}=\vec{\nabla}\times\vec{\psi}$, where
$\vec{\psi}=[0\;0\;\psi]{{}^{T}}$.

The governing equations are solved in the following manner. The stream
function equation, Eq. (62) is solved first by using CD2 scheme for
discretization and the BiCGSTAB iterative method 50 for eventual solution of
the discrete Poisson equation. After this step, the vorticity at the
boundaries is computed from its definition: $\omega=-\nabla^{2}\psi$. Next,
vorticity at new time instant is updated by solving Eq. (63), wherein the CD2
Lax-Wendroff scheme applied to convection terms is employed. The Lax-Wendroff
scheme for the 2D VTE is given in Eq. (64) for reference. The process is
repeated until the simulation time reaches a maximum user specified value or
when the flow reaches a steady state. In performing these steps, one should
also note that a constant value, Dirichlet boundary condition is imposed on
$\psi$ ($\psi=\psi_{0}$) in order to satisfy the non-penetrative condition at
the solid walls.

$\begin{split}\omega_{ij}^{n+1}=\omega_{ij}^{n}+\Delta
t\left(-u\frac{\partial\omega}{\partial x}-v\frac{\partial\omega}{\partial
y}\right)_{ij}^{n}+\frac{\Delta t^{2}}{2}&\left[\left(u\frac{\partial
u}{\partial x}+v\frac{\partial u}{\partial
y}\right)\frac{\partial\omega}{\partial
x}+u^{2}\frac{\partial^{2}\omega}{\partial x^{2}}+\left(u\frac{\partial
v}{\partial x}+v\frac{\partial v}{\partial
y}\right)\frac{\partial\omega}{\partial y}\right.\\\\[6.45831pt]
&\left.+2uv\frac{\partial^{2}\omega}{\partial x\partial
y}+v^{2}\frac{\partial^{2}\omega}{\partial
y^{2}}\right]_{ij}^{n}+\frac{1}{Re}\nabla^{2}\omega_{ij}^{n}\end{split}$ (64)

###  6.2 Numerical Solution at LDC at Post-critical Re=10,000

The LDC problem is solved here for a post-critical Re of 10,000 using the
explicit CD based Lax-Wendroff scheme. The present case is chosen as it is an
excellent test case for benchmarking the accuracy of numerical schemes/codes.
Accurate solutions for this Re must display transient, triangular vortex
structure during the flow evolution 48, 49. Hence, its capture serves as a
direct confirmation of accuracy for unsteady flows. The numerical solution is
obtained for a uniform, equispaced grid of $1112\times 1112$ with a timestep
of $\Delta t=8.1\times 10^{-5}$. These simulation parameters result in
$Pe_{x}=Pe_{y}=0.01$ and CFL number based on upper lid velocity $N_{cx}=0.09$,
respectively and correspond to the optimal conditions determined for the
scheme.

In Fig. 9, vorticity contours are plotted at the indicated time instants
showing the complex vortex dynamics of the flow. At early times, one notes the
formation of a pentagonal vortical structure in the core which subsequently
evolves into a triangular structure due to the shearing action. The triangular
vortex is observed which is surrounded by gyrating satellites which rotate
along with it. The triangular vortex structure shrinks in size at much later
times thus displaying its transient nature. This evolution provides an
excellent match with the observations in 48, 49 which was used for high
accuracy combined compact schemes. The present results, therefore,
demonstrates the potential of the Lax-Wendroff scheme for solving unsteady
flows and further consolidates the utility of GSA in determining the optimal
limits for simulation parameters for realizing a desired level of accuracy.

Figure 9: Evolution of vorticity field showing the transient, triangular core
vortex for the LDC problem for Re=10,000 computed by the explicit central
difference based Lax-Wendroff method applied to convection terms of 2D
incompressible NSE.

##  7 Summary and Conclusions

In the present paper, the Lax-Wendroff (LW) method based on explicit central
differences is comprehensively analyzed using global spectral analysis (GSA)
for the model 1D and 2D linear convection-diffusion equations (CDE) in order
to obtain optimal simulation parameters for performing accurate implicit LES
(ILES). The LW method is one of the earliest approaches in developing second
order in time DRP schemes via the modified differential equation (MDE), where
higher order time derivatives are converted to spatial derivatives using the
governing equation. While the LW scheme was analyzed recently using GSA,
important questions regarding its variants, analysis for higher than one
dimensional problems and optimal parameters for accuracy are addressed here in
brief. The optimal simulation parameters presented here ensure accurate
simulation of CDE and thereby Navier-Stokes equations, ILES to be specific.
This is achieved by a rigorous quantification of the numerical properties viz.
the resolution, numerical diffusion, dispersion and signal propagation speed
for the LW method for 2D CDE.

A consequence of the application of the GSA for the 1D and 2D CDE is that the
diffusion and solution propagation properties for numerical schemes become
space-time dependent. Hence, the numerical schemes have corresponding
wavenumber dependent coefficients such as numerical diffusion coefficient
$\alpha_{num}$, numerical phase speed $c_{num}$ and numerical group velocity
$v_{g,num}$, which determine the evolution of the solution. To accurately
solve the governing equations, these properties must match their physical
counterparts, and thereby these define dispersion relation preserving (DRP)
schemes.

Formulating the MDE for the 1D CDE, two variants of the LW method are
developed based on the treatment of convection and convection-diffusion terms
i) as applied only to the convection term (3rd and 4th derivatives are set to
zero in Eq. (25)) and ii) applied to all the terms given in Eq. (25),
respectively. The first approach is popularized and practiced by many
practitioners without adequate numerical analysis, while the second approach
is the full LW scheme in the present context. One of the interesting
observations noted in the application of the full LW procedure for 1D CDE is
the inclusion of a third and fourth order derivative terms along with an added
second order derivative term. The latter is noted as the reason for the LW
method to be inconsistent for the 1D convection equation. Presence of the
third and fourth derivatives introduce the effect of dispersion and weak anti-
diffusion when using the LW method with the latter moderating the added
numerical diffusion to the already existing physical term.

GSA shows marginal benefits offered by the full LW scheme over the other
variant. Analysis of two values of Peclet numbers $0.01$ and $0.02$, shows the
full LW method to possess an increased stability region over the other
variant. The analysis also reveals parameter combinations in the
$(N_{c},kh)$-plane for both schemes where the overall numerical diffusion is
equal to the physical diffusion despite the added second and fourth order
numerical diffusion terms. The Lax-Wendroff method applied to convection term
alone shows better performance because of the near vertical contour line
${\alpha_{num}\over\alpha}=1$, offering more resolution over the full scheme.
Furthermore, the full LW scheme is computationally more expensive due to the
$3^{rd}$ and $4^{th}$ order derivative terms. This provides the necessary
justification for adopting the LW method based on its application to the
convection term. Hence, the LW scheme applied to convection only terms is
analyzed for 2D CDE in order to determine optimal parameters for ILES and
Navier-Stokes simulations.

The LW method is evaluated for solution of 2D CDE using GSA due to a one-to-
one correspondence between the model equation and the Navier-Stokes equations
44. An additional cross-derivative term is noted for the scheme which is
reported for the first time and whose effect is to cause asymmetry in the
numerical properties. Rigorous quantification/ evaluation of the DRP region is
performed to determine the best numerical parameters viz. a time step/ grid
resolution for accurate ILES. This is evaluated by prescribing tolerance
limits of $10^{-4}$ and $10^{-6}$ in $|1-|{G_{num}\over G_{phys}}||$ i.e.
$|1-|{G_{num}\over G_{phys}}||\leq\epsilon$, with $\epsilon$ as the tolerance.
The latter tolerance limit can be considered as representative of ILES and
unresolved DNS scenarios and such cases are termed as fine simulations.
Optimal CFL values are determined for representative Pe values for which the
scheme performs the best as noted from Figs. 5 and 6. In determining these
limits, only the error in numerical amplification factor is considered.
Analysis of the complete sources of error for the 2D CDE given by Eq. (61)
also shows contribution from the dispersion property. A complete assessment of
the performance of the method with respect to these errors are presented in
Tables 1 and 2. It is also noted that as CFL or Peclet number increases, the
scheme becomes more stable.

The optimal parameters obtained for the LW scheme are corroborated by solving
2D Navier-Stokes equations for the square lid driven cavity problem for post-
critical Reynolds number of $10,000$. For the post-critical case optimal CFL
conditions are used as the time evolution is important. The simulation shows a
very good agreement with the benchmark results for the unsteady case capturing
the transient triangular vortex which is considered as a stringent test for
validating high accuracy solution methods. The present work demonstrates
accurate analyses of the LW method as corroborated by Navier-Stokes simulation
thus highlighting the potential of GSA approach.

## Acknowledgments

### Author contributions

All authors have contributed equally to the data generation, analysis, writing
and preparation of the manuscript.

### Financial disclosure

None reported.

### Conflict of interest

The authors declare no potential conflict of interests.

## Supporting information

None available.

## References

  * 1 Courant, R., Friedrich, K. O. and Lewy, H. On the partial differential equations of mathematical physics. IBM J. , 215-234 (1967)
  * 2 Ames, W. F. Numerical Methods for Partial Differential Equation, 2ndS edn., Academic Press, New York, USA (1977)
  * 3 Sengupta, T. K. High Accuracy Computing Method: Fluid Flows and Wave Phenomena, Cambridge Univ. Press, New York, USA (2013)
  * 4 Richardson, L. F. Weather Prediction by Numerical Process, Cambridge Univ. Press, Cambridge, UK (1922) (Reprinted by Dover (1965))
  * 5 von Neumann, J. and Richtmyer, R.D. On the numerical solution of partial differential equations of parabolic type. Los Alamos Rept., Series A LA-657, 1–17 (1947).
  * 6 Charney, J. G., Fjörtoft, R. and von Neumann, J. Numerical integration of the barotropic vorticity equation. Tellus, 2(4), 237-254 (1950)
  * 7 Morton, K.W. and Mayers, D.F. Numerical Solution of Partial Differential Equations. $2^{nd}$ Ed., Cambridge Univ. Press, U.K. (2005).
  * 8 Zingg, D. W. Comparison of high-accuracy finite-difference schemes for linear wave propagation. SIAM J. Sci. Comp. 22(2), 476-502 (2000)
  * 9 Haltiner, G. J. and Williams, R. T. Numerical Prediction and Dynamic Meteorology. $2^{nd}$ edn., John Wiley $\&$ Sons., New York, USA (1980)
  * 10 Durran, D. R. Numerical Methods for Wave Equations in Geophysical Fluid Dynamics, Springer Verlag, New York, USA (1999)
  * 11 Zingg, D. W., Lomax, H. and Jurgens, H. High accuracy finite difference schemes for linear wave propagation. SIAM. J. Sci. Comp., 17(2), 328–346 (1996).
  * 12 Hirsch, C. Numerical Computation of Internal and External Flows. Vols. I and II. Computational Methods for Inviscid and Viscous Flows. Wiley, Chichester, UK (1990)
  * 13 Sengupta, T.K. and Dipankar, A., A comparative study of time advancement methods for solving Navier–Stokes equation. J. Sci. Comp., 21(2), 225–250 (2004).
  * 14 Sengupta, T.K., Dipankar, A. and Sagaut, P. Error dynamics: beyond von Neumann analysis. J. Comp. Phys. 226, 1211–1218 (2007).
  * 15 Milne, W.E. Numerical Solution of Differential Equations. Wiley, New York, USA (1953).
  * 16 Du Fort, E.C. and Frankel, S.P. Stability conditions in the numerical treatment of parabolic differential equations. Math. Table NRC, Washington, 7B5, 135–153 (1953).
  * 17 Warming, R.F. and Hyett, B.J. The modified equation approach to the stability and accuracy analysis of finite difference methods, J. Comp. Phys. 14, 159–179 (1974).
  * 18 Li, J. and Yang, Z. The von Neumann analysis and modified equation approach for finite difference schemes, App. Math. Comp., 225, 610-620 (2013).
  * 19 Garabedian, P.R. Estimation of the relaxation factor of small mesh sizes. Math. Tables Aids Comput., 10, 1983–185 (1956).
  * 20 Harten, A., Hyman, J.M. and Lax, P.D., On the finite difference approximations and entropy conditions for shocks. Commun. Pure Appl. Math. 29, 297-322 (1976).
  * 21 Griffiths, D.F. and Sanz-Serna, J.M. On the scope of the method of modified equations. SIAM J. Sci. Stat. Comput., 7, 994–1008 (1986).
  * 22 Lax, P.D. and Wendroff, B. System of conservations laws. Comm. Pure App. Maths., 13, 217–237 (1960).
  * 23 Winnicki, I., Jasinski, J. and Pietrek, S., New approach to Lax-Wendroff modified differential equation for linear and nonlinear advection. Numer. Methods Partial Differential Eq., 1–30 (2019).
  * 24 Shokin, Y.I. The Method of Differential Approximation. Springer–Verlag, Berlin, Germany (1983).
  * 25 Yanenko, N.N., Fedotova, Z.I., Tusheva, L.A. and Shokin, Yu.I. Classification of difference schemes of gas dynamics by the method of differential approximation–I. Comput. Fluids, 11(3), 187–206 (1983).
  * 26 Vichnevetsky, R. and Bowles, J. B. Fourier Analysis of Numerical Approximations of Hyperbolic Equations. SIAM Stud. App. Math., 5, Philadelphia, USA (1982).
  * 27 Lele, S.K. Compact finite difference schemes with spectral like resolution. J. Comput. Phys., 103, 16–42 (1992).
  * 28 Sengupta, T. K., Ganeriwal, G. and De, S. Analysis of central and upwind compact schemes. J. Comput. Phys., 192(2), 677–694 (2003).
  * 29 Sengupta, T. K. Fundamentals of Computational Fluid Dynamics. Universities Press, Hyderabad, India (2004).
  * 30 Carpenter, M.H., Gottlieb, D. and Abarbanel, S. The stability of numerical boundary treatments for compact high-order finite difference schemes. J. Computational Physics, 108, 272-295 (1993).
  * 31 Hu, F. Q., Hussaini, M. Y. and Manthey, J. L. Low-dissipation and low-dispersion Runge-Kutta schemes for computational acoustics. J. Computational Physics, 124, 177-191 (1996).
  * 32 Zhong, X. High-order finite difference schemes for numerical simulation of hypersonic boundary-layer transition. J. Computational Physics, 144, 622-709 (1998).
  * 33 Trefethen, L.N. Group velocity in finite difference schemes. SIAM Review, 24(2), 113–136 (1982).
  * 34 LeVeque, R.J. Finite Difference Methods for Ordinary and Partial Differential Equations. Steady State and Time-Dependent Problems. SIAM, Philadelphia, USA (2007).
  * 35 Strikwerda, J.C. Finite Difference Schemes and Partial Differential Equations. $2^{nd}$ Ed., SIAM, Philadelphia, USA (2004).
  * 36 Tam, C.K.W. and Webb, J.C. Dispersion-relation-preserving finite difference schemes for computational acoustics. J. Comput. Phys., 107, 262–281 (1993).
  * 37 Sengupta, T.K. Sengupta, A. and Saurabh, K. Global spectral analysis of multi-level time integration schemes: Numerical properties for error analysis. Appl. Math. Computation, 304, 41–57 (2017).
  * 38 Sengupta, T.K., Sagaut, P., Sengupta, A. and Saurabh, K. Global spectral analysis of three-time level integration schemes: Focusing phenomenon, Comput. Fluids, 157, 182–195 (2017).
  * 39 Wang, G., Duchaine, F., Papadogiannis, D., Moreau, S. and Gicquel, L. An overset grid method for large eddy simulation of turbomachinery stages. J. Comput. Phys., 274, 333–355 (2014)
  * 40 Brunet, V., Croner, E., Minot, A., de Laborderie, J., Lippinois, E., Richard, S., Boussuge, J.-F., Dombard, J., Duchaine, F., Gicquel, L. et al. Comparison of various CFD codes for LES simulations of turbomachinery: from inviscid vortex convection to multi-stage compressor. ASME Turbo Expo, Jun 2018, Oslo, Norway. 10.1115/GT2018-75523. hal-02003192 (2018)
  * 41 Schoenfeld, T. and Rudgyard, M. Steady and unsteady flow simulations using the hybrid flow solver AVBP. AIAA J., 37(11), 1378–1385 (1999)
  * 42 Rochette, B., Collin-Bastiani, F., Gicquel, L., Vermorle, O., Veynante, D. and Poinsot, T. Influence of chemical schemes, numerical method and dynamic turbulent combustion modelling on LES of premixed turbulent flames. Combustion and Flame, 191, 417–430 (2018)
  * 43 Suman, V. K., Sengupta, T. K., Durga Prasad, C. J., Surya Mohan, K. and Sanwalia, D. Spectral analysis of finite difference schemes for convection diffusion equation. Comput. Fluids, 150, 95–114 (2017).
  * 44 Suman, V. K., Sengupta, T. K. and Mathur, J. S. Effects of numerical anti-diffusion in closed unsteady flows governed by two-dimensional Navier-Stokes equation. Comput. Fluids, 201, 104479 (2020)
  * 45 Sengupta, S., Sengupta, T. K., Puttam, J. K. and Suman, V. K. Global Spectral Analysis for Convection-Diffusion-Reaction Equation in One and Two-dimensions: Effects of Numerical Anti-Diffusion and Dispersion. J. Comput. Phys. 109310 (2020)
  * 46 Sharma, N., Sengupta, A., Rajpoot, M.K., Samuel, R.J. and Sengupta, T.K. Hybrid sixth order spatial discretization scheme for non-uniform Cartesian grids. Comput. Fluids, 157, 208-231 (2017)
  * 47 Sengupta, T. K. and Suman, V. K., Focusing phenomenon in numerical solution of two-dimensional Navier-Stokes equation. In the Proc. of High-Performance Computing of Big Data for Turbulence and Combustion, (Eds. Pirozzoli & Sengupta). Springer, Cham, 1–29 (2019)
  * 48 Sengupta, T. K., Lakshmanan, V. and Vijay, V. V. S. N.  A new combined stable and dispersion relation preserving compact scheme for non-periodic problems.  J. Comput. Phys., 228(8), 3048–3071 (2009)
  * 49 Sengupta, T. K., Vijay, V. V. S. N. and Bhaumik, S.  Further improvement and analysis of CCD scheme: Dissipation discretization and de-aliasing properties.  J. Comput. Phys., 228(17), 6150–6168 (2009)
  * 50 Van der Vorst, H. A. Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of non-symmetric linear systems. SIAM J. Sci. Stat. Comput., 12, 631–644 (1992)
  * 51 Sengupta, T.K., Suman, V.K., Sundaram, P. and Sengupta, A. Analysis of Pseudo-spectral Methods Used for Numerical Simulations of Turbulence. WSEAS Trans. Comp. Res., 10, 9–24 (2022)
  * 52 Sengupta, S., Sreejith, N.A., Mohanamuraly, P., Staffelbach, G., and Gicquel, L. Global spectral analysis of the Lax–Wendroff-central difference scheme applied to Convection–Diffusion equation Comp. Fluids, 242, 105508 (2022)
  * 53 Lou, S., Chao, Y., Ma, L-B and Jiang, Z-H The Flux Reconstruction Method with Lax–Wendroff Type Temporal Discretization for Hyperbolic Conservation Laws J. Sci. Comp., 82(42) (2020)
  * 54 Burger, R., Sudarshan Kumar, K., and Zorio, D. Approximate Lax–Wendroff discontinuous Galerkin methods for hyperbolic conservation laws Comp. Math. Appl., 74(6), 1288-1310 (2017)
  * 55 Sengupta, .T K., and Bhole, A. Error dynamics of diffusion equation: Effects of numerical diffusion and dispersive diffusion J. Comp. Phys., 266, 240-251 (2014)