Stephen Wolframstephenwolfram.com
Publications by Stephen Wolfram * Articles * Cellular Automata * Cellular Automaton Fluids: Basic Theory (1986)
Cellular Automaton Fluids: Basic Theory (1986)


2. Macroscopic Equations for a Sample Model

2.1. Structure of the Model

The model is based on a regular two-dimensional lattice of hexagonal cells, as illustrated in Fig. 1. The site at the center of each cell is connected to its six neighbors by links corresponding to the unit vectors through given by

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 . 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 . 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 , as would triples of particles 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.

2.2. Basis for Kinetic Theory

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 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 configurations, 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. 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. (2) 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 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. Kinetic Equations

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 gives the probability of finding a particle with velocity at position and time , averaged over all other features of the configuration (e.g., Ref. 23).

Two processes lead to changes in 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 with velocity move at the next time step to the adjacent cell at position . As a result, the distribution function evolves according to

For large lattices and long time intervals, position and time may be approximated by continuous variables. One may define, for example, scaled variables and , where . In terms of these scaled variables, the difference equation (2.3.1) becomes

In deriving macroscopic transport equations, this must be converted to a differential equation. Carrying out a Taylor expansion, one obtains

If all variations in the are assumed small, and certainly less than , it suffices to keep only first-order terms in . In this way one obtains the basic transport equation

This has the form of a collisionless Boltzmann transport equation for (e.g., Ref. 25). It implies, as expected, that is unaffected by particle motion in a spatially uniform system.

Collisions can, however, change even in a uniform system, and their effect can be complicated. Consider, for example, collisions that cause particles in directions and to scatter in directions and . The rate for such collisions is determined by the probability that particles in directions and occur together in a particular cell. This probability is defined as the joint two-particle distribution function . The collisions deplete the population of particles in direction at a rate . Microscopic reversibility guarantees the existence of an inverse process, which increases the population of particles in direction at a rate given in this case by . Notice that in a model where there can be at most one particle on each link, the scattering to directions and in a particular cell can occur only if no particles are already present on these links. The distribution function 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 will for now be summarized by a simple ``collision term'' , which in general depends on two-particle and higher order distribution functions. (In the models considered here, is always entirely local, and cannot depend directly on, for example, derivatives of distribution functions.) In terms of , the kinetic equation (2.3.3) extended to include collisions becomes

With the appropriate form for , this is an exact statistical equation for (at least to first order in ).

But the equation is not in general sufficient to determine . It gives the time evolution of in terms of the two-particle and higher order distribution functions that appear in . 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, of which Eq. (2.3.5) is the first level.

The Boltzmann transport equation approximates (2.3.5) by assuming that 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 is thus approximated as . 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 is an approximation, and corrections can arise, as discussed in Section 2.5. 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.

2.4. Conservation Laws

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

while the momentum density , where is the average fluid velocity, is given by

The conservation of these quantities places important constraints on the behavior of the .

In a uniform system , so that Eq. (2.3.5) becomes

and Eqs. (2.4.1) and (2.4.2) imply

Using the kinetic equation (2.3.5), Eq. (2.4.4) implies

With the second term in the form , Eq. (2.4.6) can be written exactly in terms of macroscopic quantities as

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

Defining the momentum flux density tensor

Eq. (2.4.8) becomes

No simple macroscopic result for 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 . 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.

Assuming , Eq. (2.4.6) to second order becomes

Writing the term in the form

this term is seen to vanish for any which 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.

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

The second term vanishes if 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).

2.5. Chapman-Enskog Expansion

If there is local equilibrium, as discussed in Section 2.2, then the microscopic distribution functions should depend, on average, only on the macroscopic parameters and and their derivatives. In general, this dependence may be very complicated. But in hydrodynamic processes, and vary only slowly with position and time. In addition, in the subsonic limit, .

With these assumptions, one may approximate the by a series or Chapman-Enskog expansion in the macroscopic variables. To the order required for standard hydrodynamic phenomena, the possible terms are

where the 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 first-order 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 from the vectors , and .

The relation

where here and , and and are space indices, has been used in Eq. (2.5.1) to choose the forms of the and terms so as to satisfy the constraints (2.4.1) and (2.4.2), independent of the values of the coefficients and . In terms of (2.5.1), Eq. (2.4.1) yields immediately

while (2.4.2) gives

The specific values of and 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 , all the are given by

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

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:

When , can be evaluated in the approximation (2.5.1) using the relations

and

The result is

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

where

The form (2.5.10) for 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 . These corrections must be compared with other effects included in Eq. (2.5.11). The rescaling implies that spatial gradient terms in the Chapman-Enskog expansion can be of the same order as the correction terms in Eq. (2.4.13). When the term in the Chapman-Enskog expansion (2.5.1) for the is substituted into the last term of Eq. (2.4.13), it gives a contribution

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

2.6. Navier-Stokes Equation

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

where is pressure, and 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 proportional to be traceless.

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 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 must also be neglected.

The term proportional to remains, but can be combined with the 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 a cellular automaton fluid, considered here, bulk viscosity is given by

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

where is the kinematic viscosity. An approximate method of evaluating 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

which is not in general equal to 1. In continuum fluids, the covariant derivative usually has the form 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.

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

the equation for has coefficient 1 for the 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

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 and . Correlation function methods indicate, however, that additional effects yield logarithmically divergent contributions to in two dimensions (e.g., Ref. 29). The full viscous term in this case may in fact be of the rough form .

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 . 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 , 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.

In the transonic flow regime , 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 and higher in the Chapman-Enskog expansion (2.5.1).

In the homogeneous fluid approximation , one may take

The constraints (2.4.1) and (2.4.2) imply

where 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,

The term in Eq. (2.7.6) has the isotropic form given in Eq. (2.5.11). The term is, however, anisotropic.

To obtain an isotropic term, one must generalize the model, as discussed in Section 3. One possibility is to allow vectors corresponding to corners of an -sided polygon with . In this case, the continuum equation deduced from the Chapman-Enskog expansion (2.7.1) becomes

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 in the Chapman-Enskog expansion (2.7.1). The possible fourth-order terms are given by contractions of with products of or . They yield a piece in the Chapman-Enskog expansion of the form

where Eq. (2.4.1) implies the constraints (for )

The resulting continuum equations may be written in terms of vectors formed by contractions of and . The complete result is

where, on the right-hand side, the first group of terms are all , while the second group are . Further corrections involve higher derivative terms, such as .

For a channel flow with , , the time-independent terms in Eq. (2.7.11) have an component

and zero component.

previous  l  next