Practical Application of Lattice-Gas and Lattice Boltzmann Methods to Dispersion Problems

Abstract

We compare lattice-gas and lattice Boltzmann techniques for the quantitative study of dispersion in natural systems. The principal advantage of lattice gases is speed (computational efficiency), and the ability to dilute boundary errors by using large numbers of sites; however, with large numbers of sites, viscous equilibration time can severely limit the practicality of the method. The principal advantage of lattice Boltzmann is the tunability of the transport coefficients, allowing simulations in narrow channels; however, errors in boundary conditions appear to be more serious for dispersion phenomena, than for the widely studied bulk flow (single-component) problems, and these errors may prevent the use of very small channels. Despite these problems, agreement with experimental results is often quite good, and the experimental error for dispersion measurements will often swamp the errors inherent in the numerical methods. For modeling double-diffusive fingering (salt fingering), the lattice Boltzmann method is clearly superior.

Introduction

When lattice-gas automata (LGA) were introduced by Frisch, Hasslacher and Pomeau (1986) as a means to solve the Navier-Stokes equations, the technique was first envisioned as a "numerical wind tunnel" to study high-speed, complex flows bounded by geometrically simple solid barriers. However, Rothman (1988) noted that LGA were powerful tools for modeling "slow" flows (Reynolds number < 100) in complex porous and geological media.  The method was quickly extended to immiscible liquids (Rothman and Keller, 1988). Early models were two-dimensional, or restricted to single-component flows; but the advance of supercomputing, and the development of the lattice Boltzmann (LB) technique, allowed full 3D calculation of flow in porous media, with realistic representations of immiscible liquids in petroleum-bearing sandstones (Martys and Chen, 1996).  The adoption of the Bhatnagar-Gross-Krook (BGK) collision operator (Qian et al., 1992) simplified LB calculations, and now multi-component LB calculations of complex media are practical on single-CPU computers.

If a tracer is injected into a flowing fluid, it is dispersed -- that is, spread and stretched by a combination of flow and molecular diffusion (Aris, 1956).  For example, the tracer might be smoke in air, or a dye or solute in water. When influenced by buoyancy, multiple diffusion coefficients, or constraint in convoluted channels, the process can be remarkably complex.  LGA and LB techniques are well-suited to modeling dispersion, particularly in the intermediate Peclet number range, where neither diffusion nor fluid flow dominates the spreading of solute.  Baudet et al. (1989) were apparently the first to apply LGA techniques to dispersion in planar channels, and more complex geometries were considered by Gutfraind et al. (1995) and Perea-Reeves and Stockman (1997).  The principal advantage of LGA and LB is that they handle complicated boundary and initial conditions very well, free from many "gridding" and stability constraints that plague conventional numerical methods.

In this note, we examine several applications of LGA and LB to dispersion problems in the natural sciences, and compare the relative advantages of each technique.  We focus on oft-neglected, practical aspects of the calculations: errors associated with boundary conditions, accuracy required for meaningful comparisons with experimental data, computer coding, and the factors that dictate problem size and run time.  All the examples are 2D or quasi-3D, and can be run on single-CPU personal computers or workstations.

Background

LGA methods.  LGA are collections of discrete particles constrained to move on fixed, geometric lattices. Figure 1 shows a "microscopic" view of the stages in a single LGA time step, for particles near a solid wall; figure 2 follows the process through many time steps. There is a translation or "streaming" step in which particles move from site to site on the grid, then a collision step.  Both figures depict the 2D hexagonal "FHP" lattice (Frisch et al., 1986), which has only 6 possible particle velocities, and one particle speed.  Collisions in the fluid conserve momentum and particle number; collisions with walls are typically "bounce-back", that is, particle directions are simply reversed.  In modeling multi-component systems, the particles are often described as having "colors" that follow the particles and give rise to diffusion.  In computer memory, the particles are represented by individual bits; the translation and collision steps are carried out entirely with simple integer operations.

Figure 1.  Three stages of a single time step for a 2D hexagonal lattice near a solid wall (brick pattern).  This is a "microscopic" view, containing 11 particles over ca. 15 fluid sites and 5 wall sites.  A typical simulation uses 106 sites and 2·106 particles.  Collisions in the fluid conserve momentum and particle number.  At the walls, particles undergo bounce-back collisions, which simply reverse particle direction. Click on thumbnail for a larger image.

mpeg movie [571 kB]

Figure 2.  A movie of particle motions in 8x8 sites near a solid object (denoted by green dots).  This simulation was part of a much larger calculation with two components, represented by the blue- and red-colored arrows, so particles can leave and enter from outside the viewfield.  You may wish to single-step through the movie.

Momentum and particle averages for these simple systems, taken over many sites or many time steps, obey the equations:

        (1)

and

                                                               (2)

where t is time,  is the fluid velocity (u_x , u_y , u_z),  is the density in particles/site, P is pressure,  is a body force (e.g. gravity),  is the kinematic viscosity, C is the concentration of a tracer dispersed through the system, and D_m is the molecular diffusion coefficient.  Equation (1) is very similar to the Navier-Stokes equation, and equation (2) is the advection-diffusion equation.  For steady-state flow, the  term can be adsorbed into the viscosity, pressure and force terms, making equation (1) effectively equivalent to Navier-Stokes, subject to constraints discussed below.  The inclusion of simple rules for nearest-neighbor interactions (Rothman and Keller, 1988) gives rise to immiscibility that obeys the LaPlace equation for bubble pressure, and Langmuir sorption and buoyancy are added easily.

The following mpeg movies illustrate the complexity of problems that can be handled by LGA; you may wish to set your mpeg viewer to "loop" mode.

Apart from generating "pretty pictures" of physical phenomena, LGA can also be used for quantitative tests and predictions.  They are particularly useful for investigating the effects of buoyancy and scaling on dispersion in small-scale laboratory tests, as aides to experimental design and interpretation. While LB methods have nearly supplanted LGA in the literature, LGA are still useful for a surprising number of applications, and have the advantages of rigorous particle conservation and high computational speed.

As noted, LGA calculations require "averaging" to recover the macroscopic equations (1) and (2).  In general, the averaging process poses a much more severe constraint on velocity field calculations than on dispersion modeling.  For example, we found it necessary to average over 10⁵ steps to resolve recirculation zones at Reynolds numbers (Re) of 50, in the modeling of flow through rough fractures (Brown et al., 1995).  With such long averaging times, the LGA technique is often restricted to slowly varying or steady-state flows.  Even with tight coding (e.g., assembly language), the averaging process can be the most time-consuming part of a flow calculation; this limitation is rarely mentioned, and must temper enthusiastic reports of LGA computational performance.  On the other extreme, Perea-Reeves and Stockman (1997) found concentration averages over just 4 to 64 steps were adequate for dispersion measurements at a Peclet number (Pe) of 50.

Lattice Boltzmann methods.  The lattice Boltzmann technique is described in detail by Martys and Chen (1996). Here we consider only LB methods with the simple BGK collision operator (Qian et al., 1992).  The principal difference between LGA and LB is that the latter uses probabilistic particle distributions, rather than discrete particles; thus it is a floating-point method. The LB equation is written:

         (3)

where the f_i represent the particle distribution functions for each vector direction i, t is time,  is an (x,y,z) position in the lattice,  is the velocity in the ith direction,  f_i^eq  is the equilibrium particle distribution, and  is a relaxation parameter.  The left side of the equation is analogous to the "translation" stage in LGA, and the right to the "collision" stage. For example, in the two-dimensional "D2Q9" model of Qian et al., there are 9 velocities () on a square lattice: one has speed=0 and corresponds to a "rest" particle; four have speed=1 and are at 0, 90, 180 and 270 degrees; and four have speed= at 45, 135, 225 and 315 degrees (figure 7).

Figure 7.  Schematic 2D lattice Boltzmann calculation on a square lattice, after the translation step. Shown are 6 fluid sites and 3 wall sites (the wall is shown with a brick pattern).  The rest particle is denoted by f0. Click on the thumbnail to view a larger image.

Through careful choice of the equilibrium distribution, LB momentum and concentration obey the Navier-Stokes and advection-diffusion equations, and are Galilean-invariant. The equilibrium distribution is a simple function of the particle speed  and density  at each site, e.g. for the Qian et al. (1992) D2Q9 model:

Where Q_i = 1/9 for i = 1,3,5 and 7, and Q_i = 1/36 for i = 2,4,6 and 8.

Other distributions are given by Noble et al. (1996). The determines the viscosity and diffusion coefficients of the system; in the Qian et al. model, for example, the kinematic viscosity  = ( - 0.5) / 3.

There are several methods for including dispersion in an LB calculation (Martys and Chen, 1996; Flekkøy, 1993).  The Flekkøy method is perhaps the simplest; a carrier fluid following equations 3 and 4 produces a velocity field, and tracers (held in separate arrays) are advected by the carrier.  The concentrations of tracers at each (x,y,z) can be used to calculate a buoyancy force that feeds back into the  in equations 4a-4d (Martys and Chen, 1996).  Each tracer is characterized by its own _s, and the molecular diffusion coefficient of each tracer is given by D_m = (_s - 0.5) / 3 (for the square lattice in figure 7).

Because the LB method is floating point-based, it is inherently much less noisy than the LGA method; in general, there is no need to take momentum and particle averages to obtain meaningful velocities or concentrations at each (x,y,z).  Furthermore, the diffusion coefficients, which depend on , can be made much smaller than the LGA diffusion coefficients.  The consequences of these differences are apparent in this example of the Rayleigh-Taylor instability for miscible fluids (mpeg movie, 162 kB), which can be compared with figure 4, the LGA calculation for immiscible liquids.  Each frame in the LB calculation is taken from a single time step, with no averaging.  With only 128x128 sites, fine whorls of "smoke" are resolved, and the interfaces are less random-looking than in the LGA example.  With LGA, it is hard to get a diffusion coefficient D_m < 0.3, so it is hard to model miscible Rayleigh-Taylor or fingering displacements without making the simulation very large (Holme and Rothman, 1992).

Relative efficiency.   Figure 8 plots the MUPs, or millions (of site) updates per second, achievable on single-CPU computers against the year each CPU became widely available.  While LGA methods are still substantially faster on a per-site basis, the gap between the methods is narrowing. Although LGA use much less memory per site, they have a rather high ratio of memory to arithmetic operations, so the speed of LGA calculations is more bound by memory throughput, typically in the translation step. For LB on single-CPU systems, the speed of floating-point computations is generally the limitation; with our LB algorithms, the memory-intensive translation step rarely consumes more than 10 to 20% of the total computation time. As we will shortly show, there are other considerations that can make LB calculations effectively much faster than LGA. On parallel machines like the CM-5, the LB translation step must be handled carefully, else it will become the rate-limiting factor in performance (Noble et al., 1995, 1996).  The wide variation in CM-5 throughput results partly from the range of boundary conditions implemented with each algorithm. Modern single-CPU computers have begun to overlap the CM-5 range, but we expect the advances in single-CPU technology will soon be reflected in performance gains for parallel systems.

Figure 8.  Increase in performance of LB and LGA with time.  LGA algorithms are multi-site coded in x86 assembly language; LB algorithms are in C or HPF.  For the 21164a (DEC alpha), HWS refers to our general version of the code, and RJH refers to Robert J. Harley's specialized version. Click on thumbnail for larger image. The dashed red lines refer to 64-node (256 vector processor) CM-5 systems; upper value is estimated from Muders (1995) for single-phase, D2Q9 flow (Muders actually used a D3Q15 model), and lower is from Noble et al. (1995).

Coding considerations: LGA.  The two basic approaches to LGA coding have been referred to as "multi-spin" (in analogy with Ising models) and "multi-site" (Kohring, 1991; Brosa and Stauffer, 1989). Both approaches represent the absence or presence of a particle, with a particular velocity at a site (x,y,z), with a single bit; 0 for absent, 1 for present. The multi-spin method packs all bits corresponding to a given lattice direction into a single array. Thus for the FHP model (figure 1), there are 6 arrays corresponding to the 0°, 60°, 120°, 180°, 240° and 300° velocities, plus an array for the rest particle (if used), and an array to denote the presence of solids at the site. Within each multi-spin array, each machine word holds the particle information for 32 or 64 consecutive sites. The multi-spin method greatly speeds up the translation step, at the expense of making the collision step more complicated. While very high LGA update speeds can be achieved with the "multi-spin" coding technique, such models are generally restricted to single-component flow. In contrast, the multi-site method packs all the velocity bits for a site together in adjacent memory locations; with the FHP model, typically 6 velocity bits plus one solid bit (representing absence or presence of a solid wall) are packed into a byte, and 4 to 8 bytes packed per machine word. With multi-site, the translation step is more complicated, involving arithmetic ANDs and ORs to isolate and translate the individual velocity bits, operating on 4 to 8 sites in parallel. The multi-site collision step is simple, typically accomplished with look-up tables. Multiple components are most easily handled with separate arrays for "red" vs. "blue" vs. "green" components, though the size of the collision tables can increase rapidly as components are added.

We used the multi-site method, and hand-coded the translation, collision, immiscibility and averaging steps in 80x86 assembly language for speed.  To keep the look-up tables compact, we used a 6-velocity FHP model.  The diffusion algorithm was based on McNamara (1990), with a few added collisions to reduce the diffusion coefficient; extensive testing showed that the "molecular" diffusion coefficient did not change with flow speeds up to 0.2 lu/ts (lattice units per time step).  The immiscibilty algorithm was modified after Rothman and Keller (1988). On a 90 MHz Pentium®, this method achieves 4 million site updates per second for one-component flow, 2 million/second for two-component flow, and 1.x million/second with immiscibility. As a practical comparison, for the viscous fingering example (mpeg, 1916 kB), a 90 MHz Pentium® achieved ~30% of the site update speed reported by Lutsko et al. (1992) for a 400 node Transputer network. However, on single-processor 80x86 systems, the multi-site method is somewhat memory-bound, and is "running out of steam" with each advance in CPU technology.  The current generation of microprocessors (such as the Intel MMX®, Sun Sparclite®, and DEC 21164pc®) may obtain some speed increment through the use of "multi-media" instructions. Typically these instructions alias the floating point registers to represent 8 packed integer bytes; all bytes in a register can be acted on in parallel, not just by the traditional AND and OR instructions, but by integer multiplication, division and addition as well. It is doubtful, though, that such instructions would improve performance by more than a factor of two.

Coding considerations: LB.  We coded our LB algorithms in C, since C compilers are widely available, generally inexpensive, and are typically supplied with advanced development environments.  The simple dynamic allocation features of C/C++ are particularly useful, though dynamic allocation is now available in FORTRAN 90 as well.  In our 2D LB codes, we have declared arrays in two ways. The first is:

float ****F;

with the array elements F[s][i][y][x], indexed so that s selects the component of a multicomponent system, i selects the vector direction (i = 0 to 8 for the D2Q9 model), and y and x are the Cartesian spatial positions. The advantage of this order is that it allows much of the translation step to be accomplished with pointer swaps and efficient block moves. In the collision step, the F[s][i][y][x] for a given s, y and x are copied to a local array f[i] for floating-point manipulations; this copy step removes nearly all the aliasing problems of C (relative to FORTRAN).  The second method is:

typedef float site[NVELOC];     /* NVELOC = 9 for 2D */

site ***F;

which avoids one level of indirection, and is presumably more efficient, holding all vectors in adjacent memory locations.  However, this approach makes the translation stage more complicated, and requires buffering (that is, copying data to temporary buffers to void overwriting existing sites).  Our timings found the second method slightly slower than the first.

It is possible to achieve much higher speeds for LB algorithms, for somewhat specialized circumstances, by careful attention to processor architecture. Robert Harley of INRIA (denoted RJH on figure 8) modified our code to achieve 912000 site updates per second for 2D, two component flow on a 500 MHz DEC alpha 21164a processor, with the combination of Ziegler and no-flux boundary conditions described below. We estimate Harley's code could achieve ca. 1.8 to 2.0 MUPs for single component flow. This speed was obtained by copying all the F[][][][] variables for a given x,y, and for each component, to local register variables; by factoring algebraic expressions to reduce the number of additions and multiplications; and by inlining critical functions (e.g functions describing wall collisions). Some of this speed would be lost in moving to 3D, where it would be hard to hold all variables in registers, or if one used more complicated boundary collisions and more than two components.  However, this example shows that substantial speed gains can be obtained without resorting to assembly language coding or complicated logic.

The most striking contrast between LGA and LB coding is the simplicity of the latter. While the basic translation and collision steps of LGA are very simple, code for application of an unbiased body force, buoyancy, Hele-Shaw drag forces, and sorption can be extremely complicated in LGA. For example, the application of unbiased body force required approximately 100 lines of C code for LGA, and 1 line for LB.

Galilean Invariance and LGA.   A well-known problem with conventional LGA is the lack of Galilean invariance. The LGA bulk flow equation (1) contains the  term, which is exactly 1 in the Navier-Stokes equation.  For the FHP automaton,  = , so that .  When 1, the advection of vorticity does not obey the Galilean transforms.  For single component flow, Galilean invariance can be restored by adsorbing the  term in the velocity, pressure and force terms of equation (1), such that , etc.  However, this method decouples the Navier-Stokes equation from the advection-diffusion equation (2), such that solutes do not advect at the same rate as bulk flow.

For many dispersion problems, there is a better method to regain the Navier-Stokes equation.  Often we are interested in measuring dispersion of a passive tracer in a system with steady-state bulk flow; i.e for the case where . If 0 < < 3, the  term can be adsorbed in the viscosity, pressure and force, such that , etc.  This latter approach does not decouple the Navier-Stokes and advection-dispersion equations.

However, sometimes one must model dispersion in a system that is not at steady-state flow; for example, a system affected by buoyancy.  One recourse is to set =3, so that =0.  That is, we force the inertial term in equation (1) to zero, modeling Stokesian flow.  For many natural systems, Stokesian flow is quite reasonable; diffusion coefficients of solutes in water are typically 3 orders of magnitude lower than the viscosity of water, so Pe = 100 corresponds to Re = 0.1.  In addition, we often must compare against analytical solutions derived on the assumption of no inertial forces.

Accurate Boundary Conditions.  We may distinguish three types of boundary errors associated with dispersion calculations. The first applies to most numerical methods, and involves the choice of inlet conditions to represent, with a minimal number of sites or nodes, a process that properly occurs on a much larger scale. This first type of error was discussed by Stockman et al. (1997), and is also implicit in the discussion of viscous equilibrium below. The second concerns "pixelation" error, or the manner of representing smooth surfaces with the LB or LGA grid; a similar problem occurs in finite difference calculations, and elegant LB solutions have been suggested by Noble et al. (1995) and Ginzbourg and d'Humières (1996). The third, discussed in this section, concerns the manner of representing particle collisions at solid-fluid boundaries, and is mainly restricted to LB and LGA.

In the early days of LGA, the bounce-back collision was considered sufficiently accurate to describe the interface between a solid wall and a fluid.  The acceptance of bounce-back was easy with LGA, because the efficiency of LGA admitted use of wide channels, which diluted the error, and the noisiness of the method made it hard to discern errors in fluid velocity near a wall. With the advent of noiseless but computationally-intensive LB methods, and the corresponding need to decrease channel width, bounce-back was examined more carefully. Cornubert et al. (1991) showed for flow parallel to a wall, bounce-back placed the effective (zero-velocity) wall position 1/2 row into the fluid from the wall solid nodes. It is now known that if the wall position is simply "reinterpreted" to the 1/2 position, bulk-flow velocity fields calculated with bounce-back can be quite accurate, at least for slow flows (Chen et al., 1996). However, with the exception of the study by McNamara et al. (1995), little has been done with boundaries for tracer transport and dispersion, for which the lack of accuracy is potentially more serious.

To understand the boundary dilemma, we examine the LB calculation after the translation step.  At this stage, particle vectors on the solid-fluid interface point into the wall, but not out from the wall (figure 7). The bounce-back method is a simple solution to this problem; when the LB or LGA collision algorithm detects a wall node, the particle directions are reversed, sending the particles back into the fluid. This method seemed very natural for LGA, and has been used in many LB codes.

Simple application of the "reinterpreted" bounce-back to tracers causes an obvious problem; though the effective wall position for zero flow is 1/2 lu out, tracer can still get "behind" the 1/2 position, and reside on the real wall. There is a second shortcoming of the bounce-back method: it does not allow the fluid at the wall to inherit the diffusion and viscosity coefficients of the bulk fluid; in LB, bounce-back skips the relaxation process, with its characteristic . If one wishes to use the relaxation process, bounce-back has to be replaced with a method that specifies the "unknown" particle distributions that point into the fluid. Figure 9 shows two simple methods for bulk flow and dispersion. We emphasize that better solutions are available, especially for flow (Noble et al., 1995, Ginzbourg and d'Humières, 1996); we use these simple methods to illustrate the problem. The top part of the figure shows the Ziegler (1993) method; a similar approach was used by Martys and Chen (1996). For bulk flow, this method places a zero-velocity boundary on or near the wall. The Ziegler method is inappropriate for dispersion of Flekkøy-type tracers, since it smears any distribution of tracer tangential to the wall. At the bottom of the figure is a simple method to preserve concentration gradients along a "no-flux" wall; this second method is applied only to the tracers. With both, the contrived particle distribution is sent through the normal collision-relaxation process (equation 3).

Figure 9. Simple boundary conditions near a solid wall, after the translation step; initial known particle distributions are pink, and contrived particle distributions are orange and green. Top: Ziegler (1993) condition for fluid velocity (carrier fluid); the out-going particle distributions in directions 2,3 and 4 are inverses of the incoming distributions 6,7 and 8, respectively; 5=1= weighted average(6,7,8).  Bottom: no-flux condition for tracers; particles are allowed to stream along walls, and 2,3 and 4 are mirror reflections of 8,7 and 6, respectively.  Both subsequently go through the relaxation step.

If we place a slug of solute in a 5 lu-wide channel, with the true D_m = 0.05, and no flow (Pe=0), we can measure the apparent diffusion coefficient by the relaxation method. Under these conditions, the reinterpreted bounce-back gives an apparent diffusion coefficient ~98% of the true value, whereas the method shown in figure 9 gives back ~99.9999% of the true value. However, for Pe>0, the error can be substantially worse.

Figure 10 and figure 11 show how the choice of boundary method affects the accuracy of a simple dispersion calculation at Pe ~ 20. The basic geometry of the "experiment" is shown in figure 10: a slug of solute is injected into a flow between two parallel plates, and is stretched laterally by both advection and dispersion (similar to the previous Poiseiulle flow [mpeg] example). The distribution of solute, projected onto the x-axis, quickly becomes Gaussian, and is characterized by its standard deviation .

The effective diffusion coefficient D* (also called the dispersion coefficient) in the numerical "experiment" can be measured by the method of moments:

D* = 0.5 · dm₂ / dt                                                                             (5)

where m₂ is the second moment of the solute distribution -- that is, the variance, or ².

According to the Taylor-Aris theory (Aris, 1956), the real solution (after several diffusion times, described below) should be:

D*/D_m = 1 + Pe² / 210,                                                                     (6)

where Pe = W·U / D_m , W is the channel width, and U is the average flow speed through the channel.

Figure 11 plots the discrepancy between these two equations, i.e.:

[(D*/D_m)_measured - (D*/D_m)_Taylor] / (D*/D_m)_Taylor ,                        (7)

as a function of the channel width W. For this analysis, we kept flow speed U constant at either 0.05 or 0.01 l.u./timestep, and varied D_m to maintain a constant Pe = 20. We allowed the dispersion to complete 4 characteristic diffusion times (W ² / D_m); plots of m₂ vs. timestep were linear with correlation coefficient R ² = 0.999999 to 0.9999999999. At low W, the bounce-back method is actually more accurate, possibly because the sources of error cancel. However, the accuracy of the bounce-back method apparently improves first-order in W, while the no-flux method shows ~ second-order improvement in the range considered.

	Figure 10. Geometry for Taylor-Aris dispersion. Click on thumbnail for larger image.
	Figure 11. Fractional error (equation 7) in the dispersion coefficient for flow between parallel plates (figure 10), for the "reinterpreted" bounce-back method, and for the combination of no-flux and Ziegler methods (figure 9), as a function of the channel width W. Click on thumbnail for larger image.

One way to deal with this uncertainty is to make the channels wide, diluting the effect of the walls; here LGA seem to have an advantage, since the method is roughly 30 times faster than LB per node. However, for many dispersion problems, there are serious penalties for making the system larger. Unless the system is very simple, we will generally not know the equilibrium velocity distribution a priori; the normal approach is to start with a quiescent fluid, and apply a body force to accelerate the particles, until a steady-state flow distribution is reached. The number of calculation timesteps required for this "viscous equilibrium" is W ²/.  In a 2D calculation, the clock time per time step also goes as W ², so the net clock time to achieve viscous equilibrium goes as W ⁴.  For calculations that follow a slug of solute through a channel, the channel length L must accommodate several characteristic diffusion lengths U·(W ²/ D_m), so the total number of sites must be proportional to Pe·W ², ultimately resulting in a similar W ⁴ clock time constraint. Therefore simply making the channels 10 times wider may exact a severe penalty.

Ironically, the dispersion boundary problem is often less severe for flow through "rough", convoluted, or alveolated channels. For these more complicated geometries, much of the interface between the fluid and solid occurs in stagnant pockets on the side of the channel, where local Pe ~ 0, and the error in the measured diffusion coefficient is small.

Achievable Peclet numbers.  The LGA approach is constrained by the limited choice in diffusion coefficients. In general, we wish to match the Peclet number Pe for an experiment, typically defined as W·U / D_m, where W is some characteristic width (for example, the width of a channel or fracture), U is a fluid speed averaged across the channel, and D_m is the molecular diffusion coefficient. For LGA diffusion, we use a modification of the McNamara (1990) method, which gives a D_m of 0.277 lu²/ts at = 3.0 particle/site. Some variation in D_m can be obtained by varying , but for practical values of , the range is a factor of ~2. Thus, to vary Pe, one is really restricted to varying W or U.  For the 2-D hexagonal LGA, the practical U range is from about 0.001 to 0.1 lu/ts; lower values require inordinately long averaging times to overcome statistical noise, and higher U (higher Mach numbers) cause errors due to Galilean non-invariance (Hayot, 1987). With the LB method, D_m as low as 10^-4 can be achieved (Flekkøy, 1993) with 32-bit floating point, simply by varying the relaxation parameter _s  for the tracer species. However, in our experience, it is not always possible to vary independently the U, W, and D_m in LB systems, and maintain numerical stability. In the fracture-mixing problem to be described shortly, it was necessary to keep U ~0.02 lu/ts for W = 40, = 0.033 and D_m = 0.00033 lu²/ts. For the fracture-mixing geometry, higher fluid speeds produced negative concentrations that grew and eventually "infected" the entire simulation.  Flekkøy noticed similar, transient negative concentrations at the start of his Hele-Shaw simulations.

Examples

Background for examples 1 and 2 can be found in Stockman et al. (1997), Li (1995), and Perea-Reeves and Stockman (1997); here we focus on practical aspects of the calculations not covered in prior publications, and new LB results. Example 3 reflects ongoing research, and the laboratory experiments are described by Cooper et al. (1997).

Example 1: Mixing at Fracture Intersections.  The fracture mixing problem (figure 12) illustrates the limitations of both methods, though in this case, LB seems a clear winner.  Our intent is to determine the extent of solute mixing when two streams of fluid meet at a fracture junction; the answer is significant for modeling the advance of gaseous or dissolved pollutants through fractured rock.

Figure 12. Geometry for fracture mixing problem (Click on thumbnail for larger image). Fluid containing "black" tracer enters via channel 1; fluid with "white" tracer enters via channel 2. The two streams partially mix by diffusion, so some "gray" fluid exits via channels 3 and 4. The analysis domain is (2·L + W) · (2·L + W).

Typically, fractures in rocks are roughly 0.0001 to 0.1 cm wide, but the junctions are 10¹ to 10³ cm apart. Large-scale transport codes usually treat the fractured system as a network, applying an empirical mixing rule at the junctions, but they do not explicitly model the flow and solute distribution within each fracture. There are two end member assumptions about the mixing rule: (A) complete mixing, appropriate for low Pe; and (B) streamline routing, appropriate for high Pe (for this problem, the Pe is traditionally defined as W·U / (D_m·) ). The process is described by the mixing ratio:

M_r = J₃ / (J₃ + J₄)  ~  C₃ / (C₃ + C₄)                                        (8)

where the solute fluxes (J) and concentrations (C) are measured in the outlet legs.

For the "even flow" case, where the average flow velocities U into channels 1 and 2 are equal, assumption (A) predicts half the solute goes down each of legs 3 and 4 (M_r = 0.5), and assumption (B) predicts it all exits through leg 4 (M_r = 0). The problem is how one deals with "intermediate" Pe, where neither method is accurate, or how one determines the Pe where the transition from one method to the next is appropriate. Here we will consider modeling the problem with both LGA and LB, emphasizing calculations with "high" Pe.

Figure 13 compares the LB and LGA calculations for M_r, as a function of Pe, with experiments performed by Chunhong Li (1995). The agreement is reasonably good, especially considering that the location of the M_r = 0.25 cross-over point was previously uncertain to about two orders of magnitude in Pe. The scatter in the experimental data is notable, and indicates that LB or LGA calculations need not be extremely accurate to provide useful models for this physical process. Experiments of this nature, which require careful control of pump rates at very low fluid velocities, and non-intrusive measurement of solute concentration in minute quantities of fluid, are extremely difficult; the scatter in Li's data is remarkably small given these difficulties.

It is worthwhile to compare the calculation times required by the fracture mixing problem, for LB and LGA. The highest Pe number reachable with LGA was 221; the calculation required channels 572 lu wide, and took ~70 hours on a 200 MHz Pentium® Pro computer. The LB equivalent used channels 40 lu wide and took 1-2 hours. The wide channels were necessary for LGA, because the D_m was fixed at 0.277; consequently, the viscous equilibrium time became the main determinant of run time. Perhaps the only advantages of LGA for this problem are: the LGA method is well-understood, stable, and can be tested for exact conservation of particles and momentum; and it is relatively easy to force Stokesian flow (zero inertial term) by setting =3.0, so the results can be compared with other Stokesian calculations.

Figure 13. The solute mixing ratio (equation 8) for the geometry in figure 12 (Click on the thumbnail for a larger image). The experiments performed by Li measured the mixing of NaCl and KCl solutions in a plexiglas "fracture junction".

Example 2: Dispersion in Alveolated Channels.  An example that shows LGA in a reasonably good light, is the modeling of gas dispersion in "alveolated" or pocketed channels.  Dispersion experiments, meant to approximate behavior in the human lung, were performed by Tsuda et al. (1991); the basic geometry is shown in figure 14 and figure 15.  For a lung to be effective, the Pe for gaseous dispersion must be modest; this is equivalent to saying that there is an optimum breathing rate.  Tsuda et al. emphasized 0.1 < Pe < 100, and this range is achievable by LGA. The experiments used a quasi-2D geometry, with two parallel aluminum plates faced with a hexagonal honeycomb pattern to represent the alveoli; tracer gas containing several mole percent helium or sulfur hexafluoride was introduced at one end.  The effective diffusion coefficient (equation 5) describes the lateral (x-direction) spreading of a slug of tracer gas as it is pushed between the pocketed plates. Because we wish to keep the fins of the honeycomb thin to match the experiment, and because D*/D_m is typically proportional to W² at higher Pe, any uncertainty in the location of the boundaries is undesirable; thus we used channels ~100 to 200 lu wide to keep the relative error to a few %. In addition, the numerical experiment must be followed for several characteristic distances, x_c = U·W²/D_m; at the higher Pe numbers, this requirement translates into channels up to 200 lu by 10⁴ lu. The LGA calculation required 6 MB of memory, including arrays to hold running averages of particle velocity and concentration. The LB equivalents used narrower channels, but took approximately the same amount of calculation time.

	Figure 14. Dispersion of a slug of solute in flow through an alveolated channel (this is a small portion of the geometry; at high Pe, channels lengths were 100 times the channel width). The channel consists of two horizontal, planar plates, with small vertical septa protruding into the flow. This is a snapshot of the mpeg below. Click on thumbnail for larger image.
mpeg movie (146 kB) of above	Figure 15. A movie of dispersion in a small portion of an alveolated channel. The colors scale to the maximum solute concentration in the view, so as the slug leaves the right side, the solute lagging in the side pockets "lights up" and becomes more visible. This is an LB calculation, and is less "noisy" than the similar LGA calculation in figure 6.

Figure 16 compares the LGA and LB calculations with experimental results and the "stagnant pocket theory" derived by Tsuda et al. (1991) :

D*/D_m = 1/(1 + ß) + (Pe/(1 + ß))² · ((35ß² + 16ß + 2),                    (9)

where ß is the ratio of pocket volume to central channel volume. In figure 16, we use ß=0.73 to 0.75; this choice is motivated by the need to keep the central channel relatively wide, to minimize the uncertainty in location of the solid boundaries, yet keep the total number of sites relatively small for a comparison with LB. For Pe > 1, the pockets enhance the dispersion coefficient by up to a factor of 6 compared to simple Taylor-Aris dispersion in a channel with planar walls of similar separation. The cause of the increase is illustrated in the mpeg of the dispersion process; some solute gets trapped in the pockets, where the x-direction flow velocity is small, and thus lags behind the main solute pulse and increases the variance of the solute distribution (m₂). However, for very low Pe < 1, the pockets diminish the dispersion coefficient by up to a factor of 2, because when the flow is very slow, the walls of the alveoli inhibit lateral spreading. The hardest problems to model are those at high (>20) and low (<1) Pe. At high Pe, the number of sites has to be large to keep the solute slug from running off the end; at low Pe, we must worry about solute diffusing off the ends, especially since we are looking at a subtle effect (D*/D_m tends to be only slightly different from 1); with the LGA method, it is very hard to obtain accurate measurements of the very low required particle speeds. Give these difficulties, the agreement between LB and LGA is quite good. The slight disagreement with experiments at higher Pe partly reflects the difference in the way we and Tsuda defined the channel velocity used for the calculation of Pe. We used the average channel speed, in the spirit of Aris' original derivation, while Tsuda et al. used the experimentally measurable bulk flow speed. We transformed the latter to our definition, but the transformation required assumptions about the velocity field. Once again, we note the uncertainties in controlling and monitoring diffusion experiments can dwarf the errors in LB or LGA calculations.

Figure 16. Comparison of LGA and LB calculations with experimental data from Tsuda et al. (1991), the stagnant pocket equation (9) derived by those authors, and the Taylor-Aris equation (6) in a channel with the same width. Click on thumbnail for a larger image.

The real strength of LGA and LB lies in the ability to model more subtle experimental artifacts, such as the effect of buoyancy (figure 6, mpeg, 159 kB), which are unapproachable with conventional numerical methods. Reejhsinghani et al. (1966) showed that remarkably small density differences, between the solute slug and the carrier fluid, could large errors in the measured dispersion coefficient, even for horizontal tubes. Since the tracer gases used by Tsuda et al. were several % denser or lighter than the carrier gas, it is reasonable to ask if the results would be affected by buoyancy. From our LGA simulations, we estimated that the maximum effect for the geometry used by Tsuda et al. would be < 6% for a thin slug of solute, and < 25% for a step change source, which is not large compared to other experimental uncertainties. It seems that with the alveolated channel, the rough surfaces actually ameliorate the effects of buoyancy contrast.

Example 3: Double-Diffusive Fingering. Double-diffusive fingering is responsible for a variety of natural phenomena, and affects the transport of materials and heat in metals casting (Chen and Lu, 1992) and the evolution of stars (Grossman and Taam, 1996), and the transport of solutes around salt domes, saline basins, sea-floor vents and hydrothermal systems, and contaminated soils (Sarkar and Phillips, 1992). Perhaps the best known example of the phenomenon is "salt-fingering", which occurs when a layer of warm salt water overlies colder, fresher water as in the Mediterranean outflow (Holyer, 1984).  The juxtaposition is initially gravitationally stable; however, heat diffuses approximately 100 times faster than salt, and as the fresh water warms near the interface, it becomes more buoyant and begins to "finger" upward through the salty water. Modeling this phenomenon in the lab is difficult, since it is hard to make a "no-flux" heat boundary. Consequently, researchers often use layered salt and sugar solutions instead, with the salt solution on the bottom; the salt diffusion coefficient is about three times higher than the sucrose value. Often the solutions are placed in a thin Hele-Shaw cell, so the progress of the experiment is easily observed and the composition of the solutions can be monitored by absorption spectroscopy; in addition, the Hele-Shaw apparatus approximates the behavior of solutes in fractures or porous media, and the cell can be tipped about a horizontal axis to vary the effective permeability.

Figure 17 shows the progress of a laboratory salt-sugar experiment with false color representing the concentration of NaCl; each frame is approximately 2 cm wide. This is a very small window into the cell, which is 31 cm high, 17 cm wide, and 0.05 cm thick. The "red" salt solution (initially 0.18696 molal, 1.00618 g/cm³ at 20°C) is placed in the lower half of the cell, then a thin plexiglas divider strip is placed in the slot above the solution (the slot is visible slightly above the midpoint of each frame). Then the sugar solution (initially 0.0404 molal, 1.003548 g/cm³ at 20°C) is poured above the strip and allowed to settle. When the strip is withdrawn, there is an unavoidable oscillation of the surface; it is expected that the magnitude of the initial oscillation affects the onset, rate, and type of fingering, but there is some difficulty in characterizing the magnitude of the disturbance for modeling purposes. The wake of the strip is probably turbulent, so the solutes become somewhat mixed; furthermore, since the peaks and troughs of the oscillation are gravitationally unstable, they initially flatten rather rapidly, which may cause a Taylor-Aris spreading of the interface. While we can't match these initial conditions exactly with the LB method, we can attempt to resolve the importance of each factor. Figure 17 begins after the strip has been withdrawn and the initial, rapid mixing has occurred.

	Figure 17. Double-diffusive fingering experiment (Click on thumbnail for larger image). Each frame is ca. 2 cm wide, representing a small window into the larger Hele-Shaw cell. Red represents the NaCl solution, and blue the sucrose (sugar) solution; initially the NaCl solution is denser. Total time between 1st and last frame is ca. 56 minutes.
	Figure 18. Double-diffusive fingering lattice Boltzmann calculation (Click on thumbnail for larger image). Each frame is ca. 2 cm wide, representing a small window into the larger Hele-Shaw cell. Colors same as above.
Movie (mpeg, 157 kB)	Figure 19. LB calculation depicted in figure 18, as mpeg movie.

The experiment is characterized by solutal Rayleigh numbers; for the jth component:

Ra_j = B_j · C_j · (g·sinØ · (h²/12)·L / (D_mj ·)                               (10)

where C_j is the concentration, B_j = d / dC_j is the molar expansivity of the jth component, g is the gravitational acceleration, Ø is the angle of the Hele-Shaw cell with respect to horizontal, h is the cell plate separation, and L is a characteristic length scale for the initial "disturbance" of the interface (theoretically 0 for a perfectly sharp and smooth interface). It may be desirable to match the Prandtl number Pr =/ D_m , though since the flow is generally slow and Stokesian, Pr is largely ignored in much of the literature for double-diffusive convection.

Our LB model for this process is modified after methods developed by Flekkøy's (1993) and Flekkøy et al. (1995), which in turn use the Hele-Shaw drag term tested by Holme and Rothman (1992). The LB method is 2D, but emulates 3D by specifying a fictitious cell thickness h, which is perpendicular to the page in figure 18. The drag force is then:

 .                                                          (11)

It is somewhat difficult to match the Ra_j for the experiments with LB. We could make h and L large, but that implies using a large number of sites, thus long run times. To obtain high Ra_j without using very large numbers of sites, we can use either large g·sinØ or small D_mj ·. We can't make the diffusion coefficients arbitrarily small (especially in a 32-bit model), so adjusting the "gravity" is the principal option. If the viscosity is made small, we risk not matching the Pr (though the significance of this match is not clear at low Reynolds number Re). Though we wish to model essentially incompressible fluids, the LB fluids are numerically compressible, so very high gravity will cause artifacts by inducing non-physical density gradients along the model cell. Ultimately, it is barely possible to match Ra_j = 5000 for a 256x64 site automaton, if the carrier fluid is water and the solutes have D_m ca. 10^-5 cm²/sec.

The effect of the LB "compressibility" is lessened somewhat by the manner of calculating the buoyancy force. The density at a given (x,y) is calculated via:

                          (12)

where the concentration differences are calculated relative to a reference value, corresponding to ₀. If we choose the reference carefully, and calculate buoyancy proportional to  -₀ , the total force across the automaton is lessened, even though the relative buoyancy is maintained.

Ironically, it would really be much easier to model salt-heat diffusion with LB than salt-sugar diffusion. When the diffusion coefficient is much larger for one component than the other, it is most important to match the D_m / for the faster diffusing component. For the salt-sugar experiment, the D_m / ~ 10^-3, which is hard to achieve in a 32-bit LB; for true salt fingering, D_m / ~ 3·10^-2.

While this quasi-3D LB is adequate for most Hele-Shaw salt-sugar experiments, it fails when the fingers begin to move rapidly. In a Hele-Shaw cell, the Pe defined relative to the cell thickness is  h ·U / D_m . From the Taylor-Aris equation, we see that for Pe > (210)^{0. 5} ~ 14.5, the effective diffusion coefficient D* will be more than twice the D_m ; for our experiments, this condition occurs when the mean finger speed exceeds ~ 3·10^-3 cm/csec. From figure 17, the average finger speed appears to be ~1.4·10^-3, perilously close to the point where advective dispersion, from the narrow plate gap, begins to dominate plain molecular diffusion.

We have experimented with approximate methods to model the effect of Taylor-Aris dispersion in our quasi-3D LB, by using an effective D*_j  that is a function of fluid velocity. We note that the D*_j in a moving finger is anisotropic; a method for calculating such an anisotropic D*_j has been suggested by Maillot and Main (1996). We have tested this method, but it has the constraint that D*_j /D_mj must be < 6 for the D2Q9 lattice, and so far it appears to introduce numerical instability. Eventually it will be necessary to use a true 3D LB to model higher-speed finger motions.

Though Taylor-Aris dispersion can affect the spreading of fingers, it probably has a much more profound effect on the initial condition of the experiments and LB calculation. It is apparent from figure 19 (mpeg) that the initial step disturbance of the interface relaxes very quickly, since it is gravitationally unstable; during this relaxation stage, the fluid speeds are actually much greater than any achieved through fingering itself. Consequently, Taylor-Aris dispersion is unavoidable at the early stage, and greatly smears out the interface. Thus the initial, sharp interface used in our LB calculations might appropriately be replaced with a smeared concentration distribution.

Could this experiment have been modeled with LGA?  The answer depends on what we believe to be the most important dimensionless parameters, and how willing we are to make significant changes to the algorithms.  With conventional LGA, it is hard to make the viscosity very different from the diffusion coefficient, so matching the experimental D_m / is unlikely; however, if we force Stokesian flow by setting =3, this point may be moot.  With the McNamara-type diffusion coefficient, and some selectivity in the rules, it is reasonable to design a system where D_m1 / D_m2 = 3; it is not easy to design a system where D_m1 / D_m2 = 100, so modeling real salt fingering is probably out of the question with LGA. The minimum LGA viscosity will be on the order 0.1 lu/ts, and the diffusion coefficients will be ~0.25 and 0.75, so the denominator in the solutal Rayleigh number will be ~10³ times larger than we can achieve with the LB model. The easiest way to compensate for this difference would be to make the h²·L term in the numerator 10³ times larger; ultimately, this means that we must increase the x and y LGA spatial dimensions by a factor of ten. The time scale, 12·L·/ (h²·g), would be approximately the same for both LB and LGA methods, since the LGA will accept about a ten times smaller g per site per time step.  The LGA would require about 10² times as many sites, but the LGA codes run about 30 times as fast as the LB codes per site, so in summary the LGA calculation might take ~3 times as long. However, the real limitation of the LGA method would be the lack of tunabilty; the D_m1 / D_m2 for the model would be a function of  alone, and would be fixed for the Stokesian case (=3 in the 6-velocity model). Developing the 3-color collision tables and debugging the Hele-Shaw perturbation model could be especially tedious, and there would be no simple way to add a Taylor-Aris function to dynamically change the diffusion coefficient as a function of fluid speed.

Conclusions

We have shown LB and LGA can be used for quantitative dispersion models, and facilitate the analysis of laboratory experiments. The uncertainties in laboratory dispersion studies are often larger than the uncertainties inherent in modeling the systems with discrete velocity methods. While a great deal of effort has gone into developing accurate LB boundary conditions for bulk flow, boundaries have a more serious effect on the accuracy of LB dispersion calculations. However, the variability in initial and boundary conditions of lab experiments (illustrated with the double-diffusive fingering problem) may swamp the errors in LB and LGA boundaries.

It is apparent that LB methods will supplant LGA for experiments that require tunable transport coefficients, particularly when the system is characterized by multiple diffusion coefficients that vary by orders of magnitude. The principal advantages of LGA are the low memory use and high speed per site; when it is necessary to match the exact geometry of an experiment (by using many sites to define solid-fluid boundaries), or it is necessary to minimize boundary errors by "dilution", these are legitimate advantages. However, the rapid increase in the speed of floating point hardware, the drop in computer memory prices, and the advent of accurate boundary methods will make LB methods even more attractive in the near future.

Acknowledgments

This work was supported by the U.S. Dept. of Energy, Office of Basic Energy Sciences, Geosciences Research Program under contract # DE-AC04-94AL85000, and by the Sandia LDRD Office. Additional support was provided by the National Science Foundation, Grant No. EAR-9405837.

We thank Joanne Fredrich, Ruaidhri O'Connor, John Wilson and Bob Glass for support and helpful discussions, Eirik Flekkøy and Dan Rothman for advice on anisotropic Hele-Shaw models, Robert J. Harley for his optimization strategies, and David Noble for advice on boundary conditions. Akira Tsuda generously made his dispersion data available for analysis, and provided us with details of the experimental set-up.

References

Aris, R. (1956)  On the dispersion of a solute in a fluid flowing through a tube, Proc. Royal Soc. 235A, 67-77.

Baudet, C.; Hulin, J.P.; Lallemand, P. and d’Humières, D. (1989)  Lattice-gas automata: a model for the simulation of dispersion phenomena, Phys. Fluids A 1, 507-512.

Brosa, U. and Stauffer, D. (1989) Vectorized multisite coding for hydrodynamic cellular automata, Jour. Stat. Phys. 57, 399-403.

Brown, S.R.; Stockman, H.W. and Reeves, S.J. (1995)  Applicability of the Reynolds equation for modeling fluid flow between rough surfaces, Geophys. Res. Lett. 22, 2537-2540.

Chen, F. and Lu, J.W. (1992) Onset of salt-finger convection in anisotropic and inhomogeneous porous media, Int. Jour. Heat Mass Transfer 35, 3451-3464.

Chen, S.J.; Martínez, D. and Mei, R. (1996)  On boundary conditions in lattice Boltzmann methods, Phys. Fluids 8, 2527-2535.

Cooper, C.; Glass, R.J. and Tyler, S. (1997)  Experimental investigation of the stability boundary for double-diffusive finger convection in a Hele-Shaw cell, Water Resources Res. 33, 517-526.

Cornubert, R.; d'Humières, D. and Levermore, D. (1991)  A Knudsen layer theory for lattice gases, Physica D 47, 241-259.

Flekkøy, E.G. (1993)  Lattice Bhatnagar-Gross-Krook models for miscible fluids, Phys. Rev. E 47(6), 4247-4257.

Flekkøy, E.G.; Oxaal, U.; Feder, J. and Jøssang, T. (1995)  Hydrodynamic dispersion at stagnation points: simulations and experiments, Phys. Rev. E 52, 4952-4962.

Frisch, U.; Hasslacher, B. and Pomeau, Y. (1986)  Lattice-gas automata for the Navier-Stokes equation, Phys. Rev. Lett. 56, 1505-1508.

Ginzbourg, L. and d'Humières, D. (1996)  Local second-order boundary methods for lattice Boltzmann models, Jour. Stat. Phys. 84, 927-971.

Grossman S.A. and Taam, R.E. (1996)  Double-diffusive mixing-length theory, semiconvection and massive star evolution, Monthly Notices Royal Astro. Soc. 283(#4) 1165-1178.

Gutfraind, R.; Ippolito, I. and Hansen, A. (1995)  Study of tracer dispersion in self-affine fractures using lattice-gas automata, Phys. Fluids 7, 1938-1948.

Hayot, F. (1987) The effect of Galilean non-invariance in lattice gas automaton one-dimensional flow, Complex Systems 1, 753-761.

Holme, R. and Rothman, D.H. (1992)  Lattice-gas and lattice-Boltzmann models of miscible fluids, Jour. Stat. Phys. 68, 409-429.

Holyer, J.Y. (1984) The stability of long, steady, two-dimensional salt fingers, Jour. Fluid Mech. 147, 169-185.

Kohring, G.A. (1991) Parallelization of short- and long-range cellular automata on scalar, vector, SIMD and MIMD machines, Inter. Jour. Modern Phys. 2, 755-772.

Li, C. (1995)  Low Peclet Number Mixing Behavior at Fracture Junctions, Ph.D. thesis, 172 pp., New Mexico Institute of Mining and Technology.

Lutsko, J.F.; Boon, J.P. and Somers, J.A. (1992)  Lattice gas automata simulations of viscous fingering in a porous medium, in: T.M.M. Verheggen (ed.) Numerical Methods for the Simulation of Multi-Phase and Complex Flow, Springer-Verlag, Berlin, 124-135.

Maillot, B. and Main, I.G. (1996)  A lattice BGK model for the diffusion of pore fluid pressure, including anisotropy, heterogeneity, and gravity effects, Geophys. Res. Lett. 23, 13-16.

Martys, N. and Chen, H. (1996)  Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method, Phys. Rev. E  53, 743-750.

McNamara, G.R. (1990)  Diffusion in a lattice gas automaton, Europhys. Lett. 12, 329-334.

McNamara, G.R.; Garcia, A.L. and Alder, B.J. (1995)  Stabilization of thermal lattice Boltzmann methods, Jour. Stat. Phys. 81, 395-408.

Muders, D. (1995)  Three-Dimensional Parallel Lattice Boltzmann Hydrodynamic Simulations of Turbulent Flows in Interstellar Dark Clouds, Ph.D. thesis, 81 pp., Universität Bonn. (available on www)

Noble, D.R.; Georgiadis, J.G. and Buckius, R.O. (1995)  Direct assessment of lattice Boltzmann hydrodynamics and boundary conditions for recirculating flows, Jour. Stat. Phys. 81, 17-33.

Noble, D.R.; Georgiadis, J.G. and Buckius, R.O. (1996)  Comparison of accuracy for lattice Boltzmann and finite difference simulations of steady viscous flow, Inter. Jour. Numerical Meth. Fluids 23, 1-18.

Perea-Reeves, S.J. and Stockman H.W. (1997)  A lattice-gas study of dispersion in alveolated channels, in press Chem. Eng. Sci.

Qian, Y.H.; d’Humières, D. and Lallemand, P. (1992)  Lattice BGK models for Navier-Stokes equation, Europhys. Lett. 17(6 BIS), 479-484.

Reejhsinghani, N.S.; Gill, W.N. and Barduhn, A.J. (1966)  Part III. Experiments in horizontal tubes including observations on natural convection effects, A.I.Ch.E. Jour. 12, 916-921.

Rothman, D.H. (1988)  Cellular-automaton fluids: a model for flow in porous media, Geophysics 53, 509-518.

Rothman, D.H. and Keller, J.M. (1988)  Immiscible cellular-automaton fluids, Jour. Stat. Phys. 52, 1119-1127.

Sarkar, A. and Phillips, O.M. (1992) Effects of horizontal gradients on thermohaline instabilities in a thick porous layer, Phys. Fluids A 4, 1165-1175.

Stockman, H.W.; Li, C. and Wilson, J.L. (1997)  A lattice-gas and lattice Boltzmann study of mixing at continuous fracture junctions: importance of boundary conditions, in press Geophys. Res. Lett.

Tsuda, A.; Federspiel, W.J.; Grant, P.A. Jr. and Fredberg, J.J. (1991)  Axial dispersion of inert species in alveolated channels, Chem. Eng. Sci. 46, 1419-1426.

Ziegler, D.P. (1993)  Boundary conditions for lattice Boltzmann simulations. Jour. Stat. Phys. 71, 1171-1177.

Modern Unix systems typically come configured with mpeg viewers. The SGI (ca. 1996) mpeg viewers will handle all the movies in the article, though they are sluggish on 100 MHz R4000 systems.

Microsoft Windows® and Apple Macintosh® systems are typically not supplied with mpeg viewers. Here are some popular viewers, and associated problems and strengths relevant to this article (comments by H.W. Stockman). At the time of this writing (April 1997), these programs were widely available for download via the www. After downloading and installing the viewer, it is usually necessary to register the program as a helper application in the browser (e.g., via the : Options | General Preferences | Helpers dialogue in Netscape®.)

These comments do not constitute endorsements of any programs.

For MS Windows ®:

vmpegwin.exe: The "lite" version 1.7 is freeware, and can handle the largest mpegs in this report (2.1 MB). The viewer has good single-step, loop and size control, and is very fast. However, the way the program handles successive invocations is unusual; it will not automatically play a second mpeg after it is invoked by the browser-- it must be shut down between invocations. The shutdown can be handled automatically if the /close argument is added to the command line (when the program is registered in the browser as a "helper" application).

mpegplay.exe: The shareware version will not play the largest files in this document. The viewer ocassionally trips up on some 24-bit color files, but will handle successive invocations easily. Allows loop and single-step.

InterVU: Freeware versions are available for both Mac and Windows (PC). The viewer is highly integrated with the browser and will play very large files, but the freeware version has no easy loop or single-step capability.

NetTOOB: The tested version did not correctly play many of the larger mpegs, and oddly distorted others.

Microsoft mplayer (with active movie): The tested version tended to make colors blocky and would not play many of the larger files.

For the Macintosh ®:

I don't use Macintosh computers very often, so I can't offer authoritative information. Maynard Handley's Sparkle plays mpegs, is widely available on the net, and seems very capable for turning mpegs into QuickTime. InterVU (described above) has a free version for the Macintosh, and MacZilla is another popular Mac mpeg/avi/Quicktime viewer.