Annu. Rev. Fluid Mech. 1998. 30:329–64 Copyright c 1998 by Annual Reviews Inc. All rights reserved
LATTICE BOLTZMANN METHOD FOR FLUID FLOWS
Shiyi Chen1,2 and Gary D. Doolen2
r />Research Division, T. J. Watson Research Center, P.O. Box 218, Yorktown Heights, NY 10598; 2Theoretical Division and Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, NM 87545; e-mail: email@example.com
KEY WORDS: lattice Boltzmann method, mesoscopic approach, ?uid ?ow simulation
We present an overview of the lattice Boltzmann method (LBM), a parallel and ef?cient algorithm for simulating single-phase and multiphase ?uid ?ows and for incorporating additional physical complexities. The LBM is especially useful for modeling complicated boundary conditions and multiphase interfaces. Recent extensions of this method are described, including simulations of ?uid turbulence, suspension ?ows, and reaction diffusion systems.
In recent years, the lattice Boltzmann method (LBM) has developed into an alternative and promising numerical scheme for simulating ?uid ?ows and modeling physics in ?uids. The scheme is particularly successful in ?uid ?ow applications involving interfacial dynamics and complex boundaries. Unlike conventional numerical schemes based on discretizations of macroscopic continuum equations, the lattice Boltzmann method is based on microscopic models and mesoscopic kinetic equations. The fundamental idea of the LBM is to construct simpli?ed kinetic models that incorporate the essential physics of microscopic or mesoscopic processes so that the macroscopic averaged properties obey the desired macroscopic equations. The basic premise for using these simpli?ed kinetic-type methods for macroscopic ?uid ?ows is that the macroscopic dynamics of a ?uid is the result of the collective behavior of many microscopic particles in the system and that the macroscopic dynamics is not sensitive to the underlying details in microscopic physics (Kadanoff 1986). By developing a simpli?ed version of the kinetic equation, one avoids solving 329 0066-4189/98/0115-0329$08.00
CHEN & DOOLEN
complicated kinetic equations such as the full Boltzmann equation, and one avoids following each particle as in molecular dynamics simulations. Even though the LBM is based on a particle picture, its principal focus is the averaged macroscopic behavior. The kinetic equation provides many of the advantages of molecular dynamics, including clear physical pictures, easy implementation of boundary conditions, and fully parallel algorithms. Because of the availability of very fast and massively parallel machines, there is a current trend to use codes that can exploit the intrinsic features of parallelism. The LBM ful?lls these requirements in a straightforward manner. The kinetic nature of the LBM introduces three important features that distinguish it from other numerical methods. First, the convection operator (or streaming process) of the LBM in phase space (or velocity space) is linear. This feature is borrowed from kinetic theory and contrasts with the nonlinear convection terms in other approaches that use a macroscopic representation. Simple convection combined with a relaxation process (or collision operator) allows the recovery of the nonlinear macroscopic advection through multi-scale expansions. Second, the incompressible Navier-Stokes (NS) equations can be obtained in the nearly incompressible limit of the LBM. The pressure of the LBM is calculated using an equation of state. In contrast, in the direct numerical simulation of the incompressible NS equations, the pressure satis?es a Poisson equation with velocity strains acting as sources. Solving this equation for the pressure often produces numerical dif?culties requiring special treatment, such as iteration or relaxation. Third, the LBM utilizes a minimal set of velocities in phase space. In the traditional kinetic theory with the MaxwellBoltzmann equilibrium distribution, the phase space is a complete functional space. The averaging process involves information from the whole velocity phase space. Because only one or two speeds and a few moving directions are used in LBM, the transformation relating the microscopic distribution function and macroscopic quantities is greatly simpli?ed and consists of simple arithmetic calculations. The LBM originated from lattice gas (LG) automata, a discrete particle kinetics utilizing a discrete lattice and discrete time. The LBM can also be viewed as a special ?nite difference scheme for the kinetic equation of the discrete-velocity distribution function. The idea of using the simpli?ed kinetic equation with a single-particle speed to simulate ?uid ?ows was employed by Broadwell (Broadwell 1964) for studying shock structures. In fact, one can view the Broadwell model as a simple one-dimensional lattice Boltzmann equation. Multispeed discrete particle velocities models have also been used for studying shock-wave structures (Inamuro & Sturtevant 1990). In all these models, although the particle velocity in the distribution function was discretized, space and time were continuous. The full discrete particle velocity model, where
LATTICE BOLTZMANN METHOD
space and time are also discretized on a square lattice, was proposed by Hardy et al (1976) for studying transport properties of ?uids. In the seminal work of the lattice gas automaton method for two-dimensional hydrodynamics, Frisch et al (1986) recognized the importance of the symmetry of the lattice for the recovery of the Navier-Stokes equation; for the ?rst time they obtained the correct Navier-Stokes equation starting from the lattice gas automata on a hexagonal lattice. The central ideas in the papers contemporary with the FHP paper include the cellular automaton model (Wolfram 1986) and the 3-D model using the four-dimensional face-centered-hyper-cubic (FCHC) lattice (d’Humi` res e et al 1986). The lattice gas automaton is constructed as a simpli?ed, ?ctitious molecular dynamic in which space, time, and particle velocities are all discrete. From this perspective, the lattice gas method is often called lattice gas cellular automata. In general, a lattice gas automaton consists of a regular lattice with particles residing on the nodes. A set of Boolean variables n i (x, t)(i = 1, · · · , M) describing the particle occupation is de?ned, where M is the number of directions of the particle velocities at each node. The evolution equation of the LG automata is as follows: n i (x + ei , t + 1) = n i (x, t) +
i (n(x, t)),
(i = 0, 1 · · ·, M), (1)
where ei are the local particle velocities. Starting from an initial state, the con?guration of particles at each time step evolves in two sequential sub-steps, (a) streaming, in which each particle moves to the nearest node in the direction of its velocity, and (b) collision, which occurs when particles arriving at a node interact and change their velocity directions according to scattering rules. For simplicity, the exclusion principle (no more than one particle being allowed at a given time and node with a given velocity) is imposed for memory ef?ciency and leads to a Fermi-Dirac local equilibrium distribution (Frisch et al 1987). The main feature of the LBM is to replace the particle occupation variables, n i (Boolean variables), in Equation 1 by single-particle distribution functions (real variables) f i = n i and neglect individual particle motion and particle-particle correlations in the kinetic equations (McNamara & Zanetti 1988), where denotes an ensemble average. This procedure eliminates statistical noise in the LBM. In the LBM, the primitive variables are the averaged particle distributions, which are mesoscopic variables. Because the kinetic form is still the same as the LG automata, the advantages of locality in the kinetic approach are retained. The locality is essential to parallelism. An important simpli?cation of the LBM was made by Higuera & Jim? nez e (1989) who linearized the collision operator by assuming that the distribution is close to the local equilibrium state. An enhanced collision operator approach which is linearly stable was proposed by Higuera et al (1989). A particular
CHEN & DOOLEN
simple linearized version of the collision operator makes use of a relaxation time towards the local equilibrium using a single time relaxation. The relaxation term is known as the Bhatnagar-Gross-Krook (BGK) collision operator (Bhatnagar et al 1954) and has been independently suggested by several authors (Qian 1990, Chen et al 1991). In this lattice BGK (LBGK) model, the local equilibrium distribution is chosen to recover the Navier-Stokes macroscopic equations (Qian et al 1992, Chen et al 1992). Use of the lattice BGK model makes the computations more ef?cient and allows ?exibility of the transport coef?cients.
LATTICE BOLTZMANN EQUATIONS LBE: An Extension of LG Automata
There are several ways to obtain the lattice Boltzmann equation (LBE) from either discrete velocity models or the Boltzmann kinetic equation. There are also several ways to derive the macroscopic Navier-Stokes equations from the LBE. Because the LBM is a derivative of the LG method, we will introduce the LBE beginning from a discrete kinetic equation for the particle distribution function, which is similar to the kinetic equation in the LG automata in Equation 1: f i (x + ei x, t + t) = f i (x, t) +
f (x, t)),
(i = 0, 1 · · ·, M),
where f i is the particle velocity distribution function along the ith direction; i = i ( f (x, t)) is the collision operator which represents the rate of change of f i resulting from collision. t and x are time and space increments, respectively. When x/ t = |ei |, Equations 1 and 2 have the same discretizations. i depends only on the local distribution function. In the LBM, space is discretized in a way that is consistent with the kinetic equation, i.e. the coordinates of the nearest neighbor points around x are x + ei . The density ρ and momentum density ρu are de?ned as particle velocity moments of the distribution function, f i , ρ=
f i ei ,
M where i ≡ i=1, . i is required to satisfy conservation of total mass and total momentum at each lattice: i i
If only the physics in the long-wave–length and low-frequency limit are of interest, the lattice spacing x and the time increment t in Equation 2 can be regarded as small parameters of the same order ε. Performing a Taylor
LATTICE BOLTZMANN METHOD
expansion in time and space, we obtain the following continuum form of the kinetic equation accurate to second order in ε: ? fi + ei · ?t fi + ε 1 ei ei : 2 f i + ei · 1 ? 2 fi ? fi + ?t 2 ?t 2 =
To derive the macroscopic hydrodynamic equation, we employ the ChapmanEnskog expansion, which is essentially a formal multiscaling expansion (Frisch et al 1987), ? ? ? =ε + ε2 , ?t ?t1 ?t2 ? ? =ε . ?x ? x1
The above formula assumes that the diffusion time scale t2 is much slower than the convection time scale t1 . Likewise, the one-particle distribution function f i can be expanded formally about the local equilibrium distribution function eq fi , fi = fi + ε fi
eq fi eq (neq)
Here depends on the local macroscopic variables (ρ and ρu) and should satisfy the following constraints: fi
i (neq) eq
f i ei = ρu.
fi = f i(1) + ε f i(2) + O(ε2 ) is the nonequilibrium distribution function, which has the following constraints: f i(k) = 0,
f i(k) ei = 0,
for both k = 1 and k = 2. Inserting f i into the collision operator
the Taylor expansion gives:
f) = +
f eq ) + ?
i ( f ) (1) fj ?fj eq
f eq )
f j(2) +
? 2 i ( f eq ) (1) (1) f j fk ? f j ? fk
+ O(ε 3 ).
From Equation 5, we note that when ε → 0, we have: to a linearized collision operator,
f eq ) = 0. This leads (10)
Mi j eq fj ? fj , ε
i where Mi j ≡ ? ?(f fj ) is the collision matrix (Higuera & Jim? nez 1989), which e determines the scattering rate between directions i and j. For a given lattice,
CHEN & DOOLEN
Mi j only depends on the angle between directions i and j and has a limited set of values. For mass and momentum conservation collision, Mi j satisfy the following constraints (Benzi et al 1992):
Mi j = 0,
ei Mi j = 0.
If we further assume that the local particle distribution relaxes to an equilibrium state at a single rate τ , 1 (12) Mi j = ? δi j , τ we arrive at the lattice BGK collision term (Bhatnagar et al 1954), 1 neq 1 i = ? fi = ? f (1) + ε f i(2) , (13) ε τ ετ i and the LBGK equation: fi ? fi . (14) τ The BGK collision term (Qian 1990, Chen et al 1991, Qian et al 1992, Chen et al 1992) was used previously in the full Boltzmann simulation for studying shock formation (Chu 1965) and was used recently for a shock-capturing using ?nite-volume methods (Xu & Prendergast). From Equation 5 one obtains the following equations: f i (x + ei , t + 1) = f i (x, t) ? ? fi + ei · ?t1
eq eq 1 fi eq
f i(1) , τ 1 f i(1) + ei ei : 2
to order ε0 and ? (1) ? eq f + f + ei · ?t1 i ?t2 i + ei ·
? eq 1 ? 2 eq 1 f + f = f i(2) , 2 ?t1 i 2 ?t1 i τ
to order ε 1 . Using Equation 15 and some algebra, we can rewrite the ?rst order equation as 2 ? f i(1) + 1? ?t2 τ ? f i(1) + ei · ?t1
(1) 1 fi
f i(2) . τ
From Equations 15 and 17 we obtain the following mass and momentum equations: ?ρ + · ρu = 0, (18) ?t
LATTICE BOLTZMANN METHOD
?ρu + ?t
which are accurate to second order in ε for Equation 2. Here the momentum ?ux tensor has the form:
(ei )α (ei )β f i + 1 ?
f i(1) ,
and (ei )α is the component of the velocity vector ei in the α-coordinate direction. To specify the detailed form of αβ , the lattice structure and the corresponding equilibrium distribution have to be speci?ed. For simplicity and without loss of generality, we consider here the two-dimensional square lattice with nine velocities: ei = (cos(π/2(i ? 1), sin(π/2(i ? 1)) for i = 1, 3, 5, 7, √ ei = 2(cos(π/2(i ? 1) + π/4, sin(π/2(i ? 1) + π/4) for i = 2, 4, 6, 8; e0 = 0 corresponds to a zero-speed velocity. The requirement for using the nine-velocity model, instead of the simpler ?ve-velocity square lattice, comes from the consideration of lattice symmetry: the LBE cannot recover the correct Navier-Stokes equations unless suf?cient lattice symmetry is guaranteed (Frisch et al 1986). Note that the Navier-Stokes equations has a second-order nonlinearity. The general form of the equilibrium distribution function can be written up to O(u 2 ) (Chen et al 1992): fi
= ρ a + bei · u + c(ei · u)2 + du 2 ,
where a, b, c and d are lattice constants. This expansion is valid only for small velocities, or small Mach numbers u/Cs , where Cs is the sound speed. Using the constraints in Equation 7, the coef?cients in Equation 21 can be obtained analytically (Qian et al 1992): 9 3 eq f i = ρwi 1 + 3ei · u + (ei · u)2 ? u 2 , 2 2 (22)
with w0 = 4/9, w1 = w3 = w5 = w7 = 1/9, and w2 = w4 = w6 = w8 = 1/36. Inserting the above formula into Equation 20 we have,
(ei )α (ei )β f i 1? 1 2τ
= pδαβ + ρu α u β , (23)
α (ρuβ )
(ei )α (ei )β f i1 = ν(
β (ρuα )),
where p = ρ/3 is the pressure, which gives a constant sound speed, Cs = √ 1/ 3, and ν = (2τ ? 1)/6 is the kinematic viscosity.
CHEN & DOOLEN
The resulting momentum equation is ?uα (24) + β · uα uβ = ? α p + ν β · ( α ρuβ + β ρuα ), ?t which is exactly the same as the Navier-Stokes equations if the density variation δρ is small enough (Qian & Orszag 1993). ρ
LBE: Approximation to the Continuum Boltzmann Equation
Although the LBE evolved from its Boolean counterpart, the LG automaton method, it has been shown recently by two groups independently (He & Luo 1997, Abe 1997) that the LBE can be obtained from the continuum Boltzmann equation for discrete velocities by using a small Mach number expansion. In these derivations, the starting point is the Boltzmann BGK equation (Bhatnagar et al 1954, Chu 1965, Xu & Prendergast 1994): 1 ?g + ξ · g = ? (g ? g eq ), (25) ?t ετ where g ≡ g(x, ξ, t) is the single-particle distribution function in continuum phase space (x, ξ), and g eq is the Maxwell-Boltzmann equilibrium distribution function: ρ 3 g eq ≡ exp ? (ξ ? u)2 . (26) (2π/3) D/2 2 Here D is the spatial dimension, and for simplicity, the particle velocity ξ and √ the ?uid √ velocity u have been normalized by 3RT , giving a sound speed of cs = 1/ 3. T is the temperature. The macroscopic ?uid variables are the velocity moments of the distribution function g: 1 (27) (ξ ? u)2 gdξ, ρ = gdξ, ρu = ξgdξ, ρε = 2 where ε = D T is the internal energy. Assuming that the ?uid velocity in Equa2 tion 26 is a small parameter (compared with the sound speed), the equilibrium distribution g eq up to O(u 2 ) has the following form (Koelman 1991): ρ 3 3 9 exp ? ξ 2 1 + 3(ξ · u) + (ξ · u)2 ? u 2 . (28) (2π/3) D/2 2 2 2 For discrete velocity models, only a small set of particle velocities ei = ξi (i = 1, . . . , M), and their distribution functions at these velocities, gi (x, t) = g (x, ei , t), are used; and the kinetic evolution in Equation 25 will only require the solution of gi . The ?rst two de?nitions in Equation 27 can be approximated using the discrete velocities in a Gaussian-type quadrature: g eq = ρ(x, t) =
Wi gi (x, t),
ρu(x, t) =
Wi ei gi (x, t).
LATTICE BOLTZMANN METHOD
This de?nes an effective distribution function f i (x, t): f i (x, t) = Wi gi (x, t), which satis?es the simple conservation relations in Equation 3. Given ei , Wi is constant, and f i satis?es the same equation as g: ? fi + ei · ?t
1 eq f ? fi , ετ i
with f i = wi ρ[1 + 3(ei · u) + 9 (ei · u)2 ? 3 u 2 ] and wi = Wi /(2π/3) D/2 2 2 exp(? 3 ei2 ). 2 To obtain the weight wi , He & Luo (1997) use a third-order Hermite formula to approximate the integrals in Equation 27. Abe (1997) assumes wi has a simple truncated functional form based on ei . For the nine-velocity square lattice, both papers found that w0 = 4/9, wi (i = 1, 3, 5, 7) = 1/9 and wi (i = 2, 4, 6, 8) = 1/36, which give the same equilibrium distribution function as the original lattice Boltzmann model in Equation 22. If in Equation 31 the time derivative is replaced by a ?rst order time difference, and the ?rst order upwind discretization for the convective term ei · f i is used and a downwind collision term (x ? ei , t) for (x, t) is used (Cao et al 1997), then we obtain the ?nite difference equation for f i : f i (x, t + t) = f i (x, t) ? α f i (x, t) ? f i (x ? ? xei , t) β eq (32) f i (x ? xei , t) ? f i (x ? xei , t) , τ where α = t|ei |/ x, β = t/ε, and t and x are the time step and the grid step, respectively. Choosing α = 1 and β = 1, Equation 32 becomes the standard LBE as shown in Equation 2. From the discretization process, we note Equation 32 only has ?rst order convergence in space and time to Equation 31. However, as shown by Sterling and Chen (1995), because Equation 32 has a Lagrangian nature in its spatial discretization, the discretization error has a special form which can be included in the viscous term, resulting in second-order accuracy both in space and time.
Boundary Conditions in the LBM
Wall boundary conditions in the LBM were originally taken from the LG method. For example, a particle distribution function bounce-back scheme (Wolfram 1986, Lavall? e et al 1991) was used at walls to obtain no-slip vee locity conditions. By the so-called bounce-back scheme, we mean that when a particle distribution streams to a wall node, the particle distribution scatters back to the node it came from. The easy implementation of this no-slip velocity condition by the bounce-back boundary scheme supports the idea that the
CHEN & DOOLEN
LBM is ideal for simulating ?uid ?ows in complicated geometries, such as ?ow through porous media. For a node near a boundary, some of its neighboring nodes lie outside the ?ow domain. Therefore the distribution functions at these no-slip nodes are not uniquely de?ned. The bounce-back scheme is a simple way to ?x these unknown distributions on the wall node. On the other hand, it was found that the bounce-back condition is only ?rst-order in numerical accuracy at the boundaries (Cornubert et al 1991, Ziegler 1993, Ginzbourg & Adler 1994). This degrades the LBM, because numerical accuracy of the LBM in Equation 2 for the interior mesh points is second-order. He et al (1997) con?rmed this result by analyzing the slip velocity near the wall node for Poiseuille ?ow. To improve the numerical accuracy of the LBM, other boundary treatments have been proposed. Skordos (1993) suggested including velocity gradients in the equilibrium distribution function at the wall nodes. Noble et al (1995) proposed using hydrodynamic boundary conditions on no-slip walls by enforcing a pressure constraint. Inamuro et al (1995) recognized that a slip velocity near wall nodes could be induced by the bounce-back scheme and proposed to use a counter slip velocity to cancel that effect. Maier et al (1996) modi?ed the bounce-back condition to nullify net momentum tangent to the wall and to preserve momentum normal to the wall. Zou and He (1997) extended the bounce-back condition for the nonequilibrium portion of the distribution. Ziegler (1993) noticed that if the boundary was shifted into ?uid by one half mesh unit, i.e. placing the nonslip condition between nodes, then the bounceback scheme will give second-order accuracy. Simulations demonstrated that these heuristic models yield good results for ?uid ?ows around simple wall boundaries. The above boundary treatments were also analyzed by Zou et al (1995) and He et al (1997) by solving the lattice BGK equation (Equation 2). It appears, however, that the extension of these simple assumptions to arbitrary boundary conditions is dif?cult. Chen et al (1996) viewed the LBM as a special ?nite difference scheme of the kinetic equation (Equation 31). They adopted staggered mesh discretization from traditional ?nite difference methods and proposed using a second-order extrapolation scheme of distributions in the ?ow to obtain the unknown particle distribution functions. Their extrapolation scheme is simple and can be extended to include velocity, temperature, and pressure (Maier et al 1996, Zou & He 1997) boundary conditions and their derivatives. In this treatment, boundary conditions can be assigned to grid nodes or to any locations in space for little computational effort. Numerical simulations including time-dependent Couette ?ow, a lid-driven square-cavity ?ow, and ?ow over a column of cylinders were carried out, showing good agreement with analytical solutions and ?nite difference solutions.
LATTICE BOLTZMANN METHOD
Numerical Ef?ciency, Stability, and Nonuniform Grid
The discrete velocity equation (Equation 31) is a hyperbolic equation that approximates the Navier-Stokes equation in the nearly-incompressible limit (Majda 1984). The numerical accuracy depends on Mach number (Reider & Sterling 1995). The advective velocity ei in Equation 31 is a constant vector, in contrast to the spatial dependent velocity in compressible hydrodynamic equations, which prevents the LBM from solving nonlinear Riemann problems. From a numerical analysis point of view, the LBM, like other kinetic equations, is a relaxation method (Cao et al 1997) for the macroscopic equations, which has much in common with the explicit “penalty” or “pseudocompressibility method” (Sterling & Chen 1996). This view was used by Ancona (1994) to generalize the LBM to include fully Lagrangian methods for the solution of partial differential equations. The advective velocity ei in Equation 31 is constant in contrast to the spatial dependent velocity in compressible hydrodynamic equations. The kinetic-relaxation method for solving a hyperbolic conservation system was proposed by Jin and Xin (1995). This approach uses the relaxation approach to model the nonlinear ?ux terms in the macroscopic equations and thus it does not require nonlinear Riemann solvers either. Using this relaxation method, Jin and Katsoulakis (1997) calculated curvature-dependent front propagation. In principle, the LBM and the kinetic-relaxation method are very much alike. The kinetic-relaxation was developed mainly for shock capture in Euler systems, whereas the LBM is more focused on viscous complex ?ows in the nearly-incompressible limit. Nadiga and Pullin (1994) proposed a simulation scheme for kinetic discrete-velocity gases based on local thermodynamic equilibrium. Their method seems more general and able to obtain high order numerical accuracy. Their ?nite volume technique was further developed (Nadiga 1995) to solve the compressible Euler equation by allowing the discrete velocities to adapt to the local hydrodynamic state. Elton et al (1995) studied issues of convergence, consistency, stability, and numerical ef?ciency for lattice Boltzmann models for viscous Burgers’ equation and advection-diffusion system. The numerical ef?ciency of the LBM was studied by several authors. Succi et al (1991) noted that the number N f po of ?oating point calculations for the LBM would be ?r N D , while N f po ? (25 log2 N )N D for the pseudospectral method, where N is the number of lattice points along each direction in D dimensions. r is 150 for the two-dimensional LBM using the full matrix Mi j in Equation 10 and 40 for the LBGK model. For 3-D low Reynolds number simulations, an LBM performance of 2.5 times faster than the pseudo-spectral method was reported (Chen et al 1992). The scalability of the LBM was studied by Noble et al (1996). They found a better-than-linear scalability for the LBM
CHEN & DOOLEN
when the number of processors increases while simultaneously increasing size of the computational domain. For a 2-D problem with 21762 lattice points, they achieved an overall 13.9 giga?ops using 512 processors of CM-5 machines. In the continuum Boltzmann equation, the equilibrium distribution function given by Equation 26 corresponds to the maximum entropy state. Thus, any initial state will evolve towards a state of higher entropy. This result is known as Boltzmann’s H-theorem, which ensures an increase of entropy and ensures stability. An H-theorem has been derived also for the LG method (Frisch et al 1987). In the LBM, however, only a small number of discrete velocities is used and usually the equilibrium distribution is truncated to O(u 2 ) (see Equation 21). Consequently, one cannot simultaneously guarantee an H-theorem and allow the correct form of the macroscopic equations. Therefore, the LBM, although of a kinetic nature, is subject to numerical instability. Furthermore, we know also that the traditional LBM equation (Equation 2) is a ?nite-difference solution of the discrete-velocity Boltzmann equation in Equation 31 and the discretization in Equation 32 with α = 1 marginally satis?es the CourantFriedrichs-Lewy condition. A von Neumann linearized stability analysis of the LBM was carried out by Sterling & Chen (1996). In this analysis, they expanded f i as f i (x, t) = f i(0) + f i (x, t), where the global equilibrium populations f i(0) were assumed to be constants that depend only on a constant mean density and constant velocity. Taylor-expanding Equation 2, they arrived at the following linearized equation: f i (x + ei t, t + t) = Ji j f j (x, t), (33)
where Ji j is the Jacobian matrix corresponding to the coef?cient of the linear term in Equation 2 and does not depend on location or time. Spatial dependence of the stability analysis was carried out by taking the Fourier transformation of Equation 33 and solving the eigenvalue problem for diag[exp(?ik · ei t)]Ji j . For k = 0, they found the linear stability condition: τ ≤ 1/2 which is consistent with the positivity of viscosity. Detailed stability analysis were carried out for the hexagonal 2-D LBM model (Chen et al 1992), the square 2-D LBM model and the 3-D 15-velocity model (Qian et al 1992). They concluded that a stable LBM requires that the mean ?ow velocity be below a maximum velocity that is a function of several parameters, including the sound speed Cs , the relaxation time τ , and the wave number. The numerical stability can also be improved if one starts from the continuum kinetic equation (Equation 31) (Ancona 1994, McNamara et al 1995, Cao et al 1997) using different ?nite difference discretizations. Until recently, one of the limitations to the numerical ef?ciency in the LBM, as compared with other CFD methods, is that the discretization of Equation 31 was constrained on a special class of uniform and regular lattices. To increase
LATTICE BOLTZMANN METHOD
numerical ef?ciency and accuracy, nonuniform grid methods have been developed. Koelman (1991) suggested including spatial nonuniformity in the lattice structure. Nannelli & Succi (1992) and Amati et al (1997) developed the ?nite volume LBM (FVLBM) based on Equation 31. They de?ned a nonuniform coarse lattice whose cell typically contained several original lattice units. The ?1 evolution equation for the mean value Fi ≡ VC C f i d 3 x in the Cell C requires the evaluation of the ?ux across the boundaries of C, where VC is the volume. A piece-wise constant or a piece-wise linear interpolation was used to approximate the ?ux. This FVLBM has been used to study two-dimensional ?ow past a bluff body (Succi & Nannelli 1994) and three-dimensional turbulent channel ?ow (Amati et al 1997). The results of the 3-D channel ?ow simulation reveal that the FVLBM is about one order of magnitude faster than the traditional nonuniform CFD methods. On the other hand, while most simulation results agree well with other traditional methods, the comparison of simulation results is correspondingly less satisfactory, owing partially to the low-order interpolation schemes. Cao et al (1997) developed a nonuniform ?nite difference LBM for Equation 31 and simulated annular ?ows in polar coordinates. He et al (1996) adhered to the LBM equation (Equation 2) and developed a time-dependent interpolation scheme that increases grid density in high-shear regions. They successfully simulated ?ows in the 2-D symmetric channel with sudden expansion and ?ow around a circular cylinder (He & Doolen 1997).
LATTICE BOLTZMANN SIMULATION OF FLUID FLOWS Flows with Simple Boundaries
DRIVEN CAVITY FLOWS For two-dimensional cavity ?ows, there are many papers using traditional schemes, including ?nite-difference (Schreiber & Keller 1983, E & Liu 1996) and multigrid (Vanka 1986). The fundamental characteristics of the 2-D cavity ?ow are the emergence of a large primary vortex in the center and two secondary vortices in the lower corners. The values of the stream function and the locations of the centers of these vortices as a function of Reynolds numbers have been well studied. The lattice Boltzmann simulation of the 2-D driven cavity by Hou et al (1995) covered a wide range of Reynolds numbers from 10 to 10,000 using a 2562 lattice. They carefully compared simulation results of the stream function and the locations of the vortex centers with previous numerical simulations and demonstrated that the differences of the values of the stream function and the locations of the vortices between the LBM and other methods were less than 1%. This difference is
CHEN & DOOLEN
within the numerical uncertainty of the solutions using other numerical methods. The spatial distributions of the velocity, pressure and vorticity ?elds were also carefully studied. An error analysis of the compressibility effect from the model was carried out. Two-dimensional cavity ?ow was also studied by Miller (1995). Instead of using a uniform velocity on the top of the cavity, in Miller’s work a shear-?ow condition with certain velocity distribution was given which allowed an analytical solution of the velocity for the whole cavity. Excellent agreement between the LBM simulations and the analytical solutions for the velocity and pressure ?elds was found. The effects of the variation of the cavity aspect ratio on vortex dynamics were also studied. Hou (1995) simulated three-dimensional cubic cavity ?ow using the threedimensional 15-velocity LBM (Qian et al 1992, Chen et al 1992). 1283 lattice points were used with Re = 3,200. Flows at this Reynolds number had been extensively studied earlier, including the ?nite-difference simulation and the experimental work by Prasad & Koseff (1989). Flow structures, including the velocity and vorticity ?elds in different planes were analyzed using the LBM simulation. They compared well with previous numerical studies. The correlation of the Taylor-G¨ rtler-like vortices in the yz- and x z-planes was o reported for the ?rst time. Figure 1 displays the mean velocity pro?les from 3-D LBM, 2-D LBM simulations, and experimental work (Prasad & Koseff 1989) in the symmetry plane along the vertical and the horizontal centerlines (left), and the root-mean-square (rms) velocity pro?les (right). The agreement
? u v ? Figure 1 Mean velocity pro?les (left), U and U , in the symmetry plane along the vertical and the horizontal centerlines. Solid line is for the 3-D LBE simulation, dotted line is a 2-D simulation (Hou 1995), and circles represent experimental results (Prasad & Koseff 1989). rms velocities (right), Ur ms and Vr ms , in the symmetry plane along vertical and horizontal centerlines. Solid line is the LBE simulation, upward triangles and downward triangles are experimental results (Prasad & Koseff 1989).
LATTICE BOLTZMANN METHOD
between the 3-D LBM results and experimental results shows that the LBM is capable of simulating complex 3-D unsteady ?ows.
FLOW OVER A BACKWARD-FACING STEP The two-dimensional symmetric sudden expansion channel ?ow was studied by Luo (1997) using the LBM. The main interest in Luo’s research was to study the symmetry-breaking bifurcation of the ?ow when Reynolds number increases. In this simulation, an asymmetric initial perturbation was introduced and two different expansion boundaries, square and sinusoidal, were used. This simulation reproduced the symmetricbreaking bifurcation for the ?ow observed previously, and obtained the critical Reynolds number of 46.19. This critical Reynolds number was compared with earlier simulation and experimental results of 40.45 and 47.3, respectively. Qian et al (1996) used the ?ow over a backward-facing step as a benchmark for validating the LBM. They studied the recirculation length as a function of Reynolds number and the channel expansion ratio. The results from the LBM compared well with previous results using ?nite-difference methods, spectralelement methods, and experiments. The appearance of the second vortex and the transition leading to turbulence were also studied. FLOW AROUND A CIRCULAR CYLINDER The ?ow around a two-dimensional circular cylinder was simulated using the LBM by several groups of people. The ?ow around an octagonal cylinder was also studied (Noble et al 1996). Higuera & Succi (1989) studied ?ow patterns for Reynolds number up to 80. At Re = 52.8, they found that the ?ow became periodic after a long initial transient. For Re = 77.8, a periodic shedding ?ow emerged. They compared Strouhal number, ?ow-separation angle, and lift and drag coef?cients with previous experimental and simulation results, showing reasonable agreement. Pressure distributions around the surface of the cylinder as a function of Reynolds number were carefully studied by Wagner (1994). He concluded that the difference between his Strouhal numbers and those from other simulations was less than 3.5%. Considering the compressible nature of the LBM scheme, this result is signi?cant. It demonstrates that in the nearly incompressible limit, the LBM simulates the pressure distribution for incompressible ?uids well (Martinez et al 1994). The above two simulations were based on uniform lattices, which poorly approximate circular geometry. The no-slip boundary was enforced at the boundary mesh point which is jagged and not necessary on the cylinder. He & Doolen (1997) revisited the problem of two-dimensional ?ow around a circular cylinder based on the interpolation-supplement LBM (He et al 1996). In this simulation, the underlying lattice was square, but a spatial interpolation was used. A polar coordinate was constructed and the simulation domain consisted of a circular region. The no-slip wall boundary condition was exactly enforced
CHEN & DOOLEN
on the cylinder boundary. It is clear that the computational accuracy in this simulation is better than previous LBM simulations. The streaklines showed vortex shedding quite similar to experimental observation. The LBM calculated Strouhal numbers. Lift and drag coef?cients for ?ows at different Reynolds numbers compared well with previous simulations using other schemes. The error was within experimental uncertainty and is the same as errors among other schemes.
Flows in Complex Geometries
An attractive feature of the LBM is that the no-slip bounce-back LBM boundary condition costs little in computational time. This makes the LBM very useful for simulating ?ows in complicated geometries, such as ?ow through porous media, where wall boundaries are extremely complicated and an ef?cient scheme for handing wall-?uid interaction is essential. The fundamental problem in the simulation of ?uid ?ows through porous media is to understand and model the macroscopic transport from microscopic ?uid dynamics. A starting point for this problem is Darcy’s law, which assumes a linear relation between the pressure gradient, P/L, and the volume ?ow rate per unit area q: q = ?K /(ρ0 ν) P/L. Here K is the permeability. Previous numerical simulations, including ?nite-difference schemes (Schwarz et al 1994) and networking models (Koplik & Lasseter), were either limited to simple physics, small geometry size, or both. Lattice gas automata were also used for simulating porous ?ows and verifying Darcy’s law in simple and complicated geometries (Rothman 1988). Succi et al (1989) utilized the LBM to measure the permeability in a 3-D random media. Darcy’s law was con?rmed. Cancelliere et al (1990) studied the permeability K as a function of solid fraction η in a system of randomly positioned spheres of equal radii. The simulation covered a wide range of η (0.02 ≤ η ≤ 0.98). The values of K from the LBM compared well with the Brinkman approximation: K = K 0 [1 + 3/4η(1 ? √ 8/η ? 3)] for η < 0.2, and agreed well with the semiempirical KozenyCarman equation: K = (1 ? η)3 /(6s 2 ). Here K 0 is the effective permeability for a single sphere and s is the speci?c surface area for the system. Later, Heijs & Lowe (1995) con?rmed this result independently. In addition, they studied the validity of the Kozeny-Carman equation for soil samples where ?ow occurs only through some speci?c continuous connected pore and ?ows occurring at smaller scales are negligible. For this condition, the Kozeny-Carman equation provided a much less successful estimate of the permeability. Flows through sandstones measured using X-ray microtomography were simulated by Buckles et al (1994), Soll et al (1994), and Ferr? ol & Rothman (1995). They found that e the permeability for these sandstones, although showing large variation in space and ?ow directions, in general agreed well with experimental measurements
LATTICE BOLTZMANN METHOD
within experimental uncertainty. Ferr? ol & Rothman (1995) also studied the e effect of grid resolution on permeability. They found that the permeability strongly depends on the viscosity when the minimum channel in the pores is several grid points wide. Spaid & Phelan (1997) investigated the injection process in resin transfer molding. For this heterogeneous porous media simulation, they used the LBM simulation of full Navier-Stokes ?ows to model a ?ow around a circular cylinder. Inside the cylinder, the velocity is governed by the Brinkman equation (ν 2 u ? βu = P/ρ, where β is a model parameter). Excellent agreement between the LBM simulation and lubrication theory for cell permeabilities was reported.
Simulation of Fluid Turbulence
DIRECT NUMERICAL SIMULATION A major difference between the LBM and the LG method is that the LBM can be used for smaller viscosities. Consequently the LBM can be used for direct numerical simulation (DNS) of high Reynolds number ?uid ?ows. To validate the LBM for simulating turbulent ?ows, Martinez et al (1994) studied decaying turbulence of a shear layer using both the pseudospectral method and the LBM. The initial shear layer consisted of uniform velocity reversing sign in a very narrow region. The initial Reynolds number was 10,000. The simulation used a 5122 mesh and ran for 80 largeeddy turnover times. Martinez et al carefully compared the spatial distribution, time evolution of the stream functions, and the vorticity ?elds. Energy spectra as a function of time, small scale quantities, including enstrophy, palinstrophy, and 4th-order enstrophy as a function of time were also studied. The correlation between vorticity and stream function was calculated and compared with theoretical predictions. They concluded that the LBE method provided a solution that was “accurate” in the sense that time histories of global quantities, wavenumber spectra, and vorticity contour plots were very similar to those obtained from the spectral method. In particular, they reproduced details of the wavenumber spectra at high wavenumbers as well as the detailed structure of vortex distributions. Their conclusions are signi?cant for simulations of turbulence in which small-scale dynamics are important. In that paper they found discrepancies in the spatial locations of vortical structures at very late stages. Figure 2 (top left) compares the energy spectra from a pseudospectral method (dotted line) and the LBM (solid line) for t = 5, and displays the time evolution of the enstrophy function (bottom left). Figure 2 also depicts contour plots of the vorticity function at the same time step (right). The agreement of these two methods is excellent. Two-dimensional forced isotropic turbulence was simulated by Benzi & Succi (1990) to study the enstrophy cascade range. Twodimensional forced turbulence was also simulated by Qian et al (1995) to study
CHEN & DOOLEN
Figure 2 Velocity energy spectra versus wavenumber, k (top left). The enstrophy of the system, ω2 (bottom left), as a function of time. Right, a comparison of the vorticity distributions from a pseudospectral method and the LBM (Martinez et al 1994).
the energy inverse cascade range. They reproduced the k ?5/3 inertial range scaling, in good agreement with theoretical predictions (Kraichnan 1967). Chen et al (1992) made a similar validation of LBM by comparing results of three-dimensional isotropic turbulence using the LBM and the pseudospectral method. In that work, they simulated three-dimensional Beltrami ?ows, the decaying Taylor-Green vortex, and decaying three-dimensional turbulence. The agreement between the two methods for spatial and time distributions of velocities and vortices was good. This assessment was echoed by Trevi? o & n Higuera (1994), who studied the nonlinear stability of Kolmogorov ?ows using the pseudospectral method and the LBM at various Reynolds numbers. The LBM is a useful direct numerical simulation (DNS) tool for simulating nonhomogeneous turbulent ?ows. To validate the generalized extended selfsimilarity for anisotropic ?ows where the simple extended self-similarity is not
LATTICE BOLTZMANN METHOD
valid, Benzi et al (1996) used the LBM to simulate three-dimensional shear ?ow where a sinusoidal force in the x direction varying along the z direction was applied. This simulation generated useful velocity and vorticity information, and the scaling exponents of the velocity increment and their dependence on the shear rate were studied. Based on their analysis of this ?ow, they concluded that the generalized extended self-similarity was valid for anisotropic ?ows. Succi et al (1991) studied the bifurcation of a two-dimensional Poiseuille ?ow and obtained the amplitude of the primary and secondary bifurcated models as a function of the Reynolds number. Eggels (1996) utilized the LBM on the FCHC lattice (d’Humi` res et al 1986) as a DNS method for simulating threee dimensional channel ?ows with Re? = 180 (based on the friction velocity). Their simulation results, including mean stream-wise velocity pro?les and rms velocities as functions of the distance from the wall, were compared with those reported by Kim et al (1987), showing good agreement.
LBM MODELS FOR TURBULENT FLOWS As in other numerical methods for solving the Navier-Stokes equations, a subgrid-scale (SGS) model is required in the LBM to simulate ?ows at very high Reynolds numbers. Direct numerical simulation is impractical due to the time and memory constraints required to resolve the smallest scales (Orszag & Yakhot 1986). Hou et al (1996) directly applied the subgrid idea in the Smagorinsky model to the LBM by ?ltering the particle distribution function and its equation in Equation 2 using a standard box ?lter. To account for the nonlinear term in the equilibrium distribution (Equation 22), Hou et al (1996) replaced the relaxation time τ in Equation 13 2 by (3ν0 + Cs 2 |S|2 ) + 1/2, where ν0 = (2τ ? 1)/6 is the kinetic viscosity, is the ?lter width, |S| is the amplitude of the ?ltered large-scale strain-rate tensor, and Cs is the Smagorinsky constant. This simple Smagorinsky-type LBM subgrid model was used to simulate ?ows in a two-dimensional cavity at Reynolds numbers up to 106 using a 2562 mesh (Hou et al 1996). Threedimensional turbulent pipe ?ow was simulated by Somers (1993) using the LBM SGS model at a Reynolds number of 50,000 for mesh sizes up to 803 . Friction as a function of Reynolds number compared reasonably well with the Blasius power law for turbulent ?ows. Succi et al (1995) suggested using the relaxation idea in the LBM to solve the turbulent K ? equations. Hayot & Wagner (1996) used a dependent relaxation τ (t) in the lattice BGK equation (Equation 13), which leads to an effective turbulent viscosity. Eggels (1996) and Eggels & Somers (1995) included the turbulent stress tensor σt directly eq into the equilibrium distribution, f i in Equation 22, by adding an extra term: M σ f i = ρ[eiα eiβ σt αβ ? D tr(σt )]. Using this model, Eggels (1996) carried out a large-eddy simulation of the turbulent ?ow in a baf?ed stirred tank reactor. The impact of a mechanical impeller was modeled via a spatial and temporal
CHEN & DOOLEN
dependent force on the momentum equation. It was reported that the subgrid LBE turbulent model coupled with the impeller model produced satisfactory results. Both mean ?ow and turbulent intensities compared well with experimental data. This simulation demonstrated the potential of the LBM SGS model as a useful tool for investigating turbulent ?ows in industrial applications of practical importance.
LBM SIMULATIONS OF MULTIPHASE AND MULTICOMPONENT FLOWS
The numerical simulation of multiphase and multicomponent ?uid ?ows is an interesting and challenging problem because of dif?culties in modeling interface dynamics and the importance of related engineering applications, including ?ow through porous media, boiling dynamics, and dendrite formation. Traditional numerical schemes have been successfully used for simple interfacial boundaries (Glimm et al 1981, Brackbill et al 1992, Chang et al 1996). The LBM provides an alternative for simulating complicated multiphase and multicomponent ?uid ?ows, in particular for three-dimensional ?ows.
Method of Gunstensen et al
Gunstensen et al (1991) were the ?rst to develop the multicomponent LBM method. It was based on the two-component LG model proposed by Rothman & Keller (1988). Later, Grunau et al (1993) extended this model to allow variations of density and viscosity. In these models, red and blue particle distribution functions f i(r ) (x, t) and f i(b) (x, t) were introduced to represent two different ?uids. The total particle distribution function (or the color-blind particle distribution function) is de?ned as: f i = f i(r ) + f i(b) . The LBM equation is written for each phase: f ik (x + ei , t + 1) = f ik (x, t) +
k i k i (x, t),
were k denotes either the red or blue ?uid, and =
k 1 i
k 2 i
is the collision operator. The ?rst term in the collision operator, ( ik )1 , represents the process of relaxation to the local equilibrium similar to the LBGK model in Equation 13:
k 1 i
?1 k k(eq) f ? fi . τk i
is the local equilibrium distribution depending on the local macroHere, f i scopic variables ρ (k) and u(k) . τk is the characteristic relaxation time for species
LATTICE BOLTZMANN METHOD
k. The viscosity of each ?uid can be selected by choosing the desired τk . Conservation of mass for each phase and total momentum conservation are enforced at each node during the collision process: ρr =
f ir =
, ei ,
f ib =
f ik ei =
where ρ = ρr + ρb is the total density and ρu is the local total momentum. k(eq) The form of f i can be chosen to be similar to Equation 22. The additional collision operator ( k )2 contributes to the dynamics in the interfaces and generates a surface tension: Ak k 2 |F| (ei · F)2 /|F|2 ? 1/2 , = (38) i 2 where F is the local color gradient, de?ned as: F(x) =
ei (ρr (x + ei ) ? ρb (x + ei )).
Note that in a single-phase region of the incompressible ?uid model, F vanishes. Therefore, the second term of the collision operator ( ik )2 only contributes to interfaces and mixing regions. The parameter Ak is a free parameter, which determines the surface tension. The additional collision term in Equation 38 does not cause the phase segregation. To maintain interfaces or to separate the different phases, the LBM by Gunstensen et al follows the LG method of Rothman & Keller (1988) to force the local color momentum, j = i ( f ir ? f ib )ei , to align with the direction of the local color gradient after collision. In other words, the colored distribution functions at interfaces were redistributed to maximize ?j · F. Intuitively, this step will force colored ?uids to move toward ?uids with the same colors. The multicomponent LBM by Gunstensen et al (1991) has two drawbacks. First, the procedure of redistribution of the colored density at each node requires time-consuming calculations of local maxima. Second, the perturbation step with the redistribution of colored distribution functions causes an anisotropic surface tension that induces unphysical vortices near interfaces. D’ortona et al (1994) modi?ed the model of Gunstensen et al. In their model, the recoloring step is replaced by an evolution equation for f ik that increases computational ef?ciency.
Method of Shan & Chen
Shan & Chen (1993) and Shan & Doolen (1995) used microscopic interactions to modify the surface-tension–related collision operator for which the surface
CHEN & DOOLEN
k 2 i)
interface can be maintained automatically. In these models, ( 38 was replaced by the following formula:
k 2 i k
in Equation (40)
= ei · Fk .
Here F is an effective force on the kth phase owing to a pairwise interaction between the different phases: Fk (x) = ?
Vkk (x, x + ei )ei .
Here Vkk is an interaction pseudopotential between different phases (or components): Vkk (x, x ) = G kk (x, x )ψ k (x)ψ k (x ).
G kk (x) is the strength of the interaction; and ψ (x) is a function of density for the k phase at x which has taken the following empirical form: ψ = ρ0 [1 ? exp(?ρ/ρ0 )], where ρ0 is a constant free parameter. It was shown that this form of the effective density ψ gives a non-ideal equation of state, which separates phases or two component ?uids. Qian et al (1995) have recently demonstrated that other pseudopotential forms, including fractional potential and a van der Waals potential, will produce similar results. It should be noted that the collision operator ( ik )2 in the Shan-Chen model does not satisfy local momentum conservation. This treatment is physically plausible because the distant pairwise interactions between phases change the point-wise local momenta at the positions involved in the interactions. This feature differs from the model of Gunstensen et al where the total local momentum is conserved (see the third operation in Equation 37). It is arguable that this spurious conservation might be one reason why the model exhibits unphysical features near interfaces. For simplicity, only the nearest-neighbor interactions were involved in the Shan-Chen model, in which G kk (x, x ) is constant when x ? x = ei , and zero otherwise. G kk acts like a temperature; when G is smaller than the critical value G c (depending on the lattice structure and initial density), the ?uids separate. Theoretical calculations indicate (Shan & Chen 1994) that the surface tension σ ? M G/[2D(D + 1)], which was also veri?ed by numerical simulation. In the Shan-Chen model, the separation of ?uid phases or components is automatic (Chen 1993). This is an important improvement in numerical ef?ciency compared with the original LBM multiphase models. The Shan-Chen model also improves the isotropy of the surface tension.
Free Energy Approach
The above multiphase and multicomponent lattice Boltzmann models are based on phenomenological models of interface dynamics and are probably most
LATTICE BOLTZMANN METHOD
suitable for isothermal multicomponent ?ows. One important improvement in models using the free-energy approach (Swift et al 1995, 1996) is that the equilibrium distribution can be de?ned consistently based on thermodynamics. Consequently, the conservation of the total energy, including the surface energy, kinetic energy, and internal energy can be properly satis?ed (Nadiga & Zaleski 1996). The van der Waals formulation of quasilocal thermodynamics for a twocomponent ?uid in thermodynamic equilibrium at a ?xed temperature has the following free-energy functional: (r) = dr[ψ(T, ρ) + W ( ρ)]. (43)
The ?rst term in the integral is the bulk free-energy density, which depends on the equation of state. The second term is the free-energy contribution from density gradients and is related to the surface tension. For simple multiphase ?uid ?ows, W = κ ( ρ)2 , where κ is related to the surface tension. The full 2 pressure tensor in a nonuniform ?uid has the following form: Pαβ = pδαβ + κ ?ρ ?ρ . ? xα ? xβ (44)
The pressure is obtained from the free energy: p(r) = ρ δδρ ? (r) = p0 ? κρ 2 ρ ? κ ( ρ)2 , where p0 = ρψ (ρ) ? ψ(ρ) is the equation of state. For a 2 van der Waals equation of state, ψ = ρT In ρ/(1 ? ρb) ? aρ 2 , where a and b are free parameters. To incorporate the above formulation into the lattice Boltzmann equation, Swift et al (1995) added a term to the original equilibrium distribution in Equaeq eq eq tion 21: f i = f i + G αβ eiα eiβ . They vary the coef?cients in f i (see Equation 21) and G αβ to guarantee the conservation of mass and momentum and to eq obtain the following stress tensor condition: i f i eiα eiβ = Pαβ + ρu α u β .
Numerical Veri?cation and Applications
Two fundamental numerical tests associated with interfacial phenomena have been carried out using the multiphase and multicomponent lattice Boltzmann models. In the ?rst test, the lattice Boltzmann models were used to verify Laplace’s formula by measuring the pressure difference between the inside and the outside of a droplet: Pinner ? Pouter = σ , where Pinner and Pouter denote R the pressures inside and outside of the droplet, respectively. R is the radius of the droplet and σ is the surface tension. The simulated value of σ has been compared with theoretical predictions, and good agreement was reported (Gunstensen et al 1991, Shan and Chen 1993, Swift et al 1995). Alternatively,
CHEN & DOOLEN
the surface tension can also be measured using the mechanical de?nition σ =
(PN ? PT )dz,
where PN and PT are the normal and the tangential pressure along a ?at interface, respectively, and the integration is evaluated in a direction perpendicular to the interface. It is shown from LBM simulations (Gunstensen et al 1991, Grunau et al 1993, Swift et al 1995) that the surface tension obtained from the mechanical de?nition agrees with the Laplace formula de?nition. In the second test of LBM interfacial models, the oscillation of a capillary wave was simulated (Gunstensen et al 1991, Shan & Chen 1994, Swift et al 1995). A sine wave displacement of a given wave vector was imposed on an interface that had reached equilibrium. The resulting dispersion relation was measured and compared with the theoretical prediction (Laudau & Lifshitz 1959): ω2 = σ/ρk 3 . Good agreement was observed, validating the LBM surface tension models. Because the lattice Boltzmann multiphase and multicomponent models do not track interfaces, it is easy to simulate ?uid ?ows with very complicated interfaces, such as the domain growth in the spinodal decomposition process for binary ?uids (Alexander et al 1993a, Rybka et al 1995). In these simulations, hydrodynamics plays an important role and multiphase or multicomponent ?uids evolved from an initial mixed state owing to surface tension, viscosity, and inertial force. A convenient way to characterize the growth kinetics is to use the correlation function of the order parameter, G(r, t) ≡ φ(r )φ(0) , where φ is the order parameter, de?ned as (ρr ? ρb )/(ρr + ρb ). The Fourier transform of G(r, t) is then the structure function S(k,t). As time evolves, the structure function becomes more sharply peaked, and its maximum value moves toward small k. In a wide variety of phase segregating systems, S(k,t) has a simple form at late times: S(k, t) = R D (t)F(x). Here R(t) is the average domain size, D is the dimension, and F(x) is a universal function depending on x = k R(t). Numerical calculations of R(t) and F(x) have been carried out using the LBM (Alexander et al 1993a). It was found that R(t) scales ?t 2/3 in two dimensions and ?t in three dimensions for two-component ?uids. This agrees with theoretical predictions. F(x) was found to agree well with Porod’s law for large x, ?x D+1 , but small scales showed strong dependence of ?ow properties, which was later con?rmed by Wu et al (1995) using the Langevin equation. Osborn et al (1995) studied similar scaling problems for liquid-gas two-phase domain growth. They found that R(t) ? t 1/2 and ?t 2/3 at high and low viscosities, respectively. A similar crossover for two-component ?uids from t 1/3 to t 2/3 was also reported as time increased. The hydrodynamic effects on the
LATTICE BOLTZMANN METHOD
spinodal decomposition for three and four components in two and three dimensions in critical and off-critical quenching were also studied by Chen & Lookman (1995). Cieplak (1995) carried out a study of two-dimensional rupture dynamics under initial perturbation and compared results with molecular dynamics simulations. The coalescence process was also studied for a droplet in a background ?uid falling in a gravitational ?eld to the bottom of a container in which the bottom is a “bare” wall, a shallow liquid with the same kind as the droplet, or a deeper liquid. Interesting interfacial structures were observed. In Figure 3 we show condensation and subsequent coalescence of liquid drops in supersaturated vapor at two different times. The LBE model with interparticle interaction (Shan & Chen) was used in this simulation on a 1283 lattice. The ?uid obeys the equation of state: p = 1 ρ + 3 G(1 ? e?ρ )2 . G (chosen to 3 2 be ?0.55) is a parameter that gives the relative strength of the interaction. This equation of state gives a non-monotonic isotherm with a critical point at G = ? 4 . The initial density is constant except for a randomly distributed small 9 perturbation. One feature of the free-energy approach is that one can parameterize the free energy to simulate complicated ?uid ?ows, including lamellar ?uids. In their study, Gonnella et al (1997) used the following Ginzburg-Landau potential for the ψ function in Equation 43 based on the order parameter, φ: ψ = λ 2 φ + θ φ 4 . They included high-order gradients in the W function: W = 2 4 κ ( φ)2 + ζ ( 2 φ)2 where θ and ζ are positive. λ < 0 and k < 0 lead to 2 2 lamellar phases. The magnitude of ζ is directly connected to the elasticity of
Figure 3 Condensation and subsequent coalescence of liquid drops in supersaturated vapor. Early, left, late, right (Shan 1997).
CHEN & DOOLEN
Figure 4 Con?guration of a lamellar ?uid when a shear ?ow is applied. The shear velocities are v = 0 (a), v = 0.1 (b), and v = 0.3 (c) (Gonnella et al 1997).
the ?uid. In Figure 4, we present ?ow patterns for a lamellar ?uid under a shear induced by forcing a constant velocity v in the left and ?v in the right simulation domains, respectively. As the shear was applied, ordered lamellae quickly formed, aligned at an angle to the direction of shear. A further increase of the shear rate resulted in a decrease in the angle between the direction of the lamellae and that of the shear. This pattern directly re?ects the balance between the shear force and elastic force of the lamellar ?uids. It would be of great interest to use this model for simulating bio-?uids where elasticity is important. Lattice Boltzmann multiphase ?uid models have been extensively used to simulate multicomponent ?ow through porous media in order to understand the fundamental physics associated with enhanced oil recovery, including relative permeabilities (Vangenabeek et al 1996). The LBM is particularly useful for this problem because of its capability of handling complex geometrical boundary conditions and varying physical parameters, including viscosities (Grunau et al 1993) and wettabilities (Blake et al 1995). In the work by Buckles et al (1994), Martys & Chen (1995), and Ferreol & Rothman (1995), realistic sandstone geometries from oil ?elds were used. Very complicated ?ow patterns were observed. The numerical values of the relative permeability as a function of percent saturation of wetting ?uid agree qualitatively with experimental data.
LATTICE BOLTZMANN METHOD
Figure 5 A sandstone sample used in the multiphase LBM simulation (top left). Other panels display two-phase ?uid ?ows through the sandstone at different times (Buckles et al 1994).
Gunstensen & Rothman (1993) studied the linear and nonlinear multicomponent ?ow regimes corresponding to large and small ?ow rates, respectively. They found, for the ?rst time, that the traditional Darcy’s law must be modi?ed in the nonlinear regime because of the capillary effect. Figure 5 shows a typical pore geometry (top left). It is a portion of a sandstone sample obtained from a Mobil offshore oil reservoir. Two-phase ?uid ?ows through the sandstone are displayed at successive times (top right, bottom left, bottom right). The sandstone is transparent. The dark color indicates invading water and the grey color indicates oil. As seen in the plot, the lattice Boltzmann simulation preserves the fundamental phenomena observed in experiments that the water phase forms long ?ngers through the porous medium because of the wettability properties of the water. The LBM is becoming an increasingly popular means of modeling multiphase ?uid ?ows in porous media because of its ability to simulate the exact Navier-Stokes equation in a parallel fashion, to handle complicated geometry, and to simulate surface dynamics and wettability (or contact angles).
CHEN & DOOLEN
LATTICE BOLTZMANN SIMULATION OF PARTICLES IN FLUIDS
To simulate particles suspended in ?uids, an approximate treatment of ?uidparticle interaction must be incorporated. Much of the pioneering work and some interesting applications in this area were carried out by Ladd (1993, 1994a,b, 1997). Important variations and applications were mostly associated with Behrend (1995) and Aidun & Lu (1995). In the model of Ladd, a solid boundary was mapped onto the lattice and a set of boundary nodes, rb , was de?ned in the middle of links, whose interior points represent a particle. A no-slip boundary condition on the moving particle requires the ?uid velocity to have the same speed at the boundary nodes as the particle velocity ub which has the translational part U, and the rotational part b . Assuming that the center position of the particle is R, then, ub = U + b × (rb ? R). The distribution function f i is de?ned for grid points inside and outside the particle. To account for the momentum change when ub was not zero, Ladd proposed adding a term to the distribution function for both sides of the boundary nodes: f i (x) = f i (x) ± B(ei · ub ), (46)
where B is a coef?cient proportional to the mass density of the ?uid and depends on the detailed lattice structure; the + sign applies to boundary nodes at which the particle is moving toward the ?uid and ? for moving away from the ?uid. Equation 46 includes a mass exchange between the ?uid inside and the ?uid outside the particle. Aidun & Lu (1995) modi?ed Equation 46 by adding an extra term that forces mass conservation for the ?uid inside the particle. The details of the boundary rule and variations were studied extensively by Behrend (1995), who improved the ef?ciency of the algorithm used by Ladd while retaining similar accuracy for translational motion, although the rotational motion is less accurate. The approximations used so far to simulate moving boundaries are computationally convenient and ef?cient. In fact, the work for simulating N particles in Ladd’s scheme scales linearly with N , in contrast to ?nite-element schemes, which scale as N 2 (Hu 1996) when the exact incompressible condition is being enforced. The accuracy of the scheme was carefully and extensively studied for creeping ?ows and ?ows at ?nite Reynolds numbers (Ladd 1994b). It compared well with ?nite-difference and ?nite-element methods. The drag coef?cient of a circular particle in a 2-D channel with gravitational ?eld and the settling trajectory compared well with results from traditional numerical methods and the solution of the Stokes equation (Aidun & Lu 1996).
LATTICE BOLTZMANN METHOD
In the LBM, ?uctuations in the distribution functions are ignored. These ?uctuations are important for the study of Brownian motion and other effects. To include the ?uctuations in the ?uid necessary to drive Brownian motion, Dufty & Ernst (1996) and Ladd (1993) added stochastic terms to the distribution function. These distribution function ?uctuations are constrained to conserve mass and momentum, but they contribute a ?uctuating part to the stress tensor of the ?uid. The variance of the stress-tensor ?uctuations is related to the effective temperature and the viscosity of the ?uid via the ?uctuation-dissipation theorem. This method allows—apparently for the ?rst time—treatment of the Brownian short-time regime and the pre-Brownian time regime. Close quantitative agreement is found between experiment and the ?uctuating LBM in the short-time regime (Segre et al 1995). Ladd, using up to 32,000 suspended particles in the ?uctuation LBM, showed that “there is no evidence that the longrange hydrodynamic interactions are screened by changes in the pair correlation function at large distances.” His results leave open the explanation of the experimentally reported absence of the theoretically expected divergence of velocity ?uctuations (Ladd 1997).
SIMULATION OF HEAT TRANSFER AND REACTION-DIFFUSION
The lattice Bhatnagar-Gross-Krook (LBGK) models for thermal ?uids have been developed by several groups. To include a thermal variable, such as temperature, Alexander et al (1993b) used a two-dimensional 13-velocity model on the hexagonal lattice. In this work, the internal energy per unit mass E was de?ned through the second-order moment of the distribution function, 2 ρE = σ,i f σ,i (eσ,i ? u) /2. Here the index σ = 0, 1, and 2, indicates velocity magnitudes. To have a consistent compressible hydrodynamic equation, the collision operator σ,i was chosen to satisfy local energy conserva2 tion: σ,i eσ,i /2 = 0, and the third-order velocity dependent terms, such i 3 as (eσ,i · u) , (eσ,i · u)u 2 and u 3 , had been included in the equilibrium distribution (Equation 21). The transport coef?cients, including viscosity and heat conductivity derived from the Chapman-Enskog expansion, agreed well with numerical simulation. A Couette ?ow with a temperature gradient between two parallel planes was also simulated and the results agreed well with theoretical predictions when the temperature difference was small. Vahala et al (1995) examined this model by studying the effect of 2-D shear velocity turbulence on a steep temperature gradient pro?le. Qian (1993) developed 3-D thermal LBM models based on 21 and 25 velocities. Chen et al (1994) extended the above thermal models by including more speeds based on general square lattices in
CHEN & DOOLEN
D dimensions. They utilized the freedom in the equilibrium distribution to eliminate nonphysical high-order effects (Qian & Orszag 1993). Because the multispeed BGK collision operator constrained the momentum and the thermal modes to relax at the same speed, the Prandtl number was ?xed. The inclusion of additional properties to provide variation of the Prandtl number was studied by Chen et al (1995) and McNamara et al (1995). Two limitations in multispeed LBM thermal models severely restrict their application. First, because only a small set of velocities is used, the variation of temperature is small. Second, all existing LBM models suffer from numerical instability (McNamara et al 1995), owing to the absence of an H-theorem. It seems that the numerical instability is more severe in the thermal LBM models than in the isothermal models. Both problems can be partially resolved if one treats the temperature equation as an active scalar and solves it using the relaxation method by adding an independent distribution function (Bartoloni et al 1993). On the other hand, it is dif?cult for the active scalar approach to incorporate the correct and full dissipation function. Twodimensional Rayleigh-B? nard (RB) convection was simulated using this active e scalar scheme for studying scaling laws (Bartoloni et al 1993) and probability density functions (Massaioli et al 1993) at high Prandtl numbers. Twodimensional free-convective cavity ?ow was also simulated (Eggels & Somers 1995), and the results compared well with benchmark data. Two-D and 3-D Rayleigh-B? nard convections were carefully studied by Shan (1997) using a e passive scalar temperature equation and a Boussinesq approximation. This scalar equation was derived based on the two-component model of Shan & Chen (1993). The calculated critical Rayleigh number for the RB convection agreed well with theoretical predictions. The Nusselt number as a function of Rayleigh number for the 2-D simulation was in good agreement with previous numerical simulation using other methods. The LBM was extended by Dawson et al (1993) to describe a set of reactiondiffusion equations advected by velocities governed by the Navier-Stokes equation. They studied hexagonal patterns and stripe patterns caused by the Turing instability in the Selkov model. The effect of ?uid ?ows on the chemical reaction at solid surfaces was simulated (Chen et al 1995) to study geochemical processes, including dissolution and precipitation on rock surfaces and chimney structures at sea?oor hydrothermal vents. A similar study on the effect of nutrient diffusion and ?ow on coral morphology was carried out by Kaandorp et al (1996). Alvarez-Ramirez et al (1996) used the model of Dawson et al (1993) to study the effective diffusivity of a heterogeneous medium when the inclusion was impermeable or permeable with a different diffusivity. They found that the LBM simulation results compared well with Monte Carlo simulations and
LATTICE BOLTZMANN METHOD
agreed well with the predictions using Maxwell’s equation. Holme & Rothman (1992) carried out a study of viscous ?ngering in two-dimensional Hele-Shaw patterns with P? clet number = 1700. Creeping ?ow in a Hele-Shaw cell was e also simulated by Flekk?y et al (1996) to investigate the inertial effect at very small Reynolds number, which is sensitive to effects of hydrodynamic irreversibility. Anomalous diffusion of LBM ?uids in fractal media and Taylor hydrodynamic dispersion in a two-dimensional channel were simulated by Cali et al (1992). They found that the diffusion scaling exponents agreed well with the exact solution for a Sierpinski gasket. The measured longitudinal diffusion coef?cient as a function of ?ow speed in 2-D Taylor hydrodynamic dispersion agreed well with an asymptotic prediction. This research demonstrates that the LBM is capable of simulating the diffusion process in irregular geometries. The effects of turbulent mixing behind a wake on Turing patterns were studied by Weimar & Boon (1996). They found that the mixing caused by turbulence increased eddy diffusivity and destroyed the pattern formation. Qian et al (1995) simulated front dynamics and scaling properties for the irreversible reaction: A + B → C.
This review is intended to be an overview of the subject of the lattice Boltzmann method and its applications in ?uid mechanics. The LBM is so diverse and interdisciplinary that it is not possible to include all interesting topics. We refer readers to related review articles written by Benzi et al (1992), Rothman & Zaleski (1994), Chen et al (1995), and Qian et al (1995). We have introduced the fundamentals of the lattice Boltzmann method, including the lattice Boltzmann equation and its relation to the macroscopic Navier-Stokes equation. Some LBM methods, including multiphase LBM models and ?uid-particle interaction models, have been discussed in detail. We have demonstrated that simulation results from the lattice Boltzmann methods are in good quantitative agreement with experimental results and results from other numerical methods. We wish to emphasize the kinetic nature of the LBM. Because the LBM is a mesoscopic and dynamic description of the physics of ?uids, it can model problems wherein both macroscopic hydrodynamics and microscopic statistics are important. The LBM can be considered to be an ef?cient numerical method for computational ?uid dynamics. It is also a powerful tool for modeling new physical phenomena that are not yet easily described by macroscopic equations. From a computational point of view, the lattice Boltzmann equation is hyperbolic and can be solved locally, explicitly, and ef?ciently on parallel computers.
CHEN & DOOLEN
Not only is the scheme computationally comparable with traditional numerical methods, but it is also easy to program and to include new physics because of the simplicity of the form of the LBM equations. The lattice Boltzmann method is still undergoing development. Many models, including simulation of granular ?ows (Flekk?y & Herrmann 1993, Tan et al 1995), viscoelastic ?ows (Aharonov & Rothman 1993), magnetohydrodynamics (Martinez et al 1994), and microemulsions (Boghosian et al 1996) were recently proposed. Although innovative and promising, these existing LBM methods, including multicomponent LBM models, require additional benchmarking and veri?cation. The current lattice models for multiphase and reacting systems are most suitable for isothermal problems. The development of a reliable LBM for thermal systems will allow the simulation of heat transfer and surface phenomena simultaneously (Kato 1997). This would open many new areas of application. ACKNOWLEDGMENTS We are grateful to N Cao, H Chen, B Hasslacher, X He, S Hou, S Jin, AJC Ladd, T Lookman, L Luo, D Martinez, W Matthaeus, B Nadiga, Y Qian, X Shan, S Succi, and MR Yeomans for useful discussions.
Visit the Annual Reviews home page at http://www.AnnualReviews.org.
Literature Cited Abe T. 1997. Derivation of the lattice Boltzmann method by means of the discrete ordinate method for the Boltzmann equation. J Comp. Phys. 131:241–46 Aharonov E, Rothman DH. 1993. NonNewtonian ?ow (through porous media): a lattice Boltzmann method. Geophys. Res. Lett. 20:679–82 Aidun CK, Lu Y-N. 1995. Lattice Boltzmann simulation of solid particles suspended in ?uid. J. Stat. Phys. 81:49–61 Alexander FJ, Chen S, Grunau DW. 1993a. Hydrodynamic spinodal decomposition: growth kinetics and scaling. Phys. Rev. B 48:634– 37 Alexander FJ, Chen S, Sterling JD. 1993b. Lattice Boltzmann thermohydrodynamics. Phys. Rev. E 47:R2249–52 Alvarez-Ramirez J, Nieves-Mendoza S, Gonzalez-Trejo J. 1996. Calculation of the effective diffusivity of heterogeneous media using the lattice Boltzmann method. Phys. Rev. E 53:2298–2303 Amati G, Succi S, Benzi R. 1997. Turbulent channel ?ow simulation using a coarsegrained extension of the lattice Boltzmann method. Fluid Dyn. Res. 19:289–302 Ancona MG. 1994. Fully-Lagrangian and lattice-Boltzmann methods for solving systems of conservation equations. J. Comp. Phys. 115:107–20 Bartoloni A, Battista C, Cabasino S, Paolucci PS, Pech J, et al. 1993. LBE simulations of Rayleigh-Bernard convection on the APE100 parallel processor. Int. J. Mod. Phys. C 4:993– 1006 Behrend O. 1995. Solid-?uid boundaries in particle suspension simulations via the lattice Boltzmann method. Phys. Rev. E 52:1164– 75 Benzi R, Succi S. 1990. Two-dimensional turbulence with the lattice Boltzmann equation. J. Phys. A 23:L1–5 Benzi R, Succi S, Vergassola M. 1992. The lattice Boltzmann-equation—theory and applications. Phys. Rep. 222:145–97 Benzi R, Struglia MV, Tripiccione R. 1996. Extended self-similarity in numerical simula-
LATTICE BOLTZMANN METHOD
tions of three-dimensional anisotropic turbulence. Phys. Rev. E 53:R5565–68 Bhatnagar PL, Gross EP, Krook M. 1954. A model for collision processes in gases. I: small amplitude processes in charged and neutral one-component system. Phys. Rev. 94:511–525 Blake TD, Deconinck J, Dortona U. 1995. Models of wetting: immiscible lattice Boltzmann automata versus molecular kinetic theory. Langmuir. 11:4588–92 Boghosian BM, Coveney PV, Emerton AN. 1996. A lattice-gas model for microemulsions. Proc. Roy. Soc. London. Ser. A 452:1221–50 Brackbill JU, Kothe DB, Zemach C. 1992. A continuum method for modeling surface tension. J. Comp. Phys. 100:335–54 Broadwell JE. 1964. Study of rare?ed shear ?ow by the discrete velocity method. J. Fluid Mech. 19:401–14 Buckles J, Hazlett R, Chen S, Eggert KG, Grunau DW, Soll WE. 1994. Flow through porous media using lattice Boltzmann method. Los Alamos Sci. 22:112–21 Cali A, Succi S, Cancelliere A, Benzi R, Gramignani M. 1992. Diffusion and hydrodynamic dispersion with the lattice Boltzmann method. Phys. Rev. A 45:5771–74 Cao N, Chen S, Jin S, Martinez D. 1997. Physical symmetry and lattice symmetry in lattice Boltzmann method. Phys. Rev. E 55:R21– 24 Chang YC, Hou TY, Merriman B, Osher S. 1996. A level set formulation of Eulerian interface capturing methods for incompressible ?uid ?ows. J. Comp. Phys. 124:449–64 Chen H. 1993. Discrete Boltzmann systems and ?uid ?ows. Comp. Phys. 7:632–37 Chen H, Chen S, Matthaeus WH. 1992. Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method. Phys. Rev. A. 45:R5339–42 Chen S, Chen HD, Martinez D, Matthaeus W. 1991. Lattice Boltzmann model for simulation of magnetohydrodynamics. Phys Rev. Lett. 67:3776–79 Chen S, Dawson SP, Doolen GD, Janecky DR, Lawniczak A. 1995. Lattice methods and their applications to reacting systems. Comp. Chem. Eng. 19:617–46 Chen S, Lookman T. 1995. Growth kinetics in multicomponent ?uids. J. Stat. Phys. 81:223– 35 Chen S, Martinez D, Mei R. 1996. On boundary conditions in lattice Boltzmann methods. Phys. Fluids 8:2527–36 Chen S, Wang Z, Shan XW, Doolen GD. 1992. Lattice Boltzmann computational ?uid dynamics in three dimensions. J. Stat. Phys. 68:379–400
Chen Y, Ohashi H, Akiyama M. 1994. Thermal lattice Bhatnagar-Cross-Krook model without nonlinear deviations in macrodynamic equations. Phys. Rev. E 50:2776–83 Chen Y, Ohashi H, Akiyama M. 1995. Prandtl number of lattice Bhatnagar-Gross-Krook ?uid. Phys. Fluids 7:2280–82 Chu CK. 1965. Kinetic-Theoretic description of the formation of a shock wave. Phys. Fluids 8:12–21 Cieplak M. 1995. Rupture and coalescence in 2dimensional cellular-automata ?uids. Phys. Rev. E 51:4353–61 Cornubert R, d’Humi` res D, Levermore D. e 1991. A Knudsen layer theory for lattice gases. Physica D 47:241–59 Dawson SP, Chen S, Doolen GD. 1993. Lattice Boltzmann computations for reactiondiffusion equations. J. Chem. Phys. 98:1514– 23 d’Humi` res D, Lallemand P, Frisch U. 1986. e Lattice gas model for 3D hydrodynamics. Europhys. Lett. 2:291–97 D’ortona U, Salin D, Cieplak M, Banavar JR. 1994. Interfacial phenomena in Boltzmann cellular automata. Europhys. Lett. 28:317–22 Dufty JW, Ernst MH. 1996. Lattice BoltzmannLangevin equation. Fields Inst. Comm. 6:99– 107 E W-N, Liu J-G. 1996. Essentially compact scheme for unsteady viscous incompressible ?ows. J. Comp. Phys. 126:122–38 Eggels JGM. 1996. Direct and large-eddy simulation of turbulent ?uid ?ow using the latticeBoltzmann scheme. Int. J. Heat Fluid Flow 17:307–23 Eggels JGM, Somers JA. 1995. Numerical simulation of free convective ?ow using the lattice-Boltzmann scheme. J. Heat Fluid Flow 16:357–64 Elton BH, Levermore CD, Rodrigue H. 1995. Covergence of convective-diffusive lattice Boltzmann methods. SIAM J. Sci. Comp. 32:1327–54 Ferreol B, Rothman DH. 1995. LatticeBoltzmann simulations of ?ow-through Fontainebleau sandstone. Transp. Por. Med. 20:3–20 Flekk?y EG, Herrmann HJ. 1993. Lattice Boltzmann models for complex ?uids. Physica A. 199:1–11 Flekk?y EG, Rage T, Oxaal U, Feder J. 1996. Hydrodynamic irreversibility in creeping ?ow. Phys. Rev. Lett. 77:4170–73 Frisch U, Hasslacher B, Pomeau Y. 1986. Lattice-gas automata for the Navier-Stokes equations. Phys. Rev. Lett. 56:1505–8 Frisch U, d’Humi? res D, Hasslacher B, Lallee mand P, Pomeau Y, Rivet J-P. 1987. Lattice gas hydrodynamics in two and three dimensions. Complex Syst. 1:649–707
CHEN & DOOLEN
Hou S. 1995. Lattice Boltzmann method for incompressible viscous ?ow. PhD thesis. Kansas State Univ., Manhattan, Kansas Hou S, Zou Q, Chen S, Doolen G, Cogley AC. 1995. Simulation of cavity ?ow by the lattice Boltzmann method. J. Comp. Phys. 118:329– 47 Hou S, Sterling J, Chen S, Doolen GD. 1996. A lattice Boltzmann subgrid model for high Reynolds number ?ows. Fields Inst. Comm. 6:151–66 Hu HH. 1996. Direct simulation of ?ows of solid-liquid mixtures. Int. J. Multiphase ?ow. 22:335–52 Inamuro T, Sturtevant B. 1990. Numerical study of discrete-velocity gases. Phys. Fluids 2:2196–2203 Inamuro T, Yoshino M, Ogino F. 1995. A nonslip boundary condition for lattice Boltzmann simulations. Phys. Fluids 7:2928–30 Jin S, Katsoulakis M. Relaxation approximations to front propagation. J. Diff. Eq. In press Jin S, Xin ZP. 1995. The relaxation schemes for systems of conservation laws in arbitrary space dimensions. Comm. Pure Appl. Math. 48:235–76. Kaandorp JA, Lowe CP, Frenkel D, Sloot PMA. 1996. Effect of nutrient diffusion and ?ow on coral morphology. Phys. Rev. Lett. 77:2328– 31 Kadanoff L. 1986. On two levels. Phys. Today 39:7–9 Kato T, Kono K, Seta K, Martinez D, Chen S. 1997. Boiling two-phase ?ow by the lattice Boltzmann method Int. J. Mod. Phys. In press Kim J, Moin P, Moser R. 1987. Turbulence statistics in fully-developed channel ?ow at low Reynolds number. J. Fluid Mech. 177:133–66 Koelman JMVA. 1991. A simple lattice Boltzmann scheme for Navier-Stokes ?uid ?ow. Europhys. Lett. 15:603–7 Koplik J, Lasseter T. 1985. One- and two-phase ?ow in network models of porous media. Chem. Eng. Comm. 26:285–95 Kraichnan RH. 1967. Inertial ranges in twodimensional turbulence. Phys. Fluids 10: 1417–23 Ladd AJC. 1993. Short-time motion of colloidal particles: numerical simulation via a ?uctuating lattice-Boltzmann equation Phys. Rev. Lett. 70:1339–42 Ladd AJC. 1994a. Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part. 1. theoretical foundation. J. Fluid Mech. 271:285–309 Ladd AJC. 1994b. Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 2: Numerical results. J. Fluid Mech. 271:311–39 Ladd AJC. 1997. Sedimentation of homoge-
Ginzbourg I, Adler PM. 1994. Boundary ?ow condition analysis for the three-dimensional lattice Boltzmann model. J. Phys. II 4:191– 214 Glimm JE, Issacson E, Marchesin D, McBryan O. 1981. Front tracking for hyperbolic systems. Adv. Appl. Math. 2:91–119 Gonnella G, Orlandini E, Yeomans JM. 1997. Spinodal decomposition to a lamellar phase: effects of hydrodynamic ?ow. Phys. Rev. Lett. 78:1695–98 Grunau D, Chen S, Eggert K. 1993. A lattice Boltzmann model for multiphase ?uid ?ows. Phys. Fluids A 5:2557–62 Gunstensen AK, Rothman DH. 1993. Lattice Boltzmann studies of immiscible two phase ?ow through porous media J. Geophys. Res. 98:6431–41 Gunstensen AK, Rothman DH, Zaleski S, Zanetti G. 1991. Lattice Boltzmann model of immiscible ?uids. Phys. Rev. A. 43:4320– 27 Hardy J, de Pazzis O, Pomeau Y. 1976. Molecular dynamics of a classical lattice gas: transport properties and time correlation functions. Phys. Rev. A 13:1949–61 Hayot F, Wagner L. 1996. A non-local modi?cation of a lattice Boltzmann model. Europhys. Lett. 33:435–40 He X, Doolen GD. 1997. Lattice Boltzmann method on curvilinear coordinates system: vortex shedding behind a circular cylinder. Phys. Rev. E. In press He X, Luo L-S. 1997. A priori derivation of the lattice Boltzmann equation. Phys. Rev. E. 55:6333–36 He X, Luo L-S, Dembo M. 1996. Some progress in lattice Boltzmann method. Part I nonuniform mesh grids. J. Comp. Phys. 129:357– 63 He X, Zou Q, Luo L-S, Dembo M. 1997. Analytic solutions of simple ?ow and analysis of non-slip boundary conditions for the lattice Boltzmann BGK model. J. Stat. Phys. 87:115–36 Heijs AWJ, Lowe CP. 1995. Numerical evaluation of the permeability and the Kozeny constant for two types of porous media. Phys. Rev. E 51:4346–52 Higuera FJ, Jim? nez J. 1989. Boltzmann ape proach to lattice gas simulations. Europhys. Lett. 9:663–68 Higuera FJ, Succi S. 1989. Simulating the ?ow around a circular cylinder with a lattice Boltzmann equation. Europhys. Lett. 8:517–21 Higuera FJ, Succi S, Benzi R. 1989. Lattice gas dynamics with enhanced collisions. Europhys. Lett. 9:345–49 Holme R, Rothman DH. 1992. Lattice-gas and lattice-Boltzmann models of miscible ?uids. J. Stat. Phys. 68:409–29
LATTICE BOLTZMANN METHOD
neous suspensions of non-Brownian spheres Phys. Fluids 9:491–99 Laudau LD, Lifshitz EM. 1959. Fluid Dyn. New York: Pergamon Lavall? e P, Boon JP, Noullez A. 1991. Bounde aries in lattice gas ?ows. Physica D 47:233– 40 Luo L. 1997. Symmetry breaking of ?ow in 2D symmetric channels: simulations by lattice Boltzmann method. Int. J. Mod. Phys. C. In press Maier RS, Bernard RS, Grunau DW. 1996. Boundary conditions for the lattice Boltzmann method. Phys. Fluids 8:1788–1801 Majda A. 1984. Compressible Fluid and Systems of Conservation Laws in Several Space Variables. New York: Springer-Verlag Martys NS, Chen H. 1996. Simulation of multicomponent ?uids in complex threedimensional geometries by the lattice Boltzmann method. Phys. Rev. E 53:743–50 Martinez DO, Chen S, Matthaeus WH. 1994. Lattice Boltzmann magnetohydrodynamics. Phys. Plasmas 1:1850–67 Martinez DO, Matthaeus WH, Chen S, Montgomery DC. 1994. Comparison of spectral method and lattice Boltzmann simulations of two-dimensional hydrodynamics. Phys. Fluids 6:1285–98 Massaioli F, Benzi R, Succi S. 1993. Exponential tails in two-dimensional RayleighB? nard convection. Europhys. Lett. 21:305– e 10 McNamara GR, Garcia AL, Alder BJ. 1995. Stabilization of thermal lattice Boltzmann models. J. Stat. Phys. 81:395–408 McNamara GR, Zanetti G. 1988. Use of the Boltzmann equation to simulate lattice-gas automata. Phys. Rev. Lett. 61:2332–35 Miller W. 1995. Flow in the driven cavity calculated by the lattice Boltzmann method. Phys. Rev. E 51:3659–69 Nadiga BT. 1995. An Euler solver based on locally adaptive discrete velocities. J. Stat. Phys. 81:129–46 Nadiga BT, Pullin DI. 1994. A method for near-equilibrium discrete-velocity gas ?ows. J. Comp. Phys. 112:162–72 Nadiga BT, Zaleski S. 1996. Investigations of a two-phase ?uid model. Eur. J. Mech. B/Fluids 15:885–96 Nannelli F, Succi S. 1992. The lattice Boltzmann-equation on irregular lattices. J. Stat. Phys. 68:401–7 Noble DR, Chen S, Georgiadis JG, Buckius RO. 1995. A consistent hydrodynamic boundary condition for the lattice Boltzmann method. Phys. Fluids 7:203–9 Noble DR, Georgiadis JG, Buckius RO. 1996. Comparison of accuracy and performance for lattice Boltzmann and ?nite difference simu-
lations of steady viscous ?ow. Int. J. Numer. Meth. Fluids 23:1–18 Orszag SA, Yakhot V. 1986. Reynolds number scaling of cellular automaton hydrodynamics. Phys. Rev. Lett. 56:1691–93 Osborn WR, Orlandini E, Swift MR, Yeomans JM, Banavar JR. 1995. Lattice Boltzmann study of hydrodynamic spinodal decomposition. Phys. Rev. Lett. 75:4031–34 Prasad AK, Koseff JR. 1989. Reynolds-number and end-wall effects on a lid-driven cavity ?ow. Phys. Fluids A 1:208–18 Qian YH. 1990. Lattice gas and lattice kinetic theory applied to the Navier-Stokes equations. PhD thesis. Universit? Pierre et Marie e Curie, Paris Qian YH. 1993. Simulating thermohydrodynamics with lattice BGK models. J. Sci. Comp. 8:231–41 Qian YH, d’Humi` res D, Lallemand P. 1992. e Lattice BGK models for Navier-Stokes equation. Europhys. Lett. 17:479–84 Qian YH, Orszag SA. 1993. Lattice BGK models for the Navier-Stokes equation: nonlinear deviation in compressible regimes. Europhys. Lett. 21:255–59 Qian YH, Succi S, Massaioli F, Orszag SA. 1996. A benchmark for lattice BGK model: ?ow over a backward-facing step. Fields Inst. Comm. 6:207–15 Qian YH, Succi S, Orszag SA. 1995. Recent advances in lattice Boltzmann computing. Annu. Rev. Comp. Phys. 3:195–242 Rothman DH. 1988. Cellular-automaton ?uids: a model for ?ow in porous media. Geophys. 53:509–18 Rothman DH, Keller JM. 1988. Immiscible cellular-automaton ?uids. J. Stat. Phys. 52:1119–27 Rothman DH, Zaleski S. 1994. Lattice-gas models of phase separation: interfaces, phase transitions and multiphase ?ow. Rev. Mod. Phys. 66:1417–79 Rybka RB, Cieplak M, Salin D. 1995. Boltzmann cellular-automata studies of the spinodal decomposition. Physica A 222:105–18 Reider MB, Sterling JD. 1995. Accuracy of discrete-velocity BGK models for the simulation of the incompressible Navier-Stokes equations. Comput. Fluids 24:459–67 Schreiber R, Keller HB. 1983. Driven cavity ?ows by ef?cient numerical techniques. J. Comp. Phys. 49:310–33 Schwarz LM, Martys N, Bentz DP, Garboczi EJ, Torquato S. 1993. Cross-property relations and permeability estimation in model porous media. Phys. Rev. E 48:4584–91 Segre PN, Behrend OP, Pusey PN. 1995. Shorttime Brownian motion in colloidal suspensions: experiment and simulation. Phys. Rev. E 52:5070–83
CHEN & DOOLEN
dimensional ?ows in complex geometries with the lattice Boltzmann method. Europhys. Lett. 10:433–38 Succi S, Nannelli F. 1994. The ?nite volume formulation of the lattice Boltzmann equation. Transp. Theory Stat. Phys. 23:163–71 Tan M-L, Qian YH, Goldhirsch I, Orszag SA. 1995. Lattice-BGK approach to simulating granular ?ows. J. Stat. Phys. 81:87–103 Trevi? o C, Higuera F. 1994. Lattice Boltzmann n and spectral simulations of non-linear stability of Kolmogorov ?ows. Rev. Mex. Fis. 40:878–90 Vahala G, Pavlo P, Vahala L, Chen H. 1995. Effect of velocity shear on a strong temperature gradient—a lattice Boltzmann approach. Phys. Lett. A 202:376–82 Vangenabeek O, Rothman DH. 1996. Macroscopic manifestations of microscopic ?ows through porous-media phenomenology from simulation. Annu. Rev. Earth Planet. Sci. 24:63–87 Vanka SP. 1986. Block-implicit multigrid solution of Navier-Stokes equations in primitive variables. J. Comp. Phys. 65:138–58 Wagner L. 1994. Pressure in lattice Boltzmann simulations of ?ow around a cylinder. Phys. Fluids 6:3516–18 Weimar JR, Boon JP. 1996. Nonlinear reactions advected by a ?ow. Physica A 224:207– 15 Wolfram S. 1986. Cellular automaton ?uids. 1: Basic theory. J. Stat. Phys. 45:471–526 Wu Y, Alexander FJ, Lookman T, Chen S. 1995. Effects of hydrodynamics on phase transition kinetics in two-dimensional binary ?uids. Phys. Rev. Lett. 74:3852–55 Xu K, Prendergast KH. 1994. Numerical Navier-Stokes solutions from gas kinetic theory. J. Comp. Phys. 114:9–17 Ziegler DP. 1993. Boundary conditions for lattice Boltzmann simulations. J. Stat. Phys. 71:1171–77 Zou Q, He X. 1997. On pressure and velocity boundary conditions for the lattice Boltzmann BGK model. Phys. Fluids. 9:1591–98 Zou Q, Hou S, Doolen GD. 1995. Analytical solutions of the lattice Boltzmann BGK model. J. Stat. Phys. 81:319–34
Shan X. 1997. Simulation of Rayleigh-B? nard e convection using a lattice Boltzmann method. Phys. Rev. E 55:2780–88 Shan X, Chen H. 1993. Lattice Boltzmann model for simulating ?ows with multiple phases and components. Phys. Rev. E 47:1815–19 Shan X, Chen H. 1994. Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation. Phys. Rev. E 49:2941–48 Shan X, Doolen G. 1995. Multicomponent lattice-Boltzmann model with interparticle interaction. J. Stat. Phys. 81:379–93 Skordos PA. 1993. Initial and boundary conditions for the lattice Boltzmann method. Phys. Rev. E 48:4823–42 Soll W, Chen S, Eggert K, Grunau D, Janecky D. 1994. Applications of the latticeBoltzmann/lattice gas techniques to multi?uid ?ow in porous media. In Computational Methods in Water Resources X, ed. A Peters, pp. 991–99. Dordrecht: Academic Somers JA. 1993. Direct simulation of ?uid ?ow with cellular automata and the latticeBoltzmann equation. Applied Sci. Res. 51:127–33 Spaid MAA, Phelan FR Jr. 1997. Lattice Boltzmann methods for modeling microscale ?ow in ?brous porous media. Phys. Fluids. 9:2468–74 Sterling JD, Chen S. 1996. Stability analysis of lattice Boltzmann methods. J. Comp. Phys. 123:196–206 Swift MR, Osborn WR, Yeomans JM. 1995. Lattice Boltzmann simulation of nonideal ?uids. Phys. Rev. Lett. 75:830–33 Swift MR, Orlandini SE, Osborn WR, Yeomans JM. 1996. Lattice Boltzmann simulations of liquid-gas and binary-?uid systems. Phys. Rev. E 54:5041–52 Succi S, Amati G, Benzi R. 1995. Challenges in lattice Boltzmann computing. J. Stat. Phys. 81:5–16 Succi S, Benzi R, Higuera F. 1991. The latticeBoltzmann equation—a new tool for computational ?uid dynamics. Physica D 47:219– 30 Succi S, Foti E, Higuera F. 1989. Three-