Lattice gas automata for fluid mechanics

Physica 140A (1986) 326-335 North-Holland, Amsterdam

LATTICE GAS AUTOMATA FOR FLUID MECHANICS D. D'HUMII~RES and P. LALLEMAND CNRS, Laboratoire de Physique de l'Ecole Normale Sup~rieure, 24 rue Lhomond, 75231 Paris Cedex 05, France

A lattice gas is the representation of a gas by its restriction on the nodes of a regular lattice for discrete time steps. It was recently shown by Frisch, Hasslacher and Pomeau that such very simple models lead to the incompressible Navier-Stokes equation provided the lattice has enough symmetry and the local rules for collisions between particles obey the usual conservation laws of classical mechanics. We present here recent results of numerical simulations to illustrate the power of this new approach to fluid mechanics which may give new tools for numerical studies and build a bridge between cellular automata theory and complex physical problems.

1. Introduction

Fluid mechanics is one of the most important fields for practical applications: airplanes, weather prediction, etc., and one of the most challenging for theoretical studies. Despite a lot of theoretical work, a large part of the present knowledge comes from experiments. Fluid dynamics can be studied from two points of view. The first one consists of describing the fluid from a macroscopic point of view. In this case the equations of motion are derived as the most general equations for mass, momentum and energy conservation which satisfy the usual space invariance (translation, rotation and Galilean invariances). These equations have been known for more than one hundred years as the continuity equation,

Otp + div (pu) = O,

(la)

the Navier-Stokes equation,

POtU + p ( u . grad)u = - g r a d P + ~lAu + ~ grad (div u),

(ab)

and the heat equation; where p, u and P are the local fluid density, velocity and pressure and 7/and ~ are the shear and bulk viscosities1). Despite their apparent simplicity, the non-linear terms make these equations very difficult to use: the analytical solution is known only in a few simple cases. Generally, eqs. (1) are used under their dimensionless form in order to compare the results obtained for different fluids or geometries. Then, the isothermal flows depend on two im0378-4371/86/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

LATTICEGASAUTOMATAFOR FLUIDMECHANICS

327

portant dimensionless parameters: the Reynolds number: Re = puL/~ (where u and L are the typical velocity and length scales of the problem), which gives the influence of the non-linear terms in eqs. (1), non-stationary and turbulent flows occurring for high Reynolds numbers, and the Mach number: M = u/c (where c is the sound velocity), which gives the influence of compressible effects, low Mach numbers implying incompressible flows. The second approach, the molecular dynamics, consists of a statistical description of the fluid from its microscopic properties. The transport coefficients (sound velocity, viscosities. . . . ) are derived from the molecular properties such as the interparticle potential. However, the full description of the microscopic world allows only to simulating the behavior of a small number of particles. A first attempt to get correct microscopic properties from a simplified description of the microscopic world was done by Broadwell and Gatignol 2) when they studied a gas in which the possible molecular velocities are restricted to a finite set. Then, Hardy, de Pazzis and Pomeau (HPP) 3) tried to use the simplest possible model: a gas on a square lattice with very simple collision and propagation rules. However, this model was too simple to give realistic macroscopic behavior. More recently, the fast development of integrated circuits and computers has triggered a new interest in cellular automata, which are regular networks of simple finite state machines with local connections between ceils. The most complete study of classes of behavior in cellular automata has used analysis of machine simulations of such systems in one spatial dimension4) and it was shown that they are capable of extremely complex behavior. Conjectures have been made that there exist cellular automata capable of simulating the solution space of partial differential equations of physical interestS). No non-trivial example of this conjecture was known until Frisch, Hasslacher and Pomeau (FHP) recently found a way to cure most of the pathology of the HPP model and were able to prove that this new lattice gas evolves according the Navier-Stokes equation in the limit of large lattice size and low velocity6). We present here recent results of computer simulations of the FHP models and some of its variants. First, we give a brief description of the model. Next we present the general procedure for the computer simulations; then we describe the linear hydrodynamic regime, measure transport coefficients and simulate two dimensional flows at moderate Reynolds number: a flow at the inlet of a duct and a v o n Karman street behind a flat plate.

2. The models

We consider particles moving on a triangular lattice with unit velocity c i in direction i between a node and one of its six neighbors (i = 1. . . . . 6). At each

328

D. D'HUMII~RESAND P. LALLEMAND

time step each particle stays on a node and interacts on it with the other incoming particles according to collision laws assumed to conserve the number of particles and the total momentum on the node. Additionally, there is an exclusion principle such that no two particles with the same velocity may occupy the same node at the same time. Then the particles propagate according their new velocity. We have also used a variant of this model with "rest particles" which allows an additional particle with zero velocity at each site (i = 0 for notational purposes); these particles may be considered to have an internal energy to satisfy energy conservation. These lattice gas models can be implemented on a two dimensional cellular automata with a triangular lattice topology, where each cell has six neighbors and a seven bit state. The collision laws of the lattice gas give the transition rules for the internal states of each cell and the propagation rules of the particles give the correspondence between the internal states and the cell inputs and outputs. The original FHP model uses only six particles and the collision rules given in fig. l a (five possible collisions, model I), which is the minimal set of rules to prevent spurious mass or momentum conservation. In fact, most of the computer simulations were done with the model with resting particles and the additional rules shown in fig. lb and those of fig. la with rest particles (twenty-two possible collisions, model II). More recently we used a model with rest particles and all the possible collisions conserving mass and momentum (seventy-six possible collisions, model III). The macroscopic quantities, density p and momentum pu, are related to the macroscopic local average populations N~ of the Boolean states by 7) p = ENi,

(2)

pu = E N i e i .

i

i

before

',,

after

',, ar /

/

a

_/ X

/ --X

b

X

X I

.

~

Fig, 1. Schematic representation of some collision rules; (a) no rest particles; (b) additional rules for

rest particles.

LATTICE GAS AUTOMATAFOR FLUID MECHANICS

329

As the particles satisfy an exclusion principle, at thermodynamic equilibrium the N~ are given by Fermi-Dirac distributions and can be expanded around small U by 6'8'9) N~ = d{1 +

(p/3d)cs"

(3)

u},

where the reduced density d is the average density per direction (without rest particles p = 6d, with rest particles p = 7d). Using a Chapman-Enskog expansion up to second order terms in velocity and gradients one gets

~,(ou,~) + Eaa[g(o)ou~ut~ ] = -~,~P + nAu,, + #

~.(div u),

where the Greek subscripts refer to the x and y components and given by

g(P)

= (P - 3 ) / ( 0 - 6),

g(P) and

(4) P are

(5a)

P = 0/2,

for model I without resting particles and by

g(p)

= 7 ( 2 p - 7 ) / ( p - 7),

P =

3p/7,

(5b)

for the models II and III with rest particles. The term g(p) comes from the F e r m i - D i r a c distribution and restricts the use of lattice gas models to incompressible flows. In this case, g(p) is a constant which may be absorbed in a rescaled velocity, the resulting Reynolds number being Re = pg(p)uL/~. The values of the viscosities can be computed in the Boltzmann approximation from the collision operator1°). As noticed by H6nonn ) an additional negative term comes from the discrete nature of the lattice gas and decreases the effective viscosities. Since the definition of the density on a triangular lattice is slightly different from the usual one and since the reduced density is the natural parameter for the computer simulations, in what follows we will use reduced viscosities T/r = dTI/p and ~r = d~/p. A straightforward extension of the results of ref. 10 gives

7/r

1 12(1 -- d "3)

d 8'

~r = 0,

(6a)

for model I,

7/r

1

d

28(1 - d)3(1 - 4 d / 7 )

8'

~r

1

d

98(1 -- d ) '

28'

(6b)

330

D. D'HUMII~RESAND P. LALLEMAND

for the model II and

"Or

=

1

d

28(1 - d)(1 - 8d(1 - d ) / 7 )

8'

(6c) ~r

=

1

d

98(1 - d)(1 - 2d(1 - d))

28'

for the model III (units of time, length and mass are the time step, the link length and the particle mass, respectively).

3. Computer simulations The present work was done by simulation of the lattice gas models on an FPS-164 using lattices of order 106 nodes with a typical speed of 106 updates per second12-15). The evolution of the lattice gas is computed according to a parallel iteration in two steps. During the first one, the post-collision state is computed using either a look-up table (in this case the states are coded with eight bits, seven for the particles and one to choose between the different possible head-on collisions, and eight nodes are packed in each of the 64 bits words of the FPS) or the combination of Boolean operators giving the collision rules (in this case 64 nodes are packed in a word and seven words are used to code the different particles, the choice between the head-on collisions is made randomly at each time step). During the second step the bits coding the different particles are moved from one node to the next one. The lattice boundary conditions can be either periodic in both directions or of the " w i n d tunnel" type: particles are generated and absorbed at the lattice edges up and downstream according to eq. (3), u being the velocity of the wind tunnel. The boundary conditions on obstacles are "bounce-back" reflections: the particle velocity after the collision is the opposite of the incoming velocity. Initial flows were generated by a Monte Carlo procedure with average population N i related to the local velocity by eq. (3). Macroscopic quantities are obtained by averaging the N/ according to eqs. (2) over rectangular regions with shapes and sizes adapted to the studied flow. 3.1. Linear hydrodynamics The velocity of sound c, the shear viscosity 7/ and the bulk viscosity ~ of the lattice gas models described above have been measured, using the relaxation of an initial periodic perturbation (ult + U_L)COS(k " r) of the velocity fieldlr), k is

LATTICE GAS A U T O M A T A F O R F L U I D MECHANICS

'Ij -I

I

I

I

I

I

I

I

I

I

I

I

I

256

331

[

I

t

512

Fig. 2. Time evolution of u ± (dashed line), Ull (solid line) and p (dot-dash line) normalized by the initial values of the velocity components.

the wave vector of the perturbation, and ull and u_L are the velocity components parallel and perpendicular to the wave vector. The relaxation in time of the velocity u(r, t) and of the density perturbation 8p(r, t) are given by

u(r, t) = (ullcos(~ot + cp) exp ( - k 2 t ( , / + ~ ) / 2 p ) + U_L e x p ( - k 2 t ~ / p ) ) c o s ( k So(r, t) = (PUll/c) sin(~ot)exp(-k2t(~

• r), + ~ ) 1 2 p ) s i n ( k . r),

(7a) (7b)

with to = ck, tanq~ = k(,; + ~)/2pc. Starting from the initial conditions, at each time step the momentum and the density are averaged along lines perpendicular to the wave vector. The result is Fourier transformed to get the component of the momentum and density corresponding to k. Typical relaxation curves are given in fig. 2. From these curves, c, ~1 and ~ are measured by least squares fits of eqs. (7) to the time evolution of u± (k), ull(k ) and p(k). The measured sound velocities are isotropic and agree with theoretical values 1 / v ~ - for model I and V/3-/7 for models II and III. The measured values of the viscosities are summarized in fig. 3, along with the theoretical curves computed from eqs. (6). These measurements were obtained on 2562 lattices with periodic boundary conditions and for wavelength between 16 and 32 nodes; the size of the symbols corresponds roughly to the error bars. The measured viscosities agree with theoretical prediction for models II and III and confirm the presence of the negative lattice viscosities. Without rest particles the experimental values of ,; are

332

D. D'HUMII~RES AND P. LALLEMAND 0 . 2

.

.

.

.

,

.

.

.

.

,

.

.

.

.

,

•

0"15

•

•

•

.

.

.

.

,

•

.

.

.

o

•

o 0.1 ._~

,,"

0.05

/"

a

x

0

.

0

.

.

.

.

-"f~

'

0.1

,

,

,

0.2 Density

,

,

,

J

,

0.3 per link

•

•

'

0.4

.

.

.

.

0.5

Fig. 3. Theoretical shear (solid lines) and bulk (dashed lines) reduced viscosities as a function of the

reduced density, compared with the results of numerical simulations for different lattice gas models: original FHP model (solid squares and circles) model with rest particles and limited collision rules (open squares and triangles) and model with rest particles and all the possible collisions (oblique and triangular crosses).

too high a n d ~ is found negative. Thus, the presence of rest particles apparently i m p r o v e s the qualitative behavior of the lattice gas while decreasing the viscosity significantly. Fig. 4 shows the achievable Reynolds n u m b e r for the different m o d e l s for uL = 1; the m a x i m u m value is obtained with model I I I and is six times larger than the m a x i m u m value for the original F H P model. These results s h o w that R e y n o l d s numbers depend not only on the n u m b e r of nodes as noted in refs. 6 a n d 17 but also on the lattice gas models. Further improvements on

__d

~D

Od

0

0

0.1

0.2 Density

0.3 per link

0.4

0.5

Fig. 4. Available Reynolds number as a function of the reduced density (solid lines: theory; symbols: measured values) for the different models with the same symbols as in fig. 3.

LATTICE GAS AUTOMATA FOR F L U I D MECHANICS

333

them may increase the maximum achievable Reynolds numbers which is presently of order of 10 3. 3.2. Non-linear flow simulations The experimental points of fig. 4 are computed from the measured viscosity values and the theoretical values for g(p). In order to test the non-linear terms in the lattice gas dynamics, we studied the flow at the inlet of a 2-D duct using a 512 x 3072 latticeXS). The average density and velocity of the gas are 1.54 (d = 0.22) and 0.30 and the lattice gas model is model II. In fig. 5, we have compared the momentum obtained for the lattice gas flow to the values computed according to Slichting 18) without any adjustable parameters, using the viscosity value from the relaxation measurements and the theoretical value for g(p). The agreement between the two methods shows that the lattice gas dynamics is correctly described by eq. (4). In fig. 6 is shown a v o n Karman street behind a "bounce-back" fiat plate placed normal to a mean flow of velocity 0.51 in a 1024 × 512 wind-tunne114). The lattice gas model is model II, the length L of the plate is 120 and the density p = 1.4 (d = 0.2). The corresponding Reynolds number is 90. The flow is

Y,/a 2

i

v~ o o

o,~

~2 (a)

0,3

i 0,4

0

o3

0,2

0.3

0.4

(b)

Fig. 5. Velocity component across a channel, relative distance from the inlet (a) 0.525, (b) 5.74. The dots are obtained by the lattice gas simulations, the solid lines are calculated using the Slichting method.

334

D. D'HUMII~RES AND P. LALLEMAND

Fig. 6. Von Karman street behind a flat plate at 5000 time steps (Re - 90) with "wind-tunnel" boundary conditions; the mean flow has been subtracted.

time-dependent with a period T = 1200 time steps and a Strouhal number L / u T of approximately 0.2 close to the experimental value 0.16.

4. Conclusion Quantitative agreement between theory and simulation has been demonstrated in both the linear and non-linear regimes, for moderate Reynolds numbers. It should be noted that the introduction of obstacles into the flow is particularly simple and represents, for the cases studied, a computational overhead of a few percent. The present models are limited to Reynolds numbers of order of 103 and incompressible flows. Moreover, since local equilibrium is a function of p and pu only, such models cannot simulate thermal phenomena. However, more complicated models derived from the FHP model 8,9) may overcome most of these limitations in the near future. In this case, lattice gas simulations will be a new tool for experimental work in hydrodynamics, with the main advantage of being inherently stable.

Acknowledgements We thank U. Frisch, B. Hasslacher, M. H~non, Y. Pomeau and J.P. Rivet for helpful discussions and GRECO 70 "Experimentation num~rique" for support.

LATTICE GAS AUTOMATA FOR FLUID MECHANICS

335

References 1) L.D. Landau and F.M. Lifshitz, Fluid Mechanics (Pergamon, London, 1959); the Navier-Stokes eq. (lb) is given for two dimensional flows, for D dimensions ~ must be replaced by ~ + (D -

2)n/D. 2) I.E. Broadwell, Phys. of Fluids 7 (1964) 1243; R. Gatignol. Th~orie cin6tique des gaz h r6partition discr&e de vitesse, Lecture Notes in Physics (Springer, Berlin, 1975). 3) J. Hardy, and Y. Pomean, J. Math. Phys. 13 (1972) 1042; J. Hardy, O. de Pazzis and Y. Pomeau, J. Math. Phys. 14 (1973) 1746; Phys. Rev. A 13 (1976) 1949. 4) S. Wolfram, Rev. Mod. Phys. 55 (1984) 96. 5) Y. Pomeau, J. Phys. A 17 (1984) 415; G. Vichniac, Physica I@D (1984) 96 and refs. therein; N.H. Packard and S. Wolfram, Two dimensional celhilar automata, I.A.S. preprint, (Sept. 1984). 6) U. Frisch, B. Hasslacher and Y. Pomeau, Phys. Rev. Lett. 56 (1986) 1505. 7) Using these formulae, the density is defined as the number of particles per node; this is v:3/2 of the density which is usually defined as the number of particles per square unit. 8) D. d'Humi~res, P. Lallemand and U. Frisch, Europhys. Lett. 2 (1986) 291. 9) S. Wolfram, Cellular automaton fluids 1: basic theory, Institute for Advanced Studies preprint (1986). 10) J.P. Rivet and U. Frisch, C.R. Acad. Sci. Paris II 302 (1986) 267. 11) M. H6non, Calcul de la viscosit~ dans le r6seau triangulaire and Viscosit~ d'un r6seau, preprints Observatoire de Nice, B.P. 139, 06003 Nice Cedex, France; U. Frisch and J.P. Rivet, Lattice gas hydrodynamics: Green-Kubo formula, submitted to C.R. Acad. Sci. Paris II (1986). 12) D. d'Humi6res, Y. Pomean and P. Lallemand, C.R. Acad. Sci. Paris II 3@1 (1985) 1391. 13) D. d'Humi6res, P. Lallemand and T. Shimomura, Computer simulations of lattice gas hydrodynamics, submitted to Phys. Rev. Lett. (1985). 14) D. d'Humi6res, Y. Pomean and P. Lallemand, Two dimensional hydrodynamics calculations with a lattice gas, Innovative Numerical Methods in Engineering 241 (A Computational Mechanics Publication, Springer-Vedag, Berlin, 1986). 15) D. d'Humi6res and P. Lallemand, C.R. Acad. Sci. Paris II 302 (1985) 983. 16) R.D. Mountain, Rev. Mod. Phys. 38 (1966) 205. 17) S.A. Orszag and V. Yakhot, Phys. Rev. Lett. 56 (1986) 1691. 18) H. Slichting, Z.A.M.M. 14 (1934) 368; Boundary Layer Theory (Pergamon, London, 1955).