PHYSICAL REVIEW E 83, 026702 (2011)

Lattice-Boltzmann method combined with smoothed-profile method for particulate suspensions
Saeed Jafari,1,2,* Ryoichi Yamamoto,3 and Mohamad Rahnama2
1

Department of Petroleum engineering, Shahid Bahonar University of Kerman, Kerman, Iran2
Department of Mechanical Engineering, Shahid Bahonar University of Kerman, Kerman, Iran3
Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan(Received 21 July 2010; revised manuscript received 8 December 2010; published 9 February 2011)
We developed a simulation scheme based on the coupling of the lattice-Boltzmann method with the smoothedprofile method (SPM) to predict the dynamic behavior of colloidal dispersions. The SPM provides a coupling
scheme between continuum fluid dynamics and rigid-body dynamics through a smoothed profile of the fluidparticle interface. In this approach, the flow is computed on fixed Eulerian grids which are also used for theparticles. Owing to the use of the same grids for simulation of fluid flow and particles, this method is highly
efficient. Furthermore, an external boundary is used to impose the no-slip boundary condition at the fluid-particle
interface. In addition, the operations in the present method are local; it can be easily programmed for parallelmachines. The methodology is validated by comparing with previously published data.
DOI: 10.1103/PhysRevE.83.026702

PACS number(s): 02.70.−c, 82.70.Dd

I. INTRODUCTION

use of a finer lattice grid with more nodes at the boundaryor higher-order bounce back based on interpolations [21,22]
can reduce this fluctuation, these approaches will increase thecomputational cost. Also, a continuous bounce-back method
was developed to incorporate solid-fluid boundary conditions
on length scales smaller than the grid spacing [23,24].
Another way to simulate fluid-particle interactions is touse a combination of the immersed-boundary method (IBM)and the LBM. In this kind of method [18–20], an immersedboundary body force is introduced to capture the particle
motion. In one such method, the immersed-boundary force
density is worked out via a linear spring relation due todisplacement between the boundary point and the computational point at the solid boundary [18]. A user-defined springparameter may have a significant effect on the computational
efficiency and accuracy. Another method computes the force
density at the particle boundary by setting the velocity on theimmersed Lagrangian boundary points equal to the velocity ofthe particle [19]. In these methods a discrete δ function is used
to transform the force density from the boundary points to theflow elements. Recently, two overlapping grid systems wereused [20], a regular Eulerian grid for the fluid domain and a
Lagrangian for the particle domain, and the no-slip boundary
condition at the solid interface was applied by adding a force
density to the fluid domain to force the difference between
the fluid and the solid velocity at the boundary nodes tobe zero; the counterforce acting on the solid particles wasused to update the position and velocity of particles basedon Newtonian dynamics equations. As in previous methods,
a discrete δ function was used to transform the force densityfrom the boundary points to the Eulerian points.
To avoid complicated interpolations and transformationsof the force density, a different direct forcing method is
introduced based on the smoothed-profile method (SPM) [25].
It uses fixed Eulerian grids for the host fluid. It represents
the particles by certain smooth body forces in the NavierStokes equations instead of treating the particles as boundary
conditions to the fluid. The SPM solves a single set of
fluid dynamics equations in the entire domain including theparticle volumes, without any internal boundary conditions.
A smoothly spreading interface layer is used to represent the

The flow of a dense suspension of colloidal particles is
encountered in many industrial processing units. Proper investigation of fluid-particle interactions is crucial for obtaining
good performance. Numerical simulations can aid in extending
our understanding to more complex colloidal suspensions.
In the last two decades, the lattice-Boltzmann method(LBM) [1–4] has attracted much attention as a promising alternative for simulation of fluid flows with complex physics suchas multiphase flow and multicomponent flow [5–9], simulationof particulate suspensions and colloid hydrodynamics [10–13],
and adding thermal fluctuations to fluids [13–16]. The LBM isa method based on the solution of the Boltzmann equation on alattice with a discrete velocity field. It was shown that the basic
conservation equations of fluid flow (Navier-Stokes equations)can be recovered from the Boltzmann equation [17]. Solutionof the Boltzmann equation provides a velocity distributionfunction from which macroscopic fluid properties, such as
density, velocity, and pressure, can be obtained. Some advantages of using the LBM in computing fluid flow problems ascompared to computational fluid dynamics (CFD) are the lackof a convective term in the Boltzmann equation and the simple
pressure computation using an equation of state [2]. Moreover,the streaming-and-collision computational procedure of the
LBM, which is a local operation in computation, makes it anexcellent candidate for parallel computing.
The LBM has been used to simulate particles suspended inhost fluids [15,16,18–20]. In some of these methods [15,16],
the bounce-back rule is employed to treat the no-slip boundary
condition on the fluid-particle interface, and the particle
interface is represented by the boundary nodes, which are a
set of midpoints of the links between two fixed grids. In theseapproaches, the transfer of distribution functions is used toaccount for the momentum transfer at the fluid-particle interface. However, with these methods, the solid boundary will not
move continuously and smoothly in space; instead it will jumpfrom one midpoint to another, causing fluctuation. Although*jafari@alumni.iut.ac.ir

1539-3755/2011/83(2)/026702(7)026702-1©2011 American Physical Society

SAEED JAFARI, RYOICHI YAMAMOTO, AND MOHAMAD RAHNAMAPHYSICAL REVIEW E 83, 026702 (2011)

particle boundaries to give a transition from the rigid-body
motion to the fluid motion. By using this simple modification,
regular Cartesian coordinates can be employed for manyparticle systems with any arbitrary particle shape, rather thanboundary-fitted coordinates. The solid-fluid interface has afinite volume supported by multiple grid points. Thus, a roundparticle shape can be treated in a fixed Cartesian grid without
difficulty. We do not need to consider a lot of points on theparticle surface [19] or inside the particle [20]. Furthermore,
the computational demands for this method include sensitivity
to the number of grid points (volume of the total system);
however, it is insensitive to the number of particles [25]. The
common feature of using Cartesian mesh in the LBM and SPM
motivates us to combine them into one efficient method. The
numerical method, including the explanation of fluid-particle
interaction body force, is described in the following section.
Several examples are solved in Sec. III to show the accuracy
and the validation of the present model. Concluding remarks
are given in the last section.
FIG. 1. An example of a smoothed profile (solid line).
II. NUMERICAL METHOD
A. Fluid flow simulation

The two-dimensional LBM is used for fluid flow simulation.
In this approach, the kinetic evolution of the lattice-Boltzmann
equation with the Bhatnagar-Gross-Krook (BGK) collision
approximation is expressed as1fα (x + cα t,t + t) = fα (x,t) − fα (x,t) − fαeq (x,t) ,τ(1)
where τ is the dimensionless relaxation time, and fα is the
distribution function along the α direction. Equation (1) is adiscrete finite-difference equation and is second-order accuratewith respect to space and time. For the D2Q9 lattice model, that
is, a model of a two-dimensional problem with nine velocities,
the discrete velocities are given by c0 = (0,0); cα = (±1, 0)cand (0, ± 1)c for α = 1−4 and cα = (±1, ± 1)c for α = 5−8.
The streaming speed c is defined as x/t, where x andt are the lattice spacing and time step, respectively.
eq
The equilibrium distribution function fα for the incompressible D2Q9 model is

Navier-Stokes equation applied to the LBM to implement the
no-slip boundary condition for simulation of particulate flow.
In the smoothed-profile method, the surface of the colloid istreated not as a sharp interface having no thickness, but rather
an interface is introduced having a width comparable to thegrid spacing. The colloid density profile is defined so as tochange smoothly within this finite-thickness interface.
In order to determine regions in which a colloidal particle
exists, the density field is introduced:
ϕ(x,t) =ϕi (x,t).

Here, φi ∈ [0,1] is the density profile of the ith particle.
In the fluid regions, φi = 0, and in particle regions, φi = 1.
Between these values its value changes continuously in theinterface region of thickness ξ (Fig. 1). The functional formof φ is defined arbitrarily; specific examples are contained inRef. [25]. In the present study, the following function is used:
ϕi (x) =⎧s(R − |x − Ri |),x < −ξ/2,
⎨ 0,
+
1
,
|x| < ξ/2,
s(x) = 12 sin πxξ⎩1,x > ξ/2,

A method is developed here that is based on an external boundary force [25] at the solid-fluid interface in the(4)where R is the radius of the ith particle and Ri is the center
position of the ith particle. The density field and velocity fieldof the particles at t = tn are defined using Eq. (4):ϕ(x,tn ) =NPs(R − |x − Ri (tn )|),

(5)

i=1

ϕ (x,tn )uP (x,tn ) =NPs(R − |x − Ri (tn )|)i=1

× [Vi (tn ) +B. Fluid-solid interaction

(3)

i=1

fαeq = ωα ρ[1 + 3(cα · u)/c2 + 9(cα · u)2 /2c4 − 3u2 /2c2 ],(2)where ω0 = 4/9, ωα = 1/9 for i = 1,2,3,4 and ωα = 1/36 for
i = 5,6,7,8, and u is the local fluid velocity. The macroscopic
variables for the fluid mass density and fluid momentum are
defined by ρ = 8α=1 fα and ρ u = 8α=1 cα fα . By using
the Chapman-Enskog procedure, the incompressible NavierStokes equations can be obtained from Eq. (1) with a shear
viscosityν = cs2 t(τ − 1/2), where cs is the speed of soundand cs2 = c2 /3. More details regarding the derivation of
the Navier-Stokes equations from the LBM can be foundelsewhere [1,2].

NPi (tn ) × {x − Ri (tn )}](6)

where φuP is the velocity field of the particles. In the SPM,
time evolution is computed for {Ri ,Vi ,i } (i = 1, . . . , NP ),which are the position, translational velocity, and angular

026702-2

LATTICE-BOLTZMANN METHOD COMBINED WITH . . .

PHYSICAL REVIEW E 83, 026702 (2011)
-UW/2

UW/2

FIG. 3. A neutrally buoyant cylinder in simple shear flow.

By adding an additional term to the collision function, thefluid-solid interaction force from (8) is included in the latticeBoltzmann equation asfα (x + cα t,tn + t) = fα (x,tn ) −+1
fα (x,tn ) − fαeq (x,tn )τωα t[fH (x,tn ) · cα ].
cs2

(9)

As the particle nodes coincide with the fluid nodes, theparticle-fluid interaction force is completely local and thisforce is zero outside the particles’ domain.
The remaining steps are to update the velocity and angularvelocity of particles. The force FiH and torque NiH exerted by
the host fluid are obtained from the momentum conservation
law. From momentum exchange during the time interval t,
we obtain
tn +ttn

FH
i ds =
tn +tFIG. 2. Streamlines of the flow around a circular cylinder at
Re = 20 and 40.

velocity of the ith particle, respectively, and NP is the numberof particles. Let the position, velocity, and angular velocity of
the colloidal particles {Ri ,Vi ,i } (i = 1, . . . , NP ) be given at
the nth discretization time tn . From these variables, φ and φuP
are constructed.
The initial velocity in the fluid domain and the particles’
positions and velocities are known. The fluid nodes covered by
the particle must have the same velocity as the solid particle. A
body force is introduced over the particle inner domain to force
the fictitious fluid to satisfy the rigid-body motion constraint;
this force is zero outside the particle domain.
The fluid-solid interaction force φfP acting on the solid
particle nodes is given bytn +tϕ (x,tn )fP (x,tn )ds = ϕ(x,tn )[uP (x,tn ) − u(x,tn )].

tn

(7)

tn

fH (x,tn ) = −ϕ (x,tn )fP (x,tn )ds.
tn

(8)ρϕi (x,tn ) [u(x,tn ) − uP (x,tn )] d∀P ,(10)NHi ds =∀P[x − Ri (tn )] × ρϕi (x,tn )× [u(x,tn ) − uP (x,tn )]d∀P .

(11)Here, in order to compute FiH and NiH , the volume integral
over the particle volume, ∀P , is used, which is easier tocompute than the surface integral [26,27].
The velocity and angular velocity of the colloidal particles
are updated using this fluid drag and other forces:
Vi (tn + t) = Vi (tn ) + MP−1
i (tn + t) = i (tn ) + I−1
P

tn +t

c
ext
FH
ds,
i + Fi + Fi

tn
tn +t
tn

(12)
ext
NH
i + Nids. (13)

In Eq. (12) Fic and Fiext are the collisional force on the ith
particle and the external force on the ith particle, respectively,
and Niext in Eq. (13) is the external torque on the ith particle,
MP is the particle mass, and IP is the particle moment ofinertia tensor. The new particle positions are calculated by

In the above equation, uP (x,tn ) and u(x,tn ) are the particle
velocity and fluid velocity at t = tn and x, respectively. The
resulting force acting on the fluid boundary nodes is given bytn +t∀P

Ri (tn + t) = Ri (tn ) +tn +t

Vi ds.

(14)tn

The density field and velocity field of the particles forthe next time step are evaluated using the updated particle
positions and velocities.

026702-3

SAEED JAFARI, RYOICHI YAMAMOTO, AND MOHAMAD RAHNAMA
0.01

0.6
0.55

Present study
Feng et al. (2004)

0.5

Present study
Feng et al. (2004)
0.005

VP/UW

0.45

YP/H

PHYSICAL REVIEW E 83, 026702 (2011)

0.4

0.35

0

0.3
0.25
0.2

0

30

60

90

-0.005

120

0

25

50

T

75

100

T

FIG. 4. Nondimensional y component of the position versus
nondimensional time.

FIG. 6. Nondimensional y component of the velocity versus
nondimensional time.

Based on the above equations, the computational scheme
can be summarized as follows:(1) At t = t0 , the fluid velocity and particle velocity (translational and rotational), together with the particle position(initial position when t = t0 ) are known.
(2) The velocity field of the particles is obtained using theparticle positions and velocities [Eqs. (5) and (6)].
(3) The fluid-solid interaction force fH is computed using
Eq. (8).
(4) The flow field is solved using the LBM [Eq. (9)].
(5) The interaction force is applied to all solid boundary
nodes; the hydrodynamic force and the hydrodynamic torqueacting on the particle are calculated according to Eqs. (10) and(11).
(6) The particle velocities and positions are updated[Eqs. (12), (13), and (14)]. The computations loop back tostep 2.

cylinder, a neutrally buoyant cylinder in simple shear flow,
two circular cylinders approaching each other in a channel,
and sedimentation of two circular cylinders in a viscous fluid.

III. RESULTS

Several problems were investigated by the present method
to show its accuracy. They include flow over a circular
0.3

Up/UW

0.25

Present study
Feng et al. (2004)

A. Flow over a circular cylinder

The flow over a fixed cylinder is first investigated toshow the accuracy of the present method. Depending on the
Reynolds number, different behavior can occur. The Reynolds
number for this method is defined by
U∞ D,νwhere D, U∞ , and ν are the diameter of the cylinder, thefree stream velocity, and the dynamic viscosity of the fluid,
respectively. The simulation was performed at Re = 20 and
40 and results were compared with previous numerical data.
Figure 2 shows the computed streamlines in a domain of L ×
W = 40D × 40D with a uniform mesh size of 800 × 800. The
cylinder center is at (2L/5,W/2) to decrease the effect of the2outer boundary. The drag coefficient [CD = Fx /(0.5ρDU∞)]
and the length of the recirculation zone (LW = 2L/D) arecompared with previous published data [19] in Table I. The
present results agree well with published data. Here, in order
to compute the drag force F x , the volume integral over theparticle volume is used [Eq. (10)].
Re =

0.2

B. A neutrally buoyant cylinder in simple shear flow

0.15

A schematic diagram of a neutrally buoyant cylinder inshear flow is shown in Fig. 3. The diameter of the cylinder is D

0.1

TABLE I. Comparison of the results of the present study with theresults of the previous study of Niu et al. [19].

0.05

Reynolds number
0

0

25

50

T

75

20

100

40
FIG. 5. Nondimensional x component of the velocity versus
nondimensional time.
026702-4

MethodCD

LW

Present study
Niu et al. [19]Present study
Niu et al. [19]

2.112
2.144
1.598
1.589

1.914
1.89
4.81
4.52

LATTICE-BOLTZMANN METHOD COMBINED WITH . . .

PHYSICAL REVIEW E 83, 026702 (2011)

d

R
U/2

-U/2

FIG. 8. Two approaching circular cylinders in a channel.

previous numerical results (Fig. 4) and shows good agreement.
Figures 5 and 6 show the comparison of two components of theparticle velocity obtained in the current study and in previous
work [18] and perfect agreement is observed. In order to show
the effect of interface thickness on the results, the upwardvelocity for different interface thicknesses is shown in Fig. 7.
It is observed that using the smoothed-profile method reduces
the small fluctuations.
C. Two circular cylinders in a channel

Results of the simulation of two circular cylinders in achannel (Fig. 8) approaching each other are presented inFig. 9. The computational domain is 2L × L, and the radiusof the circular cylinder is 0.2875L. Figure 9 displays thenondimensional normal force between two cylinders versus
the nondimensional distance between them, as they approacheach other with speeds ±U/2 = ±0.01, respectively. In thisfigure h is defined by
h = d − 2R,
where d is the distance between the two cylinders, R is theradius of the cylinder, and C is the drag force on either cylinderwhen far from the other cylinder. Figure 9 shows that the abilityof this method to capture lubrication effects is enhanced byusing a finer mesh.
103(Fx-C)/(ρνU)

Present work (L=32)
Present work (L=64)
Present work (L=128)

FIG. 7. Nondimensional y component of the velocity versus
nondimensional time for different interface thicknesses (a) ξ = 0,(b) ξ = x, (c) ξ = 2x.

and it is positioned at rest halfway between the lower wall andthe center between two parallel plates. The gap between thetwo parallel plates is 4D. The upper and lower walls are moving
at velocity UW /2 in opposite directions. The computational
domain is L × H = 100D × 4D, and D is selected to be
20. The fluid Reynolds number is Re = U W H /ν = 40. The
trajectory of the particle in the Y direction is compared with

10

2

1

10 -2
10

10

-1

h/R

10

0

FIG. 9. Nondimensional normal force between two cylinders,approaching each other with speeds ±U/2 = ±0.01 versus their
nondimensional separation.

026702-5

SAEED JAFARI, RYOICHI YAMAMOTO, AND MOHAMAD RAHNAMA
D. Sedimentation of two circular cylinders in a viscous fluid

In order to show the capability of the current method, thesedimentation of two circular particles in a Newtonian fluid
is investigated. We have considered a channel 2 cm wide (inthe x direction) and 8 cm high (in the y direction). The fluid
viscosity is 0.01 g/cm s and the density is 1 g/cm3 . The particle
density is 1.01 g/cm3 and the radius is 0.1 cm. Gravity acts

FiP −P =2.4ε

PHYSICAL REVIEW E 83, 026702 (2011)

in the negative y direction. The simulation is started at t = 0
s by dropping the two particles at the center of the channel atheights of 7.2 and 6.8 cm.
To prevent the particles from overlapping witheach other or the wall, artificial forces are added.
In the current study the Lennard-Jones potential{EP P = 0.4ε[(2a/Ri,j )12 − (2a/Ri,j )6 ]} is used, where
i and j index the particles. We have0
Ri −Rj[2(2a/R)14 − (2a/Ri,j )8 ] (2a)i,j
2 ,j =1,j =i

Ri,j > 2a + ζ,
Ri,j  2a + ζ,

(15)

0
14
8 Ri −Rj w
j w=1 [2(a/Ri,j w ) − (a/Ri,j w ) ] (a)2

Ri,j > a + ζ,
Ri,j  a + ζ,

(16)

NP

and
2.4ε

j wmax

where ε = a 2 ,Ri,j = |Ri − Rj |, and Ri,j w = |Ri − Rj w |. Rj w
represents the wall position and ζ is set to one lattice unit in
the present work.
It is known that two particles dropped close to each otherin a Newtonian fluid undergo drafting, kissing, and tumbling[18,28]. The leading particle creates a wake of low pressure, inwhich the trailing particle is caught, thus making it fall fasterthan the leading one; this stage is called drafting. Then the
trailing particle with increased speed induces a kissing contact
with the leading one in which the two particles form a long
body with the center line along the stream. This state, which iscalled kissing, is unstable, and will be broken down eventually;
after that the particles tumble.
Figures 10 and 11 show the plots of instantaneous verticaland horizontal positions of the two particles. The previous
numerical results [18] are also included and their comparison
confirms the accuracy of the present work. These figures

8

6

1.75
1.5
1.25

0.75

YP (cm)Particle 1 (Present study)
Particle 2 (Present study)
Particle 1 (Feng et al. (2004))
Particle 2 (Feng et al. (2004))

0.5
0.25
0

0

1

2

3

4

5

t (sec)FIG. 11. Plot of instantaneous horizontal positions of the twoparticles.

5

IV. SUMMARY

4
3
2
1
0

1

illustrate that the results are in good agreement with the kissingof particles. It is known that the tumbling phenomenon is arealization of instability and different modes of tumbling havebeen reported in previous work [18,28].

Particle 1 (Present study)
Particle 2 (Present study)
Particle 1 (Feng et al. (2004))
Particle 2 (Feng et al. (2004))72

XP (cm)FiP −W =

0

1

2

3

4

5

t (sec)FIG. 10. Plot of instantaneous vertical positions of the twoparticles.

The LB-SP method, which combines the smoothed-profileand the lattice-Boltzmann methods, has been proposed for
simulating particle suspensions. In this method, the same
computational grid is used for flow domain and particles. The
no-slip boundary condition is imposed by using an externalboundary force in computational nodes inside the particles.
The main advantage of the present method is preservationof the merits of the LBM in simulating fluid flow. The
LB-SP method has been validated by simulating the flow
over a circular cylinder, a neutrally buoyant cylinder in simple
shear flow, two circular cylinders approaching each other in

026702-6

LATTICE-BOLTZMANN METHOD COMBINED WITH . . .

PHYSICAL REVIEW E 83, 026702 (2011)a channel, and sedimentation of two circular cylinders in a
viscous fluid. All of the results confirm the accuracy andefficiency of the present method and its ability to simulate
fluid-particle interaction problems.

[1] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics
and Beyond (Oxford University Press, Oxford, 2001).
[2] S. Chen and G. Doolen, Annu. Rev. Fluid Mech. 30, 329(1998).
[3] D. Yu, R. Mei, L. S. Luo, and W. Shyy, Prog. Aerosp. Sci. 39,
329 (2003).
[4] C. K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech. 42, 439(2010).
[5] X. Shan and H. Chen, Phys. Rev. E 47, 1815 (1993).
[6] M. R. Swift, W. R. Osborn, and J. M. Yeomans, Phys. Rev. Lett.
75, 830 (1995).
[7] T. Inamuro, T. Ogata, S. Tajima, and N. Konishi, J. Comput.
Phys. 198, 628 (2004).
[8] A. G. Yiotis, J. Psihogios, M. E. Kainourgiakis, A. Papaioannou,
and A. K. Stubos, Colloids Surf. A 300, 35 (2007).
[9] H. Huanga and X. Y. Lu, Phys. Fluids 21, 092104 (2009).
[10] M. Cates, K. Stratford, R. Adhikari, P. Stansell, J. Desplat,
I. Pagonabarraga, and A. Wagner, J. Phys. Condens. Matter 16,
3903 (2004).
[11] K. Stratford, R. Adhikari, I. Pagonabarraga, and J. Desplat,
J. Stat. Phys. 121, 163 (2005).
[12] R. W. Nash, R. Adhikari, and M. E. Cates, Phys. Rev. E 77,
026709 (2008).
[13] A. J. C. Ladd and R. Verberg, J. Stat. Phys. 104, 1191 (2001).

ACKNOWLEDGMENTS

S.J. would like to acknowledge Professor Takaji
Inamuro and Dr. Shugo Yasuda for their useful comments anddiscussions.

[14] R. Adhikari, K. Stratford, A. J. Wagner, and M. E. Cates, Eur.
Phys. Lett. 71, 473 (2005).
[15] A. J. C. Ladd, J. Fluid Mech. 271, 285 (1994).
[16] A. J. C. Ladd, J. Fluid Mech. 271, 311 (1994).
[17] Y. H. Qian, D. D’Humiéres, and P. Lallemand, Eur. Phys. Lett.
17, 479 (1992).
[18] Z. G. Feng and E. E. Michaelides, J. Comput. Phys. 195, 602(2004).
[19] X. D. Niu, C. Shu, Y. T. Chew, and Y. A. Peng, Phys. Lett. A
354, 173 (2006).
[20] J. Wu and C. K. Aidun, Int. J. Numer. Meth. Fluids. 62, 765(2010).
[21] R. Mei, L. S. Luo, and W. Shyy, J. Comput. Phys. 155, 307(1999).
[22] M. Bouzidi, D. d’Humieres, P. Lallemand, and L. S. Luo,J. Comput. Phys. 172, 704 (2001).
[23] R. Verberg and A. J. C. Ladd, Phys. Rev. Lett. 84, 2148 (2000).
[24] R. Verberg and A. J. C. Ladd, Phys. Rev. E 65, 016701 (2001).
[25] Y. Nakayama and R. Yamamoto, Phys. Rev. E 71, 036707 (2005).
[26] X. He, Q. Zou, L. S. Luo, and M. Dembo, J. Comput. Phys. 129,357 (1996).
[27] Q. Behrend, Phys. Rev. E 52, 1164 (1995).
[28] N. A. Patankar, Annual Research Briefs-2001, Center for
Turbulent Research, Stanford University, 2001 (unpublished).

026702-7