Main Text

Thermodynamics and hydrodynamics describe the overall behaviour of many systems, independent of the precise microscopic construction of each system. One can thus study thermodynamics and hydrodynamics using simple models, which are more amenable to efficient simulation, and potentially to mathematical analysis.

Cellular automata (CA) are discrete dynamical systems which give simple models for many complex physical processes [1]. This paper considers CA which can be viewed as discrete approximations to molecular dynamics. In the simplest case, each link in a regular spatial lattice carries at most one ``particle'' with unit velocity in each direction. At each time step, each particle moves one link; those arriving at a particular site then ``scatter'' according to a fixed set of rules. This discrete system is well-suited to simulation on digital computers. The state of each site is represented by a few bits, and follows simple logical rules. The rules are local, so that many sites can be updated in parallel. The simulations in this paper were performed on a Connection Machine Computer [2] which updates sites concurrently in each of 65536 Boolean processors [3].

In two dimensions, one can consider square and hexagonal (six links at ) lattices. On a square lattice [4], the only nontrivial local rule which conserves momentum and particle number takes isolated pairs of particles colliding head on to scatter in the orthogonal direction (no interaction in other cases). On a hexagonal lattice [5], such pairs may scatter in either of the other two directions, and the scattering may be affected by particles in the third direction. Four particles coming along two directions may also scatter in different directions. Finally, particles on three links separated by may scatter along the other three links. At fixed boundaries, particles may either ``bounce back'' (yielding ``no slip'' on average), or reflect ``specularly'' through .

On a microscopic scale, these rules are deterministic, reversible and discrete. But on a sufficiently large scale, a statistical description may apply, and the system may behave like a continuum fluid, with macroscopic quantities, such as hydrodynamic velocity, obtained by kinetic theory averages.

Figure 1 illustrates relaxation to ``thermodynamic equilibrium''. The system randomizes, and coarse-grained entropy increases. This macroscopic behaviour is robust, but microscopic details depend sensitively on initial conditions. Small perturbations (say of one particle) have microscopic effects over linearly-expanding regions [6]. Thus ensembles of ``nearby'' initial states usually evolve to contain widely-differing ``typical'' states. But in addition, individual ``simply-specified'' initial states can yield behaviour so complex as to seem random [7,8], as in figure 1. The dynamics thus ``encrypts'' the initial data; given only coarse-grained, partial, information, the initial simplicity cannot be recovered or recognized by computationally feasible procedures [7], and the behaviour is effectively irreversible.

Microscopic instability implies that predictions of detailed behaviour are impossible without ever more extensive knowledge of initial conditions. With complete knowledge (say from a simple specification), the behaviour can always be reproduced by explicit simulation. But if effective predictions are to be made, more efficient computational procedures should be found. The CA considered here can in fact act as universal computers [9]: with appropriate initial conditions, their evolution can implement any computation. Streams of particles corresponding to ``wires'' can meet in logical gates implemented by fixed obstructions or other streams. As a consequence, the evolution is computationally irreducible [10]; there is no general shortcut to explicit simulation. No simpler computation can reproduce all the possible phenomena.

Some overall statistical predictions can nevertheless be made. In isolation, the CA seem to relax to an equilibrium in which links are populated effectively randomly with a particular average particle density and net velocity (as in figure 1). On length scales large compared to the mean free path , the system then behaves like a continuum fluid. The effective fluid pressure is , giving a speed of sound . Despite the microscopic anisotropy of the lattice, circular sound wavefronts are obtained from point sources (so long as their wavelength is larger than the mean free path) [11].

Assuming local equilibrium, the large-scale behaviour of the CA can be approximated by average rules for collections of particles, with particular average densities and velocities. The rules are like finite difference approximations to partial differential equations, whose form can be found by a standard Chapman-Enskog expansion [12] of microscopic particle distributions in terms of macroscopic quantities. The results are analogous to those for systems [13] in which particles occur with an arbitrary continuous density at each point in space, but have only a finite set of possible velocities corresponding to the links of the lattice. The hexagonal lattice CA is then found to follow exactly the standard Navier-Stokes equations [5,14]. As usual, the parameters in the Navier-Stokes equations depend on the microscopic structure of the system. Kinetic theory suggests a kinematic viscosity [15].

[ Figure 1 ] Relaxation to ``thermodynamic equilibrium'' in the hexagonal lattice cellular automaton (CA) described in the text. Discrete particles are initially in a simple array in the centre of a site square box. The upper sequence shows the randomization of this pattern with time; the lower sequence shows the cells visited in the discrete phase space (one particle track is drawn thicker). The graph illustrates the resulting increase of coarse-grained entropy calculated from particle densities in regions of a box.

Figures 2 and 3 show hydrodynamic phenomena in the large scale behaviour of the hexagonal lattice CA. An overall flow is obtained by maintaining a difference in the numbers of left- and right-moving particles at the boundaries. Since local equilibrium is rapidly reached from almost any state, the results are insensitive to the precise arrangement used. Random boundary fluxes imitate an infinite region; a regular pattern of incoming particles nevertheless also suffices, and reflecting or cyclic boundary conditions can be used on the top and bottom edges.

The hydrodynamics of the CA is much like a standard physical fluid [16]. For low Mach numbers , the fluid is approximately incompressible, and the flows show dynamical similarity, depending only on Reynolds number (). The patterns obtained agree qualitatively with experiment [3]. At low , the flows are macroscopically stable; perturbations are dissipated into microscopic ``heat''.

[ Figure 2 ] Time evolution of hydrodynamic flow around a plate in the CA of figure 1 on a site lattice. Hydrodynamic velocities are obtained as indicated by averaging over site regions. There is an average density of 0.3 particles per link (giving a total of particles). An overall velocity is maintained by introducing an excess of particles (here in a regular pattern) on the left hand boundary.

[ Figure 3 ] Hydrodynamic flows obtained after time steps in the CA of figure 2, for various overall velocities .

As increases, periodic vortex streets are at first produced, and then vortices are shed in an irregular, turbulent, fashion. Perturbations now affect details of the flow, though not its statistical properties. The macroscopic irregularity does not depend on microscopic randomness; it occurs even if microscopically simple (say spatially and temporally periodic) initial and boundary conditions are used, as illustrated in figure 2. As at the microscopic level, it seems that the evolution corresponds to a sufficiently complex computation that its results seem random [7].

The CA discussed here should serve as a basis for practical hydrodynamic simulations. They are simple to program, readily amenable to parallel processing, able to handle complex geometries easily [17], and presumably show no unphysical instabilities. (Generalization to three dimensions is straightforward in principle [18].)

Standard finite difference methods [19] consider discrete cells of fluid described by continuous parameters. These parameters are usually represented as digital numbers with say 64 bits of precision. Most of these bits are, however, probably irrelevant in determining observable features of flow. In the CA approach, all bits are of essentially equal importance, and the number of elementary operations performed is potentially closer to the irreducible limit.

The difficulty of computation in a particular case depends on the number of cells that must be used. Below a certain dissipation length scale (in dimensions), viscosity makes physical homogeneous turbulent fluids smooth [16]. In finite difference schemes, individual cells can represent fluid regions of this size. But complete calculations with the CA considered here probably require increasing numbers of cells in each region [20]. Approximate ``turbulence models'' involving fewer cells may however be devised.

Several further extensions of the CA scheme can be considered. First, on some or all of the lattice, basic units containing say particles, rather than single particles, can be used. The properties of these units can be specified by digital numbers with bits, but exact conservation laws can still be maintained. This scheme comes closer to adaptive grid finite difference methods [19], and potentially avoids detailed computation in featureless parts of flows.

A second, related, extension introduces discrete internal degrees of freedom for each particle. These could represent different particle types, directions of discrete vortices [19], or internal energy (giving variable temperature [21]).

This paper has given further evidence that simple cellular automata can reproduce the essential features of thermodynamic and hydrodynamic behaviour. These models make contact with results in dynamical systems theory and computation theory. They should also yield efficient practical simulations, particularly on parallel-processing computers.

Cellular automata can potentially reproduce behaviour conventionally described by partial differential equations in many other systems whose intrinsic dynamics involves many degrees of freedom with no large disparity in scales.

We are grateful to U. Frisch, B. Hasslacher, Y. Pomeau and T. Shimomura for sharing their unpublished results with us, and to N. Margolus, S. Omohundro, S. Orszag, N. Packard, R. Shaw, T. Toffoli, G. Vichniac and V. Yakhot for discussions. We thank many people at Thinking Machines Corporation for their help and encouragement. The work of S.W. was supported in part by the U.S. Office of Naval Research under contract number N00014-85-K-0045.

previous  l  next