# Preconditioned Central Moment Lattice Boltzmann Method on a Rectangular
Lattice Grid for Accelerated Computations of Inhomogeneous Flows

Eman Yahia<EMAIL_ADDRESS>Kannan N. Premnath
<EMAIL_ADDRESS>Department of Mechanical Engineering
College of Engineering, Design and Computing
University of Colorado Denver
1200 Larimer street, Denver, Colorado 80217 , U.S.A.

###### Abstract

Convergence acceleration of flow simulations to their steady states at lower
Mach numbers can be achieved via preconditioning the lattice Boltzmann (LB)
schemes that alleviate the associated numerical stiffness, which have so far
been constructed on square lattices. We present a new central moment LB method
on rectangular lattice grids for efficient computations of inhomogeneous and
anisotropic flows by solving the preconditioned Navier-Stokes (PNS) equations.
Moment equilibria corrections are derived via a Chapman-Enskog analysis for
eliminating the truncation errors due to grid-anisotropy arising from the use
of the rectangular lattice and the non-Galilean invariant cubic velocity
errors resulting from an aliasing effect on the standard D2Q9 lattice for
consistently recovering the PNS equations. Such corrections depend on the
diagonal components of the velocity gradients, which are locally obtained from
the second order non-equilibrium moments and parameterized by an associated
grid aspect ratio $r$ and a preconditioning parameter $\gamma$, and the speed
of sound in the collision model is naturally adapted to $r$ via a physically
consistent strategy. We develop our approach by using a robust non-orthogonal
moment basis and the central moment equilibria are based on a matching
principle, leading to simpler expressions for the corrections for using the
rectangular grids and for representing the viscosities as functions of the
relaxation parameters, $r$ and $\gamma$, and its implementation is modular
allowing a ready extension of the existing LB schemes based on the square
lattice. Numerical simulations of inhomogeneous and anisotropic shear-driven
bounded flows using the preconditioned rectangular central moment LB method
demonstrate the accuracy and significant reductions in the numbers of steps to
reach the steady states for various sets of characteristic parameters.

###### keywords:

Lattice Boltzmann method, Rectangular lattice, Central moments,
Preconditioning , Inhomogeneous and anisotropic flows , Convergence
acceleration

††journal: Journal Name

##  1 Introduction

With the rapid development in the computational fluid dynamics (CFD)
techniques and their engineering and scientific applications, the need for
accurate and efficient simulations of fluid flows has become vitally
important. The lattice Boltzmann (LB) method has proved especially
advantageous in simulating a variety of fluid flows, including complex flows
such as multiphase and multicomponent systems, porous media flows, and
turbulence [1, 2, 3]. Its firm basis resulting from certain a discretization
of the Boltzmann equation and its ability to accommodate considerations beyond
hydrodynamics, such as the higher order kinetic moments, contributed to many
refinements of this approach. The LB method is naturally parallelizable,
flexible in adopting models from kinetic theory, and it facilitates boundary
condition implementations on Cartesian grids with relative ease, which have
paved the way for its growing number of applications. Briefly, the method
involves tracking the spatial and temporal evolution of the distribution of
the particle populations due to collisions and advection along the
characteristic discrete velocity directions referred to as a lattice. The
collision step is often represented by a model involving the relaxation to
certain equilibria, either directly involving the distribution functions [4]
or their raw moment [5], central moment [6] or cumulant [7] representations,
performed with using a single relaxation time [4] or multiple relaxation times
[5], or by a model that is compliant with certain notions of entropy [8]. The
asymptotic continuum limit of such collide-and-stream steps on a lattice
satisfying the necessary symmetry and isotropy considerations correspond to
the dynamics of the fluid flow represented by the Navier-Stokes (NS)
equations.

The LB method typically uses uniform Cartesian grids, resulting, for example,
from the choice of a square lattice in two-dimensions (2D). Real world
problems are often dominated by inhomogeneous and anisotropic flows, including
wall-bounded shear flows. For example, in turbulent boundary layers or flow
through channels or ducts, the eddy sizes are markedly different in different
coordinate directions, and which progressively increase in the direction
normal to the wall. Similar situations arise in simulating flows within
enclosures with high geometric skewness characterized by large disparities in
the length scales or the aspect ratios. Thus, the use of square/cubic grids
are to solve such problems of practical interest are associated with high
computational costs both in terms of time and memory, which can be orders of
magnitude more than efficient algorithms based on nonuniform, stretched grids.
Hence, it becomes highly important to develop more efficient approaches that
use grids which naturally accommodate the spatial variations in the flow
features. As a result, much attention has been paid to extend the LB schemes
involving non-uniform grids. While unstructured grids or schemes involving
interpolations for considering stretched grids can be used in this regard (see
e.g., the monograph [9] for a survey of such methods), it can result in
complicated implementations or introduce additional numerical dissipation
[10]. On the other hand, one of the hallmarks of the LB methods is that the
perfect lock-step advection or streaming used in their standard formulations
incurs minimal overall dissipation while maintaining the simplicity of its
implementation. For more efficiently simulating anisotropic flows while
preserving this important feature, it becomes natural to utilize rectangular
lattices rather than square lattices in 2D.

Thus, significant focus has been shown towards developing LB methods using
rectangular lattice grids during the last decade following the initial
investigation in this direction by Koelman [11]. For example, Hegele _et al_
[12], Peng _et al_ [13], Wang _et al_ [14] presented rectangular LB algorithms
based on a single relaxation time (SRT) model via an extended lattice set,
corrections to equilibrium distribution functions, and counteracting source
terms, respectively, to recover the NS equations. Moreover, rectangular LB
schemes based on raw moments using multiple relaxation times (MRT) were
developed by Bouzidi _et al_ [15] and Zhou _et al_ [16], which were analyzed
and an improved LB scheme on a rectangular grid with the necessary correction
terms that is consistent with the NS equations was constructed by Peng _et al_
[17]. However, many of these methods involved cumbersome implementations,
complicated expressions for the corrections, and numerical stability issues
when the grid aspect ratio of the rectangular lattice (defined later) is
significantly far off from unity (i.e., characterizing strong grid stretching
in one of the directions relative to the other) or for simulating flows with
relatively low viscosities or high Reynolds numbers. On the other hand,
recognizing that the use of central moments, which naturally preserves the
Galilean invariance of those moments independently supported by the lattice,
can significantly improve the stability and accuracy when compared to the use
of raw moments [6, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 31,
32, 33], we recently constructed a rectangular central moment LB method (RC-
LBM) [34], which was then further extended to three-dimensions with an
improved implementation strategy [35]. While the original central moment LB
scheme was constructed using an orthogonal moment basis [6], Geier _et al._
[7] in 2015 provided a detailed discussion on the role of the moment basis in
their development of a cumulant LB method and also constructed a variety of
collision models, including those based on raw moments, central moments and
cumulants using non-orthogonal moment basis and presented them in the various
appendices of [7]. Moreover, the numerical stability advantages of using such
non-orthogonal moment basis relative to the orthogonal moment basis were
demonstrated via a linear stability analysis in [36]. Besides, earlier studies
on cascaded LB schemes performed mathematical analysis and demonstrated
consistency to the Navier-Stokes equations using such simpler basis [18, 19].
The use of non-orthogonal central moments in algorithmic implementations in LB
schemes was later adopted in Refs. [22, 23]. Hence, in contrast to the prior
rectangular LB schemes, the RC-LBM used a natural non-orthogonal moment basis,
a matching principle to construct the equilibria involving higher order
velocity terms resulting in a simpler and significantly more robust
implementation. Also, in Ref. [35], we also explicitly demonstrated the
computational advantages of using a rectangular lattice in lieu of a square
lattice in solving inhomogeneous and anisotropic flows. Moreover, the central
moment LB method on a cuboid lattice presented in Ref. [35] is modular in
construction thereby allowing ready extension of the existing algorithms on a
cubic lattice to cuboid lattices, and provides a unified formulation with
corrections that are applicable for a wide variety of all the standard
collision models.

Since the LB methods, which are time-marching and weakly compressible flow
solvers, represent the fluid motion in the incompressible limit
asymptotically, the smaller the Mach number is used simulations, the better is
their accuracy, which also enhances their numerical stability. However, a
general issue in CFD, including those for the LB methods, for computing flows
via reducing the Mach number to relatively small values is the associated
increase in the stiffness that causes slower convergence to the steady state.
This is due to the large condition number of evolution equations resulting
from wide contrasts in the flow speeds and the acoustic speeds in such cases.
One approach to reduce the number of steps taken for convergence is to
precondition the system of flow equations, wherein such disparities between
the characteristic speeds are reduced at the cost of sacrificing the temporal
accuracy. As shown by Turkel (see e.g., [37, 38]), this has been accomplished
in the context of classical CFD methods by solving the so-called
preconditioned NS equations which involve an adjustable preconditioning
parameter. Guo _et al_ [39] introduced the first preconditioned LB scheme
using an SRT model, which was then extended to a MRT model involving raw
moments and forcing terms by Premnath _et al_ [40]. Izquierdo and Fueyo [41]
demonstrated optimal preconditioning of a MRT-LB scheme, while Meng _et al_
[42] introduced a preconditioned MRT-LB algorithm for simulations of steady
two-phase flows in porous media. More recently, Hajabdollahi and Premnath [43]
developed a cascaded central moment LB scheme for simulation of preconditioned
NS equations, which was then further improved by eliminating the non-Galilean
invariant cubic velocity errors that are dependent on the preconditioning
parameter and demonstrating significant reductions in the number of steps for
convergence of a variety of flows in Ref. [26]. Moreover, we note that,
recently, a preconditioned SRT-LB approach based on a finite-volume
discretization on unstructured grids was developed and studied by Walsh and
Boyle [44]. However, generally, prior investigations constructed
preconditioned LB algorithms on square lattices, and rectangular LB schemes
for the solution of preconditioned NS equations while maintaining the collide-
and-stream steps with perfect lock-step advection have not yet been discussed
in the literature. Development and analysis of such preconditioned LB schemes
on rectangular lattices could enable convergence acceleration of inhomogeneous
and anisotropic flows thereby further improving the computational efficiency
achieved with the use of rectangular lattice grids, which is the focus of this
paper.

In this work, we aim to construct a new preconditioned rectangular central
moment lattice Boltzmann method (referred to as the PRC-LBM in what follows).
In this regard, we employ a simpler non-orthogonal moment basis and the
central moment equilibria are constructed by matching with those of the
continuous Maxwell distribution with appropriate modifications in order to
consistently recover the preconditioned NS equations. By performing a Chapman-
Enskog analysis, we will derive corrections to the equilibria that eliminate
the truncation errors due to grid anisotropy and the non-Galilean invariant
cubic velocity terms arising due to aliasing effects on the standard D2Q9
lattice appearing in the emergent equations under the asymptotic limit when
compared to the preconditioned NS equations. The resulting corrections will be
shown to depend on the preconditioning parameter, grid aspect ratio, and the
normal components of the velocity gradient tensor, where the latter will be
expressed in terms of second-order non-equilibrium moments which will allow
their computations locally without using finite difference approximations. It
may be noted that in our previous 2D rectangular LBM [34], the transformation
matrices for mappings between the distribution functions and raw moments,
which depend on the grid aspect ratio, are constructed to separate the trace
of the second order moments from its other components that allows independent
specification of the bulk viscosity and shear viscosity, respectively. By
contrast, in this paper, following our recent work in Ref. [35], the PRC-LBM
will segregate the bulk viscosity from the shear viscosity only within the
step involving relaxation of moments to preconditioned equilibria under
collision, and the pre- and post-collision mapping matrices involve a simpler
moment basis and account for the grid aspect ratios only via certain diagonal
scaling matrices. This results in a modular implementation of the PRC-LBM, and
its formulation involves simpler expressions for the necessary corrections and
the transport coefficients. Moreover, we note while all the prior rectangular
LB schemes (e.g., [17, 34]) indicated that the speed of sound should be
adjusted to accommodate for its variations with the grid aspect ratio when
compared to the speed sound for the square lattice, they did not provide any
rationale or explicit formulas to accomplish this other than providing some
tabulated data. In this work, we provide some physical arguments to
consistently obtain the speed of sound for the rectangular lattice, with its
explicit parametrization by the grid aspect ratio. Finally, we will perform
some numerical studies to demonstrate the accuracy of the PRC-LBM and
reductions in the number of steps for convergence to the steady state for
simulations of selected cases of inhomogeneous and anisotropic flows for
different choices of the preconditioning parameter and the grid aspect ratio.

This paper is organized as follows. The following section (Sec. 2) discusses a
consistent approach for the selection of the speed of sound for lattice
schemes based on rectangular lattice grids. Next, in Sec. 3, we present a
Chapman-Enskog analysis of the preconditioned non-orthogonal moment LB
formulation on a rectangular D2Q9 lattice and derive the correction terms
necessary to eliminate the truncation errors due to grid anisotropy and non-
Galilean invariant velocity terms arising from aliasing effects. Such
corrections are shown to be parameterized by the grid aspect ratio,
preconditioning parameter and the velocity gradients, where the latter are
obtained locally from non-equilibrium moments. The construction of the
preconditioned rectangular central moment LBM for an efficient implementation
is discussed in Sec. 4, with the attendant algorithmic details of the PRC-LBM
provided in A. Then, in Sec. 5, we present numerical results for some case
studies involving anisotropic and inhomogeneous shear flows, validating the
PRC-LBM for accuracy and demonstrating convergence acceleration via
preconditioning on rectangular lattice grids for various characteristic
parameters. Moreover, comparisons between the preconditioned rectangular
central moment LBM (PRC-LBM) another formulation involving the preconditioned
rectangular raw moment LBM are made in Sec. 6. The conclusions of this work
are highlighted in Sec. 7.

##  2 Selection of Speed of Sound on Rectangular Lattice Grid for Physical
Consistency

Before discussing the preconditioned rectangular central moment LB scheme and
its analysis, we will now present a general physical consideration on the
selection of the speed of sound for rectangular lattice grids and its relation
to the sound speed of the usual square lattice. For the D2Q9 _square_ lattice
shown in Fig. 1(a), with a lattice spacing $\Delta x$ and a time step $\Delta
t$ resulting in the particle speed $c=\Delta x/\Delta t$, based on
considerations of isotropy and Galilean invariance, it is well known that its
optimal value of the speed of sound $c_{s*}$ is given by

$c_{s*}=\frac{1}{\sqrt{3}}c=\frac{1}{\sqrt{3}}\frac{\Delta x}{\Delta
t}=\frac{1}{\sqrt{3}}.$ (1)

Thus, $c_{s*}=1/\sqrt{3}$ in the usual lattice units (i.e., when $\Delta
x=\Delta t=1.0$). For the two possible arrangements of the D2Q9 _rectangular_
lattice shown in Figs. 1(b) and 1(c), we can generally define a _grid aspect
ratio_ $r$ representing the ratio of the grid spacing in the $y$ direction
with respect to that in the $x$ direction, i.e., $r=\Delta y/\Delta x$.

Figure 1: Two-dimensional, nine velocities (D2Q9) lattice with different
possible arrangements based on the grid aspect ratio $r$.

Denoting the particle speeds in the $x$ and $y$ directions, respectively, as
$c_{x}=\Delta x/\Delta t=c$ and $c_{y}=\Delta y/\Delta t$, it then follows
that $c_{x}=c$ and $c_{y}=r\;c$. From these two particle speeds, one may
introduce two possible values of the speed of sound $c_{sx}$ and $c_{sy}$ in
the two coordinate directions: $c_{sx}=c_{s*}$ and $c_{sy}=r\;c_{s*}$.
However, it can be readily established from a Chapman-Enskog analysis that the
pressure field is related to the local density times the square of the speed
of sound, which reflects the equation of state of a weakly compressible
athermal fluid motion represented by the LB method. Since the pressure field
is an isotropic and a scalar quantity, the rectangular LB schemes need to use
only one among the possible values of the speed of sound to satisfy physical
consistency. Moreover, given that the speed of sound is a fraction of the
particle speed, and in order to be consistent with the Courant-Friedrichs-Lewy
condition, we prescribe that the effective speed of sound on the rectangular
lattice be chosen as the one that provides the _minimum_ value among the two
possibilities, i.e., $c_{s}=\mbox{min}(c_{sx},c_{sy})$, so that it picks the
more limiting case. In other words, our selection procedure for the speed of
sound on the rectangular lattice $c_{s}$ relative to that for the
corresponding square lattice $c_{s*}$ can be written as

$c_{s}=q\;c_{s*},\quad q=\mbox{min}(1,r).$ (2)

Thus, if $r<1$, $c_{s}=rc_{s*}$ (see Fig. 1(b)) and when $r>1$, $c_{s}=c_{s*}$
(see Fig. 1(c)). Moreover, when $r=1$, it naturally recovers the optimal value
of $1/\sqrt{3}$ used for the square lattice. Thus, Eq. (2) automatically
adapts the sound speed according to the grid aspect ratio, $r$, unlike
previous LB schemes based on the rectangular lattice (e.g., [17, 34]), where
no such expressions are provided. Also, noting that a Chapman-Enksog analysis
relates the kinematic viscosity of the fluid $\nu$, to a relaxation parameter
and the square of the speed of sound $c_{s}^{2}$ (see the next section), it
follows that $\nu$ is then parameterized by $q^{2}$, which facilitates
maintaining numerical stability self consistently as the grid aspect ratio $r$
is varied with the use of rectangular lattice grids. Finally, we mention here
that the specification of the Mach number Ma for the rectangular LB
simulations should be based on Eq. (2), i.e., for any characteristic flow
speed $U$, $\mbox{Ma}=U/c_{s}=U/(qc_{s*})=\mbox{Ma}_{*}/q$, where
$\mbox{Ma}_{*}$ is the Mach number used for the square lattice.

We will formulate our central moment LB scheme on a $D2Q9$ rectangular lattice
grid, where the particle velocity components $e_{x}$ and $e_{y}$ in the
directions $x$ and $y$ (see Fig. 1) following the definition of the grid
aspect ratio $r$ can be written as

$\displaystyle\ket{e_{x}}$ $\displaystyle=$
$\displaystyle(0,1,0,-1,0,1,-1,-1,1)^{\dagger},$ (3a)
$\displaystyle\ket{e_{y}}$ $\displaystyle=$
$\displaystyle(0,0,r,0,-r,r,r,-r,-r)^{\dagger},$ (3b)

where $\ket{\cdot}$ denotes a column vector based on the standard ‘ket’
notation and ${\dagger}$ refers to taking the transpose. We will also need the
following 9-dimensional vector in defining the moment basis in the next
section:

$\displaystyle\ket{1}=(1,1,1,1,1,1,1,1,1)^{\dagger}.$ (4)

##  3 Chapman-Enskog Analysis of Preconditioned LBE a Rectangular Lattice
Grid: Isotropy Corrections, Macroscopic Flow Equations, and Local Expressions
for Velocity Gradients

In the following, we will construct a preconditioned rectangular central
moment LBM on the D2Q9 rectangular lattice with grid aspect ratio-adapted
speed of sound to solve the following preconditioned NS equations:

$\displaystyle\partial_{t}\rho+\bm{\nabla}\cdot\left(\rho\bm{u}\right)=0,$
(5a)
$\displaystyle\partial_{t}\left(\rho\bm{u}\right)+\bm{\nabla}\cdot\left(\cfrac{\rho\bm{u}\bm{u}}{\gamma}\right)=-\cfrac{1}{\gamma}\bm{\nabla}{p}+\cfrac{1}{\gamma}\;\bm{\nabla}\cdot\bm{\mathsf{\tau}}+\cfrac{\bm{F}}{\gamma},$
(5b)

where $\gamma$ is the preconditioned parameter used to achieve convergence
acceleration to the steady state, $\bm{u}$ and $\rho$ are the fluid velocity
and density, respectively, $p=\gamma c_{s}^{2}\rho$ is the pressure field and
$\bm{\mathsf{\tau}}$ is the viscous stress tensor
$\bm{\mathsf{\tau}}=\rho\nu(\bm{\nabla}\bm{u}+(\bm{\nabla}\bm{u})^{\dagger})$,
and $\bm{F}$ is the body force.

###  3.1 Moment basis, and definitions of central moments and raw moments

In this regard, as in our previous work [34, 35], we employ a linearly
independent set of non-orthogonal basis vectors for moments, by noting that
they are chosen to especially allow for the separation of the isotropic parts
from the non-isotropic parts of the second order moments for independent
specification of the transport coefficients (i.e., the shear and bulk
viscosities). For the D2Q9 lattice, such basis vectors are defined using a
combination of the monomials of the type $\ket{e_{x}^{m}e_{y}^{n}}$, where $m$
and $n$ are integer exponents, as follows:

$\bm{\mathsf{T}}=\Big{[}\ket{T_{0}},\ket{T_{1}},\ket{T_{2}},\ldots,\ket{T_{8}}\Big{]}^{{\dagger}},$
(6)

where

$\displaystyle\ket{T_{0}}=\ket{1},\qquad\ket{T_{1}}=\ket{e_{x}},\qquad\ket{T_{2}}=\ket{e_{y}},\qquad\ket{T_{3}}=\ket{e_{x}^{2}+e_{y}^{2}},\qquad\ket{T_{4}}=\ket{e_{x}^{2}-e_{y}^{2}},$
$\displaystyle\ket{T_{5}}=\ket{{e_{x}}{e_{y}}},\qquad\ket{T_{6}}=\ket{e_{x}^{2}{e_{y}}},\qquad\ket{T_{7}}=\ket{{e_{x}}e_{y}^{2}},\qquad\ket{T_{8}}=\ket{e_{x}^{2}e_{y}^{2}}.$
(7)

Then, defining the sets of discrete distribution functions $\mathbf{f}$, their
equilibria $\mathbf{f}^{eq}$, and the source terms $\mathbf{S}$, which
represent the effect of the body force $\bm{F}=(F_{x},F_{y})$ on the fluid
motion, respectively, as

$\displaystyle\mathbf{f}=\left(f_{0},f_{1},f_{2},\ldots,f_{8}\right)^{{\dagger}},\quad\mathbf{{f}}^{eq}=\left({f}_{0}^{eq},{f}_{1}^{eq},{f}_{2}^{eq},\ldots,{f}_{8}^{eq}\right)^{{\dagger}},\quad$
$\displaystyle\mathbf{S}=\left({S}_{0},{S}_{1},{S}_{2},\ldots,{S}_{8}\right)^{{\dagger}},$
(8a)

we can then express their _raw moments_ of order ($m+n$), $k_{mn}^{\prime}$,
$k_{mn}^{eq\prime}$, and $\sigma_{mn}^{\prime}$, respectively,

$\displaystyle k_{mn}^{\prime}=\sum_{\alpha=0}^{8}f_{\alpha}\;e_{\alpha
x}^{m}e_{\alpha y}^{n},$ (9a) $\displaystyle
k_{mn}^{eq\prime}=\sum_{\alpha=0}^{8}f_{\alpha}^{eq}\;e_{\alpha
x}^{m}e_{\alpha y}^{n},$ (9b)
$\displaystyle\sigma_{mn}^{\prime}=\sum_{\alpha=0}^{8}S_{\alpha}\;e_{\alpha
x}^{m}e_{\alpha y}^{n}.$ (9c)

Similarly, we can write the corresponding _central moments_ $k_{mn}$,
$k_{mn}^{eq}$ and $\sigma_{mn}$, respectively, by subtracting the particle
velocities ($e_{\alpha x},e_{\alpha y}$) by the fluid velocity ($u_{x},u_{y}$)
as follows:

$\displaystyle k_{mn}=\sum_{\alpha=0}^{8}f_{\alpha}\;(e_{\alpha
x}-u_{x})^{m}(e_{\alpha y}-u_{y})^{n},$ (10a) $\displaystyle
k_{mn}^{eq}=\sum_{\alpha=0}^{8}f_{\alpha}^{eq}\;(e_{\alpha
x}-u_{x})^{m}(e_{\alpha y}-u_{y})^{n}$ (10b)
$\displaystyle\sigma_{mn}=\sum_{\alpha=0}^{8}S_{\alpha}\;(e_{\alpha
x}-u_{x})^{m}(e_{\alpha y}-u_{y})^{n}.$ (10c)

For convenience, we collect the various raw moments supported by the lattice
set in view of the moment basis given in Eq. (3.1) in the form of the
following 9-dimensional vectors as

$\displaystyle\mathbf{n}=\left(k_{00}^{\prime},k_{10}^{\prime},k_{01}^{\prime},k_{20}^{\prime}+k_{02}^{\prime},k_{20}^{\prime}-k_{02}^{\prime},k_{11}^{\prime},k_{21}^{\prime},k_{12}^{\prime},k_{22}^{\prime}\right)^{{\dagger}},$
(11a)
$\displaystyle\mathbf{n}^{eq}=\left(k_{00}^{eq\prime},k_{10}^{eq\prime},k_{01}^{eq\prime},k_{20}^{eq\prime}+k_{02}^{eq\prime},k_{20}^{eq\prime}-k_{02}^{eq\prime},k_{11}^{eq\prime},k_{21}^{eq\prime},k_{12}^{eq\prime},k_{22}^{eq\prime}\right)^{{\dagger}},$
(11b)
$\displaystyle\mathbf{\Psi}=\Big{(}\sigma_{00}^{\prime},\sigma_{10}^{\prime},\sigma_{01}^{\prime},\sigma_{20}^{\prime}+\sigma_{02}^{\prime},\sigma_{20}^{\prime}-\sigma_{02}^{\prime},\sigma_{11}^{\prime},\sigma_{21}^{\prime},\sigma_{12}^{\prime},\sigma_{22}^{\prime}\Big{)}^{{\dagger}},$
(11c)

Then, the mappings between the various raw moments and the distribution
functions can be compactly expressed via the matrix $\bm{\mathsf{T}}$ as

$\displaystyle\mathbf{n}=\bm{\mathsf{T}}\mathbf{f},\quad\mathbf{n}^{eq}=\bm{\mathsf{T}}\mathbf{f^{eq}},\quad\mathbf{\Psi}=\bm{\mathsf{T}}\mathbf{S}.$
(12)

Here, it should be mentioned that we use the combinations of the second order
moment basis $\ket{e_{x}^{2}+e_{y}^{2}}$ and $\ket{e_{x}^{2}-e_{y}^{2}}$ to
retain the flexibility of an independent specification of the bulk viscosity
and shear viscosity, and the matrix $\bm{\mathsf{T}}$ as formulated above then
facilitates in the demonstration of the consistency of our approach with the
preconditioned NS equations and in the derivation of the required attendant
corrections in the reminder of this section. However, in the actual
implementation of the algorithm in the next section (see Sec. 4), we introduce
the effects equivalent to the independent evolution of the moments related to
$\ket{e_{x}^{2}+e_{y}^{2}}$ and $\ket{e_{x}^{2}-e_{y}^{2}}$ only within the
sub-step involving the relaxations under collision and not for performing the
mappings between the distribution functions and moments.

###  3.2 Preconditioned Lattice Boltzmann Equation

Next, it is important to note that the use of the rectangular lattice would
result in an anisotropic form of the viscous stress tensor dependent on the
grid aspect ratio $r$. Such spurious effects along with the truncation errors
arising from the non-Galilean invariant aliasing effects on the D2Q9 lattice
dependent on the cubic velocity terms and the preconditioning parameter need
to be eliminated via certain counteracting corrections, which appear in the
evolution of the non-equilibrium part of the second order moments (see Ref
[34]). Since by construction, the non-equilibrium second order central moments
are identical to those of raw moments, it suffices to perform a consistency
analysis and derive the necessary correction terms based on the simpler
preconditioned rectangular raw moment MRT formulation of the lattice Boltzmann
equation (MRT-LBE) written in a compact matrix-vector form given by

$\mathbf{f}(\bm{x}+\mathbf{e}\Delta t,t+\Delta
t)-\mathbf{f}(\bm{x},t)=\bm{\mathsf{T^{-1}}}\Big{[}\bm{\mathsf{\Lambda}}\;\left(\;\mathbf{n}^{eq}-\mathbf{n}\;\right)+\left(\bm{\mathsf{I}}-\frac{\Lambda}{2}\right)\mathbf{\Psi}\Delta
t\Big{]},$ (13)

where $\bm{\mathsf{I}}$ is an identity matrix of dimension $9\times 9$ and
$\bm{\mathsf{\Lambda}}$ is a diagonal matrix holding the relaxation parameters
given by

$\bm{\mathsf{\Lambda}}=\mbox{diag}\;\big{(}0,0,0,\omega_{3},\omega_{4},\omega_{5},\omega_{6},\omega_{7},\omega_{8}\big{)}.$
(14)

The solution of this LBE (Eq. (13)) yields the distribution functions
$f_{\alpha}=f_{\alpha}(\bm{x},t+\Delta t)$, whose leading order moments then
provide the hydrodynamic fields as

$\rho=\sum_{\alpha=0}^{8}f_{\alpha},\qquad\rho\bm{u}=\sum_{\alpha=0}^{8}f_{\alpha}\bm{e}_{\alpha}+\frac{\bm{F}}{2\gamma}\Delta
t,\qquad p=\gamma c_{s}^{2}\rho.$ (15)

The key issue here is the specification of the moment equilibria components
appearing in $\mathbf{n}^{eq}$ in Eq. (13) so that the preconditioned NS
equations (Eq. (5)) can be recovered consistently on a rectangular lattice
grid.

###  3.3 Preconditioned Equilibria and Sources: Raw Moments and Central
Moments

In this regard, our starting point is the matching of the components of the
discrete raw moment equilibria supported by the D2Q9 lattice with those
following from the continuous Maxwell distribution, where the speed of sound
is based on Eq. (2). Then, we account for the modifications needed for
recovering the preconditioned NS equations, which was obtained in an earlier
analysis performed in Ref. [26] in the case of the _square_ lattice. Using
these as our initial formulation, we can then write the leading terms of the
raw moment equilibria as

$\displaystyle k_{00}^{eq\prime}=\rho,\quad\quad k_{10}^{eq\prime}=\rho
u_{x},\quad\quad k_{01}^{eq\prime}=\rho u_{y},$ $\displaystyle
k_{20}^{eq\prime}=q^{2}c_{s*}^{2}\rho+\frac{\rho
u_{x}^{2}}{\gamma},\quad\quad\qquad
k_{02}^{eq\prime}=q^{2}c_{s*}^{2}\rho+\frac{\rho
u_{y}^{2}}{\gamma},\quad\quad\quad k_{11}^{eq\prime}=\frac{\rho
u_{x}u_{y}}{\gamma},$ $\displaystyle
k_{21}^{eq\prime}=\rho\left(q^{2}c_{s*}^{2}+\frac{u_{x}^{2}}{\gamma^{2}}\right)u_{y},\quad\quad
k_{12}^{eq\prime}=\rho\left(q^{2}c_{s*}^{2}+\frac{u_{y}^{2}}{\gamma^{2}}\right)u_{x},\quad
k_{22}^{eq\prime}=\rho q^{4}c_{s*}^{4}+\rho
q^{2}c_{s*}^{2}\left(u_{x}^{2}+u_{y}^{2}\right)+\rho u_{x}^{2}u_{y}^{2}.$ (16)

The expressions in Eq. (3.3) need to be corrected further to consistently
recover the preconditioned NS equations on _rectangular_ lattice grids, which
will be accomplished later in this section. Moreover, the raw moment
equilibria of the source terms also need to be scaled appropriately by the
preconditioning parameter $\gamma$ as follows [40, 26]:

$\displaystyle\sigma_{00}^{\prime}=0,\quad\sigma_{10}^{\prime}=\frac{F_{x}}{\gamma},\quad\sigma_{01}^{\prime}=\frac{F_{y}}{\gamma},\qquad\sigma_{20}^{\prime}=\frac{2F_{x}u_{x}}{\gamma^{2}},\quad\sigma_{02}^{\prime}=\frac{2F_{y}u_{y}}{\gamma^{2}},\quad\sigma_{11}^{\prime}=\frac{\left(F_{x}u_{y}+F_{y}u_{x}\right)}{\gamma^{2}},$
(17)

and $\hat{\sigma}_{mn}^{\prime}=0,$ if $(m+n)\geq 3$.

Then, the countable set of preconditioned discrete central moment equilibria
on the D2Q9 lattice can be obtained from the corresponding the raw moment
equilibria given in Eq. (3.3) via the binomial transformations as

$\displaystyle k_{00}^{eq}=\rho,\qquad k_{10}^{eq}=0,\qquad k_{01}^{eq}=0,$
$\displaystyle
k_{20}^{eq}=q^{2}c_{s*}^{2}\rho+\left(\frac{1}{\gamma}-1\right)\rho
u_{x}^{2},\qquad
k_{02}^{eq}=q^{2}c_{s*}^{2}\rho+\left(\frac{1}{\gamma}-1\right)\rho
u_{y}^{2},\qquad k_{11}^{eq}=\left(\frac{1}{\gamma}-1\right)\rho u_{x}u_{y},,$
$\displaystyle
k_{21}^{eq}=\left(\frac{1}{\gamma^{2}}-\frac{3}{\gamma}+2\right)\rho
u_{x}^{2}u_{y},\qquad
k_{12}^{eq}=\left(\frac{1}{\gamma^{2}}-\frac{3}{\gamma}+2\right)\rho
u_{x}u_{y}^{2},\qquad k_{22}^{eq}=q^{4}c_{s*}^{4}\rho,$ (18)

Note that since the fourth order component of the equilibrium central moment
$k_{22}^{eq}$ does not appear at the leading order in Chapman-Enskog analysis
of the preconditioned NS equations, for simplicity, we set it as
$k_{22}^{eq}=q^{4}c_{s*}^{4}\rho$ following our previous work [34]. Moreover,
similarly the central moment components of the source terms follow from the
corresponding raw moments Eq. (17) using binomial expansions as

$\displaystyle\sigma_{00}=0,\quad\sigma_{10}=\frac{F_{x}}{\gamma},\quad\sigma_{01}=\frac{F_{y}}{\gamma},$
$\displaystyle\sigma_{20}=\left(\frac{1}{\gamma^{2}}-\frac{1}{\gamma}\right)2F_{x}u_{x},\quad\sigma_{02}=\left(\frac{1}{\gamma^{2}}-\frac{1}{\gamma}\right)2F_{y}u_{y},\quad\sigma_{11}=\left(\frac{1}{\gamma^{2}}-\frac{1}{\gamma}\right)\left(F_{x}u_{y}+F_{y}u_{x}\right),$
(19)

and $\sigma_{mn}=0,$ if $(m+n)\geq 3$. An alternative approach to implementing
body forces in central moment LB schemes has been proposed by Fei and Luo [24,
25]. It has been used for various applications, including thermal flows and
multiphase flows, and involves including the effect of the body force on the
higher order central moments. By contrast, the forcing scheme given above and
others such as in [7, 45] involve the effect of body forces up to the second
order moments, and are constructed to recover the hydrodynamics (Navier-Stokes
equations) as prescribed by the Chapman-Enskog analysis.

###  3.4 Chapman-Enskog Analysis: Identification of Truncation Errors due to
Grid Anisotropy, Preconditioning and Non-Galilean Invariance from Aliasing
Effects on the D2Q9 Rectangular Lattice

Next, we will perform a Chapman-Enskog (C-E) analysis [46] in order to
determine the truncation errors arising from grid anisotropy with the use of
the rectangular lattice and the non-Galilean invariant (GI) cubic velocity
terms due to the aliasing effects manifesting as a result of the discreteness
of the D2Q9 lattice. This would be carried out following the approach taken in
our previous works [19, 26, 34]. First, expanding the moments about their
equilibria and the time derivative by means of a multiple time expansion, we
write

$\mathbf{n}=\sum_{j=0}^{\infty}\epsilon^{j}\mathbf{n}^{(j)},\quad\partial_{t}=\sum_{j=0}^{\infty}{\epsilon}^{j}{\partial_{t_{j}}},$
(20)

where $\epsilon=\Delta t$ represents the perturbation parameter serving in
what follows to delineating the terms of different orders. Substituting the
above equation in Eq. (13) and rewriting its left side via a Taylor series
expansion and converting the resulting expression in terms of moments using
$\mathbf{f}=\bm{\mathsf{T}}^{-1}\mathbf{n}$, we obtain the evolution equations
of the moments of different orders of $\epsilon$, i.e., $O(\epsilon^{k})$,
where $k=0,1$, and $2$ as

$\displaystyle O(\epsilon^{0}):\mathbf{n}^{(0)}=\mathbf{n}^{eq},$ (21a)
$\displaystyle
O(\epsilon^{1}):\left(\partial_{t_{0}}+\bm{E}_{i}\partial_{i}\right)\mathbf{n}^{(0)}=-\bm{\mathsf{\Lambda}}\;\mathbf{n}^{(1)}+\mathbf{\Psi},$
(21b) $\displaystyle
O(\epsilon^{2}):\partial_{t_{1}}\;\mathbf{n}^{(0)}+\left(\partial_{t_{0}}+\bm{E}_{i}\partial_{i}\right)\;\left[\bm{\mathsf{I}}-\frac{\bm{\mathsf{\Lambda}}}{2}\right]\mathbf{n}^{(1)}=-\bm{\mathsf{\Lambda}}\;\mathbf{n}^{(2)},$
(21c)

where
$\bm{\mathsf{E}}_{i}=\bm{\mathsf{T}}\;(\mathbf{e}_{i}\;\bm{\mathsf{I}})\bm{\mathsf{T}}^{-1}$
and $\mathbf{e}_{i}=\ket{e_{i}}$, $i\in(x,y)$. Then, substituting the raw
moments and the source terms shown in Eqs. (3.3) and (17), respectively, into
Eq. (21b), the relevant moment system $O(\epsilon)$ which are relevant in
recovering the preconditioned hydrodynamics can be written as

$\displaystyle\partial_{t_{0}}\rho+\partial_{x}\rho u_{x}+\partial_{y}\rho
u_{y}=0,$ (22a) $\displaystyle\partial_{t_{0}}\rho u_{x}+\partial_{x}(\rho
q^{2}c_{s*}^{2}+\rho u_{x}^{2}/\gamma)+\partial_{y}(\rho
u_{x}u_{y}/\gamma)=F_{x}/\gamma,$ (22b) $\displaystyle\partial_{t_{0}}\rho
u_{y}+\partial_{x}(\rho u_{x}u_{y}/\gamma)+\partial_{y}(\rho
q^{2}c_{s*}^{2}+\rho u_{y}^{2}/\gamma)=F_{y}/\gamma,$ (22c)
$\displaystyle\partial_{t_{0}}\left[2\rho
q^{2}c_{s*}^{2}+\rho(u_{x}^{2}+u_{y}^{2})/\gamma\right]+\partial_{x}\left[(1+q^{2}c_{s*}^{2})\rho
u_{x}+\rho
u_{x}u_{y}^{2}/\gamma^{2}\right]+\partial_{y}\left[(r^{2}+q^{2}c_{s*}^{2})\rho
u_{y}+\rho u_{x}^{2}u_{y}/\gamma^{2}\right]=$
$\displaystyle-\omega_{3}\;n_{3}^{(1)}+2\left(F_{x}u_{x}+F_{y}u_{y}\right)/\gamma^{2},$
(22d)
$\displaystyle\partial_{t_{0}}\left[\rho(u_{x}^{2}-u_{y}^{2})/\gamma\right]+\partial_{x}\left[(1-q^{2}c_{s*}^{2})\rho
u_{x}-\rho
u_{x}u_{y}^{2}/\gamma^{2}\right]+\partial_{y}\left[(-r^{2}+q^{2}c_{s*}^{2})\rho
u_{y}+\rho u_{x}^{2}u_{y}/\gamma^{2}\right]=$
$\displaystyle-\omega_{4}\;n_{4}^{(1)}+2\left(F_{x}u_{x}-F_{y}u_{y}\right)/\gamma^{2},$
(22e) $\displaystyle\partial_{t_{0}}(\rho
u_{x}u_{y}/\gamma)+\partial_{x}\left[q^{2}c_{s*}^{2}\rho u_{y}+\rho
u_{x}^{2}u_{y}/\gamma^{2}\right]+\partial_{y}\left[q^{2}c_{s*}^{2}\rho
u_{x}+\rho u_{x}u_{y}^{2}/\gamma^{2}\right]=$
$\displaystyle-\omega_{5}\;n_{5}^{(1)}+\left(F_{x}u_{y}+F_{y}u_{x}\right)/\gamma^{2}.$
(22f)

Similarly, the $O(\epsilon^{2})$ evolution equations for the conserved moments
at the slower time scale $t_{1}$ reads from Eq. (21c) as

$\displaystyle\partial_{t_{1}}\rho=0,$ (23a)
$\displaystyle\partial_{t_{1}}\left(\rho
u_{x}\right)+\partial_{x}\left[\dfrac{1}{2}\left(1-\dfrac{\omega_{3}}{2}\right)n_{3}^{(1)}+\dfrac{1}{2}\left(1-\dfrac{\omega_{4}}{2}\right)n_{4}^{(1)}\right]+\partial_{y}\left[\left(1-\dfrac{\omega_{5}}{2}\right)n_{5}^{(1)}\right]=0,$
(23b) $\displaystyle\partial_{t_{1}}\left(\rho
u_{y}\right)+\partial_{x}\left[\left(1-\dfrac{\omega_{5}}{2}\right)n_{5}^{(1)}\right]+\partial_{y}\left[\dfrac{1}{2}\left(1-\dfrac{\omega_{3}}{2}\right)n_{3}^{(1)}-\dfrac{1}{2}\left(1-\dfrac{\omega_{4}}{2}\right)n_{4}^{(1)}\right]=0.$
(23c)

The above Eqs. (23a)-(23c) depend on the components of the non-equilibrium
moments $n_{3}^{(1)}$, $n_{4}^{(1)}$ and $n_{5}^{(1)}$, which can be obtained
from Eqs. (22)-(22). Hence, rewriting Eqs. (22)-(22) to express the non-
equilibrium moments as

$\displaystyle
n_{3}^{(1)}=\frac{1}{\omega_{3}}\Big{\\{}-\partial_{t_{0}}\left[2q^{2}c_{s*}^{2}\rho+\rho
u_{x}^{2}/\gamma+\rho
u_{y}^{2}/\gamma\right]-\partial_{x}\left[(1+q^{2}c_{s*}^{2})\rho
u_{x}\right]-\partial_{x}(\rho
u_{x}u_{y}^{2}/\gamma^{2})-\partial_{y}\left[(r^{2}+q^{2}c_{s*}^{2})\rho
u_{y}\right]$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\quad-\partial_{y}(\rho
u_{x}^{2}u_{y}/\gamma^{2})+2(F_{x}u_{x}+F_{y}u_{y})/\gamma^{2}\Big{\\}},$
(24a) $\displaystyle
n_{4}^{(1)}=\frac{1}{\omega_{4}}\Big{\\{}-\partial_{t_{0}}[(\rho
u_{x}^{2}-\rho u_{y}^{2})/\gamma]-\partial_{x}\left[(1-q^{2}c_{s*}^{2})\rho
u_{x}\right]+\partial_{x}(\rho
u_{x}u_{y}^{2}/\gamma^{2})+\partial_{y}\left[(r^{2}-q^{2}c_{s*}^{2})\rho
u_{y}\right]$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\quad-\partial_{y}(\rho
u_{x}^{2}u_{y}/\gamma^{2})+2(F_{x}u_{x}-F_{y}u_{y})/\gamma^{2}\Big{\\}},$
(24b) $\displaystyle
n_{5}^{(1)}=\frac{1}{\omega_{5}}\Big{\\{}-\partial_{t_{0}}(\rho
u_{x}u_{y}/\gamma)-\partial_{x}(q^{2}c_{s*}^{2}\rho u_{y})-\partial_{x}(\rho
u_{x}^{2}u_{y}/\gamma^{2})-\partial_{y}(q^{2}c_{s*}^{2}\rho
u_{x})-\partial_{y}(\rho
u_{x}u_{y}^{2}/\gamma^{2})+(F_{x}u_{y}+F_{x}u_{y})/\gamma^{2}\Big{\\}}.$ (24c)

Clearly, the second order non-equilibrium moments $n_{3}^{(1)}$, $n_{4}^{(1)}$
and $n_{5}^{(1)}$ involve terms related to the non-GI cubic velocity errors,
whose prefactors are dependent on the preconditioning parameter $\gamma$ and
the grid-anisotropy error terms dependent on the grid aspect ratio $r$, in
addition to those are related to the physically consistent terms that
contribute towards the viscous stress tensor. Denoting the truncation errors
related to the grid anisotropy as $E_{js}$ and the non-GI cubic velocity terms
as $E_{jg}$ for $j=3,4$ and $5$, and after replacing the time derivatives
appearing in Eqs. (24a)-(24c) in terms of the spatial derivatives of the
conserved moments via Eqs. (22a)-(22c), we can then simplify the resulting
equations by retaining terms up to $O(u_{i}^{3})$ (see Refs. [26] and [34] for
details). Then, the final expressions for the second order non-equilibrium
moment components on the rectangular lattice resulting for our LB formulation
can be written as follows:

$\displaystyle n_{3}^{(1)}$ $\displaystyle=$
$\displaystyle-\frac{2q^{2}c_{s*}^{2}}{\omega_{3}}\rho\bm{\nabla}\cdot\bm{u}+E_{3g}+E_{3s},$
(25a) $\displaystyle n_{4}^{(1)}$ $\displaystyle=$
$\displaystyle-\frac{2q^{2}c_{s*}^{2}}{\omega_{4}}\rho\left(\partial_{x}u_{x}-\partial_{y}u_{y}\right)+E_{4g}+E_{4s},$
(25b) $\displaystyle n_{5}^{(1)}$ $\displaystyle=$
$\displaystyle-\frac{q^{2}c_{s*}^{2}\rho}{\omega_{5}}\left(\partial_{x}u_{y}+\partial_{y}u_{x}\right)+E_{5g},$
(25c)

where the expressions for the truncation errors due to grid anisotropy
$E_{3s}$ and $E_{4s}$, and the non-GI cubic velocity aliasing errors $E_{3g}$,
$E_{4g}$,a nd $E_{5g}$ read as

$\displaystyle E_{3s}$ $\displaystyle=$
$\displaystyle\cfrac{1}{\omega_{3}}(3q^{2}c_{s*}^{2}-1)\rho\partial_{x}u_{x}+\cfrac{1}{\omega_{3}}(3q^{2}c_{s*}^{2}-r^{2})\rho\partial_{y}u_{y},$
(26a) $\displaystyle E_{4s}$ $\displaystyle=$
$\displaystyle\dfrac{1}{\omega_{4}}(3q^{2}c_{s*}^{2}-1)\rho\partial_{x}u_{x}-\dfrac{1}{\omega_{4}}(3q^{2}c_{s*}^{2}-r^{2})\rho\partial_{y}u_{y},$
(26b)

and

$\displaystyle E_{3g}$ $\displaystyle=$
$\displaystyle\cfrac{1}{\omega_{3}}\left[(2/\gamma+1)q^{2}c_{s*}^{2}-1\right]u_{x}\partial_{x}\rho+\cfrac{1}{\omega_{3}}\left[(2/\gamma+1)q^{2}c_{s*}^{2}-r^{2}\right]u_{y}\partial_{y}\rho+\cfrac{M_{3}}{\omega_{3}}\;\partial_{x}u_{x}+\cfrac{N_{3}}{\omega_{3}}\;\partial_{y}u_{y},$
(27a) $\displaystyle E_{4g}$ $\displaystyle=$
$\displaystyle\dfrac{1}{\omega_{4}}\left[(2/\gamma+1)q^{2}c_{s*}^{2}-1\right]u_{x}\partial_{x}\rho-\dfrac{1}{\omega_{4}}\left[(2/\gamma+1)q^{2}c_{s*}^{2}-r^{2}\right]u_{y}\partial_{y}\rho+\dfrac{M_{4}}{\omega_{4}}\partial_{x}u_{x}+\dfrac{N_{4}}{\omega_{4}}\partial_{y}u_{y},$
(27b) $\displaystyle E_{5g}$ $\displaystyle=$
$\displaystyle\frac{1}{\omega_{5}}(1/\gamma-1)q^{2}c_{s*}^{2}\left(u_{x}\partial_{y}\rho+u_{y}\partial_{x}\rho\right)+\frac{1}{\omega_{5}}(1/\gamma^{2}-1/\gamma)\rho
u_{x}u_{y}\left(\partial_{x}u_{x}+\partial_{y}u_{y}\right),$ (27c)

where the prefactors $M_{3}$, $N_{3}$, $M_{4}$ and $N_{4}$ appearing in Eqs.
(27a) and (27b) can be expressed as

$\displaystyle M_{3}$ $\displaystyle=$
$\displaystyle\rho\left[(4/\gamma^{2}-1/\gamma)u_{x}^{2}+(1/\gamma^{2}-1/\gamma)u_{y}^{2}\right],$
(28a) $\displaystyle N_{3}$ $\displaystyle=$
$\displaystyle\rho\left[(4/\gamma^{2}-1/\gamma)u_{y}^{2}+(1/\gamma^{2}-1/\gamma)u_{x}^{2}\right],$
(28b) $\displaystyle M_{4}$ $\displaystyle=$
$\displaystyle\rho\left[(4/\gamma^{2}-1/\gamma)u_{x}^{2}-(1/\gamma^{2}-1/\gamma)u_{y}^{2}\right],$
(28c) $\displaystyle N_{4}$ $\displaystyle=$
$\displaystyle\rho\left[-(4/\gamma^{2}-1/\gamma)u_{y}^{2}+(1/\gamma^{2}-1/\gamma)u_{x}^{2}\right].$
(28d)

Recognizing $q=\mbox{min}\\{r,1\\}$, it is evident that the various truncation
errors given above are dependent on the preconditioning parameter $\gamma$ and
the grid aspect ratio $r$, which need to be eliminated.

###  3.5 Corrections via Extended Moment Equilibria for Elimination of
Truncation Errors due to Grid Anisotropy, Preconditioning, and Aliasing
Effects

In this regard, we propose an extended moment equilibria
$\mathbf{n}^{\mathit{eq,eff}}$

$\displaystyle\mathbf{n}^{\mathit{eq,eff}}=\mathbf{n}^{eq}+\Delta
t\mathbf{n}^{eq(1)},$ (29)

where $\mathbf{n}^{eq(1)}$ are the corrections made to the base moment
equilibria $\mathbf{n}^{eq}$ introduced in Eqs. (11b) and (3.3). As shown in
Eqs. (25), (26), and (27), the truncation errors exist in the evolution of the
second order moments $n_{3}$, $n_{4}$ and $n_{5}$, and involve spatial
derivatives of the velocity field and the density field. Recognizing this
fact, for the purpose of consistently recovering the preconditioned NS
equations in a rectangular lattice grid, we thus write the following
expressions for the corrections to the moment equilibria, where it suffices to
introduce them to only the second order components:

$\displaystyle
n^{eq(1)}_{j}=\begin{cases}\theta_{3x}\partial_{x}u_{x}+\theta_{3y}\partial_{y}u_{y}+\lambda_{3x}\partial_{x}\rho+\lambda_{3y}\partial_{y}\rho&\quad
j=3\\\
\theta_{4x}\partial_{x}u_{x}-\theta_{4y}\partial_{y}u_{y}+\lambda_{4y}\partial_{x}\rho+\lambda_{4y}\partial_{y}\rho&\quad
j=4\\\
\theta_{5x}\partial_{x}u_{x}+\theta_{5y}\partial_{y}u_{y}+\lambda_{5y}\partial_{x}\rho+\lambda_{5y}\partial_{y}\rho&\quad
j=5\\\ 0&\quad\mbox{otherwise},\\\ \end{cases}$ (30)

Here, $\theta_{jx}$, $\theta_{jy}$, $\lambda_{jx}$, $\lambda_{jy}$, where
$j=3,4,$ and $5$, are the unknown coefficients, which will be determined by
carrying out a modified C-E expansion that includes the extended moment
equilibria given in Eq. (29). Thus, replacing the expansions appearing in Eq.
(20) with

$\displaystyle\mathbf{n}$ $\displaystyle=$
$\displaystyle\mathbf{n}^{\mathit{eq,eff}}+\epsilon\mathbf{n}^{(1)}+\epsilon^{2}\mathbf{n}^{(2)}+\ldots=\mathbf{n}^{(0)}+\underline{\epsilon\mathbf{n}^{eq(1)}}+\epsilon\mathbf{n}^{(1)}+\epsilon^{2}\mathbf{n}^{(2)}+\ldots$
$\displaystyle\partial_{t}$ $\displaystyle=$
$\displaystyle\partial_{t_{0}}+\epsilon\partial_{t_{1}}+\epsilon^{2}\partial_{t_{2}}+\ldots,$
(31)

and performing the same steps that follow Eq. (20) with using Eq. (3.5), then
the evolution of the moment systems at various orders of $\epsilon$ in Eq.
(21) modify to the following by accounting for the presence of the corrections
$\mathbf{n}^{eq(1)}$:

$\displaystyle O(\epsilon^{0}):\mathbf{n}^{(0)}=\mathbf{n}^{\mathit{eq}},$
(32a) $\displaystyle
O(\epsilon^{1}):\left({\partial_{t_{0}}}+\bm{\mathsf{E}}_{i}\partial_{i}\right)\mathbf{n}^{(0)}=-\bm{\mathsf{\Lambda}}\left[\mathbf{n}^{(1)}-\underline{\mathbf{n}^{eq(1)}}\right]+\mathbf{\Psi},$
(32b) $\displaystyle
O(\epsilon^{2}):{\partial_{t_{1}}}\mathbf{n}^{(0)}+\left({\partial_{t_{0}}}+\bm{\mathsf{E}}_{i}\partial_{i}\right)\left[\left(\bm{\mathsf{I}}-\frac{\bm{\mathsf{\Lambda}}}{2}\right)\mathbf{n}^{(1)}\right]+\left({\partial_{t_{0}}}+\bm{\mathsf{E}}_{i}\partial_{i}\right)\underline{\left[\frac{\bm{\mathsf{\Lambda}}}{2}\mathbf{n}^{eq(1)}\right]}=-\bm{\mathsf{\Lambda}}\mathbf{n}^{(2)}.$
(32c)

In view of the derivation given in the previous section and the changes
appearing in Eq. (32b) relative to Eq. (21b), the second order non-equilibrium
moments for the rectangular lattice with preconditioning in Eq. (25) modify to

$\displaystyle n_{3}^{(1)}$ $\displaystyle=$
$\displaystyle-\frac{2q^{2}c_{s*}^{2}}{\omega_{3}}\rho\bm{\nabla}\cdot\bm{u}+E_{3g}+E_{3s}+\underline{n_{3}^{eq(1)}},$
(33a) $\displaystyle n_{4}^{(1)}$ $\displaystyle=$
$\displaystyle-\frac{2q^{2}c_{s*}^{2}}{\omega_{4}}\rho\left(\partial_{x}u_{x}-\partial_{y}u_{y}\right)+E_{4g}+E_{4s}+\underline{n_{4}^{eq(1)}},$
(33b) $\displaystyle n_{5}^{(1)}$ $\displaystyle=$
$\displaystyle-\frac{q^{2}c_{s*}^{2}\rho}{\omega_{5}}\left(\partial_{x}u_{y}+\partial_{y}u_{x}\right)+E_{5g}+\underline{n_{5}^{eq(1)}},$
(33c)

where the error terms $E_{3g}$, $E_{3s}$, $E_{4g}$, $E_{4s}$, and $E_{5g}$ are
given in the previous section in Eqs. (26) and (27).

Now, in order to derive explicit formulas for the corrections $n_{3}^{eq(1)}$,
$n_{4}^{eq(1)}$ and $n_{5}^{eq(1)}$, we need certain constraint relationships
between them and the error terms. These follow from combining Eq. (32b) and
$\epsilon$ times Eq. (32c) and using the expressions for the non-equilibrium
moments in Eq. (33) along with
$\partial_{t}=\partial_{t_{0}}+\epsilon\partial_{t_{1}}$ to obtain the
effective evolution equations for the conserved moments, which would include
both the truncation error terms identified earlier and the unknown corrections
whose combined effects are set to zero so that the evolution equations
correspond to the preconditioned NS equations. See e.g., Refs. [26, 34, 35]
for details of these steps. Writing the truncation error terms compactly in
the form a vector $\mathbf{\Xi}$

$\mathbf{\Xi}=\left(\varphi_{0},\varphi_{1},\varphi_{2},\dots,\varphi_{8}\right)^{{\dagger}},$
(34)

where

$\displaystyle\varphi_{j}=\begin{cases}E_{3s}+E_{3g}&\quad j=3\\\
E_{4s}+E_{4g}&\quad j=4\\\ E_{5g}&\quad j=5\\\
0&\quad\mbox{otherwise},\end{cases}$ (35)

then the necessary constraint equation between the vector of moment
corrections $\mathbf{n}^{eq(1)}$ identified whose functional forms with
unknown coefficients are given in Eq. (30) and the above vector holding the
truncation errors $\mathbf{\Xi}$ (see Eqs. (34) and (35)) reads as

$\mathbf{n}^{eq(1)}+\left(\bm{\mathsf{I}}-\frac{\bm{\mathsf{\Lambda}}}{2}\right)\mathbf{\Xi}=0,$
(36)

which, in component form, becomes

$n_{j}^{eq(1)}+\left(1-\frac{\omega_{j}}{2}\right)(E_{js}+E_{jg})=0,\quad
j=3,4,5.$ (37)

Evaluating Eq. (37) and using Eqs. (30), (34) and (35) for $j=3,4$ and $5$,
respectively, we get

$\displaystyle\theta_{3x}\partial_{x}u_{x}+\theta_{3y}\partial_{y}u_{y}+\lambda_{3x}\partial_{x}\rho+\lambda_{3y}\partial_{y}\rho=-\left(1-\dfrac{\omega_{3}}{2}\right)E_{3g}-\left(1-\dfrac{\omega_{3}}{2}\right)E_{3s}$
$\displaystyle=-\left(\dfrac{1}{\omega_{3}}-\dfrac{1}{2}\right)\left[M_{3}+(q^{2}c_{s*}^{2}-1)\rho\right]\partial_{x}u_{x}-\left(\dfrac{1}{\omega_{3}}-\dfrac{1}{2}\right)\left[N_{3}+(q^{2}c_{s*}^{2}-r^{2})\rho\right]\partial_{y}u_{y}$
$\displaystyle-\left(\dfrac{1}{\omega_{3}}-\dfrac{1}{2}\right)\left[\left(\dfrac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-1\right]u_{x}\partial_{x}\rho-\left(\dfrac{1}{\omega_{3}}-\dfrac{1}{2}\right)\left[\left(\dfrac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-r^{2}\right]u_{y}\partial_{y}\rho.$
$\displaystyle\theta_{4x}\partial_{x}u_{x}-\theta_{4y}\partial_{y}u_{y}+\lambda_{4x}\partial_{x}\rho+\lambda_{4y}\partial_{y}\rho=-\left(1-\dfrac{\omega_{4}}{2}\right)E_{4g}-\left(1-\dfrac{\omega_{4}}{2}\right)E_{4s}$
$\displaystyle=-\left(\dfrac{1}{\omega_{4}}-\dfrac{1}{2}\right)\left[M_{4}+(q^{2}c_{s*}^{2}-1)\rho\right]\partial_{x}u_{x}-\left(\dfrac{1}{\omega_{4}}-\dfrac{1}{2}\right)\left[M_{4}+(q^{2}c_{s*}^{2}-r^{2})\rho\right]\partial_{y}u_{y}$
$\displaystyle-\left(\dfrac{1}{\omega_{4}}-\dfrac{1}{2}\right)\left[\left(\dfrac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-1\right]u_{x}\partial_{x}\rho+\left(\dfrac{1}{\omega_{4}}-\dfrac{1}{2}\right)\left[\left(\dfrac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-r^{2}\right]u_{y}\partial_{y}\rho,$
$\displaystyle\theta_{5x}\partial_{x}u_{x}+\theta_{5y}\partial_{y}u_{y}+\lambda_{5x}\partial_{x}\rho+\lambda_{5y}\partial_{y}\rho=-\left(1-\dfrac{\omega_{5}}{2}\right)E_{5g}$
$\displaystyle=-\left(\dfrac{1}{\omega_{5}}-\dfrac{1}{2}\right)\left(\dfrac{1}{\gamma^{2}}-\dfrac{1}{\gamma}\right)\rho
u_{x}u_{y}(\partial_{x}u_{x}+\partial_{y}u_{y})-\left(\dfrac{1}{\omega_{5}}-\dfrac{1}{2}\right)\left(\dfrac{1}{\gamma}-1\right)q^{2}c_{s*}^{2}\left(u_{x}\partial_{x}\rho+u_{y}\partial_{y}\rho\right),$

Comparing the terms involving the spatial gradients of the same type of
quantity in each side of the above three equations, we finally get the
coefficients for the correction terms in the second order moment equilibria as

$\displaystyle\theta_{3x}$ $\displaystyle=$
$\displaystyle-\Big{[}M_{3}+\left(3q^{2}c_{s*}^{2}-1\right)\rho\Big{]}\left(\frac{1}{\omega_{3}}-\frac{1}{2}\right),$
(38a) $\displaystyle\theta_{3y}$ $\displaystyle=$
$\displaystyle-\Big{[}N_{3}+\left(3q^{2}c_{s*}^{2}-r^{2}\right)\rho\Big{]}\left(\frac{1}{\omega_{3}}-\frac{1}{2}\right),$
(38b) $\displaystyle\lambda_{3x}$ $\displaystyle=$
$\displaystyle-\left[\left(\frac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-1\right]\left(\frac{1}{\omega_{3}}-\frac{1}{2}\right)u_{x},$
(38c) $\displaystyle\lambda_{3y}$ $\displaystyle=$
$\displaystyle-\left[\left(\frac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-r^{2}\right]\left(\frac{1}{\omega_{3}}-\frac{1}{2}\right)u_{y},$
(38d) $\displaystyle\theta_{4x}$ $\displaystyle=$
$\displaystyle-\Big{[}M_{4}+\left(3q^{2}c_{s*}^{2}-1\right)\rho\Big{]}\left(\frac{1}{\omega_{4}}-\frac{1}{2}\right),$
(39a) $\displaystyle\theta_{4y}$ $\displaystyle=$
$\displaystyle+\Big{[}N_{4}-\left(3q^{2}c_{s*}^{2}-r^{2}\right)\rho\Big{]}\left(\frac{1}{\omega_{4}}-\frac{1}{2}\right),$
(39b) $\displaystyle\lambda_{4x}$ $\displaystyle=$
$\displaystyle-\left[\left(\frac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-1\right]\left(\frac{1}{\omega_{4}}-\frac{1}{2}\right)u_{x},$
(39c) $\displaystyle\lambda_{4y}$ $\displaystyle=$
$\displaystyle+\left[\left(\frac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-r^{2}\right]\left(\frac{1}{\omega_{4}}-\frac{1}{2}\right)u_{y},$
(39d) $\displaystyle\theta_{5x}$ $\displaystyle=$
$\displaystyle-\left(\frac{1}{\gamma^{2}}-\frac{1}{\gamma}\right)\left(\frac{1}{\omega_{5}}-\frac{1}{2}\right)\rho
u_{x}u_{y},$ (40a) $\displaystyle\theta_{5y}$ $\displaystyle=$
$\displaystyle-\left(\frac{1}{\gamma^{2}}-\frac{1}{\gamma}\right)\left(\frac{1}{\omega_{5}}-\frac{1}{2}\right)\rho
u_{x}u_{y},$ (40b) $\displaystyle\lambda_{5x}$ $\displaystyle=$
$\displaystyle-\left(\frac{1}{\gamma}-1\right)\left(\frac{1}{\omega_{5}}-\frac{1}{2}\right)q^{2}c_{s*}^{2}u_{y},$
(40c) $\displaystyle\lambda_{5y}$ $\displaystyle=$
$\displaystyle-\left(\frac{1}{\gamma}-1\right)\left(\frac{1}{\omega_{5}}-\frac{1}{2}\right)q^{2}c_{s*}^{2}u_{x}.$
(40d)

These expressions (Eqs. (38)-(40)) together with Eqs. (29) and (30) are among
the main results of this work that contribute towards formulating a new
preconditioned LB approach on a rectangular lattice grid. The above choices
for the moment equilibria corrections, which depend on both the grid aspect
ratio $r$ and the preconditioning parameter $\gamma$, ensures that the
resulting algorithm using a rectangular lattice represents the preconditioned
NS equations with the shear viscosity $\nu$ and bulk viscosity $\xi$
satisfying the following relationships among the various model parameters:

$\nu=\gamma q^{2}c_{s*}^{2}\left(\frac{1}{\omega_{j}}-\frac{1}{2}\right)\Delta
t,\quad j=4,5,\;\;\quad\xi=\gamma
q^{2}c_{s*}^{2}\left(\frac{1}{\omega_{3}}-\frac{1}{2}\right)\Delta t,$ (41)

where the optimal value of $c_{s*}^{2}$ is $1/3$, and the emergent pressure
field $p$ is given by $p=\gamma q^{2}c_{s*}^{2}\rho$. We emphasize here that
the simple expressions given in Eq. (41) self-consistently parameterize the
transport coefficients in terms of $q$ which is given in Eq. (2) and maintains
desired numerical stability in LB simulations using rectangular lattice grids.
Unlike in previous works (see e.g., [13, 17, 47]), there is no need to rely on
trial and error to adjust the speed of sound when a rectangular lattice is
used and the grid aspect ratio is varied to any desired value.

###  3.6 Strain rate tensor components based on non-equilibrium moments

We will now show the diagonal components of the strain rate tensor
$\partial_{x}u_{x}$ and $\partial_{y}u_{y}$, which appear in the moment
equilibria corrections given in Eqs. (29), (30) and Eqs. (38)-(40) can be
computed locally via second-order non-equilibrium moments. First, using Eqs.
(33a) and (33b), and simplifying via Eq. (37), we obtain

$\displaystyle n_{3}^{(1)}$ $\displaystyle=$
$\displaystyle-\cfrac{2q^{2}c_{s*}^{2}}{\omega_{3}}\rho\left(\partial_{x}u_{x}+\partial_{y}u_{y}\right)+\cfrac{\omega_{3}}{2}E_{3g}+\cfrac{\omega_{3}}{2}E_{3s},$
$\displaystyle n_{4}^{(1)}$ $\displaystyle=$
$\displaystyle-\dfrac{2q^{2}c_{s*}^{2}}{\omega_{4}}\rho\left(\partial_{x}u_{x}-\partial_{y}u_{y}\right)+\dfrac{\omega_{4}}{2}E_{4g}+\dfrac{\omega_{4}}{2}E_{4s}$

Then, substituting for $E_{3s}$, $E_{4s}$, $E_{3g}$, and $E_{4g}$ using Eqs.
(26) and (27) in the last two equations and rearranging them leads to

$\displaystyle\left[-\dfrac{2q^{2}c_{s*}^{2}}{\omega_{3}}\rho+\dfrac{M_{3}}{2}+\dfrac{1}{2}(3q^{2}c_{s*}^{2}-1)\rho\right]\partial_{x}u_{x}+\left[-\dfrac{2q^{2}c_{s*}^{2}}{\omega_{3}}\rho+\dfrac{N_{3}}{2}+\dfrac{1}{2}(3q^{2}c_{s*}^{2}-r^{2})\rho\right]\partial_{y}u_{y}$
$\displaystyle=n_{3}^{(1)}-\dfrac{1}{2}\left[\left(\dfrac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-1\right]u_{x}\partial_{x}\rho-\dfrac{1}{2}\left[\left(\dfrac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-r^{2}\right]u_{y}\partial_{y}\rho.$
(42)
$\displaystyle\left[-\dfrac{2q^{2}c_{s*}^{2}}{\omega_{4}}\rho+\dfrac{M_{4}}{2}+\dfrac{1}{2}(3q^{2}c_{s*}^{2}-1)\rho\right]\partial_{x}u_{x}+\left[+\dfrac{2q^{2}c_{s*}^{2}}{\omega_{4}}\rho+\dfrac{N_{4}}{2}-\dfrac{1}{2}(3q^{2}c_{s*}^{2}-r^{2})\rho\right]\partial_{y}u_{y}$
$\displaystyle=n_{4}^{(1)}-\dfrac{1}{2}\left[\left(\dfrac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-1\right]u_{x}\partial_{x}\rho+\dfrac{1}{2}\left[\left(\dfrac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-r^{2}\right]u_{y}\partial_{y}\rho.$
(43)

Based on Eqs. (3.6) and (3.6), the required local expressions for the diagonal
components of the strain rate tensor $\partial_{x}u_{x}$ and
$\partial_{y}u_{y}$ can be obtained. In this regard, we first introduce the
following intermediate variables

$\displaystyle A$ $\displaystyle=$
$\displaystyle\cfrac{1}{2}\left[\left(\frac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-1\right]u_{x},\qquad
B=\cfrac{1}{2}\left[\left(\frac{2}{\gamma}+1\right)q^{2}c_{s*}^{2}-r^{2}\right]u_{y},$
(44a) $\displaystyle e_{3\rho}$ $\displaystyle=$
$\displaystyle-A\partial_{x}\rho-B\partial_{y}\rho,\qquad\qquad\qquad
e_{4\rho}=-A\partial_{x}\rho+B\partial_{y}\rho,$ (44b)

where the density gradients $\partial_{x}\rho$ and $\partial_{y}\rho$ may be
obtained via an isotropic finite-difference scheme. The non-equilibrium
moments $n_{3}^{(1)}$ and $n_{4}^{(1)}$ appearing in Eqs. (3.6) and (3.6) can
be computed using either raw moments or central moments as

$\displaystyle n_{3}^{(1)}$ $\displaystyle=$
$\displaystyle\left(k_{20}^{\prime}+k_{02}^{\prime}\right)-\left(k_{20}^{eq\prime}+k_{02}^{eq\prime}\right)=\left(k_{20}+k_{02}\right)-\left(k_{20}^{eq}+k_{02}^{eq}\right)$
$\displaystyle=$
$\displaystyle\left(k_{20}+k_{02}\right)-2q^{2}c_{s*}^{2}\rho-\left(\frac{1}{\gamma}-1\right)(u_{x}^{2}+u_{y}^{2})\rho,$
$\displaystyle n_{4}^{(1)}$ $\displaystyle=$
$\displaystyle\left(k_{20}^{\prime}-k_{02}^{\prime}\right)-\left(k_{20}^{eq\prime}-k_{02}^{eq\prime}\right)=\left(k_{20}-k_{02}\right)-\left(k_{20}^{eq}-k_{02}^{eq}\right)$
$\displaystyle=$
$\displaystyle\left(k_{20}-k_{02}\right)-\left(\frac{1}{\gamma}-1\right)(u_{x}^{2}-u_{y}^{2})\rho.$

Based on these considerations, we can then identify the right sides and the
left sides of Eqs. (3.6) and (3.6), respectively, conveniently by further
introducing the following additional intermediate variables

$\displaystyle R_{3}$ $\displaystyle=$ $\displaystyle
n_{3}^{(1)}+e_{3\rho}=k_{20}+k_{02}-2q^{2}c_{s*}^{2}\rho-\left(\frac{1}{\gamma}-1\right)(u_{x}^{2}+u_{y}^{2})\rho+e_{3\rho},$
(45a) $\displaystyle R_{4}$ $\displaystyle=$ $\displaystyle
n_{4}^{(1)}+e_{4\rho}=k_{20}-k_{02}-\left(\frac{1}{\gamma}-1\right)(u_{x}^{2}-u_{y}^{2})\rho+e_{4\rho},$
(45b)

and

$\displaystyle
C_{3x}=\left[-\frac{2q^{2}c_{s*}^{2}}{\omega_{3}}\rho+\frac{M_{3}}{2}+\frac{1}{2}(3q^{2}c_{s*}^{2}-1)\rho\right],\quad
C_{3y}=\left[-\frac{2q^{2}c_{s*}^{2}}{\omega_{3}}\rho+\frac{N_{3}}{2}+\frac{1}{2}(3q^{2}c_{s*}^{2}-r^{2})\rho\right],$
(46a) $\displaystyle
C_{4x}=\left[-\frac{2q^{2}c_{s*}^{2}}{\omega_{4}}\rho+\frac{M_{4}}{2}+\frac{1}{2}(3q^{2}c_{s*}^{2}-1)\rho\right],\quad
C_{4y}=\left[+\frac{2q^{2}c_{s*}^{2}}{\omega_{4}}\rho+\frac{N_{4}}{2}-\frac{1}{2}(3q^{2}c_{s*}^{2}-r^{2})\rho\right],$
(46b)

where $M_{3}$, $N_{3}$, $M_{4}$ and $N_{4}$ are given in Eq. (28). Then, Eqs.
(3.6) and (3.6) can be more compactly written as

$\displaystyle C_{3x}\partial_{x}u_{x}+C_{3y}\partial_{y}u_{y}$
$\displaystyle=$ $\displaystyle R_{3},$ (47a) $\displaystyle
C_{4x}\partial_{x}u_{x}+C_{4y}\partial_{y}u_{y}$ $\displaystyle=$
$\displaystyle R_{4}.$ (47b)

Solving the last two equations, we finally get the required local expressions
for the diagonal parts of the strain rate tensor as follows:

$\partial_{x}u_{x}=\frac{\left[C_{4y}R_{3}-C_{3y}R_{4}\right]}{\left[C_{3x}C_{4y}-C_{4x}C_{3y}\right]},\qquad\qquad\partial_{y}u_{y}=\frac{1}{C_{3y}}\left[R_{3}-C_{3x}\partial_{x}u_{x}\right].$
(48)

For completeness, we note that a similar relation for the off-diagonal
component ($\partial_{x}u_{y}+\partial_{y}u_{x}$) follows from combining Eqs.
(33c) and (37) and then simplifying via Eq. (27).

##  4 Preconditioned Rectangular Central Moment Lattice Boltzmann Method (PRC-
LBM)

In this section, we will present a robust and efficient implementation of a LB
algorithm on rectangular lattice grids for solving preconditioned NS equations
using and extending the results of the C-E analysis performed in the last
section. In this regard, we note that the effect of the moment basis
$\bm{\mathsf{T}}$ as defined in Eqs. (6) and (3.1) will be equivalently
utilized in our implementation in a more modular fashion so the LB schemes
based on the square lattice can be readily extended for utilizing rectangular
lattice grids along with preconditioning and the necessary corrections. This
involves using a simpler re-defined moment basis in conjunction with diagonal
scaling matrices based on the grid aspect ratio for performing the pre- and
post-collision transformations between the raw moments and distribution
functions, and segregate the evolution of the trace of the diagonal parts of
the second order moments from the others for achieving independent variations
of bulk and shear viscosities [7], along with accounting for the corrections
to eliminate the grid-anisotropy and non-GI truncation errors, only within the
collision step under moment relaxations (see Ref. [35]). In other words, the
linear combinations of moments as required are considered only for performing
the collision step and not for any mappings. This represents an improvement
over the implementation discussed over all the previous 2D rectangular LB
schemes for the solution of the NS equations, including our recent work [34],
and is consistent with our more recent 3D formulation [35]), but extended here
with a preconditioning strategy for convergence acceleration. Similar approach
based on a natural independent moment set without involving the mixed moments
has also been used in previous work [48, 49] by using a block diagonal
relaxation matrix [18] in the context of LB formulations using a square
lattice.

###  4.1 Reformulation of the Preconditioned Rectangular Raw Moment LBE

Thus, we first introduce a moment basis $\bm{\mathsf{Q}}$, which unlike
$\bm{\mathsf{T}}$ in Eqs. (6) and (3.1), does not contain any combinations of
the basis vectors, but only a set of bare basis vectors for the D2Q9 lattice:

$\displaystyle\bm{\mathsf{Q}}=\Big{[}\;\ket{1},\ket{e_{x}},\ket{e_{y}},\ket{e_{x}^{2}},\ket{e_{y}^{2}},\ket{{e_{x}}{e_{y}}},\ket{e_{x}^{2}e_{y}},\ket{e_{x}e_{y}^{2}},\ket{e_{x}^{2}e_{y}^{2}}\;\Big{]}^{{\dagger}},$
(49)

where $\ket{e_{x}}$, $\ket{e_{y}}$ and $\ket{1}$ are given in Eqs.(3a)-(3b)
and (4), respectively. Hence, $\bm{\mathsf{Q}}$ depends on the grid aspect
ratio $r$. We can relate this moment basis for a _rectangular_ lattice
$\bm{\mathsf{Q}}$ to an equivalent moment basis for a _square_ lattice
$\bm{\mathsf{P}}$ given by

$\displaystyle\bm{\mathsf{P}}=\Big{[}\;\ket{1},\ket{\bar{e}_{x}},\ket{\bar{e}_{y}},\ket{\bar{e}_{x}^{2}},\ket{\bar{e}_{y}^{2}},\ket{{\bar{e}_{x}}{\bar{e}_{y}}},\ket{\bar{e}_{x}^{2}\bar{e}_{y}},\ket{\bar{e}_{x}\bar{e}_{y}^{2}},\ket{\bar{e}_{x}^{2}\bar{e}_{y}^{2}}\;\Big{]}^{{\dagger}},$
(50)

where the particle velocity components of the square lattice
$\ket{\bar{e}_{x}}$ and $\ket{\bar{e}_{y}}$ are given as

$\displaystyle\ket{\bar{e}_{x}}$ $\displaystyle=$
$\displaystyle(0,1,0,-1,0,1,-1,-1,1)^{\dagger},$
$\displaystyle\ket{\bar{e}_{y}}$ $\displaystyle=$
$\displaystyle(0,0,1,0,-1,1,1,-1,-1)^{\dagger}.$

Evidently, the two moment basis matrices $\bm{\mathsf{Q}}$ and
$\bm{\mathsf{P}}$ can be readily related via a diagonal scaling matrix
$\bm{\mathsf{S}}$

$\bm{\mathsf{Q}}=\bm{\mathsf{S}}\bm{\mathsf{P}},$ (52)

where $\bm{\mathsf{S}}$ reads as

$\bm{\mathsf{S}}=\mbox{diag}{\begin{bmatrix}\;1&1&r&1&r^{2}&r&r&r^{2}&r^{2}\;\end{bmatrix}}.$
(53)

Importantly, from Eq. (52), the matrix inverse of $\bm{\mathsf{Q}}$ follows
directly from the inverse of $\bm{\mathsf{P}}$ for the square lattice, which
is quite straightforward to perform, and the inverse of the scaling matrix
$\bm{\mathsf{S}}$ using

$\bm{\mathsf{Q}}^{-1}=\bm{\mathsf{P}}^{-1}\bm{\mathsf{S}}^{-1},$ (54)

where $\bm{\mathsf{S}}^{-1}$ is obtained from Eq. (53), which being a diagonal
matrix, by simply taking the reciprocal of each of its elements, i.e.,

$\bm{\mathsf{S}}^{-1}=\mbox{diag}{\begin{bmatrix}\;1&1&r^{-1}&1&r^{-2}&r^{-1}&r^{-1}&r^{-2}&r^{-2}\;\end{bmatrix}}.$
(55)

In other words, $\bm{\mathsf{Q}}^{-1}$ for the rectangular lattice is easy to
perform knowing the corresponding $\bm{\mathsf{P}}^{-1}$ by appropriate
scalings of the latter’s elements based on the grid aspect ratio $r$. By
contrast, since $\bm{\mathsf{T}}$ is defined using combinations of the basis
vectors, its inverse $\bm{\mathsf{T}}^{-1}$, involves cumbersome expressions
with parameterizations based on $r$. This fact confers a significant advantage
of using $\bm{\mathsf{Q}}$ (and its inverse) rather than $\bm{\mathsf{T}}$ in
performing mappings between moments and distribution functions [35] and is
thus adopted in designing our LB algorithm in what follows. However, as
mentioned earlier, the effect of such combinations should still be accounted
for in the evolution of the moments, which we accomplish by formally
introducing a matrix $\bm{\mathsf{B}}$ in

$\bm{\mathsf{T}}=\bm{\mathsf{B}}\bm{\mathsf{Q}}.$ (56)

Thus, $\bm{\mathsf{B}}$ expresses the combinations of the moments (for the
second order components $\ket{e_{x}^{2}+e_{y}^{2}}$ and
$\ket{e_{x}^{2}-e_{y}^{2}}$ in the case of the D2Q9 lattice), which will be
effectively introduced in the LB scheme in the evolution of the corresponding
combinations of moments under collision and not for mappings. For this
purpose, using the moment basis defined by $\bm{\mathsf{Q}}$, we can then
define a set of bare moments $\mathbf{m}$ from the distribution functions
$\mathbf{f}$ (and vice versa) using

$\displaystyle\mathbf{m}=\bm{\mathsf{Q}}\mathbf{f},\qquad\mathbf{f}={\bm{\mathsf{Q}}}^{-1}\mathbf{m},$
(57)

where $\mathbf{m}$ is given by

$\mathbf{m}=\left(k_{00}^{\prime},k_{10}^{\prime},k_{01}^{\prime},k_{20}^{\prime},k_{02}^{\prime},k_{11}^{\prime},k_{21}^{\prime},k_{12}^{\prime},k_{22}^{\prime}\right)^{{\dagger}},$
(58)

and similarly for the sets of raw moment equilibria and the source terms,
respectively, via $\mathbf{m}^{eq}=\bm{\mathsf{Q}}\mathbf{f}^{eq}$ and
$\mathbf{\Phi}=\bm{\mathsf{Q}}\mathbf{S}$ as

$\displaystyle\mathbf{m}^{eq}$ $\displaystyle=$
$\displaystyle\left(k_{00}^{eq\prime},k_{10}^{eq\prime},k_{01}^{eq\prime},k_{20}^{eq\prime},k_{02}^{eq\prime},k_{11}^{eq\prime},k_{21}^{eq\prime},k_{12}^{eq\prime},k_{22}^{eq\prime}\right)^{{\dagger}},$
(59a) $\displaystyle\mathbf{\Phi}$ $\displaystyle=$
$\displaystyle\left(\sigma_{00}^{\prime},\sigma_{10}^{\prime},\sigma_{01}^{\prime},\sigma_{20}^{\prime},\sigma_{02}^{\prime},\sigma_{11}^{\prime},\sigma_{21}^{\prime},\sigma_{12}^{\prime},\sigma_{22}^{\prime}\right)^{{\dagger}}.$
(59b)

These represent the simpler bare moments versions of those given in Eq. (11)
for the combined moments of various quantities (which followed from Eq. (12)).

From the above developments by exploiting the properties of the various
matrices introduced and rearranging, the preconditioned raw moment based MRT-
LBE in Eq. (13) can be rewritten in the following equivalent form (see Ref.
[35] for details):

$\mathbf{f}(\bm{x}+\mathbf{e}\Delta t,t+\Delta
t)=\bm{\mathsf{P}}^{-1}\bm{\mathsf{S}}^{-1}\Big{[}\mathbf{m}+\bm{\mathsf{B}}^{-1}\bm{\mathsf{\Lambda}}\;\left(\;\bm{\mathsf{B}}\mathbf{m}^{eq}-\bm{\mathsf{B}}\mathbf{m}\;\right)+\bm{\mathsf{B}}^{-1}\left(\bm{\mathsf{I}}-\frac{\bm{\mathsf{\Lambda}}}{2}\right)\bm{\mathsf{B}}\mathbf{\Phi}\Delta
t\Big{]}.$ (60)

This equation (Eq. (60)) can be more conveniently represented by splitting
them in the form of the following sequence of sub-steps that are amenable for
implementation:

$\displaystyle\mathbf{m}$ $\displaystyle=$
$\displaystyle\bm{\mathsf{S}}\bm{\mathsf{P}}\mathbf{f},$
$\displaystyle\tilde{\mathbf{m}}$ $\displaystyle=$
$\displaystyle\mathbf{m}+\bm{\mathsf{B}}^{-1}\left\\{\bm{\mathsf{\Lambda}}\;\left(\;\bm{\mathsf{B}}\mathbf{m}^{eq}-\bm{\mathsf{B}}\mathbf{m}\;\right)+\left(\bm{\mathsf{I}}-\frac{\bm{\mathsf{\Lambda}}}{2}\right)\bm{\mathsf{B}}\mathbf{\Phi}\Delta
t\right\\},$ $\displaystyle\tilde{\mathbf{f}}(\bm{x},t)$ $\displaystyle=$
$\displaystyle\bm{\mathsf{P}}^{-1}\bm{\mathsf{S}}^{-1}\tilde{\mathbf{m}},$
$\displaystyle\mathbf{f}(\bm{x}+\mathbf{e}\Delta t,t+\Delta t)$
$\displaystyle=$ $\displaystyle\tilde{\mathbf{f}}(\bm{x},t).$ (61)

Here, we emphasize that $\bm{\mathsf{P}}$ and $\bm{\mathsf{P}}^{-1}$ perform
transformations between the distribution functions and raw moments in a way as
done for the usual _square_ lattice using the non-orthogonal moment basis,
$\bm{\mathsf{S}}$ and $\bm{\mathsf{S}}^{-1}$ reflect the simple scalings of
the raw moments by factors based on grid aspect ratio and the order of the
moment before and after collision, respectively, and $\bm{\mathsf{B}}$ and
$\bm{\mathsf{B}}^{-1}$ represent combining moments prior to their relaxations
under collision with the addition of the source terms, and their subsequent
segregation, respectively. Equation (61) expresses the preconditioned
rectangular LBM based on _raw moments_. As such, this scheme, by including
both the numerical enhancement features, viz., preconditioning and rectangular
lattice grids together, even in the context of raw moments, is new and
suitable for implementation. Nevertheless, a number of prior studies (see
e.g., [21, 28, 29, 50, 31, 31, 32, 33]), including those involving
rectangular/cuboid lattices [34, 35], have demonstrated that constructing LB
schemes involving the relaxations of _central moments_ under collision offer
significant improvements in numerical stability over those based on raw
moments. Hence, in this work we will only implement and perform a numerical
study on the generalization of the above developments to central moments,
viz., the preconditioned rectangular central moment LBE, which will be
discussed next and followed by a summary of its algorithmic steps.

###  4.2 Formulation of the Preconditioned Rectangular Central Moment LBE

For this purpose, we will utilize the independently supported bare central
moments defined in Eq. (10) for the D2Q9 lattice and collect them in the form
of the following vectors:

$\displaystyle\mathbf{m}^{c}$ $\displaystyle=$
$\displaystyle\left(k_{00},k_{10},k_{01},k_{20},k_{02},k_{11},k_{21},k_{12},k_{22}\right)^{{\dagger}},$
(62a) $\displaystyle\mathbf{m}^{c,eq}$ $\displaystyle=$
$\displaystyle\left(k_{00}^{eq},k_{10}^{eq},k_{01}^{eq},k_{20}^{eq},k_{02}^{eq},k_{11}^{eq},k_{21}^{eq},k_{12}^{eq},k_{22}^{eq}\right)^{{\dagger}},$
(62b) $\displaystyle\mathbf{\Phi}^{c}$ $\displaystyle=$
$\displaystyle\left(\sigma_{00},\sigma_{10},\sigma_{01},\sigma_{20},\sigma_{02},\sigma_{11},\sigma_{21},\sigma_{12},\sigma_{22}\right)^{{\dagger}},$
(62c)

Now, the raw moments defined in Eq. (9) can be related to the central moments
in Eq. (10) via straightforward binomial expansions involving the former in
combinations with monomials of the fluid velocity components at different
order (of the form $u_{x}^{p}u_{y}^{q}$). Thus, the mappings from the raw
moments to central moments (and vice versa) can be formally expressed as

$\displaystyle\mathbf{m}^{c}=\bm{\mathsf{F}}\mathbf{m},\qquad\mathbf{m}={\bm{\mathsf{F}}}^{-1}\mathbf{m}^{c},$
(63)

where $\bm{\mathsf{F}}$ is referred to as the frame transformation matrix
reflecting the binomial transforms of moments at different orders supported by
the D2Q9 lattice. Such a formulation to represent the transformations between
the raw moments and central moments for the complete set supported by the
lattice in the form of a shift matrix was first introduced by Fei and Luo in
[24, 25]. It is given by

$\bm{\mathsf{F}}=\begin{bmatrix}1&0&0&0&0&0&0&0&0\\\\[4.0pt]
-u_{x}&1&0&0&0&0&0&0&0\\\\[4.0pt] -u_{y}&0&1&0&0&0&0&0&0\\\\[4.0pt]
u_{x}^{2}+u_{y}^{2}&-2u_{x}&-2u_{y}&1&0&0&0&0&0\\\\[4.0pt]
u_{x}^{2}-u_{y}^{2}&-2u_{x}&2u_{y}&0&1&0&0&0&0\\\\[4.0pt]
u_{x}u_{y}&-u_{y}&-u_{x}&0&0&1&0&0&0\\\\[4.0pt]
-u_{x}^{2}u_{y}&2u_{x}u_{y}&u_{x}^{2}&-\frac{1}{2}u_{y}&-\frac{1}{2}u_{y}&-2u_{x}&1&0&0\\\\[4.0pt]
-u_{x}u_{y}^{2}&u_{y}^{2}&2u_{x}u_{y}&-\frac{1}{2}u_{x}&\frac{1}{2}u_{x}&-2u_{y}&0&1&0\\\\[4.0pt]
u_{x}^{2}u_{y}^{2}&-2u_{x}u_{y}^{2}&-2u_{x}^{2}u_{y}&\frac{1}{2}(u_{x}^{2}+u_{y}^{2})&\frac{1}{2}(u_{y}^{2}-u_{x}^{2})&4u_{x}u_{y}&-2u_{y}&-2u_{x}&1\\\
.\end{bmatrix}$ (64)

Also, as noted in Ref. [34], its inverse ${\bm{\mathsf{F}}}^{-1}$ can be read
off directly from the elements of
${\bm{\mathsf{F}}}={\bm{\mathsf{F}}}(u_{x},u_{y})$ with minor changes by
exploiting the following property that exists for such transforms:
${\bm{\mathsf{F}}}^{-1}={\bm{\mathsf{F}}}(-u_{x},-u_{y})$. Thus, and thus
naturally both of them are lower triangular matrices. Then, by an analogy with
Eq. (60), we can write the following preconditioned rectangular central moment
LBE by involving relaxations of central moments $\mathbf{m}^{c}$ under
collision (rather than raw moments $\mathbf{m}$) and including the additional
transforms between them and the raw moments (via $\bm{\mathsf{F}}$ and
$\bm{\mathsf{F}}^{-1}$) [35]:

$\mathbf{f}(\bm{x}+\mathbf{e}\Delta t,t+\Delta
t)=\bm{\mathsf{P}}^{-1}\bm{\mathsf{S}}^{-1}\bm{\mathsf{F}}^{-1}\Big{[}\mathbf{m}^{c}+\bm{\mathsf{B}}^{-1}\bm{\mathsf{\Lambda}}\;\left(\;\bm{\mathsf{B}}\mathbf{m}^{eq,c}-\bm{\mathsf{B}}\mathbf{m}^{c}\;\right)+\bm{\mathsf{B}}^{-1}\left(\bm{\mathsf{I}}-\frac{\bm{\mathsf{\Lambda}}}{2}\right)\bm{\mathsf{B}}\mathbf{\Phi}^{c}\Delta
t\Big{]}.$ (65)

It may be noted that this Eq. (65) has some similarities with those presented
by Luo and collaborators [24, 25], where they presented the cascaded central
moment LB method in a generalized MRT framework on a square lattice. Our
notations follow from those presented in an earlier work of the second author
[19] for the frame transformation matrix $\bm{\mathsf{F}}$. This matrix as
shown in Eq. (64) is identical to the shift matrix used in [24, 25].
Nevertheless, there are some key differences: Eq. (65) also involves the
forward and inverse scaling transformations (via diagonal matrices) related to
the grid aspect ratio $r$ to accommodate the use of a rectangular lattice in a
modular fashion. Moreover, for the collision step, Refs. [24, 25] combine the
relaxation parameters of the second order moments similar to that in [18]. By
contrast, here, to execute the collision step, the second order moments are
combined prior to collision, which are then relaxed at independent rates to
their equilibria (with appropriate corrections based on $r$ and $\gamma$) and
then segregated post collision. Such a strategy for performing collision was
presented by Geier _et al._ [7], and present work can be considered as an
extension of such an approach for performing flow simulations on rectangular
lattices with preconditioning. Also, it should be noted that not all collision
models admit interpretations based on matrices. For example, the highly
sophisticated and nonlinear cumulant LB scheme cannot be represented in the
form of matrices. Hence, the algorithms (including the special cases for raw
moments and central moments) presented in [7] are given only as series of
substeps and no matrices are utilized in this regard. Thus, to maintain
generality of our approach, in what follows, we will represent our PRC-LBM in
the form of a sequence of operations that conveys the essence of Eq. (65).
This last equation (Eq. (65)) can then be more conveniently split up into
various sub-steps, which then results into the following preconditioned
rectangular central moment LBM or PRC-LBM:

$\displaystyle\mathbf{m}^{c}$ $\displaystyle=$
$\displaystyle\bm{\mathsf{F}}\bm{\mathsf{S}}\bm{\mathsf{P}}\mathbf{f},$
$\displaystyle\tilde{\mathbf{m}}^{c}$ $\displaystyle=$
$\displaystyle\mathbf{m}^{c}+\bm{\mathsf{B}}^{-1}\left\\{\bm{\mathsf{\Lambda}}\;\left(\;\bm{\mathsf{B}}\mathbf{m}^{c,eq}-\bm{\mathsf{B}}\mathbf{m}^{c}\;\right)+\left(\bm{\mathsf{I}}-\frac{\bm{\mathsf{\Lambda}}}{2}\right)\bm{\mathsf{B}}\mathbf{\Phi}^{c}\Delta
t\right\\},$ $\displaystyle\tilde{\mathbf{f}}(\bm{x},t)$ $\displaystyle=$
$\displaystyle\bm{\mathsf{P}}^{-1}\bm{\mathsf{S}}^{-1}\bm{\mathsf{F}}^{-1}\tilde{\mathbf{m}}^{c},$
$\displaystyle\mathbf{f}(\bm{x}+\mathbf{e}\Delta t,t+\Delta t)$
$\displaystyle=$ $\displaystyle\tilde{\mathbf{f}}(\bm{x},t).$ (66)

This PRC-LBM is a numerically more robust formulation than its raw moment
counterpart given earlier in Eq. (61). The algorithmic details of the PRC-LBM
to facilitate its implementation are discussed in A.

##  5 Results and discussion

We will now discuss some case studies based on the PRC-LB algorithm for
simulations of shear flows at various characteristic parameters that show its
numerical validation against certain benchmark problems and the significant
advantages of combining the rectangular lattice grid and preconditioning over
the LB scheme based on the square lattice and without preconditioning. In this
regard, as noted at the end of A, in what follows, the no-slip boundary
condition for the moving walls, which generate shear flows, are accounted for
via the momentum augmented half-way bounce back approach and including the
parametrization for the grid aspect ratio for the rectangular lattice given in
our recent work [34].

###  5.1 2D Shear Flows in Lid-driven Square Cavity using PRC-LBM: Validation

First, we will assess the accuracy of the PRC-LBM for the simulation of the
classical flow within a _square_ cavity of side $H$, whose top surface moves
at a constant velocity $U$ setting up flow patterns that depend on the
Reynolds number given by $\mbox{Re}=UH/\nu$. We performed simulations at
Reynolds numbers of $\mbox{Re}=100,1000$ and $3200$ at a fixed Mach number
$\mbox{Ma}=0.05$. A rectangular lattice with grid aspect ratio of $r=0.5$
using a grid resolution of $N_{x}\times N_{y}=200\times 400$ is employed by
setting the preconditioning parameter $\gamma=0.1$ in our algorithm given in
A. The numerical results of the horizontal and vertical components of the
velocity profiles along the centerlines of the cavity predicted by the PRC-LBM
at the above three choices of Re are compared against the benchmark numerical
solutions given by Ghia _et al._ [51] in Fig. 2. It is evident that the PRC-
LBM results are in very good agreement with the benchmark data.

(a) $\mbox{Re}=100$

(b) $\mbox{Re}=1000$

(c) $\mbox{Re}=3200$

Figure 2: The components of the velocity profiles $u(y)$ and $v(x)$ along the
vertical and horizontal centerlines of a square cavity, i.e., $x=H/2$ and
$y=H/2$ respectively, computed using the PRC-LBM on a rectangular lattice grid
of aspect ratio of $r=0.5$ with the preconditioning parameter $\gamma=0.1$ at
Ma$=0.05$ for Reynolds numbers of (a) Re$=100$, (b) Re$=1000$, (c) Re$=3200$
and compared with the benchmark numerical solutions of Ref. [51] (symbols).

###  5.2 2D Shear Flows in Lid-driven Shallow and Deep Cavities using PRC-LBM:
Validation and Convergence Acceleration

Next, we will demonstrate the accuracy and computational advantages of using
the PRC-LBM for computing anisotropic and inhomogeneous shear flows inside
_rectangular_ cavities of length $L$ and height $H$, and characterized by the
geometric aspect ratio $\mbox{AR}=H/L$. As shown in Fig 3, we consider two
cases: (a) shallow cavity with aspect ratio $\mbox{AR}<1$ and (b) deep cavity
with aspect ratio $\mbox{AR}>1$. In each case, the flow is set up by the
motion of the upper lid with a velocity $U$ in the positive $x$ direction
which generates vortices that are different in size and shape based on the
Reynolds number specified by $\mbox{Re}=UH/\nu$. The confinement effect
characterized by the aspect ratio AR results in the characteristic flow scales
or the spatial gradients in velocities that can be different in different
coordinate directions, which can be more naturally and efficiently resolved by
using a rectangular lattice grid. In our previous work, we illustrated the
benefits of using the rectangular lattice over that based on the square
lattice for simulating such flows within shallow cavities with the use of
fewer grid nodes for the former when compared to the latter [34]. In the
current study, we aim to show further improvements of utilizing
preconditioning with rectangular lattice grids for convergence acceleration of
flows to their steady states, resulting in dramatic cumulative advantages of
simulating such flows using the square lattice and without preconditioning.

(a) Shallow cavity ($\mbox{AR}<1$)

(b) Deep cavity ($\mbox{AR}>1$)

Figure 3: Schematic arrangements of the flows inside a 2D (a) shallow and (b)
deep cavities of dimensions $L\times H$.

In this regard, first, simulations of flow inside a shallow rectangular cavity
of aspect ratio $\mbox{AR}=0.25$ at a Reynolds number $\mbox{Re}=100$ and Mach
number $\mbox{Ma}=0.06$ are performed using rectangular lattice with grid
aspect ratio $r=\Delta y/\Delta x=0.2$. If $N_{x}$ and $N_{y}$ are the number
of grid nodes resolving the cavity in $x$ and $y$ directions, respectively,
the grid spacings in the respective directions satisfy $\Delta x=L/N_{x}$ and
$\Delta y=H/N_{y}$, or $r=\mbox{AR}N_{x}/N_{y}$. Thus,
$N_{x}=(r/\mbox{AR})N_{y}$ in the case of the rectangular lattice and
$N_{x}=(1/\mbox{AR})N_{y}$ for the square lattice, where $r=1$. Choosing
$N_{y}=125$, this results in $100\times 125$ as the total number of grid nodes
for the rectangular lattice, while taking even somewhat smaller value of
$N_{y}=100$, however, requires a total $400\times 100$ grid nodes for the
square lattice case, which is significantly more when compared to the former.
Thus, if the computed results in each case are in agrement with one another,
this, by itself, is a saving in the memory storage and computational cost by a
factor of over $3$ in using the rectangular lattice. Then, the effect of
different levels of preconditioning in PRC-LBM is studied by considering
$\gamma=1.0,0.5,0.1$ and $0.05$. Figures 4a and 4b show comparisons of the
centerline velocity components $u$ and $v$ computed using PRC-LBM with $r=0.2$
and $100\times 125$ for different $\gamma$ with the results obtained using the
square lattice ($r=1$) with $400\times 100$ and without preconditioning
($\gamma=1$). It can be seen the PRC-LBM results are in remarkably good
agreement for the entire range of the choice of $\gamma$, when compared to the
corresponding velocity profiles for the square lattice case.

(a) $\mbox{AR}=0.25$, Re=100, $r=0.2$

(b) $\mbox{AR}=0.25$, Re=100, $r=0.2$

Figure 4: The velocity profiles along the centerlines of a shallow rectangular
cavity of aspect ratio $\mbox{AR}=0.25$ at a Reynolds number $\mbox{Re}=100$
and Mach number $\mbox{Ma}=0.06$ computed using the PRC-LBM with a grid
resolution $100\times 125$ using the rectangular lattice of grid aspect ratio
$r=0.2$ with different levels of preconditioning, i.e., $\gamma=1.0,0.5,0.1$
and $0.05$, and compared with the results obtained using a square lattice
($r=1.0$) with a grid resolution of $400\times 100$ at $\gamma=1.0$. (a) $u$
component of the velocity along the vertical centerline, and (b) $v$ component
of the velocity along the vertical centerline.

Then, in order to verify the benefits of utilizing the preconditioning
procedure with the rectangular lattice, we study the convergence histories to
the steady state in using the PRC-LBM at $r=0.2$ for two values of the Mach
number and different $\gamma$ in each case. Figure 5a shows the convergence
histories at $\mbox{Ma}=0.06$ with $\gamma=1,0.5,0.1$, and $0.08$ and Fig. 5b
shows that at smaller Mach number $\mbox{Ma}=0.01$ with $\gamma=1,0.5,0.1$,
and $0.05$, where the residual global error of the $u$ velocity component is
estimated under the second norm as $||u(t+20)-u(t)||_{2}$. At
$\mbox{Ma}=0.06$, it takes about $610,000$ steps to reach the steady state
without preconditioning ($\gamma=1$), while the PRC-LBM with $\gamma=0.08$
requires only about $24,000$ steps to reach similar residual error as the
previous case, leading to a dramatic reduction in the number of steps for
convergence by a factor of about $25$ in this case. On the other hand, at a
lower $\mbox{Ma}=0.01$, the PRC-LBM takes about $5,000,000$ without
preconditioning, while only about $97,000$ with preconditioning (using
$\gamma=0.05$), with an even larger improvement corresponding to a reduction
factor of about $51$. Clearly, at lower Mach numbers, the disparities between
the flow speed and the sound speed are larger, and the associated higher
stiffness is respectively alleviated to a greater degree with preconditioning,
which is consistent with previous investigations on square lattice grids (see
e.g., [39, 43]). Noting that we have already reduced the computational costs
by involving the rectangular lattice grids when compared to that using the
square lattice, preconditioning the rectangular central moment LBM provides a
further, i.e., cumulative improvement in solving steady state flow problems
more efficiently. However, it should be noted that while using smaller values
of $\gamma$ does favor faster convergence speed, its smallest possible value
is limited by the numerical stability considerations (see e.g., [39, 40, 41]).
The optimal value of the level of preconditioning is a compromise between
convergence rate and stability. Typically, the minimum possible $\gamma$ is
found to be proportional to the Mach number used in simulations.

(a) $\mbox{Ma}=0.06$, $r=0.2$

(b) $\mbox{Ma}=0.01$, $r=0.2$

Figure 5: Convergence histories to the steady state for simulations of flows
within a shallow rectangular cavity of aspect ratio $\mbox{AR}=0.25$ at a
Reynolds number Re$=100$ using the PRC-LBM with a grid resolution of
$100\times 125$ using a rectangular lattice of grid aspect ratio $r=0.2$ with
different values of the preconditioned parameter $\gamma$ at (a) Mach number
Ma$=0.06$, and (b) Mach number Ma$=0.01$.

Next, we simulate the flow inside a deep cavity ($H>L$) as shown in Fig. 3b by
considering $\mbox{AR}=2$ at $\mbox{Re}=100$ and $\mbox{Ma}=0.06$ computed
using the PRC-LBM using $r=1.6$. Choosing $N_{x}=100$ for $r=1.6$, this leads
to $N_{y}=125$ for the rectangular lattice case, while for the square lattice
case, with $N_{x}=100$, we need $N_{y}=200$, i.e., more number of grid nodes
in the $y$ direction. Results shown in Fig. 6 present comparisons of the
centerline profiles of the components of the velocity across the deep cavity
computed using the PRC-LBM with $\gamma=1.0,0.5,0.1$ and $0.05$ using
$100\times 125$ rectangular grids with $r=1.6$ against those based on the
square lattice ($r=1$) with $100\times 200$ grid nodes and $\gamma=1.0$. As in
the shallow cavity case, it can be seen that the preconditioned central moment
LBM for all possible choice of $\gamma$ and with fewer number of grid nodes
delivers solutions that are in very good agreement with those based on the
square lattice.

(a) $AR=2$, $\mbox{Re}=100$, $r=1.6$

(b) $AR=2$, Re=100, $r=1.6$

Figure 6: The velocity profiles along the centerlines of a deep rectangular
cavity of aspect ratio $\mbox{AR}=2$ at a Reynolds number $\mbox{Re}=100$ and
Mach number $\mbox{Ma}=0.06$ computed using the PRC-LBM with a grid resolution
$100\times 125$ using the rectangular lattice of grid aspect ratio $r=1.6$
with different levels of preconditioning, i.e., $\gamma=1.0,0.5,0.1$ and
$0.05$, and compared with the results obtained using a square lattice
($r=1.0$) with a grid resolution of $100\times 200$ at $\gamma=1.0$. (a) $u$
component of the velocity along the vertical centerline, and (b) $v$ component
of the velocity along the vertical centerline.

Moreover, the convergence histories presented in Fig. 7 for deep cavity flow
simulations at $\mbox{AR}=2$ and Re$=100$ using the PRC-LBM with $r=1.6$ with
various levels of preconditioning again show a significant reduction in the
number of steps for convergence – for example, by a factor of $14$ with
$\gamma=0.08$ for $\mbox{Ma}=0.06$, and a factor of $23$ with $\gamma=0.05$
for $\mbox{Ma}=0.02$ when compared to the corresponding cases without
preconditioning. For the latter case with $\mbox{Ma}=0.02$, it may be noted
that the choice $\gamma=0.05$ is almost at the threshold of its smallest
possible value dictated by stability considerations. As a result, it may be
leading to a transient effect with the convergence history temporarily
crossing over that for the case $\gamma=0.1$ when the residual error is
O($10^{-12}$). Nevertheless, when the residual error drops further down to
O($10^{-13}$) or smaller, as expected, the convergence histories show that the
case $\gamma=0.05$ result in faster convergence to the steady state when
compared to $\gamma=0.1$.

(a) $\mbox{Ma}=0.06$, $r=1.6$

(b) $\mbox{Ma}=0.02$, $r=1.6$

Figure 7: Convergence histories to the steady state for simulations of flows
within a deep rectangular cavity of aspect ratio AR$=2$ at a Reynolds number
Re$=100$ using the PRC-LBM with a grid resolution of $100\times 125$ using a
rectangular lattice of grid aspect ratio $r=1.6$ with different values of the
preconditioned parameter $\gamma$ at (a) Mach number Ma$=0.06$, and (b) Mach
number Ma$=0.02$.

As a last case study, we perform an investigation on the efficacy of the PRC-
LBM for simulations of flows inside a deep cavity ($\mbox{AR}=2$) at a higher
Reynolds number of $\mbox{Re}=1000$ at $\mbox{Ma}=0.06$ using a grid
resolution of $N_{x}\times N_{y}=100\times 150$ corresponding to the grid
aspect ratio of $r=1.33$. The results of the centerline velocity component
profiles obtained using the PRC-LBM at $\gamma=1.0,0.75,0.5$ and $0.1$ are
reported in Fig. 8 and compared against the square lattice results using
$100\times 200$ grid nodes at $\gamma=1.0$. Evidently, the rectangular LB
scheme with all possible choices of $\gamma$ yields solutions that are again
in very good agreement with those based on the square lattice. Moreover, these
results are further corroborated by the plots of the streamline contours at
two Reynolds numbers $\mbox{Re}=100$ and $\mbox{Re}=1000$ simulated using the
square lattice ($r=1$) and the rectangular lattice ($r=1.6$ for
$\mbox{Re}=100$ and $r=1.33$ for $\mbox{Re}=1000$) and presented in Fig. 9,
which show the ability of the PRC-LBM to compute the flow patterns accurately.

(a) $AR=2$, $\mbox{Re}=1000$, $r=1.33$

(b) $AR=2$, $\mbox{Re}=1000$, $r=1.33$

Figure 8: The velocity profiles along the centerlines of a deep rectangular
cavity of aspect ratio $\mbox{AR}=2$ at a Reynolds number $\mbox{Re}=1000$ and
Mach number $\mbox{Ma}=0.06$ computed using the PRC-LBM with a grid resolution
$100\times 150$ using the rectangular lattice of grid aspect ratio $r=1.33$
with different levels of preconditioning, i.e., $\gamma=1.0,0.75,0.5$ and
$0.1$, and compared with the results obtained using a square lattice ($r=1.0$)
with a grid resolution of $100\times 200$ at $\gamma=1.0$. (a) $u$ component
of the velocity along the vertical centerline, and (b) $v$ component of the
velocity along the vertical centerline.

(a) $\mbox{Re}=100$, $AR=2$

(b) $\mbox{Re}=1000$, $AR=2$

Figure 9: Streamline contours of the flow field in a 2D rectangular deep
cavity of aspect ratio $\mbox{AR}=2$ computed using the PRC-LBM with
$\gamma=0.1$ on a rectangular lattice at (a) Re=100 using $r=1.6$, and (b)
Re=1000 using $r=1.33$ and, in each case, compared with the results of the
non-preconditioned LBM using the square lattice ($\gamma=1.0$ and $r=1$).

Finally, the convergence histories presented in Fig. 10 for deep cavity flow
simulations using a rectangular lattice grid ($r=2$) at $\mbox{Re}=1000$ and
$\mbox{Ma}=0.06$ show that it takes about $11,154,000$ steps to reach the
steady state without preconditioning, while the PRC-LBM with preconditioning
($\gamma=0.1$) requires significantly fewer steps of about $583,000$ to reach
similar residual errors, delivering an improvement by a factor of about $19$.
In general, the higher $\mbox{Re}=1000$ case takes longer to ready steady
state when compared to the case at lower $\mbox{Re}=100$ (shown in Fig. 7)
since the former is set up by reducing the fluid viscosity when compared to
the latter, resulting in a slower diffusion of momentum and its convergence.
However, thanks to preconditioning, even at higher Re, the rectangular central
moment LBM is able to achieve substantial savings in the overall computational
effort.

Figure 10: Convergence histories to the steady states for simulations of flows
within a deep rectangular cavity of aspect ratio $\mbox{AR}=2$ at a Reynolds
number $\mbox{Re}=1000$ using the PRC-LBM with a grid resolution of $100\times
150$ using a rectangular lattice of grid aspect ratio $r=1.33$ with different
values of the preconditioned parameter $\gamma$ at Mach number
$\mbox{Ma}=0.06$.

##  6 Comparisons between Preconditioned Rectangular LB Formulations based on
Raw Moments and Central Moments

As noted in the introduction, no other formulation other than the current work
that combines both preconditioning and rectangular lattice grids are available
in the literature. Nevertheless, it should be noted that the derivation
presented in Sec. 3 can be utilized to construct a preconditioned rectangular
_raw_ moment LBM (referred to as the PRNR-LBM in what follows) by performing
the collision step in terms of relaxations involving the raw moments using the
corrections to the equilibria based on $\gamma$ and $r$ given in Sec. 3.5. The
implementation of this strategy is summarized in Eq. (61). As such this PRNR-
LBM can be regarded as a special case of the PRC-LBM based on central moments.
Here, we note that the single-relaxation-time (SRT) formulations are well-
known to have serious deficiencies when compared other collision models. For
example, they are significantly less stable in simulating flows at relatively
low viscosities or large Reynolds numbers, even when compared to the
approaches based on raw moments, and it is quite cumbersome to combine the use
of both preconditioning and rectangular lattice grids. Given these issues,
they are not given further considerations in this work, and we will focus on
attention making comparisons between the PRC-LBM and PRNR-LBM in what follows.

First, we point out that provided a given collision model yields numerically
stable flow simulations, the convergence rate to the steady state is primarily
influenced by the choice of the preconditioning parameter, and the use of two
different collision models under otherwise similar conditions, such as the
same choice of model parameters, is expected to result in a similar
convergence speed. Figure 11 presents convergence histories for the PRC-LBM
and PRNR-LBM for the simulations of flow with a shallow cavity with aspect
ratio $\mbox{AR}=0.25$, $r=0.5$, $\mbox{Re}=100$, and $\mbox{Ma}=0.06$ for two
different values of $\gamma$ (equal to $0.1$ and $0.5$). For these choices,
both the collision models result in the same convergence rate to the steady
state for a fixed preconditioning parameter.

Figure 11: Convergence histories to the steady states for simulations of flows
within a shallow rectangular cavity of aspect ratio $\mbox{AR}=0.25$ at a
Reynolds number $\mbox{Re}=100$ using the PRC-LBM and PRNR-LBM with a grid
resolution of $100\times 50$ using a rectangular lattice of grid aspect ratio
$r=0.5$ with two different values of the preconditioned parameter $\gamma$ at
Mach number $\mbox{Ma}=0.06$.

On the other hand, we will now demonstrate the advantage of performing the
collision step in a frame of reference based on the local fluid velocity in
the PRC-LBM when compared to using the PRNR-LBM involving the rest or the
lattice frame of reference is related to improving robustness or numerical
stability of computations. In this regard, we investigate the maximum Reynolds
number achieved by PRC-LBM and PRNR-LBM for simulating flows within a shallow
cavity with an aspect ratio of $\mbox{AR}=0.25$ with a rectangular grid using
$r=0.25$ by maintaining the lid velocity at a relative small constant value of
$U=0.02$ and reducing the shear viscosity of the fluid to the smallest
possible value for which the computations remain numerically stable. The shear
viscosity is varied by changing the relaxation parameters associated with the
second order moments and the relaxation parameters for the higher order
moments are set to unity for simplicity for both the collision models. Two
different grid resolutions of $100\times 100$, $200\times 200$ are considered,
and in each case, three different choices of the preconditioning parameter
$\gamma=1.0$, $0.5$ and $0.2$ are used. The results of these stability tests
are presented in Fig. 12. Clearly, even for the relatively low lid velocity
considered, the preconditioned rectangular LBM using central moments is
consistently more stable when compared to the preconditioned rectangular LBM
using raw moments, especially at smaller $\gamma$. These results extend those
presented in Ref. [34] by including preconditioning effects. Further
improvements in numerical stability are possible when simulating shear flows
with larger characteristic velocities. Moreover, as emphasized by various
studies involving central moments in LBM (see e.g., [29, 31]), the additional
computational overhead of using central moments when compared to that of using
raw moments is relatively small, by about 15%, but with the benefit of stable
simulations at significantly higher Reynolds numbers using the former when
compared to the latter.

(a) $N_{x}\times N_{y}=100\times 100$, $AR=0.25$

(b) $N_{x}\times N_{y}=200\times 200$, $AR=0.25$

Figure 12: Comparison of the maximum Reynolds number Re for numerical
stability of PRC-LBM and PRNR-LBM for simulations of flows within a shallow
rectangular cavity of aspect ratio $\mbox{AR}=0.25$ with a fixed lid-velocity
of $0.02$ at different mesh resolutions with a grid aspect ratio of $r=0.25$
for three different choices of the preconditioning parameter ($\gamma=1.0$,
$0.5$ and $0.2$).

##  7 Summary and Conclusions

In this paper, we have developed a new LB algorithm based on central moments
and using involving a preconditioning strategy and a rectangular lattice grid,
viz., the PRC-LBM, for efficient simulations of inhomogeneous and anisotropic
flows. By including a preconditioning parameter $\gamma$ in its moment
equilibria and augmenting its second order components via corrections that
eliminate the anisotropy effects associated with the rectangular lattice
characterized by the grid aspect ratio $r$ and the non-Galilean invariant
velocity terms due to the aliasing effects on the D2Q9 lattice, it can
consistently represent the preconditioned Navier-Stokes equations. Such
corrections to the equilibria, obtained via a Chapman-Enskog analysis, are
related to the diagonal components of the velocity gradient tensor, which are
expressed in terms of the non-equilibrium moments to facilitate local
computations and their coefficients simultaneously depend on both $\gamma$ and
$r$. In the construction of the PRC-LBM, unlike the prior rectangular LB
formulations, we have used the natural, non-orthogonal moment basis with a
physically consistent parametrization of the speed of sound based on $r$, and
with a $\gamma$-adjusted equilibria obtained through a matching principle
based on the continuous Maxwell distribution. These result in simpler
correction terms for using rectangular lattice grids in conjunction with
preconditioning, and which, while used here in terms of central moments for
robustness, have general applicability for other collision models. Moreover,
our algorithmic implementation is modular in nature with a clear
interpretation based on special matrices, which also naturally extends to
three dimensions using a cuboid lattice in solving the preconditioned Navier-
Stokes equations. The PRC-LBM simulations of benchmark shear-driven flows
within shallow and deep cavities with significant geometric anisotropy at
various Reynolds numbers and Mach numbers and with different values of $r$ and
$\gamma$, validate the method for accuracy and show improvements in stability
when compared to another implementation based on raw moments. Furthermore, we
demonstrate significant reductions in the computational cost with the use of
PRC-LBM via convergence acceleration to the steady states and reduced memory
storage when compared to the corresponding LB scheme using the square lattice
and without involving preconditioning.

As final concluding remarks, we note the following. Given the ubiquity of
inhomogeneous and anisotropic flows, classical schemes for CFD invariably use
stretched grids that adapt to the local flow conditions and rarely utilize
uniform square or cubic grids. Thus, the use of square/cubic lattices in LB
algorithms is far from optimal in the use of overall computational resources
in performing such flow simulations efficiently. Previous efforts in
developing LB schemes based on rectangular lattices (e.g., [15, 16, 13, 17,
47]), which were generally based on orthogonal moment basis and non-optimal
equilibria, were significantly restrictive in terms of complexity of
implementation, limited stability ranges and cumbersome approach involved in
choosing the various model parameters. All these issues have been circumvented
in our present rectangular LB formulation involving a non-orthogonal moment
basis and the construction of equilibria based on a matching principle. It is
modular in construction in the sense that an LB code for the square lattice
can be readily extended to the rectangular lattice with some minor efforts
based on the former by making few simple changes. As shown in the appendix on
the algorithmic implementation, such changes include the use of pre-collision
and post-collision grid aspect ratio-based scalings of the raw moments and the
use of extended second order moment equilibria adjusted suitably based on the
grid aspect ratio to recover the Navier-Stokes equations. These simple
additional changes to the existing square lattice-based LB codes to
accommodate the rectangular grids have been numerically demonstrated to result
in savings by an order of magnitude or more in terms of computational cost and
memory when compared to that based on the square lattice (consistent with
Refs. [34, 35]). Furthermore, the use of preconditioning with rectangular
lattice grids resulted in significant convergence acceleration to the steady
states leading to further reduction in the overall computational efforts.
Hence, the PRC-LBM represents an efficient approach for flow simulations.

## Acknowledgements

The first author (EY) thanks the Department of Mechanical Engineering at the
University of Colorado Denver for financial support. The second author would
like to acknowledge the support of the US National Science Foundation (NSF)
under Grant CBET-1705630. The second author (KNP) would like to also thank the
NSF for support of the development of a computer cluster infrastructure
through the project “CC* Compute: Accelerating Science and Education by Campus
and Grid Computing” under Award 2019089.

##  Appendix A Algorithmic Details of the PRC-LBM

The PRC-LB algorithm consists of the following sub-steps during a time step
$\Delta t$, where the distribution function $f_{\alpha}(\bm{x},t+\Delta t)$ is
updated from $f_{\alpha}(\bm{x},t)$ for $\alpha=0,1,2,\ldots,8$ on a
rectangular lattice grid according to Eq. (66) to simulate the preconditioned
NS equations:

  * 1.

Compute pre-collision raw moments
We perform

$\mathbf{m}=\bm{\mathsf{P}}\mathbf{f}$

on the components of $\mathbf{f}(\bm{x},t)$, so that the pre-collision raw
moment components $k_{mn}^{\prime}$ of $\mathbf{m}$ for the usual square
lattice are obtained, where the mapping $\bm{\mathsf{P}}$ follows from Eq.
(50). In other words,

$\displaystyle k^{{}^{\prime}}_{\scriptscriptstyle 00}$ $\displaystyle=$
$\displaystyle f_{0}+f_{1}+f_{2}+f_{3}+f_{4}+f_{5}+f_{6}+f_{7}+f_{8},$
$\displaystyle k^{{}^{\prime}}_{\scriptscriptstyle 10}$ $\displaystyle=$
$\displaystyle f_{1}-f_{3}+f_{5}-f_{6}-f_{7}+f_{8},$ $\displaystyle
k^{{}^{\prime}}_{\scriptscriptstyle 01}$ $\displaystyle=$ $\displaystyle
f_{2}-f_{4}+f_{5}+f_{6}-f_{7}-f_{8},$ $\displaystyle
k^{{}^{\prime}}_{\scriptscriptstyle 20}$ $\displaystyle=$ $\displaystyle
f_{1}+f_{3}+f_{5}+f_{6}+f_{7}+f_{8},$ $\displaystyle
k^{{}^{\prime}}_{\scriptscriptstyle 02}$ $\displaystyle=$ $\displaystyle
f_{2}+f_{4}+f_{5}+f_{6}+f_{7}+f_{8},$ $\displaystyle
k^{{}^{\prime}}_{\scriptscriptstyle 11}$ $\displaystyle=$ $\displaystyle
f_{5}-f_{6}+f_{7}-f_{8},$ $\displaystyle k^{{}^{\prime}}_{\scriptscriptstyle
21}$ $\displaystyle=$ $\displaystyle f_{5}+f_{6}-f_{7}-f_{8},$ $\displaystyle
k^{{}^{\prime}}_{\scriptscriptstyle 12}$ $\displaystyle=$ $\displaystyle
f_{5}-f_{6}-f_{7}+f_{8},$ $\displaystyle k^{{}^{\prime}}_{\scriptscriptstyle
22}$ $\displaystyle=$ $\displaystyle f_{5}+f_{6}+f_{7}+f_{8}.$

  * 2.

Scale pre-collision raw moments for rectangular lattice
Using Eq. (53) for the scaling matrix $\bm{\mathsf{S}}$, implement

$\mathbf{m}\leftarrow\bm{\mathsf{S}}\mathbf{m},$

which yields the pre-collision raw moments for the _rectangular_ lattice
$k_{mn}^{\prime}$ from those computed for the square lattice above. Thus,

$\displaystyle k^{{}^{\prime}}_{\scriptscriptstyle
00}=k^{{}^{\prime}}_{\scriptscriptstyle 00},\quad
k^{{}^{\prime}}_{\scriptscriptstyle 10}=k^{{}^{\prime}}_{\scriptscriptstyle
10},\quad k^{{}^{\prime}}_{\scriptscriptstyle
01}=rk^{{}^{\prime}}_{\scriptscriptstyle 01},$ $\displaystyle
k^{{}^{\prime}}_{\scriptscriptstyle 20}=k^{{}^{\prime}}_{\scriptscriptstyle
20},\quad k^{{}^{\prime}}_{\scriptscriptstyle
02}=r^{2}k^{{}^{\prime}}_{\scriptscriptstyle 02},\quad
k^{{}^{\prime}}_{\scriptscriptstyle 11}=rk^{{}^{\prime}}_{\scriptscriptstyle
11},\quad k^{{}^{\prime}}_{\scriptscriptstyle
21}=rk^{{}^{\prime}}_{\scriptscriptstyle 21},\quad
k^{{}^{\prime}}_{\scriptscriptstyle
12}=r^{2}k^{{}^{\prime}}_{\scriptscriptstyle 12},\quad
k^{{}^{\prime}}_{\scriptscriptstyle
22}=r^{2}k^{{}^{\prime}}_{\scriptscriptstyle 22},$

which involve a scale factor of $r^{n}$ for the moment $k_{mn}^{\prime}$.

  * 3.

Compute pre-collision central moments
Transform the pre-collision raw moments $k^{{}^{\prime}}_{\scriptscriptstyle
mn}$ in $\mathbf{m}$ into pre-collision central moments $k_{mn}$ in
$\mathbf{m}^{c}$ via

$\mathbf{m}^{c}=\bm{\mathsf{F}}\mathbf{m},$

where $\bm{\mathsf{F}}$ is provided in Eq. (64). Thus, we get

$\displaystyle k_{00}$ $\displaystyle=$ $\displaystyle k_{00}^{\prime},$
$\displaystyle k_{10}$ $\displaystyle=$ $\displaystyle
k_{10}^{\prime}-u_{x}k_{00}^{\prime},$ $\displaystyle k_{01}$ $\displaystyle=$
$\displaystyle k_{01}^{\prime}-u_{y}k_{00}^{\prime},$ $\displaystyle k_{20}$
$\displaystyle=$ $\displaystyle
k_{20}^{\prime}-2u_{x}k_{10}^{\prime}+u_{x}^{2}k_{00}^{\prime},$
$\displaystyle k_{02}$ $\displaystyle=$ $\displaystyle
k_{02}^{\prime}-2u_{y}k_{01}^{\prime}+u_{y}^{2}k_{00}^{\prime},$
$\displaystyle k_{11}$ $\displaystyle=$ $\displaystyle
k_{11}^{\prime}-u_{y}k_{10}^{\prime}-u_{x}k_{01}^{\prime}+u_{x}u_{y}k_{00}^{\prime},$
$\displaystyle k_{21}$ $\displaystyle=$ $\displaystyle
k_{21}^{\prime}-2u_{x}k_{11}^{\prime}+u_{x}^{2}k_{01}^{\prime}-u_{y}k_{20}^{\prime}+2u_{x}u_{y}k_{10}^{\prime}-u_{x}^{2}u_{y}k_{00}^{\prime},$
$\displaystyle k_{12}$ $\displaystyle=$ $\displaystyle
k_{12}^{\prime}-2u_{y}k_{11}^{\prime}+u_{y}^{2}k_{10}^{\prime}-u_{x}k_{02}^{\prime}+2u_{x}u_{y}k_{01}^{\prime}-u_{x}u_{y}^{2}k_{00}^{\prime},$
$\displaystyle k_{22}$ $\displaystyle=$ $\displaystyle
k_{22}^{\prime}-2u_{x}k_{12}^{\prime}+u_{x}^{2}k_{02}^{\prime}-2u_{y}k_{21}^{\prime}+4u_{x}u_{y}k_{11}^{\prime}-2u_{x}^{2}u_{y}k_{01}^{\prime}+u_{y}^{2}k_{20}^{\prime}-2u_{x}u_{y}^{2}k_{10}^{\prime}+u_{x}^{2}u_{y}^{2}k_{00}^{\prime}.$

  * 4.

Compute post-collision central moments: Relaxation under collision using
preconditioned extended equilibria and source terms, and corrections
Perform the relaxations of central moments to their preconditioned equilibria,
including the source terms for the body force, using the different relaxation
rates given in $\bm{\mathsf{\Lambda}}$, and with corrections to eliminate the
grid-anisotropy and non-GI truncation errors parameterized by the grid aspect
ratio $r$ and the preconditioning parameter $\gamma$ reflecting the sub-step

$\tilde{\mathbf{m}}^{c}=\mathbf{m}^{c}+\bm{\mathsf{B}}^{-1}\left\\{\bm{\mathsf{\Lambda}}\;\left(\;\bm{\mathsf{B}}\mathbf{m}^{c,eq}-\bm{\mathsf{B}}\mathbf{m}^{c}\;\right)+\left(\bm{\mathsf{I}}-\frac{\bm{\mathsf{\Lambda}}}{2}\right)\bm{\mathsf{B}}\mathbf{\Phi}^{c}\Delta
t\right\\},$

where combining certain moments for their independent evolutions under
collision and their subsequent segregation following collision are formally
shown via the operators $\bm{\mathsf{B}}$ and $\bm{\mathsf{B}}^{-1}$,
respectively.

Thus, for the D2Q9 lattice, applying the operator $\bm{\mathsf{B}}$ implies
combining the second order diagonal components of moments as

$k_{2s}=\left(k_{20}+k_{02}\right),\quad k_{2d}=\left(k_{20}-k_{02}\right),$

which will relax to their equilibria at their own relaxation rates in what
follows rather than with using individual components in this regard. From Sec.
(2), the preconditioned central moment equilibria for the rectangular D2Q9
lattice, including the necessary correction terms, can be written as

$\displaystyle k_{00}^{eq}$ $\displaystyle=$
$\displaystyle\rho,\quad\quad\quad\quad k_{10}^{eq}=0,\quad\quad\quad\quad
k_{01}^{eq}=0,$ $\displaystyle k_{2s}^{eq}$ $\displaystyle=$ $\displaystyle
k_{20}^{eq}+k_{02}^{eq}=2q^{2}c_{s*}^{2}\rho+\left(\frac{1}{\gamma}-1\right)\rho\left(u_{x}^{2}+u_{y}^{2}\right)+\left({\theta}_{bx}\partial_{x}u_{x}+{\theta}_{by}\partial_{y}u_{y}+{\lambda}_{bx}\partial_{x}\rho+{\lambda}_{by}\partial_{y}\rho\right)\Delta
t,$ $\displaystyle k_{2d}^{eq}$ $\displaystyle=$ $\displaystyle
k_{20}^{eq}-k_{02}^{eq}=\;\;\qquad\quad\quad\left(\frac{1}{\gamma}-1\right)\rho\left(u_{x}^{2}-u_{y}^{2}\right)+\left(\theta_{sx}\partial_{x}u_{x}-\theta_{sy}\partial_{y}u_{y}+\lambda_{sx}\partial_{x}\rho+\lambda_{sy}\partial_{y}\rho\right)\Delta
t,$ $\displaystyle k_{11}^{eq}$ $\displaystyle=$
$\displaystyle\left(\frac{1}{\gamma}-1\right)\rho
u_{x}u_{y}+\left(\tilde{{\theta}}_{sx}\partial_{x}u_{x}+\tilde{\theta}_{sy}\partial_{y}u_{y}+\tilde{\lambda}_{sx}\partial_{x}\rho+\tilde{\lambda}_{sy}\partial_{y}\rho\right)\Delta
t,$ $\displaystyle k_{21}^{eq}$ $\displaystyle=$
$\displaystyle\left(\dfrac{1}{\gamma^{2}}-\dfrac{3}{\gamma}+2\right)\rho
u_{x}^{2}u_{y},\qquad
k_{12}^{eq}=\left(\dfrac{1}{\gamma^{2}}-\dfrac{3}{\gamma}+2\right)\rho
u_{x}u_{y}^{2},\qquad k_{22}^{eq}=q^{4}c_{s*}^{4}\rho.$

Here, and in the following, for better clarity, we use subscript ‘b’ for the
coefficients associated with the corrections for the trace of the diagonal
components of the second order moments (for bulk viscosity $\xi$), and the
subscript ‘s’ for the corrections for those other moments related to the shear
viscosity $\nu$, rather than the numerical subscripts used in Sec. 2. Then,
rewriting Eqs. (38), (39), and (40), respectively, the coefficients appearing
the above extended moment equilibria related to the corrections can be written
as

$\displaystyle\theta_{bx}=-\Big{[}P_{\gamma}u_{x}^{2}+Q_{\gamma}u_{y}^{2}+\left(3q^{2}c_{s*}^{2}-1\right)\Big{]}\rho\left(\frac{1}{\omega_{\xi}}-\frac{1}{2}\right),$
$\displaystyle\theta_{by}=-\Big{[}P_{\gamma}u_{y}^{2}+Q_{\gamma}u_{x}^{2}+\left(3q^{2}c_{s*}^{2}-r^{2}\right)\Big{]}\rho\left(\frac{1}{\omega_{\xi}}-\frac{1}{2}\right),$
$\displaystyle\lambda_{bx}=-\left[U_{\gamma}q^{2}c_{s*}^{2}-1\right]\left(\frac{1}{\omega_{\xi}}-\frac{1}{2}\right)u_{x},$
$\displaystyle\lambda_{3y}=-\left[U_{\gamma}q^{2}c_{s*}^{2}-r^{2}\right]\left(\frac{1}{\omega_{\xi}}-\frac{1}{2}\right)u_{y},$
$\displaystyle\theta_{sx}=-\Big{[}P_{\gamma}u_{x}^{2}-Q_{\gamma}u_{y}^{2}+\left(3q^{2}c_{s*}^{2}-1\right)\Big{]}\rho\left(\frac{1}{\omega_{\nu}}-\frac{1}{2}\right),$
$\displaystyle\theta_{sy}=+\Big{[}-P_{\gamma}u_{y}^{2}+Q_{\gamma}u_{x}^{2}-\left(3q^{2}c_{s*}^{2}-r^{2}\right)\Big{]}\rho\left(\frac{1}{\omega_{\nu}}-\frac{1}{2}\right),$
$\displaystyle\lambda_{sx}=-\left[U_{\gamma}q^{2}c_{s*}^{2}-1\right]\left(\frac{1}{\omega_{\nu}}-\frac{1}{2}\right)u_{x},$
$\displaystyle\lambda_{3y}=+\left[U_{\gamma}q^{2}c_{s*}^{2}-r^{2}\right]\left(\frac{1}{\omega_{\nu}}-\frac{1}{2}\right)u_{y},$

and

$\displaystyle\tilde{\theta}_{sx}=-\left(\frac{1}{\gamma^{2}}-\frac{1}{\gamma}\right)\left(\frac{1}{\omega_{\nu}}-\frac{1}{2}\right)\rho
u_{x}u_{y},$
$\displaystyle\tilde{\theta}_{sy}=-\left(\frac{1}{\gamma^{2}}-\frac{1}{\gamma}\right)\left(\frac{1}{\omega_{\nu}}-\frac{1}{2}\right)\rho
u_{x}u_{y},$
$\displaystyle\tilde{\lambda}_{sx}=-\left(\frac{1}{\gamma}-1\right)\left(\frac{1}{\omega_{\nu}}-\frac{1}{2}\right)q^{2}c_{s*}^{2}u_{y},$
$\displaystyle\tilde{\lambda}_{sy}=-\left(\frac{1}{\gamma}-1\right)\left(\frac{1}{\omega_{\nu}}-\frac{1}{2}\right)q^{2}c_{s*}^{2}u_{x},$

where

$P_{\gamma}=\frac{4}{\gamma^{2}}-\frac{1}{\gamma},\qquad
Q_{\gamma}=\frac{1}{\gamma^{2}}-\frac{1}{\gamma},\qquad
U_{\gamma}=\frac{2}{\gamma}+1,$

and the relaxation parameters $\omega_{\xi}$ and $\omega_{\nu}$ determine the
bulk and shear viscosities, respectively, which are shown at the end of this
sub-step. The spatial derivatives of the density $\partial_{x}\rho$ and
$\partial_{y}\rho$ appearing in the above moment correction terms are obtained
from an isotropic second order finite difference scheme, while the spatial
derivatives of the velocity field $\partial_{x}u_{x}$ and $\partial_{y}u_{y}$
in such terms are computed locally using non-equilibrium moments, by rewriting
Eqs. (44), (45), (46), and (48), as follows. First, writing

$\displaystyle
A=\cfrac{1}{2}\left(U_{\gamma}q^{2}c_{s*}^{2}-1\right)u_{x},\qquad
B=\cfrac{1}{2}\left(U_{\gamma}q^{2}c_{s*}^{2}-r^{2}\right)u_{y}$
$\displaystyle e_{b\rho}=-A\partial_{x}\rho-B\partial_{y}\rho,\qquad
e_{s\rho}=-A\partial_{x}\rho+B\partial_{y}\rho,$

and subsequently defining the following intermediate quantities

$R_{2s}=\left({k_{20}}+{k_{02}}\right)-2q^{2}c_{s*}^{2}\rho+e_{b\rho},\qquad
R_{2d}=\left({k_{20}}-{k_{02}}\right)-e_{s\rho},$

and

$\displaystyle C_{bx}$ $\displaystyle=$
$\displaystyle\left[-\frac{2q^{2}c_{s*}^{2}}{\omega_{\xi}}+\frac{P_{\gamma}u_{x}^{2}+Q_{\gamma}u_{y}^{2}}{2}+\frac{(3q^{2}c_{s*}^{2}-1)}{2}\right]\rho,\quad
C_{by}=\left[-\frac{2q^{2}c_{s*}^{2}}{\omega_{\xi}}+\frac{P_{\gamma}u_{y}^{2}+Q_{\gamma}u_{x}^{2}}{2}+\frac{(3q^{2}c_{s*}^{2}-r^{2})}{2}\right]\rho,$
$\displaystyle C_{sx}$ $\displaystyle=$
$\displaystyle\left[-\frac{2q^{2}c_{s*}^{2}}{\omega_{\nu}}+\frac{P_{\gamma}u_{x}^{2}-Q_{\gamma}u_{y}^{2}}{2}+\frac{(3q^{2}c_{s*}^{2}-1)}{2}\right]\rho,\quad
C_{sy}=\left[\frac{2q^{2}c_{s*}^{2}}{\omega_{\nu}}+\frac{-P_{\gamma}u_{y}^{2}+Q_{\gamma}u_{x}^{2}}{2}-\frac{(3q^{2}c_{s*}^{2}-r^{2})}{2}\right]\rho,$

the required local expressions for the derivatives of the velocity field are
then given by

$\partial_{x}u_{x}=\cfrac{\left(C_{sy}R_{2s}-C_{by}R_{2d}\right)}{\left(C_{bx}C_{sy}-C_{sx}C_{by}\right)},\qquad\partial_{y}u_{y}=\cfrac{1}{C_{by}}\left(R_{2s}-C_{bx}\partial_{x}u_{x}\right).$

With the above specifications, the post-collision central moments resulting
from the relaxations of the various central moments to their central moment
equilibria under collision and augmented by the effect of the source terms can
now be written as follows:

$\displaystyle\tilde{k}_{00}$ $\displaystyle=$ $\displaystyle k_{00},$
$\displaystyle\tilde{k}_{10}$ $\displaystyle=$ $\displaystyle
k_{10}+\omega_{1}(k_{10}^{eq}-k_{10})+(1-\omega_{1}/2)\sigma_{10}\Delta t,$
$\displaystyle\tilde{k}_{01}$ $\displaystyle=$ $\displaystyle
k_{01}+\omega_{1}(k_{01}^{eq}-k_{01})+(1-\omega_{1}/2)\sigma_{01}\Delta t,$
$\displaystyle\tilde{k}_{2s}$ $\displaystyle=$ $\displaystyle
k_{2s}+\omega_{\xi}(k_{2s}^{eq}-k_{2s})+(1-\omega_{\xi}/2)\sigma_{2s}\Delta
t,$ $\displaystyle\tilde{k}_{2d}$ $\displaystyle=$ $\displaystyle
k_{2d}+\omega_{\nu}(k_{2d}^{eq}-k_{2d})+(1-\omega_{\nu}/2)\sigma_{2d}\Delta
t,$ $\displaystyle\tilde{k}_{11}$ $\displaystyle=$ $\displaystyle
k_{11}+\omega_{\nu}(k_{11}^{eq}-k_{11})+(1-\omega_{\nu}/2)\sigma_{11}\Delta
t,$ $\displaystyle\tilde{k}_{21}$ $\displaystyle=$ $\displaystyle
k_{21}+\omega_{21}(k_{21}^{eq}-k_{21}),$ $\displaystyle\tilde{k}_{12}$
$\displaystyle=$ $\displaystyle k_{12}+\omega_{12}(k_{12}^{eq}-k_{12}),$
$\displaystyle\tilde{k}_{22}$ $\displaystyle=$ $\displaystyle
k_{22}+\omega_{22}(k_{22}^{eq}-k_{22}),$

where $\sigma_{2s}=\sigma_{20}+\sigma_{02}$ and
$\sigma_{2d}=\sigma_{20}-\sigma_{02}$, and the central moments of the source
terms relevant for recovering the preconditioned NS equations with a body
force, $\sigma_{10}$, $\sigma_{01}$, $\sigma_{20}$, $\sigma_{02}$, and
$\sigma_{11}$ are given in Eq. (3.3). Here, the relaxation parameters
$\omega_{\nu}$ and $\omega_{\xi}$ are related to the shear viscosity $\nu$ and
bulk viscosity $\xi$, respectively, through the following expressions:

$\nu=\gamma\;q^{2}c_{s*}^{2}\left(\frac{1}{\omega_{\nu}}-\frac{1}{2}\right)\Delta
t,\qquad\xi=\gamma
q^{2}c_{s*}^{2}\left(\frac{1}{\omega_{\xi}}-\frac{1}{2}\right)\Delta t.$

Notice that these transport coefficients are functions of the preconditioning
parameter $\gamma$ and the grid aspect ratio $r$ (via $q$). In this work, we
set the relaxation parameter associated with the bulk viscosity as well as
those for the other, especially the higher order, moments to unity, i.e.,
$\omega_{\xi}=\omega_{1}=\omega_{21}=\omega_{12}=\omega_{22}=1.0$. To complete
this sub-step for collision, we now segregate the post-collision combined
central moments $\tilde{k}_{2s}$ and $\tilde{k}_{2d}$ into the bare central
moments $\tilde{k}_{20}$ and $\tilde{k}_{02}$ (reflecting the application of
the inverse operator $\bm{\mathsf{B}}^{-1}$) via

$\tilde{k}_{20}=\frac{1}{2}(\tilde{k}_{2s}+\tilde{k}_{2d}),\qquad\tilde{k}_{02}=\frac{1}{2}(\tilde{k}_{20}-\tilde{k}_{02}).$

As a result, all the post-collision bare central moments supported by the
rectangular D2Q9 lattice is now computed.

  * 5.

Compute post-collision raw moments
Evaluating

$\tilde{\mathbf{m}}=\bm{\mathsf{F}}^{-1}\tilde{\mathbf{m}}^{c},$

the post-collision raw moments $\tilde{k}_{mn}^{\prime}$ can be obtained from
the corresponding central moments $\tilde{k}_{mn}$ computed in the previous
sub-step, where $\bm{\mathsf{F}}^{-1}=\bm{\mathsf{F}}(-u_{x},-u_{y})$, with
(see Eq. (64) for $\bm{\mathsf{F}}$). Thus, we have

$\displaystyle\tilde{k}_{00}=\tilde{k}_{00}^{\prime},$
$\displaystyle\tilde{k}_{10}=\tilde{k}_{10}^{\prime}+u_{x}\tilde{k}_{00}^{\prime},$
$\displaystyle\tilde{k}_{01}=\tilde{k}_{01}^{\prime}+u_{y}\tilde{k}_{00}^{\prime},$
$\displaystyle\tilde{k}_{20}=\tilde{k}_{20}^{\prime}+2u_{x}\tilde{k}_{10}^{\prime}+u_{x}^{2}\tilde{k}_{00}^{\prime},$
$\displaystyle\tilde{k}_{02}=\tilde{k}_{02}^{\prime}+2u_{y}\tilde{k}_{01}^{\prime}+u_{y}^{2}\tilde{k}_{00}^{\prime},$
$\displaystyle\tilde{k}_{11}=\tilde{k}_{11}^{\prime}+u_{y}\tilde{k}_{10}^{\prime}+u_{x}\tilde{k}_{01}^{\prime}+u_{x}u_{y}\tilde{k}_{00}^{\prime},$
$\displaystyle\tilde{k}_{21}=\tilde{k}_{21}^{\prime}+2u_{x}\tilde{k}_{11}^{\prime}+u_{x}^{2}\tilde{k}_{01}^{\prime}+u_{y}\tilde{k}_{20}^{\prime}+2u_{x}u_{y}\tilde{k}_{10}^{\prime}+u_{x}^{2}u_{y}\tilde{k}_{00}^{\prime},$
$\displaystyle\tilde{k}_{12}=\tilde{k}_{12}^{\prime}+2u_{y}\tilde{k}_{11}^{\prime}+u_{y}^{2}\tilde{k}_{10}^{\prime}+u_{x}\tilde{k}_{02}^{\prime}+2u_{x}u_{y}\tilde{k}_{01}^{\prime}+u_{x}u_{y}^{2}\tilde{k}_{00}^{\prime},$
$\displaystyle\tilde{k}_{22}=\tilde{k}_{22}^{\prime}+2u_{x}\tilde{k}_{12}^{\prime}+u_{x}^{2}\tilde{k}_{02}^{\prime}+2u_{y}\tilde{k}_{21}^{\prime}+4u_{x}u_{y}\tilde{k}_{11}^{\prime}+2u_{x}^{2}u_{y}\tilde{k}_{01}^{\prime}+u_{y}^{2}\tilde{k}_{20}^{\prime}+2u_{x}u_{y}^{2}\tilde{k}_{10}^{\prime}+u_{x}^{2}u_{y}^{2}\tilde{k}_{00}^{\prime}.$

  * 6.

Apply inverse scaling of post-collision raw moments for rectangular lattice
Perform

$\tilde{\mathbf{m}}\leftarrow\bm{\mathsf{S}}^{-1}\tilde{\mathbf{m}}$

using Eq. (55) for the inverse scaling $\bm{\mathsf{S}}^{-1}$, which involves
applying an inverse scale factor $r^{-n}$ for the raw moment
$\tilde{k}_{mn}^{\prime}$ computed in the previous sub-step so that we have

$\displaystyle\tilde{k}_{00}^{\prime}=\tilde{k}_{00}^{\prime},\quad\tilde{k}_{10}^{\prime}=\tilde{k}_{10}^{\prime},\quad\tilde{k}_{01}^{\prime}=\frac{1}{r}\tilde{k}_{01}^{\prime},$
$\displaystyle\tilde{k}_{20}^{\prime}=\tilde{k}_{20}^{\prime},\quad\tilde{k}_{02}^{\prime}=\frac{1}{r^{2}}\tilde{k}_{02}^{\prime},\quad\tilde{k}_{11}^{\prime}=\frac{1}{r}\tilde{k}_{11}^{\prime},\quad\tilde{k}_{21}^{\prime}=\frac{1}{r}\tilde{k}_{21}^{\prime},\quad\tilde{k}_{12}^{\prime}=\frac{1}{r^{2}}\tilde{k}_{12}^{\prime},\quad\tilde{k}_{22}^{\prime}=\frac{1}{r^{2}}\tilde{k}_{22}^{\prime}.$

This enables a more efficient transformation of the raw moments to
distribution functions involving the inverse of simpler moment basis
$\bm{\mathsf{P}}$ (see Eq. 50) for the square lattice, i.e.,
$\bm{\mathsf{P}}^{-1}$ in the next sub-step.

  * 7.

Compute post-collision distribution functions
Invoking the inverse mapping

$\tilde{\mathbf{f}}=\bm{\mathsf{P}}^{-1}\tilde{\mathbf{m}},$

we then obtain the post-collision distribution functions on the rectangular
D2Q9 lattice as follows:

$\displaystyle\tilde{f}_{0}$ $\displaystyle=$
$\displaystyle\left(\tilde{k}_{00}^{\prime}-\tilde{k}_{20}^{\prime}-\tilde{k}_{02}^{\prime}+\tilde{k}_{22}^{\prime}\right),$
$\displaystyle\tilde{f}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\tilde{k}_{10}^{\prime}+\tilde{k}_{20}^{\prime}-\tilde{k}_{12}^{\prime}-\tilde{k}_{22}^{\prime}\right),$
$\displaystyle\tilde{f}_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\tilde{k}_{01}^{\prime}+\tilde{k}_{02}^{\prime}-\tilde{k}_{21}^{\prime}-\tilde{k}_{22}^{\prime}\right),$
$\displaystyle\tilde{f}_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(-\tilde{k}_{10}^{\prime}+\tilde{k}_{20}^{\prime}+\tilde{k}_{12}^{\prime}-\tilde{k}_{22}^{\prime}\right),$
$\displaystyle\tilde{f}_{4}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(-\tilde{k}_{01}^{\prime}+\tilde{k}_{02}^{\prime}+\tilde{k}_{21}^{\prime}-\tilde{k}_{22}^{\prime}\right),$
$\displaystyle\tilde{f}_{5}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(\tilde{k}_{11}^{\prime}+\tilde{k}_{21}^{\prime}+\tilde{k}_{12}^{\prime}+\tilde{k}_{22}^{\prime}\right),$
$\displaystyle\tilde{f}_{6}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(-\tilde{k}_{11}^{\prime}+\tilde{k}_{21}^{\prime}-\tilde{k}_{12}^{\prime}+\tilde{k}_{22}^{\prime}\right),$
$\displaystyle\tilde{f}_{7}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(\tilde{k}_{11}^{\prime}-\tilde{k}_{21}^{\prime}-\tilde{k}_{12}^{\prime}+\tilde{k}_{22}^{\prime}\right),$
$\displaystyle\tilde{f}_{8}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(-\tilde{k}_{11}^{\prime}-\tilde{k}_{21}^{\prime}+\tilde{k}_{12}^{\prime}+\tilde{k}_{22}^{\prime}\right).$

  * 8.

Stream distribution functions along particle characteristics
Performing perfect shift advection on the neighboring lattice nodes via

$f_{\alpha}(\bm{x},t+\Delta
t)=\widetilde{f}_{\alpha}(\bm{x}-\bm{e}_{\alpha}\Delta t,t),$

we obtain the updated distribution functions $f_{\alpha}$ at the end of time
step $t+\Delta t$.

  * 9.

Update hydrodynamic fields
From the distribution functions obtained in the previous sub-step, update the
fluid density and velocities, as well as the pressure field via

$\rho=\sum_{\alpha=0}^{8}f_{\alpha},\quad\rho\bm{u}=\sum_{\alpha=0}^{8}f_{\alpha}\bm{e}_{\alpha}+\frac{1}{2\gamma}\bm{F}\Delta
t,\quad p=\gamma q^{2}c_{s*}^{2}\rho.$

It may be noted that for implementing moving wall no-slip boundary conditions
for simulating shear flows on rectangular grids, the extension of the
momentum-augmented half-way bounce back scheme with parametrization by the
grid aspect ratio $r$ presented in our previous work [34] can be used in the
present LB formulation by noting that the speed of sound $c_{s}$ should be
specified consistently using $c_{s}=qc_{s*}$, where $q=\mbox{min}(r,1)$.

## References

  * [1] G. R. McNamara, G. Zanetti, Use of the Boltzmann equation to simulate lattice-gas automata, Physical review letters 61 (20) (1988) 2332.
  * [2] R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation: theory and applications, Physics Reports 222 (3) (1992) 145–197.
  * [3] P. Lallemand, L.-S. Luo, M. Krafczyk, W.-A. Yong, The lattice Boltzmann method for nearly incompressible flows, Journal of Computational Physics 431 (2021) 109713.
  * [4] Y.-H. Qian, D. d’Humières, P. Lallemand, Lattice BGK models for navier-stokes equation, EPL (Europhysics Letters) 17 (6) (1992) 479.
  * [5] D. d’Humieres, Multiple–relaxation–time lattice Boltzmann models in three dimensions, Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 360 (1792) (2002) 437–451.
  * [6] M. Geier, A. Greiner, J. G. Korvink, Cascaded digital lattice Boltzmann automata for high Reynolds number flow, Physical Review E 73 (6) (2006) 066705\.
  * [7] M. Geier, M. Schönherr, A. Pasquali, M. Krafczyk, The cumulant lattice Boltzmann equation in three dimensions: Theory and validation, Computers & Mathematics with Applications 70 (4) (2015) 507–547.
  * [8] I. V. Karlin, A. Ferrante, H. C. Öttinger, Perfect entropy functions of the lattice Boltzmann method, Europhys. Lett. 47 (2) (1999) 182.
  * [9] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. M. Viggen, The lattice Boltzmann method, Springer International Publishing 10 (2017) 978–3.
  * [10] P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, galilean invariance, and stability, Physical Review E 61 (6) (2000) 6546.
  * [11] J. Koelman, A simple lattice Boltzmann scheme for navier-stokes fluid flow, EPL (Europhysics Letters) 15 (6) (1991) 603.
  * [12] L. A. Hegele Jr, K. Mattila, P. C. Philippi, Rectangular lattice-boltzmann schemes with bgk-collision operator, Journal of Scientific Computing 56 (2) (2013) 230–242.
  * [13] C. Peng, Z. Guo, L.-P. Wang, A lattice-BGK model for the Navier-Stokes equations based on a rectangular grid, Computers & Mathematics with Applications (2016).
  * [14] Z. Wang, J. Zhang, Simulating anisotropic flows with isotropic lattice models via coordinate and velocity transformation, International Journal of Modern Physics C 30 (10) (2019) 1941001.
  * [15] M. Bouzidi, D. d’Humières, P. Lallemand, L.-S. Luo, Lattice Boltzmann equation on a two-dimensional rectangular grid, Journal of Computational Physics 172 (2) (2001) 704–717.
  * [16] J. G. Zhou, MRT rectangular lattice Boltzmann method, International Journal of Modern Physics C 23 (05) (2012) 1250040.
  * [17] C. Peng, H. Min, Z. Guo, L.-P. Wang, A hydrodynamically-consistent MRT lattice Boltzmann model on a 2d rectangular grid, Journal of Computational Physics 326 (2016) 893–912.
  * [18] P. Asinari, Generalized local equilibrium in the cascaded lattice Boltzmann method, Physical Review E 78 (1) (2008) 016701.
  * [19] K. N. Premnath, S. Banerjee, Incorporating forcing terms in cascaded lattice boltzmann approach by method of central moments, Physical Review E 80 (3) (2009) 036702.