(J.-P. Rivet, J. P. Boon) Lattice Gas Hydrodynamic
(J.-P. Rivet, J. P. Boon) Lattice Gas Hydrodynamic
(J.-P. Rivet, J. P. Boon) Lattice Gas Hydrodynamic
EDITORS
G. M. Zaslavsky, R. Z. Sagdeev,
D. A. Usikov and A. A. Chernikov
1. Weak chaos and quasi-regular patterns
J. L. McCauley
2. Chaos, dynamics and fractals: an
algorithmic approach to deterministic
chaos
K. Nakamura
3. Quantum chaos: a new paradigm for
nonlinear dynamics
C. Beck and F. Schlogl
4. Thermodynamics of chaotic systems:
an introduction
P. Meakin
5. Fractals, scaling and growth far from
equilibrium
R. Badii and A. Politi
6. Complexity - hierarchical structures
and scaling in physics
H. Kantz and T. Schreiber
7. Nonlinear time series analysis
T. Bohr, M. H. Jensen, G. Paladin and A. Vulpiani
8. Dynamical systems approach to turbulence
P. Gaspard
9. Chaos, scattering and statistical mechanics
E. Scholl
10. Nonlinear Spatio-Temporal Dynamics and
Chaos in Semiconductors
J.-P. Rivet and J. P. Boon
11. Lattice Gas Hydrodynamics
Lattice Gas
Hydrodynamics
J.-P. Rivet
Observatoire de la Cote d'i4zur, France
and J. P. Boon
Universite de Bruxelles, Belgium
CAMBRIDGE
UNIVERSITY PRESS
CAMBRIDGE UNIVERSITY PRESS
Cambridge, New York, Melbourne, Madrid, Cape Town, Singapore, Sao Paulo
Published in the United States of America by Cambridge University Press, New York
www.cambridge.org
Information on this title: www.cambridge.org/9780521419444
© Cambridge University Press 2001
Preface xiii
vn
viii Contents
... Feynman told us to explain it like this: We have noticed in nature that the
behavior of a fluid depends very little on the nature of the individual particles in
that fluid. [...] We have therefore taken advantage of this fact to invent a type
of imaginary particle that is especially simple for us to simulate. This particle is
a perfect ball bearing that can move at a single speed in one of six directions.
The flow of these particles on a large enough scale is very similar to the flow of
natural fluids.
W.D. Hillis, Physics Today, February 1989
The story of lattice gas automata started around 1985 when pioneering studies
established theoretically and computationally the feasibility of simulating fluid
dynamics via a microscopic approach based on a new paradigm : a fictitious
oversimplified micro-world is constructed as an automaton universe based not
on a realistic description of interacting particles (as in molecular dynamics),
but merely on the laws of symmetry and of invariance of macroscopic physics.
Imagine point-like particles residing on a regular lattice where they move from
node to node and undergo collisions when their trajectories meet at the same
node. The remarkable fact is that, if the collisions occur according to some
simple logical rules and if the lattice has the proper symmetry, this automaton
shows global behavior very similar to that of real fluids. Then, to observe the
fluid-like behavior of the automaton, we look at the lattice from a distance and,
to measure its properties, we perform measurements over averaged quantities
to obtain for instance the fluid velocity field or to visualize its stream lines (see
Figure 1).
At the beginning, most of the effort was invested into lattice gas hydrody-
namics through theoretical developments recovering Navier-Stokes equations,
xin
Preface
Figure 1 Beyond some threshold value of the Reynolds number, the flow
behind an obstacle produces a non-stationary wake, the Benard-von Karman
vortex street illustrated in (a) by a picture taken by Prandtl and Tietjens in
1934; (b) shows the result of a lattice gas simulation (performed on the
dedicated machine 'RAP-T, see Chapter 10) under conditions similar to those
of the laboratory experiment. In both (a) and (b), the Reynolds number is of
the order of 200. ((a) In An Introduction to Fluid Dynamics, G.K. Batchelor,
Cambridge University Press, 1967; (b) J.P. Rivet, 1987.)
Preface xv
models for relatively high Reynolds numbers and transition to turbulence, multi-
species and multi-phase models for low Reynolds number complex flows (with
or without surface tension), boundary layers, lattice Boltzmann approach to
turbulence, thermal lattice gases, reactive systems. Despite an apparent diver-
sity, lattice gas automata developed for generic or specific physical problems
all share a common basis, conceptually and operationally. The reason can be
found in the following considerations.
Often systems with a large number of degrees of freedom exhibit macroscopic
behavior where the details of the microscopic dynamics are relatively unimpor-
tant. This feature opens the way to a generic description: systems with different
microscopic characteristics can be described at the macroscopic scale with a
generic set of equations where the specific nature of the constituents of the
system enter only a restricted number of coefficients at the level of constitutive
relations. Since the macroscopic description does not provide any content to
these coefficients, their specificity can be (at least partly) eliminated when the
equations are properly recast into non-dimensional form. Then the global behav-
ior of the system depends on a limited number of universal control parameters
where the microscopic nature of the elementary constituents does not appear
explicitly. The macroscopic level of analysis therefore provides a phenomenolog-
ical description where (i) the complexity of identifying the microscopic degrees
of freedom is bypassed, and (ii) many different physical phenomena are recog-
nized to belong to a limited number of classes. Classical fluid mechanics offers
a striking example through Reynolds's dynamical similarity law (1883). Never-
theless, it is obvious that in this process of reasoning the connection between
the phenomenology and the underlying microscopic mechanisms has been lost.
Consequently, the macroscopic level of description which uses average quantities
prevents analysis of, for instance, how large-scale phenomena can be triggered
by local and/or transient deviations from averages, that is, fluctuations. This
is one of the objects of the statistical mechanical approach which establishes
the microscopic basis of the properties of many-body systems. But we have
observed that it is not necessary to have a complete knowledge of the details
of the microscopic interactions to understand how macroscopic phenomena
emerge. It will appear as logical conduct to develop a microscopic approach to
macroscopic phenomena by modeling the microscopic dynamics by means of
a simplified description provided the basic requirements of fundamental physics
are correctly incorporated, essentially - although further requirements may be
in order - the conservation laws and the symmetry properties.
Our description uses a class of simple discrete models, where space and time
are incremented by finite elementary units, as in the physical picture that we
presented where point particles move step by step on a regular lattice. A
logical question to ask is: why are we interested in such discrete models? The
fact is that despite their simplicity at small scale (at the microscopic level),
xvi Preface
these discrete systems can show extremely complex behavior at larger scales
(the mesoscopic and macroscopic levels), and this global behavior can be very
similar to what is observed in natural phenomena (see Figure 1). Now many
of these natural phenomena are very complicated (for instance the diversity
of forms that a fluid can exhibit in the turbulent regime), and often the best
one can do is to offer a phenomenological analysis of what is observed. Yet
the objective is to try to understand how things come about by providing
an analysis in terms of the elementary processes that are responsible for the
emergence of large-scale phenomena, a formidable task if one starts with a
realistic description of all the basic interactions that govern the underlying
microscopic dynamics.
The philosophy that constitutes the conceptual basis of the lattice gas au-
tomaton (LGA) approach is that a microscopically simple system which can
exhibit complex macroscopic behavior in accordance with observed phenomena
should contain, at the elementary level, the essentials which are responsible for
the emergence of complexity, and thereby can teach us something about the
basic mechanisms from which complexity builds up.
As a matter of fact, when LGAs were initially developed around 1985, the
goal was to produce simplified models that could provide new computational
accessibility to complicated problems in fluid dynamics such as 2-D and 3-D
turbulence. While earlier studies (in the seventies) had considered even simpler
lattice models for basic problems in statistical mechanics, subsequently a variety
of models were constructed for more involved physical applications: hydrody-
namic instabilities, flows in porous media, phase transitions, reactive systems,
etc. All these models share a common structure where point particles move
along the links of a regular lattice and interact on the nodes through collisions.
The key point is that these simple basic ingredients (along with some symmetry
requirement for the lattice structure) suffice for the emergence of complex global
behavior from elementary processes. Furthermore, the lattice gas automaton ex-
hibits two other important features: (i) it usually resides on a large lattice, and
so possesses a large number of degrees of freedom; (ii) its microscopic Boolean
nature (which we shall describe in detail in the next two chapters) combined with
the (generally) stochastic rules which govern its microscopic dynamics, result in
intrinsic fluctuations. Because of these spontaneous fluctuations and because of
its large number of degrees of freedom, the lattice gas can be considered as a
'reservoir of thermal excitations' in much the same way as an actual fluid (see
Figure 2).
The principal goal of this book is to carry two messages. The first is to
show how an automaton universe with simple microscopic dynamics can exhibit
macroscopic behavior in accordance with the phenomenological description of
classical physics. The second is to establish that the correlations of the lattice gas
intrinsic fluctuations capture the essentials of actual fluctuations in real fluids.
Preface
The two aspects are of course closely connected. Therefore our emphasis is on
the statistical mechanics of lattice gas automata as a microscopic approach to
the description of fluid systems, and in this sense we proceed very much along
the lines of the general philosophy of classical statistical mechanics.
After this simple presentation of the lattice gas picture and of the ideas
underlying the philosophy of the lattice gas automaton approach, we introduce
the basic notions of cellular automata and of lattice gas automata from the
physical point of view and from the mathematical point of view, and we
discuss the complementary aspects of the two viewpoints (Chapter 1). With
these basic elements, we construct the general formalism of the microscopic
dynamics of the lattice gas automaton leading to the microdynamic equations
and the corresponding microscopic properties (Chapter 2). The formalism is then
illustrated through the most used LGA models which are discussed in Chapter 3.
At this point we have all the material necessary to establish the basic statistical
mechanics of the lattice gas leading to the most important mathematical object
of our subject: the lattice Boltzmann equation, and its consequences, the H-
theorem and the equilibrium state (Chapter 4). A discussion of the equilibrium
properties of the lattice gas follows naturally providing interesting analogies
with the static correlation functions of classical statistical mechanics.
The macroscopic dynamics of the lattice gas is developed following two com-
plementary paths. In Chapter 5, we proceed according the classical Chapman-
Enskog method to derive the lattice gas hydrodynamic equations and the ex-
pressions for the transport coefficients. These results show how the connection
can be established between the lattice gas macrodynamics and the equations
of classical fluid dynamics; we also discuss the problem of Galilean invariance,
which calls for special attention. The hydrodynamic regimes which emerge from
the limits of the macrodynamic equations are explored in Chapter 9; in par-
ticular, in the incompressible fluid limit, we recover exactly the Navier-Stokes
equations of classical fluid dynamics.
When we consider small deviations from local equilibrium, we can use the
linearized Boltzmann equation to describe the dynamics of fluctuations, and
from the natural separation into fast and slow variables, we obtain the linearized
hydrodynamic limit (Chapter 6), which sets the stage for the investigation
of hydrodynamic fluctuations whose correlation functions are computed in
Chapter 7. In the hydrodynamic limit, we retrieve the results of the Landau-
Placzek theory, and we show that the lattice gas dynamic structure factor exhibits
the same structure as the power spectrum of density fluctuations in real fluids.
We then consider the lattice gas dynamics beyond the hydrodynamic regime -
that is at shorter wavelengths - where the LGA approach offers the interesting
possibility of exploring all the modes of the system up to the kinetic regime.
These results are also extended to the two-component system (the 'colored' lattice
gas) to incorporate mass diffusion. Larger deviations for systems arbitrarily far
from equilibrium are considered in Chapter 8 where, using a projector technique,
we derive the full non-linear hydrodynamic equations and the general form of
the Green-Kubo expression for the transport coefficients.
Numerical simulations are an important contribution to our subject: they
Preface xix
provide tests of validity of the LGA theory and, in cases such as those illustrated
in Figures 1 and 2, the results of lattice gas simulations can be compared with
laboratory experimental results; but most importantly they can also offer access
- theoretically or numerically - to regimes which cannot be easily reached by
other methods. Chapter 10 presents examples of lattice gas simulations which
are illustrated and discussed for various hydrodynamic flows, along with the
essentials of the procedure for computer implementation. Finally, we close with
an epilogue in the form of a guide for further reading (Chapter 11) indicating
topics and areas of research where lattice gas automata offer, technically or
conceptually, a particularly interesting and useful approach.1
The ideas expressed in this book and the material presented here reflect
our own choices and views. The materialization of these ideas has been influ-
enced in various ways by the interaction, the stimulation and the discussions
we had over the past fifteen years with our colleagues and friends. To all of
them we express our appreciation, with special thanks to David Dab, Matthieu
Ernst, Uriel Frisch, Patrick Grosfils, David Hanon, Michel Henon, Dominique
d'Humieres, Ray Kapral, Pierre Lallemand, Alain Noullez, Dan Rothman, Al-
berto Suarez, Sauro Succi, Olivier Tribel, Jorg Weimar, and Stephane Zaleski.
JPR acknowledges support by the Centre National de la Recherche Scientifique
(CNRS, France) and JPB acknowledges support by the Fonds National de la
Recherche Scientifique (FNRS, Belgium).
1
As much as possible, we avoided disrupting the course of the main text by deferring
mathematical details and technical developments to the appendix.
Chapter 1
Basic ideas
When the lattice gas was introduced in statistical physics around 1985 (see Frisch,
Hasslacher and Pomeau, 1986), it was originally constructed as a physical model
for hydrodynamics. In fact, the concept of the lattice gas is as much a physical
concept - and we shall indeed start with intuitive physical ideas - as it is a
mathematical concept, as a more formal definition can also be given. We first
present the point of view of the physicist (Section 1.1), then we describe the
lattice gas automaton from the mathematical viewpoint (Section 1.2), and in
Section 1.3 we discuss the two aspects.
conservations are necessary, but not sufficient conditions in order that a lattice
gas exhibit physically realistic collective large-scale dynamics: an appropriate
crystallographic space group for the lattice2 is another necessary condition for
macroscopic dynamics with correct continuous invariances. These points will be
emphasized in later chapters.
Another important feature is the 'exclusion principle': two absolutely indis-
tinguishable particles cannot reside simultaneously on the same node. However,
two particles with different velocities can coexist at a node as well as parti-
cles with the same velocity as far as they can be distinguished by some other
property (mass, color, label, etc.). A consequence of the exclusion principle is
the following: provided there exists only a finite number b of distinguishable
kinds of particles, the occupation state of a node is completely determined by
b Boolean variables, each variable encoding the presence or the absence of a
particle of each kind on the node. This feature is of considerable importance
for the implementation of lattice gases on special-purpose machines as well as
on general purpose computers.
So the physicist would view a lattice gas as a universe of fictitious zero-
dimensional particles undergoing one-dimensional displacements on the links
and colliding on the nodes of a regular Bravais lattice embedded in a D-
dimensional space. But the Boolean encoding of the state of a lattice gas leads
naturally to a mathematically more formal point of view.
We shall first introduce the concepts of 'finite automata' and 'cellular automata'.
alphabet. The above definition prescribes a deterministic evolution rule for the
internal state of the automaton. Nevertheless, a non-deterministic automaton
also fits the above definition if its non-deterministic evolution rule can be
considered as a deterministic rule depending on a finite number of additional
input symbols belonging to finite alphabets, produced by 'external' (pseudo-)
random generators.
Since the internal state of a finite automaton can change, the automaton
undergoes some kind of evolution, and there must be an underlying notion of
'before' and 'after': a 'past' and a 'future'. However, these primitive notions do
not really imply a temporal structure for the automaton, since the concept of
'time interval between events' is not included in the definition.
(i) The individual automata are tied geometrically to the nodes of a regular
Bravais lattice embedded in a Euclidean space of dimension D. This relation
between individual automata and lattice nodes yields a one-to-one mapping
between the set of all the individual automata of the cellular automaton
and a connected part of the Bravais lattice. Because of this one-to-one
mapping, the individual automata can also be called 'nodes'. The nodes
are labeled by their position vector r*; the star subscript indicates that the
vector takes only discrete values.
(ii) The elementary evolution process of the cellular automaton is repeated at
regularly spaced discrete times separated by a time increment At (the time
step), which is usually taken equal to unity by a proper choice of the time
scale.
(hi) The number of possible internal states of any individual automaton is
2b, where b is an integer. Any possible internal state of an individual
automaton can thus be represented by b Boolean variables, which will
be called 'channels'3 because they correspond to 'communication channels'
between nodes, as we shall see below. The b channels on any node are
labeled by an integer i ranging from 1 to b or from 0 to b — 1. Since all
the individual automata are taken identical, the number b of channels per
node and their labeling are node-independent.
(iv) The elementary evolution process of the cellular automaton that occurs at
each time step is a sequence of two distinct phases:
(a) The collision phase: during this phase, the internal state of each indi-
vidual automaton evolves according to a purely local, node-independent
rule which can be non-reversible and non-deterministic. More precisely,
the post-collision state of an individual automaton depends only on its
own pre-collision state and on afinitenumber of random input symbols
for non-deterministic models (see comments on local versus non-local
collisions in Section 1.3).
(b) The propagation phase: during this phase, the information present in the
automata's 'memory' (the post-collision state) is shifted from node to
node without loss or gain, according to a node-independent, reversible
and deterministic rule. In practice, each channel index i between 0
and b — 1 is associated to a unique lattice vector c,-, such that the
Boolean value of channel i on any node r* is transferred during the
propagation phase to channel i of node r* + c*. The vectors c,- are called
the 'velocity vectors', and must verify the following conditions: (i) the
set of velocity vectors remains globally invariant under the action of the
3
Channels are sometimes called 'cells', but the word 'cell' may be misleading because of
possible confusion with the crystallographic notion of (primitive) cell of a lattice.
1.3 Comments 5
1.3 Comments
1.3.1 The velocity vectors
We first justify the name 'velocity vectors'. If the time scale is chosen in such
a way that the time step At is equal to unity, then the information (particle)
present in channel i at node r* goes to node r* + ct in one unit of time. The
CjS can thus be viewed as the microscopic velocity vectors of fictitious particles
occupying channels i.
All velocity vectors are not necessarily distinct: lattice gas models with mul-
tiple channels can be designed to model, for instance, gas mixtures or reaction-
diffusion systems where distinguishable particles representing different 'chemical
species' have the same kinematics (i.e. the same propagation rules). In addition,
some of the velocity vectors can be zero: this is the case for models with rest
particles which stay on a node during the propagation phase.
The velocity vectors, which must have the same local symmetries as the
lattice, are also constrained to include D generating vectors of the Bravais
lattice. This guarantees that any two nodes of the lattice can be connected by
a finite chain of velocity vectors, and thus, that the lattice cannot separate into
sub-lattices with independent evolutions. Note that for a given Bravais lattice,
there may exist several sets of D lattice vectors that generate the whole lattice,
and several sets of velocity vectors compatible with the symmetries of the lattice.
To illustrate this point, take for instance the three-dimensional body-centered
lattice ('bcc' lattice) with unit lattice constant (length of the spatial periodicity).
This Bravais lattice can be generated by the set £f\ containing the three lattice
an
vectors {\>—\,—\),{—\i\i—\)> d (i> i> \)- The smallest set of velocity vectors
that can be built around Sf\ contains eight vectors which can be constructed
by applying all isometries of the local invariance group of the bcc lattice, to
all three vectors of Sf\. These eight vectors with modulus *f have the form
(^15^25^35). Here, (^1,62,^3) are any of the eight possible combinations of'+'
4
The point group or local invariance group of a lattice is the subgroup of the space
group of the lattice which leaves individual nodes unchanged; see Ashcroft and Mermin
(1976).
5
This definition should not be considered as universal: with the same basic ideas, one can
relax, modify or add constraints to construct different definitions of lattice gases.
6 1 Basic ideas
and '—'. Now, the bcc lattice can also be generated by the set &\ containing the
three vectors (1,0,0), (0,1,0) and (5,5,5). The smallest set of velocity vectors
that can be built around ^2 contains fourteen vectors that include the eight
vectors with modulus ^ , plus six vectors with unit modulus: (±1,0,0), (0, ±1,0)
and (0,0, ±1). This example shows that with one single Bravais lattice, one can
construct several lattice gas models with different propagation rules and thus
different complexities. It motivates the introduction of the concept of 'minimal
model' for a given Bravais lattice. A minimal model built on a given Bravais
lattice is a model which has the smallest possible number b of channels per
node compatible with the local symmetries of the lattice. There may be several
minimal models associated to the same Bravais lattice, since they can differ
by their collision rules, but also by the 'physical' properties attached to the
channels.
correct provided the collision phase can 'break' particles with multiple mass in
such a way that it produces no spurious conservations.
Note that n distinguishable particles require n Boolean quantities, whereas n
indistinguishable particles only require p ~ log 2 n Boolean quantities, provided n
is of the form 2P — 1. The problem of distinguishable versus indistinguishable
particles occurs, for example, for some three-dimensional lattice gas models (see
Section 3.7).6
Indeed, for the latter, the set of nearest neighbors is different for 'odd' and 'even'
nodes. No node-independent connection scheme involving nearest neighbors can
be designed on the honeycomb lattice.
One can classify exhaustively and enumerate all possible Bravais lattices in
given dimensions, according to their symmetry groups: there exist only five
different Bravais lattices in two dimensions, and fourteen in three dimensions,
as first observed by Auguste Bravais in 1845 (see Cauchy, 1851). Notice that the
constraint of using a Bravais lattice is not too severely restrictive, since several
lattice gas models can be built on one single Bravais lattice by changing the
connection scheme and the collision rules. Moreover, discrete models within the
spirit of lattice gases, but involving non-regular or non-Bravais or quasi-lattices
can exist, but will not be considered here, as we keep with the definition of
lattice gases given in Section 1.2.3.
may depend on the phase where the iteration sequence is stopped (see Chapter 4).
Care must be taken in using proper averaging procedure.
8
A different point of view, discrediting the concept of cellular automata for lattice gases,
is presented by Henon (Henon, 1989).
Chapter 2
10
2.1 Basic concepts and notation 11
2.1.3 Observables
In order to establish the connection between purely abstract automata and
concrete gas models, physical properties (mass, momentum, energy, etc.) must
be associated with each channel of each node of the lattice. We emphasize
that such physical quantities (observables) are not associated with the particles
2.1 Basic concepts and notation 13
themselves, but with the channels: particles loose their identity during collisions,
whereas channels do not. Nevertheless, for simplicity and intuition, we shall use
terms such as 'the mass of particle f or 'the momentum of particle f rather
than 'the mass associated with channel number i if occupied' or 'the momentum
associated with channel number i if occupied'. Mass, momentum and energy are
examples of observables with physical content; yet, the concept of observable, as
we define it below, is in principle not restricted to quantities with an equivalent
in the physical world.
It is systematically assumed that the amount of a given observable assigned
to a channel depends only on the channel index, and not on the position of the
node: if occupied, channel i of any node always carries the same amount of the
observable.
Under the above assumption, an observable is defined as a collection of
b values (qi9 i = l,...,fo), indicating the amount of the underlying observable
quantity which is present in the channel labeled by the index i9 when this channel
is occupied. The qfi can be scalars, vectors or more generally tensors of arbitrary
order.
The 'value of the observable measured at node r*' is the total amount of the
observable quantity present at node r*: X^=i 4ini(r*)- It is called the 'microscopic
density per node' of the observable or simply its 'microscopic density', and is
denoted by N*(r*) in the most general case:
1=1
2.1.3.1 Example 1
The collection g = ( l , l , l , . . . , l ) o f b quantities equal to 1 is a scalar observable
called the 'particle number observable', as its value measured at node r* is
J2 n,-(r*) which is clearly the total number of particles present at node r*.
2.1.3.2 Example 2
Denoting by m,- the individual mass assigned to any particle present in channel i,
the collection q = (mi,..., nib) is a scalar observable called the 'mass observable'.
Its value ^ m ;fy( r *)> measured at node r*, is the total mass present at node r*.
If the individual masses m,-, (i = 1,..., b) are all equal to a single value m (which
can be taken equal to 1 by a proper choice of the mass scale), then the mass
observable is equal to the particle number observable. This is the most frequent
case.
2.1.3.3 Example 3
Denoting by pia((x = 1,...,D) the components of the individual momentum as-
signed to any particle present in channel i, the collection qa = (pi a ,... ,p&a), (a =
1,...,D) is a vector observable called the 'momentum observable'. Its value
qa = 5^Pia»i(r*), 2 measured at node r*, is the total momentum present at node
r*. For many models, the individual momentum p, of channel i is chosen equal
to niiCi, which is physically consistent with the interpretation of the C;S as mi-
croscopic velocities. This choice also guarantees that the microscopic density of
momentum is just the microscopic flux of mass.
2.1.3.4 Example 4
Denoting by e\ the individual energy assigned to any particle present in channel i,
the collection q = (e\,..., e&) is a scalar observable called the 'energy observable'.
Its value J2eini(r*)> measured at node r*, is the total energy present at node r*.
If particles carry only kinetic energy, it is reasonable to take e\ equal to p?/(2mj).
does not explicitly3 depend on r*, then the generalized observable is said to be
'homogeneous'. When the value of the observable measured at node r* only
depends on thefy-sat this very node r*, then the generalized observable is said
to be 'local'.
The following simple examples should provide a more physical picture of
generalized observables.
This generalized observable is linear and homogeneous, but non-local, since its
value, measured at node r*, depends on the Boolean field at nodes r* + r*', with
x* e «/. The resulting value can be understood as the population of channel
number 1, space-averaged over the domain J surrounding node r*.
This generalized observable is linear and local, but non-homogeneous, since its
value depends explicitly on r* through the scalar potential U(u).
The cellular automaton rule governing the evolution4 of a lattice gas can be
viewed as a logical equation expressing the Boolean field at time U + 1 from
the Boolean field at time U. When written in algebraic form, the resulting
equation is called the 'microdynamic equation' of the lattice gas, as it governs
the microscopic evolution (microdynamics) of the lattice (Frisch, d'Humieres,
Hasslacher, Lallemand, Pomeau and Rivet, 1987). The microdynamic equation
is an exact equation which contains the complete dynamics of the lattice gas.
3
The value of q depends implicitly on the space variable r*, through the dependence in
the Boolean field.
4
The time scale is chosen such that the time increment of the cellular automaton is one.
16 2 Microdynamics: general formalism
For the non-deterministic case, the transition between times U and U + l involves
random choices at each node. The microdynamic equation then reads:
This definition of the propagation operator raises the problem of finite size
lattices. Indeed, if the subset <£ of the Bravais lattice under consideration is
finite, the node r* + c, may be outside j£f, even if the node r* is inside.
One solution to this problem is to introduce a periodic wrapping in such a
way that any particle driven out of S£ be re-injected elsewhere in if. More
precisely, we need a finite domain if' of the underlying Bravais lattice which
2.2 The microdynamic equation 17
contains all nodes of <£ and which can perfectly tile the full Bravais lattice
by successive lattice-preserving translations. So there must exist a subgroup of
the lattice-preserving translations group such that the nodes 5 in $£' and their
images by the elements of this subgroup cover the full infinite Bravais lattice
without any gap or overlap (see Figure 2.1 (a)). It is then possible to wrap $£'
on itself and to 'connect' opposite boundaries in such a way that j§?' become a
topological torus (see Figure 2.1(b)). In terms of particle motions, this wrapping
leads to periodic boundary conditions6 at the boundaries of if'.
Periodic
boundary conditions
For a given lattice gas model and a given domain j£f, there can exist several
different periodic boundary conditions. Indeed, the same domain <£ can lead to
several perfect tilings with different tile shapes, even if one restricts the choice
to tiles which contain a minimal number of nodes (see Figure 2.2(a)). Moreover,
the same tile shape can lead to several different tilings with different translation
subgroups (see Figure 2.2(b)).
We emphasize that this periodic wrapping is applied only to the propagation
operation; it does not imply that the macroscopic variables satisfy periodic con-
ditions in the usual sense. Indeed, special collision operators can be applied to
the boundary nodes in such a way that the boundaries of <£ behave macro-
scopically as solid walls, or as mass and/or momentum injectors for example.
These special collision operators can be applied as well to the boundaries of
$£ to produce particular macroscopic boundary conditions, and/or to localized
parts inside j£?, to simulate for instance fixed solid obstacles of almost any
arbitrary shape. Some of the most useful special collision rules are described
and discussed in Section 2.4.
Figure 2.2 The same domain <£ can lead to several different tilings with (a)
different tile shapes and/or (b) different translation subgroups. The
two-dimensional triangular Bravais lattice is used to illustrate the point.
2.2 The microdynamic equation 19
transition probabilities, the matrix elements are real numbers between zero and
one, and satisfy the normalization condition:
A(s^) = l9 Vsey. (2.9)
s'ey
Under the above assumption, one can write the collision operator in an explicit
form that involves the collision matrix indirectly:
where Sj = 1 — Sj and H;- = 1 — rij. The rijS and n ; s are taken at node r*.
The £(s^>s') s are the elements of a 2b x 2 b square matrix of Bernoulli random
variables (values 0 or 1) whose averages are the A(s-+sf)s, and such that
selects, among all possible pairs of Boolean states (5, s'\ the one for which s is
the pre-collision state of node r*, and s' is the actually realized post-collision
state. Equation (2.10) assigns to ^(n)(r*) the value of this selected post-collision
state s'.
The above formalism is applicable to both deterministic and non-deterministic
models. For deterministic models, the matrix elements A(s^>sf) are either 0 or 1,
and the £(s->s')s are deterministic and thus simply equal to the A(s^>s')s.
20 2 Microdynamics: general formalism
(s,s')ey2 ;=1
where the rijS on the r.h.s. are implicitly taken at node r* and at time U. An
alternative form, strictly equivalent to (2.12), reads:
b
The equivalence between (2.12) and (2.13) follows immediately from the identity:
b
(s,s')ey2 7=1
which is a consequence of (2.11). The virtue of the alternative form (2.13) is that
it distinguishes 'effective collisions' with s' =/= s from 'passive collisions' which
leave the Boolean state unchanged.
The microdynamic equation plays for lattice gases the role of Hamilton's
equations for a system of interacting particles; it involves no approximation and
constitutes the starting point of all further theoretical developments.
As for real gases, the macroscopic (large-scale) dynamics of lattice gases depends
crucially on their microscopic properties. We now define some of the most
important microscopic characteristics that a given lattice gas model may or may
not possess. The consequences of these microscopic characteristic properties on
the large-scale macroscopic behavior will be discussed in Chapters 4 and 5.
sey sey
2.3 Microscopic properties of a lattice gas 21
which is the most commonly used expression for the semi-detailed balance
condition.
The semi-detailed balance will be crucial in deriving the macroscopic theory
of lattice gases. In particular, it affects deeply the nature of equilibrium states
and the validity of the Boltzmann approximation.
2.3.2 Duality
For any Boolean state s, the dual state s is obtained by exchanging particles and
'holes' (occupied and unoccupied channels) or Os and Is in s. If s is considered
as a collection of b logical variables, then the dual state is s, the logical bitwise
negation of s.
A lattice gas model is said to be 'self-dual' if the dynamics of holes is, at least
statistically, identical to the dynamics of particles. The propagation operator
moves information without any loss or gain, and is thus self-dual. Consequently,
a necessary and sufficient condition for a model to be self-dual is that the
collision matrix A(s^sf) satisfy the relation:
V(s, s') e y2, A(s->sf) = A(s->sf). (2.18)
5, s') e y\
This equation, which must hold for all discrete times U and for all nodes in jSf,
expresses that the amount of the collisional invariant q present at time U at node
r* is conserved by the collision operator, and distributed over the b neighboring
nodes r*+c; by the propagation operator.8 Equation (2.20) displaying collisional
invariance will be explicitly used as a starting point for the macroscopic theory
of lattice gases.
For instance, the triplet formed by any collisional invariant, the whole lattice
J&f, and any integer k is a static geometrical invariant, since the propagation
operator only redistributes the total amount of the observable initially present in
JSf, without any loss or gain. Of course, this only holds when periodic conditions
are applied (with no boundaries elsewhere).
2.3.3.3 Moving geometrical invariant
In order to define a moving geometrical invariant, we need a triplet (^,«/,T),
where q is an observable, J> a sub-lattice of JSf, and T a lattice-preserving
translation which does not leave J globally invariant. It is then assumed that,
for any initial configuration of the lattice, the amount of observable present in
J be transferred completely to T(«/). Re-expressed in terms of the Boolean field,
this reads:
b b
2^(r*>U + 1
)=EE^(r*,**). (2.22)
2.3.3.4 Conserved quantities
A conserved quantity is simply a static geometrical invariant with k = 1. It is a
doublet formed by an observable q and a sub-lattice J> of if, such that the value
of the observable summed over all nodes in J remains unchanged under the
action of the full evolution operator, for any initial configuration of the lattice.
A collisional invariant, summed over the whole lattice is always a conserved
quantity (except if special conditions are applied, e.g. local injection or removal
of particles, biased velocity field).
2.3.4 G-invariance
The property of G-invariance is quite intuitive. It characterizes lattice gas mod-
els whose microscopic structure (propagation rules, collision rules, observables)
is compatible with the local symmetries of their underlying Bravais lattice.
G-invariant lattice gases are those whose physical structure has the same sym-
metries as the Bravais lattice.
2.3.4.1 Definition
A lattice gas model is said to be G-invariant whenever, for any isometry g in the
crystallographic point group (local invariance group) G of the Bravais lattice,
there exists a unique index permutation g acting on {1,2,...,b} such that:
24 2 Microdynamics: general formalism
g t e ) = 4g(/)> i = 1,...,&,
c15 CM
c21
c2
C22
\ L/ / c,
1\ clf / c, c
C}
c,
Cu
—^ —^ 9
\1
c4 c17
< ^
> c 2l
/ t
\ C2<
c23 c18 Cn cn
Figure 2.3 (a) A lattice gas model (the HPP model) with only one class of
channels, (b) A more complex model with five classes.
element of the orthogonal group Op in dimension D, that is, the group of all
distance-preserving linear transformations acting on R2* (see Tung, 1985). There
exists an equivalent and more intuitive criterion to characterize the isotropy of
a tensor: a tensor of order n is isotropic if, and only if, its contraction with
n arbitrary 'test' vectors of R^ can be expressed exclusively in terms of inner
products of these vectors. With this criterion, it is obvious that odd-order tensors
cannot be isotropic, unless they are null. It is also easy to convince oneself that
a second order tensor T aia2 is isotropic if, and only if, it is a multiple of the unit
tensor 5aiai (Kronecker symbol). In the same way, a fourth order isotropic tensor
T^WWA must be a linear combination of 5^3^, <5aia3(5a2a4, and <5aiaA2a3- More
generally, an nth rank isotropic tensor (with n even) must be a linear combination
of all products of the form <5aff(1)a<r(2)... <5a<y(fI_1)a<r(B), where a is any permutation on
(1,...,*).
If, in addition, the tensor is fully symmetric (i.e. invariant under any per-
mutation of its indices), then there exists a second isotropy criterion which is
equivalent to the previous one, but technically much simpler to use: A fully
symmetric nth order tensor T (n) is isotropic if, and only if, its n-fold contraction
with one single arbitrary test vector x of R D is proportional to the nth power
of the modulus ||x|| of the test vector, that is, if, and only if, there exists a real
number a such that:
sense, but its local geometry (the c/s) may however yield crystallographic isotropy
up to some order.
Crystallographic tensors built on the c/ vectors of lattice gas models are
fully symmetric by construction. Hence, the second isotropy criterion provides
a simple algebraic method to check the crystallographic isotropy properties of
a given lattice gas model. Indeed, consider the nth order crystallographic tensor
j{n) = ^ = 1 c i a i ...Cian of a lattice gas model; if n is odd, the isotropy test is
simple: isotropy holds if, and only if, T^ is zero. If n is even, then the full
contraction of T(w) with an arbitrary vector x reduces to X^=i(c* ' XT which is a
polynomial function with even order of the coordinates (xi,..., xp) of the vector
x. The test for crystallographic isotropy is then to check whether this function
is proportional to (]Ca=i x«)2. Such algebraic manipulations, involving only the
use of multinomial expansions, can be carried out by hand at least for orders
lower than or equal to 6 or 8. (The use of a symbolic manipulator can simplify
the task.)
2.3.6 Irreducibility
Irreducibility is a symmetry-related property which is not as natural as the
notions defined previously. However, some of the calculations presented in
Chapters 4 and following, are dramatically simplified if this property is verified.
We thus wish to spend some time to define and visualize clearly this notion, and
to explain its consequences.10
2.3.6.1 Definition
A lattice gas model is said to be 'irreducible' if, and only if, for any channel index
i= 1,..., ft, the invariant subspace £; of the subgroup Gt of all lattice-preserving
isometries that also preserve c,-, is linearly generated by c,-. In other words, any
vector u preserved by all the isometries in G which also leave c,- unchanged,
must be a multiple of c,-.
We shall try to justify the name 'irreducibility' given to this property. Let us
consider the subgroup Gt and the hyper-plane 11/ of dimension D — 1, orthogonal
to c/. We construct the set G\ collecting the restrictions to 11/ of all elements in G/.
The irreducibility property simply states that G\ is an irreducible representation
of Gt on n*.11
2.3.6.2 Examples
To illustrate the notion of irreducibility for a lattice gas, we use once more
the simple HPP and extended HPP models presented in Section 2.3.4. Let us
10
On first reading, the reader can omit the details of the present section and just keep in
mind its conclusions and consequences.
11
For basics on irreducible representations of groups, the reader is referred to text books
on group theory, for example Eliott and Dawber (1979), or Tung (1985).
2.4 Special rules 29
consider first the simple HPP model (see Figure 2.3(a) for notations). We pick
channel 1 corresponding to the velocity vector ci = (1,0). The subgroup G\ of G
leaving ci unchanged contains only the identity (Ro) and the mirror symmetry
Sy. The only vectors in R 2 that remain invariant under the action of G\ are
clearly proportional to ci. Since all four vectors ci to C4 are equivalent, what
is true for ci is true for all c* vectors. This completes the proof that the simple
HPP model possesses the property of irreducibility. Consider now the extended
HPP model defined in Figure 2.3(b). We pick channel 13 corresponding to the
velocity vector C13 = (2,1). The subgroup G13 of G leaving C13 unchanged only
contains the identity (Ro)- All vectors in ]RD are invariant under the action
of G13, not only those proportional to C13. This shows that the extended HPP
model does not possess the irreducibility property.
If the part of lattice i f actually occupied by the lattice gas were infinite in
all directions, then the cellular automaton rule described by the microdynamic
12
In this context, 'symmetry' means that Ti0Lp = T^ a , where Titxp are the Cartesian
components of the tensor Tj.
13
For simplicity, we denote by the same symbol g the isometry acting on vectors and on
second order tensors.
30 2 Microdynamics: general formalism
axes, at least for specular reflections. If solid obstacle conditions are applied to
the boundaries of if, for example to simulate wall-bounded flows, then some
care must be taken to avoid leakage. The layer of boundary nodes where the
obstacle collision operator is applied must be sufficiently thick to guarantee that
during the propagation phase, no particle be sent directly outside the layer of
boundary nodes. This point is discussed in Chapter 10.
There exist several ways to realize particle reflection, leading to several dif-
ferent reflection operators, and to different macroscopic effects. In the next
paragraphs, we describe some of the simplest reflections: the 'bounce-back',
'specular', and 'diffusive' reflection operators.
M7I/,
WW/,
Diffusive
Px t=> -PXr^
p
y P
Figure 2.4 The figure illustrates, in a two-dimensional case, the effect on the
total momentum of (a) the bounce-back reflection, (b) the specular reflection
and (c) the diffusive reflection. The solid arrows represent total momenta at
surface nodes. They do not represent the individual velocity vectors c,- as in
the figures of Chapter 3. Px and Py denote respectively the parallel and
orthogonal components of the total momentum P.
32 2 Microdynamics: general formalism
conserves the total mass and energy observables, but the momentum observable
is reversed.
The bounce-back reflection mode is by far the simplest one and does not
depend on the orientation of the surface. The easiest realization is to apply
the central symmetry to the state of obstacle nodes. A quite efficient numerical
strategy to implement this central symmetry is described in Chapter 10. At
the macroscopic level, the global effect of bounce-back reflections is that all
components of the mean velocity are canceled at the surface of the obstacle,
where 'no-slip' boundary conditions are thus realized.
14
A binary symmetry axis is such that a 180°rotation around the axis leaves the Bravais
lattice globally invariant.
2.5 Comments 33
2.5 Comments
We have presented the framework and the tools for the description and classi-
fication of lattice gases at the microscopic level. These microdynamic concepts
and properties have a crucial importance for the theory to be developed in later
chapters. The material presented so far may be perceived as rather abstract.
Therefore, before engaging in the systematic study of the statistical properties of
lattice gases at equilibrium and out of equilibrium, we devote the next chapter
to the description of various examples of lattice gas models in order to illustrate
the microdynamic formalism.15
15
The reader familiar with the microdynamic formalism and with lattice gas models may
want to move directly to Chapter 4.
Chapter 3
Historically, the first lattice gas model was introduced in the early seventies by
Hardy, de Pazzis and Pomeau (1973) with motivations focusing on fundamental
34
3.1 The HPP model 35
aspects of statistical physics (see also Hardy et al, 1972, 1976 and 1977). Their
model (HPP) is based on the two-dimensional square lattice with unit lattice
constant (see Figure 3.1).
The HPP model is a 4-bit model, where each node has b = 4 channels,
corresponding to the four directions of the square lattice (east, north, west,
south) labeled by integers from 1 to 4. The c,-s have unit modulus and their
Cartesian components are:
These vectors connect a node to its four nearest neighbors. So, the model's
connectivity B and the number b of channels per node are both equal to
the coordination number of the two-dimensional square lattice, namely four.
Furthermore, the four connection vectors are identical to the four velocity
vectors C*, which are all different and non-zero (see Section 2.1.1).
The masses m,- of particles are equal, and by a suitable choice of the mass
scale, their common value may be taken as one unit mass.1 The momentum p,
of each particle is equal to c,-, which is physically consistent with unit mass and
unit time step. The kinetic energy of each particle is taken equal to \, again
with physical consistency. The choice of unitary quantities for lattice constant,
time step, and particle mass completely defines a natural system of units for
microscopic mechanical quantities.
As prescribed by the general definition of lattice gases (see Section 1.2.3),
1
-I 1 l_ -4 _
l
1 1 1 1 1 1k 1
\ 1 J 1 1
1 ^1 1 T
1 1 1 1 ,. 1
-\ — 1— + -\ — i 1-
1
Figure 3.2 Collision rules for the HPP model, reduced by symmetry. The full
set of collisions can be reconstructed from the above reduced set by the
application of all the isometries of ^. The number T above the open arrow
indicates that the collisional change of configuration occurs with probability 1.
3.1 The HPP model 37
model by taking all the £(s^s')s equal to the ^l(s->s')s (since the model is
deterministic). It reads:
for i = 1,...,4:
- ^+1^+2^+3
+
Here, n; stands for (1 — nt) and the subscript i (channel index) is cyclic (modulo
4). All variables on the right hand side are taken implicitly at time U and node
r*.
Because of the simplicity of the HPP microdynamics, (3.2) can be derived
merely by intuition: the three terms on the right hand side simply express that
the collision phase can either:
We note that (3.2) is compatible with the exclusion principle, since it guaran-
tees that the nts at time U + 1 cannot be anything other than 0 or 1, provided
the n,s at time U are 0 or 1.
For simplicity, we have chosen to write (3.2) for a microscopic evolution law
where the collision phase is performed before the propagation phase. If the
reverse order is chosen, the microdynamical equation takes almost the same
form, except that each of the n,s and n^ on the right hand side must then be
taken at a different node (namely at nodes r* — c,-), instead of being all taken
at the same node r*. This minor difference is generally unimportant, since the
order of the two successive phases becomes macroscopically irrelevant when a
large number of steps are executed sequentially.2
The microdynamical equation can also be written in logical form; it suffices
to consider the n^s as logical variables with values T R U E ' or 'FALSE' (instead
of algebraic variables with values T or '0'), to obtain:
o\
1 (3.4)
m Py 0
V-i/
In addition, the HPP model has several static and moving geometrical
invariants (see Section 2.3.3), such as the x- (or y-) component of the
momentum, summed over any line of lattice nodes with equal x- (or y-)
coordinates (d'Humieres, Qian and Lallemand, 1989, 1990).
(v) The model is clearly G-invariant (invariant under the crystallographic group
of the lattice, see Section 2.3.4), and all four c,- vectors clearly belong to the
same class, as defined in Section 2.3.4.
(vi) The HPP lattice gas model has third order crystallographic isotropy (see
Section 2.3.5). As a justification, Table 3.1 gives the crystallographic tensors
of order 1 to 4, contracted with a vector x with coordinates (xi,X2)-
Clearly, the contraction of the fourth order tensor is not proportional to
Since all C;S have the same modulus 1 and all the particles have the same
mass, the kinetic energy observable is simply proportional to the mass observable.
States with the same mass will also have the same kinetic energy, so, the energy
conservation and the mass conservation are equivalent. As a consequence, the
HPP model ignores thermal effects.
The degree of crystallographic isotropy of the HPP model (namely three)
is not sufficient to produce isotropic large-scale dynamics as described by the
Navier-Stokes equations for a Newtonian fluid.3 This requires at least fourth
order crystallographic isotropy (see Chapter 5). It is actually possible to construct
models with fourth order isotropy on the square lattice. However, such models
are not minimal, in the sense of the concept introduced in Section 1.3.1, and
require a larger set of c* vectors: in addition to the HPP c,s, which must now
be endowed with multiplicity four (4 channels for each of the c;s), four vectors
must be added* connecting each node to its next-nearest neighbors (1 channel
per vector). Such 'modified' HPP models require b = 20 channels per node, and
are not homokinetic. About ten years after the introduction of the HPP model,
the 'anisotropy disease' has been cured in a less expensive way by models based
on the triangular lattice, which are discussed in the next sections.
sponding to the six directions of the triangular lattice, labeled by integers from
1 to 6. The CjS have unit modulus and their Cartesian components are:
(3.5)
C4
=(~o 1
These vectors connect each node to its six nearest neighbors. Thus, the model's
connectivity B and the number b of channels per node are both equal to
the coordination number of the two-dimensional triangular lattice, namely six.
Furthermore, the six connection vectors are identical to the six velocity vectors
Cj, the latter being all different and non-zero (see Section 2.1.1).
The FHP-1 model shares with the HPP model the following properties: The
masses m* of all the particles are equal, and by a suitable choice of the mass
scale, the common value is taken as unity. The momentum p,- of each particle
is equal to c,-; this is physically consistent with unit mass and unit time steps.
The kinetic energy of each particle is taken equal to \, again for physical
consistency. The option that lattice constant, time step and particle mass be
equal to unity completely defines a natural system of units for microscopic
mechanical quantities.
The propagation phase in the FHP-1 model proceeds in exactly the same way
as for the HPP model. An essential difference appears in the collision phase:
while the collisions are deterministic for the HPP model, some stochasticity
enters the FHP microdynamics. Indeed, two particles coming from opposite
directions (for example 1 and 4) undergo a binary collision with an output state
rotated by +60° or —60° (as shown in Figure 3.4(a)), with probabilities p and
1 — p respectively. Another new feature is the introduction of deterministic triple
collisions: three particles coming from three directions forming 120° angles
between each other, will be deflected by 60° (as shown in Figure 3.4(b)). All
other states remain unchanged. Note that, among the 26 = 64 possible states,
the five states shown in Figure 3.4 are the only ones that can undergo effective
collisions. The collisional efficiency of FHP-1 is 7.81% (see Section 3.1 for
definition).
fori=l,...,6:
. + 1) = nt
(3.6)
Here, ni stands for (1 — n,) and the index i is cyclic (modulo 6). All variables
on the right hand side are taken implicitly at time U and node r*. £(+)(r*, £*)
--*—X--
\ / \ / \ /
Figure 3.4 Collision rules for the FHP-1 model, reduced by symmetry. The
full set of collisions can be reconstructed from the above reduced set by the
application of all the isometries of (§. The numbers above the open arrows are
the transition probabilities. The most commonly used value p = 1/2 has been
chosen here, since it is the only value compatible with the G-invariance.
42 3 Microdynamics: various examples
/1 \ / 1 \ / 0 \
1 2 4
~2 ] 2T
m: !. ,, p x ':
Px , py : . . (3.7)
-1 0
\J
1 2 -4
\-4/
The geometrical invariants of FHP models are discussed by d'Humieres,
Qian and Lallemand (1989, 1990).
(v) The FHP-1 model is G-invariant only when p is equal to 1 — p, that is, if
p = 1/2. Indeed, consider the three states s<°> = (100100), s^ = (010010)
and s ( - ) = (001001). The mirror symmetry gx with respect to the x-axis
induces on the single-node phase-space y a, transformation (SX that leaves
s(0) unchanged and permutes s (+) and s^~\ Since &x changes the collision
s(0) _• s(+) j n t 0 th e collision s ^ —• s^ and reciprocally, the probabilities of
these two collisions must be equal in order that the model be G-invariant.
3.2 TheFHP-1 model 43
In addition, all six velocity vectors clearly belong to the same class, since
any of these six vectors can be changed into any other by a rotation of
some integer multiple of TC/3.
(vi) The FHP-1 model has the fifth order crystallographic isotropy (see Sec-
tion 2.3.5). As a justification, Table 3.2 gives the crystallographic tensors
of order 1 to 6, contracted with a test vector x of coordinates (xi,*2)-
Clearly, the contraction of the fourth order tensor is proportional to ||x|| 4 ,
and the fifth order tensor is zero. Hence, by the criterion of Section 2.3.5,
crystallographic isotropy holds up to the fifth order for the FHP lattice.
(vii) The FHP-1 model is irreducible (see Section 2.3.6). Indeed, the only isom-
etry in G that preserves ci, for example, is the mirror symmetry with
respect to the x-direction, whose invariant sub-space only contains vectors
proportional to ci. The same argument holds for C2 to C6.
As for the HPP model, the energy observable is proportional to the mass
observable and its conservation is a consequence of the mass conservation. So,
the FHP-1 model has large-scale behavior ignoring thermal effects.
The basic reason for incorporating triple collisions in the FHP-1 model follows
from the fact that without triple collisions the model would have an additional
collisional invariant, namely ( 1 , - 1 , 1 , - 1 , 1 , - 1 ) , which is physically irrelevant.
This spurious invariant would strongly affect the macroscopic dynamics of the
model and would make it unphysical.
The degree of crystallographic isotropy of FHP-1 is larger than four, which is
compatible with the possibility of having realistic large-scale dynamics governed
in certain limits by the Navier-Stokes equation, provided the model be G-
invariant, that is, for p = 1/2, which is the most commonly used value. Other
values produce chirality.
The FHP-2 model is a variant of the FHP-1 model that includes the possibility
of one rest particle per node, in addition to the six moving particles of FHP-1.
Each node then has b = 7 channels, corresponding to particles moving along
the six directions of the triangular lattice and to the rest particle. The channels
corresponding to moving particles are labeled by integers from 1 to 6, and the
channel corresponding to the rest particle is labeled 0.
The six connection vectors are the same as for FHP-1. They do not include
the vector Co which is zero. As a consequence, the connectivity B, which is still
equal to the coordination number of the triangular lattice (that is, 6), is now
different from the number b = 1 of channels per node. The six vectors labeled
by i = 1 to 6 have equal unit moduli, and the vector Co corresponding to the
rest particle obviously has zero modulus.
The masses m* of all particles are equal to one and the momentum p, of each
particle is equal to c,. Thus, the momentum of the rest particle is zero. The
kinetic energy is taken equal to \ for moving particles, and to 0 for the rest
particle. As for the HPP and FHP-1 models, the above choices are physically
consistent with unit particle mass and unit time step.
The propagation phase for moving particles is the same as for FHP-1, and
does not affect the rest particles. The collision rules of the FHP-2 model are
similar to the collision rules of FHP-1 with four additional collisions coupling
moving and rest particles: (i) a moving particle arriving at a node occupied
only by a rest particle produces a pair of moving particles with angles +60°
and —60° measured from the direction of the incoming particle (Figure 3.5.e),
(ii) the reverse of the former (Figure 3.5.f), (iii) the collisions of FHP-1 with
a passive rest particle left unchanged by the collision (Figure 3.5.C and 3.5.d).
All possible collisions of FHP-2 that cannot be deduced from one another
by a lattice-preserving isometry are summarized in Figure 3.5. Only 22 among
the 27 = 128 possible states can undergo effective collisions. The collisional
efficiency of FHP-2 is 17.2% (see Section 3.1 for definition).
FHP-1, except that there is an additional equation for no, and a few more terms
on the right hand side. The result reads:
for i = 1, . . . , 6 :
n,-(r*+c,-,f* + l) = fit
+
+
/\ 7\
\ / \/
- -K - - > - -
\/ \ / \ \ / x / \
A** A /\
Figure 3.5 Collision rules for the FHP-2 model, reduced by symmetry. The
full set of collisions can be reconstructed from the above reduced set by the
application of all the isometries of ^ . Black dots represent rest particles.
46 3 Microdynamics: various examples
for z = 0:
n0
Here, the index i is cyclic (modulo 6), in such a way that index i = 7 is equivalent
to index 1 (and not 0); other notations are as for FHP-1 (see Section 3.2.1).
Note that in (3.8) five terms are independent of no, which reflects that the
binary and triple collisions of FHP-1 occur indifferently with or without a
passive rest particle. The strict application of (2.13) would have required two
pairs of random Bernoulli variables (6(H, 6^1) and (6^1, £p\), for binary head-
on collisions with and without passive rest particle respectively. In fact, since
these two pairs of random variables have the same statistics, they are taken
equal: this leads to the above form for the microdynamical equation of FHP-2.
1\ / 0 \
1 1 0
l
1 2
1
m 1 , Px : 2 (3.10)
1 -1 0
1
1 2
1/ V5/
(v) The model is G-invariant only if p = 1/2.
The velocity vectors can belong to two different classes: The first class,
conventionally labeled by / = 0, contains only Co, the null velocity vector of
the rest particle; the second class, labeled by / = 1, contains the non-zero
velocity vectors ci to C6-
(vi) The model has fifth order crystallographic isotropy, for the same reason
as FHP-1, since the rest particles, which have zero velocity, do not modify
the crystallographic tensors of any order.
3.3 The FHP-2 model 47
(vii) The model is irreducible. Indeed, the only isometry in G that preserves
ci, for example, is the mirror symmetry with respect to the x-direction,
whose invariant subspace only contains vectors proportional to ci. The
same argument holds for C2 to C6. For Co which belongs to Class 7 = 0
however, the reasoning is as simple, albeit different: all the isometries in G
leave Co unchanged. The only vector left unchanged by all the isometries
in the point group G of the triangular lattice is of course the null vector,
which is Co.
A further variant of the 7-bit FHP-2 model is FHP-3 where the collision rules
are designed to include as many collisions as possible, under the constraint
of having the same collisional invariants as with FHP-2. Figure 3.6 shows the
collision rules of FHP-3. As for Figures 3.2, 3.4 and 3.5, the lattice-preserving
isometries are used to reduce the number of collisions shown. In addition to the
isometries, the self-duality property, that holds for FHP-3 but not for FHP-1
and FHP-2, has also been used to further reduce the number of configurations
in the figure. A complete collision table could be reconstructed from Figure 3.6,
by applying duality combined with the isometries of G. When this reconstruction
is performed, it can be shown that 76 among the 27 = 128 possible states can
undergo effective collisions. The collisional efficiency of FHP-3 is 59.4 % (see
Section 3.1 for definition).
\/ \/
/ \ / \ / \ x / \ /
_\/_A/_ CrN -:4r-*--
A A
x / \ / \ / \ /
\ /\ /\ / 1 \ / \^\ /
~/ jif^^\ I ^ ~p TV \
/ \ / \ / \
V;
rX—V-
/\ / \
X / \ /
• * "
-iL—^~-
\ i/4
1/4 -V- ~^~1r~^
7\~~
x / \ /
^ -^
-xrV- ~ A ~ A""
\|2 -;\ / \ /
-?-*—
•£.*7\~- 7^ \ / \ /\ /
1/2 _^_J*__\^
Figure 3.6 Collision rules for the FHP-3 model, reduced by symmetry. The
full set of collisions can be reconstructed from the above reduced set by the
application of all the isometries of ^, combined, when necessary, with the
duality operation (exchanges between particles and holes). Black dots
represent rest particles.
50 3 Microdynamics: various examples
(iv) The only collisional invariants are the mass and the two components of
the momentum, as for FHP-2 and FHP-1.
(v) The model is G-invariant if, and only if, p = 1/2.
The class structure is the same as for FHP-2.
(vi) The model has fifth order crystallographic isotropy.
(vii) The model is irreducible, for the same reason as FHP-2.
The models presented so far involve identical particles with identical masses, that
cannot be differentiated by any physical property but their velocity. Physically,
these models involve a single 'chemical species', and all the resulting dynamics,
at any scale, is purely mechanical. In order to study diffusion phenomena,
the FHP-3 model must be modified to include two different chemical species
(Bernardin and Sero-Guillaume, 1990; McNamara, 1990; Noullez, 1990; Hanon
and Boon, 1997). The version presented here is the simplest one: the two
species are chemically passive and mechanically identical (same mass); there
is no chemical reaction between species (the collisions conserve the number
of particles of each species), and the mechanical behavior of the particles is
not affected by the specificity of the species. Consequently, the macroscopic
motions of the mixture are identical to those of a single component fluid such
as FHP-3, and each species acts as a passively diffusing scalar with respect to
the mixture. To get a better picture, we can associate colors (say red and blue)
to the two species, and consider the CFHP ('colored' FHP) model as an FHP-3
model whose particles are painted red or blue, without modifying their motions
and collisions. Indeed, the FHP-3 collision rules apply to the particles arriving
at a node, independently of their color. Then, the color attribute is attached
3.5 The 'colored' FHP model (CFHP) 51
randomly to the outgoing particles, with the constraint to conserve the number
of red and blue particles.
There exist several variations on the theme of colored lattice gases, e.g.
incorporation of reactive processes (see Boon et a/., 1996) or surface tension
effects (see Rothmann and Zaleski, 1994). They will not be considered here (see
Chapter 11).
The encoding of the CFHP model raises a little problem. Indeed, each channel
can be a priori in three different states: empty, occupied by a red particle, or
occupied by a blue particle. This situation does not seem compatible with a
binary encoding. One way to recover the Boolean formalism is to consider that
each of the seven channels of the FHP-3 model is now a double channel: one
'sub-channel' for the red particles, and one 'sub-channel' for the blue ones.
Fourteen Boolean quantities per lattice node are now required, as illustrated in
Figure 3.7. The (sub-)channels 0 to 6 encode the presence of red particles, and
the (sub-)channels 7 to 13, the presence of blue particles. The microdynamical
equation of the CFHP model is designed to avoid the 'forbidden' situations
where both red and blue sub-channels are occupied simultaneously. Of course,
initially, the lattice configuration has to be preset so that each pair of sub-
channels contains at most one particle (either red or blue). This rule expresses
an exclusion principle per pair of sub-channels.
y y y y y y
/ \Y
V
A
A AyA& \7 A A A
Y
A V
A
Vc V
'4
;\/ \ /
/ \ A
y
A
V
4
A
\/ A A \ / \/ \
A A A A A A
Figure 3.7 The 14 channels and velocity vectors c, for the colored-FHP model
(CFHP). 'Red' particles reside on channels i = 0 to i = 6, and 'blue' particles
on channels i = 7 to i = 13. The velocity vectors of channels 0 (red) and 7
(blue) are null and are not labeled. They are represented by the two concentric
circles.
52 3 Microdynamics: various examples
This evolution rule clearly guarantees that if the initial state does not have more
that one particle per pair of subchannels (red or blue, but not both), then so
will the final state. From now on, we will assume this is the case. Thus, the
single-node phase space y is restricted to the 3 7 'authorized' states, among the
2 14 a priori possible states of a model with b = 14.
(vi) The model has fifth order crystallographic isotropy, as FHP-3 and for the
same reasons,
(vii) The model is irreducible, for the same reason as FHP-3.
The two-dimensional HPP and FHP models do not contain thermal effects,
because there is no energy conservation independently of mass conservation
(the energy observable is proportional to the mass observable). So these models
are restricted to the description of non-thermal fluids.
We now introduce a lattice gas model with non-trivial thermodynamics:
the '19-bit model', constructed by Grosfils, Boon and Lallemand (1992); for
consistency in acronymic notation, we shall refer to this model as the 'GBL
model'. The main feature of the model is that it is a multi-speed lattice gas
residing on the two-dimensional triangular lattice. Figure 3.8 shows how the
velocity vectors connect any lattice node to (i) itself (one rest particle per node;
||Co|| = 0), (ii) its six nearest neighbors (||CJ|| = 1, i = 1,...,6), (iii) its six next-
nearest neighbors (||Cj|| = ^/3, i = 13,...,18), (iv) its six next-to-next-nearest
neighbors (||c;|| = 2, i = 7,...,12).
The number of channels b is equal to 19, whereas the connection number B
is 18, since there is one rest particle per node. Both numbers are independent
of the coordination number of the triangular lattice (that is, 6).
A7
Ci Y ^ ^7 y y v
Y Y
A
AA/ A/ v \/c'\/ N\ A
y V V v\ A
C / Y
A A V A/ A
/V\i
Y \ / \
c x
/ [7v\/ y
\ A A 11
Cl8
V« A A
The CjS have moduli 0, 1, ^/3 and 2 and their Cartesian components are:
co =
0 '
^ 'l\ ^_(~\
c4 = I . , c5 = % ) , c6 =
C 7 = l A ) , C8=( R \ , C9 = ! / - ! - (3-H)
1
ClO = I I, Cn = I _ ), C12 =
n R
-V5
C13 = | t/i I , CM = I j | ], C1S =
C C
16 = I J I' 17 = I ^ I. C« =
\ 2 / \ 2 /
The masses m, of the particles are equal to unity; their momenta are equal
to the CjS, and their kinetic energy is equal to ^ ||c/1|2. It is clear that the kinetic
energy observable (0, 5,..., \,2,...,2, §,..., \) is not proportional to the mass
observable (1,...,1). In contrast to FHP-2 and FHP-3, rest particles carry no
internal energy.5 Collisions can couple the different velocity modulus levels, and
conserve the kinetic energy observable. In other words, there exist collisions
which conserve the mass, momentum and kinetic energy observables, without
conserving individually the population of each velocity modulus level. Figure 3.9
shows some of these 'population-mixing' collisions, which do not exist in the
FHP models. The GBL model also includes all FHP collisions at the different
velocity modulus levels. The complete set of collision rules involves too many
possible cases to be shown explicitly. However, the general rule can be defined
in a synthetic way as follows. We define a family of Boolean states as a subset
of the single-node phase space y that contains all Boolean states with the same
total mass, the same total momentum, and the same total kinetic energy. The
set of all these families defines a partition of the single-node phase space, since
any possible Boolean state belongs to one and only one family. The collision of
the GBL model can be defined by a collision matrix A(s^s') such that:
An exhaustive count shows that 517 750 among the 219 = 524288 possible states
can undergo effective collisions. The collisional efficiency of GBL is 94.4%
(which is different from 517 750/524288 because ineffective collisions can occur
with probabilities different from 0 or 1). Note that an increase in collisional effi-
ciency (raising its value up to 98.9%) can be gained by imposing that ineffective
collisions are only admitted for Boolean states belonging to families with only
one element. The probability for an effective collision 5 —> s' is then A(s^s') =
(card(?F) — \)~ , if both s and sf belong to a family with more than one element.
x
V - \/_ v
^\ A A \ A 7
\/ \ / \ / \ / \/ \ /
A A
/ 0 \ / 0 \ /0\
1 Cly l
2
1 Cly 2
m : . ,, px ::
Px *:. », p« '•:
Py . ,, ee :: . . (3.12)
3
1 2
(v) The model is G-invariant. The proof is quite similar to the proof of self-
duality. Any isometry g of G induces a relation between families with the
same number of elements. The G-invariance follows immediately.
The class structure within the c, vectors is more complex than for FHP
models. The class 7 = 0 only contains Co, the (null) velocity vector of the
rest particle. The class 1 = 1 contains the velocity vectors ci to C6, which
have unit modulus. The class 1=2 contains vectors Cj to C12, with modulus
equal to 2. The class 7 = 3 contains vectors C13 to cig, with modulus y/3.
(vi) The model exhibits fifth order crystallographic isotropy, since its set of
velocity vectors can be viewed as the superposition of a zero vector plus
three groups of six vectors, each group being obtained from the six vectors
of the FHP model by a similitude operation (rotation and rescaling). The
degree of crystallographic isotropy is thus the same as for FHP models.
(vii) The GBL is irreducible. Indeed, the only isometry in G that preserves ci
is the mirror symmetry with respect to the x-direction, whose invariant
subspace only contains vectors proportional to ci. The same argument
holds for vectors C2 to cig. For Co, the same argument as for FHP-2 holds.
the two monoclinic lattices (simple and body-centered), the triclinic lattice and
the trigonal lattice, do not even satisfy the second order crystallographic isotropy.
The three-dimensional hexagonal Bravais lattice can serve as a basis for minimal
models which satisfy the second order crystallographic isotropy only if its aspect
ratio is equal to \/3/2, but even in this case, the fourth order condition is not
satisfied.
So, strategies must be developed to overcome this difficulty, and three-
dimensional lattice gas models can indeed be constructed with fully isotropic
large-scale dynamics. These strategies are discussed in the next three sections.
6
This is not in conflict with the exclusion principle, since particles with equal velocities
are assumed to be distinguishable.
3.7 Three-dimensional models 59
and
Clearly, the right hand side of (3.13) is proportional to ||x||2 for any choice
of Ci and £2, but the right hand side of (3.14) is proportional to ||x||4 only if
^Ci + 2^2 = §Ci, that is, for d = 2^2- The easiest way to satisfy this condition
is to choose £1 = 2 and £2 = 1, which leads to b = 8£i + 6^2 = 22 channels per
node.
The same idea can be applied to the simple cubic lattice, with propagation to
nearest and next-nearest neighbors with suitable weights. It yields a model with
B = 18 connections (6 to nearest neighbors and 12 to next-nearest ones) and
b = 24 channels per node (each of the six links to nearest neighbors is a double
channel).
(i) The four symmetries Si, S2, S3, and S4 with respect to the four hyper-planes
x\ = 0, X2 = 0, X3 = 0, and X4 = 0 (coordinate reversals),
(ii) The six symmetries Pn, P13, Pu, P23, P24, and P34 with respect to the
hyper-planes x\—X2 = 0, x\ — X3 = 0, xi — X4 = 0, X2 — X3 = 0, xi — X4 = 0,
and X3 — X4 = 0 (coordinate exchanges),
(iii) The symmetries Zi and E2 with respect to the hyper-planes X1+X2—X3—X4 =
0, and xi — X2 — X3 — X4 = 0 respectively.
This set of symmetries generates the point group G of the FCHC lattice, but it
is not a minimal set since only five of them, for example (Si,Pi2,Pi3,Pi4,£i) are
sufficient to generate G. A complete discussion on the symmetry group of the
FCHC lattice is given by Henon (1987a,b).
The minimal set of velocity vectors that can be built on the four generating
vectors is:
±1\ (-
±1 0
C5,...,8 =
0
o/ Vo (3.15)
o\
±1
f°
1
0
±1 » Cl7,...,20 =
o C21,.,24= + 1
0 / \+l
Here, the numbering of the C;S is arbitrarily based on the lexicographic order
of the alphabet (+1,-1,0). Note that for operational convenience in numerical
simulations, another order may be chosen in such a way that c;+i2 = — c* for
i= 1,..., 12. The reason for this choice is related to technical details concerning
the practical implementation of bounce-back reflections on solid obstacles (see
Chapter 10).
A three-dimensional lattice gas model can be constructed from the four-
dimensional FCHC lattice by the following dimension-reducing procedure:
(i) We consider the portion of the four-dimensional lattice limited by the two
hyper-planes X4 = +1 and X4 = — 1. This (hyper-)slice of lattice corresponds
7
The acronym FCHC is slightly ambiguous; in the Appendix, Section A.2, we show that
the lattice should be named 'two-dimensional face-centered body-centered hyper-cubic'
lattice. For simplicity and convenience, we nevertheless use FCHC as commonly
accepted.
3.7 Three-dimensional models 61
to one lattice period in the fourth direction, and tiles exactly the FCHC
lattice by successive translations of vector (0,0,0,2).
(ii) We impose periodic conditions in the fourth direction; in other words, we
cancel the fourth component of the c,-s and project all the nodes in this
portion onto the hyper-plane X4 = 0.
Figure 3.10 illustrates this procedure for a simple case: the construction of a
one-dimensional lattice from a two-dimensional square lattice.
Figure 3.11 shows how the three-dimensional (pseudo-four-dimensional) model
so constructed is merely the three-dimensional cubic lattice, with double chan-
nels connecting nearest neighbors and single channels connecting next-nearest
neighbors. There are still b = 24 channels per node in three dimensions, but six
of them are double channels (B = 18).
The three-dimensional model constructed by projection of a four-dimensional
lattice can also be viewed as a multiple-link model based on the simple cubic
lattice with propagation to nearest and next-nearest neighbors with multiplicity
2 and 1 respectively. While it is easier to describe this model as a multiple-link
lattice gas model (Rivet, 1987a), but because of its lack of numerical efficiency,
the isometric algorithm was not further exploited.
form would involve more lines than the present book contains; the only useful
expression is the general form given by (2.12) or (2.13).
(i) The class structure may depend on the number of rest particles. However,
the 24 velocity vectors of the moving particles always belong to the same
class. Indeed, let us take ci (see Equation (3.15) for the definition of the C;
vectors). It is clear that C2 to C4 belong to the same class as ci, since the
application of combinations of the symmetries Si and S2 to ci generates
these three vectors. The application of P23 to vectors ci to C4 generates
vectors C5 to eg, and so on. Thus, any vector ci to C24 can be changed into
any other vector c* by an isometry in G, which proves that all 24 vectors
belong to the same class.
(ii) The mass and the four components of the momentum are the only linearly
independent collisional invariants. The energy observable is proportional
to the mass observable and so is trivially conserved for models without
The reader may anticipate the following question: in order that three-
dimensional lattice gas large-scale dynamics be governed by the incompressible
fluid Navier-Stokes equations, the collisional invariants must be the mass, the
energy, the three components of the momentum, and no other quantity ; but what
about the fourth component of the momentum which is also conserved in FCHC
models ? In fact, this additional invariant is not spurious, as it corresponds to an
additional macroscopic variable with its own large-scale dynamics, but which
does not perturb the large-scale dynamics of the three-dimensional relevant
macroscopic quantities. Indeed, consider a fully four-dimensional FCHC model
with large-scale collective dynamics governed by the four-dimensional Navier-
Stokes equations. Now, when cutting a (hyper-)slice in the four-dimensional
lattice and wrapping it onto itself in the fourth dimension (see Section 3.7.3),
the macroscopic variables no longer depend on the fourth component X4 of
5 0 isotropic
6 12(5
l^ll 6) anisotropic
3.7 Three-dimensional models 67
68
4.1 The Liouville description 69
4.1.1 Macrostates
Most common laboratory probes (thermometers, pressure sensors, velocimeters,
etc.) yield time- and space-averaged measures, simply because they have a finite
macroscopic size and a finite response time, that is, they operate on a length
scale and a time scale which are large compared to the characteristic microscopic
lengths and times. For theoretical developments however, it is logical and
convenient, following Gibbs, to consider ensemble-averages (Huang, 1963). This
procedure requires the ergodic hypothesis assuming that ensemble-averages
correctly describe macroscopic observables which are in practice time- and
space-averages.
Formally, we define a 'macrostate' of the lattice as a function acting on the
phase space F, such that it associates to each microstate S e F a real number
between 0 and 1:
*: F - [0,1]
V }
1. (4.2)
SGT
4. l .2 Ensemble - averages
Having defined a probability measure in the phase space, we can now define
ensemble-averaged quantities. Consider an arbitrary generalized observable (see
Section 2.1.4), that is, a physical quantity whose measure at node r* at time
U is an arbitrary function q(n,u) of the Boolean field n, and possibly of the
70 4 Equilibrium statistical mechanics
*)P(S). (4.3)
SeT
SeT i=l
which reduces after summation inversion to:
b
, (4.4)
Although formally different from the Liouville equation for continuous systems,
Equation (4.6) plays an equivalent role for lattice gases: it governs the time-
evolution of the probability measure in phase space T. We therefore call (4.6)
the lattice Liouville equation'.
We now address the problem of finding the equation governing the time-
evolution of the ensemble-averaged populations of each channel at each node.
propagation phase shuffles the particles in such a way that their correlations
vanish at each node, so that the pre-collision (post-propagation) state of each
node be totally factorized. Obviously, right after a collision, particles are in a
strongly correlated (post-collision) state precisely because of the effect of the
collisions, but we assume that during the subsequent propagation phase, these
correlations are completely damped out.
The Boltzmann approximation, for lattice gases as for real gases, is a strong
hypothesis whose validity is expected to hold at sufficiently low gas density.
For lattice gases however, the approximation is rather well verified even at
moderate densities, at least for models satisfying semi-detailed balance (like the
GBL model and the FCHC-3 model; see Chapter 7).
/,(r
s,s'ey
where the rij s on the r.h.s. are implicitly taken at time U and node r*. For clarity,
these arguments are not written explicitly, and will be omitted systematically
when no ambiguity results. We also use the more compact notation ft(r*,t*)
rather than /(c,-,r*, £*), which is formally closer to the classical notation /(p,r,t),
keeping in mind that /;(r*,£*) and /(c;,r*,f*) denote the same object: (w,-(r*,**)).
To arrive at (4.7), we have used the fact that the £(s-»s')s in (2.12) are random
Bernoulli variables, statistically independent from the Boolean field, and whose
averages are the ^4(s-»s')s.
Equation (4.7) is an exact equation but, as it stands, there is no way to express
explicitly the r.h.s. in terms of the fts. Therefore, we now use the Boltzmann
approximation (see Section 4.2.1), that is, we make the assumption that the
pre-collision configuration is totally uncorrelated. With this assumption, the
ensemble-average of the product on the r.h.s. of (4.7) can be factorized into the
product of the ensemble-averages, to obtain:
fj', (4.10)
s,s'ey 7=1
/ =klog2n-H,
where k is an arbitrary positive constant that will be chosen equal to one if H is
to be expressed in bits. We shall not reproduce Shannon's proof, but rather use
a simple argument to convince ourselves that this measure is 'reasonable' and
76 4 Equilibrium statistical mechanics
consistent with the intuitive and natural results of the deterministic case. Indeed,
suppose that an incomplete measure reveals that x can take indifferently any
of the values x^\..., x(m) (m < n), and no other value. This partial information
is equivalent to the statistical information that the probabilities p ( 1 ) ,...,p ( m ) are
equal to 1/m, and p(m+1\...,p^ are zero. With this set of probabilities, (4.11)
yields H = log 2 m, which is consistent with the result presented above for the
deterministic case.
YeT
I =log22q + ^0>(Y)log20>(Y) (4.12)
Yer
= q-H.
Suppose now that the only available knowledge is not the full statistics of
X, but the individual probabilities pi, i = l,...,q that x\ take the value 1.
According to (4.11), the entropy and the quantity of information associated
with the knowledge of the individual variable xio are, respectively:
H
k = ~ I Pio lo g2 Pio + (! - Pio) loS2 (* ~ Pio) 1»
L J
(4.13)
since x;0 is a Boolean variable which can take two values: 1 or 0, with respective
probabilities pt and (1 —pi). The entropy and the quantity of information associ-
ated with the knowledge of all the individual probabilities are thus, respectively:
q (4.14)
I* =^2l+ [Pi loS2 Pi + Pi l°g2 P]
where pt stands for 1 — p,-. Since the knowledge of the individual statistics does
not account for possible correlations between the variables, it is logical to expect
4.3 The if-theorem 77
that it contains less information than the knowledge of the full statistics of X.
The proof for this quite intuitive statement is given in the Appendix, Section A.3.
We now apply these concepts to lattice gases.
sey
Notice that we now use natural logarithms rather than logarithms to the base 2.
(The difference is a non-significant multiplicative constant that we shall ignore.)
To simplify the notation, we shall omit the arguments r*, U in p(r*, £*, s).
The effect of the collision phase on the local probability distribution function
p is a linear transformation involving the collisional transition probabilities
v4(s->s') (see Chapter 2), that is, the post-collision local probability distribution
function pf can be expressed as:
We now use a basic property of convex functions: for any convex function 0,
the following inequality holds:
¥ , 6 r ,
where the sums implicitly extend over all 5 in 7. If we assume that the semi-
detailed balance condition ^2sey A(s^s') = 1 is satisfied, then, after summation
over all s', (4.17) reads:
( )
s'ey \sey / s'ey sey
Using (4.16) and the normalization condition J2S'ey ^(s-^O = 1> w e obtain:
s'ey sey
78 4 Equilibrium statistical mechanics
Equation (4.19), applied to the convex function (j)(x) = —xlnx, implies that:
(4.20)
s'ey sey
This result proves that under the assumption of semi-detailed balance the
collision phase increases the local entropy of the lattice gas, and therefore
decreases the quantity of statistical information we have about the state of the
lattice node.
Note that without the semi-detailed balance assumption the collision phase
can decrease the local entropy of the node. As a simple example, consider
a lattice gas model with a collision rule such that the post-collision state of
any collision is always the same state so, whatever the pre-collision state. This
unphysical model does not satisfy the semi-detailed balance condition. Suppose
that, initially, this lattice gas is in a macrostate such that all configurations are
equally probable, namely: p(s) = 2~b, for all s in y. Then, the local entropy
initially has its maximum value b x In 2. After the collision phase, the lattice is
in a completely determined configuration, where all the nodes are in state so-
Thus, the local entropy after the collision phase has decreased to its minimum
value: zero. We now make the Boltzmann approximation, i.e. we neglect the
correlations between particles entering a collision. Then, the local entropy at the
node before collision is equal to the sum of the entropies in each channel:
b
- £ [/,In/,, +/,ln/J = -5>(s)lnp(s). (4.21)
1=1 sey
Now, there is no reason that the post-collision state be uncorrelated; so, from
the results of Section 4.3.1.3, we have, after the collision:
b
Combination of (4.20), (4.21) and (4.22) leads to the local form of the if-theorem:
b b
cannot decrease under the effect of the collision phase. This quantity is called the
local entropy.
4.3 The if-theorem 79
We now introduce the 'global entropy' H which is the sum of the local
entropies h(u), extended over all lattice nodes r*. If the lattice is either infinite
or finite with periodic boundary conditions, the propagation phase only moves
information from node to node without loss or gain. It thus has no effect on
the global entropy. As a consequence, the local if-theorem implies that the
global entropy cannot decrease under the action of the full evolution operator
(collision and propagation).
Global if-theorem: Under the Boltzmann approximation and with the semi-detailed
balance condition, the global quantity
, **) +fi(r*9t*)]nfi(u9t*)]
x+eS£ i=\
cannot decrease under the action of the lattice gas evolution rule, at least if the
lattice is either infinite or finite with periodic boundary conditions. This quantity
is called the global entropy.
This theorem, proved under the Boltzmann approximation, implies that the
evolution of the statistics of a lattice gas is irreversible.
Let us illustrate the if-theorem with a simple example. We consider a lattice
gas residing on a finite lattice with periodic boundary conditions and prepared in
a homogeneous macrostate. In the homogeneous macrostate, the global entropy
H(u) is merely Jf times the local entropy /i(r*,£*) which no longer depends on
the position variable r*. So the statement that H cannot decrease when semi-
detailed balance is verified, at least within the Boltzmann approximation, also
applies to h. Suppose now that we prepare the lattice gas in a non-equilibrium
macrostate, for instance, by biasing in a systematic way one of the collision
rules; then, at some later time, we restore the usual rules (satisfying detailed
or semi-detailed balance), and in the course of time, we perform successive
measurements which give us an estimate of the entropy. We observe, as shown
in Figure 4.1, an increase in the numerical estimate of h(U) (see Tribel and Boon,
1997) towards an asymptotic plateau corresponding to the equilibrium value as
predicted by Boltzmann's theory. Now, in Figure 4.1, we also observe that on
top of the global increase of h(u) (in accordance with Boltzmann's theory) there
are fluctuations where h can decrease locally; the reason is that the Boltzmann
decorrelation hypothesis is only an approximation.
Along the same lines of reasoning, but without the Boltzmann approximation,
one can also prove a weaker form of the H-theorem:
Liouville if-theorem: With the semi-detailed balance condition, the quantity
80 4 Equilibrium statistical mechanics
cannot decrease under the action of the lattice gas evolution rule, at least if the
lattice is either infinite or finite with periodic boundary conditions.
Now that we have established the if-theorem, the next step in the statisti-
cal mechanical description of the lattice gas is to investigate the existence of
equilibrium situations, that is, of macrostates for which all statistically relevant
quantities are stationary.
Hardy, de Pazzis and Pomeau (1973) showed that the HPP lattice gas model
has very simple statistical equilibrium solutions that they call 'invariant states'.
These states have homogeneous (node-independent) properties, and they show
no equal time (static) correlations between channels of any pair, whether on
the same node or not. Such global (homogeneous) equilibrium macrostates are
analogous to the Maxwellian equilibrium states of real gases. We shall prove
that the equilibria, for a quite general class of lattice gases (including the
HPP model), are not characterized by the Maxwell-Boltzmann distribution of
1.960
1.955 -
1.950 -
1.945 -
1.940 -
1.935 -
1.930
20 40 60 80 100
U
Figure 4.1 The local entropy h{U) = - ^ = 1 [/*(£*) ln/i(r*) + /,(£*) In/,(£*)] as
a function of time, for an FHP-1 lattice gas model on a lattice with 256 x 256
nodes, u is in the natural microscopic time unit: the automaton time step.
(From Tribel and Boon, 1997.)
4.4 Global equilibrium macrostates 81
where the /,-s should be taken for the moment as b real numbers between 0 and
1, and where overlined quantities ft and S,-(r*) denote the complements: 1 —/,,
and 1 — S/(r*), respectively. Note that the distribution (4.24) can be interpreted
as the probability distribution of a collection of b x Jf independent random
Bernoulli variables with mean values /,-, i = 1,..., b which depend on the channel
index i, but not on the node position r*. The expression (4.24) is obviously the
most general form for a fully factorized and spatially homogeneous equilibrium
measure.
Notice that the identity:
b
ft = J2 S^ I I I I fjS^f)fjSs^\ Vi = 1,..., b, Vr* e if, (4.25)
ser x*'ese ;=i
which follows straightforwardly from simple algebra, means that the probabili-
ties /; in (4.24) are the ensemble-averaged populations (n,-) of each channel.
82 4 Equilibrium statistical mechanics
To proceed, we must assume that the lattice if is either infinite or finite with
periodic boundary conditions (see Chapter 2). Under this assumption, all the
nodes of the lattice are equivalent, and the free streaming operator only transfers
information from node to node without any loss or gain. Since the equilibrium
solutions are supposed to be spatially uniform (the /,-s do not depend on the
node position r*), it becomes clear that 9(3f(S')) = &(S'), for any microstate
Sf in F. If we now look for fully factorized solutions of the form given by (4.24),
we must solve the following set of equations:
i=l sey
f j=i
The proof of equivalence does not involve any restricting hypothesis on the
lattice gas model. Furthermore, the equivalence is rather intuitive and natural
since (4.28) can be viewed as a local (single-node) version of (4.27). The physical
content of (4.28) is that the single-node equilibrium probability distribution
Ilf=i fiSift is kft unchanged by the application of the collision operator. Indeed,
the l.h.s. of (4.28) is the probability of finding the node in state sf before collision,
and the r.h.s. can be understood as the probability of finding the node in state
s' after collision.
At this point, it is not obvious whether the set (4.28) has a solution, since
it imposes 2b conditions to only b variables (the /,-s). But we shall show that
(4.28) reduces to at most (and in fact less than) b independent relations. This
'little miracle' requires the assumption of the semi-detailed balance property
(Chapter 2). With this assumption, (4.28) is strictly equivalent to the following
set of b equations:
Proving the equivalence between (4.28), (4.29) and (4.30) is not a trivial matter;
the proof involves a trick introduced in the context of gas models with discrete
velocities (Gatignol, 1975), and is given in the Appendix, Section A.5.
Before proceeding any further, it is useful to discuss the physical content of
(4.29) and (4.30). Equation (4.29) expresses that for a global equilibrium, the
ensemble-averaged population of any channel is left unchanged by collisions.
This interpretation becomes clear when we use the identity (see proof in the
Appendix, Section A.6):
S
ft = £ ^ M ) II //'/A Vi = 1,..., i, (4.31)
s,s'ey j=l
which yields:
/*=/r =
(4.34)
\K=1 /
1
Here v>FD(x) = (1+e*)" is the Fermi-Dirac distribution. The fact that we obtain
a Fermi-Dirac distribution instead of the Maxwell-Boltzmann distribution of
continuous gases, is a consequence of the exclusion principle. Equation (4.34)
describes a family of homogeneous and fully factorized equilibrium solutions
of the lattice Liouville equation (4.6), depending on d arbitrary parameters A M,
where S is the number of linearly independent collisional invariants of the
model.
The parameterization by the l M s comes out naturally in passing from (4.30)
to (4.33). In Section 4.4.3, these arbitrary coefficients will be interpreted as
Lagrange multipliers in the global entropy maximization procedure.
Two important features must be emphasized:
Since (4.36) is identical to (4.29), it follows that the lattice Boltzmann approach
is equivalent to the lattice Liouville approach if semi-detailed balance is satisfied;
consequently, the equilibrium solutions of the lattice Boltzmann equation are
Fermi-Dirac distributions.
1=1
where /,• does not depend on the position variable r* because of the homogeneity
of the equilibrium state.
In order that the macrostate be an equilibrium solution, the global entropy
must be maximum, under constraints imposed by conservation. Indeed, consider
a basis (q^K\ k = 1,..., d) of the vector space of collisional invariants. The quan-
tities (q^) = ^2q\ fi are constant, because the q^ s are collisional invariants.
86 4 Equilibrium statistical mechanics
(4.38)
where the equilibrium mean populations ffq) are given by the Fermi-Dirac
distribution (4.34).
We define the macroscopic variables p M as the macroscopic densities per node
of the collisional invariants q^K\ that is, the microscopic densities per node p* M ,
6
Section 4.5 contains material which is important for lattice gas practice, but it is rather
technical and may be omitted at first reading.
4.5 Natural parameterization of equilibria 87
ensemble-averaged over the probability distribution ^ (eq) (see Section 2.1.3 for
the notion of density of an observable):
This definition of the macroscopic variables leads to the following set of S non-
algebraic relations between the macroscopic variables p M and the Lagrange
multipliers A M :
If we are able to solve this set of equations, i.e. if we can write the A Ms as
functions of the p M s , then we can express the equilibrium distribution (4.34) in
terms of the macroscopic variables p M .
Unfortunately, the general solution to (4.39) cannot be cast in an explicit
form, and we have no exact expression for the equilibrium distribution in
terms of macroscopic variables. However, approximations can be obtained by
perturbation methods. Indeed, suppose that we have an exact solution ( J ^ K =
1,...,(5) for (4.39), corresponding to some particular values ( p ^ , K = 1,...,<5)
of the macroscopic variables. Then it is usually possible to find an approximate
solution (AM, K = 1,...,(5), corresponding to macroscopic variables (p M , K =
1,...,8) close to (pM, K = 1,...,<?).
We now sketch the asymptotic expansion leading to this approximate solution.
Let rj (<C 1) be an expansion parameter such that the perturbed macroscopic
variables read:
(The subscript between parentheses refers to the order in the expansion scheme,
whereas the superscript between square brackets refers to the index of the
associated collisional invariant.)
We look for solutions of (4.39) with the following form:
l (4.41)
the Fermi-Dirac distribution and of its first and second derivatives respectively,
computed for the unperturbed solution ( J ^ K = 1,...,<5):
\K=1 /
Using the explicit form 1/(1 + expx) for \pFD, it is easy to show that:
ft = xpFD I V j A ^ g } I
\K=1 /
= hi
3
^W1 (4-44)
(
K=l
^ \K=1 )
i=l,...,b.
It then follows that (4.39), to first order in rj, leads to the following linear set
forOlM, K = 1 , . . . , < 5 ) :
3
£Mfc*UM = pM, (4.45)
K' = l
This matrix is invertible (see proof in the Appendix, Section A.7); so the first
order linear set (4.45) has a unique solution (/ifj^, K = 1,..., 6), which is explicitly
given in terms of the macroscopic variables.
The second order of the expansion in powers of rj leads to a similar linear set
for the X^s, given in terms of the Afj^s:
EM M W = -\E^f»(E
4.5 Natural parameterization of equilibria 89
This set also has a unique solution. The structure is similar for higher orders,
but the r.h.s. grows more and more cumbersome.
The next step is to re-introduce the expressions of the perturbations X^ and
2^ into (4.44) to obtain the desired approximate expression to the average
populations /,-.
This expansion is quite general, but rather formal. We now turn to examples so
that these notions become more practical. We first apply the general expansion
scheme to a class of simple lattice gases which contains most models with fluid-
like non-thermal macroscopic properties (Section 4.5.1). Then we deal with the
more complex case of thermal models (Section 4.5.2).
(4.48)
where the Greek index a refers to Cartesian coordinates and ranges from 1
to D.10
Property (v) - all particles have same mass - leads to considerably lighter
algebra. The formalism, however, extends very simply to the case where only
8
Property (iv) excludes models FCHC-6, -7, and -8 (see Table 3.3).
9
Property (v) excludes models FCHC-5, -7, and -8 (see Table 3.3), because they have rest
particles with mass 1 and 2.
10
Property (vi) excludes the thermal GBL model which has an additional independent
collisional invariant (the kinetic energy, see Section 3.6.2), and will be treated in
Section 4.5.2.
90 4 Equilibrium statistical mechanics
moving particles have the same mass, while the mass of rest particles can take
different values (for instance, for models FCHC-5, -7, and -8).
We now use property (ii) (crystallographic isotropy up to order 2) to define
a characteristic number, £2, which will be useful in the subsequent calculations.
Property (ii) implies that the crystallographic tensor J2ici« cip *s isotropic (see
Section 2.3.5 for a definition of the notion of isotropic tensors); so there exists
a number £2 such that:
i= l
h (4.50)
[a]
p = jSk/,- (a = 1,...,/)).
We choose to study nearly equally distributed equilibria, that is, global equi-
libria for which the average population of all channels is close to a common
value. The basic unperturbed equilibrium is a macrostate in which all channels
have equal probability to be occupied: /<>,• = d, Vi = l,...,fe. The parameter
d (0 < d < 1) can be interpreted as the 'particle density per channel'.
These equilibria are also called 'low-speed equilibria', because the mean ve-
locity of particles at any node is small, as soon as the populations of all the
channels on that node are almost equal. For instance, if there are approximately
as many upward moving particles as downward moving particles, the mean
velocity is obviously small.
The values of the macroscopic variables corresponding to the unperturbed
(fully equally distributed) equilibrium are: 11
(4.51)
11
Remember that the subscript between parentheses refers to the order in the expansion
scheme, whereas the superscript between square brackets refers to the index of the
associated collisional invariant.
4.5 Natural parameterization of equilibria 91
For these values of the macroscopic variables, it is easy to verify that the
following Lagrange multipliers satisfy (4.39):
AM=0, a =1,...,£>.
We look for nearly equally distributed equilibrium states for which the mass
density is still p, but the momentum density is now j = 0 + rj]{1)9 where rj is an
expansion (small) parameter. The perturbations (p^\ K = 1,...,<5) are:
P,™ = 0,
(4.53)
From the unperturbed values (4.52) of the Lagrange multipliers, we can compute
the values of the expansion coefficients /o,- (using (4.42) and (4.48)), and fu and
f2i (using (4.43)):
hi = d,
fu = -d(l-d), (4.54)
As a consequence of the fact that all particles have equal mass, these quantities
are independent of the channel index i.
It is now straightforward to compute the 5 x 5 real matrix M defined in (4.46).
We find:
/b 0 . . 0\
0 i2 . . 0
(4.55)
\0 0
Since M^'^ is diagonal, the solution of the first order problem (4.45) is
simply:
AJJ = 0,
i / (4-56)
1 a
[] ^
case of thermal models. For non-thermal models, the algebra is similar, but
simpler.
The resulting values for the second order corrections 1 ^ to the Lagrange
multipliers are:
l 2 d
; [0] = t
« W(\-Wb& (457)
The final expression for the perturbed mean populations /,- can be obtained
from (4.44), using the expressions (4.56) and (4.57) for the first and second order
perturbations 1^ and X^ to the Lagrange multipliers. The resulting expression
reads, in terms of the macroscopic variables p = bd and j : 1 2
U
~ b
Ql
(4.58)
(4.59)
(4.61)
12
Remember that Einstein's convention is used for implicit summation over repeated
Greek indices.
13
This factor and the 'non-Galilean factor' g(p) in Chapter 5 differ by a multiplicative
constant.
4.5 Natural parameterization of equilibria 93
which gives:
/(c) ~ const. 1 +
1 / m
cacp - ^&ap) + #(w3)] • (4.62)
2 VfeBT
'H\
. . (4.63)
This set of properties is almost the same as in Section 4.5.1, except that fourth
order crystallographic isotropy is now required, and that energy conservation
is taken into account. The GBL model for thermo-hydrodynamics is a typical
example of this class of LGAs.
We use property (ii) (crystallographic isotropy up to order 4) to define two
additional characteristic numbers, £4 and ^ (£2 was defined in Section 4.5.1).
an
Indeed, property (ii) implies that the crystallographic tensors ^c^c^ d
c
J2i iotCif]CiyCid are isotropic (see Section 2.3.5 for a definition of isotropic tensors).
Moreover, two velocity vectors corresponding to two channels belonging to the
same class (see Section 2.3.4) are connected by an isometry of the crystallographic
point group G of the underlying lattice. So all velocity vectors in a given class
have the same modulus, and the quantity ^- depends only on the class index
/ of Cj. As a consequence, the fourth order tensor Yli(cV^)ci^c^ciyc^ ls a
according to property (ii). Since the three crystallographic tensors defined above
are isotropic, there exist three numbers, £2, £4 and ^ , such that:
(464)
!«-•'>
b 2
+ Say^ +
Model name D b 6
HPP 2 4 2 — —
FHP-1 2 6 3 0 0
3
FHP-2 and -3 2 7 3
165 795
GBL ^ i. 24
T9~ ~36l
FCHC-1 to -4 4 24 12 0 0
4.5 Natural parameterization of equilibria 95
^• =1 q\ q\K = 0 when K =£ K!'. The simplest and most physical choice is:
\Cba/
(4.65)
= ic2
The collisional invariants q^ and g M are the mass and the momentum of
the particles, respectively; however, g [D+1] is not the kinetic energy, but a
combination of the kinetic energy and of the mass, and is also a collisional
invariant which we call the 'effective energy'. It can be verified that the definition
of the microscopic effective energies e\ guarantees the desired orthogonality
property.
The macroscopic variables associated with these collisional invariants are
respectively the 'mass density per node' p, the 'momentum density per node'
j , and the 'effective energy density per node' w, (a linear combination of the
kinetic energy and mass densities, which we shall call 'energy' for short):
b
ja = p W = ]T Ciafu (a = 1,..., D\ (4.66)
i=i
We study nearly equally distributed equilibria, and we proceed along the same
lines as in the previous section. Here the values of the macroscopic variables for
the unperturbed equilibrium are:
= 0, a = l,...,D, (4.67)
pV+1] =
MO)
96 4 Equilibrium statistical mechanics
= 0, a = l,...,D, (4-68)
We now look for nearly equally distributed equilibrium states with mass density
p, momentum density j = 0 + rj\{l), and energy density w = 0 + rjw{l), where rj is
an expansion (small) parameter. The perturbations (pf^, K = 1,...,<5) are now:
= j m a , a=l,...,D, (4.69)
The values of the expansion coefficients /o,-, / n and fn are the same as given
in (4.54), and the d x S real matrix M defined in (4.46) reads:
(b 0 . . 0 0
0 . 0 0
&
:
= - d ( l - d) (4.70)
0 0 . . £2 0
Vo 0 . . 0
The choice of an orthogonal basis for the space of coUisional invariants has the
beneficial consequence to render M diagonal. In addition, all diagonal elements
are non-zero, since the coUisional invariants g M are linearly independent. The
particular choice for the definition of £4 in (4.64) is now seen to be justified as
we obtain simple expressions for the coefficients of M. The solution to the first
order problem (4.45) is simply:
« = 0.
d(l-d) U'
The computation of the second order perturbations /lf*] requires the evaluation
of the r.h.s. of the second order equation (4.47); the algebra is given in the
4.5 Natural parameterization of equilibria 97
l 2d
AM = ~ (-L f + - U 2 ) ,
M = 2
rf)2(^W(1)J(1)^' a=1
'--" D ' <472)
!d / 2 :2 ,
With (4.71) and (4.72), the final expression for the perturbed mean populations
fi is obtained in terms of the macroscopic variables p = bd, j and w:
P
f -
1 1
+ -rJaCia + T-wet
S2 54
(4.73)
1
« +- w e
This is the type of expression for the mean populations given in terms of
the macroscopic variables (mass, momentum and energy densities) that we
shall use to compute the fluxes of the same conserved quantities, in order to
obtain the macrodynamic equations of the lattice gas. This result is also of
practical interest for fluid dynamics simulations. Indeed (as in computational
fluid dynamics (CFD) methods) the system must be initialized with values of
the numerical variables which represent the desired physical initial situation, e.g.
describing a field of given mass density, velocity and temperature. To initialize
the Boolean field in a lattice gas simulation, we must decide the random
occupation level (0 or 1) of each channel at each node, with the constraint that
the resulting macroscopic fields match the desired initial physical situation. An
efficient way is to use the approximate equilibrium mean populations in (4.74)
as the probabilities for the random distribution of 0s and Is; this procedure
offers the advantage that the initial state of the Boolean field is close to a (local)
equilibrium state, so that spurious transient regimes are avoided (or at least
considerably shortened).
Lattice gases with single (unitary) velocity modulus (such as the FHP and FCHC
models) are well suited for the investigation of a large class of hydrodynamic
phenomena. They are by nature non-thermal models: mass and momentum
are the conserved physical quantities, and energy conservation is simply a
consequence of mass conservation. So the notion of temperature is irrelevant
in lattice gases with single non-zero velocity modulus. For a thermal lattice
gas, obviously we must have a multi-speed model. While this requirement is
a necessary condition (as fulfilled e.g. by the GBL model), the definition of a
temperature in the lattice gas does not follow straightforwardly, and in fact it is
a rather subtle matter to establish a lattice gas thermodynamics.15
In classical statistical thermodynamics, temperature T is introduced via the
canonical distribution function where exp(—jgjf) (with J f the total energy)
contains the quantity j3 = (/CBT)" 1 (with /CB the Boltzmann constant). Can we
legitimately expect a similar relation in the thermal lattice gas? One of the
complications arises from the exclusion principle (as in the Fermi gas) which,
in the lattice gas, imposes an upper bound to the energy. Consider for instance
the equipartition:
p J d V p e / J dPe\ (4.75)
15
This problem was considered by Ernst (1991), Cercignani (1994), and Grosfils (1994).
4.6 Statistical thermodynamics 99
whose analogue for the general class of models with a finite number of discrete
velocities (including the class of lattice gases) would read:
^ ^ / - ^ , (4.76)
where the sum is taken over all velocity channels. To yield a standard definition
of the temperature, the above expression should be satisfied for any value of /?,
which is clearly not the case when i is discrete and finite.
Let us consider a volume element in the lattice gas at thermodynamic equi-
librium; the size of the volume element - the region <£ of the lattice - is large
compared to vo, the elementary volume associated to a node. So we shall use the
grand canonical ensemble. We start from the probability density in phase space,
(4.38):
b
= n Y[ftq)Si{r*y;q)Si{t*\ (4.77)
where ffq) is the Fermi-Dirac distribution (4.34):
s \ r / $
(4.78)
and with the normalization: X^SGr^(eq)(^) = 1> which is easily verified with
(4.77) since the Sts take the value 0 or 1. We rewrite (4.77) as:
where
(4.80)
J l
' K=\
where et is the energy of channel i, and a and j8 are the conjugate thermodynamic
variables to be determined.16 Then we can rewrite ^{eq)(S) as:
b
0
*"^ 5 '^, (4.82)
16
We assume that there are no spurious invariants, and the variable conjugate to the
momentum does not appear since the system is at global equilibrium.
100 4 Equilibrium statistical mechanics
We can commute the sum and the products provided we take exactly the S/S that
are compatible with the state S eT; since here Si = 0 or 1, H is the Fermi-Dirac
grand partition function given by:
b
r*ei? 1=1
It will be convenient to define:
Q = log HFD
r*,i
53z«-""), (4.85)
where we have taken into account that the Jf nodes in the volume element
V (i.e. in if) have identical energy levels et. From the knowledge of the grand
partition function, we can obtain the basic thermodynamic quantities of the
system. Let us first compute the average number of particles in V'P
Since (N) = J^Yli (ni(r*))> the average density per channel is given by:
17
Note that it is not essential to indicate explicitly that the partial derivatives in the
thermodynamic expressions are to be taken at constant V since, in the lattice gas,
expansion or compression cannot be performed as in classical thermodynamics; here ££
is fixed and V is necessarily constant.
18
For simplicity, we shall omit the superscript in /| eq) .
4.6 Statistical thermodynamics 101
ft = e*-pei
We find:
and
and using (4.85), (4.86), and (4.89), it is easy to show that (4.91) can be written
as:
which is consistent with the classical thermodynamic expression for the entropy
provided a = P\x\ the consistency criterion determines \i as the chemical potential
(and consequently z = Qxp(pfi) as the fugacity), and we have:
= 0+B-i\(?Q\
explicit relation can be obtained for p(p, T) in the lattice gas. Indeed, using
(4.87), let us rewrite Equation (4.90) as:
= ze~Pet9 (4.95)
li
p
We have thus retrieved Boltzmann statistics in the low density limit. Notice
however that the status of the partition function Yli e~^6i *s different in lattice
gases where - despite the fact that there are no non-local interactions - particles
subject to the exclusion principle are not independent, and (4.85)) does not lead
to the classical result:
The reason is that in classical statistical mechanics one sums over indistinguish-
able particles (hence the factor 1/iV! in (4.97)), whereas in lattice gases one sums
over the lattice sites.
Another difference with classical statistical thermodynamics is the equation
of state of the ideal gas which we are now in a position to compute from the
above results. Using (4.90):
we obtain:
m= rlY,^YtfiT- (4-100)
m=l i
For the reason discussed above, the pressure of the 'ideal lattice gas' is larger
than the intuitively expected value which would read:
^2Tf- (4-104)
Equation (4.103) is analogous to the virial expansion of classical statistical
mechanics. However, the lattice gas is not a Newtonian system :20 the notion of
interaction potential is precluded in the classical mechanical sense, 21 and there
is no virial theorem; therefore (4.103) is not stricto sensu a virial expansion.
Nevertheless, we shall see in the next section that the coefficient B^ can be
expressed in terms of a correlation function and takes a form equivalent to the
expression of the classical second virial coefficient. In fact, at this point, we
can notice an interesting analogy. For lattice gases with single non-zero velocity
modulus, /j- e q ) = d, and (4.104) can also be written as pBi = \ d ; for the hard
disk gas, one has pBi = \ Nna2/V, where a is the hard disk diameter (see e.g.
McQuarrie, 1976). So in both cases, the result is one half the relative average
'volume' where there is one particle of the gas.
The kinetic pressure PK, defined as the density of momentum transfer:
(4105)
i
§^'
i
with D the space dimension, is different from the thermodynamic pressure p,
(4.94). To illustrate the point, let us define:
i
20
T h e lattice gas dynamics is n o t governed by Newton's equations of motion.
21
Models have been developed where non-local interactions are incorporated in the
d y n a m i c s ; one can then define a quantity which is identified as the potential of mean
force of the lattice gas (Tribel a n d Boon, 1995).
104 4 Equilibrium statistical mechanics
where the summation is taken over all velocity channels; for non-thermal models
(with |c| = 1), (4.106) reduces to c2s = 1/D, and (4.105) becomes:
Now setting jS"1 = c2s in (4.90), we find for the thermodynamic pressure:
In the low density limit (p —> 0), the thermodynamic pressure (p ~ /J" 1 p, see
(4.103)) is equal to the kinetic pressure (see (4.107)), and
(4109)
»-?5-S
is the compressibility coefficient.22 In this limit, and with P~{ = c2, one retrieves
the law of ideal gases p = p/J" 1 .
The difference between thermodynamic and kinetic pressures in thermal lattice
gases is illustrated for the GBL model (see Section 3.6) in Figure 4.2; the kinetic
pressure for the same model is shown as a function of 6 = exp(—jS/2) in
Figure 4.3.
Some authors (Chen, Zhao and Doolen, 1989, and Molvig et a/., 1989) have
used the expression PK = P ^ B ^ K , which defines a kinetic temperature TR. AS
suggested by Figure 4.3 the kinetic temperature is a complicated function of
p and of [I, and must be distinguished from the thermodynamic temperature
T = (/CBJS)"1. To discuss this point we consider the low density limit where
the mean occupation of the energy channels et is small and is given by (4.95);
the condition z -> 0 is realized for the thermodynamic state p -> 0 at P or
15.0
22
The notation cs used here will be justified in Chapter 9 where we shall see that cs is the
speed of sound propagation, which is consistent with (4.109).
4.7 Static correlation functions 105
0.1
0.0 0.2
106 4 Equilibrium statistical mechanics
where (5rcz(r*) = fy(r*) — /*, or equivalently by the static structure factor S(k),
the space-Fourier transform of (4.110):
pS(k) =
r*,r*'
r*,r*' ij
= log(exp{xn,}> = 2 ^ ^ - , (4.112)
with Qi = log(l + ea~Pei) = log(l + ea~Pe{c)). The relations between the first
cumulants and moments of interest here are (Abramowitz and Stegun, 1965):
*?] = (n,) = U
K? = (to-fif),
K?] = (to-fif),
*?] = (to - ft)4) - Hto - ft)2)2 • (4.H4)
Note that, for single speed models, (n,-) = /,- = d9 in which case the subscripts
of the KS are redundant. With (4.113) and (4.114), some straightforward algebra
yields the results: 23
({8m)3) = /,.(l-/0(l-2/,),
4
((3m) ) = /,(1 - ft)(l - 3/,- + 3ff). (4.115)
23
Higher order cumulants can be obtained from a general expression (see Appendix A. 11),
but in practice we shall not need explicit expressions beyond the fourth moment. Notice
also that (n,-) = ft and {{Snt)2) = ft(l — ft) are classical results for the mean and the
variance of the coin-tossing experiment: n,- = 0 or 1 (see e.g. Papoulis, 1965).
4.7 Static correlation functions 107
(See also (4.42) and (4.43).) It then follows from (4.111) and (4.114) that the
static structure factor is given by:
Ps(k) = j>< 2) = E ^ 1
- fi)- ( 4 - 117 )
i i
For LGAs with single velocity modulus (/,- = d), (4.117) reduces to:
> (4.119)
c i i
B2 = \p-l{\-S(k)), (4.120)
f
where:
24
Rigorously pS(k) = £ ] . icf\l — d(k)), which shows t h a t exactly a t |k| = 0 (i.e. at infinite
wavelength) 5 ( k ) = 0 (Forster, 1975). But for all practical purposes the factor (1 - <5(k))
can b e ignored.
25
L G A s with non-local interactions exhibit a k - d e p e n d e n t static structure factor reflecting
the effect of the interaction r a n g e (Tribel a n d Boon, 1995).
26
In classical statistical mechanics, Bi is expressed in terms of the radial distribution
function g(r), which is related to the static structure factor (see H u a n g , 1963, a n d Boon
and Yip, 1980): B2 = - \ /0°° dh(g(r) - 1) = \ p" 1 (1 - lim^ 0 S(k)).
108 4 Equilibrium statistical mechanics
we find:
?) (4124)
-
So the isothermal compressibility is given by:
27
For lattice gases whose microscopic dynamics includes non-local interactions (Tribel and
Boon, 1995), the limit k - • 0 must also be taken in (4.126).
Chapter 5
109
110 5 Macrodynamics: Chapman-Enskog method
Consider a lattice gas with d collisional invariants q^K\ K = 1,...,S, and sat-
isfying the semi-detailed balance condition so as to be in accordance with
the requirements for the existence of factorized global equilibrium states (see
Chapter 4).
Its macrostate at time U = 0 is described by a probability distribution function
3P (see Section 4.1.1) with the following properties:
The macroscopic fields p M (r*) vary slowly in space, which means that the
typical distance over which the macroscopic fields vary significantly is
much larger than the microscopic unit length, namely the mesh size.
2
Of course, the situation is quite different in the presence of an external force (e.g.
pressure gradient, gravity) which maintains the system out of equilibrium.
5.2 The multi-scale expansion for macrodynamics 111
where ^ (0) is just ^(loc), and ^ (1) and ^ (2) are the first and second order deviations
from local equilibrium.
The average populations /,• can be expanded in the same way:
SeT
SeT
SeT
112 5 Macrodynamics: Chapman-Enskog method
That the actual solution /,• and its lowest order approximation / / 0 ) must lead
to the same macroscopic fields implies that the first and second order corrections
must not contribute to the macroscopic variables. Thus, the corrections / / ^ and
/ / 2 ) must verify:
and (54)
Ser i=l
a n d
(5.5)
= 0, K = 1,...,5.
SeT i=l
ri=er, (5.6)
where r is the continuous version of the discrete microscopic space variable r*.
Rescaling r into ri amounts to choosing e"1 as the new length unit, instead of
the lattice constant which is the natural length unit for microscopic phenomena.
When the macroscopic variables are inhomogeneous, with gradients of order
e"1, then, by physical intuition (based on our experience with real fluids), we
anticipate two generic classes of macroscopic phenomena, superimposed on the
underlying microscopic evolution: (i) non-linear and pressure effects such as
acoustic propagation, and (ii) linear diffusive effects like viscous damping.
For a real Newtonian fluid, with given sound speed and viscosity (see e.g.
Landau and Lifschitz, 1987), non-linear and pressure effects involve first order
space derivatives. So, for inhomogeneous macroscopic fields with length scales
of order e"1, the typical time scales over which non-linear and pressure effects
occur should be of order e~l. This is particularly clear for sound waves: a
5.2 The multi-scale expansion for macrodynamics 113
On the other hand, diffusive effects involve second order space derivatives.
Thus, the inhomogeneities with typical length scales of order e~{ undergo slow
damping over time scales of order e~2. This is the case for sound attenuation in
viscous fluids (viscous damping of sound waves).
From these physical arguments, we anticipate similar macroscopic dynamical
processes in lattice gases and we introduce two auxiliary time variables t\ and
t2:
At this point t\ and t2 are not two different time variables, but simply two
expressions of the same time variable t with different scalings. To illustrate the
argument, imagine a lattice gas whose (microscopic) time step would be one
second. Then the second is the natural time unit for microscopic phenomena.
Suppose now that the lattice gas is in an inhomogeneous macrostate with a
space scale separation parameter e of 1/60. Then, t\ is just the time expressed
in minutes and t2 the time expressed in hours. For this lattice gas, the acoustic
phenomena would evolve on time scales of a few minutes, and thus are better
described in terms of a time variable expressed in minutes; diffusive effects
would only be noticeable after hours, and therefore should be analyzed with a
time variable expressed in hours.
Note that t\ and t2 are not independent variables: they are connected through
their relations with t. But it is usual in multi-scale methods (see Bender and
Orszag, 1978) to consider any function F(r, t) of space and time as a function
F{*ututi) of ri, t\ and t2, as if t\ and £2 were independent variables. 3 The time
and space derivatives of an arbitrary function F(r, t) can thus be expressed as
follows in terms of the derivatives of F(ri, t\, t2):
dF dF a n dJ
dF dF 2dF
-x-
dr
= e dfx
r- ^~
dt
= €~^-
dti
+ e dt
^--
2
All the tools necessary for the multi-scale expansion are now ready; it remains
to introduce the master equation to which the expansion will be applied.
3
Strictly (from the mathematical point of view) the two functions should be labeled
differently, but, as no ambiguity results, we keep the same notation.
114 5 Macrodynamics: Chapman-Enskog method
b (5.9)
Si-sdZis^Hn/in/L
(s,s')ey2 ;=1
Here, nj denotes the Boolean population after collision. Note that the Boolean
populations n\ and n, in the second equation are all taken at node r* and
time U. This formulation which separates the collision phase and the propagation
phase, is strictly equivalent to (2.13), but conceptually easier to handle for the
forthcoming averaging operations.
Ensemble-averaging the first equation in (5.9) with the probability distribution
8P of the corresponding macrostate gives:
/,(r*+c I -,f* + l) = /;(r*,t*), (5.10)
which merely expresses that the propagation (or free-streaming) phase just
moves information, without loss or gain, from one node to another node. The
second equation in (5.9), after ensemble-averaging over £P and over all possible
choices for non-deterministic collisions, yields:
Indeed, the random variables £(s^sf) involved in the collision phase (see Sec-
tion 2.2.3) are statistically independent from the Boolean field. So, the averaging
over all possible choices for non-deterministic collisions simply replaces the
random variables £(s->s') by their mean values A(s-*sf) (the collisional transition
probabilities). The brackets (• • •) denote ensemble-averaging over the probability
distribution ^ , that is, summing over all possible realizations S of the Boolean
field n, with each term weighted by its probability ^(S,U)\
b
Ai(<P, r*, u) = ^ 9{S) ] T (s[ - Si)A(s^sf) J J S,W Sj(up.
SeT (s,s')ey j=l
Since the product Y{SJ{T*) JSJ{T*) J is zero unless S(u) = s, the averaged collision
S S
Without any further hypothesis, the averaged collision term A,- depends a
priori on the full probability distribution ^ , and not just on the mean popula-
tions /j(r*,f*). Indeed, the average of the product II/=i n/ 7 '"/ 7 ' ^n (5-9) cannot
be factorized into the product of the averages, since the populations of different
channels at the same node may be correlated. Consequently, a closed form for
the averaged microdynamic equation cannot be obtained without the Boltzmann
approximation (see Chapter 4), which assumes that particles entering a collision
are uncorrelated; then the average of the product Il;=i njSjn/j c a n be approx-
imated by the product of the averages. With this ansatz, we obtain the lattice
Boltzmann equation as the averaged microdynamic equation (see Chapter 4).
Here we postpone the introduction of the Boltzmann approximation and keep
the full expression (5.13) for the averaged collision term.
Aj(^, r*,£*) is clearly linear in 3P\ it is thus a linear function of 2bmjr real
variables, the probabilities 0>(S, U) of each of the 2hmjr possible configurations
S of the full lattice gas at time U. We rewrite (5.13) as:
\ (5.14)
SeT
To derive the macrodynamics of the lattice gas we thus start from the following
equation for the averaged microdynamics:
Obviously the average collision term A, vanishes for any global equilibrium
)
. It is less obvious, but nonetheless true, that it also vanishes for any local
equilibrium ^(loc). Indeed, the collision process in a lattice gas, as defined in
Section 1.2.3, is purely local. Therefore, the value Aj(^,r*) of the averaged
collision term at node r* and time U depends only on the single node probability
distribution for node r* at time £*, and not on the full information contained
in 3P. On the other hand, the definition of a local equilibrium (see Section 5.1)
implies that for each node r* in j£f, there exists a global equilibrium ^ (eq) such
that the single node distribution extracted from ^ (eq) is identical to the single
node distribution extracted from ^(loc) at r*. Since all global equilibria make A,
vanish, it is clear that A,- also vanishes for local equilibria:
It is important to note that the conservation laws imply some constraints for
A,-. Indeed (as in Section 4.4.1) consider a basis (q^K\ K = 1,...,<5) of the vector
subspace of collisional invariants. Using (5.13) for Aj(^,r*), and Equation (2.19)
116 5 Macrodynamics: Chapman-Enskog method
which means that the rank of the b x 2bmjr matrix LiS is smaller than fo, and is
at best equal to b — S, 6 being the number of independent collisional invariants.
In other words, the subspace of collisional invariants is included in the kernel
of the adjoint of LiS. This property will play an important role in the derivation
of the macrodynamics of the lattice gas.
When the scale separation parameter e is small (i.e. the averaged populations
vary little from one node to the next), we can validly approximate the averaged
population /j(r* + c^<, U + 1) in Equation (5.16), by its Taylor expansion around
(5.20)
+ (9(e3).
Here, /,- and its derivatives are taken implicitly at node r* and time U (unless
otherwise specified). To obtain this result, we have used Equation (5.8), which
5.2 The multi-scale expansion for macrodynamics 117
connects the derivatives with respect to the lattice space and time variables
(r*, £*), and the derivatives with respect to the rescaled (macroscopic) space and
time variables (ri,ti,^).
We now introduce the fact that the average populations /, are close to their
local equilibrium values; so they can also be expanded in powers of e - see (5.2)
- with a zeroth order term corresponding to the local equilibrium. Combining
the two expansions, the l.h.s. of (5.16) becomes:
(5.21)
Ser
) (5-22)
SeT
33
+ (9(e )).
Using (5.21) and (5.22) as the respective expansions of the l.h.s. and r.h.s. of the
averaged microdynamic equation (5.16), we obtain, to the first order in e:
SeT
and to second order:
We shall now derive the first and second order macrodynamical equations of
the lattice gas starting from Equations (5.23) and (5.24) which follow from
the averaged microdynamic equation (5.16) expanded in powers of the scale
separation parameter e.
118 5 Macrodynamics: Chapman-Enskog method
Under this form, the solvability condition (5.26) for the first order problem (5.23)
appears clearly as a set of macroscopic continuity equations connecting densities
p M and fluxes <DM of the microscopically conserved observables (collisional
invariants) g M .
We now need to evaluate these fluxes and to express them as functions of
the macroscopic densities p M , so that the continuity equations (5.26) become a
closed set of evolution equations for the macroscopic variables, i.e. the densities
4
See Section 2.1.3 for the definition of macroscopic densities and fluxes, and Section 4.1.2
for the notion of ensemble-averaging.
5.3 First order macrodynamics 119
f.(0) _ (L . _ i / c . . _* .
0 C2 C4
b(b-2p)
T2 J*J
2p{b-p) S2 (5.27)
X \SapSyS + SayS + d d h
Ac?
^ w T )
X (<><xB(>yd
aj + Say SRS +
A
120 5 Macrodynamics: Chapman-Enskog method
These coefficients depend only on the geometric structure of the lattice gas and
not on the collision rules.5
We now use (5.27) to obtain the averaged mass, momentum and effective
energy fluxes O [0] , O M and O [D+1] , involved in the first order macrodynamic
equations (5.25). The details of the calculations are given in the Appendix,
Section A.9.
(i) The mass flux is just the momentum density; its components are:
b
Qf = ^ <W,-(0) = Jp = f>up, P = h...,D, (5.29)
where
and
i=\
^ 2£4 . (5.33)
+ J
where
5
See Section 4.5.2 and Table 4.1 for the values of these coefficients for some simple
models.
5.3 First order macrodynamics 121
The density-dependent coefficient g(p) that appears in both (5.30) and (5.32)
is called the 'non-Galilean factor', or, more precisely, the 'non-Galilean factor
for momentum'. 6 This non-Galilean factor is very important both for theory and
for applications of lattice gases. Indeed, Equation (5.30) is not Galilean invariant
unless g(p) is exactly equal to 1, which a priori is not the case. However, as
we shall see in Chapter 9, Galilean invariance can be recovered, at least for
non-thermal models, in incompressible regimes where the density p is nearly
constant: g(p) is then a constant factor that can be eliminated by proper time
rescaling. The factor g\p) in Equation (5.33) is the 'non-Galilean factor for
energy' because it plays, for the energy, the same role as the non-Galilean factor
g(p) in the momentum equation.
Somewhat in the same way as transport coefficients are expressed in terms of
time-correlation functions, these non-Galilean factors may also be given a mi-
croscopic content: they can be expressed in terms of static correlation functions
(for details see Brito, Bussemaker and Ernst, 1992). For LGAs satisfying detailed
or semi-detailed balance - so possessing a universal equilibrium distribution -
the non-Galilean factors are given in terms of the correlation functions K^
and 7c|3) discussed in Section 4.7; so, in some way, they can be interpreted as a
thermodynamic property of the lattice gas.
We now return to the first order continuity equations (5.26) which, with the
above results, we can write in the standard Eulerian form:
09 (5.35)
The mass continuity equation is valid to any order in j , whereas the momentum
and effective energy continuity equations are approximations neglecting terms
(9(rj3). The reason is that the mass flux is also the momentum density, but
the momentum and effective energy fluxes are functions of the macroscopic
variables, which can be obtained only through an asymptotic expansion in
powers of j and w.
ser
+ (9(erj2).
It is relatively easy to convince oneself that the solution ^ (1) (S) of the first order
problem must be a linear combination of the gradients diyj$ + d^jy and d^w.
The first order corrections f^ to the mean populations, that result from ^ (1) (S),
must thus also be a linear combination of these gradients. The most general
form for this combination reads:
fP = M ^ ^ l + M + NiAyWy (5.39)
where N,- is a set of i-indexed D -dimensional vectors, and where M,- is a
set of i-indexed D -dimensional second order tensors that may depend on the
macroscopic variables (p,j 2 , w). Furthermore, as the second order tensors M,- are
contracted with the symmetric expression diyjs + d^jy, they can be assumed,
without loss of generality, to be symmetrical with respect to exchanges of their
first and second Greek indices.
So far, we have made no hypothesis on the value of these vectors and tensors,
and they can be regarded as a set of unknown quantities depending on the
macroscopic variables. However, we can be more specific if the lattice gas is
G-invariant and verifies the irreducibility property (see Sections 2.3.4 and 2.3.6).
Then, we can apply the simplification procedure described in Section 2.3.6.3:
the vectors N,- must be proportional to the velocity vectors c,-, with equal
proportionality coefficients for all channels belonging to the same class :7
where / is the index of the class to which channel i belongs. Along the same
lines of reasoning, the second order tensors M, are linear combinations of c,- ® c,-
and (5, the Kronecker tensor, with coefficients which depend only on the class
of the corresponding channel:
This simplification is possible because the vectors and tensors M,- and N,- are
such that for all isometry g in G and for all i = 1,..., b, we have g(M,-) = Mg(,-)
and g(N,-) = Ng(,-) (for the definitions of g, G and g, see Section 2.3.4).
The resulting form for f^ reads:
d
fP = (Minces+M\i)dy3) «»+ dM
iy diyw.
that ft^ does not contribute to the momentum variable imposes that:
b
i5=0, Vy,5, (5.44)
and that f^ does not contribute to the effective energy requires that:
b
^2(M(i)ciyciS + M\i)5yS)ei = 0, V 7,8. (5.45)
The // (1) s contribute neither to the effective energy density nor to the mass
density, thus, they do not contribute to the kinetic energy density either, since
the latter is a linear combination of effective energy and mass densities. So, the
following relation holds between M and Mr:
= 0, Vy,<S. (5.46)
1=1
This relation, taken for 7 = 6 and summed over 7 gives the following
124 5 Macrodynamics: Chapman-Enskog method
simpler form:
i=\
(5.48)
The lattice gas model is assumed to verify the fourth order crystallographic
isotropy property (see Section 2.3.5). Thus, there exist two coefficients v(c) and
v^c\ depending on the macroscopic variables (p,j 2 ,w), such that:
ipCiyCis = —v (c)
+8as8py),
and (5-49)
b
= - v ( c ) 5ap.
The superscript '(<?)' is used to indicate that these coefficients come from the
first order perturbation f^ to the mean populations, which are solutions of the
first order problem, and thus depend on the collision rules. They are not purely
geometrical coefficients like £2, £4 and £5.
5.3 First order macrodynamics 125
and (5.50)
1=1
= 0. (5.51)
With this result, one can express the first order contribution to the momentum
flux tensor (5.48) in terms of the coefficient v(c^:
1=1
( . (5.52)
\ /
which yields immediately:
Thanks to (5.51) which connects v(c) and v^, the first order contribution to
the momentum flux tensor takes the simple form (5.53), involving only the
symmetric traceless part of the momentum density gradient tensor d j .
In a similar manner, the first order contribution J2i eictpfi^ t 0 the effective
energy flux is given by:
(5-54)
Observing that e\ depends only on the class to which channel i belongs, and
invoking again the crystallographic isotropy property, we see that there exists a
coefficient ^c\ which depend on the macroscopic variables (p,j2,w), such that:
(5.55)
126 5 Macrodynamics: Chapman-Enskog method
The resulting expression for the first order contribution to the effective energy
flux vector reads:
b
We emphasize again that the actual values of v(c) and £(c) are still unknown. At
this point, all we can say - without the Boltzmann approximation - is that these
coefficients do exist, and that the first order contributions to the momentum
and effective energy fluxes are of the form (5.53) and (5.56). In Section 5.6, we
shall use the Boltzmann approximation to compute these coefficients explicitly.
Their numerical value can also be 'measured' by computer simulations of wave
damping (see Chapter 10).
The second order problem (5.24) has exactly the same mathematical structure
as the first order problem (5.23), except that the r.h.s. is more complicated.
The analysis will also involve solvability conditions, that is, the r.h.s. must be
orthogonal to all vectors in the kernel of the adjoint of the linear operator L. s
(Fredholm's criterion).
i=\
wy (5.57)
1=1 i=l
for all K = 1,...,<5, that is, for all collisional invariants. Since we treat the
case of single-species thermal lattice gases, the linear invariants are the mass,
the momentum and the effective energy. The algebraic manipulations leading
from (5.57) to equations (5.58), (5.59) and (5.61) hereafter are given in the
Appendix, Section A. 10. In these manipulations, we use the first order macrody-
namical equations, (5.35) to (5.37), and the fact that the first order perturbations
5.5 The macrodynamic equations 127
dt2p = 0, (5.58)
(1)
which shows that there is no mass diffusion. The first order perturbations /; do
not appear in the second order solvability condition (5.58). Indeed, by definition,
the /i (1) s contribute neither to the mass density, nor to the mass flux, since the
latter is also the momentum density (see Equation (5.4)).
The case of the momentum invariant is more difficult, because the first order
perturbations /; (1) now contribute to the second order solvability condition
through the momentum flux J2iciocctpfi^ ( s e e Equation (5.53)). The resulting
solvability condition is:
where the coefficient v, which a priori depends on (p,j 2 , w), amounts to:
) (5.61)
where the coefficient £ is:
We can now merge the results of the first order problem, (5.35), (5.36) and (5.37),
and of the second order problem, (5.58), (5.59) and (5.61), to obtain the macro-
dynamical equations governing the time-evolution of the macroscopic variables
128 5 Macrodynamics: Chapman-Enskog method
dp(&jp) = dpUdpw)
V J V 7
(5.65)
(9{e2ri2)
The non-Galilean factors g and g', defined in (5.31) and (5.34) respectively,
depend on the mass density p. The pressure P is defined in (5.32), and depends
on all the macroscopic parameters (p,j2,w). The transport coefficients v and £
are defined in (5.60) and (5.62). Their values also depend a priori on all the
macroscopic parameters (p,j2, w). However, any dependence of v or £ on j 2 or w
would lead to terms of order ofe2rj2 or higher in the macrodynamic equations. To
be consistent with the order of approximation in the macrodynamic equations,
these terms will be neglected; from now on, we shall use the short-hand notation
v and £ for v(p,0,0) and £(p,0,0).
Equations (5.63), (5.64) and (5.65) give the general form of the macrodynamic
equations governing the long-time, long-wavelength behavior of the lattice gas.8
These macrodynamic equations show close resemblance with the equations
of continuous fluid dynamics (such as the Navier-Stokes equation). However,
there is a difference: the non-Galilean factors g(p) and g'(p), in (5.64) and (5.65)
respectively, do not appear in the usual equations of continuous fluid mechanics.
The lack of Galilean invariance in the fluid-like behavior of the lattice gas is
a direct consequence of the discrete nature of the microscopic velocity space,
which itself results from the lattice structure.9 In Chapter 9, we shall show
how Galilean invariance can be recovered in hydrodynamic regimes such as
incompressible flow and linear acoustic flow.
8
In Chapter 8, we give an alternative derivation of the lattice gas hydrodynamic
equations using the lattice Boltzmann approximation.
9
A factor also appears in front of the advection term in the hydrodynamic equations of
discrete kinetic models (Gatignol, 1975) where position space is continuous and velocity
space is discrete.
5.6 Transport coefficients within the Boltzmann approximation 129
jj ^ w \ (5.68)
7=1
where Q,7 is the Boltzmann collision matrix, that is, the Boltzmann collision
operator linearized around the uniform, zero-speed, zero-effective energy equi-
librium populations /; (0) = d = p/b:
% = (5.69)
One may wonder why the linearization is performed around the zero-speed and
zero-effective energy uniform equilibrium /,• = d, rather than around /, = f^°\
as it logically should. The reason is that the difference between the Boltzmann
10
The Boltzmann equation was discussed in detail in Section 4.2.2; we should keep in
mind, it is an approximate equation: strictly, in (5.67) we should write ' ^ ' rather than '='.
130 5 Macrodynamics: Chapman-Enskog method
s
i i ) A ( s ^ s l SJ dP l
\ ~ d)bP'> ( 5 - 7 °)
s,s'ey
1=1
Expression (5.70) provides the basis of a procedure for the computation of the
values of the matrix coefficients Q,7, once the collision rules A(s^s') of the lattice
gas model are given.
The Boltzmann collision matrix is not regular, because of the existence of
conserved quantities. Yet, Equation (5.68) has a family of solutions provided
the first order solvability conditions, (5.35) to (5.37), are satisfied. The desired
solution // (1) will be the only element in this family that also satisfies the
condition to have zero contribution to the mass, momentum and effective
energy (see Equation (5.4) and the discussion preceding (5.4)).
Using the solvability conditions together with the expression (5.27) for the
local equilibria fp\ one can re-express the r.h.s. of (5.68) in terms of the
gradients of the macroscopic variables p, j and w as:
E diyj5 + dl5jy
C
O f (1) ^ (, n ^
j j = y 5
Ti ^ ~ D ' 2
(5.71)
(9(erj2).
The standard procedure for solving (5.71) is to introduce Q [~1], the 'pseudo-
inverse' of the Boltzmann collision matrix. More precisely, consider the restric-
tion Q of Boltzmann's collision matrix Q, to the subspace of R b , orthogonal to
its kernel. This restricted operator is regular and its inverse operator Q" 1 can
be defined. The 'pseudo-inverse' Q [11 of Q is defined as an operator which acts
as the null operator on the kernel of Q, and as Q" 1 , on the subspace of ]Rb,
orthogonal to the kernel.
11
The linearized Boltzmann equation is discussed in detail in Section 6.1.
5.6 Transport coefficients within the Boltzmann approximation 131
£2
l
. STct-lU 2
^\ z ( 5 - 72 )
+ (9(erj2).
(5 ?3)
2 . . -
and
(5 74
] „ ., ->
Considering Equations (5.52) and (5.56) as definitions for v(c) and £(c), the
collisional contributions to the transport coefficients v and C, we can use the
relations (5.73) and (5.74) to obtain the desired expressions for v(c) and £(c)*
1 b
1
and
4i
The resulting values for v(c) and C(c) depend on the density per channel d through
Q[~l]; these are approximate results, since we used the Boltzmann approximation
in the computation. They are important theoretical results, but in practice, it is
more convenient - and also more accurate - to 'measure' the transport coeffi-
cients using a wave relaxation technique by lattice gas automaton simulations
(see Chapter 10).
132 5 Macrodynamics: Chapman-Enskog method
Non-thermal models (see Section 4.5.1) differ from thermal models by the
conservation laws: non-thermal models conserve mass and momentum, but do
not conserve a linearly independent kinetic energy observable. Examples are the
FHP-II and FHP-III models. Consequently, the macrodynamics of non-thermal
models is governed by only D + 1 macroscopic continuity equations, while
there are D + 2 equations for thermal models. However, it would be erroneous
to think that the macrodynamics of non-thermal models simply degenerates
from the macrodynamics of thermal models, by 'omitting' the energy continuity
equation and the contributions of the effective energy density w to the mass
and momentum continuity equations. Indeed, there are substantial differences.
The first difference already appears at the stage of the nearly equally distributed
equilibrium distributions, which, for non-thermal models, read:
=
~h ~^~ ~p~CifX^a
b b - 2p , &c \ . . (5 77)
(see Equation (4.58) in Section 4.5.1). This expression cannot be obtained from
the equivalent expression (4.73) for thermal models by setting w = 0. Indeed, the
non-linear jajp term for non-thermal models involves the second order tensor
(ciaCia — ^<5aj?)> whereas the same term for thermal models involves the tensor
(cfaCia — T^ajs)- These tensors are different in the most general case, for example
for models with rest particles. This slight difference in the equilibrium distri-
bution produces modifications in the first and second order macrodynamical
equations.
(5.78)
(5.80)
The coefficient A, which is equal to 1 for thermal models (see Equation (5.32)),
is given by:
, , b-2p b /4<JU D2
g(p)= b_ £>(2) + 2)V7 2 " + T
The second order macroscopic equations for the mass and momentum densi-
ties read:
dt2p = 0, (5.82)
and
and
* > ^ (5.85)
and v' are given in (5.60) and (5.62), respectively. The values of v and V also
depend a priori on the macroscopic quantities (p,j 2 ). However, any dependence
of v or v' on j 2 would lead to terms at least of order (e2rj2) in the macrodynamic
equations. For consistency with the order of approximation in the macrodynamic
equations, these terms are neglected; for simplicity, we shall use the short-hand
notation v and v' for v(p,0) and v'(p,0).
(5.88)
2
+ (9(erj ).
and
b
1 / 1 __ r „ \
- Dv^j. (5.90)
For homokinetic models (see Section 5.7.1), the bulk viscosity coefficient van-
ishes, and the Boltzmann expression for the shear viscosity coefficient is identical
to the shear viscosity of thermal lattice gases.
5.8 Comments
Linearized hydrodynamics
137
138 6 Linearized hydrodynamics
the question whether the intrinsic fluctuations in LGA capture the essentials of
actual fluctuations in real fluids as described by linearized hydrodynamics.
7=1
with the notation /,-(r*;t*) = /(r*,c,-;t*) and where we have used the property
that the action of the collision operator on f^ is zero. The b x b matrix Q i;
denotes the linearized collision operator (5.69):
ssf
ss'
_ ^(eq)\-s; _ /i _ f(eq)\-:
X jj ) v1 j ; ;
2
The equilibrium state /(eq) is a global equilibrium /(eq)(c). Fluctuations can also be
related to external (forced) perturbations when a constraint is imposed to the system;
then /(eq) must be taken as local equilibrium /(eq)(r*,c).
6.1 The linearized Boltzmann equation 139
> i
(/ ( e q ) )(l -
Ksi — s
i) siys —* $ )
n //-(eq)\sfc/-i
k \J k ' V
_
1—\ 5
ss' ^~fj )
where the last line on the r.h.s. vanishes because of the equilibrium condition
(4.36).
For non-thermal LGAs (models with single non-zero velocity modulus) the
equilibrium distribution function (6.3) is independent of the velocity: / (eq ) = d,
the average particle density per channel. Then the linearized collision operator
reduces to (see 5.70)
ss'
2
ss'
For multi-speed models (like thermal LGAs) this symmetry property does not
hold. Nevertheless a symmetrical expression can be obtained (Ernst and Das,
1992) by considering the modified operator Slijxfp where KJ^ = / J q ( l — /J q )
is related to the static correlation function S(k) (see Section 4.7). Defining the
factorized equilibrium distribution function:
b
(which is a- and /^-dependent) and using the expression obtained above for Q I; ,
we have:
ss'
( f (eq)
s/
- ) n 7Z7(
l
k \ h
where, with the explicit expression for the equilibrium distribution function (6.3),
Here p(s) = ^2k Sk and e(s) = ^2k e(ck) Sk are respectively the number of particles
in configuration s and the corresponding energy; these quantities are conserved
140 6 Linearized hydrodynamics
P(s) = J2Sk =
k k
ss'
(A\a\B) =
ij
In the present chapter we are interested in describing the space and time
dynamics of fluctuating quantities, such as the number density, the current
density, and the energy density, in the lattice gas. The relevant objects for
evaluating the dynamics of these quantities are their autocorrelation functions,
which involve the value of the fluctuating quantity at some arbitrary initial
time, provided the system is in a stationary state. Therefore the computation
of the correlation functions can be formulated as an initial value problem, and
the projection operator technique is quite appropriate (Zwanzig, 1965; Mori,
1965a, 1965b). It provides a convenient separation of the dynamical variables
into slow variables and fast variables: the latter correspond to the kinetic modes
and they decay practically instantaneously in comparison with those variables
which vary slowly in space and time. These slow variables can be identified
when the constants of motion of the system are known; we shall see that they
correspond to the hydrodynamic modes.
3
The weight K\ is p dependent.
6.2 Slow and fast variables 141
fcjhint), (6.17)
8e(u;t*) = \
We shall describe the long-time behavior of the lattice gas in terms of </>;(r;t).
Long times are large compared to the microscopic time At over which collision
and propagation take place. So far we have omitted to write explicitly the time
interval At because it is just the automaton unitary time, but here it is important
to do so for the clarity of mathematical manipulations; then the discrete time
variable is U = nAt, where n is an integer. So the linearized Boltzmann equation
(6.4) now reads:
* + At) - <fc(r*;**)
j
its physical interpretation will be given subsequently. It is easy to check that the
operator ^ , as defined by (6.21), has the properties of a projector, i.e.
Applying this definition to 0,-(r*;£*), we project $(k,c;s) onto "To to obtain the
Laplace-Fourier transformed hydrodynamic variables, i.e.
4
^0(k,c;s) = ^ K)iV a 2 a (k;s), (6.25)
a=l
with
Ax(k;s) = dp(k;s) = J2K?)kKc;s)9 (6.26)
i
4
The factor exp(s At) on the l.h.s. of (6.30) is a consequence of discreteness. Note that by
taking the limit
/ e-stdt<t>{t) =
Jo
6.2 Slow and fast variables 143
This result exhibits an appropriate form for applying the projectors 9 and 2, to
separate the equation into a 'hydrodynamic equation':
(6.33)
) (6.37)
equation:
4
sAt
e Aa(k;s) - ^2 -S?«*4r(k;s)-
<7=1 a=\
= e^'AtA^kiO) + / a(k;s). (6.39)
where (i) the second term on the l.h.s. represents the coupling between the
hydrodynamic variables, and the matrix elements J£?ao- are given by:
)NG (6.42)
The expansion of the other quantities in (6.39) is not so trivial, and they will
now be considered separately.
~ (aa\(l-ik'cAt-±(k'cAt)2)\aa)N
At2ff
~ S^ - ikiAt<aa|J<)Na - fc/UJjl-ylJZ1)Nff. (6.46)
Here k\ (I = {x,y}) denotes the components of the wave vector k,6 and |J a ) the
current associate to |a a ), with:
(
0 ikx iky 0
0 ikx
0 iky
%% o
0 ikx(c2s-c2T) iky(c2s-c2T) 0
= Lm, (6.47)
where the constants c2 and c\ are given by:
4 - i^.
we obtain:
(6.53)
a=l
this quantity is called the subtracted current, and is orthogonal to the hydrody-
namic eigenvectors.
We now expand the resolvent operator:
One can easily verify that J£X2|<£) = A|0) because J projects any vector
onto the subspace orthogonal to the subspace spanned by the vectors with zero
eigenvalue; so in practice we can substitute £l£ll by Cl. On the other hand,
we notice from Equation (6.34) that R~{ = (exp(s At) — 1 <£) operates on the
projected vector J \(f>), and since I21 \cj>) = 11(/>), we may therefore, in (6.55),
replace operationally J 2 by the identity operator. With these simplifications, we
have
with JR given by (6.56).7 We now treat the resolvent by applying the operator
identity
l
A=1B-1B{A-B)^ (658)
we obtain
l l
R = —^ —^H(K)R, (6.60)
w
sAt-Q sAt-Q
7
Expression (6.57) has a structure similar to the memory function in continuous theory
(see e.g. Boon and Yip, 1980, Chapter 3).
6.3 The hydrodynamic limit 147
(6.62)
where Kj$ is proportional to kn, that is to lowest order in k:
Kaa ~ *:£> = -ktkm (At)2 <J< |(1 + Q) sAt{_a \Ja) Na. (6.63)
This memory function is well behaved in the limit s —> 0 and for all values of
fc, provided the non-zero eigenvalues of Cl be sufficiently different from zero at
k = 0. This 'correct' long-time behavior is a consequence of the combination
of the operator (s At — Q)"1 with the projector =2 which prevents singularities
(~ s"1) that would result from the action of Cl on the hydrodynamic variables
(which have zero eigenvalues at k = 0).
where a n (k; 0), (n = 5,...,fo),corresponds to the fast modes. The negative eigen-
values (distinct from the zero eigenvalues, n = 1,...,4) are responsible for the
fast decay and eventual vanishing of J a (k;t*). 9 So we may safely assume that in
the hydrodynamic limit:
lim lim/ a (k;t*) = 0. (6.67)
££•00
-•00 / 0
(1 + sAt + ^s2At2)Ax(k;s)
- ( 1 + sAt+ ^s2At2)AtAa(k;0)
J2 Lxa(k)AtAa(k;s)
9
The decay rate is model dependent; this point and the related question of the separation
of eigenvalues at low k will be examined in the next chapter in the light of specific LGA
model computations.
6.3 The hydrodynamic limit 149
A/-2 \fl
— \J m) + (V\ \Jm)
2 sAt — Q
IffAff{k;s) = 0. (6.71)
Again we divide this equation by At, and take the limit sAt <C 1, to obtain:
52 a (k;s) + -s
+ fe/fcm
- Q
where we have omitted the term on the r.h.s. of (6.70) because when inserted in
(6.72), it will give ^sAtA^k',0) which vanishes in the limit sAt < 1. Iteration of
(6.74) gives:
J ^ N c A ^ (6.75)
where
- <W). (6.76)
Substitution of these results into (6.72) yields:
*»*m 5 1 *„,*;«)]a,(k;s) = ^ a (k;0), (6.77)
with the memory function:
(6.79)
^ ^ | ) , ( 6 . 8 0 )
z J
s— Q
These results, (6.77) with (6.80), conclude the derivation of the linearized hydro-
dynamic equations.
The memory function expresses the effect of the past history of the dynamics
related to the fast variables. But here we are interested in the dynamics of the
slow variables Aff(k; U) which relax over times {(9(k2)) very long compared to the
characteristic decay rate of Oa(7(k;T) ~ exp(Ai), where X denotes the non-zero
(negative) eigenvalues of the collision operator. Consequently, in (6.81), we may
consider, to good approximation, the following simplifications: (i) Aff(k;t*) does
not vary significantly over the summation time, and can therefore be taken out
of the sum; and (ii) Oa<T(k;r) decays sufficiently rapidly that it is practically zero
for times larger than T ~ (9(X~l), and the summation can be extended to oo. As
a result the convolution (6.81) becomes Aa(T Aff(k;t*) with:
T=0
=E H=0
11
The limit 5 —> 0 in Laplace space corresponds to the limit t —> oo in real time.
6.5 Comments 151
6.5 Comments
set by the observer in the range of sufficiently large wavelengths. In the next
chapter we shall see that this picture, as provided by a simplified system such
as the lattice gas automaton, is in accordance with the behavior of spontaneous
fluctuations as observed in real fluids.
Chapter 7
Hydrodynamic fluctuations
The object of the present chapter is small deviations from local equilibrium
which are triggered by spontaneous fluctuations. In real fluids these fluctuations
which temporarily disturb the system from local equilibrium are such that a
fluid at global equilibrium can be viewed as a reservoir of excitations extending
over a broad range of wavelengths and frequencies from the hydrodynamic
scale down to the range of the intermolecular potential. Non-intrusive scat-
tering techniques are used to probe these fluctuations at the molecular level
(neutron scattering spectroscopy) and at the level of collective excitations (light
scattering spectroscopy) (Boon and Yip, 1980). The quantity measured by these
scattering methods is the power spectrum of density fluctuations, i.e. the dy-
namic structure factor S(k,co) which is the space- and time-Fourier transform
of the correlation function of the density fluctuations. The spectral function
iS(k, co) is important because it provides insight into the dynamical behavior
of spontaneous fluctuations (or forced fluctuations in non-equilibrium systems).
Whereas the fluctuations extend continuously from the molecular level to the
hydrodynamic scale, there are experimental and theoretical limitations to the
ranges where they can be probed and computed. Indeed, no theory provides
a fully explicit analytical description of space-time dynamics establishing the
bridge between kinetic theory and hydrodynamic theory. Scattering techniques
have limited ranges of wavelengths over which fluctuation correlations can be
probed. Numerical computational techniques can realize molecular dynamics
simulations (Boon and Yip, 1980) covering in principle the whole desired range,
but in practice there are computation time and memory requirement limitations
153
154 7 Hydrodynamic fluctuations
The dynamic structure factor is defined as the double Fourier transform with
respect to space and time of the van Hove function G(r, £•), the correlation
function of density fluctuations (5p(r*,£*) around the equilibrium state (Boon
and Yip, 1980):
00
pS(k,co) = Yl Y, e-'°"*-'kl*(5p(T*,\u\)Sp(0,0))
r* t*=—oo
y1 J2 k,0)), (7.1)
t*=—oo
where
h(k) = J2 exp(-zk • r*)/z(r*) (7.2)
The average is taken over an equilibrium ensemble and the equilibrium distri-
bution has the form of the Fermi-Dirac distribution (see Chapter 4):
• (7.4)
where drii(k, t*) is the spatial Fourier transform of <5ftj(r*, £*), and KJ^ is the equal
time correlation function of the drift (see Section 4.7). With (7.3), the dynamic
structure factor (7.1) is expressed in terms of the propagator (7.5) as:
(7.6)
and the static structure factor S(k) - the Fourier transform of the equal time
density correlation function, discussed in Section 4.7 - is obtained from (7.5)
taken at t = 0:
we want to study fluctuations which disturb the gas locally from its equilibrium
state. We proceed as in Chapter 6 by expanding the collision term AJ5) around
the stationary distribution / , (7.4):
<5n,(k,U + 1)
156 7 Hydrodynamic fluctuations
which expresses the fluctuation dnt in terms of its value at the previous time step.
Solving Equation (7.11) by iteration, we obtain Srii(k, t) from the knowledge of
its initial time value:
, (7.12)
ij
which is valid for U > 0. Here [exp(—ik • c)]; ; = <5;;-exp(—ik • c,-) is a diagonal
matrix. Substituting (7.12) in (7.5) and noting from (7.5) that r i ; (k,0) = 5iJ9 we
obtain the expression for the propagator:
t=0
r In,
[ e s + i k c _ j[ _ Q J .j* ' W-1^
where the second equality follows from X^JLo*" = (1 — x) —1 , |x| < 1. After
summation over (i9j)9 we rewrite (7.14) equivalently as:
ijl
then we make use of the relation J \ QjKj = 0 (which follows from mass
conservation) to obtain:
E iVM*? = E 42) + E V
The dynamic structure factor iSf(k,co) is now readily obtained from (7.6) by
setting s = ico in (7.16):
pS(k,co) = R ^ £
t=—oo ij
Y ^ ™ + pS(k); (7.19)
this is the spectrum for the ballistic regime. On the other hand, when one is
interested in the hydrodynamic regime, i.e. when the system is probed with a
wavelength large compared to the microscopic scale (the inter-particle distance
in real fluids; the lattice unit length in LGAs) and large compared to the mean
free path between collisions (this condition prevails for dilute fluids), then one
can compute S(k,co) analytically from (7.18) in the limit of weak spatial and
temporal variations (|k| —> 0, co -> 0).1 The hydrodynamic structure factor is
the subject of Section 7.4.
ijKfaniCj) = 0, (7.20)
are the collisional invariants.3 We will also need the left eigenvectors with zero
eigenvalue of the non-symmetric matrix O. They follow from the conservation
laws for the total number of particles, the total energy and the total momentum:
£,-an(c0A<*)({n7-}) = 0,or
^ a n ( C l ) Q l 7 = 0. (7.22)
i
where \\p^) is the b-component vector |i/^)/ = K^xp^Ct) (i = 1,...,&), and z/i(k)
denotes the associated eigenvalue. The zero eigenvectors an(c) of O represent
the right eigenvectors xp^ in the long wavelength limit (k —• 0). Denoting the
eigenvectors by:
M V ^ = <V (7.29)
3
(7.20) follows directly from the observation that for a stationary distribution with
Lagrange multipliers bn (and also with bn + Sbn), A^B\f(bn)) = 0.
7.3 The hydrodynamic modes 159
With the help of these eigenfunctions, we make the following spectral decom-
position of the Boltzmann propagator (7.13):
With the notation introduced above, we rewrite the spectral function F(k, co)
as:
(fl>), (7-31)
At this point we can verify that the first swm rw/e is satisfied, that is the
integrated intensity of S(k,co) (at fixed k-value) yields the static structure factor
iS(k) (7.7). This is easily seen as follows: we evaluate the integral (by contour
integration):
and we use the relations (7.7), (7.24) and (7.25) to obtain (p\p) = pS(k) and
hence:
Each mode contributes a spectral line in (7.34), and their total weight factors
add up to unity.4
The dominant contributions to F(k, co) come from small values of the de-
nominators in (7.32) which are obtained for the eigenvalues corresponding to
the hydrodynamic modes since in the hydrodynamic regime these eigenvalues
z^(k) are either (9(k) or (9(k2) and k(= |k|) is small. So the hydrodynamic modes
with z^(fc) —• 0 as k —> 0 are slow modes (that is they exhibit slow decay in
space and time). In the long wavelength limit, there are also kinetic modes,
with Rez /i (0) < 0, which give terms with fast exponential decay in (7.30) and
4
At finite |k|-values, 'weight factors' are not necessarily positive, because the matrices are
complex, but not Hermitian.
160 7 Hydrodynamic fluctuations
In the long wavelength limit the hydrodynamic modes ty^(k) and eigenvalues
z^(k) can be determined by a Taylor series expansion:
<°>> = 0
f ^ ^ f i (7.37)
with c\ = k • c, where k is the unit vector along the k-direction. The general
solution of the zeroth order equation is an arbitrary linear combination of the
collisional invariants (7.24):
(7.38)
m\ci+zM\an)Bn=O9 (7.40)
s p I _L
s (s\s) 0 0 0
P 0 z^(p\p)
"n \r\r/ (pic?)
\F\^I / 0
" (1 A\\
0 (p\cj) zM(ci\ci) 0
± 0 0 0 zi
where we have used (7.39) and the relations (p\cf) = (p\p) and (c/|c/) =
The eigenvalues z^ to first order, the hydrodynamic modes xp^ to zeroth order,
and the currents j ^ = (c/ -f z^)\p^ are then given respectively by:
with TXX = \{c2x — c\) and a = + ; here s labels the entropy (or heat) mode, _L
the shear (or transverse momentum) mode, and a the two sound modes.
The solution to the second equation of (7.37) is:
(7.44)
The first term on the right hand side belongs to the orthogonal complement
of the null space of H. The coefficients B^ remain undetermined. By inserting
(7.44) into the third equation of (7.37) and multiplying it on the left by (t/^l
(with X = s,_L and <r), we obtain:
0 =
162 7 Hydrodynamic fluctuations
For A = fi9 the first term on the r.h.s. vanishes and we obtain the eigenvalues to
second order:
1
z(2) ^ v=-(T I - + -IT )(c \C r
Q 2
Here D T is the thermal diffusivity, v the kinematic viscosity, and F the sound
damping coefficient. To further simplify the expression for F, we use the relation:
which follows from (7.39), (7.42) and (7.43). For lattice gas models with sufficient
symmetry (like the triangular lattice with hexagonal symmetry, see Chapter 3),
a fourth rank tensor is isotropic, yielding the equalities:
we have:
i i (7.52)
B,,_, = (v-r)/[2(7C 8 ]
where the asterisk on the summation sign indicates restriction to the hydrody-
namic modes. The coefficients JV\ in (7.32) are evaluated in Section 7.3.2 for
smallfc,and yield:
Combination of (7.31), (7.55) and (7.56) yields the dynamic structure factor
in the Landau-Placzek approximation:
S(k,co) (l^±\ 2DTk2 ly> Tk2
S(k) ~ \ y ) w + (DTk2)2
2
y^(co±
^ csk)2 + {Tk2)
+ 1 [F + (y - 1)DT] *
y ccs *-f (at
(7.59)
the same result was obtained for the lattice gas in Section 7.3.I.7
In order to analyze the dynamic structure factor S(k,co) over the full spectral
range of wavelengths and thereby to assess the limits of validity of the Landau-
Placzek theory, we must know both the eigenvalues, zAt(k) = Rez/z(k) + zlmz/x(k),
and the eigenfunctions of the propagator (7.13). At fixed wavelength k, each
eigenmode contributes, according to (7.31), a spectral line with a maximum
located approximately at Imz^(k) (^ 0 if the mode is propagating) and a width
7
Notice that the small-fc limit and the co-integration cannot be interchanged, because the
co-integral does not converge uniformly near co = 0.
7.5 The eigenvalue spectrum 165
2.0
7.5 The eigenvalue spectrum 167
the transport coefficients k-dependent We see indeed, in Figure 7.1, that the
discrepancy between the macroscopic relaxation rates (dashed lines) and the
exact ones (full lines) becomes larger and larger as k increases, and that the
speed of sound cs(k) exhibits dispersion. The exact diffusivities, as we shall see
in Section 7.6, can be described in terms of generalized k-dependent transport
coefficients: D^k) = Rez^(k)//c 2 = (v(k),r(k),D T (k)}.
However, we observe that, in the range 0.6 <, k <, 1.25, there still is a clear
separation between hydrodynamic and kinetic real eigenvalues; so the contri-
butions of the fast modes may be ignored in the analysis of the hydrodynamic
spectrum. This defines the domain of generalized hydrodynamics where the gen-
eral structure of the hydrodynamic spectrum should be preserved, but where
important discrepancies will be observed between the Landau-Placzek theory
and the Boltzmann theory for S(k,co). Closer inspection of Figure 7.1 leads to
the following observation: for k parallel to a basis vector, say lx, the matrix
in (7.31) is invariant under the reflection cy <-• — c y. Consequently the matrix
in (7.30) can be decomposed in two subspaces with even and odd parity in cy.
Therefore S(k,co) in (7.31) only couples to the even subspace containing sound
fa = a = +) and heat (// = s) modes. It does not couple to the shear mode
(fi =_L) with odd Cy-parity. Therefore, when considering S(k,co), the domain of
generalized hydrodynamics extends to about k ~ 1.5, which corresponds to a
wavelength of four lattice units.
In Section 7.6, we shall test these concepts of generalized hydrodynamics by
substituting the k-dependent parameters c s(k), F(k) and Dj(k) (as obtained from
the numerical evaluation of the eigenvalues) into the Landau-Placzek expression
(7.58) and compare the result with the Boltzmann prediction (7.31) which does
not involve any long wavelength or slow mode approximation.
(Das, Bussemaker and Ernst, 1993), where the collision matrix in (7.31) is only
a small perturbation to the propagation term which constitutes the dominant
contribution to the spectral density (see Equation (7.19)).
high density is shown in Figure 7.2. According to the discussion in Section 7.5.1
(see also Figure 7.1), this spectrum is measured at a fc-value well inside the
hydrodynamic regime where the predictions of the Boltzmann theory (7.18) and
of the Landau-Placzek theory (7.58) coincide, and are in excellent agreement
with the 'experimental' spectrum as Figure 7.2 shows. The spectrum consists of
a doublet of spectral lines located at co = ±csk, corresponding to the acoustic
modes, and a peak centered at co = 0, which arises from fluctuations at constant
pressure and corresponds to the thermal diffusivity mode (for non-thermal LGA
models, this central peak is absent). Thus we observe that in the small fe-domain
the spectral density of lattice gas fluctuations is very well reproduced by the
Landau-Placzek theory.8 Most important is the observation that the dynamical
structure factor of the lattice gas exhibits the typical structure and line-shapes
of the Rayleigh-Brillouin spectrum measured in real fluids (Boon and Yip,
1980); it was precisely to explain how spontaneous hydrodynamic fluctuations
in continuous fluids give rise to this type of spectrum that the Landau-Placzek
-2.0
Figure 7.2 The LGA dynamic structure factor in the hydrodynamic regime.
S(k,co) as a function of the frequency co for the GBL lattice gas at high
density and low k. Comparison between the simulation results (solid line) and
the theoretical predictions: Landau-Placzek theory, Equation (7.58), (dashed
curve), and Boltzmann theory, Equation (7.18), (dotted curve).
Thermodynamic state: p = 6.0; e = 6.7. Automaton universe: L x L — (512)2.
k = 20/co ^ 0.245 reciprocal lattice units, co is expressed in units of 2n/T0,
where To is the total number of time steps, and S(k,a>)/S(k) is expressed in
reciprocal co units.
The slight discrepancy in the height of the doublet peaks (see Figure 7.2) is probably
due to the weaker resolution at small /c-values of the peaks where fluctuations in the
simulation data are most pronounced.
170 7 Hydrodynamic fluctuations
theory was initially developed. The comparative inspection of Figure 7.2 and
Figure 7.3 is instructive as the latter shows the power spectrum measured by
light scattering spectroscopy in a real fluid at k ~ 2 x 105 cm" 1 , i.e. a wavelength
of a few thousands angstroms.
When the value of the wavenumber increases, the Landau-Placzek theory
progressively fails to reproduce the measured spectra, because, as we discussed
in Section 7.5.2, linearized hydrodynamics with constant transport coefficients
becomes invalid at short wavelengths. For the GBL lattice gas, at k = 5Ofco — 0.6
one enters the regime of generalized hydrodynamics where the Boltzmann theory
gives very good account of the simulation data.
When k increases further, the deviations from classical hydrodynamics become
larger, as exemplified in Figure 7.4, which shows a comparison between the
measured and computed dynamical structure factors at k = 136fco — 1-67 for
the GBL model. The figure shows that the Boltzmann prediction remains valid
down to quite short wavelengths (here A ~ 4 lattice units). Here we are at
the border between generalized hydrodynamics and the kinetic regime (see
Figure 7.1): the Landau-Placzek theory with constant transport coefficients is
expected to fail completely, but when combined with the fe-dependent sound
speed and fc-dependent transport coefficients, as computed from the eigenvalue
spectrum, considerable improvement is obtained, as illustrated in Figure 7.4.
Although it is clear that hydrodynamics becomes invalid in this high k domain,
it is remarkable that the general structure of the hydrodynamic power spectrum
persists down to such short wavelengths. Similar observation can be made from
the spectra measured by neutron scattering spectroscopy in actual fluids, e.g.
liquid neon at k ~ 107 cm" 1 , i.e. a wavelength shorter than 100 angstroms (Bell
et al, 1975).
!- 5 |
1 1 1
' 1I ' I ' I ' I ' iI ' I1 '' I ' | Figure 7.3 The power
spectrum of liquid argon
in its normal liquid range
measured by light
I Q \-- II -I
- scattering spectroscopy
(Fleury and Boon, 1969).
3' [ ^ 22.883
- GHz
883 GHz 1
JjJ 22.883
- 883 GHz
GHz ^ o) is given in GHz and
S{k,co) in arbitrary units.
'}
0.5 -
\H
0.0
-4.0 -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0
co (GHz)
7.6 Power spectrum 171
-0.4
172 7 Hydrodynamic fluctuations
see Figure 7.6. In this almost collisionless regime, correlation effects through
the periodic boundary conditions might be quite substantial. Nevertheless the
Boltzmann theory yields predictions which follow more or less the general trends
of the spectral features. In fact, from the eigenvalue spectrum and a numerical
analysis of the eigenfunctions, one can see that 5(k,co) is strongly coupled to
six modes with Rez^(k) ~ 2, four of which are propagating, which explains
qualitatively why the 'Brillouin lines' in Figure 7.6 have a complex structure
resulting from the superposition of several modes.
4.U 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Figure 7.6 The dynamic
structure factor of the
GBL lattice gas at low
'.'/// \\\\;;; --
3.0 ;' :A J-i It \ - density in the kinetic
regime. Same conditions
and specifications as for
Figure 7.5 except
^ 2.0 -
k = 52/co — 0.64 reciprocal
\
lattice units.
1.0
-3.0
; J. . .k J -2.0 -1.0 0.0 1.0 2.0 3.0
0)
9
At very high frequencies, the transport functions can also become co-dependent as the
system exhibits viscoelastic behavior when the response becomes non-local in time.
7.6 Power spectrum 173
(open circles) is compared with Imz+(/c) = ±cs(k)k (solid line) computed from
the lattice Boltzmann equation; the classical behavior is shown by the dotted
line, the thermodynamic speed of sound cs being calculated from (7.39). The
measured dispersion in the sound speed is seen to be remarkably well accounted
for (up to k ~ 1.5) by the Boltzmann theory for cs(k) confirming the validity of
generalized hydrodynamics.
From the expression for the classical hydrodynamic spectrum (7.58), we can
compute the ratio of the integrated intensity of the central Rayleigh peak to
those of the Brillouin peaks, and obtain the Landau-Placzek ratio, which is
a thermodynamic quantity equal to (y — 1). For the lattice gas, the measured
value yexp — 1-34 is in very good agreement with the theoretical value 1.32,
calculated from (7.51): y = c\/c\. In the regime of generalized hydrodynamics,
the Landau-Placzek ratio then becomes fc-dependent, but it is very hard to mea-
sure quantitatively the dispersion effects on y, and thereby make a comparison
with theoretical predictions, because of the difficulty of evaluating the integrated
intensities of the spectral peaks in a quantitative manner at high /c-values (as it
is in real fluids).
In the long wavelength limit, the value of the transport coefficients can be
obtained from the line-widths of the spectral lines of S(k,co) since the Landau-
Placzek theory predicts three Lorentzians of width DTk2 (for the central Rayleigh
peak) and Tk2 (for the two shifted Brillouin peaks). In Figure 7.8, the simulation
data for Rayleigh width (black dots) and for the Brillouin width (open circles)
are shown as functions of k2. The dashed line represents the classical values
Djk2 and F/c2 and corresponds in fact to two coinciding lines because here
(by mere chance) D T ^ F whose values are given by the Boltzmann transport
coefficients in the limit k —• 0, derived in (7.46). The solid curves D T(k)k2 and
T(k)k2 computed from the Boltzmann eigenvalues (see Section 7.5) indicate that
high
density (°pen
circles), classical
hydrodynamic value
cs x k (dotted line), and
Boltzmann theory
prediction (solid curve).
Same units as in
Figure 7.5.
174 7 Hydrodynamic fluctuations
the fluid, the concentration or species density, and which will couple to the other
thermodynamic and mechanical variables.
The FHP model has been generalized to study diffusive phenomena in binary
fluids using 'macroscopic' experiments (Bernardin and Sero-Guillaume, 1990;
McNamara, 1990; Noullez, 1990). Typically the measurement would be per-
formed in a lattice gas composed of red and blue particles (where the color
is a passive property used to distinguish species which otherwise do not differ
from one another) starting from an initial condition with a step function color
density profile. The observer would then measure the evolution of the density
profile, and evaluate the diffusion coefficient by fitting the 'experimental' profile,
obtained after some time, to the solution of the diffusion equation subject to
the appropriate boundary conditions (McNamara, 1990; Noullez, 1990).
Here we shall consider a more microscopic measure along the lines of the
spectral analysis developed in this chapter. But the mixture of two real fluids
exhibits a power spectrum where the central peak is not a simple Lorentzian,
even in the long wavelength limit (Boon and Yip, 1980): it has a spectral structure
where it is difficult to separate the contributions from entropy fluctuations and
from concentration fluctuations which, in general, are not decoupled. 12 In this
respect, the lattice gas offers an interesting approach to probe diffusion in a
simple fluid mixture, because we can consider a non-thermal LGA and define an
observable - the weighted difference of the species densities - whose fluctuation
correlations yield the diffusive mode independently of the other modes, so that
the corresponding power spectrum provides a measure of diffusion dynamics
solely (Hanon and Boon, 1997).
with d the average density per channel, and Xr the concentration of red particles.
Note that:
b
1 1
fo M(1 Zr) =
^ + ^rCd + ^1UC = p (7 62)
'
is given by:
oo
where (5n,-(k, £*) is the spatial Fourier transform of 5n,(r*, t*) and Kj = /j e q '(l —
/j e q ) ). Using (7.66), we rewrite (7.64) as:
oo b/2 b/2
r y (M,)> 9 , (7.67)
and the static structure factor (the Fourier transform of the equal-time van Hove
13
Remember that V is also interpreted as the volume of the lattice i f universe, and here
the lattice is finite and has periodic boundary conditions.
7.7 Diffusion and correlations 177
or
This expression for the dynamic structure factor is similar to the result given
in (7.17) and (7.18), and is exact within the Boltzmann approximation, but
again the explicit analytical inversion of the b x b matrix in (7.72) cannot
be performed in all generality. However, perturbation methods can be used
to compute analytically S red(k,co) in the hydrodynamic limit: |k| — • 0 and
co ~ $(|k|), $(|k| 2 ); beyond the long-wavelength, long-time domain, one has
recourse to numerical evaluation of (7.72) to compute the Boltzmann power
spectrum. Since the developments proceed essentially as for the single component
lattice gas, we shall not present the detailed analysis 14 and will focus on the
main results.
factor of the 'colored' lattice gas which reads (Hanon and Boon, 1997) :15
2 KT + Kb co2 + (Dk2)2
b K2 ^ / r/c 2 r/c cs/c ± i
2 * r + Kb ^ V (« ± c*k) 2 + (r/c 2 ) 2 cs (co + csk)2
(7.73)
Here ?cr (red) and Kb (blue) are given by:
Kr = dXr(l -dXr), Kb = d(l- Xr )[l - d(l - &)] , (7.74)
cs is the speed of sound (here cs = y 3/7), F is the sound damping coefficient
(see (7.52), but without the second term since the colored LGA is non-thermal),
and D is the diffusion coefficient. We recognize in (7.73), the typical structure
of the hydrodynamic structure factor (compare with (7.58)): at fixed value of
k, the spectrum consists of a Brillouin-doublet centered around +/ccs, and of
a central peak characterizing color diffusion. It is straightforward to verify
that for Xr = 1, (7.73) reduces to the usual Landau-Placzek expression of the
hydrodynamic spectrum for the non-thermal single-component fluid.
Now if we consider the fluctuations of the observable defined by:
d = p/pmax, the reduced density (which is also the average density per channel).
Accordingly the hydrodynamic regime is defined by fc/f <C 1, the generalized
hydrodynamic regime (Boltzmann regime) by k/f < 1, and the kinetic regime
by faff > 1.
We shall discuss the power spectra obtained from automaton simulation
data and will compare the results to the analytical Landau-Placzek expressions
and to the predictions of the lattice Boltzmann theory. For the latter, one
uses the eigenvalue spectrum of the propagator T which can be evaluated
numerically over the complete fc-domain, so extending the computation of
the power spectrum to the region of k values where analytical evaluation
can no longer be performed. Figure 7.9 shows a typical eigenvalue spectrum
as computed numerically, where the 14 modes of the CFHP model can be
distinguished.
We observe that in the range k < 0.4, corresponding to wavelengths X > 15/o
(with A) the lattice unit length), the four slow modes are well separated from
2.0
180 7 Hydrodynamic fluctuations
the kinetic modes, and their behavior is correctly given by the hydrodynamic
expressions. We therefore expect that in this domain, the Landau-Placzek theory
should provide a correct description of the dynamic structure factor. This
is indeed confirmed by the comparison between the results obtained from
simulation data and theoretical predictions as shown in Figure 7.10, where
we also notice that the hydrodynamic spectrum is indistinguishable from the
Boltzmann spectrum.
From the hydrodynamic spectrum one can directly measure the diffusion
coefficient D(d,xr); the expression (7.76) of Sdiff(k,co) is a single Lorentzian
with a half-width Aco = Dk2, so that plotting Aco as a function of fe2, as in
Figure 7.1 l(a), a linear least-squares fit yields a slope whose value gives a direct
measure of the diffusion coefficient. This measure is in agreement with the
predictions of the lattice Boltzmann theory up to reduced densities d ~ 0.25
as Figure 7.11(b) illustrates. Deviations for larger values of the density are not
surprising since the Boltzmann theory does not take into account correlations
which are known to become important at moderate and high densities.
As k increases from 0.4 to 1.4, there is still a distinct scale separation between
slow and fast modes (see Figure 7.9(a)), but the eigenvalues of the slow modes
depart significantly from the hydrodynamic prediction, indicating the breakdown
of the local response hypothesis: the transport coefficients become /c-dependent.
As a result, the hydrodynamic theory no longer describes the power spectrum
correctly,16 but the complete Boltzmann spectrum is in rather good agreement
with the simulation results, as Figure 7.12(a) shows. At higher densities there
is a discrepancy between the Boltzmann spectral density and the experimental
power spectrum (see Figure 7.12(b)) which reflects the failure of the Boltzmann
3.0
2.5
2 0
1.5
1.0
0.5
0.0
0.0 0.2 0.4 0.6 0.8 1.0
d
16
Note that even at rather short wavelength (A ~ lO^o) the experimental data can still be
reasonably fit with a Landau-Placzek spectral function if the transport coefficients are
parameterized, indicating that a hydrodynamic type description holds qualitatively down
to quite short wavelengths, provided the transport coefficients are renormalized.
182 7 Hydrodynamic fluctuations
theory to evaluate correctly the transport coefficients in the high density region
where correlations are important.
With this discussion we close our investigation of the hydrodynamic fluctu-
ations in the lattice gas (with one and two species). We shall return briefly to
the two-component case in Chapter 10 when we discuss two-fluid simulations.
But we now turn to the macroscopic dynamics of the LGA; after showing in
this chapter that the lattice gas exhibits spontaneous fluctuations which are
compatible with those observed in real fluids, we go back to the important
question: are the equations describing the macroscopic behavior of the lattice
gas compatible with the Navier-Stokes equations of classical fluid dynamics?
This question was addressed in Chapter 5 and will now be revisited with a
different complementary approach.
0.05
Chapter 8
1
Most of the derivations of the LGA hydrodynamic equations given in the literature
(Frisch et al, 1986, 1987; Ernst, 1990a; Rothman and Zaleski, 1994) include
non-linearities in the Euler part of the equations, but the dissipative term is calculated
by linearizing around global equilibrium, thus limiting the validity of the results to a
regime in which linear response and simple fluctuation-dissipation relations hold. In the
development given by Wolfram (1986), the first few non-linear corrections are included,
but this derivation remains at a phenomenological level, as the non-linear coefficients
are assumed to be correctly obtained by direct Taylor expansion.
183
184 8 Macrodynamics: projectors approach
can be used to analyze phenomena taking place in systems arbitrarily far from
equilibrium, for instance in thermal LGAs under large temperature gradient.
In order to derive the hydrodynamic equations, we make use of the Boltzmann
hypothesis (see Section 4.4.2) that particles entering a collision are uncorrelated.
The average outcome of the collision is assumed to be a non-equilibrium
distribution which is close to the local equilibrium one. A local equilibrium
distribution has the same analytical form as an equilibrium one, but with space-
and time-dependent thermodynamic potentials (temperature, chemical potential,
etc.). The non-equilibrium configuration can then be split into a local equilibrium
component, which depends only on the local values of the set of densities of
conserved quantities of the system, and a non-local equilibrium contribution, a
'projection' onto the space orthogonal to the subspace spanned by the constants
of motion. This procedure is intimately related to the projection operator
methods developed to derive the hydrodynamic equations from the microscopic
dynamics in continuous fluids (Zwanzig, 1960; Mori, 1960; Procaccia, Ronis
and Oppenheim, 1979).2 We shall effectively use a simple projector that neglects
the coupling of the linear hydrodynamic variables (mass, momentum, energy
densities) to bilinear and higher order modes. With the treatment used here -
and perhaps at the cost of a somewhat formal analysis - we shall obtain a
general formulation of the lattice gas hydrodynamics valid for non-equilibrium
states arbitrarily far from global homogeneous equilibrium. We shall also make
contact with the linearized equations and show that the dissipative coefficients
are given in terms of time-correlation functions in the form of the Green-
Kubo expression of the transport coefficients. The present analysis is within the
logical continuation of Chapters 6 and 7; in Chapter 5 we give an alternative,
less formal, derivation of the lattice gas hydrodynamic equations using the
Chapman-Enskog method. Notice that, whatever the method - be it for discrete
systems or continuous fluids - there is always a common, unavoidable major
step in the analysis: the multiple scales hypothesis, by which one anticipates that
there is a natural space- and time-separation between hydrodynamic phenomena
and microscopic processes.
8.1 Preliminaries
The observables in the automaton are defined as averages over the stochastic
dynamics and over an ensemble of initial conditions. The averaging is denoted by
angular brackets, (.. .) N E , where the subindex stresses the fact that the ensemble
is a non-equilibrium one. The collection of local quantities preserved under
collision are labeled by the Greek index cp running from cp = 1 to cp = iVhydro,
where A/hydro denotes the number of hydrodynamic variables. This set should
contain the quantities conserved in the dynamics of the physical system modeled
by the automaton. The index cp = 1 corresponds to the number of particles,
which is conserved by the collisions.
Before the collision step, the configuration of a node r* e <£ is given by :4
and:
The last equality follows from the fact that the value of the hydrodynamic
variables is preserved by the collision step. The hydrodynamic variables include
some or all of the following: particle number, momentum, and, in thermal
automata, kinetic energy per node. For example, in non-thermal automata whose
collision rules preserve the number of particles and the momentum per node,
the hydrodynamic variables are the number density, p(r, t), and the momentum
density, pu(r, t), which can be collected in a two-dimensional vector:
-' = ( c j • (8J)
where £{s}{s'} is equal to 1 when, starting from a configuration {s}, the collision
produces configuration {s'}, and is 0 otherwise. For stochastic automata, ^ is a
random Boolean variable. In this case the average post-collisional configuration,
starting from a given pre-collisional configuration {ft(r, £)} is:
JiA(s^)5({n(T9t)},{s}), (8.10)
where A(s^sf) is the matrix whose entries are the probabilities of having a
configuration {s'} as the outcome of a collision starting from a configuration
{s}. The quantity S ({n}, {5}) is equal to one if the configurations {s} and
{ft} are equal, and zero if they are different. A representation of this Kronecker
delta in configuration space is (see Section 2.2.4):
nskk(l-nkf-^. (8.11)
k=l
an arbitrary distribution of the Fermi-Dirac form for the single particle distri-
bution (see Section 4.4):
ft = T - \. Al , (8-13)
1 + e x p { 0 * A t\
without correlations between fluctuations in different channels, is invariant under
the collision step. The quantities (j) are the thermodynamic fields conjugate to
5
We ignore effects associated with the presence of spurious global invariants (see Brito,
Ernst and Kirkpatrick, 1991).
8.1 Preliminaries 187
<p=l
(4(r,0)LE = 7^ fjLE(l
1
A
-'
= <A(r,t))m. (8.16)
Equation (8.16) defines the thermodynamic potentials implicitly in terms of the
local values of the densities of conserved quantities: (j)(r,t) = <>
/ ((A(r, t))). For
conserved quantities, since the local equilibrium ensemble has been defined in
such a way that (A(r, t))LE = {A(r9 f))NE, we may omit the label indicating the
ensemble over which the average is taken.
The lattice Boltzmann equation (see Section 4.2.2):
constitutes the starting point for the derivation of the macroscopic equations that
describe the evolution of the hydrodynamic fields. Multiplying Equation (8.17)
by At, summing over the index i and using the displacement operator defined
by:
/,(r + c , t) = exp {c; • ^ } /,(r, t), (8.18)
where the scalar product denoted by • is a contraction of the spatial indices, we
obtain:
b
ex
U(r, r + 1)> = Yl P {~c< • Vf} ffa ^ ^f (8.19)
1=1
188 8 Macrodynamics: projectors approach
The microscopic evolution equations (8.21) are finite difference equations con-
taining the full dynamics of the automaton able to describe its behavior at all
scales. The microscopic scales are given by the automaton lattice spacing and, for
time, by the lapse between two consecutive propagation steps. There is another
characteristic microscopic time TLE, of the order of a few collisions, which corre-
sponds to the relaxation to local equilibrium. In order to derive hydrodynamic-
type equations from the evolution equations (8.21), we introduce the parameter
e, which represents the ratio between characteristic microscopic and hydrody-
namic length scales. This ratio is assumed to be small, a fact which allows us to
carry out a multiple scale analysis (see also Chapter 5). Thus, we anticipate the
scales of physical interest in the solution of (8.21), and introduce the following
variables: a space variable ri = er, and two time variables, t\ = e t, which is
of order e° in theregime where Euler equations are valid, and ti = e2t, which
is of order e° in thedissipative regime. We further assume that deviations from
local equilibrium at each node are small. This corresponds to the Chapman-
Enskog picture of a non-equilibrium state (Chapman and Cowling, 1970), where
the deviations from local equilibrium are on the order of the gradients of the
conserved quantities, i.e. of order e. The substitutions (introduced in Chapter 5):
= 2 ^ 1 :(M(M))LE)
-c f • Vr } - 1) At SfKr, t\ (8.24)
i=l
8.2 Multiple scales analysis 189
with:
It = CiAi, (8.26)
b b
and:
b b
Equations (8.23) are the Euler equations describing the streaming part of the
dynamics, and Equations (8.24) contain the dissipative contributions.
We note that we can rewrite the second term on the l.h.s. of Equation (8.24)
by making use of the first order equation (8.23):
and:
with:
We have assumed that 5fu and consequently Sf'i9 is small (of order e); this
assumption will be discussed below.
The matrix ifLE(r, t) has the form (see Appendix A.12):
[
^ r , r) [(5/;(r - c;-, t-l)
7=1
(8.42)
where e~dt is the time-displacement operator. We rewrite Equation (8.42) as:
b b
E
Sffat) = Yl ^if(rj) f) (r,t) + ^ ^ ( r , t ) 3 / ; ( r , t - l ) , (8.43)
7=1 7=1
where the matrix operators ^~LE and y L E are given by the expressions:
<n/(M) = ^ / ( r , 0 exp{-c r V r }, (8.44)
5^(r,t) = ^ ( r , t ) (exp{-c r V r }exp(-^) - l ) . (8.45)
Equation (8.43) is iterated to yield:
Sf'i(T,t) =
E E ( l l ^ LE ( r ^-^)l ^5f(r,t-T)/P(r,t-T)
y,/ T=0 I T'=0 J fy
b ( t-\ \
The requirement that 5/-(r, t) be of order e, and the vanishing of the 'random
force' term depending on the initial conditions, <5/j(r, 0), are fulfilled by choosing
an initial local-equilibrium ensemble in which the values of the thermodynamic
variables vary significantly only over a hydrodynamic length scale. In the case
192 8 Macrodynamics: projectors approach
that this were not so, there should be an initial transient regime where the
hydrodynamic equations are not valid. The duration of this regime is on the
order of TLE, the time necessary to relax to local equilibrium.
We notice that:
T>TLE, (8.48)
which is a consequence of the fact that all the eigenvalues of J^ E (r, t) are strictly
smaller than one in absolute value. So, given that the kernel in Equation (8.46)
is non-zero only on a microscopic time scale TLE ~ (9(e°), and that during this
time the local equilibrium distribution does not change significantly (since it
depends only on slow variables, whose variation is on a time scale (9{e~1)),
Equation (8.46) can be approximated by:
t b
Within the same level of approximation we can further replace the upper limit
in the summation over T by infinity. Hence, the last term in (8.24) becomes:
4 £ {[ dh) ff(t,t)
1
[
oo
£ c-4- £{[
ij TT=1
=1
3
ffE(r,t) = Knr,t)A^ (A\A)-X {r,t), (8.51)
d(A(r,t))
8.4 Linear response and Green-Kubo coefficients 193
proved in Appendix A.12 (see Equation (A.58) with n = 1). Notice that the
quantity on the r.h.s. of (8.50):
dt(A(r,t)) + V r -(J(r,0))LEL
(8.53)
The prime in the summation indicates that the first term of the sum (T = 0) is
multiplied by 1/2. This factor appears as a consequence of the discreteness of
time in the automaton. We have also replaced the matrix J? LE(r, t) by J^LE(r, £),
since the quantity upon which it acts contains no projections onto the conserved
quantities.
Equation (8.53) is the most general form6 of the hydrodynamic equations in
the lattice Boltzmann approximation; the expression of the dissipative term (on
the r.h.s. of Equation (8.53)) is valid even far from global equilibrium. So these
equations are the fundamental equations describing the macroscopic dynamics
of a lattice gas automaton in the hydrodynamic regime; in particular, they can
be used to explore non-equilibrium phenomena such as far from equilibrium
steady states and hydrodynamic instabilities in lattice gas automata.
Now we shall show that linearization of Equation (8.53) yields the usual linear
response results, including the Green-Kubo formulae for the transport coeffi-
cients. We first recover the equations of linear hydrodynamics (see Chapters 6
and 7) by carrying out an expansion of (8.53) around the global equilibrium
state characterized by the values of the hydrodynamic variables:
6
Notice one restriction: the validity of Equation (8.53) is limited by the fact that we have
neglected the coupling of the evolution of the single-particle distribution to ^-particle
distributions.
194 8 Macrodynamics: projectors approach
Making use of Equation (A.58) (see Section A.ll), we expand the single particle
distribution ffE(r, t) around equilibrium:
xcrVr/PM)
x c r V r (Kj
^ ) ^ 4(r, 0)
In the second equality of Equation (8.58), we have used the definition of the
equilibrium projectors (Zanetti, 1989; Ernst, 1990a; Ernst and Dufty, 1990)
(^ + 1 = 1) acting on an arbitrary vector {Bj}bj=1:
(8.59)
The projected quantities have the following properties (see Appendix A. 13):
(8.61)
(A\2B)(eq) = ( 4 | B ) ( e q ) = 0. (8.62)
8.4 Linear response and Green-Kubo coefficients 195
A | j) = (A | 1 J) (eq) = 0, (8.63)
7=1
we can write:
T=0 ^ 4J
= £' T=0
= £'
T=0 rr'
(8.68)
{SniiTrfSnfaO))^ (8.69)
decays to zero sufficiently fast for the sum in (8.68) to be finite (Schmitz and
Dufty, 1990). The result is that the dissipative coefficients exhibit a Green-Kubo
form (Rivet, 1987b; Ernst, 1990b) analogous to the continuous theory expression
(Kubo, 1958; Zwanzig, 1965).
This concludes the derivation of the non-linear hydrodynamic equations
for lattice gas automata satisfying the semi-detailed balance condition, and of
their subsequent linearization yielding correlation function expressions for the
dissipative coefficients.
(dmir^Srij^O)) (8.70)
the solutions of the two equations, and the result for the normalized velocity
autocorrelation function \p(t) reads (Ernst, 1991):
where bd/vo is the number density per elementary unit volume of the lattice (e.g.
v
o = \/3/2 m * n e triangular lattice), and Ds denotes the self-diffusion coefficient.
This result is, in fact, the same as in continuous fluid theory (Ernst, Hauge
and van Leeuwen, 1971), apart from the factor (1 — d) which comes from the
exclusion principle.
10°
101 r
io- 2 -N
101-3 =-
101-5 -
101-6 =-
101-7 -
10*
10° 101 102 103 104
t
Hydrodynamic regimes
= 0, (9.1)
The macroscopic variables are the ensemble-averages of the conserved quantities (see
Section 4.5).
199
200 9 Hydrodynamic regimes
and
* + dp [gjaupj + daP =
^ +
pip j 2) = — p — — ^—j 2
(9.4)
with
The geometrical coefficients ^2 and £4 are defined in Section 4.5.2, and their
values are given for various models in Table 4.1.
The expansion parameter e in (9.2) is the space-scale separation parameter;
the hydrodynamic limit corresponds to e < 1, the condition under which the
macrodynamical equations were derived. The expansion parameter rj represents
the order of magnitude of the macroscopic momentum flux j , that is, the (small)
'deviation' from zero-velocity basic equilibrium (see Section 4.5).
where po is the constant and uniform mean density, and the velocity fluctuation
9.1 The acoustic limit 201
u(t, r) (with zero mean velocity). When p(t, r) and u(f, r) are small, we can ignore
the non-linear terms (in dp and u) in Equations (9.1) and (9.2) to obtain:
dp
6 3 D-
+ -yda— = vdpd pua + (V + —
where v and v' stand implicitly for v(po) and v'(po) respectively. Using the
standard 'nabla' notation, the preceding set of equations can be rewritten as:
(9.7)
_
^ V ( ^ ) = vAu + (V + ^-jr^v) V(V • u),
b po \ D /
where we clearly recognize the standard equations of viscously damped acoustics.
Here a comment is in order about Galilean invariance. The acoustic equa-
tions (9.6) or (9.7) are Galilean invariant while the full lattice gas macrodynam-
ical equations are not (because of the non-Galilean factor g(p)). In the acoustic
limit, Galilean invariance is recovered because the non-linear term dp(g(p)upj a)
is absent in (9.6) and (9.7).
Proceeding one step further in simplification, we consider the inviscid flow,
i.e. we consider short (yet macroscopic) time scales where viscous damping is
negligible. The dynamics is then governed by the standard Eulerian acoustic
equations:
Po
. , (9-8)
This set of equations has sound wave solutions 2 of the following form:
(9.10)
2
These sound waves are called 'longitudinal waves' because the fluid velocity u is
everywhere parallel to the wave vector.
202 9 Hydrodynamic regimes
Ma=-. (9.11)
From (9.9), we observe that the relative density fluctuation Sp/po is of the same
order of magnitude as the Mach number.
On longer time scales, the viscous terms become relevant, and Equations (9.7)
exhibit solutions describing sound wave propagation with an amplitude decaying
exponentially slowly in time:
D 2 '
and the phase shift <f> is given by:
ock
sin© = —,
Cs
Po
These are not propagating waves (there is no phase rotation); their amplitude
decays exponentially over long time scales. The shear wave solution is not
of direct interest for lattice gas simulations of fluid dynamics, but appears to
be very useful to obtain an independent measure of the kinematic viscosity
coefficient v. From the damping coefficient a of sound waves and the damping
coefficient v of shear waves, we can then compute the shear and bulk viscosities
separately (see Chapter 10).
9.2 The incompressible limit 203
dtp + u • V(p) = 0.
The continuity equation and the momentum equation for real fluids then be-
come:
V-u = 0,
p (9-14)
dtu + u V(u) = - V ( —) + vAu.
Po
The first equation is the 'incompressibility condition', and the second is the
'incompressible Navier-Stokes' equation (see for example Batchelor (1967),
page 174).
For 'real' fluids, the incompressible Navier-Stokes equation can be seen as
a low Mach number approximation of the fully compressible Navier-Stokes
equations. 3 The main feature of incompressible regimes is that the density can
be taken as approximatively equal to a constant po, except in the pressure term
which must fluctuate so that the two equations in (9.14) be compatible. This
amounts to saying that pressure fluctuations are slaved to the fluid velocity field
so that the divergence of u remains zero.
The same kind of limit can be taken for the lattice gas; the macrodynamical
equations (9.1) and (9.2) then reduce to:
dpup = 0,
Stua + g(po)dp(uaup) H d aP =
v J
Po
or, equivalently:
V • u = 0,
p (9.15)
dtu + g(po)u • V(u) = - V ( — ) + v(p o)Au.
Po
These equations are not identical to the continuous fluid incompressible Navier-
Stokes equations (9.15) because of the non-Galilean factor g(po) which is con-
stant but (in general) not equal to one. We now show that there exists a simple
rescaling by which we recover the standard incompressible Navier-Stokes equa-
tions.
3
There are some tricky points in the Mach number expansion (see e.g. Majda, 1984).
204 9 Hydrodynamic regimes
Let us denote by UQ and by ^o the typical fluid velocity and the typical space
scale of the problem (expressed in lattice units) respectively. For example, con-
sider the flow around a fixed solid obstacle in a wind tunnel-like configuration:
the typical velocity wo would be the up-stream mean velocity in the simulated
wind tunnel, and the typical length £§ would be some transverse size of the
obstacle. With a length scale and a velocity scale, we can build the typical time
scale:
to =
g(po)uo
This time scale is equivalent to the 'circulation time' commonly used in con-
tinuous fluid dynamics. The presence of the non-Galilean factor g(po) in this
definition is necessary to render the rescaled equations of motion Galilean-
invariant.
With space, time and velocity scales, we can define the following dimensionless
variables:
V • U = 0,
(9.17)
dTU + U • V(U) = - V ( n ) + —AU,
Re
where the operators V and A are taken with respect to the rescaled space variable
R. The parameter:
U
Re = ^ ^ . (9.18)
v(po)
will be called the 'Reynolds number', because it plays, in lattice gases, the same
role as the Reynolds number for real fluids.
The equations (9.17) are identical to the dimensionless form of the 'real
world' incompressible Navier-Stokes equations (9.14). They involve only one
dimensionless control parameter: the Reynolds number. The rescaling (9.16) is
sufficient to transform the incompressible macrodynamical equations of non-
thermal lattice gases into the classical incompressible Navier-Stokes equations
for real Newtonian fluids, which are Galilean-invariant.
9.3 Comments 205
9.3 Comments
dtp + V • u = 0,
One of our main objectives has been to show that single-species non-thermal
lattice gases can exhibit large-scale collective behavior governed by the same
continuous, isotropic and Galilean-invariant equations as real Newtonian fluids.
This is true despite the intrinsically Boolean, spatially discrete, anisotropic and
non-Galilean invariant structure of lattice gases. Moreover, in the past 10 years,
further lattice gas models have been designed to incorporate more complicated
physical features such as reactive processes, magneto-hydrodynamic phenomena
or surface tension (see Section 11.4 in Chapter 11).
On one hand, there has been considerable effort in basic research to under-
stand the subtleties of the statistical mechanics of lattice gases and on the other
hand intense work has been accomplished to take advantage of the similarities
between lattice gases and real fluids in order to simulate fluid motions with
simple and easily implemented lattice gas algorithms. Indeed, because of their
fully Boolean cellular automaton structure, lattice gases are excellent candidates
for efficient implementations on both dedicated and general purpose computers
with serial, vectorial, parallel or even massively parallel architecture. In addition,
various physical effects can be added at low cost. For example, the presence
in a flow of a rigid fixed obstacle is extremely easy to take into account: it
just requires replacing the standard collision rule by a bounce-back rule (see
Section 2.4.1) on all nodes covered by the obstacle. Modifying the shape or
the position of the obstacle is almost immediate, and no mesh modification is
necessary.
Besides technicalities, the lattice gas methods for fluid dynamics offer an-
other fundamental advantage: the simulation of flows where bifurcations are
207
208 10 Lattice gas simulations
suspected. A bifurcation occurs when a given regime loses its stability in favor
of another regime, as soon as a dimensionless control parameter exceeds a
critical value. The new regime originates from small perturbations which grow
instead of being damped. This growth saturates because of non-linear effects,
and leads to the new regime. These small perturbations operate some kind
of 'natural selection' by eliminating unstable regimes. In natural fluid systems,
thermal noise plays this role. In standard (non-microscopic) computational fluid
dynamics methods, the unstable regime can survive artificially, and an externally
imposed disturbance is sometimes necessary to trigger or at least to fasten the
onset of the new regime. Since the lattice gas has fundamentally a microscopic
structure and possesses spontaneous microscopic fluctuations (see Chapter 7),
the macroscopic variables are naturally noisy: they are averages of naturally
fluctuating microscopic quantities. Thus no external perturbation is needed to
trigger bifurcations. Examples of lattice gas simulations displaying bifurcations
are described in Section 10.5.
Since this chapter is dedicated to the application of the lattice gas method
to simulate fluid dynamics, we start with a rapid overview of lattice gas im-
plementation strategies and of frequently encountered programming difficulties.
We then give some notions on how to perform physically relevant numerical
simulations offluidflows.Lattice gas numerical experiments are then described
to illustrate these concepts.
The possibility of building dedicated machines, that is, electronic devices spe-
cially designed to run lattice gas algorithms was considered as early as 1984.
Several machines were constructed and successfully used, like the RAP family
of D. d'Humieres and A. Cloucqueur (ENS, Paris), and the Cellular Automaton
Machines (CAM) of T. Toffoli and N. Margolus (MIT).
We briefly describe one of the earliest of these dedicated machines: the
RAP-1 prototype (October 1985).1 The acronym 'RAP' stands for 'Reseau
d'Automates Programmables' (Programmable Automata Array). The RAP-1
machine is a desk-top, relatively versatile, dedicated hardware, which can handle
two-dimensional lattice gas models, with up to 16 channels per node. It is not,
strictly speaking, a computer, since it contains no microprocessor. It is essentially
built with:
The RAP-1 machine can update the Boolean configuration of 512 x 256 nodes
with 16 channels each, at a rate of 50 times per second, so that the instantaneous
state of the whole lattice can be displayed at the European standard video
frequency. Figure 10.1 shows the RAP-1 facility, with its external personal
computer and its video display monitor.
The automaton behavior of the RAP-1 machine can be programmed from
the external personal computer, through a set of three tables:
Figure 10.1 The RAP-1 facility. On the left is the external personal computer
required to program and run the RAP machine. On the right is the RAP-1
machine itself with its display monitor.
210 10 Lattice gas simulations
per node, residing on a piece of lattice with Jf nodes. The storage problem is the
following: how to dispatch the b x Jf bits necessary to encode the instantaneous
Boolean configuration into the available memory words, so as to minimize
memory wasting and to optimize the update speed of the automaton. Two
strategies are commonly used; each strategy has advantages and drawbacks,
that we discuss hereafter. To be simple, suppose we want to implement a lattice
gas model such as HPP (see Section 3.1) which requires 4 channels per node,
on a (non-realistic) computer with 4 bits per memory word. The essentials of
the discussion hold for any number of channels per node, and for computers
with any number of bits per word. We denote by Jf the total number of lattice
nodes that the lattice gas simulation code must manage.
The first strategy is the 'channel-wise storage': each 4-bit word stores the same
channel of 4 different nodes. Thus, 4 arrays of memory words are necessary to
store the whole Boolean field, one for each channel, and each array must contain
Jf 14 words. Figure 10.2(a) illustrates this storage mode. If the number of node
Jf is not an integer multiple of the number of bits per memory word, then
the number of words per array must be Jf 14 + 1, if Jf '/4 denotes the integer
division. The last word of each array is then only partially used. This strategy
has minor memory wasting.
The second strategy is the 'node-wise storage', which, in some way, is 'or-
thogonal' to the the channel-wise storage strategy. Each 4-bit word stores the
4 different channels of the same node. For a lattice with Jf nodes, Jf words
"3 "3 %
1 1 1 1
5 6 S
Nodel 11,(1) n 2 d) n 3 d) n 4 d) Nodel 11,(1) n 2 (D n 3 d) n 4 d)
Node 2 n,(2) n 2 (2) n3(2) n4(2) Node 2 n,(2) n2(2) n3(2) n4(2)
Node 3 11,(3) n 2 (3) n3(3) n4(3) Node 3 11,(3) n2(3) n3(3) n4(3)
Node 4 n,(4) n2(4) n3(4) n4(4) Node 4 n,(4) n2(4) n3(4) n4(4)
Figure 10.2 The two main storage modes of a Boolean lattice gas in the memory words
of a computer. The number of bits per memory words is taken to be 4 for simplicity.
The solid rectangular boxes symbolize the 4-bit memory words, (a) Bit-wise storage;
(b) node-wise storage.
212 10 Lattice gas simulations
are necessary to store the whole Boolean field. This storage mode is illustrated
in Figure 10.2(b). If the number of bits per memory word is smaller than the
number of channels per node, then two or more words may be necessary to store
the Boolean state of one single node. If, on the contrary, the number of bits per
memory word is larger than b, then the Boolean state of two or more nodes can
be stored in one single memory word. There may be some waste of memory if
the number of bits per memory word is not close to an integer multiple of b.
However, the wasted bits can be used to store auxiliary information (obstacle,
body force, etc.).
Neither of these two storage modes can be preferred a priori; the right choice
depends on several factors. The channel-wise storage wastes very little memory.
The propagation phase is conceptually simple, since it reduces to systematic
data shifts within each of the b arrays of words. If the collision rule is simple
enough to be coded with a non-prohibitive number of logical operations, then
the collision phase is also very efficient, since the logical operations (AND, OR,
NOT ...) act on all the bits of their data words at the same time, so that a
single 'turn off the loop' performs the collision on 64 nodes on a computer with
64-bit words for example. However, when the collision is too complex to be
efficiently coded with combined logical operations,2 one must resort to a look-
up table strategy which operates as follows: the b bits of a single node must be
picked from b different memory words, and gathered into a single auxiliary word
encoding the pre-collision state of the node under consideration. This word is
then sent as an address to the collision look-up table, which delivers the post-
collision state. Then, the converse operation must be performed: the b bits of
the output word must be dispatched into the b arrays of words. These gathering
and dispatching operations can be time-consuming, and the counterpart of an
efficient propagation phase is the time spent in the collision phase.
The node-wise storage can be less memory-efficient, since the number of bits
per words may not be close to an integer multiple of the number b of channels
per node. The propagation phase is less straightforward than in the channel-wise
storage. Indeed, the post-propagation state of one node must be built from b
bits extracted from b neighboring nodes, and collected in a single word. The
result of this operation must be stored in a buffer copy of the whole lattice
configuration. Indeed, if the post-propagation state of the node were stored back
in its address in the original array storing the pre-propagation state, then the
post-collision state of the neighboring nodes would be computed with already
propagated data. This is a consequence of the non-locality of the propagation
phase. The collision phase however is much simpler, at least for look-up table
strategies, since each word (or fraction of word) of the buffer copy can directly
serve as an address for the look-up table. The results furnished by the look-up
2
This is the case for most models, except the simplest HPP and FHP models.
10.2 Lattice gas algorithms on general purpose computers 213
table are then stored back in the original array, so that they are ready for the
next propagation phase.
To sum up the discussion, the channel-wise storage is in general preferable for
models with simple collisions and/or if the available amount of memory is the
bottle-neck of the computer to be used; the channel-wise storage offers better
efficiency for models with complicated collision schemes, and/or on a computer
where the real bottle-neck is the CPU-time.
10.2.3 Obstacles
When the node-wise storage mode is selected, there exists a simple and efficient
way to implement bounce-back boundary conditions (see Section 2.4.1). Suppose
that the b microscopic velocity vectors Cj are labeled such that Ci+b/2 = — c,-. If
the Boolean states are stored node-wise in computer words, so that the most
significant bit of each word stores the occupation level of channel 1, and so
on in sequence up to channel b, then the bounce-back collision is simply a
permutation of the b/2 most significant bits with the next b/2 following bits.
This requires only two masking operations, two shifts, and an 'OR' operation.
In addition, it is useless to apply this bounce-back rule in the bulk of an
obstacle, provided all obstacle nodes are initialized with zero particle density. It is
10.3 Essential features of a lattice gas simulation code 215
sufficient to apply the bounce-back rule to the 'surface' of the obstacle only, that
is, to the obstacle nodes that can be occupied because of particle propagation.
Indeed, the bounce-back rule sends back (to where it came from) any particle
arriving at the surface, so preventing any fluid particle 'contaminating' the
bulk of the obstacle. The thickness of the 'surface' layer depends on the set
of microscopic velocity vectors: if some of these vectors connect a node to its
second or third nearest neighbors, then the layer must be two or three times
thicker.
In the previous section, we addressed the main problems and choices involved in
designing the hard core of a lattice gas simulation code: the automaton updating
machinery. However, such a code cannot be reduced, even conceptually, to its
hard core. Many auxiliary software tools are needed to properly operate the
code and to extract qualitative physically relevant data. In addition to pure
programmer skill, these auxiliary tools sometimes require deep physical insight.
These tools can be roughly classified into three main categories: initialization
tools, raw physical data extraction, and post-processing.
10.3.1 Initialization
Like any computational fluid dynamics method, the lattice gas method requires
one to decide what should be the initial macroscopic fields (uniform, harmonic,
etc.). Furthermore, for lattice gas methods, there is an additional task before
starting the time-evolution: the initial microscopic Booleanfieldmust be set so as
to correspond to the desired initial macroscopic fields. This involves a random
initialization of the Jt;(r*) with 0s or Is, with the minimal constraint that the
resulting averaged macroscopic fields be actually the desired initial macroscopic
fields. Often there is an additional constraint: the random initialization of the
tt;(r*) must be realized with probabilities that correspond to (local) equilibrium
distributions, such as those discussed in Chapter 4. This eliminates (or short-
ens) the transient regime during which local equilibrium would have to settle.
Therefore the average populations fi(f*) corresponding to the desired initial
macroscopic fields must be computed according to the relations obtained in
Section 4.5. Then, at each node and for each channel, a pseudo-random number
generator must be used to assign to n^r*) the value 0 or 1, with a probability
/i(r*). A convenient way to do this is to use a random number generator deliv-
ering values cp with a uniform distribution between 0.0 and 1.0 (most built-in
random number generators do so); the integer part of cp +/,(r*), is then either
0 or 1, with respective probabilities /j(r*) and 1 — /*(r*).
216 10 Lattice gas simulations
If the initial macroscopic fields are uniform, the probabilities /j(r*) do not
depend on r*. Then they can be pre-computed once, for all the nodes (instead
of once for each node), leading to considerable speed-up.
Note also that if CPU time is a crucial issue - especially for non-homogeneous
initial fields that require computation of the probabilities /,-(r*) for each node
- a simplified procedure can be used to obtain the //(r*), usually with no
dramatic physical consequences: one can use approximation formulae for the
probabilities /;(r*), such as those obtained in Sections 4.5.1 and 4.5.2, but taking
into account the first order correction terms only. This simplifies the computation
considerably.
larger the cells, the smaller the residual noise, and the poorer the effective space-
and/or time-resolution. So there should be a compromise between the loss in
resolution and the loss in accuracy due to the residual noise. This averaging
procedure can also be viewed as a 'poor man's' filtering process, where the
microscopic data are convoluted with a square function that is equal to one in a
small space or space-time domain, and zero everywhere else. This interpretation
gives a natural generalization of the space- and/or time-averaging: the square
function can be replaced by any smooth function whose Fourier transform de-
cays smoothly at high space- and/or time-frequencies. However, the compromise
between noise reduction and effective resolution still remains: setting the cut-off
frequency too low will efficiently reduce the noise, but may reject physically
interesting small structures; conversely, setting the cut-off frequency too high
will produce prohibitively noisy fields.
10.3.3 Post-processing
Like standard floating point computational fluid dynamics methods, lattice gas
simulations often require serious processing for drawing physically significant
results from raw data, i.e. the averaged macroscopic fields. This may involve
physical aspects such as computing macroscopic fields different from those
delivered naturally by the lattice gas simulation. For example, one may wish to
compute the pressure field and the vorticity field from the density and velocity
fields. This phase may also involve graphical aspects, such as choosing the
right way to display the results to highlight a particular phenomenon, which is
of crucial importance, especially for three-dimensional simulations. Since these
problems are common to all computational fluid dynamics methods, we shall
not embark on a detailed discussion here. We shall focus attention on post-
processing problems involving specific aspects of the lattice gas method, such as
drag and lift coefficients measurements, and Strouhal number evaluations.
CO..)
where p is the mass density, u is the mean flow velocity, and S is the area of
transverse section of the obstacle, that is, the area of projection of the obstacle
onto a plane orthogonal to the mean flow velocity.
218 10 Lattice gas simulations
The transverse component of the force, expressed in the same way, yields the
'lift coefficient':
CL = y % - . (10.2)
These quantities are of obvious interest, and can be rather accurately extracted
from a lattice gas simulation (see Rivet, 1993). First, one has to compute the
momentum variation Ap during one time step, due to the collisions on the
obstacle. This must be done at each time step during the simulation, and
the values must be saved to be subsequently post-processed. Then one has to
compute the area of transverse section of the obstacle, in natural lattice units,
and the average mass density per unit volume, in natural lattice units. The volume
of the unit lattice cell must be taken into account if it differs from 1 (in natural
lattice units). The drag and lift coefficients are obtained from:
f^' (10.3)
g(p)\pu2S
and
CL = , f/\. (10.4)
g{p)\pu2S
The factor g(p) appears to express the momentum variation per natural lattice
time unit (time step) as a force, that is, a momentum variation per natural
physical time unit (see Section 9.2). It was found that the values of the drag
coefficients computed for spheres and cylinders with the lattice gas method fall
within a small percentage of experimental values.
St = ^ , (10.5)
where L is the typical size of the obstacle, and u is the mean flow velocity.
The Strouhal number - and how it varies as a function of the Reynolds
number - may be a precious indicator of what happens in the wake (Konig,
Eisenlohr and Eckelmann, 1990). This is the case for example in the wake of
a cylinder, where the three-dimensional structure of the wake may be complex.
The Strouhal number can be computed by 'lattice gas hot-wire anemometry',
i.e. by implementing the numerical analogue of a 'hot-wire velocity probe', a
laboratory device that delivers the local fluid velocity as a function of time (see
10.4 Measurement of basic lattice gas properties 219
Bonetti and Boon, 1989). In the lattice gas simulation, the 'numerical' hot-wire
probe is an averaging cell, over which the averaged velocity is computed at
each time step during the simulation, and stored for further post-processing.
The resulting (discrete) signal is very noisy, but a low frequency mode clearly
dominates. The next step is to get the best possible estimate of the frequency of
this dominant mode, which is done with standard signal processing algorithms
such as modulation and filtering. The frequency so obtained yields the Strouhal
number through the formula:
^ (10.6)
where L and u are expressed in natural lattice units. The lattice gas-specific non-
Galilean factor g(p) appears as a scaling factor between the physical time unit
and the natural lattice time unit (see Section 9.2). Strouhal numbers obtained
from LGA simulations of the wake behind a cylinder are also found to be in
agreement (within a small percentage) with experimental values.
The damping coefficient F is obtained from a least squares fit of the logarithm of
the amplitude as a function of time. Knowing v(p) from shear wave experiments,
and F from sound wave experiments, the value of the bulk viscosity v'(p) follows
straightforwardly.
The sound velocity cannot be accurately extracted by fitting the phase of the
Fourier mode to a linear function. Indeed, the phase speed in the viscous case is
not cs but \/cs2 — T2k2 (see Equation (9.12)). One convenient way to obtain cs
10.4 Measurement of basic lattice gas properties 221
Table 10.1 Non-Galilean factor g(p) and sound speed cs(p) for the
FCHC-3 lattice gas model, d is the density per channel, p = d x b is the
density per node. gth is the theoretical value of g(p) (see Equation (9.3)),
gexp is the value obtained by shear wave numerical simulations, gerr is the
root mean square estimated accuracy of the simulation. The same
subscripts are used for the sound speed cs(p) (see Equation (9.10) for the
theoretical value).
Table 10.2 Kinematic and bulk viscosity coefficients v(p) and v'(p) for
the FCHC-3 lattice gas model. The convention for the meaning of the
subscripts is the same as in Table 10.1.
strategy with one single table. So the FCHC-3 lattice gas is only approximately
G-invariant. Simulations are performed over a domain of 256 x 64 x 64 lattice
nodes (see Dubrulle et al, 1990). Tables 10.1 and 10.2 give the theoretical values,
the measured values and the estimated accuracies, for g(p) and cs(p), and for
the coefficients v(p) and v\p) of the FCHC-3 lattice gas model.
For further illustrations of lattice gas simulations, see Section 11.7 in Chapter 11.
10.5 Examples of lattice gas simulations 223
the second species may serve as a tracer for better visualization of the vortices.
Figure 10.3 shows the simulation of the Kelvin-Helmholtz instability after 1200
time steps. The upper layer flows to the left, and is more concentrated in one of
the two species, whereas the lower layer flows to the right and is more concen-
trated in the other species. The gray levels label the concentration differences.
The two layers are initially separated by a flat interface, and after 1000 time
steps the Kelvin-Helmholtz instability is clearly visible. The patterns in Fig-
ure 10.3 can be compared (qualitatively) to those in Figure 10.4 showing (real)
clouds in two atmospheric layers undergoing a Kelvin-Helmholtz instability.
Figure 10.5 illustrates a simulation with a slightly different collision rule where
the two species are subject to a reactive process according to a majority rule:
if there are more particles of species A on a node, then the particles of the less
abundant species B are transformed into particles of species A. Consequently,
the gradients are sharper, and the dispersive effect of particle diffusion is re-
duced; so one can run the simulation for longer times (4000 time steps) and
avoid fuzzy contours.
Figure 10.4 Photograph of real clouds between two atmospheric layers with
differential motion: the Kelvin-Helmholtz instability occurring between these
layers is made visible by the clouds. (Drazin and Reid, 1981.)
224 10 Lattice gas simulations
role of the initial seed: it is 'sticky' in the sense that any moving particle arriving
at one of the neighboring sites (next to the seed) stops and stays, and becomes
sticky too. In the course of time, an aggregate of fixed sticky particles forms and
grows with an apparently fractal structure. Figure 10.6 showing such a fractal
aggregate was obtained after just a few seconds of simulation on RAP-1.
Figure 10.8 Gray-level coding of the vorticity field corresponding to the velocity field in
Fig. 10.7.
10.5 Examples of lattice gas simulations 227
vortices) appear in black. The square 'pixels' clearly visible in Figure 10.8 are
the averaging cells. Purposely, no smoothing was performed in this picture,
because it would have hidden the coarse discretization of the physical space,
imposed by the necessity of space-averages. As explained in Section 10.3.2, the
size of the averaging cells results from a compromise between the accepted noise
level, and the desired physical resolution: larger cells would have led to less noisy
values of the velocity field, at the expense of a poorer space-resolution; smaller
cells would have led to a better resolution, with an increased noise level. With
both space- and time-averaging, one would have a gain in 'signal-to-noise' ratio
for a given space-resolution, but at the cost of a loss in time-resolution.
FCHC-8 model. Figure 10.9 illustrates the geometry of the simulation. The
solid arrow labeled 'I/' shows the direction of the incoming flow. The boundary
conditions are periodic in the transverse direction X and in the span-wise
direction Y. In the stream-wise direction Z, the inlet and outlet conditions are
as described in the two-dimensional simulation of Section 10.5.3. The upstream
velocity is 0.19 in lattice units, and the corresponding Mach number and
Reynolds number are Ma = 0.29 and Re = 73, respectively. The Mach number
has a low enough value to maintain the density fluctuations within 1.5% of the
mean value d = 0.42 (for more details, see Rivet, 1993).
Figure 10.10(a) is a snapshot taken 16000 time steps (74.4 circulation times)
after the impulsive start of the flow. It shows a two-dimensional cross-section
through the fluid velocity field, in a plane perpendicular to the cylinder axis. A
Benard-von Karman vortex street is clearly visible in the wake of the cylinder,
as in the two-dimensional case. The vortices are evidenced in Figure 10.10(b)
with a gray-level encoding of the span-wise (Y) component of the vorticity, in
the same plane. The macroscopic velocity field is obtained by space-averaging
over 8 x 8 x 8 cubic cells, and the vorticity field is computed a posteriori from the
macroscopic velocity field, during the post-processing phase; the square pixels
are the cross-sections of the cubic space-averaging cells. No smoothing procedure
has been applied to the picture for the reason discussed in Section 10.5.3.
Any cross-section by a plane perpendicular to the cylinder would show the
same kind of structure: oscillating vortex streets with alternately positive and
negative vortices detaching from the obstacle. But are all these oscillating vortex
streets phase-synchronized? In other words, are the vortex rolls straight and
parallel to the cylinder, as we would naively expect? To answer this question,
we make cuts by planes parallel both to the cylinder and to the mean velocity
direction. Figures 10.1 l(a), 10.11(b) and 10.1 l(c) are three snapshots taken
respectively 10000, 16000 and 24000 time steps (46.4, 74.4 and 111.7 circulation
times) after the impulsive start of the flow. The pictures give a finite-distance
perspective view of the simulation domain, seen from a direction roughly at
55° from the X axis. They show a two-dimensional section in a plane tangent
to the cylinder and parallel to the stream-wise (Z) direction. In this plane, the
velocity field is represented so as to make the span-wise structure of the vortex
rolls clearly visible; the transverse (X) component of the macroscopic velocity
Figure 10.11
Three-dimensional
incompressible flow
around a cylinder in
wind-tunnel geometry,
simulated with the
FCHC-8 lattice gas model
on a CRAY-2 computer.
The X component of the
velocity is displayed by
gray levels in a plane
tangent to the cylinder
and parallel to the
stream-wise (Z) direction.
(a) 10000 time steps (46.4
circulation times) after
the impulsive start; a
sequence of wavy vortex
rolls is visible in the wake.
(b) 16000 time steps (74.4
circulation times) after
the impulsive start; a
dislocation appears in the
vortex pattern, (c) 24000
time steps (111.7
circulation times) after
the impulsive start; the
dislocation has
disappeared and the
vortex rolls are now tilted.
V
10.5 Examples of lattice gas simulations 231
Figure 10.12
Two-dimensional lattice
gas simulation of a
two-phase fluid with
surface tension in a
porous medium. The
interface is initially flat;
an instability develops
and produces the typical
fingering patterns of the
less viscous fluid (in
black) penetrating the
most viscous fluid (in
gray). The four snapshots
are taken after 30000 (a),
50000 (b), 100000 (c),
and 112 500 (d) time
steps.
Chapter 11
Since 1985, lattice gas automata have become a widely and actively explored
field. Academic groups and industrial laboratories around the world have in-
vested considerable effort in their research activities to drive the subject in
various new and promising directions. As a result, the literature on lattice gases
and related topics has grown so rapidly and has become so voluminous that
an exhaustive list and a detailed review of all relevant lattice gas publications
would practically make a book by itself.1
Here we select research areas where lattice gases have played an important
role, and, for each of them, we quote articles considered as representative,
historically or presently. Our review is by no means complete and our choices
are certainly selective; unavoidably we haven't done justice to those whose work
escaped our selection. Nevertheless we have attempted to cover a range of works
expanding over various aspects of the subject in the available literature. Our
goal will be reached if the reader finds this chapter a helpful tool for exploring
the subject beyond the scope of this book.
1
Systematic explorations of the literature were performed by G.D. Doolen and by
D. d'Humieres from 1986 to 1995; alphabetical lists of references (often with abstracts)
can be found in special issues of Physica D, 47, pp. 299-337, 1991 (for the period
1986-1990), J. Stat. Phys., 68, pp. 611-669, 1992 (for the period 1990-1992), and Fields
Inst. Comm., 6, pp. 275-346, 1996 (for the period 1992-1995).
233
234 11 Guide for further reading
The promising results of the lattice gas method for fluid dynamics simulations
in two dimensions stimulated the search for three-dimensional models. Unfor-
tunately, the 'miracle' that renders the macroscopic behavior of FHP models
isotropic does not happen so naturally in three dimensions. Several solutions
were proposed, including the 'FCHC models' that reside on a four-dimensional
lattice (see Section 3.7.3). But these models have at least 24 channels per node,
and the design of collision rules soon becomes a subtle optimization puzzle.
Indeed, the collision rules must satisfy the correct physical properties (conserva-
tion laws, G-invariance, possibly semi-detailed balance, etc.), and in addition -
for high Reynolds number flow - they should yield the smallest possible value
for the kinematic viscosity. Moreover, they must be numerically efficient, and
shouldn't of course exceed the computer's technical limitations. These questions
are discussed in detail in the following papers where the interested reader will
also find elaborations on collision rules strategy.
& d'Humieres, D., Lallemand, P. and Frisch, U., 'Lattice gas models for 3-D
hydrodynamics', Europhys. Lett, 2, pp. 291-297, 1986.
& Henon, M., Isometric collision rules for the 4-D FCHC lattice gas',
Complex Systems, 1, pp. 475-494, 1987.
¥> Henon, M., 'Optimization of collision rules in the FCHC Lattice Gas and
addition of rest particles', in the proceedings of the workshop (Discrete
Kinetic Theory, Lattice Gas Dynamics and Foundations of Hydrodynamics',
Torino (Italy), 20-24 September 1988, ed. R. Monaco, World Scientific,
Singapore, 1989.
& Somers, J.A. and Rem, P.C., 'The construction of efficient collision tables
for fluid flow computations with cellular automata', in Cellular Automata
and Modeling of Complex Physical Systems, pp. 161-177,
ed. P. Manneville, Springer-Verlag, Berlin, 1989.
Initially, rather than tools for computer simulations, lattice gas models were
considered as 'test-benches' for theoretical ideas and methods. So it is not
surprising that the literature on lattice gas theory is so rich, since many aspects
of the subject have been explored by various approaches. In what follows we
11.3 Theoretical analyses 237
offer the reader our selection of references which we have organized in five
categories.
There have been several attempts to extend lattice gas automata to the
quantum domain. Models have been proposed for the Dirac equation and for
the Schrodinger equation. Also motivated by the observation that a natural
architecture for nanoscale quantum computation is that of a quantum cellular
automaton, authors have developed theoretical analyses and algorithmic schemes
for quantum lattice gas automata. Here are some references on a subject which,
although still in its infancy, looks promising.
• Dubrulle, B., 'Method of computation of the Reynolds number for the two
models of lattice gas involving violation of semi-detailed balance', Complex
Systems, 2, pp. 577-609, 1988.
* Dubrulle, B., Frisch, U., Henon, M. and Rivet, J.R, 'Low-viscosity lattice
gases', J. Stat. Phys., 59, pp. 1187-1226, 1990.
240 11 Guide for further reading
• Lavallee, P., Boon, J.P. and Noullez, A., 'Boundaries in lattice gas flows',
Physica D, 47, pp. 233-240, 1991.
• Cornubert, R., d'Humieres, D. and Levermore, D., A Knudsen layer
theory for lattice gases', Physica D, 47, pp. 241-259, 1991.
11.4 Models with particular features 241
11.4.5 Thermo-hydrodynamics
-: Chopard, B. and Droz, M., 'Cellular automata model for heat conduction
in fluid', Phys. Rev. Lett. A, 126, pp. 476-480, 1988.
; Luo, L.S., Chen, H., Chen, S., Doolen, G.D. and Lee, Y.C., 'Generalized
hydrodynamic transport in lattice-gas automata', Phys. Rev. A, 43,
pp. 7097-7100, 1991.
: Chen, S., Chen, H., Doolen, G.D., Gutman, S. and Lee, Y.C., A lattice
gas model for thermodynamics', J. Stat. Phys., 62, p. 1121, 1991.
: Ernst, M.H. and Das, S.R, 'Thermal cellular automata', J. Stat. Phys., 66,
pp. 465-483, 1991.
: Grosfils, R, Boon, J.R and Lallemand, R, 'Spontaneous fluctuation
correlations in thermal lattice gas automata', Phys. Rev. Lett., 68,
pp. 1077-1080, 1992.
: Das, S.R and Ernst, M.H., 'Thermal transport properties in a square
lattice gas', Physica A, 187, pp. '191-209, 1992.
-; Mora, P., 'The lattice Boltzmann phononic lattice solid', J. Stat. Phys.,
68(3/4), pp. 591-610, 1992.
-: Maillot, B., 'Modeles semi-microscopiques d'ondes elastiques, Ph.D. Thesis,
Institut de Physique du Globe, Paris (France), June 1994.
The linearized Boltzmann collision operator (see Section 11.5) can be simplified
by reducing its matrix to a diagonal form with one single non-zero eigenvalue.
This method proposed initially in 1954 by Bhatnagar, Gross and Krook amounts
to replacing the full spectrum of eigenvalues by one single relaxation time. The
method applies straightforwardly to the lattice Boltzmann scheme and is quite
useful in particular for tuning, within certain limits, the kinematic viscosity to
lower values than can be achieved by the LGA method. The procedure known
as the lattice BGK method (LBGK) is discussed from the point of the theory
as well as for various applications in:
: Bhatnagar, P., Gross, E.P. and Krook, M.K., A model for collision
processes in gases', Phys. Rev., 94, pp. 511-525, 1954.
: Qian, Y.H., d'Humieres, D. and Lallemand, P., 'Lattice BGK models for
Navier-Stokes equation', Europhys. Lett., 17, pp. 479-484, 1992.
-: Qian, Y.H. and Orszag, S.A., 'Non-linear correction to Navier-Stokes
equation derived from lattice BGK models', Europhys. Lett., 21,
pp. 255-259, 1992.
: Succi, S., d'Humieres, D., Qian, Y.H. and Orszag, S.A., 'On the small-scale
dynamical behavior of lattice BGK and lattice Boltzmann schemes', J. Sci.
Comp., 8, pp. 219-230, 1993.
u Qian, Y.H., 'Simulating thermohydrodynamics with lattice BGK models',
J. Sci. Comp., 8, pp. 231-242, 1993.
• Flekkoy, E.G., 'Lattice Bhatnagar-Gross-Krook models for miscible
fluids', Phys. Rev. E, 47(6), pp. 4247-4257, 1993.
Because of the binary and discrete nature of cellular automata, and because
of the global updating of the LGA, lattice gases are very good candidates
for easy and efficient implementation on general purpose massively parallel
(super)computers, as well as on dedicated computers. The next two sub-sections
give references to articles on dedicated hardware implementations, and general
purpose computer implementations and simulations.
-: Somers, J.A. and Rem, P.C., 'Flow computation with lattice gases', Appl.
Sci. Res., 48, pp. 391-435, 1991.
v Henon, M. and Scholl H., 'Lattice-gas simulation of a nontransverse
large-scale instability for a modified Kolmogorov flow', Phys. Rev. A, 43,
pp. 5365-5366, 1991.
-; Somers, J.A. and Rem, P.C., 'Obtaining numerical results from the 3D
FCHC-lattice gas', Lecture Notes in Physics, 398, pp. 59-78,
Springer-Verlag, Berlin, 1992.
-; Henon, M., 'Implementation of the FCHC lattice gas model on the
Connection Machine', J. Stat. Phys., 68(3/4), pp. 353-378, 1992.
-: Rivet, J.P., 'Spontaneous symmetry breaking in the 3-D wake of a long
cylinder simulated by the lattice gas method and drag coefficient
measurements', Appl. Sci. Res., 51, pp. 123-126, 1993.
We close with a list of books and proceedings volumes published since 1986;
we also include general review articles on lattice gases and related topics.
; Wolfram, S. (editor), Collection of contributions to the workshop on
Large Non-Linear Systems, Santa Fe, (NM, USA), 27-29 October 1986;
Published as a special issue of Complex Systems, 1(4), 1987.
-: Frisch, U., 'Une nouvelle strategie pour Vhydrodynamique: les reseaux
d'automates, J. Astr. Fr., 32, pp. 17-20, 1988.
-: Boon, J.P., 'Physique, complexity et automates', Bull. Soc. Fr. Phys.,
October 1989.
; Monaco, R. (editor), Proceedings of the workshop on 'Discrete Kinetic
Theory, Lattice Gas Dynamics and Foundations of Hydrodynamics, Torino
(Italy), 20-24 September 1988, World-Scientific, Singapore, 1989.
-; Doolen, G.D. (editor), Lattice Gas Methods for Partial Differential
Equations, Addison-Wesley, 1990.
-; Boon, J.R, 'Lattice gas automata: a new approach to the simulation of
complex flows', in Microscopic Simulations of Complex Flows, ed.
M. Mareschal, Plenum Press, New York, 1990.
• Doolen, G.D. (editor), Lattice Gas Methods for Partial Differential
Equations: Theory, Applications and Hardware, special issue of Physica D,
47(1/2), 1991. (Proceedings of the NATO Advanced Research Workshop
on Lattice Gas Methods for Partial Differential Equations: Theory,
Applications and Hardware, Los Alamos National Laboratory (NM,
USA), 6-8 September 1989).
11.8 Books and review articles 249
Chen, S., Doolen, G.D. and Matthaeus, W.H., 'Lattice gas automata for
simple and complex fluids', J. Stat. Phys., 64(5/6), pp. 1133-1162, 1991.
Ernst, M.H., 'Statistical mechanics of cellular automata fluids', in the
proceedings of the NATO ASI on Liquids, Freezing and the Glass
Transition, Les Houches (France), July 3-28, 1989, eds. J.P. Hansen,
D. Levesque and J. Zinn-Justin, North Holland, Amsterdam, 1991.
Boon, J.P. (editor), special issue of J. Stat. Phys., 68(3/4), 1992.
(Proceedings of the NATO Advanced Research Workshop on Lattice Gas
Automata Theory, Implementation, and Simulation, Observatory of Nice
(France), 25-28 June 1991).
Boon, J.P., Frisch, U. and d'Humieres, D., 'L'hydrodynamique modelisee
sur reseau, La Recherche, 253, pp. 390-399, April 1993.
Rothman, D.H. and Zaleski, S., Lattice-Gas Cellular Automata,
Cambridge University Press, Cambridge (UK), 1997.
Chopard, B. and Droz, M., Cellular Automata Modeling of Physical
Systems, Cambridge University Press, Cambridge (UK), 1998.
Wolf-Gladrow, D., An Introduction to Lattice-Gas Cellular Automata and
Lattice Boltzmann Models, Springer-Verlag, Berlin, 1999.
Succi, S., The Lattice Boltzmann Equation for Fluid Dynamics and Beyond,
Oxford University Press, (to be published).
Appendix
Mathematical details
Proof
Consider a non-local collision model with b channels per node. The collision
phase is such that the post-collision state n^-(r*) of any node labeled by its
position vector r* is a function J^ of the pre-collision state n,-(r* — e^) of B nodes
(the node itself plus B — 1 other nodes), with position vectors r* — e&:
In the most general case, the non-local collision process at a given node can
also involve the state of the node itself. So, one of the e^ s (say ei for example)
must be the null vector. In addition, the number B of 'collision vectors' e& is
finite, and the e^ s are globally invariant under the action of the point-group of
the lattice.
The propagation phase involves the set of b propagation vectors c* which
are also globally invariant under the point-group of the lattice. The final (post-
propagation) state nj(r*)(i=iv.^) of node r* is given by:
We define a new model with the same lattice structure, but with b = b x B
250
A.2 251
channels per node. As before, each node is labeled by its position vector r*, but
now each channel is labeled by a composite label ? = (i, k) with i = 1,..., b and
fc = l , . . . , B .
The collision phase of the new model is such that the post-collision state
n'fj^r*) of any node r* is given by:
The new collision process is now purely local, since the right hand side of (A. 3)
involves only the state n^k) of node r*.
The propagation phase of the new model involves the propagation vectors
c(l/c) = a + e/c, with i and k ranging from 1 to b and from 1 to B, respectively.
The final (post-propagation) state h^(u)(i=it^b;k=u...9B) of node r* is given as
usual by:
It is clear that if the initial state %/c)(r*, U = 0) of the new model is set
such that %/c)(r*, U = 0) = 7i;(r* — t^U = 0), then, at any further discrete time
£•, the state n^k)(u,t*) will verify W(i,i)(r*,£*) = ni(r*,t*). In other words, the
exact dynamics of the initial model with non-local collisions is 'contained' in
the dynamics of the new model with purely local collisions. This shows that a
model with non-local collisions can be re-coded as a model with more channels
per node, but with local collisions. •
Proof
The four-dimensional FCHC lattice is the set of points of R 4 whose coordinates
are signed integers with even sum (see Section 3.7). The basic simple hyper-
cubic lattice on which the FCHC lattice is constructed is the set of all points
with coordinates (2/?i,2w2,2ft3,2ft4) where (^1,^2,^3,^4) are signed integers. Con-
sider the elementary hyper-cube limited by 16 lattice points with Cartesian
252 Appendix
coordinates:
0
0
\2j
/2\ f2\ /2\ ( \ r\
( \ ( \
2 2
P 2 2 2
P 2
2
\2/
2
\2/
'\2P ' 0 ' 2 ' 2 ' 0
\ 2 / \o/ \ 2 / \2/
2
\0/
The center of this hyper-cube is the point (1,1,1,1), which clearly belongs to
the FCHC lattice since the sum of its coordinates is even. So, the FCHC lattice
is 'body-centered'.
This hyper-cube has 8 (hyper-)faces which are obtained by setting successively
each of the 4 coordinates equal to either 0 or 2. These hyper-faces are three-
dimensional cubes. For instance, the hyper-face which verifies X4 = 0 is limited
by the 8 following lattice points:
2
0\ /2\ /o\
0 0 2
^°\0 /22\ /2\
0
/2^/o2 (A.6)
2
0 ' 0 ' 0
' p 'w 2 •2 •> 2
5
0/ \0 \o/ \o
The center of this hyper-face is the point (1,1,1,0) which does not belong to
the FCHC lattice. Thanks to the invariances of the lattice, it is easy to convince
oneself that none of the 8 hyper-faces are centered. So, the FCHC lattice is not
'face-centered', at least if 'face' is understood as three-dimensional 'hyper-face'.
The hyper-cube has 24 2D-faces which are obtained by setting successively
each pair of coordinates equal to either (0,0), (2,0), (0,2) or (2,2). For instance,
the 2D-face which verifies (^3,^4) = (0,0) is limited by the 4 following points:
/2\
0 0
(A.7)
w 0 0
\oj
The center of this 2D-face is the point (1,1,0,0) which clearly belongs to the
FCHC lattice. Using again the invariances of the FCHC lattice, it is easy to
show that all 24 2D-faces are centered. Thus the FCHC lattice is also '2D-face-
centered'. •
A.3 253
Proof
We compute the difference S between the quantity of information contained
in the knowledge of all the individual statistics and the knowledge of the full
statistics (see Section 4.3.1 for notation):
8 = r-i
q
1=1
with p. = (1 — p^ We now use the fact that for a Bernoulli variable xt, the
probability p\ that x\ be equal to 1 is also the averaged value of xf, namely:
and ft =
YeT YeT
Now considering the inequality In a < a — 1, which holds for any strictly positive
value of a, the difference 5 verifies the inequality:
(5<0.
and
7=1 sey
f
7=1
Proof
We first factorize the r.h.s. of (A.9) as follows: we single out a particular node
ro in i^, and denote by if (o) the set of all nodes of if except i*o. The r.h.s. of
(A.9) can then be rewritten as:
E E ^s'^mu'f/1 n Mt^iiftf?' ,
sey SeT, 7=1 r*ej?(0) 7=1
S()
iS and 5' taken at node r*
sey
n z^^n^//7-
sey 7=1
A.5 255
(A,2)
b
Y^ Wi - Si)A(s^s') J ] fffj> =0, Vi = 1,..., b, (A. 14)
s,s'ey ;=1
Z)log (T^J) ti
are strictly equivalent when the semi-detailed balance condi-
tion is satisfied (see Section 4.4.1)
Proof
• We prove that (A.13) ^> (A.14):
We use the normalization condition (2.9) to write (A.13) under the following
equivalent form:
On the l.h.s., both s and sr are dummy variables, which we can permute to
obtain:
We now use a trick introduced by Gatignol in the context of gas models with
discrete velocities (Gatignol, 1975): we multiply both sides by log// and sum
over I Equation (A.19) becomes:
b b
sey sey
an
and we multiply both sides by n^=i fj d sum over
1
We exclude the pathological case /& = 1.
A.5 257
s,s'ey \;=1
\ w; -wl+wi\=0.
l
(A.23)
where all products are taken from j = 1 to b. The expression between braces is
always negative or zero. Indeed, denoting by x the argument of the logarithm,
this expression has the form log x — x + 1 which is zero for x = 1 and whose
first derivative (= 1/x — 1) is positive for x < 1, null for x = 1 and negative
for x > 1. Thus, logx — x + 1 must be negative for any value of x ^ 1. As
a consequence, (A.23) implies that the expression between braces must be zero
whenever A(s^>sf) is non-zero. In other words, Il/=i Tj a n d Ily=i Tj m u s t be
equal for any states s and sr with a non-zero collisional transition probability
A(s^sf), that is:
(
5 6 \
4^) n fW=A(S-S') n f
;=i 7=1
Summing both sides over all s in y and using the semi-detailed balance condition
(2.17) gives (A.13). •
258 Appendix
Proof
We use the normalization condition (2.9) to write the l.h.s. of (A.27) under the
following form:
sey 7=1
The sum over s can be restricted to the Boolean states 5 with s* = 1. The l.h.s.
of (A.27) then reads:
sey, 7=i,
The sum is obviously equal to one, since it is taken over all possible Boolean
states of the b — 1 remaining channels (channel i excluded). •
Proof
We denote by E the canonical real vector space of dimension b, and by /
the S -dimensional vector subspace of collisional invariants. The set of vectors
(q[K\ K = 1,...,<5) is a basis of 7. The vectors in E will be denoted by Roman
letters rather than bold faced symbols, which we use for vectors in the physical
space IRD.
We introduce the symmetric quadratic form acting on E x £, which associates
to any pair of vectors (x,y) € E x E, the real number:
x
' y =
A.8 259
P P — °KKr> K
> K
— 1,...,0.
We introduce the real d x d square matrix R which transforms the q^ s into the
p M s. Its matrix elements R^ are such that:
A.8 Computation of the r.h.s. of (4.47) for the class of models de-
fined in Section 4.5.2
\f4t^r'])\ (A-29)
i=i V=i /
where the first order corrections k^ to the Lagrange multipliers are, according
to (4.71):
(A30)
260 Appendix
We recall the definitions of the characteristic coefficients {2, £4 and & (see
Equation (4.64)) which depend on the geometric structure of the lattice gas:
1=1
b
;
)
We also recall that the collisional invariants q^ of the class of models under
consideration are (see Equation (4.65)):
(A.32)
V e . = iicc 2 __ Mi
2H 2b
The expansion coefficient fii is d(l — d)(l — 2d), according to (4.54). Using the
definitions (A.31) and (A.32), one can prove the following useful identities:
**
t=
e, = ^6, 2^ ciaCijiei = [— + — J
1=1 1=1
Starting from (A.29), and taking into account (A.30), we can rewrite S^ as:
S[K] 1 2 J
1
(A.34)
We first examine the case K = 0. The summation over i can be performed using
the explicit expression (A.31) for the crystallographic tensor J2C^CW- Because
crystallographic tensors of odd order are null, the second term in (A.34) does
not contribute, and the expression for S^ reads:
(A 35)
-
For the case K = a, a = 1,...,D, the algebra is quite similar. The first and
third terms in (A.34) do not contribute since for q^ = c^ they involve odd
order crystallographic tensors. The second term is computed with the second
identity in (A.33). This gives:
Jf = S[K\ (A.38)
K' = l
262 Appendix
where the r.h.s. terms S M are given by (A.35), (A.36) and (A.37). Since the
matrix M^K'^ is diagonal (see Section 4.5.2), the second order perturbations /l|*']
follow immediately:
a = 1
>-'D> (A39)
By definition (see Section 2.1.3), momentum and effective energy fluxes are:
b
1=1
(A.40)
where the /i (0) s are the equilibrium mean populations derived in Section 4.5.2
for the nearly equally distributed states of single-species thermal models:
f.(0) _ P
~ b
1 1
-rJyCiy + T-wei
Q2 C4
b(b-2p)
2p(b — p) ^2 ^ (A 41)
2
2/ 2
equally distributed equilibrium state with zero momentum density and effective
energy density.
The key-ingredients to compute the fluxes are the following relations:
1=1
)x
^ =l
\J><xpdy8 + <
( A - 42 )
1=1
• E<
These relations are direct consequences of the definition (4.64) of the geometrical
parameters £2, £4 and & (see Section 4.5.2). To these relations, we must also
add the fact that the sum over i of any product containing an odd number of
components of c, is zero (see Section 2.3.5).
With these relations, the computation of the momentum flux tensor is rel-
atively straightforward: all terms in (A.41) with odd powers of j vanish, and
using the relations (A.42) yields the desired expression (5.30):
and
For the energy flux, all terms in (A.41) with even powers of j vanish, and one
obtains (5.33):
with
264 Appendix
We recall the second order solvability conditions (5.57) for the Chapman-Enskog
expansion:
for all K = 1,..., (5, that is, for all collisional invariants.
Before analyzing separately the cases of the mass, momentum, and effective
energy invariants corresponding respectively to K = 0, K = 1,..., D and K = D+l,
one can simplify (A.43) without any loss of generality.
First, the term J2i Q^ft^ vanishes for all K, since it represents the contribution
of the first order corrections /, (1) to the collisional invariants' densities, and this
contribution is zero by the definition of f^ (see Equation (5.4) and the preceding
discussion).
Second, one can use the first order solvability condition (5.25):
to simplify the last two terms in (A.43). The resulting simplified second order
solvability condition reads:
o = dt2p^ + di?J2^^i
t '"' b ( A - 44 )
for all K = 1,...,(5, that is, for all collisional invariants. In (A.44), p^ is the
density Yli^fi^ °ftri e collisional invariant q^K\ and O M is the corresponding
two terms in (A.44) vanish. In addition, the contribution ]T^ c^f^ of /, (1) to
the mass flux (i.e. momentum density) vanishes. This leads to the second order
solvability condition (5.58) for the mass invariant:
Kv = o.
To handle the case of the momentum invariant, we need the expression (5.30)
of the zeroth order momentum flux J2t cia.CiPfP\ that we truncate at the second
order in r\:
{0)
( ) h «, P = 1,..., D.
We also need the expression (5.53) for the first order momentum flux ^ . c^c^f^:
Finally, we must compute J2i CiaCipCiyfi^ up to the first order in rj. This requires
the second identity in (A.42), that is, the very definition of £4. The resulting
expression is:
b
. C.RC. f . ( o > _
(c)
0 = dt2ja - dlP ( v
We then use the first order solvability conditions (5.35) and (5.37) for the mass
266 Appendix
or in compact form:
0 = 3 t J a -5 1 ^(v ( c ) (3 i a ^ + 5 1 ^ a - - 3 i y
+
2) (it ^ ) (5lJ/? + dlfij* ~ l
+ (9{n2).
This gives the second order solvability condition (5.59) for the momentum
invariant:
where the coefficient v, which a priori depends on (p,j 2 ,w), is given by:
v v
2D(D + 2)
The case of the effective energy invariant is somewhat simpler. We need
expression (5.33) of the zeroth order effective energy flux ^ . etCipf^0\ that we
truncate at the second order in rj:
and we use expression (5.56) for the first order effective energy flux:
b
A. 10 267
i=l
When all these results are incorporated in (A.44), the second order solvability
condition for the effective energy invariant becomes:
/ \
0 = dhw-c
Using the first order solvability condition (5.36) for the momentum, we re-express
the time-derivatives of jp in terms of space-derivatives, and we obtain:
or in compact form:
0 = dt2w-dJ^d^
This leads to the second order solvability condition (5.61) for the effective energy
invariant:
268 Appendix
with:
with:
...B^))Jf9t). (A.47)
The first few cumulants (see also Chapter 4) are given explicitly by the expres-
sions (Abramowitz and Stegun, 1965):
4\v,t) = f\*{T,t),
ic{4)(r,t) = ( K E ( r , 0 ) 4 ) - 3((KE(r,0)2)2
= /^(r, t) (1 - /fE(r, t)) (l - 6/,LE(r, t) + 6 (/fE(r, r)
(A.48)
(Risken, 1984):
/, LE (M) 1 0 0
ff*(r,t) ft*(r,t) 1 0
ft*(t,t) /,LE(r,t) /,LE(M)
(A.49)
m \ )
m=l \ /
which is obtained from the observation that the generating function for the
cumulants:
oo n
Ai (A 54)
- -
Consider this relation for n = 1:
which implies:
K n) r f K
d{A(Ttt)) ' (') = ;"+1V>f) ^ * UI4>^ (*,t). (A.58)
where:
^ A^f^t) = Ap (A.63)
The proof makes use of the fact that the average collision matrix connects only
those configurations with the same value for the conserved quantities. Hence:
-E
A.12 271
MM
J
J *»> {s}{s>}
= Aj • (A.64)
7=1
=E ^ 7=1
4- E
{s}{^}
(1-/^,0) E sM(s^s')F{LsE}(r,0
{S}{5'}
E
= ^ (r,t)4-. (A.66)
The other right eigenvectors of J?LE{r,t), which constitute the set of kinetic
eigenvectors, are:
where the scalar product has been defined in section A.ll. This property can be
272 Appendix
E —*
TC
v v /
I
7
n - (n\
7
V
/
> I V
v I /
>
v— l
T E •
v
V 1
l
T E X
J
^ ^,f»ll
J
Uj b 7=1 \i=l /
^=1
) L E (r,0. (A.70)
^ E ( r , t) Kf(r, t) Aj = 0, (A.73)
7=1
7=1
We define the projection operator ^ , the projector onto the set of constants of
motion, by its action on an arbitrary vector {£/(r, t)ffE(r, t)}f=1:
{PB)j(T9t)f?(T9t)
b
— > IA^(Y t\ A * / AI A\ (T t\ * >d RI(Y t\ v f^(r t\
/ ^ / v*5 ^/ _ / ' \±i|£i/LE V' / _/-*-'/V 5 / J I v J /
The complementary projector i> = 1 — & projects onto the set orthogonal to the
constants of motion, i.e.:
1
(r, t) * 4/) *z(r, t)f\*(r91)
(4(^)(r,0) LE = ^ ^ ( ^ ( M l / ^ M ) = 0, (A.78)
;=i
that is, the quantity (IB) (r, t) has no projections onto the set of constants of
motion. The quantity that appears on the r.h.s. of Equation (8.50) is precisely
of this form, with B/(r, t)ff(r, t) = cr Vri /^(r, t).
b
= ] T exp {-C/ • Vr} ^ (dnk(r, % - \)bntf 9 0)>gJ), (A.79)
with:
Teq)
JJ M«or /'i_^q)) (1"Sl
K
j k=l
(A.80)
Equation (A.79) is iterated to yield:
{-Cl- • Vr}
(A.81)
274 Appendix
with:
= ^exp{crVr}. (A.82)
The equal time correlation function is:
rr'
rr'
Therefore:
00
T=0 rr'
00
T=0 ij
(A.85)
References
Abramowitz, M. & Stegun, LA. 1965. Boon, J.P., editor 1992. Lattice Gas
Handbook of Mathematical Functions. Automata Theory, Implementation and
Dover, New York. Simulation (special issue of J. Stat. Phys.
Alder, BJ. & Wainwright, T.E. 1970. Decay WW))- Plenum Press, New York
of the velocity autocorrelation function, (Proceedings of a NATO Advanced
Phys. Rev. A 1(1), 18-21. Research Workshop, Nice (France),
Appert, C. & Zaleski, S. 1990. Lattice gas
Boon
with a liquid-gas transition, Phys. Rev. > J P > D a b > D > Kapral, R.
Lett. 64(1), 1-4. & Lawniczak, A. 1996. Lattice gas
A i_ rv XTAir o A* • x t n 1 n-»r automata for reactive systems* Phys. Rep.
Ashcroft, N.W. & Mermin, N.D. 1976. - . 55147
Solid State Physics. Holt-Saunders,
Philadelphia. Boon, J.P. & Yip, S. 1980. Molecular
Batchelor, G.K. 1967. An Introduction to Hydrodynamics. McGraw-Hill, New York
Fluid Dynamics. Cambridge University (Reprinted by Dover, New York, 1991).
Press, Cambridge. Brito, R., Bussemaker, H.J. & Ernst, M.H.
1991 A
Bell, H., Mceller-Wenghoffer, H , fluctuation formula for the
Kollmar, A., Stockmeyer. R., Springer, T. non-Galilean factor in lattice gas
& Stiller, H. 1975. Neutron Brillouin automata, J. Phys. A Math. Gen. 25,
scattering in fluid neon, Phys. Rev. A L949-L954.
11(1), 316-327. Brito, R. & Ernst, M.H. 1992. Ring kinetic
theor
Bender, CM. & Orszag, S.A. 1978. y f o r tagged-particle problems in
lattice
Advanced Mathematical Methods for gases> phys- Rev- A 46(2)> 875-887.
Scientists and Engineers. McGraw-Hill, Brito, R., Ernst, M.H. & Kirkpatrick, T.R.
Philadelphia. 1991. Staggered diffusivities in lattice
Bernardin, D. & Sero-Guillaume, O.E. 1990. gas cellular automata, J. Stat. Phys.
Lattice gases mixtures models for mass 62(1/2), 283-295.
diffusion, Eur. J. Mech. B 9, 21. Bussemaker, H., Ernst, M.H. & Dufty, J.W.
Bonetti, M. & Boon, J.P. 1989. Chaotic 1995. Generalized Boltzmann equation
dynamics in open flows : the excited jet, for lattice gas automata, J. Stat. Phys.
Phys. Rev. A 40(6), 3322-3345. 78(5/6), 1521-1554.
275
276 References
Cauchy, A. 1851. Rapport sur un memoire Ernst, M.H. 1990b. Linear response theory
presente par A. Bravais, et intitule for cellular automata fluids,
'Etudes sur la cristallographie', C. R. in Fundamental Problems in Statistical
Acad. Sci. Paris 32(8), 284-289. Mechanics VII, ed. H. van Beijeren.
North Holland
Cercignani, C. 1994. On the ^ ^ Ams terdam.
thermodynamics of a discrete velocity Ernst, M.H. 1991. Temperature and heat
gas, Transp. Theor. Stat. Phys. 23(1-3), conductivity in cellular automata fluids,
1-8. in Discrete Models of Fluid Dynamics,
186
Chapman, S. & Cowling, T.G. 1970. The ~ 1 9 7 ' e d A S ' A l v e s - W o r l d Scientific >
Mathematical Theory of Non-uniform Singapore.
Gases. Cambridge University Press, Ernst, M.H. & Das, S.P. 1992. Thermal
Cambridge. cellular automata fluids, J. Stat. Phys.
Chen, S., Zhao, K.H. & Doolen, G.D. 1989. 66(1/2), 465-483.
A lattice gas model with temperature, Ernst, M.H. & Dufty, J.W. 1990.
Physica D 37(1/3), 42-59. Hydrodynamics and time correlation
fu nctioM f c e l l u l a r
Clavin, P., d'Humieres, D., Lallemand, P. , °' automata, J. Stat.
& Pomeau, Y. 1986. Automates ^ WW), 57-86.
cellulaires pour les problemes a frontieres Ernst, M.H., Hauge, E.H. & van
libres en hydrodynamique a 2 et 3 Leeuwen, J.MJ. 1971. Asymptotic time
dimensions, C. R. Acad. Sci. Paris, behavior of correlation functions, Phys.
serie II 303, 1169-1174. Rev. A 4(5), 2055-2065.
Cornubert, R., d'Humieres, D. Fano, R.M. 1961. Transmission of
& Levermore, D. 1991. A Knudsen layer Information. Wiley, New York (p. 43).
theory for lattice gases, Physica D 47, Fleury, P.A. & Boon, J.P. 1969. Brillouin
241-259. scattering in simple liquids : argon and
Dab, D., Lawniczak, A., Boon, J.P. neon, Phys. Rev. 186, 244-254.
& Kapral, R. 1990. Cellular automaton Forster, D. 1975. Hydrodynamic
model for reactive systems, Phys. Rev. Fluctuations, Broken Symmetry, and
Lett. 64(20), 2462-2465. Correlation Functions. Benjamin,
Das, S.P., Bussemaker, HJ. & Ernst, M.H. Reading, Massachusetts.
1993. Generalized hydrodynamics and Frenkel, D. & Ernst, M.H. 1989.
dispersion relations in lattice gases, Phys. Simulation of diffusion in a
Rev. A 48(1), 245-255. two-dimensional lattice-gas cellular
Doolen, G., editor 1990. Lattice Gas automaton: a test of mode-coupling
Methods for Partial Differential theory, Phys. Rev. Lett. 63, 2165-2168.
Equations. Addison-Wesley, New York Frisch, U., Hasslacher, B. & Pomeau, Y.
(Reprints collection). 1986. Lattice gas automata for the
Drazin, P.G. & Reid, W.H. 1981. Navier-Stokes equation, Phys. Rev. Lett.
56
Hydrodynamic Stability. Cambridge ' 1505-1508.
University Press, Cambridge. Frisch, U., d'Humieres, D., Hasslacher, B.,
Dubrulle, B., Frisch, U., Henon, M. Lallemand, P., Pomeau, Y. & Rivet, J.P.
1987
& Rivet, J.P. 1990. Low viscosity lattice - L a t t i c e § a s hydrodynamics in two
and three
gases J Stat Phys 59 1187-1226 dimensions, Complex Systems
' 1, 649-707.
Eliott, J.P. & Dawber, P.G. 1979. Symmetry
in Physics. MacMillan Press Ltd., Gatignol. R. 1975. Theorie Cinetique des
London ^az a Repartition Discrete de vitesses.
Springer-Verlag, Berlin.
Ernst, M.H. 1990a. Statistical mechanics of
Grosfils
cellular automata fluids, in Liquids, ' R 1 9 9 4 ' Hydrodynamique
Freezing and the Glass Transition, 1-94, statistique des gaz sur reseau, Ph.D.,
eds. D. Levesque, J.P. Hansen & Universite Libre de Bruxelles.
J. Zinn-Justin. Elsevier Science, Grosfils, P., Boon, J.P. & Lallemand, P.
Amsterdam. 1992. Spontaneous fluctuation
References 277
correlations in thermal lattice gas van der Hoef, M.A. & Frenkel, D. 1991.
automata, Phys. Rev. Lett. 68, 1077-1080. Evidence for faster-than-r 4 decay of the
velocit
Grosfils, P., Boon, J.P., Brito, R. & Ernst, y autocorrelation function in a 2D
fluid P /
M.H. 1993. Statistical hydrodynamics of ' ^ Rev< Lett 6 6 ' 1591-1594.
lattice gas automata, Phys. Rev. E 48, Homsy, G.M. 1987. Viscous fingering in
2655-2668. porous media, Annu. Rev. Fluid Mech. 19,
Hanon, D. & Boon, J.P. 1997. Diffusion 271-311.
and correlations in a lattice gas Hopcroft, J.E. & Ullman, J.D. 1979.
automata, Phys. Rev. E 56(6), 6331-6339. Introduction to Automata Theory,
Hardy, J. & Pomeau, Y. 1972. Languages and Computation.
Thermodynamics and hydrodynamics for Addison-Wesley, New York.
a modeled fluid, J. Math. Phys. 13, Huang? K 1963 statistical Mechanics.
1042-1051. w i l e y ? N e w York.
Hardy, J., Pomeau, Y & de Pazzis, O. 1973. Huerre> p & Monkewitz? pA 1990 Local
Time evolution of a two-dimensional a n d global instabilities in spatially
model system: Invariant states and time developing flows, Annu. Rev. Fluid Mech.
correlation functions, J. Math. Phys. 14, 22 473-537
1746-1759.
XT , _ . _ . _ n _ __ 4 r _ d'Humieres, D., Lallemand, P. & Frisch, U.
Hardy, J., de Pazzis, O. & Pomeau, Y 1976. i n o , T ' A. . . ,, ~
^/' . t . ^ t .' . 1986. Lattice gas models for 3-D
Molecular dynamics
J of a classical lattice , , , ^ , r -
_ • j • hydrodynamics, Europhys. Lett. 2,
gas: Transport properties and time 291-297
correlation functions, Phys. Rev. A 13,
1949-1961. d'Humieres, D., Qian, YH. & Lallemand, P.
u , T p „ , , 1 f t _ Ayr. . 1989. Invariants in lattice gas models,
Hardy, J. & Pomeau, Y 1977. Microscopic . n 3. . , TTr .6,
J 1f • a • 4. in Proceedings of the Workshop on
model tor viscous now in two ,
,. . D, D A iiC n^r. n~c Discrete Kinetic Theory,y Lattice Gas
dimensions, Phys. Rev. A 16, 720-726. , .
Dynamics and Foundations oj
Henon, M. 1987a. Isometric collision rules Hydrodynamics, Torino (Italy), September
for the 4-D FCHC lattice gas, Complex 20-24, 1988, 102-113, ed. R. Monaco.
Systems 1, 475-494. World Scientific, Singapore.
Henon, M. 1987b. Viscosity of a lattice gas, d'Humieres, D., Lallemand, P. & Qian, YH.
Complex Systems 1, 763-789. 1990 Finding the linear invariants of
Henon, M. 1989. On the relation between lattice gases, in Computational Physics
and
lattice gases and cellular automata, Cellular Automata, eds. A. Pires,
in Proceedings of the Workshop on D.P. Landau & H. Herrman. World
Discrete Kinetic Theory, Lattice Gas Scientific, Singapore.
Dynamics and Foundations of Kinchin, A.I. 1957. Mathematical
Hydrodynamics, Torino (Italy), Foundations of Information Theory.
September 20-24, 1988, 160-161, ed. Dover? N e w York
R. Monaco. World Scientific, Singapore.
Konig, M., Eisenlohr, H. and Eckelmann, H.
Henon, M. 1992. Implementation of the 1 9 9 Q J h e flne s t r u c t u r e i n ^
FCHC lattice gas model on the Strouhal-Reynolds number relationship
connection machine, in Proceedings of rf ^ , a m i n a r w a k e rf a d r c u l a r
the NATO Advanced Research Workshop ^ ^ phys FMds A 1(% 1607_1614
on Lattice Gas Automata Theory,
Implementation, and Simulation, Nice Kubo. R. 1958. Some aspects of the
(France), June 25-28, 1991, J. Stat. Phys. statistical mechanical theory of
68(3/4), 353-377, ed. J.P. Boon. irreversible processes, in Lectures in
J TrrA/rAPX7 11T-, m™ Theoretical Physics, vol. 1, 120-140, eds.
van der Hoef, M.A. & Frenkel, D. 1990. __,„
W E _ . . o\ „ _ ' _Tri
. ., r 1 1 • - Bnttin & L.G. Dunham. Wiley,
Lone-time
& tails ol the velocity ^T ^T t
, . c . . / , New York,
autocorrelation function in two- and
three-dimensional lattice-gas cellular Lamb, H. 1945. Hydrodynamics, 6 th edition,
automata, Phys. Rev. A 41, 4211-42S5. Dover, New York.
278 References
Landau, L.D. & Lifschitz, E.M. 1987. Fluid states. IV. Far-from-equilibrium
Mechanics, 2nd edition, vol. 6, in 'A stationary states and the regression of
Course of Theoretical Physics'. Pergamon fluctuations, Phys. Rev. A 20(6),
Press, Oxford. 2533-2546.
Lutsko, J.F., Boon, J.P. & Somers, J.A. Resibois, P. & de Leener, M. 1977.
1992. Lattice gas automata simulations Classical Kinetic Theory of Fluids. John
of viscous fingering in a porous medium, Wiley, New York.
in Numerical Methods for the Simulation ~. , TT irkO/) ^ r ,, „, ,
J
. Risken, H. 1984. The Fokker-Planck
of Multi-Phase and Complex Flow, F .. c . ,, ,
;,, n f J ™ , * , , , ^ Equation. Springer- Verlag, DBerlin.
r
124-135, ed. T.M.M. Verheggen.
Springer-Verlag, Berlin. Rivet, J.P. 1987. Simulation d'ecoulements
,,., A ino. _ .L1 _. ., _. tri-dimensionnels par la methode des gaz
Majda, A. 1984. Compressible Fluid Flow , . , . _ „&
,' /•„ • • sur reseaux: premiers resultats, C. .R.
ana Systems of Conservation Laws in , . . , . TT _A
,o T^ • 11 o • ^r i Acad. Sci. Pans, sene II 305, 151-156.
Several Space Variables. Springer-Verlag,
Berlin. Rivet, J.P. 1987. Green-Kubo formalism
AT uiAyrour r>T J-* 1 nni for lattice gas hydrodynamics and
Mareschal, M. & Hohan, B.L., editors 1992. ^r _ ° ;, /
,,. . c . , .. r ~ , Monte-Carlo evaluation of shear
Microscopic Simulations of Complex .
TJ J J • r>i r»i viscosities, Complex Systems 1, 839-851.
Hydrodynamic Phenomena. Plenum, > r * •>
New-York. Rivet, J.P. 1988. Hydrodynamique par la
m 6 t h o d e des gaz sur r6seaux P h D
Molvig, K, Donis, P., Mffler.IL, ' '
Myczkowski, J. & Vichniac, G. 1989. Umversite de Nice.
Multi-species lattice-gas automata for Rivet, J.P. 1993. Spontaneous
realistic fluid dynamics, In Cellular symmetry-breaking in the 3-D wake of a
Automata and Modeling of Complex long cylinder simulated by the lattice gas
Systems, 206-231, eds. P. Manneville, R. method and drag coefficient
Bideaux, N. Boccara & G. Vichniac. measurements, Appl. Sci. Res. 51,
Springer-Verlag, Berlin. 123-126.
Mori, H. 1965a. Transport, collective Rothman, D.H. & Keller, J.M. 1988.
motion, and Brownian motion, Prog. Immiscible cellular automaton fluids,
Theor. Phys. 33, 423-^55. J. Stat. Phys. 52, 1119-1127.
Mori, H. 1965b. A continued fraction Rothman, D.H. & Zaleski, S. 1994.
representation of the time-correlation Lattice-gas models of phase separation:
function, Prog. Theor. Phys. 34, 399-416. interfaces, phase transitions, and
Naitoh, T, Ernst, M.H., van der Hoef, M.A. multi-phase flow, Rev. Mod. Phys. 66,
& Frenkel, D. 1991. Extended mode 1417-1479.
coupling and simulations in Rothman, D.H. & Zalseki, S. 1977. Lattice
cellular-automata fluids, Phys. Rev. A 44, Gas Cellular Automata, Cambridge
2484^2494. University Press, Cambridge.
McNamara, G.R. 1990. Diffusion in a Schmitz, R. & Dufty, J.W. 1990. Positivity
lattice gas automaton, Europhys. Lett. 12, o f transport coefficients for cellular
329
- automata, Phys. Rev. A 41(8), 4294-4297.
McQuarrie, D.A. 1976. Statistical Shannon, C.E. & Weaver, W. 1969. The
Mechanics. Harper and Row, New York. Mathematical Theory of Communication,
Noullez, A. 1990. Automates de Gaz sur fourth edition. The University of Illinois
Reseau: Aspects Theoriques et Press, Urbana.
Simulations, Ph.D., Universite Libre de Suarez? A & Boon> j p 1997a Non_linear
Succi, S. 2000. The Lattice Boltzmann Reynolds numbers, J. Fluid Mech. 206,
Equation for Fluid Dynamics and Beyond. 579-627.
Oxford University Press, Oxford (U.K.) Wolfram, S. 1986. Theory and Applications
(to be published). of CeiMar Automata. World Scientific,
Succi, S. Benzi, R. & Higuera, F. 1991. The Singapore.
lattice Boltzmann equation: A new tool Wolfram, S. 1986. Cellular automaton
for computational fluid mechanics, fluids x. b a s i c t h j S m Ph 45
Physica D 47 (1/2), 219-230. 471-525.
Tribel, O. & Boon, J.P. 1995. Lattice gas Zanetti, G. 1989. Hydrodynamics of lattice
automata with interaction potential, automata, Phys. Rev. A 40(3),
J. Stat. Phys. 81, 361-377. 1539-1548
Tribel, O. & Boon, J.P. 1997. Entropy and Zwanzig, R. 1960. Statistical mechanics of
correlations in lattice gas automata irreversibility, in Lectures on Theoretical
without detailed balance, Int. J. Mod. physics> w / J > 1 0 6 - 1 4 1 , eds. WE. Brittin,
Phys. C 8(4), 641-652. B w D o w n s & j . Downs. Wiley, New
Tung, W.K. 1985. Group Theory in Physics. York.
World Scientific, Singapore. Zwanzig, R. 1965. Time-correlation
Williamson, C.H.K. 1989. Oblique and functions and transport coefficients in
parallel modes of vortex shedding in the statistical mechanics, Ann. Rev. Phys.
wake of a circular cylinder at low Chem. 16, 67-102.
Author index
281
282 Author index
Ernst, M.H., 98, 139, 168, 171, 183, 184, 186, Kubo, R., 196
194, 196-198, 238, 240, 243, 249
Ladd, A.J.C., 241
Feder, J., 245 Lallemand, P., xvii, 15, 22, 34, 38, 42, 48, 53, 59,
Feynman, R., xiii 74, 133, 236, 237, 240-243, 246, 247
Flekkoy, E.G., 243, 245-247 Landau, L.D., 112
Fleury, P.A., xvii, 170 Lavallee, P., 240
Forster, D., 107 Lawniczak, A., 242
Frenkel, D., 198, 238, 241 Lee, Y.C., 243, 244
Frisch, U., 15, 34, 39, 59, 74, 133, 183, 235-237, Leener (de), M., 157
239, 247-249 Leeuwen (van), J.M.J., 197
Fu, C , 243 Levermore, D., 30, 59, 240
Lifschitz, E.M., 112
Gatignol, R., 83, 234, 256 Longo, E., 234
Ginzbourg, I., 241, 245 Luo, L.S., 243
Godreche, C , 249 Lutsko, J.L., 243
Groos, E.R, 246
Grosfils, P., xvii, 34, 53, 98, 238, 243 Maillot, B., 244
Gunstensen, A.K., 242, 245 Majda, A., 203
Gutman, S., 243 Mareschal, M, 154
Margolus, N., 208, 246, 247
Hanon, D., 175, 177, 239 Markus, M., 244
Hardy, J., 34, 35, 80, 235, 237 Martinez, D.O., 245
Hasslacher, B., 15, 34, 39, 74, 133, 235, 237 Matthaeus, W.H., 241, 244, 245, 249
Hatori, T., 244 Maxwell, J.C., 234
Hauge, E.H., 197 McNamara, G.R., 175, 247
Helmholtz (von), H., 222 McQuarrie, D.A., 103
Henon, M., 9, 60, 63, 74, 236, 237, 239, 247, 248 Mermin, N.D., 2, 5, 7, 11
Hess, B., 244 Meyer, D., 239
Higuera, R, 245 Miller, R., 104
Hillis, W.D., xiii Molvig, K, 104
Hoef (van der), M.A., 198, 238, 241 Monaco, R., 234, 248
Holian, B.L., 154 Monkewitz, P.A., 227
Homsy, G.M., 231 Montgomery, D.C., 244, 245
Hopcroft, J.E., 2, 235 Mora, P., 244, 245
Huang, K., 69, 107 Mori, H., 140, 184
Huerre, P., 227 Myczkowski, J., 104
Humieres (d'), D., 15, 22, 30, 34, 38, 42, 48, 59,
74, 133, 208, 236, 237, 240-242, 245-247, Naitoh, T., 198, 238
249 Neuman (von), J., 235
Noullez, A , 175, 237, 240
Kadanoff, L.R, 247
Kapral, R., 242 Olson, J.F., 242
Karman (von), T, 218, 226 Orszag, S.A, 113, 195, 246
Karman (von), T., xiv Oxaal, U., 245
Keller, J.M., 242
Kelvin, (Lord), 222 Papoulis, A., 106
Kinchin, A.I., 74 Pazzis (de), O., 34, 35, 80, 235, 237
Kirkpatrick, T.R., 186 Platkowski, T, 234
Kohring, G.A, 243 Pomeau, Y., 15, 34, 35, 39, 74, 80, 133, 235, 237,
Krook, M.K., 246 240-242, 247
Author index 283
285
286 Subject index
microdynamic global, 12
equation, 10, 15 local, 12
formalism, 10 single-node, 12
microscopic density phonons, 243
per channel, 13 plasmas, 244
per node, 13 point group, 5, 23
per unit volume, 13 porous medium, 231, 243
microstate, 69 position vector, 4, 11
minimal model, 6 post-collision state, 4
mixture (of fluids), 241 power spectrum, 178
mode-coupling theory, 196, 238 pre-collision state, 4, 72
molecular chaos hypothesis, 71 pressure, 120, 200
molecular dynamics, xiii, 153 kinetic, 103
moment propagation method, 198 thermodynamic, 103
momentum, 12 primitive cell, 13
density per node, 95 probability distribution function, 69
multi-scale, 110, 188 projection operator, 140, 184
multinomial expansion, 28 projectors, 272
multiple channel, 5 propagation
operator, 16
Navier-Stokes equations, xiii, 128, 203 phase, 4, 16
neutron scattering spectroscopy, 153
Newtonian fluid, 34, 112, 204, 207 quasi-lattice, 8
node, 2, 4, 11
node-wise storage, 211 Reseau d'Automates Programmables (RAP),
non-Galilean factor, 92, 121, 200 208, 247
for energy, 121 random force, 147
for momentum, 121 Rayleigh-Brillouin spectrum, 169
non-thermal models, 89, 199 reaction-diffusion, 241
normalization, 19, 69 reflection
bounce-back, 31
observable, 12 diffusive, 32
energy, 14 operator, 30
generalized, 14, 69 specular, 32
homogeneous, 15 regular lattice, 7
local, 15 rest particles, 5
mass, 14 Reynolds number, 204
momentum, 14 Reynolds's similarity law, xv
ordinary, 70
particle number, 14 Saffman-Taylor instability, 231
obstacle, 30, 240 scale
operator macroscopic, 110, 112
collision, 16 microscopic, 109, 110
linearized, 138 multi-, 110, 188
evolution, 16 separation, 109
projection, 140, 184 space- and time-, 109
propagation, 16 seismic waves, 243
self-dual, 21
passive scalar, 67, 225 semi-detailed balance, 20, 82
phase separation, 242 shear viscosity, 127, 202
phase space, 12 shear wave, 202
Subject index 289