Stability and stabilisation of the lattice Boltzmann method: Magic steps and salvation operations

1.2 MB·English

by  R. A. Brownlee

Stability and stabilisation of the
lattice Boltzmann method
Magic steps and salvation operations
R. A. Brownlee,∗ A. N. Gorban, and J. Levesley
7 Department of Mathematics, University of Leicester, Leicester LE1 7RH, UK
0
0 We revisit the classical stability versus accuracy dilemma for the lattice Boltzmann methods
2 (LBM). Our goal is a stable method of second-order accuracy for fluid dynamics based on the
n lattice Bhatnager–Gross–Krook method (LBGK).
a TheLBGKschemecanberecognisedasadiscretedynamicalsystemgeneratedbyfree-flightand
J entropic involution. In this framework the stability and accuracy analysis are more natural. We
9 find the necessary and sufficient conditions for second-order accurate fluid dynamics modelling. In
1 particular, it is proven that in order to guarantee second-order accuracy the distribution should
belong to a distinguished surface – the invariant film (up to second-order in the time step). This

] surface is thetrajectory of the(quasi)equilibrium distribution surface underfree-flight.
h The main instability mechanisms are identified. The simplest recipes for stabilisation add no
c
artificial dissipation (up to second-order) and provide second-order accuracy of the method. Two
e
other prescriptions add some artificial dissipation locally and prevent the system from loss of posi-
m
tivityandlocalblow-up. DemonstrationoftheproposedstableLBGKschemesareprovidedbythe
- numerical simulation of a 1D shock tubeand the unsteady 2D-flow around a square-cylinderup to
t
a ReynoldsnumberO(10000).
t
s
.
t
a I. INTRODUCTION making the relaxation time sufficiently small. However,
m in this low-viscosity regime, LBGK suffers from numeri-
- A lattice Boltzmann method (LBM) is a discrete ve- calinstabilitieswhichreadilymanifestthemselvesaslocal
d locity method in which a fluid is described by associat- blow-ups and spurious oscillations.
on ing, with each velocity vi, a single-particle distribution Anotherproblemisthedegreeofaccuracy. Anapprox-
c function fi = fi(x,t) which is evolved by advection and imation to the continuous-in-time kinetics is not equiva-
[ interaction on a fixed computational lattice. lent to an approximation of the macroscopic transport

The method has been proposed as a discretization of equation. The fluid dynamics appears as a singular
3
Boltzmann’s kinetic equation (for an introduction and limit of the Boltzmann or BGK equation for small τ.
v
4 historicreviewsee[53]). Furthermore,thecollisionoper- Anapproximationtothe correspondingslowmanifoldin
4 ator can be alluringly simplified, as is the case with the the distribution space is constructed by the Chapman–
4 Bhatnager–Gross–Krook (BGK) operator [6], whereby Enskogexpansion. Thisisanasymptoticexpansion,and
1 collisions are described by a single-time relaxationto lo- higher(Burnett)termscouldhavesingularities. Analter-
1 cal equilibria f∗: nativeapproachtoasymptoticexpansion(with“diffusive
6 i
scaling”insteadof“convectivescaling”intheChapman–
0 ∂fi +v f = 1(f∗ f ). (1) Enskog expansion) was developed in [37] in order to ob-
mat/ Thephysically∂retasonaib·l∇echioiceτforif−i∗isiasentropymax- tfraoinmtkhienienticcosm. pressible Navier–Stokesequationsdirectly
imizers, although other choices of equilibria are often
d- preferred [53]. The local equilibria fi∗ depend nonlin- atiIotnapscpheeamrsetthoatthteheslroewlaxhaytdiorondtyinmaemoicf mthaenoifvoeldrremlaaxy-
n early on the hydrodynamic moments (density, momen-
be quite large for small viscosity: t c2δ2/(2ν)
o tum, etc.). These moments are linear functions of fi, δ2/(2τ) (see below, in Sec. III). Somreelaexst∼ima1tets of lon∼g
c hence (1) is a nonlinear equation. For small τ, the t
: relaxationtime for LBGK at large Reynolds number are
v Chapman–Enskog approximation [14] reduces (1) to the
found earlier in [58]. So, instead of fast relaxation to
i compressibleNavier–Stokesequation[53]withkinematic
X a slow manifold in continuous-in-time kinetics, we could
viscosity ν τc2, where c is the thermal velocity for
r ∼ 1 1 meet a slow relaxation to a fluid dynamics manifold in
one degree of freedom.
a
the chain of discrete LBM steps.
The overrelaxation discretization of (1) (see, e.g., [5,
16,33, 39,40, 53])isknownasLBGK,andallowsone to Our approach is based on two ideas: the Ehrenfests’
chooseatime step δ τ. This decouplesviscosityfrom coarse-graining [22, 24, 30] and the method of differen-
t
≫
the time step, thereby suggesting that LBGK is capa- tial approximation of difference equations [34, 50]. The
ble of operating at arbitrarily high-Reynolds number by background knowledge necessary to discuss the LBM in
this manner is presented in Sect. II. In this section, we
answer the question: how to provide second-order accu-
racy of the LBM methods for fluid dynamics modelling?
∗correspondingauthor: [email protected] We prove the necessary and sufficient conditions for this
2
accuracy. It requires a special connection between the state,andM forthe macroscopicstate. The vectorM is
distribution f and the hydrodynamic variables. There a linear function of f: M =m(f).
i
is only one degree of freedom for the choice of fi, if the b. Equilibrium. For any allowable value of M an
hydrodynamic fields are given. Moreover,the LBM with “equilibrium”distributionshouldbegiven: amicroscopic
overrelaxation can provide approximation of the macro- state f∗ . It should satisfy the obvious, but important
M
scopic equation even when it does not approximate the identity of self-consistency:
continuous-in-time microscopic kinetics.
This approach suggests several sources of numerical m(f∗ )=M, (2)
M
instabilities in the LBM and allows several recipes for
stabilisation. A geometric background for this analysis
or in differential form
providesamanifoldthatisatrajectoryqofthequasiequi-
libriummanifoldduetofree-flight. Wecallthismanifold ∗ ∗
m(D f ) 1, i.e., m((D f )a) a. (3)
the invariant film (of nonequilibrium states). It was in- M M ≡ M M ≡
troducedin[25]andstudiedfurtherin[24,27,28]. Com-
Thestatef∗ isnotaproperthermodynamicequilibrium,
monto eachstabilisationrecipe is the desire to stayuni- M
but a conditional one under the constraint m(f) = M.
formly close to the aforementioned manifold (Sect. III).
Thereforewecallitaquasiequilibrium (othernames,such
In Sect. IV, in addition to two LBM accuracy tests,
as localequilibrium, conditionalequilibrium, generalised
a numerical simulation of a 1D shock tube and the
canonical state or pseudoequilibrium are also in use).
unsteady 2D-flow around a square-cylinder using the
For the quasiequilibrium f∗ , an equilibration opera-
presentstabilisedLBMarepresented. Forthelaterprob- M
∗
tion is the projection Π of the distribution f into the
lem, the simulation quantitatively validates the experi-
corresponding quasiequilibrium state: Π∗(f)=f∗ .
mentally obtained Strouhal–Reynolds relationship up to m(f)
Re = (10000). This extends previous LBM studies of In the fully physical situation with continuous veloc-
O ∗
this problem where the relationship had only been suc- ity space, the quasiequilibrium f is defined as a condi-
M
cessfully validated up to Re= (1000) [1, 4]. tional entropy maximizer by a solution of the optimisa-
O
Sect. V contains some concluding remarks as well as tion problem:
practical recommendations for LBM realisations.
Weuseoperatornotationthatallowsustopresentgen- S(f) max, m(f)=M, (4)
→
eralresultsincompactform. Theonlydefinitionwehave
torecallhereisthe(Gaˆteaux)differential: thedifferential where S(f) is an entropy functional.
of a map J(f) at a point f0 is a linear operator (DfJ)f0 The choice of entropy is ambiguous; generally, we can
defined by a rule: (DfJ)f0g = ddε(J(f0+εg))ε=0. start from a concave functional of the form
S(f)= s(f(x,v,t))f(x,v,t)dxdv (5)
II. BACKGROUND
Z
a. Microscopic and macroscopic variables. Let us with a concave function of one variable s(f). The choice
describe the main elements of the LBM construction. by default is s(f) = lnf, which gives the classical
−
The first element is a microscopic description, a single- Boltzmann–Gibbs–Shannon (BGS) entropy.
particle distribution function f(x,v), where x is the For discrete velocityspace, there existsome extra mo-
spacevector,andvisvelocity. Ifvelocityspaceisapprox- ment conditions on the equilibrium construction: in ad-
imated by a finite set v , then the distribution is ap- dition to (2) some higher moments of a discrete equilib-
i
proximated by a measu{re}with finite support, f(x,v) riumshould be the same as for the continuousone. This
f δ(v v ). In thatcase,the microscopicdescriptio≈n isnecessarytoprovidethepropermacroscopicequations
is tihei fini−te-diimensional vector-function f (x). for M. Existence of entropy for the entropic equilibrium
i
PThe second main element is the macroscopic descrip- definition (4) whilst fulfilling higher moment conditions
tion. This is a set of macroscopic vector fields that couldbeincontradiction,andaspecialchoiceofvelocity
are usually some moments of the distribution func- set may be necessary (for a very recent example of such
tion. The main example gives the hydrodynamic fields researchfor multispeed lattices see [17]). Another choice
(density–momentum–energy density): n,nu, (x) = is to refuse to deal with the entropic definition of equi-
1,v,v2/2 f(x,v)dv. But this is no{t an oEbl}igatory librium (4) and assume that there will be no perpetuum
{ }
choice. If we would like to solve by LBM methods the mobile of the second kind. This extends the possibility
GR rad equations [31, 49] or some extended thermody- for approximation, but creates some risk of nonphysical
namic equations [36], we should extend the list of mo- behavior of the model. For a detailed discussion of the
ments (but, at the same time, we should be ready to in- H-theorem for LBM we refer the readers to [54].
troduce more discrete velocities for a proper description Some of the following results depend on the entropic
of these extended moment systems). definition of equilibrium, but some do not. We always
In general, we use the notation f for the microscopic point out if results are “entropy–free”.
3
c. Free flight. In the LBM construction the other The first term in the right hand side of (9) – the
mainelementsare: thefree-flighttransformationandthe quasiequilibrium approximation – consists ofmoments of
collision. There are many models of collisions, but the df/dt computed at the quasiequilibriumpoint. For free-
free-flight equation is always the same flight, hydrodynamic fields and Maxwell equilibria this
term gives the Euler equations. The second term in-
∂f
∂t +v·∇xf =0, (6) cdleufdecets otfheinavactriioanncoef∆thf∗e d(iffforerfernetei-aflligDhftJ(c6()f,)tfhM∗isodnifftehre-
with exact solution f(x,v,t) = f(x vt,v,0), or for ential is just v x,Mfor the discrete version (7) this
− is the vector-c−olum·∇n v ). These terms always ap-
discrete velocities, i x
− ·∇
pear in the Chapman–Enskogexpansion. For free-flight,
∂fi +v f =0, (7) hydrodynamic fields and Maxwell equilibria they give
∂t i·∇x i the Navier–Stokes equations for a monatomic gas with
Prandtl number Pr=1:
f (x,t)=f (x v t,0). Free-flightconservesanyentropy
i i i
−
of the form (5). In general, we can start from any dy- ∂n ∂(nu )
i
= ,
namics. For application of the entropic formalism, this
∂t − ∂x
i
dynamicsshouldconserveentropy. Letthiskineticequa- Xi
tion be ∂(nuk) ∂(nukui) 1 ∂P
=
∂t − ∂x − m∂x
df i i k
=J (f). (8) X
dt c τ 1 ∂ ∂uk ∂ui 2
+ P + δ divu ,
ki
2m ∂x ∂x ∂x − 3
Forourconsiderations,thefree-flightequationwillbethe i i (cid:20) (cid:18) i k (cid:19)(cid:21)
X
main example of the conservative kinetics (8).
∂ ∂( u ) 1 ∂(Pu )
i i
The phase flow Θ for kinetic equation(8) is a shift in E = E
t ∂t − ∂x − m ∂x
time that transforms f(t0) into f(t0+t). For free-flight, Xi i Xi i
Θt :f(x,v)→f(x−vt,v). + τ 5kB ∂ P ∂T ,
Remark. We work with dynamical systems defined 22m2 ∂x ∂x
by partial differential equations. Strictly speaking, this Xi i (cid:18) i(cid:19)
(11)
means that the proper boundary conditions are fixed.
In order to separate the discussion of equation from a
where m is particle mass, k is Boltzmann’s constant,
boundary condition problem, let us imagine a system B
P =nk T isidealgaspressure,T iskinetictemperature,
withperiodicboundaryconditions(e.g.,onatorus),ora B
and the underlined terms are the results of the coarse-
system with equilibrium boundary conditions at infinity.
d. Ehrenfests’ solver of second-order accuracy for the graining additional to the quasiequilibrium approxima-
Navier–Stokes equations. Here we present a generali- tion.
sation of a well known result. Let us study the fol- All computations are straightforwardexercises (differ-
lowing process (an example of the Ehrenfests’ chain entialcalculusandGaussianintegralsforcomputationof
[22, 24, 30], a similar result gives the optimal predic- the moments, m, in the continuous case). More details
tion approach [18]): free-flight for time τ – equilibration of these computations are presented in [28].
p–rforceees-fls,igthhtefhorydtrimodeyτna–meiqcufiileibldrsataiopnpr–ox·i·m·.atDeutrhinegsothluis- fulTthoecdoymnapmariectvhisicsofsoitrymiunla(1t1o)tisheµm=eaτ2nn-kfBreTe-.pIattihsuthsee--
tion of the (compressible) Navier–Stokes equation with ory that gives µ = τcolnkBT, where τcol is the collision
viscosity ν τc2, where c is the thermal velocity for time (the time for the mean-free-path). According to
∼ 2 1 1
one degree of freedom. The error of one step of this ap- theseformulas,we getthe followinginterpretationofthe
proximation has the order (τ3). An exact expression coarse-grainingtime τ for this example: τ =2τcol.
for the transport equation tOhat is approximated by this For any particular choice of discrete velocity set v
i
process in the general situation (for arbitrary initial ki- andofequilibriumf∗ thecalculationcouldgivediffe{ren}t
M
netics, velocity set and for any set of moments) is: equations,butthe generalformula(9)remainsthe same.
The connection between discretization and viscosity was
dM τ
∗ alsostudiedin[51]. Letusprovethegeneralformula(9).
dt =m(Jc(fM))+ 2m((DfJc(f))fM∗ ∆fM∗ ), (9) We are looking for a macroscopic system that is ap-
where ∆f∗ is the defect of invariance of the quasiequi- proximated by the Ehrenfests’ chain. Let us look for
M macroscopic equations of the form
librium manifold:
∆fM∗ =Jc(fM∗ )−DM(fM∗ )m(Jc(fM∗ )), (10) dM =Ψ(M) (12)
dt
and is the difference between the vector-field J and its
c
projection on to the quasiequilibrium manifold. This re- with the phase flow Φ : M(t) = Φ M(0). The trans-
t t
sult is entropy-free. formation Φ should coincide with the transformation
τ
4
M m(Θ (f∗ )) up to second-order in τ. The match- because of the self-consistency identity (2), (3).
7→ τ M
ing condition is Due to the presence of the small parameter τ in J(f),
the zeroth approximation to f is the quasiequilibrium
∗ M
m(Θτ(fM))=Φτ(M) for every M and given τ. (13) approximation: fM(0) = fM∗ . Let us look for fM in the
form of a power series: f = f(0) +τf(1) + , with
Thisconditionistheequationforthemacroscopicvector M M M ···
(k)
fieldΨ(M). Thesolutionofthisequationisafunctionof m(f )=0 for k 1. From (18) we immediately find:
M ≥
τ: Ψ = Ψ(M,τ). For a sufficiently smooth microscopic
vector field Jc(f) and entropy S(f) it is easy to find the fM(1) =Jc(fM(0))−(DMfM(0))(m(Jc(fM(0))))=∆fM∗ . (19)
Taylor expansion of Ψ(M,τ) in powers of τ. Let us find
the firsttwoterms: Ψ(M,τ)=Ψ (M)+τΨ (M)+o(τ). It is very natural that the first term of the Chapman–
0 1
Up to second-order in τ the matching condition (13) is Enskogexpansionfor the modelequation(16) isjust the
defect of invariance for the quasiequilibrium (10).
τ2 The corresponding first-order in τ approximation for
∗ ∗
m(Jc(fM))τ +m((DfJc(f))fM∗ (Jc(fM))) 2 the macroscopic equations is:
τ2
=Ψ (M)τ +Ψ (M)τ2+(D Ψ (M))(Ψ (M)) . dM ∗
0 1 M 0 0 2 dt =m(Jc(fM))+τm((DfJc(f))fM∗ ∆fM∗ ). (20)
(14)
Weshouldrecallthatm(∆f∗ )=0. Thelasttermin(18)
From this condition immediately follows: M
vanishesinthemacroscopicprojectionforallorders. The
∗ onlydifferencebetween(20)and(9)isthecoefficient1/2
Ψ (M)=m(J (f )),
0 c M before τ in (9).
1 (15)
f. Decoupling of time step and viscosity: how to pro-
Ψ1(M)= m((DfJc(f))f∗ ∆f∗ ),
2 M M vide second-order accuracy? In the Ehrenfests’ chain
“free-flight – equilibration – ” the starting point of
where∆fM∗ isthe defect ofinvariance(10). Thus wefind each link is a quasiequilibriu·m··state: the chain starts
thatthemacroscopicequationinthefirstapproximation
fromf∗ ,then,afterfree-flight,equilibratesintof∗ ,
is (9). M(0) M(τ)
e. The Chapman–Enskog expansion for the gener- etc. The viscosity coefficient in (9) is proportional to τ.
alised BGK equation. Here we present the Chapman– Let us choose another starting point fMs in order to de-
couple time step and viscosity and preserve the second-
Enskog method for a class of generalised model equa-
order accuracy of approximation. We would like to get
tions. This class includes the well-known BGK kinetic
equation (9) with a chain time step δ =h. Analogously
equation, as well as many other model equations [26]. t
to (14) and (15), we obtain the macroscopic equation
As a starting point we take a formal kinetic equation
with a small parameter τ
dM
df 1 dt =m(Jc(fM∗ ))+m((DfJc(f))fM∗ ((fMs −fM∗ )+(h/2)∆fM∗ )),
∗
=J(f):=Jc(f)+ (Π (f) f). (16) (21)
dt τ − under the condition that fs f∗ = (h). The initial
ThetermΠ∗(f) f isnonlinearbecauseofthenonlinear point M − M O
dependency of Π−∗(f)=f∗ on m(f).
m(f) 1
We would like to find a reduced description valid for fs =f∗ (h τ)∆ +o(h) (22)
M M − 2 − fM∗
the macroscopic variables M. This means, at least, that
we are looking for an invariant manifold parameterised
provides the required viscosity. This is a sufficient con-
by M, f =fM, that satisfies the invariance equation: dition for the second-order accuracy of the approxima-
tion. Ofcourse,theself-consistencyidentitym(fs )=M
(D f )(m(J(f )))=J(f ). (17) M
M M M M should be valid exactly, as (2) is. This starting distribu-
tionisalinearcombinationofthe quasiequilibriumstate
The invariance equation means that the time deriva-
and the first Chapman–Enskogapproximation.
tive of f calculated through the time derivative of M
(M˙ = m(J(f ))) by the chain rule coincides with the Thenecessaryandsufficientconditionforsecond-order
M accuracy of the approximation is:
true time derivative J(f). This is the central equation
for model reduction theory and applications. The first
1
generalresultsaboutexistenceandregularityofsolutions m (DfJc(f))fM∗ (fMs −fM∗ + 2(h−τ)∆fM∗ ) =o(h)
to (17)wereobtainedbyLyapunov[45](see,e.g.,the re- (cid:18) (cid:19)
(23)
viewin[28]). Forthekineticequation(16)theinvariance
(with the self-consistency identity m(fs ) = M). This
equation has the form M
means, that the difference between left and right hand
1 sides of (22) should have zero moments and give zero
∗
(DMfM)(m(Jc(fM)))=Jc(fM)+ τ(fM −fM), (18) inputs in observable macroscopic fluxes.
5
Hence, the condition of second-order accuracy signif- and q f∗ . The quasiequilibrium manifold divides
M,0 ≡ M
icantly restricts the possible initial point for free-flight. qinto twoparts,q=q− q0 q+, whereq− = qM,t t<
∪ ∪ { |
This result is also entropy-free. 0 , q = q t > 0 , and q is the quasiequilibrium
+ M,t 0
Anyconstructionofcollisionsshouldkeepthesystem’s m}anifold: q{ = |q }= f∗ .
0 { M,0} { M}
startingfree-flightstepsnearthepointsfs givenby(22) There is an important temporal involution of the film:
M
and (23). The conditions (22) and (23) for second-order
accuracyofthetransportequationapproximationdonot IT(qM,t)=qM,−t. (26)
depend on a specific collision model, they are valid for
mostmodificationsoftheLBMandLBGKthatusefree- Dueto(22),forqM,t andagiventimestephthetrans-
flight as a main step. formation M m(Θh(qM,t)) approximates the solution
7→
Various multistep approximations give more freedom of (9) with τ = 2t+h for the initial conditions M and
of choice for the initial state. For the construction of time step h withsecond-orderaccuracyinh. Hence, due
such approximations below, the following mean viscosity to the mean viscosity lemma, the two-step transforma-
lemma is important: if the transformations Ωh : M tion
M′, i = 1,...,k, approximate the phase flow fior (9) f→or
M m(I (Θ (I (Θ (q ))))) (27)
time h (shift intime h) andτ =τ with second-orderac- T h T h M,t
i 7→
curacy in h, then the superposition ΩhΩh Ωh approx-
1 2··· k approximates the solution of (9) with τ = 0 (the Euler
imates the phase flow for (9) for time kh (shift in time
equations) for the initial conditions M and time step h
kh)fortheaverageviscosityτ = 1(τ + +τ )withthe
k 1 ··· k with second-order accuracy in h. This is true for any t,
same order of accuracy. The proof is by straightforward
hence, for any starting point on the invariant film with
multiple applications of Taylor’s formula.
g. EntropicformulaforD (f∗ ). Amongthemany the given value of M.
M M Toapproximatethe solutionof (9)with nonzeroτ, we
benefits of thermodynamics for stability analysis there
need an incomplete involution:
are some technical issues too. The differential of equilib-
riumD (f∗ ) appearsinmany expressions,for example
(3), (10M), (1M4), (15), (17) and (18). If the quasiequilib- ITβ(qM,t)=qM,−(2β−1)t. (28)
rium is defined by the solution of the optimisation prob-
For β = 1, we have I1 = I and for β = 1/2, I1/2 is
lem (4), then T T T
just the projection onto the quasiequilibrium manifold:
DMfM∗ = Df2S −fM∗1 mT m Df2S −fM∗1 mT −1. (24) tIhT1/e2f(oqlMlo,wt)in=gsΠeq∗(uqeMnc,te)g=iveqsMa,0s.ecAofntde-rosrodmereinintiitmiaelsstteeppsh,
This operator(cid:0) is c(cid:1)onstruct(cid:16)ed (cid:0)from(cid:1)the vec(cid:17)tor m, the approximationof (9) with τ =(1 β)h/β, 1/2 β 1:
− ≤ ≤
transposed vector mT and the second differential of en-
tropy. The inverse Hessian (∂2S/∂f ∂f )−1 is especially M =m((IβΘ )nq ). (29)
i j n T h M,t
simple for the BGS entropy, it is just f δ . The for-
i ij
mula (24) was first obtained in [48] (for an important To prove this statement we consider a transformation of
particular case; for further references see [28]). the second coordinate in q by IβΘ :
M,ϑn T h
h. Invariant film. All the points Θ (f∗ ) belong to
t M
a manifold that is a trajectory q of the quasiequilibrium ϑn+1 = (2β 1)(ϑn+h). (30)
− −
manifoldduetotheconservativedynamics(8)(inhydro-
This transformation has a fixed point ϑ∗ = h(2β
dynamic applications this is the free-flightdynamics (6).
We call this manifold the invariant film (of nonequilib- 1)/(2β) and ϑn = ϑ∗ +( 1)n(2β 1)nδ for so−me δ. −If
− −
rium states). It was introduced in [25] and studied fur- 1 β is small then relaxation may be very slow: ϑn
ϑ∗−+( 1)nδexp( 2n(1 β)), and relaxationrequires ≈
tthanergeinnt[2to4,q2a7t,t2h8e].pToihnetdfe∗fe,catnodfbinevloanrigasntcoet∆hefM∗int(e1r0s)ecis- 1/(2(1− β)) step−s. If ϑn−= ϑ∗+o(h) then the sequenc∼e
M M (29−)approximates(9)withτ =h 2ϑ∗ =(1 β)h/β
tionofthistangentspacewithkerm. Thisintersectionis n
− | | −
and second-orderaccuracy in the time step h. The fixed
one-dimensional. This means that the direction of ∆
fM∗ pointsq coincidewiththerestartpointsf∗ +ϑ∗∆
is selected from the tangent space to q by the condition: M,ϑ∗ ∗ M fM∗
(22)inthefirstorderinϑ = (h τ)/2,andthemiddle
derivative of M in this direction is zero.
pointsϑ∗+h/2ofthefree-fligh−tju−mpsq q
A point f on the invariant film q is naturally param- M,ϑ∗ M′,ϑ∗+h
7→
approximate the first-order Chapman–Enskog manifold
eterised by (M,t): f = q , where M = m(f) is the
value of the macroscopic vMa,rtiables, and t = t(f) is the fM∗ + τ2∆fM∗ .
For the entropic description of quasiequilibrium, we
time shift from a quasiequilibrium state: Θ−t(f) is a
can connect time with entropy and introduce entropic
quasiequilibrium state for some (other) value of M. By
coordinates. For each M and positive s from some in-
definition, the action of Θ on the second coordinate of
t
q is simple: Θ (q )=q . To the first-order in terval 0 < s < ς there exist two numbers t±(M,s)
M,t t M,τ M′,t+τ
t, (t+(M,s)>0, t−(M,s)<0) such that
∗ ∗
qM,t =fM +t∆fM∗ , (25) S(qM,t±(M,s))=S(fM)−s. (31)
6
The numbers t± coincide to the first-order: t+ = t−+ i. Does LBGK with overrelaxation collisions approx-
−
o(t−). imate the BGK equation? The BGK equation as well
We define the entropic involution I as a transforma- as its discrete velocity version (1) has a direction of fast
S
tion of q: contraction Π∗(f) f. The discrete chain (34) with β
−
closeto1hasnothingsimilar. Hence,the approximation
I (q )=q . (32) of a genuine BGK solution by an LBGK chain may be
S M,t± M,t∓
possible only if both the BGK and the LBGK chain tra-
The introduction of incomplete entropic involution Iβ is jectories belong to a slow manifold with high accuracy.
S This implies significant restrictions on initial data and
also obvious (see [24]).
onthe dynamicsofthe approximatedsolution,aswellas
EntropicinvolutionI coincideswith the temporalin-
S fast relaxation of the LBGK chain to the slow manifold.
volution I , up to second-order in the deviation from
T ∗ The usual Taylorseries based arguments from [32] are
quasiequilibrium state f Π (f). Hence, in the vicin-
ity of quasiequilibrium th−ere is no significant difference valid for h ≪ τ. If we assume h ≫ τ, (δt ≫ λ in the
notation of [32]) then Eqn. (10) of [32] transforms (in
between these operations, and all statements about the ournotation)intof(x+vh,v,t+h)=f∗ (x+vh,v,t+
temporal involution are valid for the entropic involution M
h)+ (τ) with M = m(f(x+vh,v,t+h)). That is,
with the same level of accuracy. f(x,vO,t+h)=f∗ (x,v,t+h)+ (τ). Accordingtothis
For the transfer from free-flight with temporal or en- M O
formula, f should almost be at quasiequilibrium after a
tropicinvolutiontothe standardLBGKmodelswemust
time step h τ, with some correction terms of order
transfer fromdynamics and involutiononq to the whole ≫
τ. This first-order in τ correction is, of course, the first
space of states. Instead of ITβ or ISβ the transformation term of the Chapman–Enskog expansion (19): τf(1) =
M
I0β :f 7→Π∗(f)+(2β−1)(Π∗(f)−f) (33) τn∆atfuM∗ra(lwreitshulptofossribalneaeprrporrooxfimoradteiornOo(fτthh)e).BTGhKissioslautvieorny,
especiallyinlightoftheChapman–Enskogexpansion[14,
is used. For β = 1, I01 is a mirror reflection in the 32], but it is not the LBM scheme with overrelaxation.
quasiequilibrium state Π∗(f), and for β = 1/2, I1/2 is The standardelement in the proof of second-orderac-
0
theprojectionontothequasiequilibriummanifold. If,for curacyoftheBGKequationapproximationbyanLBGK
agivenf =q ,the sequence(29)givesasecond-order chain uses the estimation of an integral: for time step h
0 M,t
in time step h approximation of (9), then the sequence we obtain from (1) the exact identity
1 t+h
M =m((IβΘ )nf ) (34) f (x+v h,t+h)= (f∗ (x) f (x,t′))dt′,
n 0 h 0 i i τ i,m(f) − i
Zt
(36)
also gives a second-order approximation to the same where f∗ (x) is the quasiequilibrium state that cor-
equationwithτ =(1 β)h/β. Thischainisthestandard i,m(f)
− responds to the hydrodynamic fields m(f(x,t′)). Then
LBGK model.
one could apply the trapezoid rule for integration to the
Entropic LBGK (ELBGK) methods [8, 24, 42, 54] dif-
right-hand side of (36). The error of the trapezoid rule
fer only in the definition of (33): for β = 1 it should
has the order (h3):
conserve entropy, and in general has the following form: O
t+h h h3
IEβ(f)=(1−β)f +βf˜, (35) Zt Q(t′)dt′ = 2(Q(t)+Q(t+h))− 12Q¨(t′),
where t′ [t,t+h] is a priori unknown point. But for
with f˜= (1 α)f +αΠ∗(f). The number α = α(f) is ∈
− the singularly perturbed system (1), the second deriva-
chosen so that a constant entropy condition is satisfied: tive of the term f∗ (x) f (x,t′) on the right hand
S(f)=S(f˜). For LBGK (33), α=2. i,m(f) − i
side of (36) could be of order 1/τ2, and the whole er-
Of course, computationof Iβ is much easier than that
0 ror estimate is (h3/τ3). This is not small for h > τ.
of ITβ, ISβ or IEβ: it is not necessary to follow exactly the For backward oOr forward in time estimates of the inte-
manifold q and to solve the nonlinear constant entropy gral (36), errors have the order (h2/τ2). Hence, for
condition equation. For an appropriate initial condition overrelaxation with h τ this rOeasoning is not appli-
from q (not sufficiently close to q0), two steps of LBGK cable. Many simple exa≫mples of quantitative and quali-
withIβ givethesamesecond-orderaccuracyas(29). But tative errors of this approximation for a singularly per-
0
alongchainofsuchstepscanleadfarfromthequasiequi- turbed system could be obtained by analysis of a simple
librium manifold and even from q. Here, we see stability systemoftwoequations: x˙ = 1(φ(y) x),y˙ =ψ(x,y)for
τ −
problems arising. For β close to 1, the one-steptransfor- various φ and ψ. There are examples of slow relaxation
mationIβΘ inthechain(34)almostconservesdistances (instead of fast), of blow-up instead of relaxation or of
0 h
between microscopic distributions, hence, we cannot ex- spurious oscillations, etc.
pectfastexponentialdecayofanymode,andthissystem Hence, one cannot state that LBGK with overrelax-
is near the boundary of Lyapunov stability. ationcollisionsapproximatessolutionsoftheBGKequa-
7
tion. Nevertheless,itcandoanotherjob: itcanapproxi- Invariant film
Free flight steps
matesolutionsofthemacroscopictransportequation. As
demonstrated within this section, the LBGK chain (34),
after some initial relaxation period, provides a second-
orderapproximationto the transportequation, if itgoes
close to the invariant film up to the order (h2) (this
O QE manifold
initial relaxation period may have the order (h2/τ)).
O
In other words, it gives the required second-order ap- Overrelaxation steps
proximation for the macroscopic transport equation un-
der some stability conditions.
FIG.1: Directionalinstability: afterseveraliterationsthetra-
III. STABILITY AND STABILISATION jectory is not tangent to the invariant film with the required
accuracy.
A. Instabilities
(cid:507)
j. Positivity loss. First of all, if f is far from the
quasiequilibrium, the state Iβ(f) (33) may be nonphysi-
0
cal. The positivity conditions (positivity of probabilities
orpopulations)maybe violated. Formulti-andinfinite-
dimensionalproblems it is necessaryto specify what one
meansbyfar. Intheprevioussection,f isthewholestate
FIG. 2: Neutral stability and one-step oscillations in a se-
whichincludesthestatesofallsitesofthelattice. Allthe
quenceof reflections. Bold dotted line– aperturbedmotion,
involution operators with classical entropies are defined
∆ – direction of neutral stability.
for lattice sites independently. Violation of positivity at
one site makes the whole state nonphysical. Hence, we
should use here the ℓ∞-norm: close states are close uni- B. Dissipative recipes for stabilisation
formly, at all sites.
k. Large deviations. The second problem is nonlin-
n. Positivity rule. There is a simple recipe for pos-
earity: for accuracy estimates we always use the as-
itivity preservation [11, 56]: to substitute nonpositive
sumption that f is sufficiently close to quasiequilibrium.
Iβ(f)(x) by the closest nonnegative state that belongs
Far from the quasiequilibrium manifold these estimates 0
to the straight line
do not work because of nonlinearity (first of all, the
quasiequilibrium distribution, f∗ , depends nonlinearly
on M and hence the projectionMoperator, Π∗, is nonlin- λf(x)+(1 λ)Π∗(f(x)) λ R (37)
− | ∈
ear). Again we need to keep the states not far from the n o
quasiequilibrium manifold. defined by the two points, f(x) and its corresponding
l. Directional instability. The third problemis a di- quasiequilibrium state. This operation is to be applied
rectional instability that can affect accuracy: the vector point-wise,at the points ofthe lattice where the positiv-
f Π∗(f) candeviate far fromthe tangentto q (Fig.1). ity is violated. The coefficient λ depends on x too. Let
−
Hence, we shouldnot only keepf close to the quasiequi- uscallthis recipe the positivity rule (Fig.3); itpreserves
librium, but also guarantee smallness of the angle be- positivityofpopulationsandprobabilities,butcanaffect
tween the direction f Π∗(f) and tangent space to q. the accuracy of approximation. The same rule is neces-
One could rely on t−he stability of this direction, but sary for ELBGK (35) when a positive “mirror state” f˜
we fail to prove this in any generalcase. The directional withthesameentropyasf doesnotexistsonthestraight
instabilitychangesthestructureofdissipationterms: the line (37).
accuracy decreases to the first-order in time step h and The positivity rule saves the existence of positive so-
significantfluctuationsofthePrandtlnumberandviscos- lutions, but affects dissipation because the result of the
ity, etc. may occur. This carries a danger even without adjusted collision is closer to quasiequilibrium. There
blow-ups;onecouldconceivablyberelyingonnonreliable is a family of methods that modify collisions at some
computational results. points by additional shift in the direction of quasiequi-
m. Direction of neutral stability. Further, there ex- librium. The positivity rule represents the minimal nec-
ists a neutral stability of all described approximations essary modification. It is reasonable to always use this
thatcausesone-steposcillations: asmallshiftoff inthe rule for LBM (as a “salvation rule”).
direction of ∆f∗ does not relax back for β = 1, and its o. Ehrenfests’ regularisation. To discuss methods
M
relaxation is slow for β 1 (for small viscosity). This withadditionaldissipation,theentropicapproachisvery
∼
effect is demonstratedfor a chainof mirror reflections in convenient. Let entropy S(f) be defined for each popu-
Fig. 2. lation vector f = (f ) (below, we use the same letter S
i
- 8
Wehaveintroducedtwoprocedures: thepositivityrule
and Ehrenfests’ regularisation. Both improve stability,
f
1
reduce nonequilibrium entropy, and, hence, nonequilib-
Positivity fixation
rium fluxes. The proper context for discussion of such
f* procedures are the flux-limiters in finite difference and
Positivity domain
finite volume methods. Here we refer to the classi-
cal flux-corrected transport (FCT) algorithm [10] that
f*+(2(cid:533)-1)(f*-f )
1 strictly maintains positivity, and to its further develop-
ments [9, 43, 55].
FIG. 3: Positivity rule in action. The motions stops at the q. Smooth limiters of nonequilibrium entropy. The
positivity boundary. positivity rule and Ehrenfests’ regularisation provide
rare, intense and localised corrections. Of course, it
is easy and also computationally cheap to organize
forlocal-in-spaceentropy,andhope thatthecontextwill
moregentletransformationswithsmoothshiftsofhigher
make this notation clear). We assume that the global
nonequilibriumstatestoequilibrium. Thefollowingregu-
entropyis asumoflocalentropiesforallsites. The local
larisationtransformationdistributesitsactionsmoothly:
nonequilibrium entropy is
∗ ∗
∗ f f +φ(∆S(f))(f f ). (40)
∆S(f)=S(f ) S(f), (38) 7→ −
−
where f∗ is the corresponding local quasiequilibrium at The choice ofthe function φ is highly ambiguous,for ex-
ample, φ=1/(1+α∆Sk) for some α>0, k >0. There
the same point.
are two significantly different choices: (i) ensemble-
The Ehrenfests’ regularisation [11, 12] provides “en-
independent φ (i.e., the value of φ depends on the lo-
tropy trimming”: we monitor local deviation of f from
the corresponding quasiequilibrium, and when ∆S(f,x) calvalue of∆S only)and (ii) ensemble-dependent φ, for
example
exceeds a pre-specified threshold value δ, perform local
Ehrenfests’ steps to the corresponding equilibrium. So 1+(∆S/(αE(∆S)))k−1/2
thatthe Ehrenfests’steps arenotallowedto degradethe φ=
1+(∆S/(αE(∆S)))k
accuracyofLBGKitispertinenttoselecttheksiteswith
highest∆S >δ. The aposterioriestimatesofaddeddis- where E(∆S) is the average value of ∆S in the compu-
sipationcouldeasilybe performedbyanalysisofentropy tational area, k 1 and α & 1. It is easy to select
production in Ehrenfests’ steps. Numerical experiments anensemble-depe≥ndentφwithcontroloftotaladditional
show (see, e.g., [11, 12] and Sect. IV) that even a small dissipation.
number of such steps drastically improves stability. r. ELBGK collisions as a smooth limiter. On the
To avoid the change of accuracy order “on average”, basis on numerical tests, the authors of [56] claim that
the number of sites with this step should be (Nδx/L) thepositivityruleprovidesthesameresults(inthesense
O
whereN isthetotalnumberofsites,δx isthe stepofthe ofstabilityandabsence/presenceofspuriousoscillations)
space discretization and L is the macroscopic character- as the ELBGK models, but ELBGK provides better ac-
isticlength. Butthis roughestimateofaccuracyinaver- curacy.
age might be destroyed by concentrations of Ehrenfests’ For the formal definition of ELBGK (35) our tests do
steps in the most nonequilibrium areas, for example, in not support claims that ELBGK erases spurious oscil-
boundary layers. In that case, instead of the total num- lations (see Sect. IV below). Similar observations for
ber of sites N in the estimate (Nδx/L) we should take Burgers equation has been reported in [7]. We under-
O
the number of sites in a specific region [59]. The effects stand this situation in the following way. The entropic
of concentration could be easily analysed a posteriori. method consists of at least three components:
p. Entropic steps for nonentropic equilibria. If the
approximate discrete equilibrium f∗ is nonentropic, we 1. entropic quasiequilibrium defined by entropy max-
can use ∆S (f) = S (f) instead of ∆S(f), where S imisation;
K K K
−
is the Kullback entropy. This entropy,
2. entropy balanced collisions (35) that have to pro-
f vide proper entropy balance;
i
S (f)= f ln , (39)
K − i f∗
i (cid:18) i (cid:19) 3. a method for the solution of the transcendental
X equation S(f)=S(f˜) to find α=α(f) in (35).
gives the physically reasonable entropic distance from
equilibrium, if the supposed continuum system has the It appears that the first two items do not affect spuri-
classical BGS entropy. In thermodynamics, the Kull- ous oscillations at all, if we solve equation for α(f) with
back entropy belongs to the family of Massieu–Planck– highaccuracy. Additionalviscositycould,potentially,be
Kramers functions (canonical or grand canonical poten- added by use of explicit analytic formulas for α(f). In
tials). Onecanuse(39)intheconstructionofEhrenfests’ order not to decrease entropy, the errors in these formu-
regularisationfor any choice of discrete equilibrium. las always increase dissipation. This can be interpreted
9
as a hidden transformation of the form (40), where the where
coefficients of φ also depend on f∗.
1
Comparedtofluxlimiters,nonequilibriumentropylim- ∆+ = (Θ (f∗ ) Π∗(Θ (f∗ ))),
fM∗ h h M − h M
iters have a great benefit: by summation of all entropy
1
cphaatinognesthweelicmaniteersstiimntartoedtuhceeainmtooutnhteosfysatdedmit.ional dissi- ∆−fM∗ =−h(Θ−h(fM∗ )−Π∗(Θ−h(fM∗ ))).
Therearenoerrorsofthefirst-orderin(41). Theforward
(∆+ ) and backward (∆− ) approximations are one or-
f∗ f∗
C. Non-dissipative recipes for stabilisation derMlessaccurate. ThecomMputationof∆± isofthesame
f∗
M
computationalcostasanLBGKstep,hence,ifweusethe
s. Microscopic error and macroscopic accuracy. restartformula (22) with central difference evaluation of
The invariantfilmqisaninvariantmanifoldforthe free- ∆f∗ (41), then the computational cost increases three
M
flight transformation and for the temporal and entropic times (approximately). Non-localityofcollisions(restart
involutions. The linear involution I0, as well as the EL- from the distribution fMs (22) with a nonlocal expres-
fB′GwKitihnvfolutfi′o=nIαE∆,tΠr∗a(nfs)f+orom(fsapΠo∗in(ft)f),∈i.qe.i,ntthoeavpecotinotr sfrioeen-flfoigrh∆t faM∗nd) slpocoaillsctohlelismioanisn:LnBoMnloicdaelaityofisexliancetalrinaenadr
f f′ is “−almosttangent”toq,−andthe distance fromf′ exact,nonlinearityislocal[53]). Onemightalsoconsider
to−q has the order ( f Π∗(f) 2). the inclusionofother finite difference representationsfor
O k − k
Hence,iftheinitialstatebelongstoq,andthedistance ∆ into explicit LBM schemes. The consequences of
f∗
M
from quasiequilibrium is small enough ( (h)), then this combination should be investigated.
∼ O
during several steps the LBGK chain will remain near q u. Coupled steps with quasiequilibrium ends. The
with deviation (h2). Moreover, because errors pro- mean viscosity lemma allows us to combine different
∼ O
duced by collisions (deviations from q) have zero macro- starting points in order to obtain the necessary macro-
scopicprojection,thecorrespondingmacroscopicerrorin scopic equations. From this lemma, it follows that
M during several steps will remain of order (h3). the following construction of two coupled steps with
O
To demonstrate this, suppose the error in f, δf, is restart from quasiequilibrium approximates the macro-
of order (hk), and m(δf) = 0, then for smooth fields scopic equation (9) with second-order accuracy in time
O
after a free-flight step an error of higher order appears step h.
in the macroscopic variables M: m(Θ (δf))= (hk+1), Let us take f∗ as the initial state with givenM, then
h O M
because m(Θ (δf)) = m((Θ 1)(δf)) and Θ 1 = evolve the state by Θ , apply the incomplete temporal
h h h h
(h). The last estimate require−s smoothness. − involutionIβ (28),againevolvebyΘ ,andfinallyproject
OThis simple statement is useful for the error analysis by Π∗ ontoTthe quasiequilibrium mahnifold:
we perform. We shall call it the lemma of higher macro-
scopic accuracy: a microscopic error of order O(hk) in- M 7→M′ =m(Π∗(Θh(ITβ(Θh(fM∗ ))))). (42)
duces, after a time step h, a macroscopic error of order
(hk+1), if the field of macroscopic fluxes is sufficiently It follows from the restart formula (22) and the mean
O viscosity lemma that this step gives a second-order in
small (here, the microscopic error means the error that
time happroximationto the shiftintime 2hfor(9)with
has zero macroscopic projection).
t. Restarts and approximation of ∆f∗ . The prob- τ = 2(1−β)h, 1/2 ≤ β ≤ 1. Now, let us replace ITβ
M by the much simpler transformation of LBGK collisions
lem of nondissipative LBM stabilisation we interpret as
a problem of appropriate restart from a point that is I0β (33):
sufficiently close to the invariant film. If h = τ and
collisions return the state to quasiequilibrium, then the M 7→M′ =m(Π∗(Θh(I0β(Θh(fM∗ ))))). (43)
state belongs to q for all time with high accuracy. For
According to the lemma of higher macroscopic accu-
h=τ, formulasfor restartingare also available: one can
6
choose between (22) and, more flexibly, (23). Neverthe-
(cid:300)
less, many questions remain. Firstly, what should one (cid:306)
take for ∆f∗ ? This vector has a straightforward differ-
ential definiMtion (10) (let us also recall that τ∆ is the (cid:518)*
f∗
M
first Chapman–Enskog nonequilibrium correction to the
distributionfunction(19)). Butnumericaldifferentiation
I(cid:69)
couldviolatetheexact-in-spacefree-flighttransformation 0 QE-manifold
and local collisions. There exists a rather accurate cen- I
0
traldifferenceapproximationof∆f∗ onthebasisoffree-
M
flight:
FIG. 4: The scheme of coupled steps (43).
1
∆fM∗ = 2(∆+fM∗ +∆−fM∗ )+O(h2), (41) racy this step (Fig. 4) also gives a second-order in time
10
h approximation to the shift in time 2h for (9) with Method A,
•
τ = 2(1 β)h, 1/2 β 1. The replacement of Iβ
by I0β int−roduces an≤error≤in f that is of order O(h2)T, f˜n+1 =f˜n+Π∗(Θh(f˜n))−Θh(f˜n), (45)
but both transformations conserve the value of macro-
scopic variables exactly. Hence (due to the lemma of Method B,
•
higher macroscopic accuracy) the resulting error of cou-
pledsteps(43)inthemacroscopicvariablesM isoforder f˜n+1 =Θ−h(Π∗(Θh(f˜n)))
(46)
Ocu(rha3c)y..Thismeansthatthemethodhassecond-orderac- +fM∗ −Π∗(Θ−h(Π∗(Θh(f˜n)))).
Letusenumeratethemacroscopicstatesin(43): M =
0 Method A is: shoot from the previous approximation,
M, M1/2 = m(Θh(fM∗ )) and M1 = M′. The shift from f˜ , by Θ , project onto quasiequilibrium, Π∗(Θ (f˜ )),
M toM approximatestheshiftintimehfor(9)with n h h n
0 1/2 andthencorrectionoff˜ bythefinalpointdisplacement,
τ = h. If we would like to model (9) with τ h, then n
≪ Π∗(Θ (f˜ )) Θ (f˜ ). The value of M does not change,
τ = h means relatively very high viscosity. The step becauhse mn (Π−∗(Θh(f˜n )))=m(Θ (f˜ )).
from M to M has to normalize viscosity to the re- h n h n
1/2 1
questedsmallvalue (compareto antidiffusion in [9, 10]). Method B is: shoot from the previous approximation,
f˜ , by Θ , project onto quasiequilibrium, shoot back-
The antidiffusion problem necessarily appears in most n h
CFD approachesto simulation flows with high-Reynolds wardsbyΘ−h,andthencorrectionofM usingquasiequi-
libria (plus the quasiequilibrium with required value of
numbers. Another famous example of such a problem
M, and minus one with current value of M).
is the filtering-defiltering problem in large eddy simula-
tion (LES) [52]. The antidiffusion in the coupled steps The initial approximationcould be f˜1/2, and n here is
is produced by physical fluxes (by free-flight) and pre- the number of iteration. Due to the lemma of higher
serves positivity. The coupled step is a transformation macroscopic accuracy, each iteration (45) or (46) in-
M M and takes time 2h. The middle point M is creasesthe orderofaccuracy(see alsothe numericaltest
0 1 1/2
an a7→uxiliary state only. in Sec. IV).
Let us enumerate the microscopic states in (43): f = The shooting method A (45) better meets the main
0
f∗ , f− = Θ (f∗ ), f0 = Π∗(f− ), f˜ = I1(f− ), LBMidea: eachchangeofmacroscopicvariableisdueto
M 1/2 h M 1/2 1/2 1/2 0 1/2 afree-flightstep(becausefree-flightinLBMisexact),all
f+ = f0 +(1 β)(f˜ f0 ), f− = Θ (f+ ) f =
1/2 1/2 − 1/2 − 1/2 1 h 1/2 1 other operations effect nonequilibrium component of the
Π∗(f−)=f∗ , where M =m(f−). Here, in the middle distribution only. The correction of M in the shooting
1 M1 1 1
of the step, we have 4 points: a free-flight shift of the method B (46) violates this requirement.
−
initial state (f ), the corresponding quasiequilibrium The idea that all macroscopic changes are projections
1/2
(f0 ), the mirror image (f˜ ) of the point f− with of free-flight plays, for the proposed LBM antidiffusion,
1/2 1/2 1/2 thesameroleasthemonotonicityconditionforFCT[10].
respect to the centre f0 , and the state (f+ ) that is
1/2 1/2 In particular, free-flight never violates positivity.
the image of f− after homothety with centre f0 and Ifweafindsolutionf˜totheantidiffusionproblemwith
1/2 1/2
coefficient 2β 1. M = M , then we can take f+ = f0 +(1 β)(f˜
For smooth−fields, the time shift Θh returns f˜1/2 to f10/2), an1d/2M1 = m(Θh(f1+/2)). 1B/u2t eve1n/2exact−solution−s
the quasiequilibrium manifold with possible error of or-
of (44) can cause stability problems: the entropy of f˜
der (h2). For entropic equilibria, the nonequilibrium
−
O couldbelessthantheentropyoff ,andblow-upcould
entropy of the state Θ (f˜ ) is of order (h4). This is 1/2
h 1/2 O appear. A palliative solution is to perform an entropic
anentropicestimateoftheaccuracyofantidiffusion: the
nonequilibrium entropy of f˜1/2 could be estimated from step: tofindαsuchthatS(f10/2+α(f˜−f10/2))=S(f1−/2),
below as C(M)h2, where C(M) > 0 does not depend then use f+ = f0 + (1 β)α(f˜ f0 ). Even for
1/2 1/2 − − 1/2
on h. The problem of antidiffusion can be stated as an nonentropic equilibria it is possible to use the Kullback
implicit stepping problem: find a point f˜such that entropy (39) for comparison of distributions with the
same value of the macroscopic variables. Moreover, the
m(f˜)=M, (Π∗ 1)(Θ (f˜))=0. (44) quadratic approximationto (39) will not violate second-
h
−
order accuracy, and does not require the solution of a
This antidiffusion problem is a proper two-point bound- transcendental equation.
aryvalue problem. Inafinite-dimensionalspacethe first The viscosity coefficient is proportional to τ and sig-
conditionincludes N independent equations(whereN is nificantly depends on the chain construction: for the se-
the number of independent macroscopic variables), the quence(29)wehaveτ =(1 β)h/β,andforthesequence
−
second allows N degrees of freedom, because the values ofsteps(43)τ =2(1 β)h. Forsmall1 β thelatergives
− −
ofthemacroscopicvariablesatthatendarenotfixedand around two times larger viscosity (and for realisation of
Θ (f˜) could be any point on the quasiequilibrium mani- the same viscosity we must take this into account).
h
fold). Shooting methods for the solution of this problem How can the coupled steps method (43) fail? The
looks quite simple: methodcollectsallthe highordererrorsintodissipation.