lanxicy.com

第一范文网 文档专家

第一范文网 文档专家

Journal of Statistical Physics, Vol. 45, Nos. 3/4, 1986

Cellular Automaton Fluids 1: Basic Theory

Stephen Wolfram ~'2

Received March, 1986," revision received August, 1986

Continuum equations are derived for the large-scale behavior of a class of cellular automaton models for fluids. The cellular automata are discrete analogues of molecular dynamics, in which particles with discrete velocities populate the links of a fixed array of sites. Kinetic equations for microscopic particle distributions are constructed. Hydrodynamic equations are then derived using the Chapman-Enskog expansion. Slightly modified Navier-Stokes equations are obtained in two and three dimensions with certain lattices. Viscosities and other transport coefficients are calculated using the Boltzmann transport equation approximation. Some corrections to the equations of motion for cellular automaton fluids beyond the Navier-Stokes order are given. KEY W O R D S : Cellular automata; derivation of hydrodynamics; molecular dynamics; kinetic theory; Navier Stokes equations.

1. I N T R O D U C T I O N

Cellular automata (e.g., Refs. 1 and 2) are arrays of discrete cells with discrete values. Yet sufficiently large cellular automata often show seemingly continuous macroscopic behavior (e.g., Refs. 1 and 3). They can thus potentially serve as models for continuum systems, such as fluids. Their underlying discreteness, however, makes them particularly suitable for digital computer simulation and for certain forms of mathematical analysis. On a microscopic level, physical fluids also consist of discrete particles. But on a large scale, they, too, seem continuous, and can be described by the partial differential equations of hydrodynamics (e.g., Ref. 4). The form of these equations is in fact quite insensitive to microscopic details.

~The Institute for Advanced Study, Princeton, New Jersey and Thinking Machines Corporation, 245 First Street, Cambridge, Massachusetts 02142. 2 Present address: Center for Complex Systems Research, and Departments of Physics, Mathematics and Computer Science, University of Illinois, 508 South Sixth Street, Champaign, Illinois 61820.

471

0022-4715/86/1100-0471505.00/0 ~2) 1986PlenumPublishingCorporation

472

Wolfram

Changes in molecular interaction laws can affect parameters such as viscosity, but do not alter the basic form of the macroscopic equations. As a result, the overall behavior of fluids can be found without accurately reproducing the details of microscopic molecular dynamics. This paper is the first in a series which considers models of fluids based on cellular automata whose microscopic rules give discrete approximations to molecular dynamics. 3 The paper uses methods from kinetic theory to show that the macroscopic behavior of certain cellular automata corresponds to the standard Navier Stokes equations for fluid flow. The next paper in the series ~6) describes computer experiments on such cellular automata, including simulations of hydrodynamic phenomena. Figure 1 shows an example of the structure of a cellular automaton fluid model. Cells in an array are connected by links carrying a bounded number of discrete "particles." The particles move in steps and "scatter" according to a fixed set of deterministic rules. In most cases, the rules are chosen so that quantities such as particle number and momentum are conserved in each collision. Macroscopic variations of such conserved quantities can then be described by continuum equations.

3 This work has m a n y precursors. A discrete model of exactly the kind considered here was discussed in Ref. 6. A version on a hexagonal lattice was introduced in Ref. 7, and further studied in Refs. 8, 9. Related models in which particles have a discrete set of possible velocities, but can have continuously variable positions and densities, were considered much earlier.~o ~4) Detailed derivations of hydrodynamic behavior do not, however, appear to have been given even in these cases (see, however, e.g., Ref. 15).

Fig. 1. Two successive microscopic configurations in the typical cellular a u t o m a t o n fluid model discussed in Section 2. Each arrow represents a discrete "particle" on a link of the hexagonal grid. C o n t i n u u m behavior is obtained from averages over large numbers of particles.

Cellular Automaton Fluids

473

Particle configurations on a microscopic scale are rapidly randomized by collisions, so that a local equilibrium is attained, described by a few statistical average quantities. (The details of this process will be discussed in a later paper.) A master equation can then be constructed to describe the evolution of average particle densities as a result of motion and collisions. Assuming slow variations with position and time, one can then write these particle densities as an expansion in terms of macroscopic quantities such as momentum density. The evolution of these quantities is determined by the original master equation. To the appropriate order in the expansion, certain cellular automaton models yield exactly the usual Navier-Stokes equations for hydrodynamics. The form of such macroscopic equations is in fact largely determined simply by symmetry properties of the underlying cellular automaton. Thus, for example, the structure of the nonlinear and viscous terms in the Navier Stokes equations depends on the possible rank three and four tensots allowed by the symmetry of the cellular automaton array. In two dimensions, a square lattice of particle velocities gives anisotropic forms for these terms. (6) A hexagonal lattice, however, has sufficient symmetry to ensure isotropyJ 7/ In three dimensions, icosahedral symmetry would guarantee isotropy, but no crystallographic lattice with such a high degree of symmetry exists. Various structures involving links beyond nearest neighbors on the lattice can instead be used. Although the overall form of the macroscopic equations can be established by quite general arguments, the specific coefficients which appear in them depend on details of the underlying model. In most cases, such transport coefficients are found from explicit simulations. But, by using a Boltzmann approximation to the master equation, it is possible to obtain some exact results for such coefficients, potentially valid in the lowdensity limit. This paper is organized as follows. Section 2 describes the derivation of kinetic and hydrodynamic equations for a particular sample cellular automaton fluid model. Section 3 generalizes these results and discusses the basic symmetry conditions necessary to obtain standard hydrodynamic behavior. Section 4 then uses the Boltzmann equation approximation to investigate microscopic behavior and obtain results for transport coefficients. Section 5 discusses a few extensions of the model. The Appendix gives an SMP program/17) used to find macroscopic equations for cellular automaton fluids.

474 2. M A C R O S C O P I C

Wolfram E Q U A T I O N S FOR A S A M P L E M O D E L

2.1. S t r u c t u r e of t h e M o d e l

The model ~71 is based on a regular two-dimensional lattice of hexagonal cells, as illustrated in Fig. 1. The site at the centre of each cell is connected to its six neighbors by links corresponding to the unit vectors e~ through e 6 given by e~ = (cos(2~a/6), sin(2~a/6)) (2.1.1)

At each time step, zero or one particles lie on each directed link. Assuming unit time steps and unit particle masses, the velocity and momentum of each particle is given simply by its link vector e~. In this model, therefore, all particles have equal kinetic energy, and have zero potential energy. The configuration of particles evolves in a sequence of discrete time steps. At each step, every particle first moves by a displacement equal to its velocity e~. Then the particles on the six links at each site are rearranged according to a definite set of rules. The rules are chosen to conserve the number and total momentum of the particles. In a typical case, pairs of particles meeting head on might scatter through 60 ~ as would triples of particles 120 ~ apart. The rules may also rearrange other configurations, such as triples of particles meeting asymmetrically. Such features are important in determining parameters such as viscosity, but do not affect the form of the macroscopic equations derived in this section. To imitate standard physical processes, the collision rules are usually chosen to be microscopically reversible. There is therefore a unique predecessor, as well as a unique successor, for each microscopic particle configuration. The rules for collisions in each cell thus correspond to a simple permutation of the possible particle arrangements. Often the rules are self-inverse. But in any case, the evolution of a complete particle configuration can be reversed by applying inverse collision rules at each site. The discrete nature of the cellular automaton model makes such precise reversal in principle possible. But the rapid randomization of microscopic particle configurations implies that very complete knowledge of the current configuration is needed. With only partial information, the evolution may be effectively irreversible. ~8'19/

2.2. Basis for Kinetic T h e o r y

Cellular automaton rules specify the precise deterministic evolution of microscopic configurations. But if continuum behavior is seen, an approximate macroscopic description must also be possible. Such a

Cellular A u t o m a t o n Fluids

475

description will typically be a statistical one, specifying not, for example, the exact configuration of particles, but merely the probabilities with which different configurations appear. A common approach is to consider ensembles in which each possible microscopic configuration occurs with a particular probability (e.g. Ref. 18). The reversibility of the microscopic dynamics ensures that the total probability for all configurations in the ensemble must remain constant with time. The probabilities for individual configurations may, however, change, as described formally by the Liouville equation. An ensemble is in "equilibrium" if the probabilities for configurations in it do not change with time. This is the case for an ensemble in which all possible configurations occur with equal probability. For cellular automata with collision rules that conserve momentum and particle number, the subsets of this ensemble that contain only those configurations with particular total values of the conserved quantities also correspond to equilibrium ensembles. If the collision rules effectively conserved absolutely no other quantities, then momentum and particle number would uniquely specify an equilibrium ensemble. This would be the case if the system were ergodic, so that starting from any initial configuration, the system would eventually visit all other microscopic configurations with the same values of the conserved quantities. The time required would, however, inevitably be exponentially long, making this largely irrelevant for practical purposes. A more useful criterion is that starting from a wide range of initial ensembles, the system evolves rapidly to ensembles whose statistical properties are determined solely from the values of conserved quantities. In this case, one could assume for statistical purposes that the ensemble reached contains all configurations with these values of the conserved quantities, and that the configurations occur with equal probabilities. This assumption then allows for the immediate construction of kinetic equations that give the average rates for processes in the cellular automaton. The actual evolution of a cellular automaton does not involve an ensemble of configurations, but rather a single, specific configuration. Statistical results may nevertheless be applicable if the behavior of this single configuration is in some sense "typical" of the ensemble. This phenomenon is in fact the basis for statistical mechanics in many different systems. One assumes that appropriate space or time averages of an individual configuration agree with averages obtained from an ensemble of different configurations. This assumption has never been firmly established in most practical cases; cellular automata may in fact be some of the clearest systems in which to investigate it. The assumption relies on the rapid randomization of microscopic con-

476

Wolfram

figurations, and is closely related to the second law of thermodynamics. At least when statistical or coarse-grained measurements are made, configurations must seem progressively more random, and must, for example, show increasing entropies. Initially ordered configurations must evolve to apparent disorder. The reversibility of the microscopic dynamics nevertheless implies that ordered initial configurations can always in principle be reconstructed from a complete knowledge of these apparently disordered states. But just as in pseudorandom sequence generators or cryptographic systems, the evolution may correspond to a sufficiently complex transformation that any regularities in the initial conditions cannot readily be discerned. One suspects in fact that no feasibly simple computation can discover such regularities from typical coarse-grained measurements. (~9'2~ As a result, the configurations of the system seem random, at least with respect to standard statistical procedures. While most configurations may show progressive randomization, some special configurations may evolve quite differently. Configurations obtained by computing time-reversed evolution from ordered states will, for example, evolve back to ordered states. Nevertheless, one suspects that the systematic construction of such "antithermodynamic" states must again require detailed computations of a complexity beyond that corresponding to standard macroscopic experimental arrangements. Randomization requires that no additional conserved quantities are present. For some simple choices of collision rules, spurious conservation laws can nevertheless be present, as discussed in Section 4.5. For most of the collision rules considered in this paper, however, rapid microscopic randomization does seem to occur. As a result, one may use a statistical ensemble description. Equilibrium ensembles in which no statistical correlations are present should provide adequate approximations for many macroscopic properties. At a microscopic level, however, the deterministic dynamics does lead to correlations in the detailed configurations of particles. 4 Such correlations are crucial in determining local properties of the system. Different levels of approximation to macroscopic behavior are obtained by ignoring correlations of different orders. Transport and hydrodynamic phenomena involve systems whose properties are not uniform in space and time. The uniform equilibrium ensembles discussed above cannot provide exact descriptions of such

4 The kinetic theory approach used in this paper concentrates on average particle distribution functions. An alternative but essentially equivalent approach concentrates on microscopic correlation functions (e.g., Refs. 21, 22).

Cellular Automaton Fluids

477

systems. Nevertheless, so long as macroscopic properties vary slowly enough, collisions should maintain approximate local equilibrium, and should make approximations based on such ensembles accurate.

2,3. K i n e t i c E q u a t i o n s

An ensemble of microscopic particle configurations can be described by a phase space distribution function which gives the probability for each complete configuration. In studying macroscopic phenomena, it is, however, convenient to consider reduced distribution functions, in which an average has been taken over most degrees of freedom in the system. Thus, for example, the one-particle distribution function fa(x, t) gives the probability of finding a particle with velocity e, at position x and time t, averaged over all other features of the configuration (e.g., Ref. 23). Two processes lead to changes in J~, with time: motion of particles from one cell to another, and interactions between particles in a given cell. A master equation can be constructed to describe these processes. In the absence of collisions, the cellular automaton rules imply that all particles in a cell at position X with velocity e, move at the next time step to the adjacent cell at position X + %. As a result, the distribution function evolves according to

L ( x + e,,, T + 1) = L ( X , T)

(2.3.1)

For large lattices and long time intervals, position and time may be approximated by continuous variables. One may define, for example, scaled variables x =6~.X and t = 6, T, where 6x, 6, < 1. In terms of these scaled variables, the difference equation (2.3.1) becomes f~(x + e.6x, t + 6,) - f . ( x , t) = 0 (2.3.2)

In deriving macroscopic transport equations, this must be converted to a differential equation. Carrying out a Taylor expansion, one obtains (24) 6 , a , f . + 6 xe. " Vf, + ~ 62 O, f o + 6~6 ,(e~ " V )OJ , +-~f ~(e. " V12f~ + O( 031=O (2.3.3) If all variations in the .fa are assumed small, and certainly less than O(1/6x, 1/6t), it suffices to keep only first-order terms in 6~, 6,. In this way one obtains the basic transport equation 8,f,(x, t) + %. Vf.(x, t) = 0 (2.3.4)

1 1 2

478

Wolfram

This has the form of a collisionless Boltzmann transport equation for fa (e.g., Ref. 25). It implies, as expected, that fa is unaffected by particle motion in a spatially uniform system. Collisions can, however, change fa even in a uniform system, and their effect can be complicated. Consider, for example, collisions that cause particles in directions el and e4 to scatter in directions e2 and %. The rate for such collisions is determined by the probability that particles in directions el and e 4 occur together in a particular cell. This probability is defined as the joint two-particle distribution function ~14~'~2~. The collisions deplete the population of particles in direction el at a rate ~2) 9 Microscopic rever~t14 sibility guarantees the existence of an inverse process, which increases the population of particles in direction el at a rate given in this case by P~). Notice that in a model where there can be at most one particle on each link, the scattering to directions e2 and e5 in a particular cell can occur only if no particles are already present on these links. The distribution function _F is constructed to include this effect, which is mathematically analogous to the Pauli exclusion principle for fermions. The details of collisions are, however, irrelevant to the derivation of macroscopic equations given in this section. As a result, the complete change due to collisions in a one-particle distribution function fu will for now be summarized by a simple "collision term" f2,, which in general depends on two-particle and higher order distribution functions. (In the models considered here, s is always entirely local, and cannot depend directly on, for example, derivatives of distribution functions.) In terms of (2,, the kinetic equation (2.3.3) extended to include collisions becomes OJ~ + % 9Vf. = f2. (2.3.5)

With the appropriate form for (2 a, this is an exact statistical equation for f~ (at least to first order in ~). But the equation is not in general sufficient to determine f , . It gives the time evolution of fa in terms of the two-particle and higher order distribution functions that appear in s a. The two-particle distribution function then in turn satisfies an equation involving three-particle and higher order distribution functions, and so on. The result is the exact BBGKY hierarchy of equations, ~23) of which Eq. (2.3.5) is the first level. The Boltzmann transport equation approximates (2.3.5) by assuming that t2~ depends only on one-particle distribution functions. In particular, one may make a "molecular chaos" assumption that all sets of particles are statistically uncorrelated before each collision, so that multiple-particle distribution functions can be written as products of one-particle ones. The distribution function ~2) is thus approximated as --14

Cellular Automaton Fluids

479

flf4(1-f2)(1-f3)(1-f5)(1-f6)-The resulting Boltzmann equations will

be used in Section 4. In this section, only the general form (2.3.5) is needed. The derivation of Eq. (2.3.5) has been discussed here in the context of a cellular automaton model in which particles are constrained to lie on the links of a fixed array. In this case, the maintenance of terms in (2.3.3) only to first order in 67, at is an approximation, and corrections can arise, as discussed in Section 2.5. (24) Equation (2.3.5) is, however, exact for a slightly different class of models, in which particles have a discrete set of possible velocities, but follow continuous trajectories with arbitrary spatial positions. Such "discrete velocity gases" have often been considered, particularly in studies of highly rarefied fluids, in which the mean distance between collisions is comparable to the overall system size. (11A4)

2.4. C o n s e r v a t i o n Laws

The one-particle distribution functions typically determine macroscopic average quantities. In particular, the total particle density n is given by ~f~=n

a

(2.4.1)

while the momentum density nu, where u is the average fluid velocity, is given by ~, % f , = nu

a

(2.4.2)

The conservation of these quantities places important constraints on the behavior of the f , . In a uniform system Vf~ = 0, so that Eq. (2.3.5) becomes c~,f, = (2a and Eqs. (2.4.1) and (2.4.2) imply y' f2a= 0

a

(2.4.3)

(2.4.4) (2.4.5)

e~Q~ = 0

a

Using the kinetic equation (2.3.5), Eq. (2.4.4) implies c3, ~ f , + ~ %" Vf~ = 0

a c/

(2.4.6)

480

Wolfram

With the second term in the form V. (Y]eaj~) , Eq. (2.4.6) can be written exactly in terms of macroscopic quantities as c3tn + V" (nu) = 0 (2.4.7)

This is the usual continuity equation, which expresses the conservation of fluid. It is a first example of a macroscopic equation for the average behavior of a cellular automaton fluid. Momentum conservation yields the slightly more complicated equation

eoL+Z eo(eo"V/o)=0

(/ u

(2.4.8)

Defining the momentum flux density tensor H U= ~ (ea)i (ea)jf~

a

(2.4.9)

Eq. (2.4.8) becomes

O,(nui) + OjH o.= 0

(2.4.10)

No simple macroscopic result for H U can, however, be obtained directly from the definitions (2.4.1) and (2.4.2). Equations (2.4.7) and (2.4.10) have been derived here from the basic transport equation (2.3.5). However, as discussed in Section 2.3, this transport equation is only an approximation, valid to first order in the lattice scale parameters ,$x, ,st.(24) Higher order versions of (2.4.7) and (2.4.10) may be derived from the original Taylor expansion (2.3.3), and in some cases, correction terms are obtained. ~24) Assuming ,$x = 6, = 6, Eq. (2.4.6) to second order becomes ~ [(c~, + e~" V) + ~ 6(0, + ea" V)2] = 0 Writing the 0(,$) term in the form 0 , ~ (~, + % - V ) f ~ + V .~. (0, + e~. V ) e a f a

a a

(2.4.11)

(2.4.12)

this term is seen to vanish for any f~ that satisfy the first-order equations (2.4.7) and (2.4.10). Lattice discretization effects thus do not affect the continuity equation (2.4.7), at least to second order.

Cellular Automaton Fluids

481

Corrections do, however, appear at this order in the momentum equation (2.4.10). To second order, Eq. (2.4.8) can be written as

~ (~3~+ e~" V) e.f~, + ~,5 0, ~ (~3,+ e.. V) e~,f.

a a

+ ~6~

a

1

[ % - V 0 , + ( e , . V ) 2] e , f , = 0

(2.4.13)

The second term vanishes i f f , satisfies the first-order equation (2.4.8). The third term, however, contains a piece trilinear in the %, which gives a correction to the momentum equation (2.4.10)./241

2.5. C h a p r n a n - E n s k o g E x p a n s i o n

If there is local equilibrium, as discussed in Section 2.2, then the microscopic distribution functions f,(x, t) should depend, on average, only on the macroscopic parameters u(x, t) and n(x, t) and their derivatives. In general, this dependence may be very complicated. But in hydrodynamic processes, u and n vary only slowly with position and time. In addition, in the subsonic limit, lu] "~ 1. With these assumptions, one may approximate the f , by a series or Chapman-Enskog expansion in the macroscopic variables. To the order required for standard hydrodynamic phenomena, the possible terms are

f~=f l + c ( l ) e . . u + c (2) ( % . u ) 2 - ~ [ u l 2

+ dv I(e,- V)(%"u)-~1 V .u]+ ...} 21

{

L

']

(2.5.1)

where the c (i) are undetermined coefficients. The first three terms here represent the change in microscopic particle densities as a consequence of changes in macroscopic fluid velocity; the fourth term accounts for firstorder dependence of the particle densities on macroscopic spatial variations in the fluid velocity. The structures of these terms can be deduced merely from the need to form scalar quantities f, from the vectors %, u, and V. The relation (%), (ea)j = M f i , :

a

(2.5.2)

where here M = 6 and d = 2, and i and j are space indices, has been used in Eq. (2.5.1) to choose the forms of the lul 2 and Vu terms so as to satisfy the

482

Wolfram

constraints (2.4.1) and (2.4.2), independent of the values of the coefficients c (2) and c(2). In terms of (2.5.1), Eq. (2.4.1) yields immediately

f = n/6

(2.5.3)

while (2.4.2) gives C(1)=2 (2.5.4)

The specific values of C (2) and c~ ) can be determined only by explicit solution of the kinetic equation (2.3.5) including collision terms. (Some approximate results for these coefficients based on the Boltzmann transport equation will be given in Section 4.) Nevertheless, the structure of macroscopic equations can be derived from (2.5.1) without knowledge of the exact values of these parameters. For a uniform equilibrium system with u = 0, all the f~ are given by

f~ = f = n/6

(2.5.5)

In this case, the momentum flux tensor (2.4.9) is equal to the pressure tensor, given, as in the standard kinetic theory of gases, by

1 P o = ~ (%), ( % ) i f = ~ n 3u

a

(2.5.6)

where the second equality follows from Eq. (2.5.2). Note that this form is spatially isotropic, despite the underlying anisotropy of the cellular automaton lattice. This result can be deduced from general symmetry considerations, as discussed in Section 3. Equation (2.5.6) gives the equation of state relating the scalar pressure to the number density of the cellular automaton fluid:

p = n/2

(2.5.7)

When u ~ 0, H~j can be evaluated in the approximation (2.5.1) using the relations

(ea)~ (%)S (eo)k =0

a

(2.5.8)

and

M (%)i (%)./(%)I, ( % ) t - - - 2) (6~J6kz+ 6ik hit + 6izfsk) d(d+

a

(2.5.9)

Cellular A u t o m a t o n Fluids

483

The result is

HiJ=2n 6u+4ncl2) ldiRJ--2

I

'

lUl2 6 / J - ~ 4

In c(~) [ aildJ 2'v~176

Substituting the result into Eq. (2.4.10), one obtains the final macroscopic equation

~,(nu)+~nc <~> (u,V)u+

4

{

Iu(V,u)-~V lul ~It 1

(2.5.11)

1

1 : - ~1 V n - -8nc,_)V2u~ - ~lz

where = : u(u-

V)(nc ~21)- ~ I

lul 2 V(nc/2,) + (u" V)(ne~ )) - ~ (V.u) V(nc~ 1) (2.5.12)

The form (2.5.10) for Hv follows exactly from the Chapman-Enskog expansion (2.5.1). But to obtain Eq. (2.5.11), one must use the momentum equation (2.4.10). Equation (2.4.13) gives corrections to this equation that arise at second order in the lattice size parameter 6. These corrections must be compared with other effects included in Eq. (2.5.11). The rescaling x = g x X implies that spatial gradient terms in the Chapman-Enskog expansion can be of the same order as the O(6x) correction terms in Eq. (2.4.13). When the e~.u term in the Chapman-Enskog expansion (2.5.1) for thef~ is substituted into the last term of Eq. (2.4.13), it gives a contribution ~2a)

1 q?= -"(~nc ~1)V211 ~__ - ~ I nV2u

(2.5.13)

to the right-hand side of Eq. (2.5.11). Note that ~F depends solely on the choice of e~, and must, for example, vary purely linearly with the particle density f

2.6. Navier-Stokes Equation

The standard Navier-Stokes equation for a continuum fluid in d dimensions can be written in the form

o,(nu) + an(u- v) u = -Vp +,TV~u + (~ +-~,t)v(v- u)

(2.6.1)

484

Wolfram

where p is pressure, and t/and ~ are, respectively, shear and bulk viscosities (e.g., Ref. 27). The coefficient /~ of the convective term is usually constrained to have value 1 by Galilean invariance. Note that the coefficient of the last term in Eq. (2.6.1) is determined by the requirement that the term in H ~ / p r o p o r t i o n a l to t/ be traceless. (27'57) The macroscopic equation (2.5.11) for the cellular automaton fluid is close to the Navier-Stokes form (2.6.1). The convective and viscous terms are present, and have the usual structure. The pressure term appears according to the equation of state (2.5.7). There are, however, a few additional terms. Terms proportional to uVn must be discounted, since they depend on features of the microscopic distribution functions beyond those included in the Chapman-Enskog expansion (2.5.1). The continuity equation (2.4.7) shows that terms proportional to u ( V ' u ) must also be neglected. The term proportional to V lul 2 remains, but can be combined with the Vn term to yield an effective pressure term which includes fluid kinetic energy contributions. The form of the viscous terms in (2.5.11) implies that for the cellular automaton fluid, considered here, the bulk viscosity is given by =0 (2.6.2)

The value of t/ is determined by the coefficient c~ ~ that appears in the microscopic distribution function (2.5.1), according to

1

=- n v = - -~ n c ~ '

(2.6.3)

where v is the kinematic viscosity. An approximate method of evaluating c~ ) is discussed in Section 4.6. The convective term in Eq. (2.5.11) has the same structure as in the Navier-Stokes equation (2.6.1), but includes a coefficient 1 f l = ~ C (2) (2.6.4)

which is not in general equal to 1. In continuum fluids, the covariant derivative usually has the form D , = c ~ t + u . V implied by Galilean invariance. The cellular automaton fluid acts, however, as a mixture of components, each with velocities %, and these components can contribute with different weights to the covariant derivatives of different quantities, leading to convective terms with different coefficients.

Cellular A u t o m a t o n Fluids

485

The usual coefficient of the convective term can be recovered in Eq. (2.6.1) and thus Eq. (2.5.11) by a simple rescaling in velocity: setting = #u (2.6.5)

the equation for fi has coefficient 1 for the (fi'V)fi term. Small perturbations from a uniform state may be represented by a linearized approximation to Eqs. (2.4.7) and (2.5.11), which has the standard sound wave equation form, with a sound speed obtained from the equation of state (2.5.7) as c = l/x/2 (2.6.6)

The form of the Navier-Stokes equation (2.6.1) is usually obtained by simple physical arguments. Detailed derivations suggest, however, that more elaborate equations may be necessary, particularly in two dimensions (e.g., Ref. 28). The Boltzmann approximation used in Section 4 yields definite values for c ~2) and c~ I. Correlation function methods indicate, however, that additional effects yield logarithmically divergent contributions to c~ ~ in two dimensions (e.g., Ref. 29). The full viscous term in this case may in fact be of the rough form V z log(V 2) u.

2.7. Higher Order Corrections

The derivation of the Navier-Stokes form (2.5.11 ) neglects all terms in the Chapman-Enskog expansion beyond those given explicitly in Eq. (2.5.1). This approximation is expected to be adequate only when Jul <{ c. Higher order corrections may be particularly significant for supersonic flows involving shocks (e.g., Ref. 30). Since the dynamics of shocks are largely determined just by conservation laws (e.g., Ref. 27), they are expected to be closely analogous in cellular automaton fluids and in standard continuum fluids. For Fuj/c > 2, however, shocks become so strong and thin that continuum descriptions of physical fluids can no longer be applied in detail (e.g., Ref. 14). The structure of shocks in such cases can apparently be found only through consideration of explicit particle dynamics. (11,14) In the transonic flow regime lul ~ c , however, continuum equations may be used, but corrections to the Navier-Stokes form may be significant. A class of such corrections can potentially be found by maintaining terms O(u 3) and higher in the Chapman-Enskog expansion (2.5.1).

822/45/3-4 9

486

Wolfram

In the homogeneous fluid approximation Vu = 0, one may take

f~ = f { 1 + c(~/ea" u + c(2~[(ea" u)2 + 0-2 lu[ 2] + c(3)[(e~" u)3 + 03

lul 2

(ea" U)]

lu14] + ...} (2.7.1) (2.7.2)

-t-C(4)[-(ea" U)4-[-0-4,1 lul2 (ea.a)2+o-4,2

The constraints (2.4.1) and (2.4.2) imply

c(l)=d a2 ~r33 d ( d + 2) 1

1

d 3 d+2

(2.7.3) (2.7.4)

~-d0-4,1+04,2 =0

(2.7.5)

where d is the space dimension, equal to two for the model of this section. Corrections to (2.5.11) can be found by substituting (2.7.1) in the kinetic equation (2.4.8). For the hexagonal lattice model, one obtains, for example,

O,(nux) + ~ nc~21(uxOxu~ + ux~?yuy + uyOyux - u/? xuv)

q- ~ r/C(4){ [(5 q- 40"4,1) .3 _ 3Ux.23 63xUx

+ [(3 + 20"4,1) u 3 + (3 + 604,1 ) U2Uy] ~ybl x [(3 + 40-4,1) U3 -}- 3UxUy ] OxUy 2 v

+[(l+2aa,~)u3+(9+6a4,1)uxuZ]c~yuy}=O

(2.7.6)

The O(u 2) term in Eq. (2.7.6) has the isotropic form given in Eq. (2.5.11). The O(u 4) term is, however, anisotropic. To obtain an isotropic O(u 4) term, one must generalize the model, as discussed in Section 3. One possibility is to allow vectors ea corresponding to corners of an M-sided polygon with M > 6. In this case, the continuum equation deduced from the Chapman-Enskog expansion (2.7.1) becomes

~?~(nu)+ ~1 nc(2 ) (u'V)u q- U(V" U)---~ V

[

'

lul 2

]

lul2]+u(u-V)lu12}=0

(2.7.7)

1 +-~nc~4" + 0 - 4 , a ) { l u l 2 l - ( u . V ) u + u ( V . u ) - V '(1

Cellular Automaton

Fluids

487

This gives a definite form for the next-order corrections to the convective part of the Navier-Stokes equation. Corrections to the viscous part can be found by including terms proportional to Vu in the Chapman-Enskog expansion (2.7.1). The possible fourth-order terms are given by contractions of uiujSkut with products of (e.)m or 5mn. They yield a piece in the Chapman-Enskog expansion of the form c~)[271(e, - u) 2 (%" V)(e, "u) + re lul 2 (%" V)(e,-u) +%(%'u)(u" V)(e,'u)+%(e,'u) 2 (V'u)+275 lul 2 (V- u)] (2.7.8)

where Eq. (2.4.1) implies the constraints (for d = 2) 271 -~- 2273 = 0 rl + 4r2 +

427 4 + 8"C 5 =

(2.7.9) (2.7.10)

0

The resulting continuum equations may be written in terms of vectors formed by contractions of uiujSkS/Um and u~Siukc?tu m. The complete result is

8r(nu)+~nc (2) ( u ' V ) u + u ( V ' u ) - ~ V

lul 2

q-~ nCt4)(1 q- O'4,1){[Ul 2 E(u. V ) u + u(V. u ) - V

1 2 = --~nc(~)V 2u

lu[ 2 ] + u ( u - V ) l u l 2 }

1 nc(~ ) [((.rl

-- 4272_1 12Z4) u(V " U )2 _ (271 -- 4Z2 -F 4"c4) -

X u{V[(u 9 V) U] -- (U" V)(V-u)} -~- 8274{[(U" V) II]" V} u

1

+ 2 (~1 + 4272)E(v lu12) 9v ] u

+ 2~1 u ~ v - (V lu12) _ u . (V2u)

1

- 4274(v 9 u) V lu12

+ 2271uEu. (V2u)] +4-c z lu12 V2u'~ JJ

(2.7.11)

488

Wolfram

where, on the right-hand side, the first group of terms are all O((Vu)2), while the second group are O(VVu). Further corrections involve higher derivative terms, such as UiOj~k~lUm. For a channel flow with ux = ax 2, uy = 0, the time-independent terms in Eq. (2.7.11) have an x component

1 ac(~ ~+ 5 a 3 X 4 C ( ~ ) ( T 8

1 -'~

2r2 + 2~4) +_1 a2x3c(2) + a4x7c(4)(1 -}- ~ 2

(2.7.12)

and zero y component.

3. S Y M M E T R Y

CONSIDERATIONS

3.1. Tensor S t r u c t u r e The form of the macroscopic equations (2.4.7) and (2.5.11) depends on few specific properties of the hexagonal lattice cellular automaton model. The most important properties relate to the symmetries of the tensors EII t2 9 l- =Y', (ea)i~'"(%)~, !n) n

a

(3.1.1)

These tensors are determined in any cellular automaton fluid model simply from the choice of the basic particle directions ea. The momentum flux tensor (2.4.9) is given in terms of them by

Ho _- J ' ~elF(z) + ~ ~ i j k ~4k -~- p(2)F[::(4)~ l e(1)[=(3),, ~ L ~ i j k l ~ k ij "

+ c~)[ E~4{aku,+ 0-E~Z)akuk])

. . . .

-~-

o-E(2)Uk/,/k] (3.1.2)

where repeated indices are summed, and to satisfy the conditions (2.4.1) and (2.4.2) 0- = -E(4)ukk//E(2)ij (3.1.3)

The basic condition for standard hydrodynamic behavior is that the tensors E (nl for n 4 4 which appear in (3.1.2) should be isotropic. From the definition (3.1.1), the tensors must always be invariant under the discrete symmetry group of the underlying cellular automaton array. What is needed is that they should in addition be invariant under the full continuous rotation group. The definition (3.1.1) implies that the E(n) must be totally symmetric in their space indices. With no further conditions, the E(n) could have (n+,a-1) independent components in d space dimensions. Symmetries in the under-

Cellular A u t o m a t o n Fluids

489

lying cellular automaton array provide constraints which can reduce the number of independent components. Tensors that are invariant under all rotations and reflections (or inversions) can have only one independent component. Such invariance is obtained with a continuous set of vectors ea uniformly distributed on the unit sphere. Invariance up to finite n can also be obtained with certain finite sets of vectors %. Isotropic tensors F ("/ obtained with sets of M vectors e~ in d space dimensions must take the form

E(2n+l)=0

(3.1.4) ,d (2n) (3.1.5)

E (2n)=

M d(d+ 2 ) " " ( d + 2 n - 2 )

where Aij(2~_- 6o

z j (~jk/= 4)

(3.1.6) (3.1.7)

b ij~ kl -I- (~ik ~jl -~- ~ il~jk

and in general A (2n) consists of a sum of all the ( 2 n - 1)!! possible products of Kronecker delta symbols of pairs of indices, given by the recursion relation

2n

> ,,,= E

9 "-

vii!/

12""t]-ltjw

I ""i2n

(3.1.8)

j=2

The form of the d (2n) can also be specified by giving their upper simplicial components (whose indices form a nonincreasing sequence). Thus, in two dimensions, A (4)= [3, 0, 1, 0, 3] (3.1.9)

where the 1 l l l , 2111, 2211, 2221, and 2222 components are given. In three dimensions, A(4)=[3,0,1,0,3,0,0,0,0,1,0,1,0,0,3] Similarly, 3<6~= [5,0, 1,0, 1,0, 53 (3.1.11)

(3.1.1o)

490

Wolfram

and A(6)=[15,0,3,0,3,0,15,0,0,0,0,0,0,3, 0,1,0,3,0,0,0,0,3,0,3,0,0,15] in two and three dimensions, respectively. For isotropic sets of vectors %, one finds from (3.1.5) (3.1.12)

1 ?

so that for d = 2

( % ' v ) 2 n = Q2 ~ Ivl 2 n -

(2n-- 1)!!

d(d+ 2 ) ' " ( d + 2 n - 2 )

Ivl 2n

(3.1.13)

1

Q2 = ~ , while for d = 3

3

Q4=x,~

5

Q6 =T6' Q8-

35 128

(3.1.14)

1

Q2~, --

2n + 1

(3.1.15)

Similarly, -1 ~

a

(e,- v)2~ e , . v = Q2,, iv12,, v

(3.1.16)

In the model of Section 2, all the particle velocities ea are fundamentaly equivalent, and so are added with equal weight in the tensor (3.1.1). In some cellular automaton fluid models, however, one may, for example, allow particle velocities % with unequal magnitudes (e.g., Ref. 31). The relevant tensors in such cases are

E! ")

1112" " "

in

= ~ w(leal2)(%)il 999(%)i~

a

(3.1.17)

where the weights w(l%l 2) are typically determined from coefficients in the Chapman-Enskog expansion.

3.2. Polygons

As a first example, consider a set of unit vectors % corresponding to the vertices of a regular M-sided polygon: 2rca e a --=( COS--~--, sin 2__~) (3.2.1)

Cellular A u t o m a t o n Fluids Table I, Conditions for the Tensors E("~ of Eq. (3.1.1) to Be Isotropic w i t h the Lattice Vectors e~ Chosen to Correspond to the Vertices of Regular M - S i d e d Polygons

491

E{2) E(31 E(4) E151

E (6)

E(7)

M>2 M>~2, M4=3 M>2, M # 4 M~>2, M # 3 , 5 M>4, M r M~>2, M r 5, 7

F o r sufficiently large M, a n y t e n s o r E (~ c o n s t r u c t e d from these % m u s t be isotropic. T a b l e I gives the c o n d i t i o n s on M necessary to o b t a i n i s o t r o p i c E In). In general, it can be s h o w n t h a t E (~) is i s o t r o p i c if a n d only if M does n o t divide a n y of integers n, n - 2, n - 4 ..... (321 Thus, for example, E (=) m u s t be i s o t r o p i c w h e n e v e r n > M. In the case M = 6, c o r r e s p o n d i n g to the h e x a g o n a l lattice c o n s i d e r e d in Section 2, the E (=/are i s o t r o p i c u p to n = 5. T h e m a c r o s c o p i c e q u a t i o n s o b t a i n e d in this case thus have the usual h y d r o d y n a m i c form. H o w e v e r , a square lattice, with M = 4, yields an a n i s o t r o p i c E (41, given by

E(4)I M 4 = 25(4) (3.2.2)

where 6 o') is the K r o n e c k e r delta s y m b o l with n indices. T h e m a c r o s c o p i c e q u a t i o n o b t a i n e d in this case is

O,(nux) + 1_nc(2)(u~O~ux_ uvc3~uv)

. ~ . .

_

~?x n

1

~

1

1 8 (axu~ - a~uy)

Ox(nc~v ~)

which does n o t have the s t a n d a r d N a v i e r Stokes form. (6/'5

Note that even the linearized equation for sound waves is anisotropic on a square lattice. The waves propagate isotropically, but are damped with an effective viscosity that varies with direction, and can be negative. (33J

492

Wolfram

On a hexagonal lattice, E (4) is isotropic, but E (6) has the component form 1 E(6)IM=6 = - ~ [33, 0, 3, 0, 9, 0, 27] (3.2.4)

which differs from the isotropic result (3.1.11). The corrections (2.7.6) to the Navier-Stokes equation are therefore anisotropic in this case.

3.3. Polyhedra

As three-dimensional examples, one can consider vectors % corresponding to the vertices of regular polyhedra. Only for the five Platonic solids are all the I%] 2 equal. Table II gives results for the isotropy of the E (n) in these cases. Only for the icosahedron and dodecahedron is E (4) found to be isotropic, so that the usual h y d r o d y n a m i c equations are obtained. As in two dimensions, the E (2") for the cube are all proportional to a single K r o n e c k e r delta symbol over all indices. In five and higher dimensions, the only regular polytopes are the simplex, and the hypercube and its dual. (~4) These give isotropic E {n) only for n < 3, and for n < 4 and n < 4, respectively. In four dimensions, there are three additional regular polytopes, (34~ specified by Schlafi symbols {3, 4, 3}, {3, 3, 5}, and {5, 3, 3}. (The elements of these lists give the n u m b e r of edges a r o u n d each vertex, face, and 3-cell, respectively.) The {3, 4, 3} polytope has 24 vertices with coordinates corresponding to permutations of ( _+ 1, + 1, 0, 0). It yields E {") that are isotropic up to n=4, The { 3 , 3 , 5 } polytope has 120 vertices corresponding to ( +_ 1, +_ 1, + 1, _+ 1), all permutations of ( _ 2, 0, 0, 0), and

Table II. Isotropy of the Tensors E('~ w i t h % Chosen as the M Vertices of Regular Polyhedra"

ea

M 4 8 6 20 12

E (2)

E {3)

E {4)

E (5)

E (6)

Tetrahedron Cube Octahedron Dodecahedron Icosahedron

(1, 1, 1), cyc: (1, -1, -1) (_+1, _+1, _+1) cyc: (_+ 1, 0, 0) (_+1, +1, + l), cyc: (0, +~b l, +~b) cyc: (0, _+~, +1)

Y Y Y Y Y

N Y Y Y Y

N N N Y Y

N Y Y Y Y

N N N N N

In the forms for ea (which are given without normalization), the notation "cyc:" indicates all cyclic permutations. (All possible combinations of signs are chosen in all cases.) q/ is the golden ratio (1 + ,~)/2 ~ 1.618.

Cellular Automaton

Fluids

493

even-signature permutations of ( _+q~, _+ 1, ~b-~, 0), where ~b= (l + ~-5)/2. The {5, 3, 3} polytope is the dual of {3, 3, 5}. Both yield E(") that are isotropic up to n = 8.

3.4. Group Theory

The structure of the E (') was found above by explicit calculations based on particular choices for the ea. The general form of the results is, however, determined solely by the symmetries of the set of ea. A finite group G of transformations leaves the ea invariant. (For the hexagonal lattice model of Section 2, it is the hexagonal group $6.) In general G is a finite subgroup of the d-dimensional rotation group O(d). The ea form the basis for a representation of G, as do their products E(n). If the representation R (n) carried by the E(n) is irreducible, then the E('~ can have only one independent component, and must be rotationally invariant. But R {n) is in general reducible. The number of irreducible representations that it contains gives the number of independent components of E(') allowed by invariance under G. This number can be found using the method of characters (e.g., Refs. 35 and 36). Each class of elements of G in a particular representation R has a character that receives a fixed contribution from each irreducible component of R. Characters for the representation R {n) of G can be found by first evaluating them for arbitrary rotations, and then specializing to the particular sets of rotations (typically through angles of the form 2~z/k) that appear in G. To find characters for arbitrary rotations, one writes the E{') as sums of completely traceless tensors U (n) which form irreducible representations of O(d) (e.g., Ref. 37): E(n)= U{") + u(n-2)--~

"''

"~ U (0)

(3.4.1)

The characters of the E(") are then sums of the characters Z(m) for the irreducible tensors U (ml. For proper rotations through an angle ~b, the Z{m) are given by (e.g., Ref. 37)

z(m)(~) =

e 2~im~

( d = 2) ( d = 3) (3.4.2)

;~(m)(~b) = sin[(2m + 1) ~b/2] sin ~b/2

The resulting characters for the representations R (') formed by the E(~) are given in Table III. The number of irreducible representations in R {') can be found as usual by evaluating the characters for each class in R (n) (e.g., Ref. 35). Con-

494 Table III.

Wolfram Characters of Transformations of Totally Symmetric Rank n Tensors E(n) in d Dimensions"

Dimension 2

2

Rank 2

4

Character 4cz - 1

(4c 2 + 2c - 1)(4c2- 2c - 1)

4c 2 + 2c (2c + 1)(2c -- 1)(4c2+ 2c -- 1)

3 3

2 4

a c=cos(~b), where ~b is the rotation angle. For improper rotations in three dimensions 7r ~ must be used. b

sider as an example the case of R (4) with G the octahedral group O. This group has classes E, 8C3, 9C2, 6C4, where E represents the identity, and Ck represents a proper rotation by q~= 2 ~ / k about a k-fold symmetry axis. The characters for these classes in the representation R (4~ can be found from Table III. Adding the results, and dividing by the total number of classes in G, one finds that R (4) contains exactly two irreducible representations of O. Rank 4 symmetric tensors can thus have up to two independent components while still being invariant under the octahedral group. (38) In general, one may consider sets of vectors ea that are invariant under any point symmetry group. Typically, the larger the group is, the smaller the number of independent components in the E~nt can be. In two dimensions, there are an infinite number of point groups, corresponding to transformations of regular polygons. There are only a finite of nontrivial additional point groups in three dimensions. The largest is the group Y of symmetries of the icosahedron (or dodecahedron). Second largest is the cubic group E. As seen in Table II, only u gaurantees isotropy of all tensors E(~/ up to n = 4 (compare Ref. 39). It should be noted, however, that such group-theoretic considerations can only give upper bounds on the number of independent components in the E(n). The actual number of independent components depends on the particular choice of the eo, and potentially on the values of weights such as those in Eq. (3.1.16).

3.5. Regular Lattices

If the vectors ea correspond to particle velocities, then the possible displacements of particles at each time step must be of the form Za kaea 9 In discrete velocity gases, particle positions are not constrained. But in a

Cellular A u t o m a t o n Fluids

495

cellular automaton model, they are usually taken to correspond to the sites of a regular lattice. Only a finite number of such "crystallographic" lattices can be constructed in any space dimension (e.g., Refs. 40 and 41). As a result, the point symmetry groups that can occur are highly constrained. In two dimensions, the most symmetrical lattices are square and hexagonal ones. In three dimensions, the most symmetrical are hexagonal and cubic. The group-theoretic arguments of Section 3.4 suffice to show that in two dimensions, hexagonal lattices must give tensors E" that are isotropic up to n = 4, and so yield standard hydrodynamic equations (2.5.11). In three dimensions, group-theoretic arguments alone fail to establish the isotropy of E~a) for hexagonal and cubic lattices. A system with icosahedral point symmetry would be gauranteed to yield an isotropic E~4), but since it is not possible to tesselate three-dimensional space with regular icosahedra, no regular lattice with such a large point symmetry group can exist. Crystallographic lattices are classified not only by point symmetries, but also by the spatial arrangement of their sites. The lattices consist of "unit cells" containing a definite arrangement of sites, which can be repeated to form a regular tesselation. In two dimensions, five distinct such Bravais lattice structures exist; in three dimensions, there are 14 (e.g., Refs. 40 and 41). Sites in these lattices can correspond directly to the sites in a cellular automaton. The links which carry particles in cellular automaton fluid models are obtained by joining pairs of sites, usually in a regular arrangement. The link vectors give the velocities e, of the particles. In the simplest cases, the links join each site to its nearest neighbors. The regularity of the lattice implies that in such cases, all the e a are of equal length, so that all particles have the same speed. For two-dimensional square and hexagonal lattices, the % with this nearest neighbor arrangement have the form (3.2.1). The results of Section 3.2 then show that with hexagonal lattices, such e a give E/"/ that are isotropic up to n = 4, and so yield the standard hydrodynamic continuum equations (2.6.1). Table IV gives the forms of E(') for the most symmetrical three-dimensional lattices with nearest neighbor choices for the %. None yield isotropic E (41 (compare Ref. 38). The hexagonal and face-centered cubic lattices, which have the largest point symmetry groups in two and three dimensions, respectively, are also the lattices that give the densest packings of circles and spheres (e.g., Ref. 42). One suspects that in more than three dimensions (compare Ref. 43) the lattices with the largest point symmetry continue to be those with the densest sphere packing. The spheres are placed on lattice sites; the

496 Table IV.

Wolfram Forms of the Tensors E*n> for the M o s t Symmetrical Three-Dimensional Bravais Lattices a

e. M

E(21

E (4)

E~6}

Primitive cubic cyc: ( + 1, 0, 0) Body-centeredcubic (_+1, ~1, ? Face-centered cubic cyc: ( _+1, + 1, 0)

6 8 12

26 {2) 26 (4) 26 (6) 86 (2) 8(z/14) 2c]14)) 8(z/16)--2zJ14,2)+16616)) 86 (2) 4(z/14)-- 6 (4]) 4(zJ(4'2)-- 136(6))

"The basic vectors e, (used here without normalization) are taken to join each site with its M nearest neighbors. 6 C~)represents the Kronecker delta symbol of n indices; A(m represents the rotationally tensor defined in Eqs. (3.1.6) (3.1.8). A {n''l is the sum of all possible products of pairs of Kronecker delta symbols with n and m indices, respectively. p o s i t i o n s of their n e a r e s t n e i g h b o r s are defined b y a V o r o n o i p o l y h e d r o n or W i g n e r - S e i t z cell. T h e d e n s e s t sphere p a c k i n g is o b t a i n e d w h e n this cell, a n d t h u s the n e a r e s t n e i g h b o r vectors % , are closest to f o r m i n g a sphere. I n d i m e n s i o n s d ~ 8 , it has b e e n f o u n d t h a t the o p t i m a l lattices for sphere p a c k i n g are t h o s e b a s e d o n the sets of r o o t vectors for a s e q u e n c e of s i m p l e Lie g r o u p s (e.g., Ref. 44). Results o n the i s o t r o p y of the t e n s o r s F {n) for these lattices are g i v e n in T a b l e V. M o r e i s o t r o p i c sets of % c a n be o b t a i n e d b y a l l o w i n g links to j o i n sites o n the lattice b e y o n d n e a r e s t n e i g h b o r s . (31) O n a s q u a r e lattice, o n e m a y , for e x a m p l e , i n c l u d e d i a g o n a l links, y i e l d i n g a set of v e c t o r s % = (0, -[-1), ( -t-1, 0), ( ~ 1 , ! l )

(3.5.1)

Table V. Sequence of Simple Lie Groups Whose Sets of Root Vectors Yield Optimal Lattices for Sphere Packing in d Dimensions ~

d 1 2 3 4 5 6

Group A~ A2 A3 D4 D5

E6

M

f/max

SU(2) SU(3) SU(4) SO(8) SO(10)

2 6 12 24 40 72

4 2 4 2 0

"These lattices may also yield maximal isotropy for the tensors E('< Results are given for the maximum even n at which the E~'1 are found to be isotropic. The root vectors are given in Ref. 45.

Cellular A u t o m a t o n Fluids

497

Including weights w(l%l 2) as in Eq. (3.1.16), this choice of e, yields E(2/= 2[w(1 ) + 2w(2)] c5(2) E(4) = 4w(2) A (4) + 2[w(1 ) - 4w(2)] c~(4) (3.5.2) (3.5.3)

If the ratio of particles on diagonal and orthogonal links can be maintained so that w(!) = 4w(2) (3.5.4)

then Eq. (3.5.3) shows that E(4) will be isotropic. This choice effectively weights the individual vectors (0, _+ 1) and (_+ 1, 0) with a factor xf2. As a result, the vectors (3.5.1) are effectively those for a regular octagon, given by Eq. (3.2.1) with M = 8 . Including all 24 e a with components [(e,)i[ ~<2 on a square lattice, one obtains E(2)= 2[w(1) + 2w(2) + 4w(4) + 10w(5) + 8w(8)] 3 (2) E (4/= 4[w(2) + 8w(5) + 16w(8)] A (4) + 2[w(1)--4w(2) + 16w(4)-- 1 4 w ( 5 ) - 64w(8)] 6/4) 4 E(6) =~ [w(2) + 20w(5) + 64w(8)] A [6) (3.5.7) + 2[w(1) - 8w(2) + 64w(4) - 70w(5) - 512w(8)] 6(6) With w(5)= w(8)= O, E(4) and E(6) are isotropic if w(2) 3 w(4) 1 (3.5.6) (3.5.5)

w(1) 8'

w(1) 32

(3.5.8)

They cannot both be isotropic if w(4) also vanishes. In three dimensions, one may consider a cubic lattice with sites at distances 1, x/2, and x ~ joined. The ea in this case contain all those for primitive, face-centered, and body-centered cubic lattices, as given in Table IV. The E(n) can then be deduced from the results of Table IV, and are given by E12) = 2[w(1 ) + 4w(2) + 4w(3)3 6 (2) El4)= 4[w(2) + 2w(3)] A/4) + 2[w(1) -- 2w(2) -- 8w(3)] 6 (a) (3.5.9) (3.5.10)

E (6) = 8w(2) A (6~+ 4Ew(2) - 4w(3)] A (4,2/+ 2[w(1) -- 26w(2) + 64w(3)] &i6) (3.5.11)

498

Wolfram

Isotropy of E(4) is obtained when w(1) = 2w(2) + 8w(3) and of E (61 when w(1) = lOw(2) = 40w(3) (3.5.13)

(3.5.12)

Notice that (3.5.12) and (3.5.13) cannot simultaneously be satisfied by any nonzero choice of weights. Nevertheless, so long as (3.5.12) holds, isotropic hydrodynamic behavior is obtained in this three-dimensional cellular automaton fluid. Isotropic E(6) can be obtained by including in addition vectors % of the form (_+2, 0, 0) (and permutations), and choosing w(2)= w(1), w(3) = ~ w(l), w(4) = ]-~ w(1) (3.5.14)

The weights in Eq. (3.1.17) give the probabilities for particles with different speeds to occur. These probabilities are determined by microscopic equilibrium conditions. They can potentially be controlled by using different collision rules on different time steps (as discussed in Section 4.9). Each set of collision rules can, for example, be arranged to yield each particle speed with a certain probability. Then the frequency with which different collision rules are used can determine the densities of particles with different speeds.

3.6. I r r e g u l a r Lattices

The general structure of cellular automaton fluid models considered here requires that particles can occur only at definite positions and with definite discrete velocities. But the possible particle positions need not necessarily correspond with the sites of a regular lattice. The directions of particle velocities should be taken from the directions of links. But the particle speeds may consistently be taken independent" of the lengths of links. As a result, one may consider constructing cellular automaton fluids on quasilattices (e.g., Ref. 46), such as that illustrated in Fig. 2. Particle velocities are taken to follow the directions of the links, but to have unit magnitude, independent of the spatial lengths of the links. Almost all intersections involve just two links, and so can support only two-particle interactions. These intersections occur at a seemingly irregular set of points, perhaps providing a more realistic model of collisions in continuum fluids.

Cellular A u t o m a t o n Fluids

499

M=3

M=5

M-7

Fig. 2. Lattices and quasilattices constructed from grids oriented in the directions of the vertices of regular M-sided polygons. An appropriate dual of the M = 5 pattern is the Penrose aperiodic tiling.

The possible e~ on regular lattices are highly constrained, as discussed in Section 3.5. But it is possible to construct quasilattices which yield any set of e~. Given a set of generator vectors go, one constructs a grid of equally spaced lines orthogonal to each of them. (47) The directions of these lines correspond to the e~. If the tangents of the angles between the g(, are rational, then these lines must eventually form a periodic pattern, corresponding to a regular lattice. But if, for example, the g, correspond to the vertices of a pentagon, then the pattern never becomes exactly periodic, and only a quasilattice is obtained. A suitable dual of the quasilattice gives in fact the standard Penrose aperiodic tiling. ~48) In three dimensions, one may form grids of planes orthogonal to generator vectors go. Possible particle positions and velocities are obtained from the lines in which these planes intersect. Continuum equations may be derived for cellular automaton fluids on quasilattices by the same methods as were used for regular lattices above. But by appropriate choices of generator vectors, three-dimensional quasilattices with effective icosahedral point symmetry may be obtained, so that isotropic fluid behavior can be obtained even with a single particle speed. Quasilattices yield an irregular array of particle positions, but allow only a limited number of possible particle velocities. An entirely random lattice would also allow arbitrary particle velocities. Momentum conservation cannot be obtained exactly with discrete collision rules on such a lattice, but may be arranged to hold on average.

500

Wolfram

4. E V A L U A T I O N OF T R A N S P O R T COEFFICIENTS 4.1. Introduction Section 2 gave a derivation of the general form of the hydrodynamic equations for a sample cellular automaton fluid model. This section considers the evaluation of the specific transport coefficients that appear in these equations. While these coefficients may readily be found by explicit simulation, as discussed in the second paper in this series, no exact mathematical procedure is known for calculating them. This section considers primarily an approximation method based on the Boltzmann transport equation. The results obtained are expected to be accurate for certain transport coefficients at low particle densities. 6 4.2. Basis for Boltzrnann Transport Equation The kinetic equation (2.3.5) gives an exact result for the evolution of the one-particle distribution function f , . But the collision term f2a in this equation depends on two-particle distribution functions, which in turn depend on higher order distribution functions, forming the BBGKY hierarchy of kinetic equations. To obtain explicit results for the f,+ one must close or truncate this hierarchy. The simplest assumption is that there are no statistical correlations between the particles participating in any collision. In this case, the multiparticle distribution functions that appear in f2o can be replaced by products of one-particle distribution functions f , , yielding an equation of the standard Boltzmann transport form, which can in principle be solved explicitly for the fa. Even if particles were uncorrelated before a collision, they must necessarily show correlations after the collision, As a result, the factorization of multiparticle distribution functions used to obtain the Boltzmann transport equation cannot formally remain consistent. At low densities, it may nevertheless in some cases provide an adequate approximation. Correlations produced by a particular collision are typically important only if the particles involved collide again before losing their correlations. At low densities, particles usually travel large distances between collisions, so that most collisions involve different sets of particles. The particles involved in one collision will typically suffer many other collisions before meeting again, so that they are unlikely to maintain correlations. At high

6 Some similar results have been obtained by a slightly different method in Ref. 49.

Cellular A u t o m a t o n Fluids

501

densities, however, the same particles often undergo many successive collisions, so that correlations can instead be amplified. In the Boltzmann transport equation approximation, correlations and deviations from equilibrium decay exponentially with time. Microscopic perturbations may, however, lead to collective, hydrodynamic, effects, which decay only as a power of time. (29) Such effects may lead to transport coefficients that are nonanalytic functions of density and other parameters, as mentioned in Section 2.6.

4.3. Construction of Boltzmann Transport Equation

This subsection describes the formulation of the Boltzmann transport equation for the sample cellular automaton fluid model discussed in Section 2. The possible classes of particle collisions in this model are illustrated in Fig. 3. The rules for different collisions within each class are related by lattice symmetries. But, as illustrated in Fig. 3, several choices of overall rules for each class are often allowed by conservation laws. In the simplest case, the same rule is chosen for a particular class of collisions at every site. But it is often convenient to allow different choices of rules at different sites. Thus, for example, there could be a checkerboard arrangement of sites on which two-body collisions lead alternately to scattering to the left and to the right. In general, one may apply a set of rules denoted by k at some fraction 7k of the sites in a cellular automaton. (A similar procedure was mentioned in Section 3.5 as a means for obtaining isotropic behavior on three-dimensional cubic lattices.) The randomness of microscopic particle configurations suggests that the 7k should serve merely to change the overall probabilities for different types of collisions. The term f2 a in the kinetic equation (2.3.5) for f , is a sum of terms representing possible collisions involving particles of type a. Each term gives the change in the number of type a particles due to a particular type of collisions, multiplied by the probability for the arrangement of particles involved in the collision to occur. In the Boltzmann equation approximation, the probability for a particular particle arrangement is taken to be a simple product of the densities fb for particles that should be present, multiplied by factors ( 1 - f c ) for particles that should be absent. The complete Boltzmann transport equation for the model of Section 2 thus becomes Otf~ + eu. Vf. = Q . (4.3.1)

822/45/3 4-I0

502

Wolfram

++ - 2

,.

-X

where + Y3s[A(1, 3, 5 ) - A(0, 2, 4)]

~-~

5

Fig. 3. Possible types of initial and final states for collisions in the cellular automaton fluid model of Section 2.

12 = [72cA(1, 4 ) + ( 7 2 - ~)2L)A(2, 5)] --72A(0, 3)

+73A[A(2, 4, 5) +A(1, 2, 5 ) - A ( 0 , 3, 5 ) - A ( 0 , 2, 3)

+ A(1, 4,5)+ A ( 1 , 2 , 4 ) - A ( O , 3, 4)-A(O, 1,3)]

"4- []~4A(1, 2, 4, 5)-- 74LA(0, 2, 3, 5 ) - (~4-- ~4L) A(0, 1, 3, 4)3 Here (4.3.2) (4.3.3)

L+il

Aa(il, i2 ..... i k ) = l - f ~ + i ~ l - f ~ + i ~

9" 1 - f ~ + i k j = l (1-f.+j) L+,~ ~

where all indices on the fb are evaluated modulo M, and in this case M = 6. Note that in Eq. (4.3.2), the index a has been dropped on both 12 and A. The Boltzmann transport equations for any cellular automaton fluid

Cellular Automaton Fluids

503

model have the overall form of Eqs. (4.3.1) and (4.3.2). In a more general case, the simple addition of constants 1.']to the indices a in the definition of A can be replaced by transformations with appropriate lattice symmetry group operations. Independent of the values of the ?k, f2~ is seen to satisfy the momentum and particle number constraints (2.4.4) and (2.4.5). In the following calculations it is often convenient to maintain arbitrary values for the ?~ so as to trace the contributions of different classes of collisions. But to obtain a form for f2, that is invariant under the complete lattice symmetry group, one must take

1

~22L = ~)2R = 2 )22

(4.3.4) (4.3.5)

1

~)4L ~- ]14R ---~"2 ~)4

4.4. Linear A p p r o x i m a t i o n to B o l t z m a n n T r a n s p o r t Equation

In studying macroscopic behavior, one assumes that the distribution functions f , differ only slightly from their equilibrium values, as in the Chapman-Enskog expansion (2.5.1). The f , may thus be approximated as f , = f ( 1 +~b,) (]~b,j ~ 1) (4.4.1)

With this approximation, the collision term f2, in the Boltzmann transport equation may be approximated by a power series expansion in the ~b : a ~ b vb

b b,c

o.b+.~bb~bc +

""

(4.4.2)

The matrix e) (1) here is analogous to the usual linearized collision operator (e.g., Ref. 26). Notice that for a cellular automaton fluid model with collisions involving at most K particles, the expansion (4.4.2) terminates at

o(~K).

Microscopic reversibility immediately implies that the tensors co(') are all completely symmetric in their indices. The conservation laws (2.4.4) and (2.4.5) yield conditions on all the ~o(') of the form Z ~'+(') - n

~ abc.. -abe.. v

(4.4.3) (4.4.4)

y' e ~,,(") - n

a~abc.. -u abe..

In the particular case of co~

the more stringent conditions c~ . v 0 b

(4.4.5)

504

Wolfram

and e ~ab'"(~)-- 0

b

(4.4.6)

also apply. In the model of Section 2, all particle types a are equivalent up to lattice symmetry transformations. As a result, r,,(,)+ 1 )bc.. is always given simply (a by a cyclic shift of .,,(') , so that the complete form of 09(n) can be deterabc.. mined from the first row r,~(') The o9~') are thus circulant tensors (e.g., Wlbc." Ref. 50), and the values of their components depend only on numerical differences between their indices, evaluated modulo M. Expansion of (4.3.2) now yields 0)(1) = f2(1 - f ) circ{ - ['/2f 2 -t- (73s + 473A)ff-I- 74f2],

72Lf 2+ (y3S+273A)ff+]~4Lf 2,

(1 --Tzc)fZ + (__73S+273A)ff+ (74--TaL)f2, -- E72f 2 + (--73s + 473A)if-I- 74f2],

7 : o f 2 + ( -- 73s + 273A ) f f +

3)4Lf2,

(4.4.7)

0,

( 1 - - ~ 2 L ) f 2 + ( 7 3 S + 2 7 3 A ) , f f + ( V a - - 7 4 L ) f 2}

where f = (1 --f). Taking for simplicity 72 = 1, 72L = 89 ~3s = 1, 73A = ~4i = one finds co~) =f2(l _ f ) 2 circ - 1, ~ (1 + f ) ' 2 ( 1 - 3 f ) , 2 f - 1,~ (l - f ) ' 2 (1 + f )

I

'

,

,

]

(4.4.8)

0 f(f--

1)

f ( 3 f - 1) 2f(f - 1) o

50(2 ) - - __l

.be -- 2 f2(1 a

__

f ) circ

1) 0 f( 3f -- 1 ) 2,/'(f - 1 ) - 2 ( f - 1 ) ( 2 f - 1) f ( 5 f - 3) f ( 3 f - 1) ( f - 1 ) ( 2 f - 1) - f ( f - 1) - 2f 2

--f(f-

1) 2 f ( 3 f - 2) ( f - 1) ( 2 f - 1)

-f(f-

- 2 ( f - 1 ) ( 2 f - 1) -f(5f-3)

-f(f0 -f(f - 1) - - f ( 5 f - 3)

f ( 3 f - 1) ( f - 1 ) ( 2 f - 1) 2 f ( 3 f - 2)

-f(f0 2f(f-

f(f-

1)

-2f 2 ( f - 1 ) ( 2 f - 1)

-f(5f-

1)

1) 1)

3) 2 f ( f - 1) 0

(4.4.9)

Cellular Automaton Fluids 4.5. A p p r o a c h t o E q u i l i b r i u m

505

In a spatially uniform system close to equilibrium, one may use a linear approximation to the Boltzmann equation (4.3.1): ~?~(fib~)= ~ c~ ,el,

b

(4.5.1)

This equation can be solved in terms of the eigenvalues and eigenvectors of (1) the matrix ~~ The circulant property of mob considerably simplifies the computations required. An M x M circulant matrix Uab can in general be written in the form(SO~

Uob=u[(a-b+l)modM]=U~lI+U1211+

... +UIM H(M-1)

(4.5.2)

where H = circ[0, 1, 0, 0,..., 0] (4.5.3)

is an M x M cyclic permutation matrix, and I is the M x M identity matrix. From this representation, it follows that all M x M circulants have the same set of right eigenvectors vc, with components given by v,.), = - - - ~ exp

1 21ci(c - 1 )(a - 1 )

,/M

M

(4.5.4)

Writing

M

[~(Z)-~- 2 UlaZa

the corresponding eigenvalues are found to be

1

(4.5.5)

)~ = F (exp I2~i(-~-- 1 ) ] )

(4.5.6)

Using these results, the eigenvectors of co(,~ for the model of Section 2 ) are found to be

506

Wolfram

vl =~-'~ (1, 1, 1, 1, 1, 1) v2 = 7 - ~ (1, or, - a * , - 1, - a , a*)

1 1

= (u

l

v3 = 7 - ~ (1, - ~ * , -o', 1, - a * , - a ) = (vs)* (4.5.7)

V4=-~--~ (1 , --1, 1, --1, 1, - 1 ) l(v2+v6)=

2 L (u165 2i where = exp(i~/3)= 89 + ix/3 ) and the corresponding eigenvalues are 21 = 0

1

2,/6 2,/6

1 ~(2 , 1, --1, - 2 , - 1 , 1)

` / L (0 , 1, 1, O, - 1 , - 1 )

22=0 23

~ 3f2 ~ 1

f) { [72(1 - f)2 + 473Af(1 - f ) + ?4f 2]

(4.5.8,

xf~I(l_f,2(Tff_?2c)+ f2(2_y4c)]}

4i 24 = -673sf~(1 - f ) ~ 25=(23)* 26=0

Combinations of the r corresponding to eigenvectors with zero eigenvalue are conserved with time according to Eq. (4.5.1). Three such combinations are associated with the conservation laws (2.4.1) and (2.4.2). v~ corresponds to ~ , r which is the total particle number density. (v2 + v6)/2 and (v2-v6)/2/ correspond, respectively, to the x and y components of the momentum density Za e~ba-

Cellular Automaton Fluids

507

The ~bo may always be written as sums of pieces proportional to each of the orthogonal eigenvectors vc of Eq. (4.5.7): ~b = ~ 0c(v~,)a a

c

(4.5.9)

The coefficients 01, (02 + 06)/2, and ( 0 2 - 0 6 ) / 2 i give the values of the conserved particle and momentum densities in this representation, and remain fixed with time. The general solution of Eq. (4.5.1) is given in terms of Eq. (4.5.9) by 0c(t) = 0~.(0) e ;'~' (4.5.10)

Equation (4.5.8) shows that for any positive choices of the 7x, all nonzero 2c have negative real parts. As a result, the associated 0, must decay exponentially with time. Only the combinations of ~b~ associated with conserved quantities survive at large times. This result supports the local equilibrium assumption used for the derivation of hydrodynamic equations in Section 2. It implies that regardless of the initial average densities ~b~, collisions bring the system to an equilibrium that depends only on the values of the macroscopic conserved quantities (2.4.1) and (2.4.2). One may thus expect to be able to describe the local state of the cellular automaton fluid on time scales large compared to J2c]-~ (2c~0) solely in terms of these macroscopic conserved quantities. [Section 4.2 nevertheless mentioned some effects not accounted for by the Boltzmann equation (4.3.l) that can slow the approach to equilibrium. ] One notable feature of the results (4.5.8) is that they imply that the final equilibrium values of the ~b~ are not affected by the choice of the parameters 72L and 74L, which determine the mixtures of two- and fourparticle collisions with different chiralities. When the rate for collisions with different chiralities are unequal, however, 23 and 25 acquire imaginary parts, which lead to damped oscillations in the ~b as a function of time. a When all the types of collisions illustrated in Fig. 3 can occur, Eq. (4.5.8) implies that momentum and particle number are indeed the only conserved quantities. If, however, only two-particle collisions are allowed, then there are additional conserved quantities. In fact, whenever symmetric three-particle collisions are absent, so that 73s=0, Eq. (4.5.8) implies that the quantities

~11/2 + i

Of =

Z

a=i

f~

(4.5.11)

where the index a is evaluated modulo M = 6, is conserved. Thus, independent of the value of 72L, the total momenta on the two sides of any line

508

Wolfram

(not along a lattice direction) through the cellular automaton must independently be conserved. If three-particle symmetric collisions are absent, the cellular automaton thus exhibits a spurious additional conservation law, which prevents the attainment of standard local equilibrium, and modifies the hydrodynamic behavior discussed in Section 2. Section 4.8 considers some general conditions which avoid such spurious conservation laws.

4.6. Equilibrium Conditions and Transport Coefficients

Section 4.5 discussed the solution of the Boltzmann transport equation for uniform cellular automaton fluids. This section considers nonuniform fluids, and gives some approximate results for transport coefficients. The Chapman-Enskog expansion (2.5.1) gives the general form for approximations to the microscopic distribution functions fa- The coefficients c (2) and c~ ) that appear in this expansion can be estimated using the Boltzmann transport equation (4.3.1) from the microscopic equilibrium condition 8,f~=0 (4.6.1) In estimating c (2), one must maintain terms in f2a to the second order in ~bb, but one can neglect spatial variation in the ~ba. As a result, the Boltzmann equation (4.3.1) becomes

(1) E c%b ~bb+ E C~ ~bb~bc= 0

b b,c

(4.6.2)

Substituting forms for the Oa from the Chapman-Enskog expansion (2.5.1), one obtains c(2) ~ co~l)(u" %)2 + [c(1)]2 ~ co~2~(u 9%)(u 9ec) = 0

b b,c

(4.6.3)

where c~1)= 2 according to Eq. (2.5.4). Using the forms for co(11 and 0 (2) determined by the expansion of Eq. (4.3.2), one finds that the two terms in (4.6.3) show exactly the same dependence on the 7k. The final result for c (2) is thus independent of the 7k, and is given by c (a) = 2(1 - 2f)/(1 - f ) (4.6.4)

In the Boltzmann equation approximation, this implies that the coefficient /~ of the n ( u ' V ) u term in the hydrodynamic equation (2.6.1) is ( 1 - 2 f ) / [ 2 ( 1 - f ) ] . Notice that, as discussed in Section 2.6, this coefficient is not in general equal to 1.

Cellular Automaton Fluids

509

The value of the coefficient C(v can be found by a slightly simpler 2) calculation, which depends only on the linear part coab of the expansion of (1) the collision term s Keeping now first-order spatial derivatives of the ~b=, one can determine C(v from the equilibrium condition 2) ~o(~)~bb= f % "V(bb

b

(4.6.5)

which yields

~, c~2)c~ab

b

b" V)(eb" u) = c(l)f(eu 9 V)(ea- u)

(4.6.6) Then

With the approximations used, Eq. (2.4.7) implies that V . u = 0 . Eq. (4.6.6) gives the result

c(~) = - 2 { 12f(1 - f ) [ 7 2 ( 1 _ f ) 2 + 473Af( 1 __f) + 74f2] } -1 (4.6.7) Using Eq. (2.6.3), this gives the kinematic viscosity of the cellular automaton fluid in the Boltzmann equation approximation as v = {1 2 f ( 1 - f ) [ 7 2 ( 1 - - f ) 2 + 4 7 3 A f ( l - - f ) Some particular values are

v=ElZf(1-f)3] 1

+ 74/2] } - '

(4.6.8)

(y2= l, v3a = 7 4 = 0 )

'

(4.6.9)

v=[12f(l_f)(l+2f_Zf2)]

(~)2=~3A = ~ 4 = 1)

For f = 1/6 one obtains in these cases v ~0.86 and v ~ 0.47, respectively, while for f = 1/3, v ~ 0.84 and n ~ 0.26.

4.7. A G e n e r a l N o n l i n e a r A p p r o x i m a t i o n

At least for homogeneous systems, Boltzmann's H theorem (e.g., Ref. 51 ) yields a general form for the equilibrium solution of the full nonlinear Boltzmann equation (4.3.1). The H function can be defined as H = Z j7 l o g ( L )

a

(4.7.1)

where

L-- f=

1 -f~

(4.7.2)

This definition is analogous to that used for Fermi-Dirac particles (e.g., Refs. 51, 52): the factors (1 - f = ) account for the exclusion of more than one particle on each link, as in Eq. (4.3.1). The microscopic reversibility of (4.3.1) implies that when the equilibrium condition 0 , H = 0 holds, all products jT~,ya2.., must be equal for all initial and final sets of particles

510

Wolfram

{al, a2 .... } that can participate in collisions. As a result, the log(jT,) must be simple linear combinations of the quantities conserved in the collisions. If only particle number and momentum are conserved, and there are no spurious conserved quantities such as (4.5.11), the )7a can always be written in the form (7'49'54) J~a= exp(-c~ - flu. %) (4.7.3)

The one-particle distribution functions thus have the usual Fermi-Dirac form f , = [1 + exp(c~ + flu. ea)] -1 (4.7.4) where c~ and fl are in general functions of the conserved quantities n and

lul ~.

For small lul 2, one may write ~ = ~ o + ~ 1 [u[2+ "", f l = / ~ o + f l l qu[2+ "'" (4.7.5)

These expansions can be substituted into Eq. (4.7.4) and the results compared with the Chapman-Enskog expansion (2.7.1). For u = 0, one finds immediately the "fugacity relation"

f exp( - Co) = - 1-f

(4.7.6) that for generating Euler

Then, from the polynomials) (l+4eX) ~-1+~ r

1

expansion

4

(related

to

~ x (1+~1 + 42)x 3 6(1 + 4) 4

2-~$$~

r

4(1 - 4) x ~

- ~ ) ( 1 - 10~ + 42)x4+ ... 24(1+ 4) 5 (4.7.7)

together with the constraints (2.7.3)-(2.7.5) one obtains (for d = 2) 2 1-f

1-2f cq = ( l _ f ) 2

fil =

0~2

1-2/+2/2

(4.7.8)

(1 _ f ) 3

(1 - 2f)(3 - 4 f + 4f 2) 16(1 _ f ) 4

Cellular Automaton Fluids

511

where it has been assumed that the e a form an isotropic set of unit vectors, satisfying Eq. (3.1.5). The complete Chapman Enskog expansion (2.7.1) then becomes d21-2fF, e a ) 2 _ d lul2 ] f~=f{l+du.e~-~ 2 ] - - 7 L tu" d 3 1 - 6 f + 6f 2 [3 ] + 6 -~----f~ k t u ' e = ) 3 - d + 2 l u l 2 ( u ' % ) 1 1 - 2f [32(1 - 12f+ 12fZ)(u 9%)4 + 48 (1 _ f ) 3 + 3 8 4 f ( 1 - f ) ]uJ 2 ( u ' e = ) 2 + 3 ( l l - 3 6 f + 3 6 f 2) lul 4] + " ' t (4.7.9)

where for the last term it has been assumed that d = 2. The result (4.6.4) for c ~2~follows immediately from this expansion. For cellular automaton fluid models with E (6) isotropic, the continuum equation (2.7.7) holds. The results for the coefficients that appear in this equation can be obtained from the approximation (4.7.9), and have the simple forms c~2) = d2(1 - 2f) (4.7.10) 2(1 - - f ) C(4)(1 -I- 04,1)

=

2(1 - 2 ) 3 3(1 _ f f )

(d=2)

(4.7.11)

These results allow an estimate of the importance of the next-order corrections to the Navier-Stokes equations included in Eq. (2.7.7). They suggest that the corrections may be important whenever lu/(1 -f)212 is not small compared to 1. The corrections can thus potentially be important both at high average velocities and high particle densities. The hexagonal lattice model of Section 2 yields a continuum equation of the form (2.7.6), with an anisotropic O(u2Vu) term. Equation (4.7.9) gives in this case ~?,(6fux) q3/(1 - 2f) 1- f (UxC3xU:r+ uxc~YuY+ uyOyUx - uy~?xUy)

+/(14(1-2f)_f)~ {[(55 - 84f + 84fZ) u2 + 3(13 + 4 f - 4 f Z ) u 2 ] Ux~?xux + 2[(1 + 12/'- lZf 2) u2 + 9(1 - 2/)2 u23 UxC~yUy

2 + 6[-(1 + 1 2 f - 12f 2) ux+ (1 - 2 f ) 2] uyOyu~

+ 3[(13+4f--4f~)u2+(13-28f+28f2)u

2] UyOxUy}=O

(4.7.12)

512

Wolfram

3f(1 - 2f) a~(6fu~)+ l - f (-u~a~u~+u~a~u~+%a~u~-u~u~)

f(1 - 2 f ) {[(35 - 3 6 f + 36f 2) 2 + 3(17 - 44f + 44f 2) u~] + 4(1 f)3 ux + 2[-(1 + 1 2 / - 123"2) u2+ 9(1 - 2f) 2 u 23 + 6 [ ( 1 + 1 2 f - 12f 2) U2x + ( l - - 2 f )

ux(3yux

UxOxUy 2] uyaxU~ upC3yUy}----0

(4.7.13)

2 2] + 3[(17-- 4 4 f + 44f 2) u~ + (17 -- 12f+ 12f 2) Uy

The O(uVu) term is as given in Eq. (2.5.11), The O(u3Vu) terms are anisotropic, and are not even invariant under exchange of x and y coordinates (n/2 rotation). For small densities f, Eqs. (4.7.12) and (4.7.13) become

~,(Ux) + ~ (ux0~Ux + Ux#yUy+ uyayux- uy~xuy)

+ ~-~ [(55u~ + 39u 2) u~Oxux+ 2(u 2 + 9u 2) u~Oyuy +6(u~

2

1 1

1

+ uy)2uyOyux+ 39(u~ + u2)UyO~Uy]= 0

(4.7.14)

a,(u~) J--1 ( - U~axU~+ U~axUy+ U~xUx - uranus)

+ ~-~ [(35u 2 + 51u 2) UxOyU~ 2(u 2 + 9u 2) u~O~uy +

+ 6(u2+u2y)uy~xux+51(u~+u2)uySyuy]=O

(4.7.15)

The results (4.7.10) and (4.7.11) follow from the Fermi-Dirac particle distribution (4.7.4). If instead an arbitrary number of particles were allowed at each site, the equilibrium particle distribution (4.7.4) would take on the Maxwell Boltzmann form fa = exp( - ~ - flu" ea) (4.7.16)

With this simpler form, more complete results for f~ as a function of n and u can be found. Results which are isotropic to all orders in u can be obtained only for an infinite set of possible particle directions, parametrized, say, by a continuous angle 0. In this case, the number and momentum densities (2.4.1) and (2.4.2) may be obtained as integrals

~-~f~f(O) dO= f

(4.7.17)

Cellular A u t o m a t o n Fluids

513

1 ~,27z

~Jo

e(O)f(O)dO=fu

(4.7.18)

With the distribution (4.7.16), these integrals become

1 :'2~

J0 e x p ( - c ~ - flu cos 0) dO = e ~Io(flu) = f

(4.7.19) (4.7.20)

~1 12~ e x p ( - c ~ - f u c o s O ) o

cosOu/udO=-e-~I1(fu)

u/u=fu

where u = Ju[, and the Iv(z) are modified Bessel functions (e.g., Ref. 53)

(z/2)v+2n

/v(z)=,,=~o~-~-~n)!, Io(0)= 1, Iv(0)=0 ( v > 0 ) (4.7.21) (4.7.22) (4.7.23)

fl

1

(1-x2)V-1/2e

ZXdx=~zz-V(2v-1)!!Iv(z)

Iv_ ~(z) - I~+ ~(z) = (2v/z) I,(z)

The rapid convergence of the series (4.7.21) means that Eqs. (4.7.19) and (4.7.20) provide highly accurate approximations even for a small number of discrete directions ea. [For example, with M = 6 , e = 0 , f = 1, and u = (1, 1), the error in Eq. (4.7.19) is less than 10 9.] For the simple distribution (4.7.16) the momentum flux density tensor (2.4.9) may be evaluated in direct analogy with Eqs. (4.7.19) and (4.7.20) as

1

uiu j

1

(4.7.24)

Using the recurrence relation (4.7.23), and substituting the results (4.7.19) and (4.7.20), this may be rewritten in the form

1

2

uiu ~

1

(4.7.25)

Combining Eqs. (4.7.19) and (4,7.20), one finds that the function f ( f u) is independent of f, and can be determined from the implicit equation

Ii(flu)

I'o(flu)

Io(fU) Io(fU)

= --u

(4.7.26)

514

Wolfram

Expanding in powers of u 2, as in Eq. (4.7.5), yields /3(u)= ~ /3,,u2n n=0 /30 = - 2 , /31 = - 1 , 5

/32 = - g , /33 =

(4.7.27) 19

24'

143

/34 = - 18---6

Equation (4.7.19) then gives ~(u,f)= Z ~n u2" n=O 3 e2 = ~ , 25 ~3 = 3-6' (4.7.28) 133 ~4 = 192

exp(-eo) =f

cq = 1,

In the limit u ~ l,/3 ~ - ~ . The above results immediately yield values for the transport coefficients c (n) in the Chapman-Enskog expansion: c (2) -- 2 (4.7.29) (4.7.30)

C(4)(1 -t- 0"4,1) = 5

2

independent of density. Equation (4.7.29) implies that the coefficient # of the convective term in the Navier-Stokes equation (2.6.1) is equal to 1/2. The deviation from the Galilean invariant result 1 is associated with the constraint of fixed speed particles. Figure 4 shows the exact result for /3(u) obtained from Eq. (4.7.26), compared with series expansions to various orders. Significant deviations from the O(u 2) "Navier-Stokes" approximation are seen for u > 0.4. For Fermi-Dirac distributions of the form (4.7.4), the integrals (4.7.19) and (4.7.20) can only be expressed as infinite sums of Bessel functions. 4.8. O t h e r M o d e l s The results obtained so far can be generalized directly to a large class of cellular automaton fluid models. In the main case considered in Section 3, particles have velocities corresponding to a set of M unit vectors %. If this set is invariant under inversion, then both e a and - e a always occur. As a result, two particles colliding head on with velocities e a and - e a can always scatter in any directions % and - % with b ~ a. One simple possibility is to choose the

Cellular Automaton Fluids exact

515

3 - ( 2 + ~)

O(u s)

2

1

O ( u 2)

0

0.0

~

0.2

I 0.4

I 0.6

I 0.8

u Fig. 4. Dependence of fl(u) from Eq. (4.7,16) on the magnitude u of the macroscopic velocity. The results are for Maxwell-Boltzmann particles with unit speeds and arbitrary directions in two dimensions. The function fl(u) appears both in the microscopic distribution function (4.7.16) and in the macroscopic momentum flux tensor (4.7.25). The result for fl(u) from an exact solution of the implicit equation (4.7.26) is given, together with results from the series expansion (4.7.27). The O(u2) result corresponds to the Navier-Stokes approximation. Deviation from the exact result is seen for u > 0.4.

rules at different sites so that each scattering direction occurs with equal probability. If only such two-particle collisions are possible (as in a lowdensity approximation), and only one particle is allowed on each link, then the Boltzmann transport equation becomes 1 O , L + % ' V f a = --fafa H ( 1 - - f ~ ) + M _ - - - -~ ~. fbf~ 1~ ( 1 - f ~ . ) ( 4 . 8 . 1 ) c#a,~ b~a,,~ c~b,~ where fa is the distribution function for particles with direction - e a . To second order in the expansion (4.4.2) this gives O,(f~ba) + e a 9 V(f~ba) = f2(1 _ f ) M - 3

[_

' Y~ (C~a+ ~a) + m-------2b~a.a (~+~6)]

+f~(~ -f)'-~

(1 - 2f) -~o~ + M-----2~,b~,~

(4.8.2)

51 6

Wolfram

where Z ' denotes summation over the triangular region in which the indices form a strictly increasing sequence. The form of the ~b for a homogeneous system can be obtained from a the general equilibrium conditions of Section 4.7. The coefficients c (n) in the Chapman-Enskog expansion are then given by Equations (4.7.10) and (4.7.11). The convective transport coefficient # in the Navier-Stokes equation (2.6.1) is thus given by da(1 - 2/) #= 8(1 - f )

(4.8.3)

The C~v cannot be obtained by the methods of Section 4.7. But from n)

Eq. (4.8.2) one may deduce immediately the linearized collision term

co(11 -_- f 2 ( 1 - f ) m ab 3

circ[x~] 2

7~a- M _ 2

(4.8.4)

Z1 =)~l +M/2 = - 1 ,

Then, in analogy with Eq. (4.6.8), the kinematic viscosity for the cellular automaton fluid is found to be M-2 v -- 2 d ( d + 2) MfZ(1 _U)M 3 (4.8.5) For an icosahedral set of %, with d = 3 and M = 12, this yields v-=- [36f2(1 _ f ) 9 ] - 1 (4.8.6)

Several generalizations may now be considered. First, one may allow not just one, but, say, up to tc particles on each link of the cellular automaton array. In the limit ~ ~ oe an arbitrary density of particles is thus allowed in each cell. The Boltzmann equation for this case is the same as (4.8.1), but with all (1 -f~,) factors omitted. The resulting transport coefficients are U=~

1

(4.8.7) (4.8.8)

M-2 v - 2 d ( d + 2) MU2

Another generalization is to allow collisions that involve more than two particles. The simplest such collisions are "composite" ones, formed by superposing collisions involving two or less particles. The presence of such collisions changes the values of transport coefficients, but cannot affect the basic properties of the model. The four-particle and asymmetric three-par-

Cellular Automaton Fluids

517

ticle collisions in the hexagonal lattice model of Section 2 are examples of composite collisions. They increase the total collision rate, and thus, for example, decrease the viscosity, but do not change the overall macroscopic behavior of the model. In general, collisions involving k particles can occur if the possible e, are such that

k k

2 ea~= ~ ebi

i-1 i=1

(4.8.9)

for some sets of incoming and outgoing particles ai and bi. Cases in which all the ai and bi are distinct may be considered "elementary" collisions. In the hexagonal lattice model of Section 2, only two-particle and symmetric three-particle collisions are elementary. No elementary three-particle collisions are possible on primitive and body-centered cubic three-dimensional lattices, or with % corresponding to the vertices of icosahedra or dodecahedra. For a face-centered cubic lattice, however, eight distinct triples of % sum to zero [-an example is (1,-1,0)+(0,1,-1)+(-1,0, l)], so that elementary three-particle collisions are possible. One feature of the hexagonal lattice model discussed in Section 2 is the existence of the conservation law (4.5.11) when elementary three-body symmetric collisions are absent. Such spurious conservation laws exist in any cellular automaton fluid model in which all particles have the same speed, and only two-particle collisions can occur. Elementary three-particle collisions provide one mechanism for avoiding these conservation laws and allowing the equilibrium of Section 4.7 to be attained.

4.9. M u l t i p l e Speed M o d e l s

A further generalization is to allow particles with velocities e a of different magnitudes. This generalization is significant not only in allowing two-particle collisions alone to avoid the spurious conservation laws of Section 4.5, but also in making it possible to obtain isotropic hydrodynamic behavior on cubic lattices, as discussed in Section 3.5. One may define a kinetic energy 1/2 [%[2 that differs for particles of different speeds. In studies of processes such as heat conduction, one must account for the conservation of total kinetic energy. In many cases, however, one considers systems in contact with a heat bath, so that energy need not be conserved in individual collisions. In a typical case, one may then take pairs of particles with speed si colliding head on to give pairs of particles with some other speed sj. In general, different collision rules may be used on different sites, typically

822/45/3-4-11

518

Wolfram

following some regular pattern, as discussed in Section 4.3. Thus, for example, collisions between speed si particles may yield speed sj particles at a fraction 7i~ j of the sites. The number mi of possible particles with speed si= (t%[2) 1/2 that can occur at each site is determined by the structure of the lattice. The collision rules at different sites may be arranged, as in Section 4.8, to yield particles of a particular speed s~ with equal probabilities in each of the me possible directions. In a homogenous system, the probability f, for a link with speed s~ to be populated should satisfy the master equation

0t)7i=2/'6~j2 j = ~ ,j~ i

(--Ti~jmiYZ-k-Tj~imjjTgj 2)

(4.9.1)

where it is assumed that 7i~j = Z Yj~ ~= 1

J ]

(4.9.2)

and f is the reduced particle density given by Eq. (4.7.2). With two speeds, F u becomes

ru=(-Yl_2m,

271~2ml

72~,m2

--72~1m2/

~

(4.9.3)

The solutions of Eq. (4.9.1) can be found in terms of the eigenvalues and corresponding eigenvectors of this matrix: 2=0: 2 = -(Y,~2ml q-'Y2~ lm2): (72~ lm2, 71~2rnl) ( - 1 , 1) (4.9.4) (4.9.5)

In the large-time limit, only the equilibrium eigenvector (4.9.4) should survive, giving a ratio of reduced particle densities

? 2 -- 71 ~ 2 m l ~? 72~1m2

(4.9.6)

For three particle speeds, one finds the equilibrium conditions )722_ ml 71 ,1~12

273 ~ 1

-t- 7 1 ~ 2 7 3 ~ 2 -}" 7 1 ~ 3 7 3 ~ 2

m 2 7 2 + 173 + 1 - f 7 2 ~ 173 ~ 2 -t- 72 + 373 ~ 1

(4.9.7)

72 Y?

ml 71~272~3~_71~372~1..1_71~372~3 m3 7 2 ~ 173 ~ 1 -[- 7 2 ~ 173 ~ 2 t - 72 ~ 373 ~ 1

Different choices for the 7k yield different equilibrium speed distributions. The probabilities f, give the weights w(s~) that appear in Eq. (3.1.16). Equation (4.9.6) shows that by choosing

72~ 1~

271 4 2

(4.9.8)

Cellular Automaton Fluids

519

one obtains a ratio of weights for the model of Eq. (3.5.1) that satisfy the condition (3.5.5) for the isotropy of E(4). [There is a small correction to equality in Eq. (4.9.8) associated with the difference between f and ~.] On a cubic lattice, one may similarly satisfy the condition (3.5.12) for the isotropy of E(4) simply by taking f3 = 0, and 72 - I ~ 271 - 2 (4.9.9)

In this way, one may obtain approximate isotropic hydrodynamic behavior on a three-dimensional cubic lattice.

4.10. T a g g e d P a r t i c l e D y n a m i c s

In the discussion above, all the particles in the cellular automaton fluid were assumed indistinguishable. This section considers the behavior of a small concentration of special "tagged" particles. The density g~ of tagged particles with direction e, satisfies an equation of the Fokker-Planck type (e.g., Ref. 26):

8,g~+(%.V)g~=O~

(4.10.l)

Assuming as in the Boltzmann equation approximation that there are no correlations between particles at different sites, the collision term of Eq. (4.10.1) may be written in the form

Oa = ~ Oabgb

b

(4.10.2)

where 0,b gives the probability that a particle that arrives at a particular site from direction eb leaves in direction e, with a ~ b. The probability is averaged over different arrangements of ordinary particles. Various deterministic rules may be chosen for collisions between ordinary and tagged particles. The simplest assumption is that on average the tagged particles take the place of any of the outgoing particles with equal probability. Conservation of the total number of tagged particles implies

g~ = gM

a

(4,10.3)

The total momentum of tagged particles is not conserved; the background of ordinary particles acts like a "heat bath" which can exchange momentum with the tagged particles through the noise term O~. Assuming a uniform background fluid, one may make an expansion for the ga of the form ga = (g + d(l)ea" Vg + -'.) (4.10.4)

520

Wolfram

The total number of tagged particles then satisfies the equation

c~tg+d(l) ~1 E(2)~i~j g

= 0

(4.10.5)

where the collision term disappears as a result of Eq. (4.10.3). With the e, chosen so that E!2) is isotropic, Eq. (4.10.5) becomes the standard equation q for self-diffusion, ~?t g = DV2g (4.10.6) with the diffusion coefficient D given by D = - 1 d(~) d (4.10.7)

The value of d (1) must be found by solving Eq. (4.10.1) for g, using the approximation (4.10.4). The equilibrium condition for Eq. (4.10.1) in this case becomes (e a 9V) g = d(1) ~, 0ab%" Vg

b

(4.10.8)

Thus - d (~) is given in this approximation by the mean free path 2 for particle scattering, so that the diffusion coefficient is given by the standard kinetic theory formula D =~2 For the hexagonal lattice model of Section 2, D = {2f2(1 _ f ) 2 [(1 _ f ) 2 + (73s q_ 4~73A)f ( 1 _ f ) + 3~4f2] }-~ (4.10.10)

1

(4.10.9)

5. S O M E E X T E N S I O N S

The simple physical basis for cellular automaton fluid models makes it comparatively straightforward for them to include many of the physical effects that occur in actual fluid experiments. Boundaries can be represented by special sites in the cellular automaton array. Collisions with boundaries conserve particle number, but not particle momentum. One possibility is to choose boundary collision rules that exactly reverse the velocities of all particles, so that particles in a layer close to the boundary have zero average momentum. This choice yields macroscopic "no slip" boundary conditions, appropriate for many solid surfaces (e.g., Ref. 27). For boundaries that consist of flat segments

Cellular Automaton Fluids

521

aligned along lattice directions, an alternative is to take particles to undergo "specular" reflection, yielding a zero average only for the transverse component of particle momentum, and giving "free slip" macroscopic boundary conditions. The roughness of surfaces may be modeled explicitly by including various combinations of these microscopic boundary conditions (corresponding, say, to different coefficients of accommodation). Arbitrarily complex solid boundaries may be modeled by appropriate arrangements of boundary cells. To model, for example, a porous medium one can, for example, use a random array of "boundary" cells with appropriate statistical properties. A net flux of fluid can be maintained by continually inserting particles on one edge with an appropriate average momentum and extracting particles on an opposite edge. The precise arrangement of the inserted particles should not affect the macroscopic properties of the system, since microscopic processes should rapidly establish a microscopically random state of local equilibrium. Large-scale inhomogeneities, perhaps representing "free stream turbulence" (e.g., Ref. 4), can be included explicitly. External pressure and density constraints, whether static or timedependent, can be modeled by randomly inserting or extracting particles so that local average particle densities correspond to the macroscopic distribution required. External forces can be modeled by randomly changing velocities of individual particles so as to impart momentum to the fluid at the required average rate. Moving boundaries can then be modeled by explicit motion of the special boundary cells, together with the inclusion of an appropriate average momentum change for particles striking the boundary. Gravitational and other force fields can also be represented in a "quantized approximation" by explicit local changes in particle velocities. Many other physical effects depend on the existence of surfaces that separate different phases of a fluid or distinct immiscible fluids. The existence of such surfaces requires collective ordering effects within the system. For some choices of parameters, no such ordering can typically occur. But as the parameters change, phase transitions may occur, allowing large correlated regions to form. Such phenomena will be studied elsewhere. (Surface tension effects have been observed in other two-dimensional cellular automata. (3})

6. D I S C U S S I O N

Partial differential equations have conventionally formed the basis for mathematical models of continuum systems such as fluids. But only in

522

Wolfram

rather simple circumstances can exact mathematical solutions to such equations be found. Most actual studies of fluid dynamics must thus be based on digital computer simulations, which use discrete approximations to the original partial differential equations (e.g., Ref. 55). Cellular automata provide an alternative approach to modeling fluids and other continuum systems. Their basic constituent cells are discrete, and ideally suited to simulation by digital computers. Yet collections of large numbers of these cells can show overall continuum behavior. This paper has given theoretical arguments that with appropriate rules for the individual cells, the overall behavior obtained should follow that described by partial differential equations for fluids. The cellular automata considered give simple idealized models for the motion and collision of microscopic particles in a fluid. As expected from the second law of thermodynamics, precise particle configurations are rapidly randomized, and may be considered to come to some form of equilibrium. In this equilibrium, it should be adequate to describe configurations merely in terms of probabilities that depend on a few macroscopic quantities, such as momentum and particle number, that are conserved in the microscopic particle interactions. Such averaged macroscopic quantities change only slowly relative to the rate of particle interactions. Partial differential equations for their behavior can be found from the transport equations for the average microscopic particle dynamics. So long as the underlying lattice is sufficiently isotropic, many cellular automata yield in the appropriate approximation the standard Navier-Stokes equations for continuum fluids. The essential features necessary for the derivation of these equations are the conservation of a few macroscopic quantities, and the randomization of all other quantities, by microscopic particle interactions. The Navier-Stokes equations follow with approximations of low fluid velocities and velocity gradients. The simplicity of the cellular automaton model in fact makes it possible to derive in addition next order corrections to these equations. The derivation of hydrodynamic behavior from microscopic dynamics has never been entirely rigorous. Cellular automata can be considered as providing a simple example in which the necessary assumptions and approximations can be studied in detail. But strong support for the conclusions comes from explicit simulations of cellular automaton fluid models and the comparison of results with those from actual experiments. The next paper in this series will present many such simulations. The cellular automaton method of this paper can potentially be applied to a wide variety of processes conventionally described by partial differential equations.

Cellular Automaton Fluids

523

One example is diffusion. At a microscopic level, diffusion arises from random particle motions. The cellular automata used above can potentially reproduce diffusion phenomena, as discussed in Section 4.10. But much simpler cellular automaton rules should suffice. The derivation of the diffusion equation requires that the number of particles be conserved. But it is not necessary for total particle momentum to be conserved. Instead, particle directions should be randomized at each site. Such randomization can potentially be achieved by very simple cellular automaton rules, such as that of Ref. 20. Thus, one may devise cellular automaton methods for the solution of the diffusion equation, (56) which in turn gives a relaxation method for solving Laplace, Poisson, and related equations. Whenever the physical basis for partial differential equations involves large numbers of particles or other components with local interactions, one can expect to derive an effective cellular automaton model. For systems such as electromagnetic or gravitational fields, such models can perhaps be obtained as analogues of lattice gauge theories.

APPENDIX:

SMP

PROGRAMS

This Appendix contains a sample SMP (17) computation of the macroscopic equations for the hexagonal lattice cellular automaton fluid model of Section 2. The SMP definitions are as follows:

/* t w o - d i m e n s i o n a l d:2 /* d e f i n e p o s i t i o n r:{x,y} u:{ux,ay} case */

and v e l o c i t y

vectors

*/

/* g e n e r a t e p o l y g o n a l set of l a t t i c e v e c t o r s ~/ <XTrig p o l y g o n [ $ n ] :: ( e : A r [ $ n , { C o s [ 2 P i $ / $ n ] , S i n [ 2 P i $/$n]}]) /* c a l c u l a t e t e r m s in n u m b e r density, suma[$x] :: E x [ S u m [ $ x , { a , l , L e n [ e ] } ] ] nterm[$f] :: s u m a [ $ f [ a ] ] u t e r m [ $ f ] :: suma[e[a] $f[a]] p i t e r m [ $ f ] :: s u m a [ e [ a ] * * e [ a ] $f[a]] momentum vector and stress t e n s o r */

/* d e f i n e vector analysis operators */ e g r a d [ $ x , $ a ] :: S u m [ e [ $ a ] [ i ] D t [ $ x , r [ i ] ] , { i , l , d } ] div[$x] :: S u m [ D t [ $ x [ i ] , r [ i ] ] , { i , l , d } ] /* t e r m s in C h a p m a n - E n s k o g e x p a n s i o n */ n : f Len[e] ce0[$a] : f cel[$a] : f e [ $ a ] . u ce2[$a] : f ( ( e [ $ a ] . u ) ^ 2 - u.u/2) ce2d[$a] : f ( e g r a d [ e [ $ a ] . u , $ a ] - div[u]/2} celist : { c e 0 , c e l , c e 2 , c e 2 d }

524

Wolfram

/* s p e c i f y commutativity of Dt[$f,$1,{$2=(Ord[$2,$1]>0),l}]

second

derivatives */ :: D t [ $ f , $ 2 , $ 1 ]

/* d e f i n e printing of derivatives */ Dt[Pr][[$1,$2]]::Fmt[{{0,0},{l,-l},{2,0}},D,$2,$1]

Dt[Pr][[$1,$2,{$3,1}]]::Fmt[{{0,0],{l,-l},{2,-l},{3,0}},m,$2,$3,$1]

The following is a transcript of an interactive SMP session:

#1[l] :

#I[2] :

<"cafluid.smp" polygon[6]

1/2 3

/* /* set

load up

1/2 3

definitions hexagonal

*/ lattice

1/2 3 3

for

*/

1/2

#0[2] :

{{i/2,----},{-i/2,----},{-I,0},{-I/2,2 2 Map[nterm, celist] /* find terms

..... },{I/2,- ..... },{i,0}} 2 2 number density expansion */ from

#I [3] :

contributions to in Chapman-Enskog

#0[3] : #I[4] : #0[4] : #I[5] : # 0 [ 5 ] :*

{6f,0,0,0} Map[uterm, {{0,0},{3f Map[piterm, celist] ux,3f celist] /* find constributions to momentum vector */

uy},{0,0),{0,0}} /* stress tensor */

{{{3f,0},{0,3f}},{{0,0},{0,0}}, 2 2 2 3f ux 3f uy 3f ux uy 3f ux uy -3f ux {{ . . . . . . . . . . . . . , . . . . . . . . },{ . . . . . . . . , . . . . . . . + 4 4 2 2 3f D uy x + . . . . . . . }, 4 3f D uy y . . . . . . . }}} 4 approximation momentum equation D uy) x */ */ 4 2 3f uy . . . . . . }}, 4

3f D ux 3f D uy 3f D ux x y y {{ . . . . . . . . . . . . . . . ,....... 4 4 3f D ux y {....... 4 4

3f D uy -3f D ux x x + ....... ,........ 4 4

+

#I[6]::

#I[7]::

Dt[f,$$]:0

;

/* m a k e /* 3f

incompressibility contributions to

Fac[Map[div,@5]]

#O[7]:*

{{0,0],{0,0],

(ux D u x + u x D u y + u y D u x - u y { x y y - ......................................... 2 -3f

,

(ux D u x - u x D u y - u y D u x - u y y x x ........................................... 2 3 f (D ux + D ux) 3 f (D uy + D uy) xx yy xx yy {.................. ,.................. }} 4 4

D uy) y },

Cellular Automaton Fluids

525

ACKN OWLEDG M ENTS

Many people have contributed in various ways to the material presented here. For general discussions I thank Uriel Frisch, Brosl Hasslacher, David Levermore, Steve Orszag, Yves Pomeau, and Victor Yakhot. For specific suggestions I thank Roger Dashen, Dominique d'Humieres, Leo Kadanoff, Paul Martin, John Milnor, Steve Omohundro, Paul Steinhardt, and Larry Yaffe. Most of the calculations described here were made possible by using the SMP general-purpose computer mathematics system/iT) I thank Thinking Machines Corporation for much encouragement and partial support of this work.

REFERENCES

1. S. Wolfram, ed., Theory and Applications of Cellular Automata (World Scientific, !986). 2. S. Wolfram, Cellular automata as models of complexity, Nature 311:419 (1984). 3. N. Packard and S. Wolfram, Two-dimensional cellular automata, J. Stat. Phys. 38:901 (1985). 4. D. J. Tritton, Physical Fluid Dynamics (Van Nostrand, 1977). 5. W. W. Wood, Computer studies on fluid systems of hard-core particles, in Fundamental Problems in Statistical Mechanics 3 E. D. G. Cohen, ed. (North-Holland, 1975). 6. J. Hardy, Y. Pomeau, and O. de Pazzis, Time evolution of a two-dimensional model system. I. Invariant states and time correlation functions, J. Math. Phys. 14:1746 (1973); J. Hardy, O. de Pazzis, and Y. Pomeau, Molecular dynamics of a classical lattice gas: Transport properties and time correlation functions, Phys. Rev. A 13:1949 (1976). 7. U. Frisch, B. Hasslacher, and Y. Pomeau, Lattice gas automata for the Navie~Stokes equation, Phys. Rev. Lett. 56:1505 (1986). 8. J. Salem and S. Wolfram, Thermodynamics and hydrodynamics with cellular automata, in Theory and Applications of Cellular Automata, S. Wolfram, ed. (World Scientific, 1986). 9. D. d'Humieres, P. Lallemand, and T. Shimomura, An experimental study of lattice gas hydrodynamics, Los Alamos preprint LA-UR-85-4051; D. d'Humieres, Y. Pomeau, and P. Lallemand, Simulation d'allees de Von Karman bidimensionnelles a l'aide d'un gaz sur reseau, C. R. Acad. Sci. Paris II 301:1391 (1985). I0. J. Broadwell, Shock structure in a simple discrete velocity gas, Phys. Fluids 7:1243 (1964). 11. H. Cabannes0 The discrete Boltzmann equation, Lecture Notes, Berkeley (1980). 12. R. Gatignol, Theorie einetique des gaz a repartition discrete de vitesse (Springer, 1975). 13. J. Hardy and Y. Pomeau, Thermodynamics and hydrodynamics for a modeled fluid, or. Math. Phys. 13:1042 (1972). 14. S. Harris, The Bohzmann Equation (Holt, Rinehart and Winston, 1971). 15. R. Caflisch and G. Papanicolaou, The fluid-dynamical limit of a nonlinear model Boltzmann equation, Commun. Pure Appl. Math. 32:589 (1979). 16. B. Nemnich and S. Wolfram, Cellular automaton fluids 2: Basic phenomenology, in preparation. 17. S. Wolfram, S M P Referenee Manual (Inference Corporation, Los Angeles, 1983); S. Wolfram, Symbolic mathematical computation, Commun. A C M 28:390 (1985). 18. A. Sommerfeld, Thermodynamics and Statistical Mechanics (Academic Press, 1955). 19. S. Wolfram, Origins of randomness in physical systems, Phys. Rev. Lett. 55:449 (1985). 20. S. Wolfram, Random sequence generation by cellular automata, Adv. Appl. Math. 7:123 (1986). 21. J. P. Boon and S. Yip, Molecular Hydrodynamics (McGraw-Hill, 1980).

526

Wolfram

22. E. M. Lifshitz and L. P. Pitaevskii, Statistical Mechanics, Part 2 (Pergamon, 1980), Chapter 9. 23. R. Liboff, The Theory of Kinetic Equations (Wiley, 1969). 24. D. Levermore, Discretization effects in the macroscopic properties of cellular automaton fluids, in preparation. 25. E. M. Lifshitz and L. P. Pitaevskii, Physical Kinetics (Pergamon, 1981). 26. P. Resibois and M. De Leener, Classical Kinetic Theory of Fluids (Wiley, 1977). 27. L. D. Landau and E. M. Lifshitz, Fluid Mechanics (Pergamon, 1959). 28. M. H. Ernst, B. Cichocki, J. R. Dorfman, J. Sharma, and H. van Beijeren, Kinetic theory of nonlinear viscous flow in two and three dimensions, J. Stat. Phys. 18:237 (1978). 29. J. R. Dorfman, Kinetic and hydrodynamic theory of time correlation functions, in Fundamental Problems in Statistical Mechanics 3, E. D. G. Cohen, ed. (North-Holland, 1975). 30. R. Courant and K. O. Friedrichs, Supersonic Flows and Shock Waves (Interscience, 1948). 31. D. Levermore, private communication. 32. J. Milnor, private communication. 33. V. Yakhot, B. Bayley, and S. Orszag, Analogy between hyperscale transport and cellular automaton fluid dynamics, Princeton University preprint (February 1986). 34. H. S. M. Coxeter, Regular Polytopes (Macmillan, t963). 35. M. Hammermesh, Group Theory (Addison=Wesley, 1962), Chapter 9. 36. L. D. Landau and E. M. Lifshitz, Quantum Mechanics (Pergamon, 1977), Chapter 12. 37. H. Boerner, Representations of Groups (North-Holland, 1970), Chapter 7. 38. L. D. Landau and E. M. Lifshitz, Theory of Elasticity (Pergamon, 1975), Section 10. 39. D. Levine et al., Elasticity and dislocations in pentagonal and icosahedral quasicrystals, Phys. Rev. Lett. 14:1520 (1985). 40. L. D. Landau and E. M. Lifshitz, Statistical Physics (Pergamon, 1978), Chapter 13. 41. B. K. Vainshtein, Modern Crystallography, (Springer, 1981), Chapter 2. 42. J. H. Conway and N. J. A. Sloane, to be published. 43. R. L. E. Schwarzenberger, N-Dimensional Crystallography (Pitman, 1980). 44. J. Milnor, Hilbert's problem 18: On crystallographic groups, fundamental domains, and on sphere packing, Proc. Symp. Pure Math. 28:491 (1976). 45. B. G. Wybourne, Classical Groups for Physicists (Wiley, 1974), p. 78; R. Slansky, Group theory for unified model building, Phys. Rep. 79:1 (1981). 46. B. Grunbaum and G. C. Shephard, Tilings and Patterns (Freeman, in press); D. Levine and P. Steinhardt, Qnasicrystals I: Definition and structure, Univ. of Pennsylvania preprint. 47. N. G. de Bruijn, Algebraic theory of Penrose's non-periodic tilings of the plane, Nedl. Akad. Wetenseh. Indag. Math. 43:39 (1981); J. Socolar, P. Steinhardt, and D. Levine, Quasicrystals with arbitrary orientational symmetry, Phys. Rev. B 32:5547 (1985). 48. R. Penrose, Pentaplexity: A class of nonperiodic tilings of the plane, Math. Intelligencer 2:32 (1979). 49. J. P. Rivet and U. Frisch, Automates sur gaz de reseau dans l'approximation de Boltzmann, C. R. Acad. Sci. Paris H 302:267 (1986). 50. P. J, Davis, Circulant Matrices (Wiley, 1979). 51. L. D. Landau and E. M. Lifshitz, Statistical Physics (Pergamon, 1978), Chapter 5. 52. E. Kolb and S. Wolfram, Baryon number generation in the early universe, Nucl. Phys. B 172:224 (1980), Appendix A. 53. I. S. Gradshteyn and I. M. Ryzhik, Table of integrals, Series and Products (Academic Press, 1965). 54. U. Frisch, private communication. 55. P. Roache, Computational Fluid Mechanics (Hermosa, Albuquerque, 1976). 56. S. Omohundro and S. Wolfram, unpublished (July 1985). 57. D. d'Humieres, private communication.

相关文章:

更多相关标签:

- CELLULAR AUTOMATA
- Cellular Automata
- 中国页岩气勘探开发(中石油内部资料.英文)
- 影响页岩气产能的关键因素(英文)
- RS232与RS485串行接口转换电路及其编程实现
- Basic Action Theory
- basic_laser_theory
- cellular-automaton-fluids-theory
- Cellular Automata-Basic Intro
- Chapter 1 Basic theory
- Evacuation Model for Delayed Flights Based on Cellular Automation
- Basic Concepts of the Theory of Sets
- Random walk theory of jamming in a cellular automaton model for traffic flow
- Basic TEFL TESOL Theory
- steady state velocity distribution analysis for a combined cellular automation traffic model