Lattice Boltzmann inverse kinetic approach for the incompressible
  Navier-Stokes equations
Lattice Boltzmann inverse kinetic approach for the incompressible Navier-Stokes
equations
Enrico Fonda1,Massimo Tessarotto1,2 and Marco Ellero3
1Dipartimento di Matematica e Informatica,
Università di Trieste, Italy
2Consorzio di Magnetofluidodinamica, Trieste, Italy
3Institute of Aerodynamics,
Technical University of Munich, Munich, Germany
(Dated: August 18, 2021)
In spite of the large number of papers appeared in the past which are devoted to the lattice
Boltzmann (LB) methods, basic aspects of the theory still remain unchallenged. An unsolved theo-
retical issue is related to the construction of a discrete kinetic theory which yields exactly the fluid
equations, i.e., is non-asymptotic (here denoted as LB inverse kinetic theory). The purpose of this
paper is theoretical and aims at developing an inverse kinetic approach of this type. In principle
infinite solutions exist to this problem but the freedom can be exploited in order to meet important
requirements. In particular, the discrete kinetic theory can be defined so that it yields exactly the
fluid equation also for arbitrary non-equilibrium (but suitably smooth) kinetic distribution func-
tions and arbitrarily close to the boundary of the fluid domain. This includes the specification
of the kinetic initial and boundary conditions which are consistent with the initial and boundary
conditions prescribed for the fluid fields. Other basic features are the arbitrariness of the ”equi-
librium” distribution function and the condition of positivity imposed on the kinetic distribution
function. The latter can be achieved by imposing a suitable entropic principle, realized by means of
a constant H-theorem. Unlike previous entropic LB methods the theorem can be obtained without
functional constraints on the class of the initial distribution functions. As a basic consequence, the
choice of the the entropy functional remains essentially arbitrary so that it can be identified with
the Gibbs-Shannon entropy. Remarkably, this property is not affected by the particular choice of
the kinetic equilibrium (to be assumed in all cases strictly positive). Hence, it applies also in the
case of polynomial equilibria, usually adopted in customary LB approaches. We provide different
possible realizations of the theory and asymptotic approximations which permit to determine the
fluid equations with prescribed accuracy. As a result, asymptotic accuracy estimates of customary
LB approaches and comparisons with the Chorin artificial compressibility method are discussed.
PACS numbers: 47.27.Ak, 47.27.eb, 47.27.ed
1 - INTRODUCTION - INVERSE KINETIC
THEORIES
Basic issues concerning the foundations classical hy-
drodynamics still remain unanswered. A remarkable as-
pect is related the construction of inverse kinetic theo-
ries (IKT) for hydrodynamic equations in which the fluid
fields are identified with suitable moments of an appropri-
ate kinetic probability distribution. The topic has been
the subject of theoretical investigations both regarding
the incompressible Navier-Stokes (NS) equations (INSE)
[1, 2, 3, 4, 5, 6] and the quantum hydrodynamic equations
associated to the Schrödinger equation [7]. The impor-
tance of the IKT-approach for classical hydrodynamics
goes beyond the academic interest. In fact, INSE rep-
resent a mixture of hyperbolic and elliptic pde’s, which
are extremely hard to study both analytically and nu-
merically. As such, their investigation represents a chal-
lenge both for mathematical analysis and for computa-
tional fluid dynamics. The discovery of IKT [1] provides,
however, a new starting point for the theoretical and nu-
merical investigation of INSE. In fact, an inverse kinetic
theory yields, by definition, an exact solver for the fluid
equations : all the fluid fields, including the fluid pres-
sure p(r, t), are uniquely prescribed in terms of suitable
momenta of the kinetic distribution function, solution
of the kinetic equation. In the case of INSE this per-
mits, in principle, to determine the evolution of the fluid
fields without solving explicitly the Navier-Stokes equa-
tion, nor the Poisson equations for the fluid pressure [6].
Previous IKT approaches [2, 3, 4, 5, 7] have been based
on continuous phase-space models. However, the inter-
esting question arises whether similar concepts can be
adopted also to the development of discrete inverse ki-
netic theories based on the lattice Boltzmann (LB) the-
ory. The goal of this investigation is to propose a novel
LB theory for INSE, based on the development of an IKT
with discrete velocities, here denoted as lattice Boltzmann
inverse kinetic theory (LB-IKT). In this paper we intend
to analyze the theoretical foundations and basic proper-
ties of the new approach useful to display its relation-
ship with previous CFD and lattice Boltzmann methods
(LBM) for incompressible isothermal fluids. In particu-
lar, we wish to prove that it delivers an inverse kinetic
http://arxiv.org/abs/0704.0339v1
theory, i.e., that it realizes an exact Navier-Stokes and
Poisson solver.
1a - Motivations: difficulties with LBM’s
Despite the significant number of theoretical and nu-
merical papers appeared in the literature in the last few
years, the lattice Boltzmann method [8, 9, 10, 11, 12,
13, 14] - among many others available in CFD - is prob-
ably the one for which a complete understanding is not
yet available. Although originated as an extension of
the lattice gas automaton [15, 16] or a special discrete
form of the Boltzmann equation [17], several aspects re-
garding the very foundation of LB theory still remain to
be clarified. Consequently, also the comparisons and ex-
act relationship between the various lattice Boltzmann
methods (LBM) and other CFD methods are made dif-
ficult or, at least, not yet well understood. Needless to
say, these comparisons are essential to assess the relative
value (based on the characteristic computational com-
plexity, accuracy and stability) of LBM and other CFD
methods. In particular the relative performance of the
numerical methods depend strongly on the characteris-
tic spatial and time discretization scales, i.e., the minimal
spatial and time scale lengths required by each numerical
method to achieve a prescribed accuracy. On the other
hand, most of the existing knowledge of the LBM’s prop-
erties originates from numerical benchmarks (see for ex-
ample [18, 19, 20]). Although these studies have demon-
strated the LBM’s accuracy in simulating fluid flows, few
comparisons are available on the relative computational
efficiency of the LBM and other CFD methods [17, 21].
The main reason [of these difficulties] is probably because
current LBM’s, rather than being exact Navier-Stokes
solvers, are at most asymptotic ones (asymptotic LBM’s),
i.e., they depend on one or more infinitesimal parame-
ters and recover INSE only in an approximate asymptotic
sense.
The motivations of this work are related to some of
the basic features of customary LB theory representing,
at the same time, assets and weaknesses. One of the
main reasons of the popularity of the LB approach lays
in its simplicity and in the fact that it provides an ap-
proximate Poisson solver, i.e., it permits to advance in
time the fluid fields without explicitly solving numeri-
cally the Poisson equation for the fluid pressure. How-
ever customary LB approaches can yield, at most, only
asymptotic approximations for the fluid fields. This is
because of two different reasons. The first one is the dif-
ficulty in the precise definition of the kinetic boundary
conditions in customary LBM’s, since sufficiently close to
the boundary the form of the distribution function pre-
scribed by the boundary conditions is not generally con-
sistent with hydrodynamic equations. The second reason
is that the kinetic description adopted implies either the
introduction of weak compressibility [8, 9, 11, 12, 13, 14]
or temperature [22] effects of the fluid or some sort of
state equation for the fluid pressure [23]. These assump-
tions, although physically plausible, appear unacceptable
from the mathematical viewpoint since they represent a
breaking of the exact fluid equations.
Moreover, in the case of very small fluid viscosity
customary LBM’s may become inefficient as a conse-
quence of the low-order approximations usually adopted
and the possible presence of the numerical instabilities
mentioned above. These accuracy limitations at low vis-
cosities can usually be overcome only by imposing severe
grid refinements and strong reductions of the size of the
time step. This has the inevitable consequence of rais-
ing significantly the level of computational complexity
in customary LBM’s (potentially much higher than that
of so-called direct solution methods), which makes them
inefficient or even potentially unsuitable for large-scale
simulations in fluids.
A fundamental issue is, therefore, related to the con-
struction of more accurate, or higher-order, LBM’s, ap-
plicable for arbitrary values of the relevant physical
(and asymptotic) parameters. However, the route which
should permit to determine them is still uncertain, since
the very existence of an underlying exact (and non-
asymptotic) discrete kinetic theory, analogous to the con-
tinuous inverse kinetic theory [2, 3], is not yet known.
According to some authors [24, 25, 26] this should be
linked to the discretization of the Boltzmann equation, or
to the possible introduction of weakly compressible and
thermal flow models. However, the first approach is not
only extremely hard to implement [27], since it is based
on the adoption of higher-order Gauss-Hermite quadra-
tures (linked to the discretization of the Boltzmann equa-
tion), but its truncations yield at most asymptotic the-
ories. Other approaches, which are based on ’ad hoc’
modifications of the fluid equations (for example, intro-
ducing compressibility and/or temperature effects [28]),
by definition cannot provide exact Navier-Stokes solvers.
Another critical issue is related to the numerical sta-
bility of LBM’s [29], usually attributed to the violation of
the condition of strict positivity (realizability condition)
for the kinetic distribution function [29, 30]. Therefore,
according to this viewpoint, a stability criterion should
be achieved by imposing the existence of an H-theorem
(for a review see [31]). In an effort to improve the ef-
ficiency of LBM numerical implementations and to cure
these instabilities, there has been recently a renewed in-
terest in the LB theory. Several approaches have been
proposed. The first one involves the adoption of entropic
LBM’s (ELBM [30, 32, 33, 34] in which the equilibrium
distribution satisfies also a maximum principle, defined
with respect to a suitably defined entropy functional.
However, usually these methods lead to non-polynomial
equilibrium distribution functions which potentially re-
sult in higher computational complexity [35] and less nu-
merical accuracy[36]. Other approaches rely on the adop-
tion of multiple relaxation times [37, 38]. However the
efficiency, of these methods is still in doubt. Therefore,
the search for new [LB] models, overcoming these limita-
tions, remains an important unsolved task.
1b - Goals of the investigation
The aim of this work is the development of an inverse
kinetic theory for the incompressible Navier-Stokes equa-
tions (INSE) which, besides realizing an exact Navier-
Stokes (and Poisson) solver, overcomes some of the lim-
itations of previous LBM’s. Unlike Refs. [2, 3], where a
continuous IKT was considered, here we construct a dis-
crete theory based on the LB velocity-space discretiza-
tion. In such a type of approach, the kinetic description
is realized by a finite number of discrete distribution func-
tions fi(r, t), for i = 0, k, each associated to a prescribed
discrete constant velocity ai and defined everywhere in
the existence domain of the fluid fields (the open set Ω×I
). The configuration space Ω is a bounded subset of the
Euclidean space R3and the time interval I is a subset of
R. The kinetic theory is obtained as in [2, 3] by introduc-
ing an inverse kinetic equation (LB-IKE) which advances
in time the distribution function and by properly defin-
ing a correspondence principle, relating a set of velocity
momenta with the relevant fluid fields.
To achieve an IKT for INSE, however, also a proper
treatment of the initial and boundary conditions, to be
satisfied by the kinetic distribution function, must be in-
cluded. In both cases, it is proven that they can be de-
fined to be exactly consistent - at the same time - both
with the hydrodynamic equations (which must hold also
arbitrarily close to the boundary of the fluid domain) and
with the prescription of the initial and Dirichlet bound-
ary conditions set for the fluid fields. Remarkably, both
the choice of the initial and equilibrium kinetic distri-
bution functions and their functional class remain essen-
tially arbitrary. In other words, provided suitable min-
imal smoothness conditions are met by the kinetic dis-
tributions function, for arbitrary initial and boundary
kinetic distribution functions, the relevant moment equa-
tions of the kinetic equation coincide identically with the
relevant fluid equations. This includes the possibility
of defining a LB-IKT in which the kinetic distribution
function is not necessarily a Galilean invariant.
This arbitrariness is reflected also in the choice of pos-
sible ”equilibrium” distribution functions, which remain
essentially free in our theory, and can be made for exam-
ple in order to achieve minimal algorithmic complexity.
A possible solution corresponds to assume polynomial-
type kinetic equilibria, as in the traditional asymptotic
LBM’s. These kinetic equilibria are well-known to be
non-Galilean invariant with respect to arbitrary finite
velocity translations. Nevertheless, as discussed in detail
in Sec.4, Subsection 4A, although the adoption of Galilei
invariant kinetic distributions is in possible, this choice
does not represent an obstacle for the formulation of a
LB-IKT. Actually Galilean invariance need to be fulfilled
only by the fluid equations. The same invariance prop-
erty must be fulfilled only by the moment equations of
the LB-IKT and not necessarily by the whole LB inverse
kinetic equation (LB-IKE).
Another significant development of the theory is the
formal introduction of an entropic principle, realized by
a constant H-theorem, in order to assure the strict pos-
itivity of the kinetic distribution function in the whole
existence domain Ω× I. The present entropic principle
departs significantly from the literature. Unlike previ-
ous entropic LBM’s it is obtained without imposing any
functional constraints on the class of the initial kinetic
distribution functions. Namely without demanding the
validity of a principle of entropy maximization (PEM,
[39]) in a true functional sense on the form of the distri-
bution function. Rather, it follows imposing a constraint
only on a suitable set of extended fluid fields, in particu-
lar the kinetic pressure p1(r, t).The latter is uniquely re-
lated to the actual fluid pressure p(r, t) via the equation
p1(r, t) = p(r, t) + Po(t), with Po(t) > 0 to be denoted
as pseudo-pressure. The constant H-theorem is therefore
obtained by suitably prescribing the function Po(t) and
implies the strict positivity. The same prescription as-
sures that the entropy results maximal with respect in the
class of the admissible kinetic pressures, i.e., it satisfies a
principle of entropy maximization. Remarkably, since
this property is not affected by the particular choice of
the kinetic equilibrium, the H-theorem applies also in the
case of polynomial equilibria. We stress that the choice
of the entropy functional remains essentially arbitrary,
since no actual physical interpretation can be attached to
it. For example, without loss of generality it can always
be identified with the Gibbs-Shannon entropy. Even pre-
scribing these additional properties, in principle infinite
solutions exist to the problem. Hence, the freedom can
be exploited to satisfy further requirements (for example,
mathematical simplicity, minimal algorithmic complex-
ity, etc.). Different possible realizations of the theory and
comparisons with other CFD approaches are considered.
The formulation of the inverse kinetic theory is also use-
ful in order to determine the precise relationship between
the LBM’s and previous CFD schemes and in particular
to obtain possible improved asymptotic LBM’s with pre-
scribed accuracy. As an application, we intend to con-
struct asymptotic models which satisfy with prescribed
accuracy the required fluid equations [INSE] and possi-
bly extend also the range of validity of traditional LBM’s.
In particular, this permits to obtain asymptotic accuracy
estimates of customary LB approaches. The scheme of
presentation is as follows. In Sec.2 the INSE problem
is recalled and the definition of the extended fluid fields
{V, p1} is presented. In Sec. 3 the basic assumptions
of previous asymptotic LBM’s are recalled. In.Sec.4 and
5 the foundations of the new inverse kinetic theory are
laid down and the integral LB inverse kinetic theory is
presented, while in Sec. 6 the entropic theorem is proven
to hold for the kinetic distribution function for properly
defined kinetic pressure. Finally, in Sec.7 various asymp-
totic approximations are obtained for the inverse kinetic
theory and comparisons are introduce with previous LB
and CFD methods and in Sec. 8 the main conclusions
are drawn.
2 - THE INSE PROBLEM
A prerequisite for the formulation of an inverse kinetic
theory [2, 3] providing a phase-space description of a clas-
sical (or quantum) fluid is the proper identification of the
complete set of fluid equations and of the related fluid
fields. For a Newtonian incompressible fluid, referred to
an arbitrary inertial reference frame, these are provided
by the incompressible Navier-Stokes equations (INSE) for
the fluid fields {ρ,V,p}
∇ ·V = 0, (1)
NV = 0, (2)
ρ(r,t) = ρo. (3)
There are supplemented by the inequalities
p(r,t) ≥ 0, (4)
ρo > 0. (5)
Equations (1)-(3) are defined in a open connected set
Ω ⊆ R3 (defined as the subset of R3 where ρ(r,t) > 0)
with boundary δΩ, while Eqs. (4) and (5) apply on its
closure Ω. Here the notation is standard. Thus, N is the
NS operator
NV ≡ρo
V +∇p+ f − µ∇2V, (6)
with D
+V · ∇ the convective derivative, f denotes
a suitably smooth volume force density acting on the fluid
element and µ ≡ νρo > 0 is the constant fluid viscosity.
In particular we shall assume that f can be represented
in the form
f = −∇Φ(r) + f1(r,t)
where we have separated the conservative ∇Φ(r) and the
non-conservative f1 parts of the force. Equations (1)-(3)
are assumed to admit a strong solution in Ω × I, with
I ⊂ R a possibly bounded time interval. By assumption
{ρ,V,p} are continuous in the closure Ω. Hence if in Ω×I,
f is at least C(1,0)(Ω×I), it follows necessarily that {V,p}
must be at least C(2,1)(Ω × I). In the sequel we shall
impose on {V,p} the initial conditions
V(r,to) = Vo(r), (7)
p(r, to) = po(r).
Furthermore, for greater mathematical simplicity, here
we shall impose Dirichlet boundary conditions on δΩ
V(·,t)|
δΩ = VW (·,t)|δΩ
p(·,t)|
δΩ = pW (·,t)|δΩ .
Eqs.(3) and (7)-(8) define the initial-boundary value
problem associated to the reduced INSE (reduced INSE
problem). It is important to stress that the previous
problem can also formulated in an equivalent way by re-
placing the fluid pressure p(r, t) with a function p1(r, t)
(denoted kinetic pressure) of the form
p1(r, t) = Po + p(r, t), (9)
where Po = Po(t) is prescribed (but arbitrary) real func-
tion of time and is at least Po(t) ∈ C
(1)(I). {V,p1} will
be denoted hereon as extended fluid fields and Po(t) will
be denoted as pseudo-pressure.
3 - ASYMPTOTIC LBM’S
3A - Basic assumptions
As is well known, all LB methods are based on a dis-
crete kinetic theory, using a so-called lattice Boltzmann
velocity discretization of phase-space (LB discretization).
This involves the definition of a kinetic distribution func-
tion f, which can only take the values belonging to a
finite discrete set {fi(r, t), i = 0, k} (discrete kinetic dis-
tribution functions). In particular, it is assumed that the
functions fi, for i = 0, k, are associated to a discrete set
of k+1 different ”velocities” {ai, i = 0, k} . Each ai is an
’a priori’ prescribed constant vector spanning the vector
space Rn (with n = 2 or 3 respectively for the treatment
of two- and three-dimensional fluid dynamics),and each
fi(r, t) is represented by a suitably smooth real function
which is defined and continuous in Ω×I and in particular
is at least C(k,j)(Ω× I) with k ≥ 3.
The crucial aspect which characterizes customary LB
approaches [8, 9, 10, 11, 12, 13, 14, 17, 40, 41] involves the
construction of kinetic models which allow a finite sound
speed in the fluid and hence are based on the assumption
of a (weak) compressibility of the same fluid. This is
realized by assuming that the evolution equation (kinetic
equation) for the discrete distributions fi(r, t) (i = 1, k),
depends at least one (or more) infinitesimal (asymptotic)
parameters (see below). Such approaches are therefore
denoted as asymptotic LBM’s. They are characterized
by a suitable set of assumptions, which typically include:
1. LB assumption #1: discrete kinetic equation and
correspondence principle: the first assumption con-
cerns the definition of an appropriate evolution
equation for each fi(r, t) which must hold (together
with all its moment equations) in the whole open
set Ω× I. In customary LB approaches it takes the
form of the so-called LB-BGK equation [13, 41, 42]
L(i)fi = Ωi(fi), (10)
where i = 0, k. Here L(i) is a suitable streaming
operator,
Ωi(fi) = −νc(fi − f
i ) (11)
(with νc ≥ 0 a constant collision frequency) is
known as BKG collision operator (after Bhatba-
gar, Gross and Krook [43]) and f
i is an ”equi-
librium” distribution to be suitably defined. In
customary LBM’s it is implicitly assumed that
the solution of Eq.(10), subject to suitable initial
and boundary conditions exists and is unique in
the functional class indicated above. In partic-
ular, usually L(i) is either identified with the fi-
nite difference streaming operator (see for example
[8, 11, 13, 42]), i.e., L(i)fi(r, t) = LFD(i)fi(r, t) ≡
[fi(r+ ai∆t, t+∆t)− fi(r, t)] or with the dif-
ferential streaming operator (see for instance [17,
40, 41])
L(i) = LD(i) ≡
+ ai ·
. (12)
Here the notation is standard. In particular, in the
case of the operator LFD(i), ∆t and c∆t ≡ Lo are
appropriate parameters which define respectively
the characteristic time- and length- scales associ-
ated to the LBM time and spatial discretizations.
A common element to all LBM’s is the assump-
tion that all relevant fluid fields can be identified,
at least in some approximate sense, with appro-
priate momenta of the discrete kinetic distribu-
tion function (correspondence principle). In par-
ticular, for neutral and isothermal incompressible
fluids, for which the fluid fields are provided re-
spectively by the velocity and pressure fluid fields
{Yj(r, t), j = 1, 4} ≡ {V(r, t), p(r, t)} , it is as-
sumed that they are identified with a suitable set
of discrete velocity momenta (for j = 1, 4)
Yj(r, t) =
i=0,k
Xji(r, t)fi(r, t), (13)
where Xji(r, t) (with i = 0, k and j = 1, k) are ap-
propriate, smooth real weight functions. In the
literature several examples of correspondence prin-
ciples are provided, a particular case being provided
by the so-called D2Q9 (V, p)-scheme [44, 45]
p(r, t) = c2
i=0,k
fi = c
i=0,k
i , (14)
V(r,t) =
i=1,k
aifi =
i=1,k
i , (15)
where k = 8 and c = min {|ai| > 0, i = 0, k} is a
characteristic parameter of the kinetic model to be
interpreted as test particle velocity. In customary
LBM’s the parameter cs =
(with D the dimen-
sion of the set Ω) is interpreted as sound speed of
the fluid. In order that the momenta (14) and (15)
recover (in some suitable approximate sense) INSE
, however, appropriate subsidiary conditions must
be met.
2. LB assumption #2: Constraints and asymptotic
conditions: these are based on the introduc-
tion of a dimensionless parameter ε, to be consid-
ered infinitesimal, in terms of which all relevant
parameters can be ordered. In particular, it is
required that the following asymptotic orderings
[17, 40, 41] apply respectively to the fluid fields
ρo,V(r, t), p(r, t), the kinematic viscosity ν = µ/ρo
and Reynolds number Re = LV/ν:
ρo,V(r, t), p(r, t) ∼ o(ε
0), (16)
[1 + o(ε)] ∼ o(εαR), (17)
Re ∼ 1/o(ε
αR), (18)
where αR ≥ 0. Here we stress that the position
for ν holds in the case of D2Q9 only, while the
generalization to 3D and other LB discretizations.
is straightforward. Furthermore, the velocity c and
collision frequency νc are ordered so that
c ∼ 1/o(εαc), (19)
νc ∼ 1/o(ε
αν ), (20)
∼ o(εα), (21)
with α ≡ αν−αc > 0; the characteristic length and
time scales, Lo ≡ c∆t and ∆t for the spatial and
time discretization are assumed to scale as
∼ o(εαL), (22)
∼ o(εαt), (23)
with αt, αL > 0. Here L and T are the (smallest)
characteristic length and time scales, respectively
for spatial and time variations of V(r, t) and p(r.t).
Imposing also that 1
results infinitesimal at least
of order
∼ o(εα)
it follows that it must be also αt − αL > 0. These
assumptions imply necessarily that the dimension-
less parameter M eff ≡ V
(Mach number) must be
ordered as
M eff ∼ O(εαc) (24)
(small Mach-number expansion).
3. LB assumption #3: Chapman-Enskog expansion -
Kinetic initial conditions, relaxation conditions: it
is assumed that the kinetic distribution function
fi(r, t) admits a convergent Chapman-Enskog ex-
pansion of the form
fi = f
i + δf
i + δ
i + .., (25)
where δ ≡ εα and the functions f
i (j ∈ N)
are assumed smooth functions of the form (multi-
scale expansion) f
i (ro, r1, r2, ..to, t1, t2, ..), where
rn = δ
nr, tn = δ
nt and n ∈ N. In typical LBM’s
the parameter δ is usually identified with ε (which
requires letting α = 1), while the Chapman-Enskog
expansion is usually required to hold at least up to
order o(δ2). In addition the initial conditions
fi(r, to) = f
i (r, to), (26)
(for i = 0, k) are imposed in the closure of the fluid
domain Ω. It is well known [46] that this position
generally (i.e., for non-stationary fluid fields), im-
plies the violation of the Chapman-Enskog expan-
sion close to t = to, since the approximate fluid
equations are recovered only letting δf
0, i.e., assuming that the kinetic distribution func-
tion has relaxed to the Chapman-Enskog form (25).
This implies a numerical error (in the evaluation of
the correct fluid fields) which can be overcome only
discarding the first few time steps in the numerical
simulation.
4. LB assumption #5: Equilibrium kinetic distribu-
tion: a possible realization for the equilibrium dis-
tributions f
i (i = 0, k) is given by a polynomial of
second degree in the fluid velocity [44]
i (r, t) = wi
[p− Φ(r)] + (27)
+wiρo
ai ·V
ai ·V
Here, without loss of generality, the case of the
D2Q9 LB discretization will be considered, with wi
and ai (for i = 0, 8) denoting prescribed dimension-
less constant weights and discrete velocities. Notice
that, by definition, f
i is not a Galilei scalar. Nev-
ertheless, it can be considered approximately in-
variant, at least with respect to low-velocity trans-
lations which do not violate the low-Mach number
assumption (24).
5. LB assumption #6: Kinetic boundary conditions:
They are specified by suitably prescribing the form
of the incoming distribution function at the bound-
ary δΩ. [47, 48, 49, 50, 51, 51, 52, 53, 54, 54, 55, 56,
57, 58, 59]. However, this position is not generally
consistent with the Chapman-Enskog solution (25)
(see related discussion in Appendix A). As a con-
sequence violations of the hydrodynamic equations
may be expected sufficiently close to the boundary,
a fact which may be only alleviated (but not com-
pletely eliminated) by adopting suitable grid refine-
ments near the boundary. An additional potential
difficulty is related to the condition of strict posi-
tivity of the kinetic distribution function [57] which
is not easily incorporated into the no-slip boundary
conditions [50, 51, 52].
3B - Computational complexity of asymptotic
LBM’s
The requirements posed by the validity of these hy-
potheses may strongly influence the computational com-
plexity of asymptotic LBM’s which is usually associated
to the total number of ”logical” operations which must
be performed during a prescribed time interval. There-
fore, a critical parameter of numerical simulation meth-
ods is their discretization time scale ∆t. This is - in turn
- related to the Courant number NC =
, where V
and Lo.denote respectively the sup of the magnitude of
the fluid velocity and the amplitudes of the spatial dis-
cretization. As is well known ”optimal” CFD simulation
methods typically allow Lo ∼ L and a definition of the
time step ∆t = ∆tOpt such that NC ∼
V ∆tOpt
∼ 1. In-
stead, for usual LBM’s satisfying the low-M eff assump-
tion (24), the Courant number is very small since it re-
sults NC = M
eff Lo
∼ O(εα)Lo
. This means that their
discretization time scale of ∆t is much smaller than ∆tOpt
and reads
∆t ∼M eff
∆tOpt. (28)
In addition, depending on the accuracy of the numeri-
cal algorithms adopted for the construction of the dis-
crete kinetic distribution function, also the ratio Lo
sults infinitesimal in the sense Lo
∼ o(εαL), with suitable
αL > 0. Finally, we stress that LB approaches based
on the adoption of the finite-difference streaming opera-
tor LFD(i) are usually only accurate to order o(∆t
2). For
them, therefore, the requirement placed by Eq.(28) might
be even stronger. This implies that traditional LBM’s
may involve a vastly larger computation time than that
afforded by more efficient numerical methods.
4 - NEW LB INVERSE KINETIC THEORY
(LB-IKT)
A basic issue in LB approaches [8, 11, 13, 42] con-
cerns the choice of the functional class of the discrete
kinetic distribution functions fi (i = 0, k) as well as the
related definition of the equilibrium discrete distribution
function f
i [which appears in the BGK collision opera-
tor; see Eq.(11)]. This refers in particular to their trans-
formation properties with respect to arbitrary Galilean
transformations, and specifically to their Galilei invari-
ance with respect to velocity translations with constant
velocity.
In statistical mechanics it is well known that the ki-
netic distribution function is usually assumed to be a
Galilean scalar. The same assumption can, in principle,
be adopted also for LB models. However, the kinetic
distribution functions fi and f
i do not necessarily re-
quire a physical interpretation of this type. In the se-
quel we show that for a discrete inverse kinetic theory it
is sufficient that fi and f
i be so defined that the mo-
ment equations coincide with the fluid equations (which
by definition are Galilei covariant). It is sufficient to de-
mand that both fi and f
i are identified with a ordinary
scalars with respect to the group of rotation in R2, while
they need not be necessarily invariant with respect to
arbitrary velocity translations. This means that fi is in-
variant only for a particular subset of inertial reference
frames. For example for a fluid which at the initial time
moves locally with constant velocity an element of this
set can be identified with the inertial frame which in the
same position is locally co-moving with the fluid.
The adoption of non-translationally invariant discrete
distributions fi is actually already well known in LBM
and results convenient for its simplicity. This means,
manifestly, that in general no obvious physical interpre-
tation can be attached to the other momenta of the dis-
crete kinetic distribution function. As a consequence,
the very definition of the concept of statistical entropy to
be associated to the f ′is is essentially arbitrary, as well as
the related principle of entropy maximization, typically
used for the determination of the equilibrium distribution
function f
i . Several authors, nevertheless, have investi-
gated the adoption of possible alternative formulations,
which are based on suitable definitions of the entropy
functional and/or the requirement of approximate or ex-
act Galilei invariance (see for example [29, 32, 62]).
4A - Foundations of LB-IKT
As previously indicated, there are several important
motivations for seeking an exact solver based on LBM.
The lack of a theory of this type represents in fact a
weak point of LB theory. Besides being a still unsolved
theoretical issue, the problem is relevant in order to de-
termine the exact relationship between the LBM’s and
traditional CFD schemes based on the direct discretiza-
tion of the Navier–Stokes equations. Following ideas re-
cently developed [2, 3, 4, 5, 7], we show that such a theory
can be formulated by means of an inverse kinetic theory
(IKT) with discrete velocities. By definition such an IKT
should yield exactly the complete set of fluid equations
and which, contrary to customary kinetic approaches in
CFD (in particular LB methods), should not depend on
asymptotic parameters. This implies that the inverse ki-
netic theory must also satisfy an exact closure condition.
As a further condition, we require that the fluid equa-
tions are fulfilled independently of the initial conditions
for the kinetic distribution function (to be properly set)
and should hold for arbitrary fluid fields. The latter re-
quirement is necessary since we must expect that the
validity of the inverse kinetic theory should not be lim-
ited to a subset of possible fluid motions nor depend on
special assumptions, like a prescribed range of Reynolds
numbers. In principle a phase-space theory, yielding an
inverse kinetic theory, may be conveniently set in terms of
a quasi-probability, denoted as kinetic distribution func-
tion, f(x, t). A particular case of interest (investigated in
Refs.[2, 3]) refers to the case in which f(x, t) can actu-
ally be identified with a phase-space probability density.
In the sequel we address both cases, showing that, to a
certain extent, in both cases the formulation of a generic
IKT can actually be treated in a similar fashion. This
requires the introduction of an appropriate set of consti-
tutive assumptions (or axioms). These concern in par-
ticular the definitions of the kinetic equation - denoted
as inverse kinetic equation (IKE) - which advances in
time f(x, t) and of the velocity momenta to be identified
with the relevant fluid fields (correspondence principle).
However, further assumptions, such as those involving
the regularity conditions for f(x, t) and the prescription
of its initial and boundary conditions must clearly be
added. The concept [of IKT] can be easily extended to
the case in which the kinetic distribution function takes
on only discrete values in velocity space. In the sequel
we consider for definiteness the case of the so-called LB
discretization, whereby - for each (r, t) ∈ Ω × I - the
kinetic distribution function is discrete, and in particu-
lar admits a finite set of discrete values fi(r, t) ∈ R, for
i = 0, k, each one corresponding to a prescribed constant
discrete velocity ai ∈ R
3 for i = 0, k.
4B - Constitutive assumptions
Let us now introduce the constitutive assumptions (ax-
ioms) set for the construction of a LB-IKT for INSE,
whose form is suggested by the analogous continuous
inverse kinetic theory [2, 3]. The axioms, define the
”generic” form of the discrete kinetic equation, its func-
tional setting, the momenta of the kinetic distribution
function and their initial and boundary conditions, are
the following ones:
Axiom I - LB–IKE and functional setting.
Let us require that the extended fluid fields {V,p1}
are strong solutions of INSE, with initial and boundary
conditions (7)-(8) and that the pseudo pressure po(t) is
an arbitrary, suitably smooth, real function. In particu-
lar we impose that the fluid fields and the volume force
belong to the minimal functional setting:
p1,ΦǫC
(2,1)(Ω× I),
VǫC(3,1)(Ω× I), (29)
(1,0)(Ω× I).
We assume that in the set Ω×I the following equation
LD(i)fi = Ωi(fi) + Si (30)
[LB inverse kinetic equation (LB-IKE)] is satisfied iden-
tically by the discrete kinetic distributions fi(r, t) for
i = 0, k. Here Ωi(fi) and LD(i) are respectively the BGK
and the differential streaming and operators [Eqs.(11)
and (12)], while Si is a source term to be defined. We
require that KB-IKE is defined in the set Ω× I, so that
Ωi(fi) and Si are at least that C
(1)(Ω × I) and contin-
uous in Ω × I. Moreover Ωi(fi), defined by Eq.(11), is
considered for generality and will be useful for compar-
isons with customary LB approaches. We remark that
the choice of the equilibrium kinetic distribution f
the BGK operator remains completely arbitrary. We
assume furthermore that in terms of fi the fluid fields
{V, p1} are determined by means of functionals of the
form MXj [fi] =
i=0,8
Xjfi (denoted as discrete velocity
momenta). For X = X1, X2 (with X1 = c
2, X2 =
these are related to the fluid fields by means of the equa-
tions (correspondence principle)
p1(r, t)− Φ(r) = c
i=0,8
fi = c
i=0,8
i , (31)
V(r,t)=
i=1,8
aifi =
i=1,8
i , (32)
where c = min {|ai| , i = 1, 8} is the test particle veloc-
ity and f
i is defined by Eq.(27) but with the kinetic
pressure p1 that replaces the fluid pressure p adopted
previously [44]. These equations are assumed to hold
identically in the set Ω × I and by assumption, fi and
i belong to the same functional class of real functions
defined so that the extended fluid fields belong to the
minimal functional setting (29). Moreover, without loss
of generality, we consider the D2Q9 LB discretization.
Axiom II - Kinetic initial and boundary conditions.
The discrete kinetic distribution function satisfies, for
i = 0, k and for all r belonging to the closure Ω, the
initial conditions
fi(r, to) = foi(r,to) (33)
where foi(r,to) (for i = 0, k) is a initial distribution func-
tion defined in such a way to satisfy in the same set the
initial conditions for the fluid fields
p1o(r) ≡ Po(to) + po(r)− Φ(r) = (34)
i=0,8
foi(r),
Vo(r) =
i=1,8
aifoi(r) . (35)
To define the analogous kinetic boundary conditions on
δΩ, let us assume that δΩ is a smooth, possibly moving,
surface. Let us introduce the velocity of the point of the
boundary determined by the position vector rw ∈ δΩ, de-
fined by Vw(rw(t), t) =
rw(t) and denote by n(rw, t)
the outward normal unit vector, orthogonal to the bound-
ary δΩ at the point rw. Let us denote by f
i (rw, t)
and f
i (rw , t) the kinetic distributions which carry the
discrete velocities ai for which there results respectively
(ai −Vw) ·n(rw , t) > 0 (outgoing-velocity distributions)
and (ai −Vw) · n(rw, t) ≤ 0 (incoming-velocity distribu-
tions) and which are identically zero otherwise. We as-
sume for definiteness that both sets, for which |ai| > 0,
are non empty (which requires that the parameter c
be suitably defined so that c > |Vw|). The bound-
ary conditions are obtained by prescribing the incom-
ing kinetic distribution f
i (rw , t), i.e., imposing (for all
(rw, t) ∈ δΩ× I)
i (rw, t) = f
oi (rw , t). (36)
Here f
oi (rw, t) are suitable functions, to be assumed
non-vanishing and defined only for incoming discrete ve-
locities for which (ai −Vw)·n(rw , t) ≤ 0. Manifestly, the
functions f
oi (rw, t) (i = 0, k) must be defined so that
the Dirichlet boundary conditions for the fluid fields are
identically fulfilled, namely there results
p1w(rw, t) = Po(t) + pw(rw, t)− Φ(r) = (37)
i=0,k
oi (rw, t) + f
i (rw, t)
Vw(rw, t) = (38)
i=1,k
oi (rw, t) + f
i (rw, t)
Here, again, the functions foi(r) and f
oi (rw, t) (for i =
0, k) must be assumed suitably smooth. A particular case
is obtained imposing identically for i = 0, k
foi(r,to) = f
i (r, to), (39)
oi (rw, t) = f
i (rw , t), (40)
where the identification with f
oi (rw, t) and f
oi (rw, t)
is intended respectively in the subsets ai ·n(rw, t) > 0 and
ai ·n(rw , t) ≤ 0. Finally, we notice that in case Neumann
boundary conditions are imposed on the fluid pressure,
Eq.(37) still holds provided pw(rw, t) is intended as a
calculated value.
Axiom III - Moment equations.
If fi(r, t), for i = 0, k, are arbitrary solutions of LB-
IKE [Eq.(30)] which satisfy Axioms I and II validity of
Axioms I and II, we assume that the moment equations of
the same LB-IKE, evaluated in terms of the moment op-
erators MXj [·] =
i=0,8
Xj ·, with j = 1, 2, coincide iden-
tically with INSE, namely that there results identically
[for all (r, t) ∈ Ω× I]
MX1 [Lifi − Ωi(fi)− Si] = ∇ ·V = 0, (41)
MX2 [Lifi − Ωi(fi)− Si] = NV = 0. (42)
Axiom IV - Source term.
The source term is required to depend on a finite num-
ber of momenta of the distribution function. It is as-
sumed that these include, at most, the extended fluid
fields {V,p1} and the kinetic tensor pressure
Π = 3
fiaiai − ρoVV. (43)
• Furthermore, we also normally require (except
for the LB-IKT described in Appendix B) that
Si(r, t) results independent of f
i (r,t), foi(r) and
fwi(rw , t) (for i = 0, k).
Although, the implications will made clear in the fol-
lowing sections, it is manifest that these axioms do not
specify uniquely the form (and functional class) of the
equilibrium kinetic distribution function f
i (r,t), nor
of the initial and boundary kinetic distribution func-
tions (33),(36). Thus, both f
i (r,t), foi(r,to) and the
related distribution they still remain in principle com-
pletely arbitrary. Nevertheless, by construction, the
initial and (Dirichlet) boundary conditions for the fluid
fields are satisfied identically. In the sequel we show that
these axioms define a (non-empty) family of parameter-
dependent LB-IKT’s, depending on two constant free pa-
rameters νc, c > 0 and one arbitrary real function Po(t).
The examples considered are reported respectively in the
following Sec. 5,6 and in the Appendix B.
5 - A POSSIBLE REALIZATION: THE
INTEGRAL LB-IKT
We now show that, for arbitrary choices of the distri-
butions fi(r,t) and f
i (r,t) which fulfill axioms I-IV, an
explicit (and non-unique) realization of the LB-IKT can
actually be obtained. We prove, in particular, that a pos-
sible realization of the discrete inverse kinetic theory, to
be denoted as integral LB-IKT, is provided by the source
Si = (44)
− ai ·
f1−µ∇
V −∇ ·Π+∇p
≡ S̃i,
where wi
is denoted as first pressure term. Holds, in
fact, the following theorem.
Theorem 1 - Integral LB-IKT
In validity of axioms I-IV the following statements
hold. For an arbitrary particular solution fi and for ar-
bitrary extended fluid fields :
A) if fi is a solution of LB-IKE [Eq.(30)] the moment
equations coincide identically with INSE in the set Ω×I;
B) the initial conditions and the (Dirichlet) boundary
conditions for the fluid fields are satisfied identically;
C) in validity of axiom IV the source term S̃i is non-
uniquely defined by Eq.(44).
Proof
A) We notice that by definition there results identically
S̃i =
aiS̃i = (46)
f−µ∇2V−∇ ·Π+∇p
On the other hand, by construction (Axiom I) fi (i =
1, k) is defined so that there results identically
i=0 Ωi =
0 and
i=0 aiΩi = 0. Hence the momenta MX1 ,MX2 of
LB-IKE deliver respectively
i=1,8
aifi = 0 (47)
i=1,8
aifi + ρoV · ∇V +∇p1 + f−µ∇
V = 0 (48)
where the fluid fields V,p1 are defined by Eqs.(31),(32).
Hence Eqs.(47) and (48) coincide respectively with the
isochoricity and Navier-Stokes equations [(1) and (2)].
As a consequence, fi is a particular solution of LB-IKE
iff the fluid fields {V,p1} are strong solutions of INSE.
B) Initial and boundary conditions for the fluid fields
are satisfied identically by construction thanks to Axiom
C) However, even prescribing νc, c > 0 and the real
function Po(t), the functional form of the equation can-
not be unique The non uniqueness of the functional form
of the source term S̃i(r, t) is assumed to be indepen-
dent of f
i (r,t) [and hence of Eq.(30)] is obvious. In
fact, let us assume that S̃i is a particular solution for
the source term which satisfies the previous axioms I-
IV. Then, it is always possible to add to Si arbitrary
terms of the form S̃i + δSi, with δSi 6= 0 which depends
only on the momenta indicated above, and gives van-
ishing contributions to the first two moment equations,
namely MXj [δSi] =
i=0,8
XjδSi = 0, with j = 1, 2. To
prove the non-uniqueness of the source term Si, it is suf-
ficient to notice that, for example, any term of the form
δSi =
F (r, t), with F (r, t) an arbitrary real
function (to be assumed, thanks to Axiom IV, a linear
function of the fluid velocity), gives vanishing contribu-
tions to the momentaMX1 ,MX2 . Hence S̃i is non-unique.
The implications of the theorem are straightforward.
First, manifestly, it holds also in the case in which the
BGK operator vanishes identically. This occurs letting
νc = 0 in the whole domain Ω × I. Hence the inverse
kinetic equation holds independently of the specific defi-
nition of f
i (r,t).
An interesting feature of the present approach lies in
the choice of the boundary condition adopted for fi(r,t),
which is different from that usually adopted in LBM’s
[see for example [14] for a review on the subject]. In par-
ticular, the choice adopted is the simplest permitting to
fulfill the Dirichlet boundary conditions [imposed on the
fluid fields]. This is obtained prescribing the functional
form of fi(r,t) on the boundary of the fluid domain (δΩ),
which is identified with a function foi(r, t).
Second, the functional class of fi(r,t), f
i (r,t) and of
foi(r, t) remains essentially arbitrary. Thus, in particu-
lar, the initial and boundary conditions, specified by the
same function foi(r, t), can be defined imposing the po-
sitions (39),(40). As further basic consequence, f
i (r,t)
and fi(r,t) need not necessarily be Galilei-invariant (in
particular they may not be invariant with respect to ve-
locity translations), although the fluid equations must
be necessarily fully Galilei-covariant. As a consequence
it is always possible to select f
i (r,t) and foi(r, t) based
on convenience and mathematical simplicity. Thus, be-
sides distributions which are Galilei invariant and sat-
isfy a principle of maximum entropy (see for example
[22, 30, 32, 34, 60, 61]), it is always possible to iden-
tify them [i.e., f
i (r,t), foi(r, t)] with a non-Galilean in-
variant polynomial distribution of the type (27) [mani-
festly, to be exactly Galilei-invariant each f
i (r,t) should
depend on velocity only via the relative velocity ui =
ai −V].
We mention that the non-uniqueness of the source term
S̃i can be exploited also by imposing that f
i (r,t) re-
sults a particular solution of the inverse kinetic equation
Eq.(30) and there results also foi(r, t) = f
i (r,t). In Ap-
pendix B we report the extension of THM.1 which is ob-
tained by identifying again f
i (r,t) with the polynomial
distribution (27).
6 - THE ENTROPIC PRINCIPLE - CONDITION
OF POSITIVITY OF THE KINETIC
DISTRIBUTION FUNCTION
A fundamental limitation of the standard LB ap-
proaches is their difficulty to attain low viscosities, due to
the appearance of numerical instabilities [14]. In numeri-
cal simulations based on customary LB approaches large
Reynolds numbers is usually achieved by increasing nu-
merical accuracy, in particular strongly reducing the time
step and the grid size of the spatial discretization (both
of which can be realized by means of numerical schemes
with adaptive time-step and using grid refinements).
Hence, the control [and possible inhibition] of numerical
instabilities is achieved at the expense of computational
efficiency. This obstacle is only partially alleviated by
approaches based on ELBM [22, 30, 32, 34, 60, 61]. Such
methods are based on the hypothesis of fulfilling an H-
theorem, i.e., of satisfying in the whole domain Ω × I
the condition of strict positivity for the discrete kinetic
distribution functions. This requirement is considered,
by several authors (see for example [26, 29, 62]), an es-
sential prerequisite to achieve numerical stability in LB
simulations. However, the numerical implementation of
ELBM typically induce a substantial complication of the
original algorithm, or require a cumbersome fine-tuning
of adjustable parameters [22, 37].
6A - The constant entropy principle and PEM
A basic aspect of the IKT’s here developed is the possi-
bility of fulfilling identically the strict positivity require-
ment by means of a suitable H-theorem which provides
also a maximum entropy principle. In particular, in this
Section, extending the results of THM.1 and 2, we intend
to prove that a constant H-theorem can be established
both for the integral and differential LB-IKT’s defined
above. The H-theorem can be reached by imposing for
the Gibbs-Shannon entropy functional the requirement
that for all t ∈ I there results
S(f) = −
i=0,8
fi ln(fi/wi) = 0, (49)
which implies that S(f) is necessarily maximal in a suit-
able functional set {f} . The result can be stated as fol-
lows:
Theorem 2 - Constant H-theorem
In validity of THM.1, let us assume that:
1) the configuration domain Ω is bounded;
2) at time to the discrete kinetic distribution functions
fi, for i = 0, 8, are all strictly positive in the set Ω.
Then the following statements hold:
A) by suitable definition of the pseudo pressure
Po(t), the Gibbs-Shannon entropy functional S(f) =
i=0,8
fi ln(fi/wi) can be set to be constant in the
whole time interval I. This holds provided the pseudo-
pressure Po(t) satisfies the differential equation
(1 + log fi) = (50)
ai · ∇fi − Ŝi
(1 + log fi) ,
where Ŝi = Si +
B) if the entropy functional S(f) =
i=0,8
fi ln(fi/wi) is constant in the whole
time interval I the discrete kinetic distribution functions
fi are all strictly positive in the whole set Ω× I;
C) an arbitrary solution of LB-IKE [Eq.(30)] which
satisfies the requirement A) is extremal in a suitable func-
tional class and maximizes the Gibbs-Shannon entropy .
Proof:
A) Invoking Eq.(30), there results
∂S(t)
[1 + log fi] = (51)
(ai · ∇fi − Si) (1 + log fi) ,
where Si is the source term, provided by Eq.(44). By
direct substitution it follows the thesis.
B) If Eq.(50) holds identically in there results ∀t ∈
I, S (t) = S (t0) , which implies the strict positivity of fi,
for all i = 0, 8.
C) Let us introduce the functional class
{f + αδf} = {fi = fi(t) + αδfi(t), i = 0, 8} , (52)
where α is a finite real parameter and the syn-
chronous variation δfi(t) is defined δfi(t) = dfi(t) ≡
∂fi(t)
dt. Introducing the synchronous variation of the en-
tropy, defined by δS (t) = ∂
, with ψ(α) =
S (f + αδf) , it follows
δS (t) = dt
∂S(t)
. (53)
Since in validity of Eq.(50) there results
∂S(t)
which in view of Eq.(53) implies also δS (t) = 0. It is im-
mediately follows that there results necessarily δ2S (t) ≤
0, i.e., S (t) is maximal. Therefore, the kinetic distribu-
tion function which satisfies IKE (Eq.(30)] is extremal in
the functional class of variations (52) and maximizes the
Gibbs-Shannon entropy functional.
6B - Implications
In view of statement B, THM.2 warrants the strict pos-
itivity of the discrete distribution functions fi (i = 0, 8)
only in the open set Ω × I, while nothing can be said
regarding their behavior on the boundary δΩ (on which
fi might locally vanish). However, since the inverse ki-
netic equation actually holds only in the open set Ω× I,
this does not affect the validity of the result. While the
precise cause of the numerical instability of LBM’s is still
unknown,the strict positivity of the distribution function
is usually considered important for the stability of the nu-
merical solution [29, 30]. It must be stressed that the nu-
merical implementation of the condition of constant en-
tropy Eq.(50) should be straightforward, without involv-
ing a significant computational overhead for LB simula-
tions. Therefore it might represent a convenient scheme
to be adopted also for customary LB methods.
7 - ASYMPTOTIC APPROXIMATIONS AND
COMPARISONS WITH PREVIOUS CFD
METHODS
A basic issue is the relationship with previous CFD nu-
merical methods, particularly asymptotic LBM’s. Here
we consider, for definiteness, only the case of the inte-
gral LB-IKT introduced in Sec.5. Another motivation is
the possibility of constructing new improved asymptotic
models, which satisfy with prescribed accuracy the re-
quired fluid equations [INSE], of extending the range of
validity of traditional LBM’s and fulfilling also the en-
tropic principle (see Sec.6). The analysis is useful in
particular to establish on rigorous grounds the consis-
tency of previous LBM’s. The connection [with previ-
ous LBM’s] can be reached by introducing appropriate
asymptotic approximations for the IKT’s, obtained by
assuming that suitable parameters which characterize the
IKT’s are infinitesimal (or infinite) (asymptotic parame-
ters). A further interesting feature is the possibility of
constructing in principle a class of new asymptotic LBM’s
with prescribed accuracy , i.e., in which the distribution
function (and the corresponding momenta) can be de-
termined with predetermined accuracy in terms of per-
turbative expansions in the relevant asymptotic parame-
ters. Besides recovering the traditional low-Mach number
LBM’s [17, 21, 40], which satisfy the isochoricity condi-
tion only in an asymptotic sense and are closely related to
the Chorin artificial compressibility method, it is possible
to obtain an improved asymptotic LBM’s which satisfy
exactly the same equation.
We first notice that the present IKT is characterized
by the arbitrary positive parameters νc, c and the initial
value Po(to), which enter respectively in the definition
of the BGK operator [see (11)], the velocity momenta
and equilibrium distribution function f
i . Both c and
Po(to) must be assumed strictly positive, while, to assure
the validity of THM.2, Po(to) must be defined so that
(for all i = 0, 8) f
i (r,to) > 0 in the closure Ω. Thanks
to THM.1.and 2 the new theory is manifestly valid for
arbitrary finite value of these parameters. This means
that they hold also assuming
o(εαν )
, (54)
o(εαc)
, (55)
Po(to) ∼ o(ε
0), (56)
where ε denotes a strictly positive real infinitesimal,
αν , αc > 0 are real parameters to be defined, while the
extended fluid fields {ρ,V, p1} and the volume force f
are all assumed independent of ε. Hence, with respect to
ε they scale
ρo,V,p1, f ∼ o(ε
0). (57)
As a result, for suitably smooth fluid fields (i.e., in va-
lidity of Axiom 1) and appropriate initial conditions for
fi(r, t), it is expected that the first requirement actually
implies in the whole set Ω× I the condition of closeness
fi(r, t) ∼= f
i (r, t) [1 + o(ε)] , consistent with the LB As-
sumption #4. To display meaningful comparisons with
previous LBM’s let us introduce the further assumption
that the fluid viscosity is small in the sense
µ ∼ o(εαµ), (58)
with αµ ≥ 1 another real parameter to be defined.
Asymptotic approximations for the corresponding LB-
IKE [Eq.(30)] can be directly recovered by introducing
appropriate asymptotic orderings for the contributions
appearing in the source term Si = S̃i. Direct inspec-
tion shows that these are provided by the (dimensional)
parameters
M effp,a ≡
, (59)
∣∣∇ ·Π−∇p
∣∣ , (60)
∣∣µ∇2V
∣∣ . (61)
The first two M effp,a and M
are here denoted respec-
tively as (first and second) pressure effective Mach num-
bers, driven respectively by the pressure time-derivative
and by the divergence of the pressure anisotropy Π−p1.
Furthermore, M
is denoted as velocity effective Mach
number. Physically relevant examples [of asymptotic
LBM’s] can be achieved by introducing suitable orderings
in terms of the single infinitesimal ε for the parameters
M effp,a ,M
.We stress that these orderings, in
principle, can be introduced without actually introducing
restrictions on the fluid fields, i.e., retaining the assump-
tion that the extended fluid fields are independent of ε.
Interesting cases are provided by the asymptotic order-
ings indicated below.
7A - Small effective Mach numbers (Meffp,a ,M
An important aspect of LB theory is the possibility
of constructing asymptotic LBM’s with prescribed accu-
racy with respect to the infinitesimal parameter ε, in the
sense that the fluid equations are satisfied at least cor-
rect up to terms of order o(εn) included, with n = 1
or 2, namely ignoring error terms of order o(εn+1) or
higher. Let us, first, consider the case in which all pa-
rametersM effp,a ,M
and M
are all infinitesimal w.r.
to ε (low-effective-Mach numbers). Since the parameters
c and νc are free, they can be defined so that that there
results c ∼ νc ∼ 1/o(ε) [which implies αc = αν = 1].
This requires
M effp,a ∼M
∼ o(ε2). (62)
If, we consider a low-viscosity fluid for which the kine-
matic viscosity ν = µ/ρo can be assumed of order ε [and
hence αµ = 1] it follows that
∼ o(ε2). (63)
Thanks to the assumptions (54)-(58) there follows ∇ ·
Π − ∇p ∼ o(ε) and µ∇2V ∼ o(ε),which implies that
the source term S̃i, ignoring corrections of order o(ε
becomes
S̃i ∼= S̃Ai [1 + o(ε)] , (64)
S̃Ai ≡ −
ai · f . (65)
It is immediate to determine the corresponding moment
equations, which read:
+∇ ·V = 0, (66)
NV = 0+ o(ε2), (67)
Formally the first equation can be interpreted as an evo-
lution equation for the kinetic pressure p1. Nevertheless,
in view of the ordering (62) it actually implies the iso-
choricity condition
∇ ·V = 0 + o(ε2). (68)
Instead, the second one [Eq.(67)]. due to the asymp-
totic approximation (63), reduces to the Euler equation.
Therefore in this case the asymptotic approximation (64)
is not adequate. To recover the correct Navier-Stokes
equation a more accurate approximation is needed, real-
ized requiring that the hydrodynamic equations are sat-
isfied correct to order o(ε3). A fist possibility is to con-
sider a more accurate approximation for the source term.
Restoring the pressure and viscous source terms in (64)
there results the asymptotic source term
S̃Bi ≡
− ai ·
f1−µ∇
, (69)
where in validity of the previous orderings
S̃i ∼= S̃Bi [1 + o(ε)] . (70)
The corresponding moment equations become therefore
∇ ·V = 0, (71)
NV = 0+ o(ε3). (72)
It is remarkable that in this case the isochoricity condi-
tion is exactly fulfilled, even if the source term is not the
exact one. For the sake of reference, it is interesting to
mention another possible small-Mach-number ordering.
This is obtained imposing for the parameters c and νc
, (73)
o(ε2)
, (74)
while requiring for ν = µ/ρo the same constraint adopted
by asymptotic LBM’s, namely Eq.(17). In this case one
can show that the moment equation (72) is actually satis-
fied correct to order o(ε3), while the isochoricity condition
is only satisfied to order o(ε2). The following theorem
can, in fact, be proven:
Theorem 3 - Low effective-Mach-numbers asymptotic
approximation
In validity of THM.1, let us invoke the following as-
sumptions:
1) LB assumptions #3 and #4 for the discrete kinetic
distributions fi ( i = 0, 8);
2) the free parameters c and νc are assumed to satisfy
the asymptotic orderings (73),(74);
3) the fluid viscosity µ is assumed of order µ ∼ o(ε)
4) the fluid viscosity µ is prescribed so that the kine-
matic viscosity ν = µ/ρo is defined in accordance to
Eq.(17);
5) the kinetic pressure p1 is assumed slowly varying in
the sense
∂ ln p1
∼ o(ε). (75)
It follows that the source term is approximated by
Eq.(64) and moment equations are provided by the
asymptotic equations:
+∇ ·V = 0 + o(ε3), (76)
NV = 0+ o(ε3), (77)
i.e., the isochoricity and NS equation are recovered re-
spectively correct to order o(ε2) and o(ε3).
Proof
First we notice that the ordering assumptions 2)-5)
require
M effp,a ∼ o(ε
3) (78)
∼ o(ε2), (79)
∼ o(ε4), (80)
which imply at least the validity of Eqs.(64)-(67). The
proof of Eqs.(76) and (77) is immediate. In both cases
it sufficient to notice that in validity of hypotheses 1)-3)
and in terms of a Chapman-Enskog perturbative solution
of Eq.(30) there results actually
−µ∇2V −∇ ·Π+∇p = O + o(ε3), (81)
and hence S̃i reduces to Eq.(64).
The predictions of THM.3 are relevant for comparisons
and to provide asymptotic accuracy estimates for previ-
ous asymptotic LBM’s [see Refs. [17, 21, 40]]. In fact,
the asymptotic moment equations (76) and (77) formally
coincide with the analogous moment equations predicted
by such theories, when the kinetic pressure p is replaced
by the fluid pressure p1 (i.e., if the function Po(t) is set
identically equal to zero). [17, 21, 40]. Nevertheless, the
accuracy of customary LBM’s depends on the properties
of the solutions of INSE. In fact, if one assumes
∂ ln p1
∼ o(ε0) (82)
the customary (V, p) asymptotic LBM [17, 21, 40] result
actually accurate only to order o(ε2). Therefore, in such
case to reach an accuracy of order o(ε3) the approxima-
tion (69) must be invoked for the source term.
The other interesting feature of Eqs.(76) and (77) is
that they provide a connection with the artificial com-
pressibility method (ACM) postulated by Chorin [63],
previously motivated merely on the grounds of an asymp-
totic LBM [21]. In fact, these coincides with the Chorin’s
pressure relaxation equation where c can be interpreted
as sound speed of the fluid. However - in a sense - this
analogy is purely formal and is only due to the neglect of
the first pressure source term in Si. It disappears alto-
gether in Eq.(71) if we adopt the more accurate asymp-
totic source term (69). A further difference is provided
by the adoption of the kinetic pressure p1 which replaces
the fluid pressure p (used in Chorin approach). We stress
that the choice of p1 here adopted, with Po(t) determined
by the entropic principle, represents an important differ-
ence, since it permits to satisfy everywhere in Ω× I the
condition of strict positivity for the discrete kinetic dis-
tribution functions.
7B - Finite pressure-Mach number Meffp,a
Another possible asymptotic ordering, usually not per-
mitted by customary asymptotic LBM’s, is the one in
which the test particle velocity is finite, namely c ∼ o(ε0),
the viscosity remains arbitrary and is taken of order
µ ∼ o(ε0) while again νc is assumed νc ∼ 1/o(ε
2) [i.e.,
αc = νc = 0, αν = 2]. In this case the pressure Mach
M effp,a number results finite, while velocity and the sec-
ond pressure Mach numbers are considered infinitesimal,
respectively of first and second order in ε, namely
M effp,a ∼ o(ε
∼ o(ε), (83)
∼ o(ε2).
To obtain the fluid equation with the prescribed accu-
racy, say of order o(ε2), it is sufficient to approximate
the source term S̃i in terms of S̃i ∼= S
1 + o(ε2)
. The
set of asymptotic moment equations coincide therefore
with Eqs.(71),(72). Again, the isochoricity condition is
exactly fulfilled, while in this case the NS equation is
accurate only to order o(ε2).
7C - Small effective pressure-Mach numbers
(Meffp,a ,M
) and finite velocity-Mach number (M
Finally, another interesting case is the one in which the
fluid viscosity µ remains finite (strongly viscous fluid),
i.e., in the sense µ ∼ o(ε0) [i.e., αµ = 0] while both
parameters c and νc are suitably large, and respectively
scale as c ∼ 1/o(ε), νc ∼ 1/o(ε
2) [i.e., αc = 1, αν = 2].
Due to assumptions (54)-(58) one obtains ∇ ·Π−∇p ∼
o(ε2) and µ∇2V ∼ o(ε0). It follows that the effective
Mach numbers scale respectively as
∼ o(ε3) (84)
M effp,a ∼ M
∼ o(ε2),
If we impose on µ also the same constraint set by Eq.(17),
the customary asymptotic LBM’s can be invoked also in
this case. However, since the first pressure and veloc-
ity Mach numbers are only second order accurate, the
NS equation is recovered to order o(ε2) only. Never-
theless, it is possible to recover with prescribed accuracy
the fluid equations (71),(72). This is obtained adopt-
ing the source term S̃i ∼= S̃Bi [see Eq.(69)]. As a basic
consequence, the isochoricity equation is satisfied exactly
(hence no meaningful analogy with Chorin’s approach
arises), while the NS equation results correct to order
o(ε3). These results provide a meaningful extension of
the customary asymptotic LBM’s. We stress that the
entropic approach here developed holds independently of
the asymptotic orderings here considered [for the param-
etersM effp,a ,M
]. Thus it can be used in all cases
to assure the strict positivity of the discrete distribution
function.
8 - CONCLUSIONS
In this paper we have presented the theoretical foun-
dations of a new phase-space model for incompressible
isothermal fluids, based on a generalization of customary
lattice Boltzmann approaches.We have shown that many
of the limitations of traditional (asymptotic) LBM’s can
be overcome. As a main result, we have proven that
the LB-IKT can be developed in such a way that it
furnishes exact Navier-Stokes and Poisson solvers, i.e.,
it is - in a proper sense - an inverse kinetic theory for
INSE. The theory exhibits several features, in particular
we have proven that the integral LB-IKT (see Sec.5):
1. determines uniquely the fluid pressure p(r, t) via
the discrete kinetic distribution function without
solving explicitly (i.e., numerically) the Poisson
equation for the fluid pressure. Although analo-
gous to traditional LBM’s, this is interesting since
it is achieved without introducing compressibility
and/or thermal effects. In particular the present
theory does not rely on a state equation for the
fluid pressure.
2. is complete, namely all fluid fields are expressed as
momenta of the distribution function and all hy-
drodynamic equations are identified with suitable
moment equations of the LB inverse kinetic equa-
tion.
3. allows arbitrary initial and boundary conditions for
the fluid fields.
4. is self-consistent : the kinetic theory holds for ar-
bitrary, suitably smooth initial conditions for the
kinetic distribution function. In other words, the
initial kinetic distribution function must remain ar-
bitrary even if a suitable set of its momenta are
prescribed at the initial time.
5. the associated the kinetic and equilibrium distri-
bution functions can always be chosen to belong to
the class of non-Galilei-invariant distributions. In
particular the equilibrium kinetic distribution can
always be identified with a polynomial of second
degree in the velocity.
6. is non-asymptotic, i.e., unlike traditional LBM’s it
does not depend on any small parameter, in partic-
ular it holds for finite Mach numbers.
7. fulfills an entropic principle, based on a constant-H
theorem. This theorem assures, at the same time,
the strict positivity of the discrete kinetic distri-
bution function and the maximization of the as-
sociated Gibbs-Shannon entropy in a properly de-
fined functional class. Remarkably the constant H-
theorem is fulfilled for arbitrary (strictly positive)
kinetic equilibria. This includes also the case of
polynomial kinetic equilibria.
A further remarkable aspect of the theory concerns the
choice of the kinetic boundary conditions to be satisfied
by the distribution function (Axiom II) and obtained by
prescribing the form of the incoming-velocity distribution
[see Eq.(36)]. Thanks to Eqs.(34),(35), this requirement
[of the LB-IKT] the boundary conditions for the fluid
fields are satisfied exactly while the fluid equations are
by construction identically fulfilled also arbitrarily close
to the boundary. This result, in a proper sense, applies
only to Dirichlet boundary conditions for the fluid fields
[see Eqs.(8)]. Nevertheless the same approach can be
in principle extended to the case of mixed or Neumann
boundary conditions for the fluid fields.
Moreover, we have shown that a useful implication of
the theory is provided by the possibility of constructing
asymptotic approximations to the inverse kinetic equa-
tion. This permits to develop a new class of asymptotic
LBM’s which satisfy INSE with prescribed accuracy, to
obtain useful comparisons with previous CFD methods
(Chorin’s ACM) and to achieve accuracy estimates for
customary asymptotic LBM’s. The main results of the
paper are represented by THM’s 1-3, which refer respec-
tively to the construction of the integral LB-IKT, to the
entropic principle and to construction of the low effective-
Mach-numbers asymptotic approximations. For the sake
of reference, also another type of LB-IKT, which admits
as exact particular solution the polynomial kinetic equi-
librium, has been pointed out (THM.1bis).
The construction of a discrete inverse kinetic theory of
this type for the incompressible Navier-Stokes equations
represents an exciting development for the phase-space
description of fluid dynamics, providing a new starting
point for theoretical and numerical investigations based
on LB theory. In our view, the route to more accurate,
higher-order LBM’s, here pointed out, will be important
in order to achieve substantial improvements in the effi-
ciency of LBM’s in the near future.
APPENDIX A
The basic argument regarding the accuracy of the
boundary conditions adopted by customary asymptotic
LBM’s is provided by Ref.[46]. In fact. let us assume that
on the boundary δΩ the incoming distribution function
i (rw , t) is prescribed according to Eqs.(33),(37) and
(38), being f
oi (rw, t) prescribed suitably smooth func-
tions which are non vanishing only only for incoming dis-
crete velocities ai for which (ai −Vw) ·n(rw, t) ≤ 0. For
definiteness, let us assume that f
oi (rw , t) ≡ f
i (rw, t)
where f
i (rw, t) denotes a suitable equilibrium distribu-
tion. It follows that suitably close to the boundary the
kinetic distribution differs from the Chapman-Enskog so-
lution (25). The numerical error can be overcome only
discarding the first few spatial grid (close to the bound-
ary) in the numerical simulation [46].
APPENDIX B
Unlike standard kinetic theory, the distinctive feature
of LB-IKT’s is the possibility of adopting a non-Galilei
invariant kinetic distribution function (i.e., non-invariant
with respect to velocity translations). Here we re-
port another example of discrete inverse kinetic theory
of this type. Let us modify Axiom IV so that to permit
that a particular solution of LB-IKE [Eq.(30)] is pro-
vided by fi = f
i . Here we identify f
i with the (non-
Galilei invariant) polynomial kinetic distribution defined
by Eq.(27) but with the kinetic pressure p1 that replace
the fluid pressure p. In this case one can prove that the
source term Si reads
Si = S
i ≡ S̃i +∆Si, (85)
where
∆Si =
(ai −V) · ∇V−
ai∇·V
· ai+
V − ai
3ai ·V
+ (86)
wiρoai · ∇
ai ·V
Here N1 ≡ N−ρo
, where N is the Navier-Stokes opera-
tor (6), namely N1 is the nonlinear operator which acting
onV yields N1V = ρoV·∇V+∇ [p1 − Φ (r)]+f1−µ∇
Hence, invoking INSE, ∆Si can also be written in the
equivalent form
∆Si =
(ai −V) · ∇V−
ai∇·V
· ai+
V − ai
3ai ·V
+ (87)
wiρoai · ∇
ai ·V
The following result holds:
Theorem 1bis - Differential LB-IKT
In validity of axioms I-IV and the assumption that
fi = f
i is a particular solution of Eq.(30), the following
statements hold:
i is a particular solution of LB-IKE [Eq.(30)]
if and only if the extended fluid fields {V,p1} are strong
solutions of INSE of class (29), with initial and boundary
conditions (7)-(8), and arbitrary pseudo pressure po(t) of
class C(1)(I).
Moreover, for an arbitrary particular solution fi and
for arbitrary extended fluid fields :
For an arbitrary particular solution fi :
B) fi is a solution of LB-IKE [Eq.(30)] if and only
if the extended fluid fields {V,p1} are arbitrary strong
solutions of INSE of class (29), with initial and boundary
conditions (7)-(8), and arbitrary pseudo pressure po(t) of
class C(1)(I);
C) the moment equations of L-B IKE coincide identi-
cally with INSE in the set Ω× I;
D) the initial conditions and the (Dirichlet) boundary
conditions for the fluid fields are satisfied identically;
E) the source term Si is uniquely defined by
Eqs.(85),(86);
Proof:
The proof of propositions A,B, C and D is analogous
to that provided in THM.1. Assuming Si = S
i , the
proof of B follows from straightforward algebra. In fact,
letting fi(r, t) = f
i (r, t) for all (r, t) ∈ Ω × I in the
LB-IKE [Eq.(30)], one finds that Eq.(30) is fulfilled iff
the fluid fields satisfy the Navier-Stokes, isochoricity and
incompressibility equations (1),(2) and (3). The proof
of proposition E can be reached in a similar way. The
uniqueness of the source term Si is an immediate conse-
quence of the uniqueness of the solutions for INSE.
ACKNOWLEDGEMENTS Useful comments and
stimulating discussions with K.R. Sreenivasan, Direc-
tor, ICTP (International Center of Theoretical Physics,
Trieste, Italy) are warmly acknowledged. Research de-
veloped in the framework of PRIN Project Fundamen-
tals of kinetic theory and applications to fluid dynamics,
magnetofluid dynamics and quantum mechanics (MIUR,
Ministry for University and Research, Italy), with the
support of the Consortium for Magnetofluid Dynamics,
Trieste, Italy.
[1] M. Ellero and M. Tessarotto, Bull. Am Phys. Soc. 45 (9),
40 (2000).
[2] M. Tessarotto and M. Ellero, RGD24 (Italy, July 10-16,
2004), AIP Conf. Proc. 762, 108 (2005).
[3] M. Ellero and M. Tessarotto, Physica A 355, 233 (2005).
[4] M. Tessarotto and M. Ellero, Physica A 373, 142 (2007);
arXiv: physics/0602140.
[5] M. Tessarotto and M. Ellero, “On the uniqueness of con-
tinuous inverse kinetic theory for incompressible fluids,”
in press on AIP Conf. Proc., RGD25 (St. Petersburg,
Russia, July 21-28, 2006); arXiv:physics/0611113.
[6] M. Tessarotto, M. Ellero, N. Aslan, M. Mond and
P. Nicolini, “An exact pressure evolution equation
for the incompressible Navier-Stokes equationsInverse”,
arXiv:physics/0612072 (2006).
[7] M. Tessarotto, M. Ellero and P. Nicolini, Phys.Rev. A
75, 012105 (2007); arXiv:quantum-ph/060691.
[8] G.R. McNamara and G. Zanetti, Phys. Rev. Lett. 61,
2332 (1988).
[9] F. Higuera, S. Succi and R. Benzi, Europhys. Lett. 9,345
(1989).
[10] S. Succi, R. Benzi, and F. Higuera, Physica D 47, 219
(1991).
[11] S. Chen, H. Chen, D. O. Martinez, and W. H. Matthaeus,
Phys.Rev. Lett. 67, 3776 (1991).
[12] R. Benzi et al., Phys. Rep. 222, 145 (1992).
[13] H. Chen, S. Chen, and W. Matthaeus, Phys. Rev. A 45,
R5339 (1992).
[14] S. Succi, The Lattice Boltzmann Equation for Fluid Dy-
namics and Beyond (Numerical Mathematics and Scien-
tific Computation), Oxford Science Publications (2001).
[15] U. Frisch, B. Hasslacher, and Y. Pomeau, Rev.Lett. 56,
1505 (1986).
[16] U. Frisch, D. d’Humieres, B. Hasslachaer, P. Lallemand,
Y. Pomeau, and J.-P. Rivet, Complex Syst. 1, 649 (1987).
[17] X. He and L-S. Luo, Phys. Rev. E 56, 6811 (1997).
[18] D. O. Martinez, W. H. Matthaeus, S. Chen, et al., Phys.
Fluids 6, 1285 (1994).
[19] S. L. Hou, Q. S. Zou, S. Y. Chen, G. Doolen, and A. C.
Cogley, J. Comput. Phys. 118, 329 (1995).
[20] X. He and G. Doolen, J. Comput. Phys. 134, 306 (1997).
http://arxiv.org/abs/physics/0602140
http://arxiv.org/abs/physics/0611113
http://arxiv.org/abs/physics/0612072
[21] X. He, G.D. Doolen, and T. Clark, J. Comp. Physics 179
(2), 439 (2002).
[22] S. Ansumali and I. V. Karlin, Phys. Rev. E 65, 056312
(2002).
[23] Y. Shi, T. S. Zhao and Z. L. Guo, Phys. Rev.E 73, 026704
2006.
[24] X. Shan and X. He, Phys. Rev. Lett. 80, 65 (1998).
[25] S. Ansumali, I.V. Karlin, and H. C. Öttinger, Euro-
phys.Lett. 63, 798 (2003).
[26] S.S. Chikatamarla and I.V. Karlin, Phys. Rev.Lett. 97,
190601 (2006).
[27] A. Bardow, I.V. Karlin, and A. A. Gusev, Europhys.
Lett. 75, 434 (2006).
[28] S. Ansumali and I.V. Karlin, Phys. Rev. Lett. 95, 260605
(2005).
[29] S. Succi, I.V. Karlin, and H. Chen, Rev. Mod. Phys. 74,
1203 (2002).
[30] B. M. Boghosian, J. Yepez, P.V. Coveney, and A. J. Wag-
ner, Proc. R. Soc. A 457, 717 (2001).
[31] M. E. McCracken and J. Abraham, Phys. Rev.E 71,
036701 (2005).
[32] I.V. Karlin and S. Succi, Phys. Rev. E 58, R4053 (1998).
[33] I.V. Karlin, A. N. Gorban, S. Succi, and V. Boffi,
Phys.Rev. Lett. 81, 6 (1998).
[34] I.V. Karlin, A. Ferrante, and H. C. Öttinger, Europhys.
Lett. 47, 182 (1999).
[35] W.-A. Yong and L.-S. Luo, Phys. Rev. E 67, 051105
(2003).
[36] P. J. Dellar, Europhys. Lett. 57, 690 (2002).
[37] P. Lallemand and L.-S. Luo, Phys. Rev. E 61, 6546
(2000).
[38] P. Lallemand and L.-S. Luo, Phys. Rev. E 68, 036706
(2003).
[39] E.T. Shannon, Phys. Rev. 106, 620 (1957).
[40] T. Abe, J. Comput. Phys. 131, 241 (1997).
[41] N. Cao, S.Chen, S.Jin, and D. Mart́ınez, Phys. Rev. E
55, R21 (1997).
[42] Y.-H. Qian, D. d’Humieres, and P. Lallemand, Euro-
phys.Lett. 17, 479 (1992).
[43] P. L. Bhatnagar, E. P. Gross, and M. Krook, Phys. Rev.
A 94,511 (1954).
[44] X. He and L.S. Luo, Phys. Rev. E 55, R6333 (1997).
[45] Q. Zou and X. He, Phys. Fluids 9, 1951 (1997).
[46] P.A.Skordos, Phys.Review E 48, 4823 (1993).
[47] R. Comubert, D. d’Humieres, and D. Levermore, Physica
D 47, 241 (1991).
[48] D.P. Ziegler, J. Stat. Phys. 71, 1171 (1993).
[49] I. Ginzbourg and P. M. Adler, J. Phys. II France 4, 191
(1994).
[50] A.J.C. Ladd, J. Fluid Mech. 271, 285 (1994).
[51] D.R. Noble, S. Chen, J.G. Georgiadis, R.O. Buckius,
Phys. Fluids 7 (l), 203 (1995).
[52] Q. Zou and H. He, Phys. Fluids 7, 1788 (1996).
[53] R.S. Maier, R.S. Bernard, and D.W. Grunau, Phys. Flu-
ids 7, 1788 (1996).
[54] S. Chen, Daniel Martınez and R.Mei, Phys. Fluids 8(9),
2528 (1996).
[55] R. Mei, L.S. Luo, and W. Shyy, J. Comput. Phys. 155,
307 (1999).
[56] M. Bouzidi, M. Firadaouss, and P. Lallemand, Phys. Flu-
ids 13, 452 (2001).
[57] S. Ansumali and I. V. Karlin, Phys. Rev. E 66, 026311
(2002).
[58] I. Ginzburg and D.d’Humieres, Phys. Rev. E 68, 066614
(2003).
[59] M. Junk and Z. Yang, Phys. Rev. E 72, 066701 (2005).
[60] S Ansumali and I. V. Karlin, Phys. Rev. E 62, 7999
(2000).
[61] S. Ansumali, I.V. Karlin, and H. C. Ottinger, Europhys.
Lett. 63, 798 (2003).
[62] B.M. Boghosian, P.J. Love, P. V. Coveney, I.V. Karlin,
S.Succi, J.Yepez, Phys.Rev E 68, 025103(R) (2003).
[63] A.J. Chorin, J.Comp.Phys. 2, 12 (1967).
