<EMAIL_ADDRESS>(K.N.Premnath)

# Inertial Frame Independent Forcing for Discrete Velocity Boltzmann Equation:
Implications for Filtered Turbulence Simulation

Kannan N. Premnath  1  and Sanjoy Banerjee  2  11affiliationmark:  Department
of Mechanical Engineering, University of Wyoming, Laramie, WY 82071, U.S.A.
22affiliationmark:  Department of Chemical Engineering, City College of New
York, City University of New York, New York, NY 10031, U.S.A.

###### Abstract

We present a systematic derivation of a model based on the central moment
lattice Boltzmann equation that rigorously maintains Galilean invariance of
forces to simulate inertial frame independent flow fields. In this regard, the
central moments, i.e. moments shifted by the local fluid velocity, of the
discrete source terms of the lattice Boltzmann equation are obtained by
matching those of the continuous full Boltzmann equation of various orders.
This results in an exact hierarchical identity between the central moments of
the source terms of a given order and the components of the central moments of
the distribution functions and sources of lower orders. The corresponding
source terms in velocity space are then obtained from an exact inverse
transformation due to a suitable choice of orthogonal basis for moments.
Furthermore, such a central moment based kinetic model is further extended by
incorporating reduced compressibility effects to represent incompressible
flow. Moreover, the description and simulation of fluid turbulence for full or
any subset of scales or their averaged behavior should remain independent of
any inertial frame of reference. Thus, based on the above formulation, a new
approach in lattice Boltzmann framework to incorporate turbulence models for
simulation of Galilean invariant statistical averaged or filtered turbulent
fluid motion is discussed.

###### keywords:

Lattice Boltzmann Method, Central Moments, Galilean Invariance, Turbulence,
Filtering

05.20.Dd,47.27.-i,47.27.E-

##  1 Introduction

Minimal kinetic models for the Boltzmann equation, i.e. lattice Boltzmann
equation formulations, are evolving towards as alternative physically-inspired
computational techniques for various fluid mechanics and other problems.
Originally developed as an improved variant of the lattice gas automata [1] to
eliminate statistical noise [2], the lattice Boltzmann method (LBM) has
undergone a series of major refinements, in terms of its underlying physical
models as well as numerical solution schemes for various applications over the
last two decades [3, 4, 5, 6]. In particular, its rigorous connection to the
kinetic theory [7, 8, 9] has resulted in a number of recent developments,
including models that are more physically consistent for multiphase [10, 11]
and multicomponent flows [12], models for non-equilibrium phenomena beyond the
Navier-Stokes-Fourier representation [13] and an asymptotic analysis approach
to establish consistency of the LBM from a numerical point of view [14].

The stream-and-collide procedure of the LBM can be considered as a
dramatically simplified discrete representation of the continuous Boltzmann
equation. Here, the streaming step represents the advection of the
distribution of particle populations along discrete directions, which are
designed from symmetry considerations, between successive collisions. Much of
the physical effects being modeled are represented in terms of the collision
step, which also significantly influences the numerical stability of the LBM.
Most of the major developments until recently were associated with the single-
relaxation-time (SRT) model [15, 16] based on the BGK approximation [17], and
enjoys its popularity owing, mainly, to its simplicity. However, it is prone
to numerical instability. Moreover, it is inadequate in its representation of
certain physical aspects, such as independently adjustable transport
properties of thermal transport and viscoelastic phenomena.

These limitations have been significantly addressed in the multiple-
relaxation-time (MRT) collision model [18]. This, in a sense, represents a
simplified form of the relaxation LBM proposed earlier [19, 20], with an
important characteristic difference in that the collision process is carried
out in moment space [21] instead of in the usual velocity space. By separating
the relaxation time scales of different moments, obtained by using a linear
Fourier stability analysis, its numerical stability can be significantly
improved [22, 23]. Furthermore, it has resulted in significant advantages over
the SRT-LBM for computation of various classes of fluid flow problems,
including multiphase systems [24, 25, 26], turbulent flows [27, 28] and
magnetohydrodynamics [29]. It may be noted that recently a different form of
MRT model based on the orthogonal Hermite polynomial projections of the
distribution functions, which is independent of any underlying lattice
structure, allowing representation of higher order non-equilibrium effects has
been proposed [30].

The stabilization of the LBM using a single relaxation time has been addressed
from a different perspective by enforcing the H-theorem locally in the
collision step [31, 32, 33, 34]. By using the attractors of the distribution
function based on the minimization of a Lyapunov-type functional, non-linear
stability of the LBM is achieved in this Entropic LBM. This approach has
recently been significantly extended to incorporate multiple relaxation times
with efficient implementation strategies [35, 36]. Furthermore, systematic
procedures for different types of higher-order LBM have been developed [37,
38, 39]. An important element is the construction of higher-order lattices
based on symmetry considerations which have been analyzed using group theory
[40, 41]. Further progress, from a numerical aspect, is that based on the
consistency analysis [14] and a notion of structural stability [42, 43] (shown
related to the Onsager-like relation in non-equilibrium thermodynamics [44]),
convergence of the LBM to the Navier-Stokes equations has rigorously been
shown [45].

On the other hand, it is important to clearly understand in what sense the
lattice Boltzmann equation (LBE), which is generally considered as a
mesoscopic approach, inherits or maintains the various physical invariance
properties of the continuous full Boltzmann equation (which it represents as a
much simplified model) and the Navier-Stokes equations (which it represents
numerically). Careful considerations of these aspects play an important role
in ensuring the general applicability of the approach for various, especially
challenging, problems. In this regard, and to put the present work in
perspective, it should be noted that the continuum mechanics description as
well as the microscopic statistical (continuous Boltzmann) description of
fluid motion generally satisfy a larger invariance group, with inertial frame
invariance being just an important special case. The most general form among
these is the so-called the principal of material frame indifference, also
known as the objectivity principle [46]. According to this, the constitutive
equations should have the same forms in _all_ frames of reference, whether
inertial or not. While this is considered as an important axiom based on which
the continuum mechanics is formulated [46], its role from continuous kinetic
theory point of view was the subject of considerable analysis for sometime
[47, 48, 49, 50, 51, 52].

The following are the main outcomes of these studies: the continuous full
Boltzmann equation (i) is material frame indifferent in a _strong approximate
sense_ , when there is a large scale separation between the collision times
and the macroscopic flow times [49, 51, 52] (thus providing a strong support
to the axiomatic principle generally used in the continuum description), (ii)
satisfies both the inertial frame or Galilean invariance as well as the
extended Galilean invariance (i.e. invariance under arbitrary translational
accelerations of the reference frame) _exactly_ [52]. Furthermore, it was
shown that the standard procedures (e.g. Chapman-Enskog expansion, Maxwellian
iteration) lead to frame dependent higher order contributions for the
constitutive equations in non-inertial frames in the continuum limit, while
the continuous kinetic theory itself can be frame independent [49]. Careful
considerations of these principles could guide in the development of more
generally applicable models and numerical schemes for complex problems. For
example, material frame indifference (point (i)) is generally used as an
important constraint for the constitutive equations for complex fluids (e.g.
beyond Newtonian constitutive description such as polymeric fluids) and in the
development of turbulence models in continuum mechanics. As mentioned above,
this property is satisfied in a strong approximate sense by the continuous
Boltzmann equation, but not necessarily by the tools that relate the
microscopic and macroscopic descriptions(point iii). This aspects are
pertinent in the construction of complex models from the continuous Boltzmann
equation (e.g. [53, 54]). In this work, however, we limit our discussion to an
association of the properties mentioned in a part of the point (ii) for the
LBE, i.e. for the exact invariance group – the Galilean or inertial frame
invariance.

In this regard, as a discrete approximation to the continuous full Boltzmann
equation, the development of the LBE consists of simplifications at different
levels. Thus, its various elements should be analyzed carefully to ascertain
and quantify as to how well it satisfies Galilean invariance. First, in
contrast to continuous kinetic theory, due to the choice of finite lattice
velocity sets and associated symmetries, it introduces linear dependencies of
higher order moments with those of lower order moments that are supported by
the lattice set [40]. Such degeneracies can in turn lead to negative
dependence of viscosity on fluid velocity. It generally causes the Galilean
invariance to be broken by the presence of terms that are cubic in velocity
for the standard lattice configurations (with symmetries of square in 2D and
cube in 3D) and also leads to numerical instability, especially at higher Mach
numbers. This issue can be alleviated by the use of extended lattice velocity
sets, which then relegates the degeneracies among moments to even higher
orders. Second, the collision step including the forcing terms of the LBE
should be carefully constructed in such a way that they recover correct
physics which is inertial frame independent, i.e. the Navier-Stokes equations.
Here, the use of independent set of _central_ moments for a chosen lattice
provides a natural approach to maintain Galilean invariance that can be
constructed by invoking elements directly from kinetic theory. This is the
main goal of the present work (see below). A rational means to more
efficiently account for both the above aspects is discussed in the last
section of this paper. And, third, the streaming step of the LBE is generally
constructed as a discrete Lagrangian process. In the standard implementation,
this couples the particle velocity and configurations spaces, which in turn,
constrains the numerical accuracy of the LBE in the representation of the
Navier-Stokes equations. As a result, the Galilean invariance of the LBE is
limited by its overall numerical accuracy. However, it is known that such
coupling between physical and lattice symmetries is not necessary in the
discretization of the streaming operator. In fact, it can be discretized using
classical schemes such as finite-difference or finite-element methods that
alleviate this issue [55, 56, 57]. Specifically, exploiting higher order
discretization and time integration schemes (e.g. [58]) for the streaming
operator could further improve the order of accuracy (and hence the Galilean
invariance) of the LBE. Furthermore, the use of implicit schemes could enhance
the computational efficiency in this regard.

Focusing on the second aspect mentioned above, a different type of collision
operator and forcing can be devised that can maintain Galilean invariance for
a chosen lattice velocity set and a discretization scheme for the streaming
step. Specifically, central moments are relaxed in a moving frame of reference
during collision step [59], originally proposed to improve numerical
stability, but emphasized here for its better physical coherence. The use of
central moments, which are obtained by shifting the particle velocity with the
local fluid velocity [60], rigorously enforces Galilean invariance. In
particular, while other previous approaches are generally Galilean invariant
for up to second-order moments, the central moment based approach provides a
higher order frame invariance as permitted by the discrete lattice velocity
set. This approach was examined based on the concept of generalized local
equilibrium [61]. In addition, to further improve physical coherence, the
attractors for the higher order central moments were constructed as products
of the lower order central moments, leading to the factorized central moment
method [62]. Recently, a new approach to incorporate source terms using
central moments in the LBM that are Galilean invariant by construction, which
are important for computation of various physical problems, was developed
[63]. The consistency of this technique to the Navier-Stokes equations was
shown by means of the Chapman-Enskog analysis [64] and its numerical accuracy
was established. Furthermore, the method was also extended in three-dimensions
for various lattice velocity sets and validated for a class of canonical
problems [65]. As clarified in [63, 65], numerical stability of the central
moment approach can be enhanced, when it is executed in a multiple relaxation
time formulation, similar to the standard or raw moment based approaches.
Interestingly, it has been shown recently that when some classical schemes for
flow simulation are made to satisfy Galilean invariance more rigorously, they
led to more robust implementations (e.g. [66, 67]).

Turbulence remains as among the most challenging classes of flows for which
considerable effort has been focused on the development of theory and
applications using the LBM. Since its roots can be traced to kinetic theory,
the LBM has been analyzed for the development of turbulence models from a
fundamental point of view [68, 69, 70]. It has been employed for computation
of Reynolds-averaged description of turbulent flows [71, 72]. Furthermore, it
has found applications for large eddy simulation (LES) using LBM formulations
with SRT [73], and MRT [74] with multiblock approach for efficient
implementation [27]. Recently, dynamic subgrid scale (SGS) models for LES were
incorporated into the LBM framework that resulted in reduced empiricism for
description at such scales [28]. Moreover, an improved inertial-range
consistent SGS model was also proposed [75]. A theoretical formulation for a
SGS model based on an approximate deconvolution method [78] that does not rely
on the common eddy-viscosity concept for application with the LBM was also
devised recently [76]. Lastly, the closure modeling issues of kinetic and
continuum turbulence effects were reconciled in a unified statistical/filtered
description using a modified kinetic equation [77]. Effectively, this allows
the use of macroscopic turbulence models involving divergence of the Reynolds
stress in the forcing term of the kinetic equation.

An important physical consideration for any description of turbulent flow is
that it should be invariant for all inertial frames of reference. In other
words, for general applicability, representation of turbulence for all or any
subset of its scales should be Galilean invariant. Thus, in particular, all
SGS models, and associated numerical schemes for turbulence computation,
should be frame invariant. An insightful analysis of various turbulence models
was carried out from this viewpoint in [79]. A method to achieve Galilean
invariance by means of certain redefinition of turbulent stresses was
discussed in [80]. A recent review on this subject is reported in [81].
Furthermore, it should be noted that concepts based on central moments have
played an important role in the development of theoretical foundation of
turbulence physics – such as for statistical turbulence theory [82] and
turbulence modeling [83].

Thus, in this paper, we develop a lattice Boltzmann equation based on central
moments for Galilean invariant representation of turbulent flows.
Specifically, it allows frame-independent incorporation of general models for
turbulent Reynolds stresses in a statistical/filter averaged formulation using
LBM for turbulence simulation. Furthermore, in a general setting, it maintains
the forces and stresses to be independent of any inertial frame of reference
and could also improve numerical stability in computations. In [63], we
developed a forcing scheme based on a particular ansatz involving the local
Maxwell distribution. Here, we develop a general forcing based on central
moments by a direct examination of the continuous full Boltzmann equation
itself, which unlike [63] could also self-consistently account for non-
equilibrium effects in higher order terms. In this regard, the central moments
of the resulting source terms of the continuous and discrete counterparts are
matched successively at different orders leading to a cascaded structure. In
essence, this approach can be considered as a Galilean invariant minimal
discrete model for the full Boltzmann equation including forcing terms. The
attractors for higher order central moments in the collision step of this
computational model is based on the factorization in terms of those at lower
orders by including such general forcing terms. In addition, we further
develop this approach with reduced compressibility effects for improved
representation of turbulent flow physics in the incompressible limit. The
forcing formulation developed here for incorporating turbulence models in a
statistical/filtered formulation can be extended to other problems, such as,
for example, Galilean invariant representation of forces or stresses in
complex fluids.

The paper is organized as follows. Section 2 briefly discusses the choice of
the moment basis employed in this paper and Sec. 3 the continuous Boltzmann
equation. In Secs. 4 and 5, continuous forms of the central moments for the
distribution functions and its local equilibrium, and sources due to force
fields, respectively, are introduced. The LBE based on central moments with
the general forcing terms is presented in Sec. 6. Various discrete central
moments are presented in Sec. 7 that also specifies a matching principle to
maintain Galilean invariance and the relationships among such moments are
provided in Sec. 8. Section 9 describes various discrete raw moments and the
derivation of the source terms in terms of the discrete particle velocity
space. In Sec. 10,we present the construction of the collision operator of the
central moment based LBM. The computational procedure of this approach is
provided in Sec. 11. The derivation is extended by considering reduced
compressibility effects in Sec. 12. Furthermore, Sec. 13 discusses the use of
attractors of the higher order central moments based on the concept of their
factorization in term of those at lower orders. A natural consequence of this
overall approach is that turbulence models can be represented for Galilean
invariant filtered turbulence simulation using the LBM, which is described in
Sec. 14. Finally, the summary and conclusions of this work are discussed in
Sec. 15.

##  2 Selection of Moment Basis

An important element in the development of the central moment based LBM is the
specification of a suitable basis for moments. In this work, to elucidate our
approach, the two-dimensional, nine velocity (D2Q9) lattice model (see Fig. 1)
is considered, for which the moment basis used in [63] is adopted. It should,
however, be noted that the procedure described henceforth is of general
nature, and can be extended for other lattice models and in three dimensions.
The particle velocity for this lattice model $\overrightarrow{e}_{\alpha}$ is
given by

$\overrightarrow{e_{\alpha}}=\left\\{\begin{array}[]{ll}{(0,0)}&{\alpha=0}\\\
{(\pm 1,0),(0,\pm 1)}&{\alpha=1,\cdots,4}\\\ {(\pm 1,\pm
1)}&{\alpha=5,\cdots,8}\end{array}\right.$ (1)

Figure 1:  Two-dimensional, nine-velocity (D2Q9) Lattice.

For convenience, we employ Dirac’s bra-ket notion to represent the basis
vectors, and Greek and Latin subscripts for particle velocity directions and
Cartesian coordinate directions, respectively. Noting that moments in the LBM
are discrete integral properties of the distribution function $f_{\alpha}$,
i.e. $\sum_{\alpha=0}^{8}e_{\alpha x}^{m}e_{\alpha y}^{n}f_{\alpha}$, where
$m$ and $n$ are integers, we begin with the following nine non-orthogonal
independent basis vectors obtained by combining monomials $e_{\alpha
x}^{m}e_{\alpha y}^{n}$ in an ascending order. That is,

$\displaystyle\ket{\rho}\equiv\ket{|\overrightarrow{e}_{\alpha}|^{0}}$
$\displaystyle=$ $\displaystyle\left(1,1,1,1,1,1,1,1,1\right)^{T},$
$\displaystyle\ket{e_{\alpha x}}$ $\displaystyle=$
$\displaystyle\left(0,1,0,-1,0,1,-1,-1,1\right)^{T},$
$\displaystyle\ket{e_{\alpha y}}$ $\displaystyle=$
$\displaystyle\left(0,0,1,0,-1,1,1,-1,-1\right)^{T},$
$\displaystyle\ket{e_{\alpha x}^{2}+e_{\alpha y}^{2}}$ $\displaystyle=$
$\displaystyle\left(0,1,1,1,1,2,2,2,2\right)^{T},$
$\displaystyle\ket{e_{\alpha x}^{2}-e_{\alpha y}^{2}}$ $\displaystyle=$
$\displaystyle\left(0,1,-1,1,-1,0,0,0,0\right)^{T},$ (2)
$\displaystyle\ket{e_{\alpha x}e_{\alpha y}}$ $\displaystyle=$
$\displaystyle\left(0,0,0,0,0,1,-1,1,-1\right)^{T},$
$\displaystyle\ket{e_{\alpha x}^{2}e_{\alpha y}}$ $\displaystyle=$
$\displaystyle\left(0,0,0,0,0,1,1,-1,-1\right)^{T},$
$\displaystyle\ket{e_{\alpha x}e_{\alpha y}^{2}}$ $\displaystyle=$
$\displaystyle\left(0,0,0,0,0,1,-1,-1,1\right)^{T},$
$\displaystyle\ket{e_{\alpha x}^{2}e_{\alpha y}^{2}}$ $\displaystyle=$
$\displaystyle\left(0,0,0,0,0,1,1,1,1\right)^{T},$

where the superscript ‘$T$’ represents the transpose operator.

For an efficient implementation, the above non-orthogonal basis set is
transformed into an equivalent orthogonal set through the Gram-Schmidt
procedure in the increasing order of the monomials of the products of the
Cartesian components of the particle velocities [63]:

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

This can be written explicitly in term of a matrix given by

$\mathcal{K}=\left[\begin{array}[]{rrrrrrrrr}1&0&0&-4&0&0&0&0&4\\\
1&1&0&-1&1&0&0&2&-2\\\ 1&0&1&-1&-1&0&2&0&-2\\\ 1&-1&0&-1&1&0&0&-2&-2\\\
1&0&-1&-1&-1&0&-2&0&-2\\\ 1&1&1&2&0&1&-1&-1&1\\\ 1&-1&1&2&0&-1&-1&1&1\\\
1&-1&-1&2&0&1&1&1&1\\\ 1&1&-1&2&0&-1&1&-1&1\\\ \end{array}\right],$ (4)

where we have used

$\mathcal{K}=\left[\ket{K_{0}},\ket{K_{1}},\ket{K_{2}},\ket{K_{3}},\ket{K_{4}},\ket{K_{5}},\ket{K_{6}},\ket{K_{7}},\ket{K_{8}}\right].$
(5)

##  3 Continuous Boltzmann Equation

We consider the two-dimensional (2D) continuous Boltzmann equation, for which
we aim to develop a Galilean invariant discrete model using the above basis
vectors. It represents the evolution of the continuous density distribution
function $f=f(x,y,\xi_{x},\xi_{y})$ in continuous phase space
$(x,y,\xi_{x},\xi_{y})$ subjected to a local force field
$\overrightarrow{F}=(F_{x},F_{y})$, whose origin could be internal or external
to the system. By definition, the averaged effects of $f$, weighted by various
powers of the continuous particle velocity $(\xi_{x},\xi_{y})$, i.e. its
moments, are considered to characterize the various physical processes
inherent in the motion of athermal fluids. In particular, the evolution of the
slow hydrodynamical processes are described by the local macroscopic fluid
density $\rho$ and fluid velocity $\overrightarrow{u}=(u_{x},u_{y})$. The
continuous Boltzmann equation may be written as [64, 84]

$\frac{\partial f}{\partial
t}+\overrightarrow{\xi}\cdot\overrightarrow{\nabla}_{\overrightarrow{x}}f+\frac{\overrightarrow{F}}{\rho}\cdot\overrightarrow{\nabla}_{\overrightarrow{\xi}}f=\Omega(f,f),$
(6)

where

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}fd\xi_{x}d\xi_{y},$
(7) $\displaystyle\rho\overrightarrow{u}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f\overrightarrow{\xi}d\xi_{x}d\xi_{y}.$
(8)

Here, $\Omega(f,f)$ is the collision term, which represents the cumulative
effect of binary collision of particles. The force fields modify the
distribution function exactly by the term
$-\frac{\overrightarrow{F}}{\rho}\cdot\overrightarrow{\nabla}_{\overrightarrow{\xi}}f$,
which is obtained by moving the last term on the left hand side of Eq. (6) to
its right to serve as a source term. It was shown by Grad (1949) [21] that
that solution of Eq. (6) can be approximated by the evolution equations for a
hierarchical set of moments. Here, we seek to obtain a dramatically
discretized version of this continuous Boltzmann equation by means of a
hierarchy of central moments, focusing, in particular, on the forcing term, to
obtain Galilean invariant representation of the dynamics of fluid motion.

##  4 Continuous Central Moments: Distribution Function and its Local
Attractor

We now consider the integral properties of the distribution function $f$ given
in terms of its central moments, i.e. those shifted by the macroscopic fluid
velocity. In particular, we define _continuous_ central moment of $f$ of order
$(m+n)$ as

$\widehat{\Pi}_{x^{m}y^{n}}=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(\xi_{x}-u_{x})^{m}(\xi_{y}-u_{y})^{n}d\xi_{x}d\xi_{y}.$
(9)

Here, and in the rest of this paper, the use of “hat” over a symbol represents
quantities in the space of moments. The distribution function for an athermal
fluid has a local equilibrium state in the _continuous_ particle velocity
space $(\xi_{x},\xi_{y})$, which is given by the Maxwellian as [84]

$f^{\mathcal{M}}\equiv
f^{\mathcal{M}}(\rho,\overrightarrow{u},\xi_{x},\xi_{y})=\frac{\rho}{2\pi
c_{s}^{2}}\exp\left[-\frac{\left(\overrightarrow{\xi}-\overrightarrow{u}\right)^{2}}{2c_{s}^{2}}\right],$
(10)

where $c_{s}^{2}=1/3$. Analogously, we can define the corresponding central
moment of the Maxwell distribution of order $(m+n)$ as

$\widehat{\Pi}^{\mathcal{M}}_{x^{m}y^{n}}=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f^{\mathcal{M}}(\xi_{x}-u_{x})^{m}(\xi_{y}-u_{y})^{n}d\xi_{x}d\xi_{y}.$
(11)

By virtue of the fact that $f^{\mathcal{M}}$ being an even function,
$\widehat{\Pi}^{\mathcal{M}}_{x^{m}y^{n}}\neq 0$ when $m$ and $n$ are even and
$\widehat{\Pi}^{\mathcal{M}}_{x^{m}y^{n}}=0$ when $m$ or $n$ odd. Here and
henceforth, the subscripts $x^{m}y^{n}$ mean $xxx\cdots m\mbox{-times}$ and
$yyy\cdots n\mbox{-times}$. Evaluation of the central moments of the
Maxwellian, to different orders of increasing powers, yields

$\displaystyle\widehat{\Pi}^{\mathcal{M}}_{0}$ $\displaystyle=$
$\displaystyle\rho,$ $\displaystyle\widehat{\Pi}^{\mathcal{M}}_{x}$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\widehat{\Pi}^{\mathcal{M}}_{y}$ $\displaystyle=$ $\displaystyle
0,$ $\displaystyle\widehat{\Pi}^{\mathcal{M}}_{xx}$ $\displaystyle=$
$\displaystyle c_{s}^{2}\rho,$ $\displaystyle\widehat{\Pi}^{\mathcal{M}}_{yy}$
$\displaystyle=$ $\displaystyle c_{s}^{2}\rho,$ (12)
$\displaystyle\widehat{\Pi}^{\mathcal{M}}_{xy}$ $\displaystyle=$
$\displaystyle 0,$ $\displaystyle\widehat{\Pi}^{\mathcal{M}}_{xxy}$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\widehat{\Pi}^{\mathcal{M}}_{xyy}$ $\displaystyle=$
$\displaystyle 0,$ $\displaystyle\widehat{\Pi}^{\mathcal{M}}_{xxyy}$
$\displaystyle=$ $\displaystyle c_{s}^{4}\rho.$

##  5 Continuous Central Moments: Forcing

In the presence of a force field $\overrightarrow{F}$, in view of Eq. (6) and
as discussed in Sec. 3, the distribution function will be exactly modified by
the source term

$\delta
f^{F}=-\frac{\overrightarrow{F}}{\rho}\cdot\overrightarrow{\nabla}_{\overrightarrow{\xi}}f.$
(13)

Now we define a corresponding _continuous_ central moment of order $(m+n)$ due
to change in the distribution function as a result of a force field as

$\widehat{\Gamma}^{F}_{x^{m}y^{n}}=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\delta
f^{F}(\xi_{x}-u_{x})^{m}(\xi_{y}-u_{y})^{n}d\xi_{x}d\xi_{y}.$ (14)

Substituting Eq. (13) in Eq. (14) and integrating by parts by making use of
the asymptotic limit assumptions
$lim_{\xi_{x}\rightarrow\pm\infty}(\xi_{x}-u_{x})^{m}f(x,y,\xi_{x},\xi_{y})=0$
and
$lim_{\xi_{y}\rightarrow\pm\infty}(\xi_{y}-u_{y})^{n}f(x,y,\xi_{x},\xi_{y})=0$,
for $m,n\geq 0$, we get

$\displaystyle\widehat{\Gamma}^{F}_{x^{m}y^{n}}$ $\displaystyle=$
$\displaystyle
m\frac{F_{x}}{\rho}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(\xi_{x}-u_{x})^{m-1}(\xi_{y}-u_{y})^{n}d\xi_{x}d\xi_{y}+$
(15) $\displaystyle
n\frac{F_{y}}{\rho}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(\xi_{x}-u_{x})^{m}(\xi_{y}-u_{y})^{n-1}d\xi_{x}d\xi_{y}.$

From the definition given in Eq. (9), Eq. (15) reduces to an _exact_ identity
between continuous central moment of the source term of a given order to the
components of the continuous central moment of the distribution function of an
order lower acted upon by a force field:

$\widehat{\Gamma}^{F}_{x^{m}y^{n}}=m\frac{F_{x}}{\rho}\widehat{\Pi}_{x^{m-1}y^{n}}+n\frac{F_{y}}{\rho}\widehat{\Pi}_{x^{m}y^{n-1}},$
(16)

and for the special case of the zeroth central moment of the source as
$\widehat{\Gamma}^{F}_{0}=0$. This is a key result based on which the rest of
the derivation follows. Thus, we can enumerate the _exact_ values of the
central moments of sources in an increasing order as

$\displaystyle\widehat{\Gamma}^{F}_{0}$ $\displaystyle=$ $\displaystyle 0,$
$\displaystyle\widehat{\Gamma}^{F}_{x}$ $\displaystyle=$
$\displaystyle\frac{F_{x}}{\rho}\widehat{\Pi}_{0},$
$\displaystyle\widehat{\Gamma}^{F}_{y}$ $\displaystyle=$
$\displaystyle\frac{F_{y}}{\rho}\widehat{\Pi}_{0},$
$\displaystyle\widehat{\Gamma}^{F}_{xx}$ $\displaystyle=$ $\displaystyle
2\frac{F_{x}}{\rho}\widehat{\Pi}_{x},$
$\displaystyle\widehat{\Gamma}^{F}_{yy}$ $\displaystyle=$ $\displaystyle
2\frac{F_{y}}{\rho}\widehat{\Pi}_{y},$ (17)
$\displaystyle\widehat{\Gamma}^{F}_{xy}$ $\displaystyle=$
$\displaystyle\frac{F_{x}}{\rho}\widehat{\Pi}_{y}+\frac{F_{y}}{\rho}\widehat{\Pi}_{x},$
$\displaystyle\widehat{\Gamma}^{F}_{xxy}$ $\displaystyle=$ $\displaystyle
2\frac{F_{x}}{\rho}\widehat{\Pi}_{xy}+\frac{F_{y}}{\rho}\widehat{\Pi}_{xx},$
$\displaystyle\widehat{\Gamma}^{F}_{xyy}$ $\displaystyle=$
$\displaystyle\frac{F_{x}}{\rho}\widehat{\Pi}_{yy}+2\frac{F_{y}}{\rho}\widehat{\Pi}_{xy},$
$\displaystyle\widehat{\Gamma}^{F}_{xxyy}$ $\displaystyle=$ $\displaystyle
2\frac{F_{x}}{\rho}\widehat{\Pi}_{xyy}+2\frac{F_{y}}{\rho}\widehat{\Pi}_{xxy}.$

Note that if we set
$\widehat{\Pi}_{x^{m}y^{n}}=\widehat{\Pi}_{x^{m}y^{n}}^{\mathcal{M}}$ in Eq.
(17), i.e. ignore non-equilibrium effects, we arrive at the the derivation
given in [63] as a special case.

##  6 Cascaded Central Moment Lattice-Boltzmann Method with Forcing Terms

Defining a _discrete_ distribution function supported by the discrete particle
velocity set $\overrightarrow{e}_{\alpha}$ as
$\mathbf{f}=\ket{f_{\alpha}}=(f_{0},f_{1},f_{2},\ldots,f_{8})^{T}$, and a
cascaded collision operator as
$\bm{\Omega}^{c}=\ket{\Omega_{\alpha}^{c}}=(\Omega_{0}^{c},\Omega_{1}^{c},\Omega_{2}^{c},\ldots,\Omega_{8}^{c})^{T}$
as well as a source term as
$\mathbf{S}=\ket{S_{\alpha}}=(S_{0},S_{1},S_{2},\ldots,S_{8})^{T}$ based on
central moments, we obtain the lattice Boltzmann equation (LBE) as a discrete
version of Eq. (6) by temporally integrating along particle characteristics as
follows [63]:

$f_{\alpha}(\overrightarrow{x}+\overrightarrow{e}_{\alpha},t+1)=f_{\alpha}(\overrightarrow{x},t)+\Omega_{{\alpha}(\overrightarrow{x},t)}^{c}+\int_{t}^{t+1}S_{{\alpha}(\overrightarrow{x}+\overrightarrow{e}_{\alpha}\theta,t+\theta)}d\theta.$
(18)

Here, the collision operator is written as

$\Omega_{\alpha}^{c}\equiv\Omega_{\alpha}^{c}(\mathbf{f},\mathbf{\widehat{g}})=(\mathcal{K}\cdot\mathbf{\widehat{g}})_{\alpha},$
(19)

where
$\mathbf{\widehat{g}}=\ket{\widehat{g}_{\alpha}}=(\widehat{g}_{0},\widehat{g}_{1},\widehat{g}_{2},\ldots,\widehat{g}_{8})^{T}$.
The hydrodynamic fields are obtained from the distribution function as

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\sum_{\alpha=0}^{8}f_{\alpha}=\braket{f_{\alpha}}{\rho},$ (20)
$\displaystyle\rho u_{i}$ $\displaystyle=$
$\displaystyle\sum_{\alpha=0}^{8}f_{\alpha}e_{\alpha
i}=\braket{f_{\alpha}}{e_{\alpha i}},i\in{x,y}.$ (21)

For improved accuracy in recovering Navier-Stokes solution, using a semi-
implicit representation for the source term, i.e. the last term in the above
equation (Eq. (18)) as
$\int_{t}^{t+1}S_{{\alpha}(\overrightarrow{x}+\overrightarrow{e}_{\alpha}\theta,t+\theta)}d\theta=\frac{1}{2}\left[S_{{\alpha}(\overrightarrow{x},t)}+S_{{\alpha}(\overrightarrow{x}+\overrightarrow{e}_{\alpha},t+1)}\right]$,
so that Eq. (18) is made effectively explicit by using the transformation
$\overline{f}_{\alpha}=f_{\alpha}-\frac{1}{2}S_{\alpha}$ to reduce it to [63]

$\overline{f}_{\alpha}(\overrightarrow{x}+\overrightarrow{e}_{\alpha},t+1)=\overline{f}_{\alpha}(\overrightarrow{x},t)+\Omega_{{\alpha}(\overrightarrow{x},t)}^{c}+S_{{\alpha}(\overrightarrow{x},t)}.$
(22)

It may be noted that, as in [59], we first represent collision as a cascaded
process in which the effect of collision on lower order central moments
successively influence those at higher orders in a cascaded manner. That is,
in general,
$\widehat{g}_{\alpha}\equiv\widehat{g}_{\alpha}(\mathbf{f},\widehat{g}_{\beta}),\beta=0,1,2,\ldots,\alpha-1$.
Furthermore, the form of the source term is derived to rigorously enforce
Galilean invariance. The explicit expressions for $S_{\alpha}$ and
$\mathbf{\widehat{g}}$ will be determined later in Secs. 9 and 10,
respectively. Since the main focus of this work is on improving the collision
(including forcing) step with features independent of inertial frames, we have
only considered the standard discretization for the streaming operator.
However, as discussed in the Introduction, other types of discretization
schemes could be considered to improve the order of accuracy.

##  7 Various Discrete Central Moments and Galilean Invariance Matching
Principle

For determining the structure of the cascaded collision operator
$\mathbf{\widehat{g}}$ and the source terms $S_{\alpha}$, we first need to
define the following _discrete_ central moments of the distribution function,
Maxwellian, and source term, respectively:

$\displaystyle\widehat{\kappa}_{x^{m}y^{n}}$ $\displaystyle=$
$\displaystyle\sum_{\alpha}f_{\alpha}(e_{\alpha x}-u_{x})^{m}(e_{\alpha
y}-u_{y})^{n}=\braket{(e_{\alpha x}-u_{x})^{m}(e_{\alpha
y}-u_{y})^{n}}{f_{\alpha}},$ (23)
$\displaystyle\widehat{\kappa}_{x^{m}y^{n}}^{\mathcal{M}}$ $\displaystyle=$
$\displaystyle\sum_{\alpha}f_{\alpha}^{\mathcal{M}}(e_{\alpha
x}-u_{x})^{m}(e_{\alpha y}-u_{y})^{n}=\braket{(e_{\alpha
x}-u_{x})^{m}(e_{\alpha y}-u_{y})^{n}}{f_{\alpha}^{\mathcal{M}}},$ (24)
$\displaystyle\widehat{\sigma}_{x^{m}y^{n}}$ $\displaystyle=$
$\displaystyle\sum_{\alpha}S_{\alpha}(e_{\alpha x}-u_{x})^{m}(e_{\alpha
y}-u_{y})^{n}=\braket{(e_{\alpha x}-u_{x})^{m}(e_{\alpha
y}-u_{y})^{n}}{S_{\alpha}}.$ (25)

where the exact expression for the discrete $f_{\alpha}^{\mathcal{M}}$ is not
yet known, but can be determined as a result of the derivation discussed
later. To maintain physical consistency at the discrete level, we now equate
the _discrete_ central moments of the distribution function, the Maxwellian
and the source terms, defined above, with their corresponding _continuous_
central moments, whose forms are known exactly. That is, according to this
matching principle

$\displaystyle\widehat{\kappa}_{x^{m}y^{n}}$ $\displaystyle=$
$\displaystyle\widehat{\Pi}_{x^{m}y^{n}},$ (26)
$\displaystyle\widehat{\kappa}_{x^{m}y^{n}}^{\mathcal{M}}$ $\displaystyle=$
$\displaystyle\widehat{\Pi}^{\mathcal{M}}_{x^{m}y^{n}},$ (27)
$\displaystyle\widehat{\sigma}_{x^{m}y^{n}}$ $\displaystyle=$
$\displaystyle\widehat{\Gamma}^{F}_{x^{m}y^{n}}.$ (28)

In particular, the discrete central moments of various orders for both the
Maxwellian and the source terms, respectively, become

$\displaystyle\widehat{\kappa}^{\mathcal{M}}_{0}$ $\displaystyle=$
$\displaystyle\rho,$ $\displaystyle\widehat{\kappa}^{\mathcal{M}}_{x}$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\widehat{\kappa}^{\mathcal{M}}_{y}$ $\displaystyle=$
$\displaystyle 0,$ $\displaystyle\widehat{\kappa}^{\mathcal{M}}_{xx}$
$\displaystyle=$ $\displaystyle c_{s}^{2}\rho,$
$\displaystyle\widehat{\kappa}^{\mathcal{M}}_{yy}$ $\displaystyle=$
$\displaystyle c_{s}^{2}\rho,$ (29)
$\displaystyle\widehat{\kappa}^{\mathcal{M}}_{xy}$ $\displaystyle=$
$\displaystyle 0,$ $\displaystyle\widehat{\kappa}^{\mathcal{M}}_{xxy}$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\widehat{\kappa}^{\mathcal{M}}_{xyy}$ $\displaystyle=$
$\displaystyle 0,$ $\displaystyle\widehat{\kappa}^{\mathcal{M}}_{xxyy}$
$\displaystyle=$ $\displaystyle c_{s}^{4}\rho,$

and

$\displaystyle\widehat{\sigma}_{0}$ $\displaystyle=$ $\displaystyle 0,$
$\displaystyle\widehat{\sigma}_{x}$ $\displaystyle=$
$\displaystyle\frac{F_{x}}{\rho}\widehat{\kappa}_{0},$
$\displaystyle\widehat{\sigma}_{y}$ $\displaystyle=$
$\displaystyle\frac{F_{y}}{\rho}\widehat{\kappa}_{0},$
$\displaystyle\widehat{\sigma}_{xx}$ $\displaystyle=$ $\displaystyle
2\frac{F_{x}}{\rho}\widehat{\kappa}_{x},$ $\displaystyle\widehat{\sigma}_{yy}$
$\displaystyle=$ $\displaystyle 2\frac{F_{y}}{\rho}\widehat{\kappa}_{y},$ (30)
$\displaystyle\widehat{\sigma}_{xy}$ $\displaystyle=$
$\displaystyle\frac{F_{x}}{\rho}\widehat{\kappa}_{y}+\frac{F_{y}}{\rho}\widehat{\kappa}_{x},$
$\displaystyle\widehat{\sigma}_{xxy}$ $\displaystyle=$ $\displaystyle
2\frac{F_{x}}{\rho}\widehat{\kappa}_{xy}+\frac{F_{y}}{\rho}\widehat{\kappa}_{xx},$
$\displaystyle\widehat{\sigma}_{xyy}$ $\displaystyle=$
$\displaystyle\frac{F_{x}}{\rho}\widehat{\kappa}_{yy}+2\frac{F_{y}}{\rho}\widehat{\kappa}_{xy},$
$\displaystyle\widehat{\sigma}_{xxyy}$ $\displaystyle=$ $\displaystyle
2\frac{F_{x}}{\rho}\widehat{\kappa}_{xyy}+2\frac{F_{y}}{\rho}\widehat{\kappa}_{xxy}.$

We also define a _discrete_ central moment in terms of the transformed
distribution function to facilitate subsequent developments as

$\widehat{\overline{\kappa}}_{x^{m}y^{n}}=\sum_{\alpha}\overline{f}_{\alpha}(e_{\alpha
x}-u_{x})^{m}(e_{\alpha y}-u_{y})^{n}=\braket{(e_{\alpha
x}-u_{x})^{m}(e_{\alpha y}-u_{y})^{n}}{\overline{f}_{\alpha}}.$ (31)

Owing to the transformation discussed in Sec. 6, it follows that

$\widehat{\overline{\kappa}}_{x^{m}y^{n}}=\widehat{\kappa}_{x^{m}y^{n}}-\frac{1}{2}\widehat{\sigma}_{x^{m}y^{n}}.$
(32)

##  8 Relation Between Various Discrete Central Moments

Equation (30) is given in terms of the discrete moments of the original
distribution function $f_{\alpha}$. However, the cascaded central moment LBM
with forcing term provides evolution in terms of transformed distribution
function $\overline{f}_{\alpha}$ (Eq. (22)). Thus, it is important to write
all the expressions in terms of the central moments of
$\overline{f}_{\alpha}$, or, equivalently,
$\widehat{\overline{\kappa}}_{x^{m}y^{n}}$. Thus, by recursive application of
Eq. (32) using Eq. (30) to successively higher orders, we get the following
exact relations up to the third-order central moments as

$\displaystyle\widehat{\kappa}_{0}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{0},$
$\displaystyle\widehat{\kappa}_{x}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{x}+\frac{1}{2}\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{0},$
$\displaystyle\widehat{\kappa}_{y}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{y}+\frac{1}{2}\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{0},$
$\displaystyle\widehat{\kappa}_{xx}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{xx}+\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{x}+\frac{1}{2}\frac{F_{x}^{2}}{\rho^{2}}\widehat{\overline{\kappa}}_{0},$
$\displaystyle\widehat{\kappa}_{yy}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{yy}+\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{y}+\frac{1}{2}\frac{F_{y}^{2}}{\rho^{2}}\widehat{\overline{\kappa}}_{0},$
$\displaystyle\widehat{\kappa}_{xy}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{xy}+\frac{1}{2}\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{y}+\frac{1}{2}\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{x}+\frac{1}{2}\frac{F_{x}F_{y}}{\rho^{2}}\widehat{\overline{\kappa}}_{0},$
$\displaystyle\widehat{\kappa}_{xxy}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{xxy}+\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{xy}+\frac{1}{2}\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{xx}+\frac{F_{x}F_{y}}{\rho^{2}}\widehat{\overline{\kappa}}_{x}+\frac{1}{2}\frac{F_{x}^{2}}{\rho^{2}}\widehat{\overline{\kappa}}_{y}+\frac{3}{4}\frac{F_{x}^{2}F_{y}}{\rho^{3}}\widehat{\overline{\kappa}}_{0},$
$\displaystyle\widehat{\kappa}_{xyy}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{xyy}+\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{xy}+\frac{1}{2}\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{yy}+\frac{F_{x}F_{y}}{\rho^{2}}\widehat{\overline{\kappa}}_{y}+\frac{1}{2}\frac{F_{y}^{2}}{\rho^{2}}\widehat{\overline{\kappa}}_{x}+\frac{3}{4}\frac{F_{x}F_{y}^{2}}{\rho^{3}}\widehat{\overline{\kappa}}_{0}.$

That is, the central moment of the distribution function of a given order can
be written as a function of the central moment of the transformed distribution
function of the same order and successively lower orders as well. This can be
further simplified by considering the three of the lowest order central
moments, i.e., conservative moments, which by definition are
$\widehat{\kappa}_{0}=\rho$, $\widehat{\kappa}_{x}=\widehat{\kappa}_{y}=0$.
This, in turn, leads to $\widehat{\overline{\kappa}}_{0}=\rho$,
$\widehat{\overline{\kappa}}_{x}=-1/2F_{x}$,
$\widehat{\overline{\kappa}}_{y}=-1/2F_{y}$. As a result, we have the
following relations for the non-conserved central moments up to third-order:

$\displaystyle\widehat{\kappa}_{xx}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{xx},$
$\displaystyle\widehat{\kappa}_{yy}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{yy},$
$\displaystyle\widehat{\kappa}_{xy}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{xy},$
$\displaystyle\widehat{\kappa}_{xxy}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{xxy}+\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{xy}+\frac{1}{2}\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{xx},$
$\displaystyle\widehat{\kappa}_{xyy}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{xyy}+\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{xy}+\frac{1}{2}\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{yy}.$

Thus, we can finally write the central moments of the source term in Eq. (30)
in terms of the central moments of the transformed distribution function as

$\displaystyle\widehat{\sigma}_{0}$ $\displaystyle=$ $\displaystyle 0,$
$\displaystyle\widehat{\sigma}_{x}$ $\displaystyle=$ $\displaystyle F_{x},$
$\displaystyle\widehat{\sigma}_{y}$ $\displaystyle=$ $\displaystyle F_{y},$
$\displaystyle\widehat{\sigma}_{xx}$ $\displaystyle=$ $\displaystyle 0,$
$\displaystyle\widehat{\sigma}_{yy}$ $\displaystyle=$ $\displaystyle 0,$ (33)
$\displaystyle\widehat{\sigma}_{xy}$ $\displaystyle=$ $\displaystyle 0,$
$\displaystyle\widehat{\sigma}_{xxy}$ $\displaystyle=$ $\displaystyle
2\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{xy}+\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{xx},$
$\displaystyle\widehat{\sigma}_{xyy}$ $\displaystyle=$
$\displaystyle\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{yy}+2\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{xy},$
$\displaystyle\widehat{\sigma}_{xxyy}$ $\displaystyle=$ $\displaystyle
2\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{xyy}+2\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{xxy}+\frac{F_{y}^{2}}{\rho^{2}}\widehat{\overline{\kappa}}_{xx}+\frac{F_{x}^{2}}{\rho^{2}}\widehat{\overline{\kappa}}_{yy}+4\frac{F_{x}F_{y}}{\rho^{2}}\widehat{\overline{\kappa}}_{xy}.$

Thus, higher order non-equilibrium effects in
$\widehat{\overline{\kappa}}_{x^{m}y^{n}}$ and non-linear effect in
$F_{x}^{p}F_{y}^{q}$ are evident for the central moments of the source terms
that are third- and higher orders. Let us now explicitly write the central
moments of the transformed discrete Maxwellian by means of Eq. (32) using Eqs.
(29) and (33) to yield

$\displaystyle\widehat{\overline{\kappa}}^{\mathcal{M}}_{0}$ $\displaystyle=$
$\displaystyle\rho,$
$\displaystyle\widehat{\overline{\kappa}}^{\mathcal{M}}_{x}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}F_{x},$
$\displaystyle\widehat{\overline{\kappa}}^{\mathcal{M}}_{y}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}F_{y},$
$\displaystyle\widehat{\overline{\kappa}}^{\mathcal{M}}_{xx}$ $\displaystyle=$
$\displaystyle c_{s}^{2}\rho,$
$\displaystyle\widehat{\overline{\kappa}}^{\mathcal{M}}_{yy}$ $\displaystyle=$
$\displaystyle c_{s}^{2}\rho,$ (34)
$\displaystyle\widehat{\overline{\kappa}}^{\mathcal{M}}_{xy}$ $\displaystyle=$
$\displaystyle 0,$
$\displaystyle\widehat{\overline{\kappa}}^{\mathcal{M}}_{xxy}$
$\displaystyle=$
$\displaystyle-\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{xy}-\frac{F_{y}}{2\rho}\widehat{\overline{\kappa}}_{xx},$
$\displaystyle\widehat{\overline{\kappa}}^{\mathcal{M}}_{xyy}$
$\displaystyle=$
$\displaystyle-\frac{F_{x}}{2\rho}\widehat{\overline{\kappa}}_{yy}-\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{xy},$
$\displaystyle\widehat{\overline{\kappa}}^{\mathcal{M}}_{xxyy}$
$\displaystyle=$ $\displaystyle
c_{s}^{4}\rho-\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{xyy}-\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{xxy}-\frac{F_{y}^{2}}{2\rho^{2}}\widehat{\overline{\kappa}}_{xx}-\frac{F_{x}^{2}}{2\rho^{2}}\widehat{\overline{\kappa}}_{yy}-2\frac{F_{x}F_{y}}{\rho^{2}}\widehat{\overline{\kappa}}_{xy}.$

The main idea in the determination of the collision operator for the cascaded
version of the central moment method is to relax the central moments of the
transformed distribution function to its corresponding local attractor,
successively at various orders as given in Eq. (34) (see Sec. 10). Before
proceeding further to do this, we first need certain quantities in the rest or
lattice frame of reference, i.e. the raw moments, in which the computations
are actually performed. These are obtained in the next section.

##  9 Various Discrete Raw Moments and Source Terms in Particle Velocity Space

The raw moments, i.e. those in the rest frame of reference, can be related to
the central moments by means of the binomial theorem [85, 86]. For any state
variable $\varphi$ supported by the discrete particle velocity set, the
transformation relation between the two reference frames is thus given by [63]

$\displaystyle\braket{(e_{\alpha x}-u_{x})^{m}(e_{\alpha
y}-u_{y})^{n}}{\varphi}$ $\displaystyle=$ $\displaystyle\braket{e_{\alpha
x}^{m}e_{\alpha y}^{n}}{\varphi}+\braket{e_{\alpha
x}^{m}\left[\sum_{j=1}^{n}C^{n}_{j}e_{\alpha
y}^{n-j}(-1)^{j}u_{y}^{j}\right]}{\varphi}+$ (35)
$\displaystyle\braket{e_{\alpha y}^{n}\left[\sum_{i=1}^{m}C^{m}_{i}e_{\alpha
x}^{m-i}(-1)^{i}u_{x}^{i}\right]}{\varphi}+$
$\displaystyle\braket{\left[\sum_{i=1}^{m}C^{m}_{i}e_{\alpha
x}^{m-i}(-1)^{i}u_{x}^{i}\right]\left[\sum_{j=1}^{n}C^{n}_{j}e_{\alpha
y}^{n-j}(-1)^{j}u_{y}^{j}\right]}{\varphi}$

where $C^{p}_{q}=p!/(q!(p-q)!)$. We now define the following notations for
depicting various _discrete raw_ moments, based on which an operational LBE
will be devised later:

$\displaystyle\widehat{\kappa}_{x^{m}y^{n}}^{{}^{\prime}}$ $\displaystyle=$
$\displaystyle\sum_{\alpha}f_{\alpha}e_{\alpha x}^{m}e_{\alpha
y}^{n}=\braket{e_{\alpha x}^{m}e_{\alpha y}^{n}}{f_{\alpha}},$ (36)
$\displaystyle\widehat{\overline{\kappa}}_{x^{m}y^{n}}^{{}^{\prime}}$
$\displaystyle=$ $\displaystyle\sum_{\alpha}\overline{f}_{\alpha}e_{\alpha
x}^{m}e_{\alpha y}^{n}=\braket{e_{\alpha x}^{m}e_{\alpha
y}^{n}}{\overline{f}_{\alpha}},$ (37)
$\displaystyle\widehat{\sigma}_{x^{m}y^{n}}^{{}^{\prime}}$ $\displaystyle=$
$\displaystyle\sum_{\alpha}S_{\alpha}e_{\alpha x}^{m}e_{\alpha
y}^{n}=\braket{e_{\alpha x}^{m}e_{\alpha y}^{n}}{S_{\alpha}}.$ (38)

Note that the superscript “prime” (′) is used to distinguish the raw moments
from the central moments that are designated without the primes. Here,
analogous to Eq. (32), the relation
$\widehat{\overline{\kappa}}_{x^{m}y^{n}}^{{}^{\prime}}=\widehat{\kappa}_{x^{m}y^{n}}^{{}^{\prime}}-\frac{1}{2}\widehat{\sigma}_{x^{m}y^{n}}^{{}^{\prime}}$
holds. Let us first find expressions for
$\widehat{\overline{\kappa}}_{x^{m}y^{n}}^{{}^{\prime}}=\braket{\overline{f}_{\alpha}}{e_{\alpha
x}^{m}e_{\alpha y}^{n}}$ to proceed further. As in [63], for convenience, we
define the following operator acting on the transformed distribution function
$\overline{f}_{\alpha}$ in this regard:

$a(\overline{f}_{\alpha_{1}}+\overline{f}_{\alpha_{3}}+\overline{f}_{\alpha_{3}}+\cdots)+b(\overline{f}_{\beta_{1}}+\overline{f}_{\beta_{2}}+\overline{f}_{\beta_{3}}+\cdots)+\cdots=\left(a\sum_{\alpha}^{A}+b\sum_{\alpha}^{B}\cdots\right)\otimes\overline{f}_{\alpha},$
(39)

where $A=\left\\{\alpha_{1},\alpha_{2},\alpha_{3},\cdots\right\\}$,
$B=\left\\{\beta_{1},\beta_{2},\beta_{3},\cdots\right\\}$,$\cdots$. For
conserved basis vectors, we write them in terms of the hydrodynamic variables
and force fields as

$\displaystyle\widehat{\overline{\kappa}}_{0}^{{}^{\prime}}$ $\displaystyle=$
$\displaystyle\braket{\overline{f}_{\alpha}}{\rho}=\sum_{\alpha=0}^{8}\overline{f}_{\alpha}=\rho,$
(40) $\displaystyle\widehat{\overline{\kappa}}_{x}^{{}^{\prime}}$
$\displaystyle=$ $\displaystyle\braket{\overline{f}_{\alpha}}{e_{\alpha
x}}=\sum_{\alpha=0}^{8}\overline{f}_{\alpha}e_{\alpha x}=\rho
u_{x}-\frac{1}{2}F_{x},$ (41)
$\displaystyle\widehat{\overline{\kappa}}_{y}^{{}^{\prime}}$ $\displaystyle=$
$\displaystyle\braket{\overline{f}_{\alpha}}{e_{\alpha
y}}=\sum_{\alpha=0}^{8}\overline{f}_{\alpha}e_{\alpha y}=\rho
u_{y}-\frac{1}{2}F_{y},$ (42)

and, for the non-conserved basis vectors, using Eq. (39) in terms of subsets
of particle velocity directions as

$\displaystyle\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}$ $\displaystyle=$
$\displaystyle\braket{\overline{f}_{\alpha}}{e_{\alpha
x}^{2}}=\sum_{\alpha=0}^{8}\overline{f}_{\alpha}e_{\alpha
x}^{2}=\left(\sum_{\alpha}^{A_{3}}\right)\otimes\overline{f}_{\alpha},$ (43)
$\displaystyle\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}}$ $\displaystyle=$
$\displaystyle\braket{\overline{f}_{\alpha}}{e_{\alpha
y}^{2}}=\sum_{\alpha=0}^{8}\overline{f}_{\alpha}e_{\alpha
y}^{2}=\left(\sum_{\alpha}^{A_{4}}\right)\otimes\overline{f}_{\alpha},$ (44)
$\displaystyle\widehat{\overline{\kappa}}_{xy}^{{}^{\prime}}$ $\displaystyle=$
$\displaystyle\braket{\overline{f}_{\alpha}}{e_{\alpha x}e_{\alpha
y}}=\sum_{\alpha=0}^{8}\overline{f}_{\alpha}e_{\alpha x}e_{\alpha
y}=\left(\sum_{\alpha}^{A_{5}}-\sum_{\alpha}^{B_{5}}\right)\otimes\overline{f}_{\alpha},$
(45) $\displaystyle\widehat{\overline{\kappa}}_{xxy}^{{}^{\prime}}$
$\displaystyle=$ $\displaystyle\braket{\overline{f}_{\alpha}}{e_{\alpha
x}^{2}e_{\alpha y}}=\sum_{\alpha=0}^{8}\overline{f}_{\alpha}e_{\alpha
x}^{2}e_{\alpha
y}=\left(\sum_{\alpha}^{A_{6}}-\sum_{\alpha}^{B_{6}}\right)\otimes\overline{f}_{\alpha},$
(46) $\displaystyle\widehat{\overline{\kappa}}_{xyy}^{{}^{\prime}}$
$\displaystyle=$ $\displaystyle\braket{\overline{f}_{\alpha}}{e_{\alpha
x}e_{\alpha y}^{2}}=\sum_{\alpha=0}^{8}\overline{f}_{\alpha}e_{\alpha
x}e_{\alpha
y}^{2}=\left(\sum_{\alpha}^{A_{7}}-\sum_{\alpha}^{B_{7}}\right)\otimes\overline{f}_{\alpha},$
(47) $\displaystyle\widehat{\overline{\kappa}}_{xxyy}^{{}^{\prime}}$
$\displaystyle=$ $\displaystyle\braket{\overline{f}_{\alpha}}{e_{\alpha
x}^{2}e_{\alpha y}^{2}}=\sum_{\alpha=0}^{8}\overline{f}_{\alpha}e_{\alpha
x}^{2}e_{\alpha
y}^{2}=\left(\sum_{\alpha}^{A_{8}}\right)\otimes\overline{f}_{\alpha},$ (48)

where

$\displaystyle A_{3}$ $\displaystyle=$
$\displaystyle\left\\{1,3,5,6,7,8\right\\},$ $\displaystyle A_{4}$
$\displaystyle=$ $\displaystyle\left\\{2,4,5,6,7,8\right\\},$ $\displaystyle
A_{5}$ $\displaystyle=$
$\displaystyle\left\\{5,7\right\\},B_{5}=\left\\{6,8\right\\},$ $\displaystyle
A_{6}$ $\displaystyle=$
$\displaystyle\left\\{5,6\right\\},B_{6}=\left\\{7,8\right\\},$ $\displaystyle
A_{7}$ $\displaystyle=$
$\displaystyle\left\\{5,8\right\\},B_{7}=\left\\{6,7\right\\},$ $\displaystyle
A_{8}$ $\displaystyle=$ $\displaystyle\left\\{5,6,7,8\right\\}.$

Now, we transform the central moments of the source terms (Eq. (30)) to the
corresponding raw moments by considering Eq. (25) and using the frame
transformation relation (Eq. (35)). This yields

$\displaystyle\widehat{\sigma}_{0}^{{}^{\prime}}=\braket{S_{\alpha}}{\rho}=0,$
(49)
$\displaystyle\widehat{\sigma}_{x}^{{}^{\prime}}=\braket{S_{\alpha}}{e_{\alpha
x}}=F_{x},$ (50)
$\displaystyle\widehat{\sigma}_{y}^{{}^{\prime}}=\braket{S_{\alpha}}{e_{\alpha
y}}=F_{y},$ (51)
$\displaystyle\widehat{\sigma}_{xx}^{{}^{\prime}}=\braket{S_{\alpha}}{e_{\alpha
x}^{2}}=2F_{x}u_{x},$ (52)
$\displaystyle\widehat{\sigma}_{yy}^{{}^{\prime}}=\braket{S_{\alpha}}{e_{\alpha
y}^{2}}=2F_{x}u_{y},$ (53)
$\displaystyle\widehat{\sigma}_{xy}^{{}^{\prime}}=\braket{S_{\alpha}}{e_{\alpha
x}e_{\alpha y}}=F_{x}u_{y}+F_{y}u_{x},$ (54)
$\displaystyle\widehat{\sigma}_{xxy}^{{}^{\prime}}=\braket{S_{\alpha}}{e_{\alpha
x}^{2}e_{\alpha
y}}=2F_{x}\left(\frac{\widehat{\overline{\kappa}}_{xy}}{\rho}+u_{x}u_{y}\right)+F_{y}\left(\frac{\widehat{\overline{\kappa}}_{xx}}{\rho}+u_{x}^{2}\right),$
(55)
$\displaystyle\widehat{\sigma}_{xyy}^{{}^{\prime}}=\braket{S_{\alpha}}{e_{\alpha
x}e_{\alpha
y}^{2}}=F_{x}\left(\frac{\widehat{\overline{\kappa}}_{yy}}{\rho}+u_{y}^{2}\right)+2F_{y}\left(\frac{\widehat{\overline{\kappa}}_{xy}}{\rho}+u_{x}u_{y}\right),$
(56)
$\displaystyle\widehat{\sigma}_{xxyy}^{{}^{\prime}}=\braket{S_{\alpha}}{e_{\alpha
x}^{2}e_{\alpha
y}^{2}}=2F_{x}u_{x}\left(\frac{\widehat{\overline{\kappa}}_{yy}}{\rho}+u_{y}^{2}\right)+2F_{y}u_{y}\left(\frac{\widehat{\overline{\kappa}}_{xx}}{\rho}+u_{x}^{2}\right)+$
$\displaystyle\frac{2F_{x}}{\rho}\widehat{\overline{\kappa}}_{xyy}+\frac{2F_{y}}{\rho}\widehat{\overline{\kappa}}_{xxy}+\frac{F_{x}^{2}}{\rho^{2}}\widehat{\overline{\kappa}}_{yy}+\frac{F_{y}^{2}}{\rho^{2}}\widehat{\overline{\kappa}}_{xx}+$
$\displaystyle
4\left[\frac{F_{x}u_{y}}{\rho}+\frac{F_{y}u_{x}}{\rho}+\frac{F_{x}F_{y}}{\rho^{2}}\right]\widehat{\overline{\kappa}}_{xy}.$
(57)

Clearly, the raw moments of source terms for third-order or higher contain
non-equilibrium and non-linear contributions. Eqs. (55)-(57) require explicit
expressions for central moments of transformed distributions such as
$\widehat{\overline{\kappa}}_{xx}$, $\widehat{\overline{\kappa}}_{yy}$,
$\widehat{\overline{\kappa}}_{xy}$, $\widehat{\overline{\kappa}}_{xxy}$ and
$\widehat{\overline{\kappa}}_{xyy}$, in terms of raw moments to facilitate
computation. They can be readily obtained in terms of raw moments from their
respective definitions and by using the binomial theorem (Eq. (35)) and
subsequent simplification as

$\displaystyle\widehat{\overline{\kappa}}_{xx}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}+F_{x}u_{x}-\rho
u_{x}^{2},$ (58) $\displaystyle\widehat{\overline{\kappa}}_{yy}$
$\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}}+F_{y}u_{y}-\rho
u_{y}^{2},$ (59) $\displaystyle\widehat{\overline{\kappa}}_{xy}$
$\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{xy}^{{}^{\prime}}+\frac{1}{2}(F_{x}u_{y}+F_{y}u_{x})-\rho
u_{x}u_{y},$ (60)

for second-order and

$\displaystyle\widehat{\overline{\kappa}}_{xyy}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{xyy}^{{}^{\prime}}-2u_{y}\widehat{\overline{\kappa}}_{xy}^{{}^{\prime}}-u_{x}\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}}-\frac{1}{2}F_{x}u_{y}^{2}-F_{y}u_{x}u_{y}+2\rho
u_{x}u_{y}^{2},$ (61) $\displaystyle\widehat{\overline{\kappa}}_{xxy}$
$\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{xxy}^{{}^{\prime}}-2u_{x}\widehat{\overline{\kappa}}_{xy}^{{}^{\prime}}-u_{y}\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}-\frac{1}{2}F_{y}u_{x}^{2}-F_{x}u_{x}u_{y}+2\rho
u_{x}^{2}u_{y},$ (62)

for third-order moments. Based on the above, we now obtain the source terms
projected to the orthogonal moment basis vectors, i.e.
$\braket{K_{\beta}}{S_{\alpha}}$, $\beta=0,1,2,\ldots,8$, which would then
provide corresponding explicit expressions in terms of the particle velocity
space. Thus, from Eqs. (5) and (49)-(57), the following projected source
moments are derived:

$\displaystyle\widehat{m}^{s}_{0}=\braket{K_{0}}{S_{\alpha}}$ $\displaystyle=$
$\displaystyle 0,$ (63)
$\displaystyle\widehat{m}^{s}_{1}=\braket{K_{1}}{S_{\alpha}}$ $\displaystyle=$
$\displaystyle F_{x},$ (64)
$\displaystyle\widehat{m}^{s}_{2}=\braket{K_{2}}{S_{\alpha}}$ $\displaystyle=$
$\displaystyle F_{y},$ (65)
$\displaystyle\widehat{m}^{s}_{3}=\braket{K_{3}}{S_{\alpha}}$ $\displaystyle=$
$\displaystyle 6(F_{x}u_{x}+F_{y}u_{y}),$ (66)
$\displaystyle\widehat{m}^{s}_{4}=\braket{K_{4}}{S_{\alpha}}$ $\displaystyle=$
$\displaystyle 2(F_{x}u_{x}-F_{y}u_{y}),$ (67)
$\displaystyle\widehat{m}^{s}_{5}=\braket{K_{5}}{S_{\alpha}}$ $\displaystyle=$
$\displaystyle(F_{x}u_{y}+F_{y}u_{x}),$ (68)
$\displaystyle\widehat{m}^{s}_{6}=\braket{K_{6}}{S_{\alpha}}$ $\displaystyle=$
$\displaystyle-6F_{x}\left(\frac{\widehat{\overline{\kappa}}_{xy}}{\rho}+u_{x}u_{y}\right)-3F_{y}\left(\frac{\widehat{\overline{\kappa}}_{xx}}{\rho}+u_{x}^{2}-\frac{2}{3}\right),$
(69) $\displaystyle\widehat{m}^{s}_{7}=\braket{K_{7}}{S_{\alpha}}$
$\displaystyle=$
$\displaystyle-3F_{x}\left(\frac{\widehat{\overline{\kappa}}_{yy}}{\rho}+u_{y}^{2}-\frac{2}{3}\right)-6F_{y}\left(\frac{\widehat{\overline{\kappa}}_{xy}}{\rho}+u_{x}u_{y}\right),$
(70) $\displaystyle\widehat{m}^{s}_{8}=\braket{K_{8}}{S_{\alpha}}$
$\displaystyle=$ $\displaystyle
18F_{x}u_{x}\left(\frac{\widehat{\overline{\kappa}}_{yy}}{\rho}+u_{y}^{2}-\frac{2}{3}\right)+18F_{y}u_{y}\left(\frac{\widehat{\overline{\kappa}}_{xx}}{\rho}+u_{x}^{2}-\frac{2}{3}\right)+$
(71) $\displaystyle
18\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{xyy}+18\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{xxy}+9\frac{F_{x}^{2}}{\rho^{2}}\widehat{\overline{\kappa}}_{yy}+9\frac{F_{y}^{2}}{\rho^{2}}\widehat{\overline{\kappa}}_{xx}+$
$\displaystyle
36\left[\frac{F_{x}u_{y}}{\rho}+\frac{F_{y}u_{x}}{\rho}+\frac{F_{x}F_{y}}{\rho^{2}}\right]\widehat{\overline{\kappa}}_{xy}.$

In Eqs. (69)-(71), $\widehat{\overline{\kappa}}_{xx}$,
$\widehat{\overline{\kappa}}_{yy}$, $\widehat{\overline{\kappa}}_{xy}$,
$\widehat{\overline{\kappa}}_{xyy}$ and $\widehat{\overline{\kappa}}_{xxy}$
can be obtained from Eqs. (58)-(62), respectively. This can be written in
matrix form as

$\displaystyle\mathcal{K}^{T}\mathbf{S}=(\mathcal{K}\cdot\mathbf{S})_{\alpha}$
$\displaystyle=$
$\displaystyle(\braket{K_{0}}{S_{\alpha}},\braket{K_{1}}{S_{\alpha}},\braket{K_{2}}{S_{\alpha}},\ldots,\braket{K_{8}}{S_{\alpha}})$
(72) $\displaystyle=$
$\displaystyle(\widehat{m}^{s}_{0},\widehat{m}^{s}_{1},\widehat{m}^{s}_{2},\ldots,\widehat{m}^{s}_{8})^{T}\equiv\ket{\widehat{m}^{s}_{\alpha}}.$

Now, by exploiting the orthogonal property of $\mathcal{K}$ [63], i.e.
$\mathcal{K}^{-1}=\mathcal{K}^{T}\cdot D^{-1}$, where the diagonal matrix is
$D=\mbox{diag}(\braket{K_{0}}{K_{0}},\braket{K_{1}}{K_{1}},\braket{K_{2}}{K_{2}},\ldots,\braket{K_{8}}{K_{8}})=\mbox{diag}(9,6,6,36,4,4,12,12,36)$,
we exactly invert Eq. (72) to finally obtain source terms in velocity space
$S_{\alpha}$ as

$\displaystyle S_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{9}\left[-\widehat{m}^{s}_{3}+\widehat{m}^{s}_{8}\right],$
(73) $\displaystyle S_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{36}\left[6\widehat{m}^{s}_{1}-\widehat{m}^{s}_{3}+9\widehat{m}^{s}_{4}+6\widehat{m}^{s}_{7}-2\widehat{m}^{s}_{8}\right],$
(74) $\displaystyle S_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{36}\left[6\widehat{m}^{s}_{2}-\widehat{m}^{s}_{3}-9\widehat{m}^{s}_{4}+6\widehat{m}^{s}_{6}-2\widehat{m}^{s}_{8}\right],$
(75) $\displaystyle S_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{36}\left[-6\widehat{m}^{s}_{1}-\widehat{m}^{s}_{3}+9\widehat{m}^{s}_{4}-6\widehat{m}^{s}_{7}-2\widehat{m}^{s}_{8}\right],$
(76) $\displaystyle S_{4}$ $\displaystyle=$
$\displaystyle\frac{1}{36}\left[-6\widehat{m}^{s}_{2}-\widehat{m}^{s}_{3}-9\widehat{m}^{s}_{4}-6\widehat{m}^{s}_{6}-2\widehat{m}^{s}_{8}\right],$
(77) $\displaystyle S_{5}$ $\displaystyle=$
$\displaystyle\frac{1}{36}\left[6\widehat{m}^{s}_{1}+6\widehat{m}^{s}_{2}+2\widehat{m}^{s}_{3}+9\widehat{m}^{s}_{5}-3\widehat{m}^{s}_{6}-3\widehat{m}^{s}_{7}+\widehat{m}^{s}_{8}\right],$
(78) $\displaystyle S_{6}$ $\displaystyle=$
$\displaystyle\frac{1}{36}\left[-6\widehat{m}^{s}_{1}+6\widehat{m}^{s}_{2}+2\widehat{m}^{s}_{3}-9\widehat{m}^{s}_{5}-3\widehat{m}^{s}_{6}+3\widehat{m}^{s}_{7}+\widehat{m}^{s}_{8}\right],$
(79) $\displaystyle S_{7}$ $\displaystyle=$
$\displaystyle\frac{1}{36}\left[-6\widehat{m}^{s}_{1}-6\widehat{m}^{s}_{2}+2\widehat{m}^{s}_{3}+9\widehat{m}^{s}_{5}+3\widehat{m}^{s}_{6}+3\widehat{m}^{s}_{7}+\widehat{m}^{s}_{8}\right],$
(80) $\displaystyle S_{8}$ $\displaystyle=$
$\displaystyle\frac{1}{36}\left[6\widehat{m}^{s}_{1}-6\widehat{m}^{s}_{2}+2\widehat{m}^{s}_{3}-9\widehat{m}^{s}_{5}+3\widehat{m}^{s}_{6}-3\widehat{m}^{s}_{7}+\widehat{m}^{s}_{8}\right].$
(81)

This is the explicit set of expressions for the source terms in velocity space
$S_{\alpha}$ given in terms of $\overrightarrow{F}$, $\overrightarrow{u}$ and
$\widehat{\overline{\kappa}}_{x^{m}y^{n}}$, with $2\leq(m+n)\leq 3$ and $0\leq
m,n\leq 2$.

Again, using the orthogonal property of $\mathcal{K}$, we can obtain the raw
moments of the collision kernel

$\sum_{\alpha}(\mathcal{K}\cdot\mathbf{\widehat{g}})_{\alpha}e_{\alpha
x}^{m}e_{\alpha y}^{n}=\sum_{\beta}\braket{K_{\beta}}{e_{\alpha
x}^{m}e_{\alpha y}^{n}}\widehat{g}_{\beta},$ (82)

which is of central importance in the subsequent derivation. Note that for
collision invariants, $\widehat{g}_{0}=\widehat{g}_{1}=\widehat{g}_{2}=0$. We
get

$\displaystyle\sum_{\alpha}(\mathcal{K}\cdot\mathbf{\widehat{g}})_{\alpha}=\sum_{\beta}\braket{K_{\beta}}{\rho}\widehat{g}_{\beta}$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\sum_{\alpha}(\mathcal{K}\cdot\mathbf{\widehat{g}})_{\alpha}e_{\alpha
x}=\sum_{\beta}\braket{K_{\beta}}{e_{\alpha x}}\widehat{g}_{\beta}$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\sum_{\alpha}(\mathcal{K}\cdot\mathbf{\widehat{g}})_{\alpha}e_{\alpha
y}=\sum_{\beta}\braket{K_{\beta}}{e_{\alpha y}}\widehat{g}_{\beta}$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\sum_{\alpha}(\mathcal{K}\cdot\mathbf{\widehat{g}})_{\alpha}e_{\alpha
x}^{2}=\sum_{\beta}\braket{K_{\beta}}{e_{\alpha x}^{2}}\widehat{g}_{\beta}$
$\displaystyle=$ $\displaystyle 6\widehat{g}_{3}+2\widehat{g}_{4},$
$\displaystyle\sum_{\alpha}(\mathcal{K}\cdot\mathbf{\widehat{g}})_{\alpha}e_{\alpha
y}^{2}=\sum_{\beta}\braket{K_{\beta}}{e_{\alpha y}^{2}}\widehat{g}_{\beta}$
$\displaystyle=$ $\displaystyle 6\widehat{g}_{3}-2\widehat{g}_{4},$ (83)
$\displaystyle\sum_{\alpha}(\mathcal{K}\cdot\mathbf{\widehat{g}})_{\alpha}e_{\alpha
x}e_{\alpha y}=\sum_{\beta}\braket{K_{\beta}}{e_{\alpha x}e_{\alpha
y}}\widehat{g}_{\beta}$ $\displaystyle=$ $\displaystyle 4\widehat{g}_{5},$
$\displaystyle\sum_{\alpha}(\mathcal{K}\cdot\mathbf{\widehat{g}})_{\alpha}e_{\alpha
x}^{2}e_{\alpha y}=\sum_{\beta}\braket{K_{\beta}}{e_{\alpha x}^{2}e_{\alpha
y}}\widehat{g}_{\beta}$ $\displaystyle=$ $\displaystyle-4\widehat{g}_{6},$
$\displaystyle\sum_{\alpha}(\mathcal{K}\cdot\mathbf{\widehat{g}})_{\alpha}e_{\alpha
x}e_{\alpha y}^{2}=\sum_{\beta}\braket{K_{\beta}}{e_{\alpha x}e_{\alpha
y}^{2}}\widehat{g}_{\beta}$ $\displaystyle=$ $\displaystyle-4\widehat{g}_{7},$
$\displaystyle\sum_{\alpha}(\mathcal{K}\cdot\mathbf{\widehat{g}})_{\alpha}e_{\alpha
x}^{2}e_{\alpha y}^{2}=\sum_{\beta}\braket{K_{\beta}}{e_{\alpha
x}^{2}e_{\alpha y}^{2}}\widehat{g}_{\beta}$ $\displaystyle=$ $\displaystyle
8\widehat{g}_{3}+4\widehat{g}_{8}.$

Finally, the LBE in Eq. (22) can be rewritten in terms of collision and
streaming steps, respectively, as

$\displaystyle\widetilde{\overline{f}}_{\alpha}(\overrightarrow{x},t)$
$\displaystyle=$
$\displaystyle\overline{f}_{\alpha}(\overrightarrow{x},t)+\Omega_{{\alpha}(\overrightarrow{x},t)}^{c}+S_{{\alpha}(\overrightarrow{x},t)},$
(84)
$\displaystyle\overline{f}_{\alpha}(\overrightarrow{x}+\overrightarrow{e}_{\alpha},t+1)$
$\displaystyle=$
$\displaystyle\widetilde{\overline{f}}_{\alpha}(\overrightarrow{x},t),$ (85)

where the symbol “tilde” ($\sim$) in the first equation refers to the post-
collision state. In terms of the transformed distribution, the hydrodynamic
fields can be computed by means of the following:

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\sum_{\alpha=0}^{8}\overline{f}_{\alpha}=\braket{\overline{f}_{\alpha}}{\rho},$
(86) $\displaystyle\rho u_{i}$ $\displaystyle=$
$\displaystyle\sum_{\alpha=0}^{8}\overline{f}_{\alpha}e_{\alpha
i}+\frac{1}{2}F_{i}=\braket{\overline{f}_{\alpha}}{e_{\alpha
i}}+\frac{1}{2}F_{i},\qquad i\in{x,y}.$ (87)

##  10 Structure of the Collision Operator: Cascaded Central Moments

Let us now arrive at the expressions for the cascaded formulation of the
collision operator using central moments in the presence of forcing terms
based on the results obtained in the last few sections. The basic procedure
can be stated as follows. Beginning from the lowest order central moments that
are non-collisional invariants (i.e. $\widehat{\overline{\kappa}}_{xx}$ and
higher), they are successively set equal to their local attractors based on
the transformed Maxwellians (Eq. (34)). This step provides tentative
expressions for $\widehat{g}_{\alpha}$ based on the equilibrium assumption. We
then modify them to allow for relaxation during collision by multiplying them
with corresponding relaxation parameters [59]. In this step, given the
cascaded nature of the collision (i.e.
$\widehat{g}_{\alpha}\equiv\widehat{g}_{\alpha}(\mathbf{f},\widehat{g}_{\beta}),\beta=0,1,2,\ldots,\alpha-1$,
or the dependence of higher order terms on those that are lower orders), care
needs to be exercised to multiply the relaxation parameters only with those
terms that are not yet in post-collision states (i.e. terms not involving
$\widehat{g}_{\beta},\beta=0,1,2,\ldots,\alpha-1$ for $\widehat{g}_{\alpha}$).
Various details involved in this procedure are given in [63]. For brevity,
here we summarize the final results which are as follows:

$\displaystyle\widehat{g}_{3}$ $\displaystyle=$
$\displaystyle\frac{\omega_{3}}{12}\left\\{\frac{2}{3}\rho+\rho(u_{x}^{2}+u_{y}^{2})-(\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}+\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}})-\frac{1}{2}(\widehat{\sigma}_{xx}^{{}^{\prime}}+\widehat{\sigma}_{yy}^{{}^{\prime}})\right\\},$
(88) $\displaystyle\widehat{g}_{4}$ $\displaystyle=$
$\displaystyle\frac{\omega_{4}}{4}\left\\{\rho(u_{x}^{2}-u_{y}^{2})-(\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}-\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}})-\frac{1}{2}(\widehat{\sigma}_{xx}^{{}^{\prime}}-\widehat{\sigma}_{yy}^{{}^{\prime}})\right\\},$
(89) $\displaystyle\widehat{g}_{5}$ $\displaystyle=$
$\displaystyle\frac{\omega_{5}}{4}\left\\{\rho
u_{x}u_{y}-\widehat{\overline{\kappa}}_{xy}^{{}^{\prime}}-\frac{1}{2}\widehat{\sigma}_{xy}^{{}^{\prime}}\right\\},$
(90) $\displaystyle\widehat{g}_{6}$ $\displaystyle=$
$\displaystyle\frac{\omega_{6}}{4}\left\\{2\rho
u_{x}^{2}u_{y}+\widehat{\overline{\kappa}}_{xxy}^{{}^{\prime}}-2u_{x}\widehat{\overline{\kappa}}_{xy}^{{}^{\prime}}-u_{y}\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}-\frac{1}{2}\widehat{\sigma}_{xxy}\right\\}-\frac{1}{2}u_{y}(3\widehat{g}_{3}+\widehat{g}_{4})$
(91) $\displaystyle-2u_{x}\widehat{g}_{5},$ $\displaystyle\widehat{g}_{7}$
$\displaystyle=$ $\displaystyle\frac{\omega_{7}}{4}\left\\{2\rho
u_{x}u_{y}^{2}+\widehat{\overline{\kappa}}_{xyy}^{{}^{\prime}}-2u_{y}\widehat{\overline{\kappa}}_{xy}^{{}^{\prime}}-u_{x}\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}}-\frac{1}{2}\widehat{\sigma}_{xyy}\right\\}-\frac{1}{2}u_{x}(3\widehat{g}_{3}-\widehat{g}_{4})$
(92) $\displaystyle-2u_{y}\widehat{g}_{5},$ $\displaystyle\widehat{g}_{8}$
$\displaystyle=$
$\displaystyle\frac{\omega_{8}}{4}\left\\{\frac{1}{9}\rho+3\rho
u_{x}^{2}u_{y}^{2}-\left[\widehat{\overline{\kappa}}_{xxyy}^{{}^{\prime}}-2u_{x}\widehat{\overline{\kappa}}_{xyy}^{{}^{\prime}}-2u_{y}\widehat{\overline{\kappa}}_{xxy}^{{}^{\prime}}+u_{x}^{2}\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}}+u_{y}^{2}\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}\right.\right.$
(93)
$\displaystyle\left.\left.+4u_{x}u_{y}\widehat{\overline{\kappa}}_{xy}^{{}^{\prime}}\right]-\frac{1}{2}\widehat{\sigma}_{xxyy}^{{}^{\prime}}\right\\}-2\widehat{g}_{3}-\frac{1}{2}u_{y}^{2}(3\widehat{g}_{3}+\widehat{g}_{4})-\frac{1}{2}u_{x}^{2}(3\widehat{g}_{3}-\widehat{g}_{4})$
$\displaystyle-4u_{x}u_{y}\widehat{g}_{5}-2u_{y}\widehat{g}_{6}-2u_{x}\widehat{g}_{7}.$

In the above, the raw moments of various orders, i.e.
$\widehat{\overline{\kappa}}_{x^{m}y^{n}}^{{}^{\prime}}$ for different $m$ and
$n$ are required, which may be obtained from Eqs. (40)-(48). Similarly, the
raw moments of sources of various orders, i.e.
$\widehat{\sigma}_{x^{m}y^{n}}^{{}^{\prime}}$ needed in the above are given in
Eqs. (49)-(57) (see Sec. 9). Here, $\omega_{\beta}$, where
$\beta=3,4,5,\ldots,8$, are the relaxation parameters, satisfying
$0<\omega_{\beta}<2$. When a multiscale Chapman-Enskog expansion [64] is
applied to this central moment LBM based on central moments, it recovers the
Navier-Stokes equation with the relaxation parameters
$\omega_{3}=\omega^{\chi}$ and $\omega_{4}=\omega_{5}=\omega^{\nu}$
controlling bulk and shear viscosities, respectively (e.g.,
$\nu=c_{s}^{2}\left(\frac{1}{\omega^{\nu}}-\frac{1}{2}\right)$) [63]. The rest
of the parameters can be adjusted to improve numerical stability.

##  11 Cascaded Central Moment Lattice Boltzmann Equation

The collision step with the addition of forcing terms (see Eq. (19) and Eq.
(84)) in the stream-and-collide procedure of the LBM, obtained by matching
those of the continuous Boltzmann equation as discussed in the previous
sections, is expanded element-wise and can be summarized as follows:

$\displaystyle\widetilde{\overline{f}}_{0}$ $\displaystyle=$
$\displaystyle\overline{f}_{0}+\left[\widehat{g}_{0}-4(\widehat{g}_{3}-\widehat{g}_{8})\right]+S_{0},$
$\displaystyle\widetilde{\overline{f}}_{1}$ $\displaystyle=$
$\displaystyle\overline{f}_{1}+\left[\widehat{g}_{0}+\widehat{g}_{1}-\widehat{g}_{3}+\widehat{g}_{4}+2(\widehat{g}_{7}-\widehat{g}_{8})\right]+S_{1},$
$\displaystyle\widetilde{\overline{f}}_{2}$ $\displaystyle=$
$\displaystyle\overline{f}_{2}+\left[\widehat{g}_{0}+\widehat{g}_{2}-\widehat{g}_{3}-\widehat{g}_{4}+2(\widehat{g}_{6}-\widehat{g}_{8})\right]+S_{2},$
$\displaystyle\widetilde{\overline{f}}_{3}$ $\displaystyle=$
$\displaystyle\overline{f}_{3}+\left[\widehat{g}_{0}-\widehat{g}_{1}-\widehat{g}_{3}+\widehat{g}_{4}-2(\widehat{g}_{7}+\widehat{g}_{8})\right]+S_{3},$
$\displaystyle\widetilde{\overline{f}}_{4}$ $\displaystyle=$
$\displaystyle\overline{f}_{4}+\left[\widehat{g}_{0}-\widehat{g}_{2}-\widehat{g}_{3}-\widehat{g}_{4}-2(\widehat{g}_{6}+\widehat{g}_{8})\right]+S_{4},$
(94) $\displaystyle\widetilde{\overline{f}}_{5}$ $\displaystyle=$
$\displaystyle\overline{f}_{5}+\left[\widehat{g}_{0}+\widehat{g}_{1}+\widehat{g}_{2}+2\widehat{g}_{3}+\widehat{g}_{5}-\widehat{g}_{6}-\widehat{g}_{7}+\widehat{g}_{8}\right]+S_{5},$
$\displaystyle\widetilde{\overline{f}}_{6}$ $\displaystyle=$
$\displaystyle\overline{f}_{6}+\left[\widehat{g}_{0}-\widehat{g}_{1}+\widehat{g}_{2}+2\widehat{g}_{3}-\widehat{g}_{5}-\widehat{g}_{6}+\widehat{g}_{7}+\widehat{g}_{8}\right]+S_{6},$
$\displaystyle\widetilde{\overline{f}}_{7}$ $\displaystyle=$
$\displaystyle\overline{f}_{7}+\left[\widehat{g}_{0}-\widehat{g}_{1}-\widehat{g}_{2}+2\widehat{g}_{3}+\widehat{g}_{5}+\widehat{g}_{6}+\widehat{g}_{7}+\widehat{g}_{8}\right]+S_{7},$
$\displaystyle\widetilde{\overline{f}}_{8}$ $\displaystyle=$
$\displaystyle\overline{f}_{8}+\left[\widehat{g}_{0}+\widehat{g}_{1}-\widehat{g}_{2}+2\widehat{g}_{3}-\widehat{g}_{5}+\widehat{g}_{6}-\widehat{g}_{7}+\widehat{g}_{8}\right]+S_{8}.$

The collision kernel $\widehat{g}_{\beta}$ needed here can be computed from
the expressions given in the previous section (see Sec. 10). The source terms
in the velocity space can be obtained from Eqs. (73)-(81) (see Sec. 9). The
remaining streaming step is carried out as usual by using the post collision
values $\widetilde{\overline{f}}_{\alpha}$ obtained from above. Once the local
distribution function is known, macroscopic fluid density and velocity fields
satisfying the Galilean invariant Navier-Stokes equations in the presence of
force fields can be obtained from Eqs. (86) and (87), respectively.

##  12 Cascaded Collision Operator with Reduced Compressibility Effects

While a main goal of this work is the introduction of a self-consistent
approach based on the continuous Boltzmann equation to incorporate non-
equilibrium effects into the central moment approach for general
applicability, it is also useful to consider its limiting cases. For example,
the incompressible limit of fluid flow corresponds to considering very small
deviations from the local equilibrium, a special case with various
applications. In particular, this would allow simple representation of
incompressible turbulence considered later in this work.

Being a kinetic approach, the lattice Boltzmann method is inherently
compressible in nature. On the other hand, when it is desired to reproduce the
“incompressible” Navier-Stokes equations as mentioned above, it is important
to reduce such compressibility effects. An approach in this regard was
introduced earlier in [87]. Here, we will extend it further in the context of
the central moment LBM in the presence of forcing. It may be noted that the
fundamental expressions for the continuous central moments for the local
equilibrium as well as the forcing given in Secs. 4 and 5, respectively, from
which their discrete counterparts are derived, remains unchanged for this
case. However, the key element to incorporate a systematic reduction of
compressibility effects lies in the following careful definition of the raw
moments of the hydrodynamic fields:

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\sum_{\alpha=0}^{8}f_{\alpha}=\sum_{\alpha=0}^{8}\overline{f}_{\alpha},$
(95) $\displaystyle\rho_{0}\overrightarrow{u}$ $\displaystyle=$
$\displaystyle\sum_{\alpha=0}^{8}f_{\alpha}\overrightarrow{e}_{\alpha}=\sum_{\alpha=0}^{8}\overline{f}_{\alpha}\overrightarrow{e}_{\alpha}+\frac{1}{2}\overrightarrow{F}.$
(96)

where $\rho=\rho_{0}+\delta\rho$. Here $\rho_{0}$ and $\delta\rho$ are the
constant reference value and fluctuations of density, respectively. That is,
in the above, contributions of density fluctuations are eliminated from first-
order moments representing the components of momentum. Thus, we get

$\displaystyle\braket{\widetilde{\overline{f}}_{\alpha}}{\rho}$
$\displaystyle=$ $\displaystyle\rho,$ (97)
$\displaystyle\braket{\widetilde{\overline{f}}_{\alpha}}{e_{\alpha x}}$
$\displaystyle=$ $\displaystyle\rho_{0}u_{x}+\frac{1}{2}F_{x},$ (98)
$\displaystyle\braket{\widetilde{\overline{f}}_{\alpha}}{e_{\alpha y}}$
$\displaystyle=$ $\displaystyle\rho_{0}u_{y}+\frac{1}{2}F_{y}.$ (99)

Using the procedure discussed in the previous sections and with the above
specialized re-definition of the conserved moments, we obtain, after some
simplification, the cascaded collision operator with reduced compressibility
effects. They are reported here in the following:

$\displaystyle\widehat{g}_{3}$ $\displaystyle=$
$\displaystyle\frac{\omega_{3}}{12}\left\\{\frac{2}{3}\rho+(\rho_{0}-\delta\rho)(u_{x}^{2}+u_{y}^{2})-(\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}+\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}})-\frac{1}{2}(\widehat{\sigma}_{xx}^{{}^{\prime}}+\widehat{\sigma}_{yy}^{{}^{\prime}})\right\\},$
(100) $\displaystyle\widehat{g}_{4}$ $\displaystyle=$
$\displaystyle\frac{\omega_{4}}{4}\left\\{(\rho_{0}-\delta\rho)(u_{x}^{2}-u_{y}^{2})-(\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}-\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}})-\frac{1}{2}(\widehat{\sigma}_{xx}^{{}^{\prime}}-\widehat{\sigma}_{yy}^{{}^{\prime}})\right\\},$
(101) $\displaystyle\widehat{g}_{5}$ $\displaystyle=$
$\displaystyle\frac{\omega_{5}}{4}\left\\{(\rho_{0}-\delta\rho)u_{x}u_{y}-\widehat{\overline{\kappa}}_{xy}^{{}^{\prime}}-\frac{1}{2}\widehat{\sigma}_{xy}^{{}^{\prime}}\right\\},$
(102) $\displaystyle\widehat{g}_{6}$ $\displaystyle=$
$\displaystyle\frac{\omega_{6}}{4}\left\\{(2\rho_{0}-\delta\rho)u_{x}^{2}u_{y}+\widehat{\overline{\kappa}}_{xxy}^{{}^{\prime}}-2u_{x}\widehat{\overline{\kappa}}_{xy}^{{}^{\prime}}-u_{y}\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}-\frac{1}{2}\widehat{\sigma}_{xxy}\right\\}-\frac{1}{2}u_{y}(3\widehat{g}_{3}+\widehat{g}_{4})$
(103) $\displaystyle-2u_{x}\widehat{g}_{5},$ $\displaystyle\widehat{g}_{7}$
$\displaystyle=$
$\displaystyle\frac{\omega_{7}}{4}\left\\{(2\rho_{0}-\delta\rho)u_{x}u_{y}^{2}+\widehat{\overline{\kappa}}_{xyy}^{{}^{\prime}}-2u_{y}\widehat{\overline{\kappa}}_{xy}^{{}^{\prime}}-u_{x}\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}}-\frac{1}{2}\widehat{\sigma}_{xyy}\right\\}-\frac{1}{2}u_{x}(3\widehat{g}_{3}-\widehat{g}_{4})$
(104) $\displaystyle-2u_{y}\widehat{g}_{5},$ $\displaystyle\widehat{g}_{8}$
$\displaystyle=$
$\displaystyle\frac{\omega_{8}}{4}\left\\{\frac{1}{9}\rho+(3\rho_{0}-\delta\rho)u_{x}^{2}u_{y}^{2}-\left[\widehat{\overline{\kappa}}_{xxyy}^{{}^{\prime}}-2u_{x}\widehat{\overline{\kappa}}_{xyy}^{{}^{\prime}}-2u_{y}\widehat{\overline{\kappa}}_{xxy}^{{}^{\prime}}+u_{x}^{2}\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}}+u_{y}^{2}\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}\right.\right.$
(105)
$\displaystyle\left.\left.+4u_{x}u_{y}\widehat{\overline{\kappa}}_{xy}^{{}^{\prime}}\right]-\frac{1}{2}\widehat{\sigma}_{xxyy}^{{}^{\prime}}\right\\}-2\widehat{g}_{3}-\frac{1}{2}u_{y}^{2}(3\widehat{g}_{3}+\widehat{g}_{4})-\frac{1}{2}u_{x}^{2}(3\widehat{g}_{3}-\widehat{g}_{4})$
$\displaystyle-4u_{x}u_{y}\widehat{g}_{5}-2u_{y}\widehat{g}_{6}-2u_{x}\widehat{g}_{7}.$

The above collision operator that selectively introduces density fluctuations
where necessary can reduce compressibility effects for a inertial frame
invariant flow field while retaining its linear acoustics. It can thus allow
for a better comparison with “incompressible” macroscopic fluid dynamic
equations, particularly for turbulent flows as discussed later.

##  13 Factorized Central Moment Model for Collision

In this section, we will derive an alternative form of the central moment LBE
with forcing terms based on a different choice of the local attractor in the
collision step for improved physical coherence. Continuous kinetic theory, as
originally initiated by Maxwell [88], features two important properties for
the local equilibrium or the Maxwell distribution – Galilean invariance and
factorization in Cartesian components of the particle velocity. As discussed
recently [62, 39], it could prove useful to inherent these properties at the
discrete particle velocity level. The use of central moments maintains
Galilean invariance by construction. Factorization property of the
distribution function implies that particle velocities are random variables.
An extension of the factorization idea beyond equilibrium as a model for
describing the relaxation process during collision was proposed to construct
local attractors [62]. Specifically, the basic postulate behind this model is
that the Cartesian products of the post-collision values of the orthogonal
central moments of lower orders that are not in equilibrium forms as the basis
for the attractors of the higher order moments. Here, we further extend this
to include source terms so that the model can incorporate force fields. Thus,
the attractors for central moments of different orders are given as

$\displaystyle\widehat{\kappa}_{x}^{at}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\kappa}}_{x}=0,$ (106)
$\displaystyle\widehat{\kappa}_{y}^{at}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\kappa}}_{y}=0,$ (107)
$\displaystyle\widehat{\kappa}_{xy}^{at}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\kappa}}_{x}\widetilde{\widehat{\kappa}}_{y}=0,$
(108) $\displaystyle\widehat{\kappa}_{xxy}^{at}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\kappa}}_{xx}\widetilde{\widehat{\kappa}}_{y}=0,$
(109) $\displaystyle\widehat{\kappa}_{xyy}^{at}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\kappa}}_{x}\widetilde{\widehat{\kappa}}_{yy}=0,$
(110) $\displaystyle\widehat{\kappa}_{xxyy}^{at}$ $\displaystyle=$
$\displaystyle\widetilde{\widehat{\kappa}}_{xx}\widetilde{\widehat{\kappa}}_{yy},$
(111)

while the second-order longitudinal central moments are obtained from
Maxwellian as given in an earlier section (see Sec. 7), i.e.
$\widehat{\kappa}_{xx}^{at}=\widehat{\kappa}_{yy}^{at}=\rho c_{s}^{2}$. In
essence, the distinguishing feature of the factorized central moment lies in
the use of modified attractors for third and higher order moments. Now, using
the following central moment identity of the post-collision state
$\widetilde{\widehat{\kappa}}_{x^{m}y^{n}}=\widehat{\kappa}_{x^{m}y^{n}}+\sum_{\beta}\braket{K_{\alpha}}{(e_{\alpha
x}-u_{x})^{m}(e_{\alpha y}-u_{y})^{n}}\widehat{g}_{\beta}$, for $m=2,n=0$ and
$m=0,n=2$, we get

$\displaystyle\widetilde{\widehat{\kappa}}_{xx}$ $\displaystyle=$
$\displaystyle\widehat{\kappa}_{xx}+(6\widehat{g}_{3}+2\widehat{g}_{4}),$
(112) $\displaystyle\widetilde{\widehat{\kappa}}_{yy}$ $\displaystyle=$
$\displaystyle\widehat{\kappa}_{yy}+(6\widehat{g}_{3}-2\widehat{g}_{4}).$
(113)

Note that it also follows that
$\widehat{\kappa}_{xx}=\widehat{\overline{\kappa}}_{xx}$ and
$\widehat{\kappa}_{yy}=\widehat{\overline{\kappa}}_{yy}$. We can then rewrite
everything in terms of transformed raw moments, i.e.
$\widehat{\overline{\kappa}}_{xx}=\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}-\rho
u_{x}^{2}+F_{x}u_{x}$ and
$\widehat{\overline{\kappa}}_{yy}=\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}}-\rho
u_{y}^{2}+F_{y}u_{y}$. These yield

$\displaystyle\widetilde{\widehat{\kappa}}_{xx}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}-\rho
u_{x}^{2}+F_{x}u_{x}+(6\widehat{g}_{3}+2\widehat{g}_{4}),$ (114)
$\displaystyle\widetilde{\widehat{\kappa}}_{yy}$ $\displaystyle=$
$\displaystyle\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}}-\rho
u_{y}^{2}+F_{y}u_{y}+(6\widehat{g}_{3}-2\widehat{g}_{4}).$ (115)

In effect, the attractor for the fourth-order moment, i.e. Eq. (111) reduces
to

$\displaystyle\widehat{\kappa}_{xxyy}^{at}$ $\displaystyle=$
$\displaystyle\left[\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}-\rho
u_{x}^{2}+F_{x}u_{x}+(6\widehat{g}_{3}+2\widehat{g}_{4})\right]\times$ (116)
$\displaystyle\left[\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}}-\rho
u_{y}^{2}+F_{y}u_{y}+(6\widehat{g}_{3}-2\widehat{g}_{4})\right].$

Now, to obtain an operational step in terms of the transformed variables, we
use the relation
$\widehat{\overline{\kappa}}_{x^{m}y^{n}}^{at}=\widehat{\kappa}_{x^{m}y^{n}}^{at}-\frac{1}{2}\widehat{\sigma}_{x^{m}y^{n}}$
to finally get the following expression for the fourth-order central moment

$\displaystyle\widehat{\overline{\kappa}}_{xxyy}^{at}$ $\displaystyle=$
$\displaystyle\left[\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}-\rho
u_{x}^{2}+F_{x}u_{x}+(6\widehat{g}_{3}+2\widehat{g}_{4})\right]\left[\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}}-\rho
u_{y}^{2}+F_{y}u_{y}+(6\widehat{g}_{3}-2\widehat{g}_{4})\right]$ (117)
$\displaystyle-\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{xyy}-\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{xxy}-\frac{F_{y}^{2}}{2\rho^{2}}\widehat{\overline{\kappa}}_{xx}-\frac{F_{x}^{2}}{2\rho^{2}}\widehat{\overline{\kappa}}_{yy}-2\frac{F_{x}F_{y}}{\rho^{2}}\widehat{\overline{\kappa}}_{xy}.$

By replacing $\widehat{\overline{\kappa}}_{x^{m}y^{n}}^{\mathcal{M}}$ with
$\widehat{\overline{\kappa}}_{x^{m}y^{n}}^{at}$ as given above, we can derive
the collision kernel with factorized central moments as attractors. It follows
that the expressions in Sec. 10 for $\widehat{g}_{\beta}$, $\beta=3,4,5,6,7$
are the same as before with the exception for $\widehat{g}_{8}$. The
expression for $\widehat{g}_{8}$ in Eq. (93) is modified such that term
$\frac{1}{9}\rho(=\widehat{\overline{\kappa}}_{xxyy}^{\mathcal{M}})$ in this
equation is now replaced by $\widehat{\overline{\kappa}}_{xxyy}^{at}$ given in
Eq. (117). In a similar vein, the above expression can be modified for reduced
compressibility effects (see Sec. 12) as

$\displaystyle\widehat{\kappa}_{xxyy}^{at}$ $\displaystyle=$
$\displaystyle\left[\widehat{\overline{\kappa}}_{xx}^{{}^{\prime}}-(\rho_{0}-\delta\rho)u_{x}^{2}+F_{x}u_{x}+(6\widehat{g}_{3}+2\widehat{g}_{4})\right]\times$
(118)
$\displaystyle\left[\widehat{\overline{\kappa}}_{yy}^{{}^{\prime}}-(\rho_{0}-\delta\rho)u_{y}^{2}+F_{y}u_{y}+(6\widehat{g}_{3}-2\widehat{g}_{4})\right]$
$\displaystyle-\frac{F_{x}}{\rho}\widehat{\overline{\kappa}}_{xyy}-\frac{F_{y}}{\rho}\widehat{\overline{\kappa}}_{xxy}-\frac{F_{y}^{2}}{2\rho^{2}}\widehat{\overline{\kappa}}_{xx}-\frac{F_{x}^{2}}{2\rho^{2}}\widehat{\overline{\kappa}}_{yy}-2\frac{F_{x}F_{y}}{\rho^{2}}\widehat{\overline{\kappa}}_{xy}$

to modify $\widehat{g}_{8}$ in Eq. (105).

##  14 Galilean Invariant Filtered Turbulence Representation using Lattice
Kinetic Framework

Based on the various elements derived in the previous sections, we are now in
a position to construct an approach for simulation of Galilean invariant
turbulent flow field by incorporating appropriate turbulence models in the
LBM. The starting point in the statistical continuum description of turbulence
is the Reynolds decomposition of the velocity field of the fluid into
‘resolved’ and ‘unresolved’ parts. The resolved part is obtained by applying
either some averaging in space or time (in the Reynolds Averaged Navier-Stokes
(RANS) context) or by applying a filter (in the LES). Application of this
decomposition to the Navier-Stokes (NS) equation leads to additional unknown
terms involving products of the unresolved fields, which are Reynolds stresses
(in RANS) or the subgrid stresses (in LES). This closure problem then becomes
the main focus of turbulence modeling. Due to the scale invariance property of
the NS equations [83], the averaged and the filtered equations, as well as the
additional stress-like closure terms have similar forms. Thus, a unified
statistical approach may be adopted for turbulence modeling. It is interesting
to note that ideas based on kinetic theory provided the original inspiration
for the Reynolds decomposition [89] as well as early works on developing
turbulence models.

The underlying motivation here is to develop a unified statistical averaged
description (for RANS) or formal spatial filtered representation (for LES) of
_inertial frame invariant_ turbulence in a kinetic approach based on the LBM
derived in earlier sections. This would also allow reconciliation of continuum
and non-continuum effects on turbulence as discussed recently [77]. The
following notation for Reynolds decomposition is adopted here. For any scalar
$\phi$, vector $\overrightarrow{v}$ and tensor $T_{ij}$, we have

$\displaystyle\phi$ $\displaystyle=$
$\displaystyle\underline{\phi}+\phi^{{}^{\prime}},\quad\mbox{with}\quad\underline{\phi^{{}^{\prime}}}=0,$
$\displaystyle\overrightarrow{v}$ $\displaystyle=$
$\displaystyle\underline{\overrightarrow{v}}+\overrightarrow{v}^{{}^{\prime}},\quad\mbox{with}\quad\underline{\overrightarrow{v}^{{}^{\prime}}}=0,$
$\displaystyle T_{ij}$ $\displaystyle=$
$\displaystyle\underline{T_{ij}}+T_{ij}^{{}^{\prime}},\quad\mbox{with}\quad\underline{T_{ij}^{{}^{\prime}}}=0,$

where $\underline{(\cdot)}$ is an operator representing either some form of
statistical average or filter to obtain the resolved part and the symbols with
primes denote the unresolved parts of the field. As discussed in [77],
application of the above decomposition directly to the continuous Boltzmann
equation (Eq. (6)) leads to certain difficulties. In particular, using
$f=\underline{f}+f^{{}^{\prime}}$ for the distribution function in Eq. (6),
which leads to a statistically averaged kinetic equation, does not provide a
clear interpretation of turbulence physics. The local collision term needs to
model all essential physics, including the non-linear and non-local momentum
transfer effects of turbulence. Moreover, the use of the averaged attractor
based on the Maxwellian $\underline{f^{\mathcal{M}}}$ within the collision
term $\underline{\Omega(f,f)}$ leads to modeling difficulties since
$\underline{\exp\left[-\frac{(\overrightarrow{\xi}-\overrightarrow{u})^{2}}{2c_{s}^{2}}\right]}\neq\exp\left[-\frac{(\overrightarrow{\xi}-\underline{\overrightarrow{u}})^{2}}{2c_{s}^{2}}\right]$.
Thus, an alternative approach is needed to coherently represent
continuum/kinetic effects on turbulence.

To circumvent these issues, a transformation for the velocities was recently
suggested [77], which is adopted here to provide Galilean invariant turbulence
representation. The key element is to clearly separate the advective
turbulence effects due to unresolved velocity field
$\overrightarrow{u}^{{}^{\prime}}$ from the dissipative collision that
represent microscopic effects. This is accomplished by an inspection of the
local Maxwellian given in terms of the microscopic particle velocity
$\overrightarrow{\xi}$ and the macroscopic fluid velocity
$\overrightarrow{u}$. That is, it consists of the term involving the peculiar
velocity

$\overrightarrow{\xi}-\overrightarrow{u}$

as its argument which should be made independent of the unresolved part of the
macroscopic fluid velocity $\overrightarrow{u}^{{}^{\prime}}$, when the
averaging operator is applied. That is,

$\overrightarrow{\xi}-(\underline{\overrightarrow{u}}+\overrightarrow{u}^{{}^{\prime}})=(\overrightarrow{\xi}-\overrightarrow{u}^{{}^{\prime}})-\underline{\overrightarrow{u}}$

should be transformed appropriately, which can be accomplished by defining a
new variable $\overrightarrow{\eta}$ as

$\overrightarrow{\eta}=\overrightarrow{\xi}-\overrightarrow{u}^{{}^{\prime}}.$
(119)

Now, the Maxwellian in the transformed peculiar velocity
$\overrightarrow{\eta}-\underline{\overrightarrow{u}}$ commutes with the
operator for averaging or filtering. That is,
$\underline{\exp\left[-\frac{(\overrightarrow{\eta}-\overrightarrow{u})^{2}}{2c_{s}^{2}}\right]}=\exp\left[-\frac{(\overrightarrow{\eta}-\underline{\overrightarrow{u}})^{2}}{2c_{s}^{2}}\right].$
This facilitates the separation of various aspects of turbulence physics
modeling. Based on this a new distribution function
$h(\overrightarrow{x},\overrightarrow{\eta},t)$ and its local Maxwellian are
defined by

$h(\overrightarrow{x},\overrightarrow{\eta},t)=f(\overrightarrow{x},\overrightarrow{\xi},t),$
(120)

and

$h^{\mathcal{M}}(\overrightarrow{\eta},\overrightarrow{u})=\frac{\rho}{2\pi
c_{s}^{2}}\exp\left[-\frac{(\overrightarrow{\eta}-\underline{\overrightarrow{u}})^{2}}{2c_{s}^{2}}\right],$
(121)

respectively. The continuous Boltzmann equation, i.e. Eq. (6) (without the
forcing term for simplicity) is then transformed into a modified kinetic
equation in terms of $\eta$ and $h$ as follows. From Eq. (120),
$\overrightarrow{\nabla}_{\eta}h=\overrightarrow{\nabla}_{\xi}f$. When
$\eta=\mbox{constant}$, we have
$(\overrightarrow{\nabla}_{x}\overrightarrow{\xi})_{\eta}=\overrightarrow{\nabla}_{x}\overrightarrow{u}^{{}^{\prime}}$
and
$(\partial_{t}\overrightarrow{\xi})_{\eta}=\partial_{t}\overrightarrow{u}^{{}^{\prime}}$.
Hence, the derivatives in new variables are

$\displaystyle\overrightarrow{\nabla}_{x}h$ $\displaystyle=$
$\displaystyle\overrightarrow{\nabla}_{x}f+(\overrightarrow{\nabla}_{x}\overrightarrow{\xi})_{\eta}\cdot\overrightarrow{\nabla}_{\xi}f=\overrightarrow{\nabla}_{x}f+(\overrightarrow{\nabla}_{x}\overrightarrow{u}^{{}^{\prime}})_{\eta}\cdot\overrightarrow{\nabla}_{\eta}h,$
$\displaystyle\partial_{t}h$ $\displaystyle=$
$\displaystyle\partial_{t}f+(\partial_{t}\overrightarrow{\xi})_{\eta}\cdot\overrightarrow{\nabla}_{\xi}f=\partial_{t}f+(\partial_{t}\overrightarrow{u}^{{}^{\prime}})_{\eta}\cdot\overrightarrow{\nabla}_{\eta}h.$

The continuous Boltzmann equation is thus modified to [77]

$\partial_{t}h+\overrightarrow{\eta}\cdot\overrightarrow{\nabla}_{x}h+\overrightarrow{u}^{{}^{\prime}}\cdot\overrightarrow{\nabla}_{x}h-\overrightarrow{a}^{{}^{\prime}}\cdot\overrightarrow{\nabla}_{\eta}h=\Omega(h,h),$
(122)

where
$\overrightarrow{a}^{{}^{\prime}}=\partial_{t}\overrightarrow{u}^{{}^{\prime}}+\eta\cdot\overrightarrow{\nabla}_{x}\overrightarrow{u}^{{}^{\prime}}+\overrightarrow{u}^{{}^{\prime}}\cdot\overrightarrow{\nabla}_{x}\overrightarrow{u}^{{}^{\prime}}$.
Considering incompressible flows, where the unresolved velocity field
satisfies $\overrightarrow{\nabla}\cdot\overrightarrow{u}^{{}^{\prime}}=0$, we
get
$\overrightarrow{u}^{{}^{\prime}}\cdot\overrightarrow{\nabla}_{x}h=\overrightarrow{\nabla}_{x}(h\overrightarrow{u}^{{}^{\prime}})$
and
$\overrightarrow{a}^{{}^{\prime}}\cdot\overrightarrow{\nabla}_{x}h=\overrightarrow{\nabla}_{x}(h\overrightarrow{a}^{{}^{\prime}})$.
As a result, Eq. (122) is further simplified to

$\partial_{t}h+\overrightarrow{\eta}\cdot\overrightarrow{\nabla}_{x}h+\overrightarrow{\nabla}_{x}(h\overrightarrow{u}^{{}^{\prime}})-\overrightarrow{\nabla}_{x}(h\overrightarrow{a}^{{}^{\prime}})=\Omega(h,h).$
(123)

Now, applying the statistical averaging or filtering operator on Eq. (123) and
using $\underline{\Omega(h,h)}=\Omega(\underline{h},\underline{h})$, we get
the new kinetic equation

$\partial_{t}\underline{h}+\overrightarrow{\eta}\cdot\overrightarrow{\nabla}_{x}\underline{h}+\overrightarrow{\nabla}_{x}(\underline{h\overrightarrow{u}^{{}^{\prime}}})-\overrightarrow{\nabla}_{x}(\underline{h\overrightarrow{a}^{{}^{\prime}}})=\Omega(\underline{h},\underline{h}).$
(124)

The averaged density and momentum can then be obtained by taking moments of
$\overline{h}$. That is,
$\underline{\rho}=\int\underline{h}d\overrightarrow{\eta},\underline{\rho\overrightarrow{u}}=\int\underline{h}\overrightarrow{\eta}d\overrightarrow{\eta}.$
Now, $\underline{h(\eta)}$ using Eq. (124) is better suited to represent
turbulence physics for the following reasons. The term involving
$\underline{h\overrightarrow{u}^{{}^{\prime}}}$ represents transport in
physical space, i.e. redistributes $\underline{h}$ to smoothen any gradients.
$\underline{h\overrightarrow{a}^{{}^{\prime}}}$ represents transport in
velocity phase space, and acts as a source/sink for energy cascade [77]. In
particular, both these quantities can be directly related to continuum based
closure models. On the other hand, the role of collision operator is then to
simply represent averaged effect of irreversible molecular collisions. The
averaged kinetic equation, Eq. (124), can be further simplified by considering
the following simple microscopic closure [77]

$\displaystyle\underline{h\overrightarrow{u}^{{}^{\prime}}}$
$\displaystyle\approx$
$\displaystyle\underline{h}\thinspace\thinspace\underline{\overrightarrow{u}^{{}^{\prime}}}=0,$
(125) $\displaystyle\underline{h\overrightarrow{a}^{{}^{\prime}}}$
$\displaystyle\approx$
$\displaystyle\underline{h}\thinspace\thinspace\underline{\overrightarrow{a}^{{}^{\prime}}}=\underline{h}\overrightarrow{\nabla}_{x}\cdot(\underline{\overrightarrow{u}^{{}^{\prime}}\overrightarrow{u}^{{}^{\prime}}})$
(126)

That is, $\underline{h}$ is uncorrelated with both
$\underline{u}^{{}^{\prime}}_{i}$ and $\underline{a}^{{}^{\prime}}_{i}$, which
reproduces the averaged momentum equations with additional Reynolds stress
terms $\underline{u_{i}^{{}^{\prime}}u_{j}^{{}^{\prime}}}$ that can be closed
by means of any conventional macroscopic turbulence models. Thus, Eq. (124)
can now be rewritten as

$\partial_{t}\underline{h}+\overrightarrow{\eta}\cdot\overrightarrow{\nabla}_{x}\underline{h}=\Omega(\underline{h},\underline{h})+\overrightarrow{\nabla}_{x}\cdot(\underline{\overrightarrow{u}^{{}^{\prime}}\overrightarrow{u}^{{}^{\prime}}})\cdot\overrightarrow{\nabla}_{\eta}\underline{h},$
(127)

which represents the evolution of the statistical averaged/filtered turbulence
field by means of the Reynolds stresses that appear as a forcing term in a
kinetic framework.

Let us now develop a Galilean invariant lattice kinetic equation, i.e. which
provides inertial frame invariant representation with respect to the
_resolved_ velocity field obtained by statistical averaging/filtering. For
brevity and to avoid the use of additional new notations, let us rewrite Eq.
(127) by replacing $\underline{h}$ by $\underline{f}$ (and $\overline{\eta}$
and $\overline{\xi}$) to make use of the developments of the previous
sections. That is,

$\partial_{t}\underline{f}+\overrightarrow{\xi}\cdot\overrightarrow{\nabla}_{x}\underline{f}=\Omega(\underline{f},\underline{f})+\overrightarrow{\nabla}_{x}\cdot(\underline{\overrightarrow{u}^{{}^{\prime}}\overrightarrow{u}^{{}^{\prime}}})\cdot\overrightarrow{\nabla}_{\xi}\underline{f}.$
(128)

from which the resolved hydrodynamic fields can be defined as follows:

$\underline{\rho}=\int\underline{f}d\overrightarrow{\xi},\qquad\underline{\rho\overrightarrow{u}}=\int\underline{f}\overrightarrow{\xi}d\overrightarrow{\xi}.$
(129)

Now, we define operator averaged continuous central moments as

$\underline{\widehat{\Pi}}_{x^{m}y^{n}}=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\underline{f}(\xi_{x}-\underline{u_{x}})^{m}(\xi_{y}-\underline{u_{y}})^{n}d\xi_{x}d\xi_{y}$
(130)

and similarly for the continuous central moments of the Maxwellian
$\underline{\widehat{\Pi}}_{x^{m}y^{n}}^{\mathcal{M}}$ based on replacing
$h^{\mathcal{M}}$ and $\overrightarrow{\eta}$ by $f^{\mathcal{M}}$ and
$\overrightarrow{\xi}$, respectively. The Cartesian components of the
unresolved turbulent Reynolds stresses may be written as

$\displaystyle\underline{a_{x}^{{}^{\prime}}}$ $\displaystyle=$
$\displaystyle-\partial_{x}(\underline{u_{x}^{{}^{\prime}}u_{x}^{{}^{\prime}}})-\partial_{y}(\underline{u_{x}^{{}^{\prime}}u_{y}^{{}^{\prime}}}),$
(131) $\displaystyle\underline{a_{y}^{{}^{\prime}}}$ $\displaystyle=$
$\displaystyle-\partial_{x}(\underline{u_{x}^{{}^{\prime}}u_{y}^{{}^{\prime}}})-\partial_{y}(\underline{u_{y}^{{}^{\prime}}u_{y}^{{}^{\prime}}}),$
(132)

where
$\underline{\overrightarrow{a}^{{}^{\prime}}}=(\underline{a_{x}^{{}^{\prime}}},\underline{a_{y}^{{}^{\prime}}})$,
from which we analogously define a source/sink continuous central moment as

$\underline{\widehat{\Gamma}^{a}_{x^{m}y^{n}}}=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\underline{\delta
f^{a^{{}^{\prime}}}}(\xi_{x}-\underline{u_{x}})^{m}(\xi_{y}-\underline{u_{y}})^{n}d\xi_{x}d\xi_{y}.$
(133)

Here, $\underline{\delta
f^{a^{{}^{\prime}}}}=-\underline{\overrightarrow{a}^{{}^{\prime}}}\cdot\overrightarrow{\nabla}_{\xi}\underline{f}$.
It readily follows from Eq. (16) that
$\underline{\widehat{\Gamma}}^{a}_{x^{m}y^{n}}$ also satisfies the following
exact identity
$\underline{\widehat{\Gamma}}^{a}_{x^{m}y^{n}}=ma_{x}^{{}^{\prime}}\underline{\widehat{\Pi}}_{x^{m-1}y^{n}}+na_{y}^{{}^{\prime}}\underline{\widehat{\Pi}}_{x^{m}y^{n-1}}$.
That is, the statistical averaged/filtered central moment of sources/sinks due
to unresolved fields of a given order is dependent on the product of the
Cartesian components of the gradients of turbulent stresses with the lower
order central moments of the averaged/filtered distribution function. The
corresponding discrete central moment LBM can be devised by considering the
following averaged representation of discrete vectors supported by the
particle velocity set:
$\mathbf{\underline{f}}=\ket{\underline{f_{\alpha}}}=(\underline{f}_{0},\underline{f}_{1},\underline{f}_{2},\ldots,\underline{f}_{8})^{T}$,
$\mathbf{\widehat{\underline{g}}}=\ket{\widehat{\underline{g}}_{\alpha}}=(\widehat{\underline{g}}_{0},\widehat{\underline{g}}_{1},\widehat{\underline{g}}_{2},\ldots,\widehat{\underline{g}}_{8})^{T}$,
$\mathbf{\underline{S}}=\ket{\underline{S_{\alpha}}}=(\underline{S}_{0},\underline{S}_{1},\underline{S}_{2},\ldots,\underline{S}_{8})^{T}$,
and
$\bm{\underline{\Omega}}^{c}\equiv\bm{\underline{\Omega}}^{c}(\underline{\mathbf{f}},\mathbf{\underline{\widehat{g}}})=(\mathcal{K}\cdot\mathbf{\underline{\widehat{g}}})=(\underline{\Omega}_{0}^{c},\underline{\Omega}_{1}^{c},\underline{\Omega}_{2}^{c},\ldots,\underline{\Omega}_{8}^{c})^{T}$,
and invoking Galilean invariance matching principle, i.e. matching the
continuous and discrete central moments of various quantities at successively
higher orders as discussed in earlier sections. In particular, the statistical
averaged/filtered discrete collision operator $\bm{\underline{\Omega}}^{c}$
can be obtained by considering reduced compressibility effects and factorized
attractors as in Sec. 13. Furthermore, the corresponding source terms in
velocity space $\mathbf{\underline{S}}$ can be constructed using the procedure
outlined in Sec. 9. The operator averaged LBE, in terms of the transformed
distribution function $\underline{\overline{f}}_{\alpha}$ for improved
accuracy, can be finally written as

$\underline{\overline{f}}_{\alpha}(\overrightarrow{x}+\overrightarrow{e}_{\alpha},t+1)=\underline{\overline{f}}_{\alpha}(\overrightarrow{x},t)+\underline{\Omega}_{{\alpha}(\overrightarrow{x},t)}^{c}+\underline{S}_{{\alpha}(\overrightarrow{x},t)},$
(134)

where
$\underline{\overline{f}}_{\alpha}=\underline{f}_{\alpha}-\frac{1}{2}\underline{S}_{\alpha}$.
Here, as before, we have adopted the standard discretization for the streaming
step (see the comment following Eq. (22)). The resolved hydrodynamic fields in
the reduced compressibility formulation can then be obtained as

$\displaystyle\underline{\rho}$ $\displaystyle=$
$\displaystyle\sum_{\alpha=0}^{8}\underline{\overline{f}}_{\alpha}=\braket{\underline{\overline{f}}_{\alpha}}{\rho},$
(135) $\displaystyle\rho_{0}\underline{u}_{i}$ $\displaystyle=$
$\displaystyle\sum_{\alpha=0}^{8}\underline{\overline{f}}_{\alpha}e_{\alpha
i}+\frac{1}{2}\underline{\rho}\thinspace\thinspace\underline{a^{{}^{\prime}}_{i}}=\braket{\underline{\overline{f}}_{\alpha}}{e_{\alpha
i}}+\frac{1}{2}\underline{\rho}\thinspace\thinspace\underline{a^{{}^{\prime}}_{i}}.\qquad
i\in{x,y}$ (136)

This provides a minimal lattice kinetic equation for incorporating turbulence
models, where the unresolved turbulent motion are inertial frame invariant
with respect to the resolved fluid motion. Here, we clarify the meaning of
this statement as follows. Unlike other areas in fluid mechanics, where models
have been developed starting from continuous kinetic theory, its role for
fluid turbulence has been more limited. This is mainly due to the fact that
kinetic theory generally considers distinct scale separation of physical
processes. On the other hand, turbulence is a flow phenomenon intrinsically
containing a continuous spectrum of scales with no scale separation. As such,
therefore, turbulence modeling developments have to rely much on phenomenology
whose mathematical forms are then constrained by invariance principles (e.g.
material frame indifference and inertial frame invariance mentioned earlier in
the introduction) and realizability considerations [81, 90]. Thus, except for
some early models such as those based on mixing length concepts and derivation
of some recent phenomenological models (e.g. [54]) based on kinetic theory,
turbulence modeling developments are generally based on macroscopic models.
The ultimate goal of our central moment approach for the filtered kinetic
equation discussed above, is, then, to simulate resolved turbulent fields
which are inertial frame independent, when an appropriate macroscopic
turbulence model for the unresolve field is used in the forcing term.

##  15 Summary and Conclusions

A discrete lattice kinetic model for the continuous Boltzmann equation,
including forcing, based on central moments is derived. The collision operator
as well as the source term of this lattice Boltzmann equation is constructed
by matching the corresponding continuous and discrete central moments
successively at various orders. The local attractor of the collision operator
is constructed to satisfy the factorization property of the Maxwellian during
relaxation process. An exact hierarchical identity of the central moment of
sources, that incorporates non-equilibrium effects, is maintained at the
discrete level. The resulting approach provides Galilean invariant
hydrodynamic fields in the presence of any external or self-consistent
internal forces in a discrete kinetic framework. It is further extended to
incorporate reduced compressibility effects for better representation of
incompressible flow, a limiting case. An important physical characteristic of
turbulent flows is that it is inertial frame independent for all or any subset
of scales. In consequence, for general applicability, all turbulence models
and their simulation approaches, should satisfy this requirement. A
statistical averaged/filtered lattice kinetic equation based on central
moments that maintains Galilean invariant representation of unresolved fluid
motion with respect to the resolved fields of turbulent flow is thus
constructed. The formalism presented here can extended to other lattice
velocity sets and in three-dimensions as well as to other physical problems
such as complex fluids.

In this regard, we make the following remark on the development of more
efficient schemes for the former aspect. Symmetry and finiteness of the
standard lattice sets lead to degeneracies of higher order moments in terms
those at lower orders that can result in frame dependent contributions to
viscous stresses. This necessitates considerations of extended lattice sets.
In this case, it is proposed that the _central moments_ relaxation (as well as
forcing) be considered _only_ up to those higher order moments that have
bearing on the physics of _hydrodynamics_ , such as stress tensors and heat
flux vectors. In turn, this would impose Galilean invariance of the
macroscopic description of the fluid motion. The relaxation of the rest of the
higher moments (including forcing) related to the fast _kinetic or ghost
modes_ can be considered in terms of the standard or _raw moments_. Here, the
form of the hierarchical identity for the sources derived in this paper for
those higher (kinetic) moments would be the same with the simple replacement
of the central moment terms by the corresponding raw moments. It is envisaged
that such mixed central/raw moment approach would exhibit interesting
mathematical structures for the resulting collision operator and the sources,
as well as being computationally more effective. This strategy is currently
under investigation.

## References

  * [1] U. Frisch, B. Hasslacher and Y. Pomeau, Lattice-Gas Automata for the Navier-Stokes Equation, Phys. Rev. Lett. 56 (1986) 1505.
  * [2] G.R. McNamara and G. Zanetti, Use of the Boltzmann Equation to Simulate Lattice-Gas Automata, Phys. Rev. Lett. 61 (1988) 2332.
  * [3] R. Benzi, S. Succi and M. Vergassola, The Lattice Boltzmann Equation: Theory and Applications, Phys. Rep. 222 (1992) 145.
  * [4] S. Chen and G. Doolen, Lattice Boltzmann Method for Fluid Flows, Annu. Rev. Fluid Mech. 8 (1998) 2527.
  * [5] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Clarendon Press, Oxford, 2001.
  * [6] D. Yu, R. Mei, L.-S. Luo and W. Shyy, Viscous Flow Computations with the Method of Lattice Boltzmann Equation, Prog. Aero. Sci. 39 (2003) 329.
  * [7] T. Abe, Derivation of the Lattice Boltzmann Method by means of the discrete ordinate method for the Boltzmann equation, J. Comp. Phys. 131 (1997) 241.
  * [8] X. He and L.S. Luo, A priori derivation of the lattice Boltzmann equation, Phys. Rev. E 55 (1997) 115.
  * [9] X. Shan and X. He, Discretization of the Velocity Space in the Solution of the Boltzmann Equation, Phys. Rev. Lett. 80 (1998) 65.
  * [10] L.-S. Luo, Theory of the lattice Boltzmann method: Lattice Boltzmann models for nonideal gases, Phys. Rev. E 62 (2000) 4982.
  * [11] X. He and G. Doolen, Thermodynamic Foundations of Kinetic Theory and Lattice Boltzmann Models for Multiphase Flows, J. Stat. Phys. 107 (2002) 1572.
  * [12] P. Asinari and L.-S. Luo, A Consistent Lattice Boltzmann Equation with Baroclinic Coupling for Mixtures, J. Comput. Phys. 227 (2008) 3878.
  * [13] X. Shan, X.F. Yuan and H. Chen, Kinetic Theory Representation of Hydrodynamics: A Way Beyond the Navier-Stokes Equation, J. Fluid Mech. 550 (2006) 413.
  * [14] M. Junk, A. Klar and L.S. Luo, Asymptotic Analysis of the Lattice Boltzmann Equation, J. Comp. Phys. 210 (2005) 676.
  * [15] Y. Qian, D. d’Humières and P. Lallemand, Lattice BGK Models for Navier-Stokes Equation, Europhys. Lett. 17 (1992) 479.
  * [16] H. Chen, S. Chen and W.H. Matthaeus, Recover of the Navier-Stokes Eqauations using the Lattice-Gas Boltzmann Method, Phys. Rev. A 45 (1992) 5339.
  * [17] P. Bhatnagar, E. Gross and M. Krook, A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems, Phys. Rev. 94 (1954) 511.
  * [18] D. d‘Humières, Generalized Lattice Boltzmann Equations, Progress in Aeronautics and Astronautics (Eds. B.D. Shigal and D.P Weaver) 159 (1992) 450.
  * [19] F. J. Higuera and J. Jiménez, Boltzmann Approach to Lattice Gas Simulations, Europhys. Lett. 9 (1989) 663.
  * [20] F. J. Higuera, S. Succi and R. Benzi, Lattice gas-dynamics with enhanced collisions, Europhys. Lett. 9 (1989) 345.
  * [21] H. Grad, On the Kinetic Theory of Rarefied Gases, Comm. Pure Appl. Math. 2 (1949) 331.
  * [22] P. Lallemand and Li-Shi Luo, Theory of the Lattice Boltzmann Method: Dispersion, Dissipation, Isotropy, Galilean Invariance, and Stability, Phys. Rev. E 61 (2000) 6546.
  * [23] D. d‘Humières, I. Ginzburg, M. Krafczyk, P. Lallemand and L.-S. Luo, Multiple-Relaxation-Time Lattice Boltzmann Models in Three Dimensions, Phil. Trans. R. Soc. Lond. A 360 (2002) 437.
  * [24] K.N. Premnath and J. Abraham, Three-dimensional Multi-Relaxation Time (MRT) Lattice-Boltzmann Models for Multiphase Flow, J. Comp. Phys. 224 (2007) 539.
  * [25] K.N. Premnath and J. Abraham, Simulations of Binary Drop Collisions with a Multiple-Relaxation-Time Lattice-Boltzmann Model, Phys. Fluids 17 (2005) 122105.
  * [26] K.N. Premnath, M.E. McCracken and J. Abraham, A Review of Lattice Boltzmann Methods for Multiphase Flows Relevant to Engine Sprays, Trans. of the SAE: J. Engines 114 (2005) 929.
  * [27] K.N. Premnath, M.J. Pattison and S. Banerjee, Generalized Lattice Boltzmann Equation with Forcing Term for Computation of Wall Bounded Turbulent Flows, Phys. Rev. E 79 (2009) 026703.
  * [28] K.N. Premnath, M.J. Pattison and S. Banerjee, Dynamic Subgrid Scale Modeling of Turbulent Flows using Lattice-Boltzmann Method, Physica A 388 (2009) 2640.
  * [29] M.J. Pattison, K.N. Premnath, N.B. Morley and M. Abdou, Progress in Lattice Boltzmann Methods for Magnetohydrodynamic Flows Relevant to Fusion Applications, Fusion Engg. Des. 83 (2008) 557.
  * [30] X. Shan and H. Chen, A General Multiple-Relaxation Boltzmann Collision Model, Int. J. Mod. Phys. C 18 (2007) 635.
  * [31] I.V. Karlin, A. Ferrante H.C.Öttinger, Perfect Entropy Functions of the Lattice Boltzmann Method, Euro. Phys. Lett 47 (1999) 182.
  * [32] B.M. Boghosian, J. Yepez, P.V. Coveney and A. Wagner, Entropic Lattice Boltzmann Methods, Proc. R. Soc. Lond. A 457 (2001) 717.
  * [33] S. Ansumali and I.V. Karlin, Single Relaxation Time Model for Entropic Lattice Boltzmann Methods, Phys. Rev. E 65 (2002) 056312.
  * [34] S. Succi, I. Karlin and H. Chen, Role of the H theorem in Lattice Boltzmann Hydrodynamic Simulations, Rev. Mod. Phys. 74 (2002) 1203.
  * [35] P. Asinari and I.V. Karlin, Generalized Maxwell State and H Theorem for Computing Fluid Flows using the Lattice Boltzmann Method, Phys. Rev. E 79 (2009) 036703.
  * [36] P. Asinari and I.V. Karlin, Quasiequilibrium Lattice Boltzmann Models with Tunable Bulk Viscosity for Enhancing Stability, Phys. Rev. E 81 (2010) 016702.
  * [37] S.S. Chikatamarla and I.V. Karlin, Lattices for the Lattice Boltzmann Method, Phys. Rev. E 79 (2009) 046701.
  * [38] F. Dubois and P. Lallemand, Towards Higher Order Lattice Boltzmann Schemes, J. Stat. Mech.: Theory and Experiment P06006 (2009) 1.
  * [39] I.V. Karlin and P. Asinari, Factorization Symmetry in the Lattice Boltzmann Method, Physica A 389 (2010) 1530.
  * [40] R. Rubinstein and L.-S. Luo, Theory of the Lattice Boltzmann Equation: Symmetry Properties of Discrete Velocity Sets, Phys. Rev. E 77 (2008) 036709.
  * [41] H. Chen and I. Goldhirsch and S.A. Orszag, Discrete Rotational Symmetry, Moment Isotropy and Higher Order Lattice Boltzmann Models, J. Sci. Comput. 34 (2008) 87.
  * [42] M.K. Banda and W.-A. Yong and A. Klar, A Stability Notion for Lattice Boltzmann Equations, SIAM J. Sci. Comput. 27 (2006) 2098.
  * [43] M. Junk and W.-A. Yong, Weighted $\mathbb{L}^{2}$-Stability of the Lattice Boltzmann Method, SIAM J. Numer. Anal. 47 (2009) 1651.
  * [44] W.-A. Yong, Onsager-like Relation for the Lattice Boltzmann Method, Comp. Math. Appl. 58 (2009) 862.
  * [45] M. Junk and Z. Yang, Convergence of Lattice Boltzmann Methods for Navier-Stokes Flows in Periodic and Bounded Domains, Numer. Math. 112 (2009) 65.
  * [46] C. Truesdell and W. Noll, The Non-Linear Field Theories of Mechanics, Springer, New York, 2004.
  * [47] I. Muller, On the Frame Dependence of Stress and Heat Flux, Arch. Rational Mech. Anal. 45 (1972) 45.
  * [48] C.-C. Wang, On the Concept of Frame-Indifference in Continuum Mechanics and in the Kinetic Theory of Gases, Arch. Rational Mech. Anal. 58 (1975) 381.
  * [49] C.G. Speziale, On Frame-Indifference and Iterative Procedures in the Kinetic Theory of Gases, Int. J. Engg. Sci. 19 (1981) 63.
  * [50] L.C. Woods, Frame-Indifferent Kinetic Theory, J. Fluid Mech. 136 (1983) 423.
  * [51] C.G. Speziale, The Effect of the Reference Frame on the Thermophysical Properties of an Ideal Gas, Int. J. Thermophys. 7 (1986) 99.
  * [52] C.G. Speziale, Comments on the “Material Frame-Indifference” Controversy, Phys. Rev. A 36 (1987) 4522.
  * [53] A.J. Wagner, Derivation of a Non-Objective Oldroyd Model from the Boltzmann Equation, arXiv:cond-mat/0105067v1 (2001).
  * [54] P. Degond and M. Lemou, Turbulence Models for Incompressible Fluids Derived from Kinetic Theory, J. Math. Fluid Mech. 4 (2002) 257.
  * [55] N. Cao, S. Chen, S. Jin, and D. Mart nez, Physical Symmetry and Lattice Symmetry in the Lattice Boltzmann Method, Phys. Rev. E 55 (1997) R21.
  * [56] V. Sofonea and R.F. Sekerka, Viscosity of Finite Difference Lattice Boltzmann Models, J. Comput. Phys. 184 (2003) 422.
  * [57] T. Lee and C.-L. Lin, An Eulerian Description of the Streaming Process in the Lattice Boltzmann Equation, J. Comput. Phys. 185 (2003) 445.
  * [58] M. Min and T. Lee, A Spectral-Element Discontinuous Galerkin Lattice Boltzmann Method for Nearly Incompressible Flows, J. Comput. Phys. 230 (2011) 245.
  * [59] M. Geier, A. Greiner and J.G. Korvink, Cascaded Digital Lattice Boltzmann Automata for High Reynolds Number Flow, Phys. Rev. E 73 (2006) 066705.
  * [60] H. Struchtrup, Macroscopic Transport Equations for Rarefied Gas Flows: Approximation Methods in Kinetic Theory, Springer, New York, 2005.
  * [61] P. Asinari, Generalized Local Equilibrium in the Cascaded Lattice Boltzmann Method, Phys. Rev. E 78 (2008) 016701.
  * [62] M. Geier, A. Greiner and J.G. Korvink, A Factorized Central Moment Lattice Boltzmann Method, Eur. Phys. J. Special Topics 171 (2009) 55.
  * [63] K.N. Premnath and S. Banerjee, Incorporating Forcing Terms in Cascaded Lattice Boltzmann Approach by Method of Central Moments, Phys. Rev. E 80 (2009) 036702.
  * [64] S. Chapman and T.G. Cowling, Mathematical Theory of Nonuniform Gases, Cambridge University Press, New York (1964).
  * [65] K.N. Premnath and S. Banerjee, On the Three Dimensional Central Moment Lattice Boltzmann Method, J. Stat. Phys. 143 (2011) 747.
  * [66] G. Scovazzi, A Discourse on Galilean Invariance and SUPG-Type Stabilization, Comp. Methods Appl. Mech. Engg. 196 (2007) 1108.
  * [67] G. Scovazzi, Galilean Invariance and Stabilized Methods for Compressible Flows, Int. J. Num. Meth. Fluids 54 (2007) 757.
  * [68] H. Chen, S. Succi and S. Orszag, Analysis of Subgrid Scale Turbulence using the Boltzmann Bhatnagar-Gross-Krook Kinetic Equation, Phys. Rev. E 59 (1998) R2527.
  * [69] S. Ansumali, I. Karlin and S. Succi, Kinetic Theory of Turbulence Modelling: Smallness Parameter, Scaling and Microscopic Derivation of Smagorinsky Model, Physica A 338 (2004) 379.
  * [70] H. Chen, S. Orszag and I. Staroselsky, Expanded Analogy between Boltzmann Kinetic Theory of Fluids and Turbulence, J. Fluid Mech. 519 (2004) 301.
  * [71] C. Teixeira, Incorporating Turbulence Models into the Lattice-Boltzmann Method, Int. J. Mod. Phys. C 9 (1998) 1159.
  * [72] H. Chen, S. Kandasamy, S. Orszag, R. Shock, S. Succi and V. Yakhot, Extended Boltzmann Kinetic Equation for Turbulent Flows, Science 301 (2003) 633.
  * [73] S. Hou, J. Sterling, S. Chen and G.D. Doolen, A Lattice Boltzmann Subgrid Model for High Reynolds Number Flows, Fields Inst. Comm. 6 (1996) 151.
  * [74] H. Yu, L.-S. and S. Girimaji, LES of Turbulent Square Jet Flow using an MRT Lattice Boltzmann Model, Comput. Fluids 35 (2006) 957.
  * [75] Y.H. Dong, P. Sagaut and S. Marie, Inertial Consistent Subgrid Model for Large-Eddy Simulation Based on the Lattice Boltzmann Method, Phys. Fluids 20 (2008) 035104.
  * [76] P. Sagaut, Toward Advanced Subgrid Models for Lattice-Boltzmann-Based Large-Eddy Simulation: Theoretical Formulations, Comp. Math. Appl. 59 (2010) 2194.
  * [77] S. Girimaji, Boltzmann Kinetic Equation for Filtered Fluid Turbulence, Phys. Rev. Lett. 99 (2007) 034501.
  * [78] S. Stolz and N.A. Adams, An Approximate Deconvolution Procedure for Large-Eddy Simulation, Phys. Fluids 11 (1999) 1699.
  * [79] C.G. Speziale, Galilean Invariance of Subgrid-Scale Stress Models in the Large-Eddy Simulation of Turbulence, J. Fluid Mech. 156 (1985) 55.
  * [80] M. Germano, A Proposal for a Redefinition of the Turbulent Stresses in the Filtered Navier-Stokes Equations, Phys. Fluids 29 (1986) 2323.
  * [81] C. Fureby and G. Tabor, Mathematical and Physical Constraints on Large-Eddy Simulations, Theoret. Comput. Fluid Dynamics 9 (1997) 85.
  * [82] A.S. Monin and A.M. Yaglom, Statistical Fluid Mechanics, Vol I: Mechanics of Turbulence, Dover Publications, 2007.
  * [83] M. Germano, Turbulence: The Filtering Approach, J. Fluid Mech. 238 (1992) 325.
  * [84] S. Harris, An Introduction to the Theory of the Boltzmann Equation, Dover Publications, New York, 1999.
  * [85] M.N. Kogan, Rarefied Gas Dynamics, Plenum Press, New York, 1969.
  * [86] I. Müller and T. Ruggeri, Extended Thermodynamics, Springer Tracts in Natural Philosophy, New York, 1993.
  * [87] X. He and L.-S. Luo, Lattice Boltzmann Model for the Incompressible Navier-Stokes Equation, J. Stat. Phys. 88 (1997) 927.
  * [88] J.C. Maxwell, Illustrations of the Dynamical Theory of Gases, Phil. Mag. 19 (1860) 19.
  * [89] O. Reynolds, On the Dynamical Theory of Incompressible Viscous Fluids and the Determination of the Criterion, Phil. Trans. Roy. Soc. London A 186 (1895) 123.
  * [90] P.A. Durbin and B. A. Pettersson Reif, Statistical Theory and Modeling for Turbulent Flows, Wiley, West Sussex, 2010.