This document discusses how simple cellular automata can capture the essential features of thermodynamics and hydrodynamics. At the microscopic level, the cellular automata behave like discrete approximations of molecular dynamics and relax towards equilibrium. On a large scale, they behave like continuum fluids and can efficiently simulate hydrodynamic phenomena. The cellular automata follow the Navier-Stokes equations and exhibit phenomena like vortex streets and turbulent behavior in a manner similar to physical fluids. They provide a simple, parallelizable approach for practical hydrodynamic simulations.
This document discusses how simple cellular automata can capture the essential features of thermodynamics and hydrodynamics. At the microscopic level, the cellular automata behave like discrete approximations of molecular dynamics and relax towards equilibrium. On a large scale, they behave like continuum fluids and can efficiently simulate hydrodynamic phenomena. The cellular automata follow the Navier-Stokes equations and exhibit phenomena like vortex streets and turbulent behavior in a manner similar to physical fluids. They provide a simple, parallelizable approach for practical hydrodynamic simulations.
This document discusses how simple cellular automata can capture the essential features of thermodynamics and hydrodynamics. At the microscopic level, the cellular automata behave like discrete approximations of molecular dynamics and relax towards equilibrium. On a large scale, they behave like continuum fluids and can efficiently simulate hydrodynamic phenomena. The cellular automata follow the Navier-Stokes equations and exhibit phenomena like vortex streets and turbulent behavior in a manner similar to physical fluids. They provide a simple, parallelizable approach for practical hydrodynamic simulations.
This document discusses how simple cellular automata can capture the essential features of thermodynamics and hydrodynamics. At the microscopic level, the cellular automata behave like discrete approximations of molecular dynamics and relax towards equilibrium. On a large scale, they behave like continuum fluids and can efficiently simulate hydrodynamic phenomena. The cellular automata follow the Navier-Stokes equations and exhibit phenomena like vortex streets and turbulent behavior in a manner similar to physical fluids. They provide a simple, parallelizable approach for practical hydrodynamic simulations.
Thermodynamics and Hydrodynamics with Cellular Automata
James B. Salem Thinking Machines Corporation, 245 First Street, Cambridge, MA 02144
and
Stephen Wolfram The Institute for Advanced Study, Princeton NJ 08540.
(November 1985)
Simple cellular automata which seem to capture the essential features of thenno- dynamics and hydrodynamics are discussed. At a microscopic level, the cellular auto- mata are discrete approximations to molecular dynamics, and show relaxation towards equilibrium. On a large scale, they behave like continuum fiuids, and suggest efficient methods for hydrodynamic simulation.
Thermodynamics and hydrodynamics describe the overall behaviour of many systems, indepen- dent 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 com- plex 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 60°) lattices. On a square lattice [4], the only nontrivial rule which consenres momentum and particle number takes isolated pairs of particles colliding head on 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 corning along two directions may also scatter in different directions. Finally, particles on three links separated by 120° may scatter along the 2
other three links. At fixed boundaries, particles may either "bounce back" (yielding "no slip" on aver- age), or reflect "specularly" through 120°. 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 recog- nized 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 predic- tions are to be made, more efficient computational procedures should be found. It is known, however, that the CA considered here can 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 compu- tation 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 p and net velocity (as in figure 1). On length scales large compared to the mean free path A., the system then behaves like a continuum fluid. The effective fluid pressure is p=p/2, giving a speed of sound c=1I..fi.. 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 quanti- ties. The results are analogous to those for systems [13] in which particles occur with an arbitrary con- tinuous 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 micros- copic structure of the system. Kinetic theory suggests a kinematic viscosity v=AJ2 [15]. Figures 2 and 3 show hydrodynamic phenomena in the large scale behaviour of the hexagonal lat- tice CA. An overall flow U 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 [14]. 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 Ma=Ulc, the fluid is approximately incompressible, and the flows show dynamical similarity, depending only on Reynolds number Re::.UUv (1.>)..). The patterns obtained agree qualitatively with experiment [3]. At low Re, the flows are macroscopically stable; perturbations are dissipated into microscopic "heat". At higher Re, vortex streets are produced, and there are signs of irregular, turbulent, behaviour. Perturbations now affect details of the flow, though not its statistical properties. The macroscopic irre- gularity does not depend on microscopic randomness; it occurs even if microscopically simple (say spa- tially and temporally periodic) initial and boundary conditions are used, as illustrated in figure 2. As at 3
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, and presumably show no unphysical instabilities. (Generalization to three dimensions is straightforward in principle [17].) Standard finite difference methods [18] 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 irrevelant in determining observable features of flow. In the CA approach, all bits are of essentially equal importance, and the number of elementary operations per- formed 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 a, viscosity makes physical fluids smooth. Each CA cell can thus potentially represent a packet of fluid on this length scale. In fully developed turbulence, a appears to scale approximately as Re-dl4 in d dimensions [16]. These estimates govern cell sizes needed in finite difference schemes. The CA discussed so far incorporate no microscopic irreversibility: macroscopic energy dissipa- tion occurs only through degradation to microscopic "heat". But in turbulence there is continual dissi- pation of -Red energy at small scales. To absorb this energy, the CA cells must be Re 1l4 (d=3) or Re l12 (d=2) times larger than a, making such CA less efficient than irreversible finite difference schemes [19]. "Turbulence models" with irreversibility in the microscopic CA rules may overcome this difficulty. Several further extensions of the CA scheme can be considered. First, on some or all of the lat- tice, basic units containing say n particles, rather than single particles, can be used. The properties of these units can be specified by digital numbers with O(logn) bits, but exact conservation laws can still be maintained. This scheme comes closer to adaptive grid finite difference methods [18], 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 [18], or internal energy (giving variable temperature [20]). This paper has given 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, particu- larly 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 free- dom with no large disparity in scales. We thank U. Frisch, B. Hasslacher, Y. Pomeau and T. Shimomura for sharing their unpublished results with us, and we thank N. Margolus, S. Omohundro, S. Orszag, N. Packard, R. Shaw, T. Toffoli, G. Vichniac and V. Yakhot for discussions. We are also grateful to 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 NOOOI4-85-K-0045.
1. See for example S. Wolfram, "Cellular automata as models of complexity", Nature 311, 419 (1984) where applications to thermodynamics and hydrodynamics were mentioned but not explored. 2. D. Hillis, The Connection Machine (MIT press, 1985). This application is discussed in S. Wol- fram, "Scientific computation with the Connection Machine", Thinking Machines Corporation report (March 1985). 3. More detailed results of theory and simulation will be given in a forthcoming series of papers. 4
4. J. Hardy, O. de Pazzis and Y. Pomeau, "Molecular dynamics of a classical lattice gas: transport properties and time correlation functions", Phys. Rev. A13, 1949 (1976). 5. U. Frisch, B. Hasslacher and Y. Pomeau, "A lattice gas automaton for the Navier-Stokes equa- tion", Los Alamos preprint LA-UR-85-3503. 6. The expansion rate gives the Lyapunov exponent as defined in N. Packard and S. Wolfram, "Two-dimensional cellular automata", 1. Stat. Phys. 38, 901 (1985). Note that the effect involves many particles, and does not arise from instability in the motion of single particles, as in the case of hard spheres with continuous position variables (e.g. O. Penrose, "Foundations of statistical mechanics", Rep. Prog. Phys. 42, 129 (1979).) 7. S. Wolfram, "Origins of randomness in physical systems", Phys. Rev. Lett. 55, 449 (1985); "Random sequence generation by cellular automata" , Adv. Appl. Math. (in press). 8. Simple patterns are obtained with very simple or symmetrical initial conditions. On a hexagonal lattice, the motion of an isolated particle in a rectangular box is described by a linear congruence relation, and is ergodic when the side lengths are not commensurate. 9. N. Margolus, "Physics-like models of computation", Physica lOD, 81 (1984). 10. S. Wolfram, "Undecidability and intractability in theoretical physics", Phys. Rev. Lett. 54, 735 (1985). 11. cf T. Toffoli, "CAM: A high-performance cellular automaton machine", Physic a lOD, 195 (1984). 12. e.g. S. Harris, The Boltzmann Equation, (Holt, Reinhart and Winston, 1971). 13. R. Gatignol, Theorie cinetique des gaz a repartition discrete de vitesse, (Springer, 1975); J. Broadwell, "Shock structure in a simple discrete velocity gas", Phys. Fluids 7, 1243 (1964). 14. On a square lattice, the total momentum in each row is separately conserved, and so cannot be convected by velocity in the orthogonal direction. Symmetric three particle collisions on a hexag- onallattice remove this spurious conservation law. 15. The rank three and four tensor coefficients of the nonlinear and viscous terms in the Navier- Stokes equations are isotropic for a hexagonal but not a square lattice. Higher order coefficients are anisotropic in both cases. In two dimensions, there can be logarithmic corrections to the Newtonian fluid approximation: these can apparently be ignored on the length scales considered, but yield a formal divergence in the viscosity (cf [3]). 16. e.g. D. J. Tritton, Physical fluid dynamics, (Van Nostrand, 1977). 17. Icosahedral symmetry yields isotropic fluid behaviour, and can be achieved to an arbitrary degree of approximation with a suitable periodic lattice (cf D. Levine et ai., " Elasticity and dislocations in pentagonal and icosahedral quasicrystals", Phys. Rev. Lett. 54, 1520 (1985); P. Bak, "Sym- metry, stability, and elastic properties of icosahedral incommensurate crystals", Phys. Rev. B32, 5764 (1985». 18. e.g. P. Roache, Computational fluid dynamics, (Hermosa, Albuquerque, 1976). 19. S. Orszag and V. Yakhot, "Reynolds number scaling of cellular automaton hydrodynamics", Princeton University Applied and Computational Math. report (November 1985). 20. In simple cases the resulting model is analogous to a deterministic microcanonical spin system (M. Creutz, "Deterministic Ising dynamics", Ann. Phys., to be published.) t=O t=100 t=200
.' , .:' : .."~...'
.' :
s 0.14
o t 1000
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 32x32 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 illus- trates the resulting increase of coarse-grained entropy Uilogz,Pi calculated from particle densities in 32x32 regions of a 256><256 box.
Figure 2. Time evolution of hydrodynamic Bow around an inclined plate in the CA of figure 1 on a 2048x2048 site lattice. Hydrodynamic velocities are obtained as indicated by averaging over 64x64 site regions. There is an average density of 0.3 particles per link (giving a total of 7xl06 particles). An overall velocity U=O.l is maintained by introducing an excess of particles (here in a regular pattern) on the left hand boundary.
