Force evaluation in the lattice Boltzmann method involving curved geometry
Renwei Mei,
*
Dazhi Yu,
†
and Wei Shyy
‡
Department of Aerospace Engineering, Mechanics and Engineering Science, University of Florida, Gainesville, Florida 32611-6250
Li-Shi Luo
§
ICASE, MS 132C, NASA Langley Research Center, Building 1152, 3 West Reid Street, Hampton, Virginia 23681-2199
共Received 26 September 2000; revised manuscript received 7 November 2001; published 3 April 2002兲
The present work investigates two approaches for force evaluation in the lattice Boltzmann equation: the
momentum-exchange method and the stress-integration method on the surface of a body. The boundary con-
dition for the particle distribution functions on curved geometries is handled with second-order accuracy based
on our recent works 关Mei et al., J. Comput. Phys. 155, 307 共1999兲; ibid. 161, 680 共2000兲兴. The stress-
integration method is computationally laborious for two-dimensional flows and in general difficult to imple-
ment for three-dimensional flows, while the momentum-exchange method is reliable, accurate, and easy to
implement for both two-dimensional and three-dimensional flows. Several test cases are selected to evaluate
the present methods, including: 共i兲 two-dimensional pressure-driven channel flow; 共ii兲 two-dimensional uni-
form flow past a column of cylinders; 共iii兲 two-dimensional flow past a cylinder asymmetrically placed in a
channel 共with vortex shedding兲; 共iv兲 three-dimensional pressure-driven flow in a circular pipe; and 共v兲 three-
dimensional flow past a sphere. The drag evaluated by using the momentum-exchange method agrees well with
the exact or other published results.
DOI: 10.1103/PhysRevE.65.041203 PACS number共s兲: 47.10.⫹g, 47.11.⫹j, 05.20.Dd
I. INTRODUCTION
A. Background of the lattice Boltzmann equation method
The method of lattice Boltzmann equation 共LBE兲 solves
the microscopic kinetic equation for particle distribution
function f (x,
,t), where
is the particle velocity, in phase
space (x,
) and time t, from which the macroscopic quanti-
ties 共flow mass density
and velocity u) are obtained
through moment integration of f(x,
,t). Because the solu-
tion procedure is explicit, easy to implement and parallelize,
the LBE method has increasingly become an attractive alter-
native computational method for solving fluid dynamics
problems in various systems 关1–4兴. The most widely used
lattice Boltzmann equation 关1–4兴 is a discretized version of
the model Boltzmann equation with a single relaxation time
approximation due to Bhatnagar, Gross, and Krook 共BGK
model兲关5兴,
t
f⫹
•“f ⫽
1
关
f⫺ f
(0)
兴
, 共1兲
where f
(0)
is the Maxwell-Boltzmann equilibrium distribu-
tion function and is the relaxation time. The mass density
and momentum density
u are the first (D⫹ 1) hydrody-
namic moments of the distribution function f and f
(0)
in D
dimensions. It can be shown that the particle velocity space
can be discretized and reduced to a very small set of discrete
velocities
兵
␣
兩
␣
⫽ 1,2,...,b
其
, and the hydrodynamic mo-
ments of f and f
(0)
as well as their fluxes can be preserved
exactly, because the moment integral can be replaced by
quadrature exactly up to a certain order in
关6–9兴. With
velocity space
properly discretized, Eq. 共1兲 reduces to a
discrete velocity model of the Boltzmann equation:
t
f
␣
⫹
␣
•“f
␣
⫽
1
关
f
␣
⫺ f
␣
(0)
兴
. 共2兲
In the above equation, f
␣
(x,t)⬅ f (x,
␣
,t) and f
␣
(0)
(x,t)
⬅ f
(0)
(x,
␣
,t) are the distribution function and the equilib-
rium distribution function of the
␣
th discrete velocity
␣
,
respectively. Equation 共2兲 is then discretized in space x and
time t into
f
␣
共
x
i
⫹ e
␣
␦
t
,t⫹
␦
t
兲
⫺ f
␣
共
x
i
,t
兲
⫽⫺
1
关
f
␣
共
x
i
,t
兲
⫺ f
␣
(eq)
共
x
i
,t
兲
兴
,
共3兲
where
⫽ /
␦
t
is the dimensionless relaxation time and e
␣
is
a discrete velocity vector. The coherent discretization of
space and time is done in such a way that
␦
x⫽ e
␣
␦
t
is always
the displacement vector from a lattice site to one of its neigh-
boring sites. The equilibrium distribution function f
␣
(eq)
(x
i
,t)
in the lattice Boltzmann equation 共3兲 is obtained by expand-
ing the Maxwell-Boltzmann distribution function in Taylor
series of u up to second order 关6,7兴, and can be expressed in
general as
f
␣
(eq)
⫽ w
␣
冋
1⫹
3
c
2
共
e
␣
•u
兲
⫹
9
2c
4
共
e
␣
•u
兲
2
⫺
3
2c
2
u
2
册
, 共4兲
where c⬅
␦
x
/
␦
t
;
␦
x
is the lattice constant of the underlying
lattice space; and coefficient w
␣
depends on the discrete ve-
locity set
兵
e
␣
其
in D spatial dimensions. In what follows, we
shall use the lattice units of
␦
x
⫽ 1 and
␦
t
⫽ 1. The Appendix
*
Electronic address: rwm@aero.ufl.edu
†
Electronic address: ydz@aero.ufl.edu
‡
Electronic address: wss@aero.ufl.edu
§
Electronic address: luo@icase.edu; http://www.icase.edu/⬃luo
PHYSICAL REVIEW E, VOLUME 65, 041203
1063-651X/2002/65共4兲/041203共14兲/$20.00 ©2002 The American Physical Society65 041203-1
provides the details of coefficient w
␣
and the discrete veloc-
ity set
兵
e
␣
其
for the two-dimensional nine-velocity model
共D2Q9兲 and the three-dimensional nineteen-velocity model
共D3Q19兲关10兴. Figure 1 shows the discrete velocity sets of
the two models. It should be pointed out that there exist other
discrete velocity sets
兵
e
␣
其
that have sufficient symmetry for
the hydrodynamics 关6,7兴. A comparative study of three three-
dimensional LBE models including the fifteen-velocity
model 共D3Q15兲, the nineteen-velocity model 共D3Q19兲, and
the twenty-seven-velocity model 共D3Q27兲, in terms of accu-
racy and computational efficiency has been conducted by
Mei et al. 关28兴. It was found that the nineteen-velocity model
共D3Q19兲 offers a better combination of computational stabil-
ity and accuracy. The D2Q9 and D3Q19 models will be used
in this study for force evaluation in two-dimensional and
three-dimensional flows, respectively. Equation 共3兲 is conven-
ienty solved in two steps:
collision: f
˜
␣
共
x
i
,t
兲
⫽ f
␣
共
x
i
,t
兲
⫺
1
关
f
␣
共
x
i
,t
兲
⫺ f
␣
(eq)
共
x
i
,t
兲
兴
,
共5a兲
streaming: f
␣
共
x
i
⫹ e
␣
␦
t
,t⫹
␦
t
兲
⫽ f
˜
␣
共
x
i
,t
兲
, 共5b兲
which is known as the LBGK scheme 关1,2兴. The collision
step is completely local and the streaming step is uniform
and requires little computational effort, which makes Eq. 共5兲
ideal for parallel implementation. The simplicity and com-
pact nature of the LBGK scheme, however, necessitate the
use of the square lattices of constant spacing (
␦
x
⫽
␦
y
), and
consequently lead to the unity of the local Courant-
Friedrichs-Lewy number, because
␦
t
⫽
␦
x
⫽ 1.
B. Boundary condition for a curved geometry
in the LBE method
Consider a part of an arbitrary curved wall geometry, as
shown in Fig. 2, where the filled small circles on the bound-
ary, x
w
, denote the intersections of the boundary with vari-
ous lattice-to-lattice links. The fraction of an intersected link
in the fluid region, ⌬, is defined by
⌬⫽
储
x
f
⫺ x
w
储
储
x
f
⫺ x
b
储
. 共6兲
Obviously the horizontal or vertical distance between x
b
and
x
w
is ⌬
␦
x
on the square lattice, and 0⭐⌬⭐1. In Eq. 共5b兲,
the value of f
˜
␣
(x
i
,t) needs to be constructed according to
the location of the boundary and the boundary conditions, if
the grid point x
i
⫽ x
b
lies beyond the boundary. In the past,
the bounce-back boundary condition has been use to deal
with a solid boundary in order to approximate the no-slip
boundary condition at the solid boundary 关11–23兴. However,
it is well understood that this bounce-back boundary condi-
tion satisfies the no-slip boundary condition with a second-
order accuracy 共for the Couette and Poiseuille flows兲 at the
location one-half lattice spacing (⌬⫽ 1/2) outside of a
boundary node where the bounce-back collision takes place;
and this is only true with simple boundaries of straight line
parallel to the lattice grid 关18–20兴. For a curved geometry,
simply placing the boundary halfway between two nodes
will alter the geometry on the grid level and degrade the
accuracy of the flow field and the force on the body at finite
and higher Reynolds number. To circumvent this difficulty,
Mei and Shyy solved Eq. 共2兲 in curvilinear coordinates using
a finite difference method to compute f
␣
关24兴.Heand
Doolen used body-fitted curvilinear coordinates with interpo-
FIG. 1. Discrete velocity set
兵
e
␣
其
. Two-dimensional nine-
velocity 共D2Q9兲 model 共top兲. Three-dimensional nineteen-velocity
共D3Q19兲 model 共bottom兲.
FIG. 2. Layout of the regularly spaced lattices and curved wall
boundary. The circles (䊊), discs (䊉 ), shaded discs (䊊), and dia-
monds (〫 ) denote fluid nodes, boundary locations (x
w
), solid
nodes that are also boundary nodes (x
b
) inside solid, and solid
nodes, respectively.
RENWEI MEI, DAZHI YU, WEI SHYY, AND LI-SHI LUO PHYSICAL REVIEW E 65 041203
041203-2
lation throughout the entire mesh, except at the boundaries
where the bounce-back boundary condition is used 关25兴.In
the recent works of Filippova and Ha
¨
nel 关26兴 and Mei et al.
关27,28兴, a second-order accurate boundary condition for
curved geometry was developed in conjunction with the use
of Cartesian grid in order to retain the advantages of the LBE
method. An interpolation scheme is employed only at the
boundaries to obtain f
˜
␣
(x
i
,t). The detailed assessment on
the impact of the boundary condition on the accuracy of the
flow field have been given in Ref. 关27兴 for some two-
dimensional flows and in Ref. 关28兴 for some three-
dimensional flows.
Because the bounce-back type boundary conditions play
an important role in lattice Boltzmann simulations, it is im-
portant for us to understand how the boundary conditions
work. First of all, one must realize that it is impossible for
any kinetic numerical scheme to impose a given velocity 共the
Dirichlet boundary condition兲 on a given grid node, because
the Knudsen layer type of phenomena 关29–31兴 would be
manifested in kinetic schemes 关18–20,32兴. For example, in
the Poiseuille and the Couette flows, the location where hy-
drodynamic boundary conditions are satisfied are one-half
grid spacing away from the boundary grids where the
bounce-back boundary conditions are imposed 关18–20兴. For
flows around an arbitrary shaped body analytical solutions
do not exist. Nevertheless, substantial evidence shows that
the bounce-back boundary conditions combined with inter-
polations, and including the one-half grid spacing correction
at boundaries, are in fact second-order accurate and thus ca-
pable of handling curved boundaries 关22,23,25,33兴. This
point is also demonstrated in the present work.
C. Force evaluation and related works
In spite of numerous improvements in the LBE method
over the last several years, one important issue that has not
been systematically studied is the accurate determination of
the fluid dynamic force involving curved boundaries. Need-
less to say, accurate evaluation of the force is crucial to the
study of fluid dynamics, especially in fluid-structure interac-
tion. Several force evaluation schemes, including
momentum-exchange 关13,15兴 and integration of surface
stress 关25,34兴, have been used to evaluate the fluid dynamic
force on a curved body in the context of the LBE method.
He and Doolen 关25兴 evaluated the force by integrating the
total stress on the surface of the cylinder and the components
of the stress tensor were obtained by taking respective veloc-
ity gradients. Even though the body-fitted grid was used, an
extrapolation was needed to obtain the stress in order to cor-
rect the half-grid effect due to the bounce-back boundary
condition. Filippova and Ha
¨
nel 关26兴 also developed a
second-order accurate boundary condition for curved bound-
aries. However, the fluid dynamics force on a circular cylin-
der asymmetrically placed in a two-dimensional channel was
obtained by integrating the pressure and deviatoric stresses
on the surface of the cylinder by extrapolating from the
nearby Cartesian grids to the solid boundary 关26,34兴. To gain
insight into the method of surface stress integration, it is
instructive to examine the variation of the pressure on the
surface of a circular cylinder at finite Reynolds number ob-
tained by using the LBE method for flow over a column of
cylinders 共see Ref. 关27兴, and Sec. III B兲. Figure 3 shows the
pressure coefficient
C
P
⫽
p⫺ p
⬁
1
2
U
2
,
on the surface obtained by using second-order extrapolation,
where p
⬁
is the far upstream pressure. Only those boundary
points, x
w
, intersected by the horizontal or vertical veloci-
ties, i.e., e
1
, e
3
, e
5
, and e
7
, are considered in the result given
by Fig. 3. If the boundary points intersected by the links in
the diagonal velocities, i.e., e
2
, e
4
, e
6
, and e
8
, are also con-
sidered, the variation of C
P
would be more noisy. The com-
ponents of the deviatoric stress tensor show a similar noisy
pattern. It is not clear how the noise in the pressure and
stresses affect the accuracy of the fluid dynamic force in the
stress-integration method. While the programing in the ex-
trapolation and integration is manageable in two-dimensional
cases, it is rather laborious in three-dimensional 共3D兲 cases.
In Fig. 3, the LBE result of C
P
(
) 共indicated by symbol ⫻ )
is compared with that obtained by using a 3D multiblock,
body-fitted coordinates, and pressure-based Navier-Stokes
solver 关35–37兴 with a much finer resolution: 201 points
around the cylinder and the smallest grid size along the ra-
dial direction dr⫽ 0.026 共relative to the radius r⫽ 1). Not
surprisingly, the result obtained by using the Navier-Stokes
solver with body-fitted grids and a much finer resolution is
smoother than the LBE result with a Cartesian grid of
coarser resolution. Nevertheless, the LBE solution still es-
sentially agrees with the Navier-Stokes solution.
Instead of the stress-integration method, Ladd used the
momentum-exchange method to compute the fluid force on a
sphere in suspension flow 关13兴. In the flow simulation using
the bounce-back boundary condition, the body is effectively
replaced by a series of stairs. Each segment on the surface
FIG. 3. Distribution of the pressure coefficient C
P
on the sur-
face of a 2D circular cylinder of radius r⫽ 6.6, and center-to-center
distance H/r⫽ 20. The stagnation point is located at
⫽ 180°. The
LBE result denoted by symbol ⫻ is obtained with
⫽ 0.6 and Re
⫽ 40. The solid line is the result obtained by using a 3D multiblock,
body-fitted grid, and pressure-based Navier-Stokes solver with a
much finer resolution.
FORCE EVALUATION IN THE LATTICE BOLTZMANN . . . PHYSICAL REVIEW E 65 041203
041203-3
has an area of unity for a cubic lattice. The force on each link
关halfway between two lattices at x
f
and x
b
⫽ (x
f
⫹ e
␣
␦
t
)in
which x
b
resides in the solid region兴 results from the
momentum-exchange 共per unit time兲 between two opposing
directions of the neighboring lattices
1
␦
t
关
e
␣
f
␣
共
x
f
兲
⫺ e
␣
¯
f
␣
¯
共
x
f
⫹ e
␣
␦
t
兲
兴
,
in which e
␣
¯
⬅⫺ e
␣
. Whereas the momentum-exchange
method is very easy to implement computationally, its appli-
cability and accuracy for a curved boundary have not been
systematically studied. To recapitulate, there are two major
problems associated with the method of surface stress inte-
gration. First, the components of the stress tensor are often
noisy on a curved surface due to limited resolution near the
body and the use of Cartesian grids. The accuracy of such a
method has not been addressed in the literature. Second, the
implementation of the extrapolation for the Cartesian com-
ponents of the stress tensor to the boundary surface and the
integration of the stresses on the surface of a three-
dimensional geometry are very laborious in comparison with
the intrinsic simplicity of the lattice Boltzmann simulations
for flow field. The problems associated with the method of
the momentum-exchange are as follows. 共a兲 The scheme was
proposed for the case with ⌬⫽ 1/2 at every boundary inter-
section x
w
. Whether this scheme can be applied to the cases
where ⌬⫽1/2, when, for example, the boundary is not
straight, needs to be investigated. 共b兲 As in the case of stress
integration method, the resolution near a solid body is often
limited and the near wall flow variables can be noisy. If one
uses the momentum-exchange method to compute the total
force, it is not clear what the adequate resolution is to obtain
reliable fluid dynamic force on a bluff body at a given 共mod-
erate兲 value of Reynolds number, say, Re⬇O(10
2
).
D. Scope of the present work
In what follows, two methods of force evaluation, i.e., the
stress-integration and the momentum-exchange methods,
will be described in detail. The shear and normal stresses on
the wall in a pressure driven channel flow will first be exam-
ined to assess the suitability of the momentum-exchange
method when ⌬⫽1/2 and analyze the errors incurred. The
results on the drag force for flow over a column of circular
cylinders using these two methods will be subsequently as-
sessed for the consistency. The drag coefficient at Re⫽ 100
will be compared with the result of Fornberg 关38兴 obtained
by using a second-order accurate finite difference scheme
with sufficient grid resolution. For flow over a cylinder
asymmetrically placed in a channel at Re⫽ 100, the unsteady
drag and lift coefficients were computed and compared with
the results in the literature. The momentum-exchange
method is further evaluated for three-dimensional fully de-
veloped pipe flow and for a uniform flow over an two-
dimensional array of spheres at finite Reynolds number. We
found that the simple momentum-exchange method for force
evaluation gives fairly reliable results for the two-
dimensional and three-dimensional flows.
II. METHODS FOR FORCE EVALUATION
IN THE LBE METHOD
A. Second-order accurate no-slip boundary condition
for curved geometry
The analysis of boundary conditions for a curved bound-
ary in the lattice Boltzmann equation is accomplished by
applying Chapman-Enskog expansion for the distribution
function at the boundary. The following approximation for
postcollision distribution function on the right-hand side of
Eq. 共5b兲 can lead to a second-order accurate no-slip bound-
ary condition 关26–28兴
f
˜
␣
¯
共
x
b
,t
兲
⫽
共
1⫺
兲
f
˜
␣
共
x
f
,t
兲
⫹
f
␣
*
共
x
b
,t
兲
⫹ 2w
␣
3
c
2
e
␣
¯
•u
w
,
共7兲
where
f
␣
*
共
x
b
,t
兲
⫽ w
␣
共
x
f
,t
兲
冋
1⫹
3
c
2
共
e
␣
•u
bf
兲
⫹
9
2c
4
共
e
␣
•u
f
兲
2
⫺
3
2c
2
u
f
2
册
⫽ f
␣
(eq)
共
x
f
,t
兲
⫹ w
␣
共
x
f
,t
兲
3
c
2
e
␣
•
共
u
bf
⫺ u
f
兲
, 共8兲
and
u
bf
⫽ u
ff
⫽ u
f
共
x
f
⫹ e
␣
¯
␦
t
,t
兲
,
⫽
共
2⌬⫺ 1
兲
共
⫺ 2
兲
,0⭐⌬⬍
1
2
,
共9a兲
u
bf
⫽
1
2⌬
共
2⌬⫺ 3
兲
u
f
⫹
3
2⌬
u
w
,
⫽
共
2⌬⫺ 1
兲
共
⫹ 1/2
兲
,
1
2
⭐⌬⬍ 1. 共9b兲
The above treatment is applicable for both the two-
dimensional and three-dimensional lattice Boltzmann mod-
els.
By substitution of Eq. 共8兲, Eq. 共7兲 becomes
f
˜
␣
¯
共
x
b
,t
兲
⫽ f
˜
␣
共
x
f
,t
兲
⫺
关
f
˜
␣
共
x
f
,t
兲
⫺ f
␣
(eq)
共
x
f
,t
兲
兴
⫹ w
␣
共
x
f
,t
兲
3
c
2
e
␣
•
共
u
bf
⫺ u
f
⫺ 2u
w
兲
. 共10兲
Thus, the above treatment of curved boundary can be
thought as a modification of the relaxation 共the viscous ef-
fect兲 near the wall 共with the relaxation parameter
), in ad-
ditional to a forcing term accounting for the momentum-
exchange effect due to the wall.
RENWEI MEI, DAZHI YU, WEI SHYY, AND LI-SHI LUO PHYSICAL REVIEW E 65 041203
041203-4
B. Force evaluation based on stress integration
He and Doolen 关25兴 evaluated the force by integrating the
total stresses on the boundary of the cylinder
⍀,
F⫽
冕
⍀
dAn
ˆ
•
兵
⫺ pI⫹
关共
“:u
兲
⫹
共
“:u
兲
T
兴
其
, 共11兲
where n
ˆ
is the unit out normal vector of the boundary
⍀.In
Ref. 关25兴, a body-fitted coordinate system together with grid
stretching was used such that a large number of grids can be
placed near the body to yield reliable velocity gradient
i
u
j
.
In general, u is not the primary variable in the LBE simula-
tions and the evaluation of u using
兺
␣
e
␣
f
␣
based on f
␣
’s
suffers a loss of accuracy due to the cancellation of two close
numbers in f
␣
’s; consequently, the evaluation of the deriva-
tive
i
u
j
will result in further degradation of the accuracy.
Filippova 关34兴 used a similar integration scheme to obtain
the dynamic force on the body for the force on a circular
cylinder 关26兴 except that the deviatoric stresses were evalu-
ated using the nonequilibrium part of the particle distribution
function 关see Eq. 共13兲 below兴. However, since the Cartesian
grid was used, the stress vectors on the surface of the body
共with arbitrary ⌬) have to be computed through an extrapo-
lation procedure based upon the information in the flow field.
This leads to further loss of accuracy for finite lattice size
␦
x
when the shear layer near the wall is not sufficiently re-
solved.
In Eq. 共11兲, the pressure p can be easily evaluated using
the equation of state p⫽ c
s
2
. For D2Q9 and D3Q19 models,
c
s
2
⫽ 1/3 so that p⫽
/3. The deviatoric stress for two-
dimensional incompressible flow
ij
⫽
共
i
u
j
⫹
j
u
i
兲
, 共12兲
can be evaluated using the nonequilibrium part of the distri-
bution function f
␣
(neq)
⫽
关
f
␣
⫺ f
␣
(eq)
兴
ij
⫽
冉
1⫺
1
2
冊
兺
␣
f
␣
(neq)
共
x,t
兲
冉
e
␣
,i
e
␣
,j
⫺
1
D
e
␣
•e
␣
␦
ij
冊
,
共13兲
where e
␣
,i
and e
␣
,j
are ith and jth Cartesian component of
the discrete velocity e
␣
, respectively. For the flow past a
circular cylinder, a separate set of surface points on the cyl-
inder can be introduced in order to carry out the numerical
integration given by Eq. 共11兲. The values of the pressure and
each of the six components of the symmetric deviatoric
stress tensor on the surface points can be obtained using a
second-order extrapolation scheme based on the values of p
and
ij
at the neighboring fluid lattices. The force exerting on
the boundary
⍀ is computed as
F⫽
冕
⍀
dA n
ˆ
•
兵
⫺ pI⫹
关共
“:u
兲
⫹
共
“:u
兲
T
兴
其
extrapolated
.
共14兲
It is worth commenting here that for the two-dimensional
flow past a cylinder, nearly half of the entire code was taken
up by the above force evaluation procedure.
C. Method based on the momentum-exchange
In order to employ the momentum-exchange method effi-
ciently, two scalar arrays, w(i,j) and w
b
(i,j) are introduced.
A value of 0 is assigned to w(i,j) for the lattice site (i,j)
that are occupied by fluid; a value of 1 is assigned to w(i,j)
for those lattice nodes inside the solid body. The array
w
b
(i,j) is set to zero everywhere except for those boundary
nodes, x
b
, where a value of 1 is assigned. For a given non-
zero velocity e
␣
,e
␣
¯
denotes the velocity in opposite direc-
tion, i.e., e
␣
¯
⬅⫺ e
␣
共see Fig. 2兲. For a given boundary node
x
b
inside the solid region with w
b
(i,j)⫽ 1 and w(i,j)⫽ 1,
the momentum-exchange with all possible neighboring fluid
nodes over a time step
␦
t
⫽ 1is
兺
␣
⫽0
e
␣
关
f
˜
␣
共
x
b
,t
兲
⫹ f
˜
␣
¯
共
x
b
⫹ e
␣
¯
␦
t
,t
兲
兴关
1⫺ w
共
x
b
⫹ e
␣
¯
␦
t
兲
兴
.
Simply summing the contribution over all boundary nodes x
b
belonging to the body, the total force 共acted by the solid body
on the fluid兲 is obtained as
F⫽
兺
all x
b
兺
␣
⫽0
e
␣
关
f
˜
␣
共
x
b
,t
兲
⫹ f
˜
␣
¯
共
x
b
⫹ e
␣
¯
␦
t
,t
兲
兴
⫻
关
1⫺ w
共
x
b
⫹ e
␣
¯
␦
t
兲
兴
. 共15兲
In the momentum-exchange method the force F is evaluated
after the collision step is carried out and the value of f
˜
␣
¯
at
boundary given by Eq. 共7兲 has been evaluated. The
momentum-exchange occurs during the subsequent stream-
ing step when f
˜
␣
¯
(x
b
,t⫹
␦
t
) and f
˜
␣
(x
f
,t⫹
␦
t
) move to x
f
and
x
b
, respectively. As mentioned in the introductory section,
the effect of variable ⌬ is not explicitly included, but it is
implicitly taken into account in the determination of
f
˜
␣
¯
(x
b
,t⫹
␦
t
). The applicability of Eq. 共15兲 will be examined
and validated.
Clearly, the force is proportional to the number of bound-
ary nodes x
b
in the above formula of F and the number of the
boundary nodes increase linearly with the size of the body in
a two-dimensional flow. However, since the force is normal-
ized by
U
2
r in the formula for C
D
in two-dimensions 关see
Eq. 共24兲兴, the drag coefficient C
D
should be independent of
r.
FIG. 4. The channel flow configuration in the LBE simulations
with an arbitrary ⌬.
FORCE EVALUATION IN THE LATTICE BOLTZMANN . . . PHYSICAL REVIEW E 65 041203
041203-5
![FIG. 5. The LBE simulations of the channel flow, withD 50.2, 1/3, 0.5, and 0.7. The pressure drop is]xp521.0310 26 in lattice units.~a! Ratio between the wall forcern]yuxuy5H evaluated by using Eq.~18!, and the exact value2txy,w](/figures/fig-5-the-lbe-simulations-of-the-channel-flow-withd-50-2-1-3-6p5s7tlf.webp)