NASA/CR-2002-211662

ICASE Report No. 2002-22

y

Force Evaluation in the Lattice Boltzmann Method

Involving Curved Geometry

Renwei Mei, Dazhi Yu, and Wei Shyy

University of Florida, Gainesville, Florida

Li-Shi Luo

ICASE, Hampton, Virginia

July 2002
The NASA STI Program Office... in Profile

Since its founding, NASA has been dedicated

to the advancement of aeronautics and space

science. The NASA Scientific and Technical

Information (STI) Program Office plays a key

part in helping NASA maintain this

important role.

The NASA STI Program Office is operated by

Langley Research Center, the lead center for

NASA's scientific and technical information.

The NASA STI Program Office provides

access to the NASA STI Database, the

largest collection of aeronautical and space

science STI in the world. The Program Office

is also NASA's institutional mechanism for

disseminating the results of its research and

development activities. These results are

published by NASA in the NASA STI Report

Series, which includes the following report

types:

TECHNICAL PUBLICATION. Reports of

completed research or a major significant

phase of research that present the results

of NASA programs and include extensive

data or theoretical analysis. Includes

compilations of significant scientific and

technical data and information deemed

to be of continuing reference value. NASA's

counterpart of peer-reviewed formal

professional papers, but having less

stringent limitations on manuscript

length and extent of graphic

presentations.

TECHNICAL MEMORANDUM.

Scientific and technical findings that are

preliminary or of specialized interest,

e.g., quick release reports, working

papers, and bibliographies that contain

minimal annotation. Does not contain

extensive analysis.

CONTRACTOR REPORT. Scientific and

technical findings by NASA-sponsored

contractors and grantees.

CONFERENCE PUBLICATIONS.

Collected papers from scientific and

technical conferences, symposia,

seminars, or other meetings sponsored or

cosponsored by NASA.

SPECIAL PUBLICATION. Scientific,

technical, or historical information from

NASA programs, projects, and missions,

often concerned with subjects having

substantial public interest.

TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific

and technical material pertinent to

NASA's mission.

Specialized services that complement the

STI Program Office's diverse offerings include

creating custom thesauri, building customized

data bases, organizing and publishing

research results.., even providing videos.

For more information about the NASA STI

Program Office, see the following:

• Access the NASA STI Program Home

Page at http://www.sti.nasa.gov

• Email your question via the Internet to

help@sti.nasa.gov

• Fax your question to the NASA STI

Help Desk at (301) 621-0134

• Telephone the NASA STI Help Desk at

(301) 621-0390

Write to:

NASA STI Help Desk

NASA Center for AeroSpace Information

7121 Standard Drive

Hanover, MD 21076-1320
NASA/CR-2002-211662

ICASE Report No. 2002-22

y

Force Evaluation in the Lattice Boltzmann Method

Involving Curved Geometry

Renwei Mei, Dazhi Yu, and Wei Shyy

University of Florida, Gainesville, Florida

Li-Shi Luo

ICASE, Hampton, Virginia

ICASE

NASA Langley Research Center

Hampton, Virginia

Operated by Universities Space Research Association

Prepared for Langley Research Center

under Contract NAS 1-97046

July 2002
Available from the following:

NASA Center for AeroSpace Information (CASI)

7121 Standard Drive

Hanover, MD 21076-1320

(301) 621-0390

National Technical Information Service (NTIS)

5285 Port Royal Road

Springfield, VA 22161-2171

(703) 487-4650
FORCEEVALUATIONIN THE LATTICE BOLTZMANNMETHOD INVOLVING

CURVEDGEOMETRN

RENWEI MEI*,DAZHIyUI,WEISHYY$, ANDLI-SHILUG§

Abstract, Tilepresentworkinvestigatestwoapproaches [brforceewfluationin thelatticeBoltzmann

equation:themomentum-exchange methodandthestress-integration methodonthesurfa.eeof a body.

Theboundaryconditionfortheparticledistributionfunctionsoncurve(]geometriesishandledwithsecond

orderaccuracybase(]onourrecentworks.Thestress-integration methodiseornputg_ionally laboriousfor

two-dimensional flowsandingeneral difficulttoimplement forthree-dimensional flows,whilethemomentumexchangemethodisreliable,accurate, andeasytoimplement tbrbothtwo-dimensional andthree-dimensional

flows.Several testcasesareselectedtoevaluatethepresentmethods, including:(i)two-dimensional pressuredrivenchannel flow;(ii)two-dimensional urfiformflowpastacolumnofcylinders;(iii)two-dimensional flow

pastacylinderasymmetricallyplacedinachannel(withvortexshedding); (iv)three--dimensional pressure--

drivenflowin acircularpipe;and(v) three-dimensional flowpastasphere.Thedragevaluatedbyusing

themomentum--exchange methodagreeswellwiththeexactorotherpublishedresults.

Key words, lattice Boltznmnn method, t'erce evaluasiou on fluid-solid interface, momentum-exchange

method, stress-integration method, boundary condition tbr curved geometries, accuracy, 3-D flows

Subject classification. Fhfid Mechanics

1. Introduction,

1.1. 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 quantities (flow mass

density p and velociw u) are obtained through moment integration of f(x, _, t). Because the solution procedure is explicit, easy to in_plerneni, and parallelize, the LBE method has increasingly become an attractive

alternative computational method for solving fluid dynamics problems in various systems [1, 2, 3, 4]. The

most widely used lattice Boltzmann equation [I, 2, 3, 4] is a discretized version of the model Boltzmann

equation with a single relaxation time approxinmt;ion due t;o Bhatnagar, Gross, and Krook (BGK model)

l

a_.f + _. Vf = _[f- f(o)], (1.1)

where f(o) is the Maxwell-Boltzmann equilibrium distribution function and _ is the relaxation time. The

mass density p and momenum_ density pu are the first (D + 1) hydrodynamic moments of the distribution

*Department of Aerospace Engineering, Mechanics & Engineering Science, University of Florida., Gainesville, PL 32611-6250

(email address: rwm_aero.ufl.edu).

1Department of Aerospace Engineering, Mechanics & Engineering Science, University of Florida, Gainesville, FL 32611-6250

(email address: ydz(C..aero.ufl.edu).

*Department of Aerospaee Engineering, Mechanics &: Engiaeering Science, University of Florida, Gainesville, FL 32611--6250

(email address: ws@aero.ufl.edu).

§ICASE, Mail Stop 132C, NASA Langley Research Center, 3 West Reid Street, Building 1152, Hampton, VA 23681-2199

(email a.ddress: luo@icase.edu). This research wa.s supported by the National Aeronautics and Space Administration under

NASA Contract No. NAS1-.97046 while the author was in residence at ICASE, NASA Langley l:_esea.reh Center, ttampt;on, VA

23681--2199.
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 {g_la = 1, 2,..., b}, and the hydrodynamic moments

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, 7, 8, 9]. With velocity space { properly discretized,

Eq. (1.1) reduces to a discrete velocity model of the Bol{:zmann equation

Otf_ + _" g/,-_ = lA[f_ ----f(0)]_. _, (11.2)

In the above equation, f,_(x, t) :---iI(x, _, t) and f_o)(x, t) :----f(o)(x, _, t) are the distribution flmction and

the equilibrium distribution funct.ion of the (tth discrete velocity {_, respectively. Equat, ion (1.2) is then

discretized in space x and time t into

where _- .....t/6t is _he dimensionless relaxation time and e_ is a discrete velocity vector. The coherent

discretization of space and time is done in such a way that (ix = e_dt is always the displa.cement vector from

a lat,tice site to one of its neighboring sites. The equilibrium distribution function fc_eq)(xi, t) in t,he ]attice

Boltzmann equation (1.3) is obtained by expanding the Maxwell--Boltzmann distribution function in Taylor

series of "a up to second order [6, 7], and can be expressed in general as

[ 3, u) + 9 , - _'tt 2] (1.4)

where c _= fixSt; 5._:is the lattice constant of the underlying lattice space; and coefficient 'w_ depends on the

discrete velocity se_ {Co,} in D spatial dimensions. In what follows, we shall use the lattice units of 5x = ]

and (it .....1. The Appendix provides the details of coefficient; w_ and the discrete velocity set; {e(_} for the twodimensional nine-velocity model (D2Q9) and the three-dimensional nineteen-velocity model (D3Q19) [10].

Figure 1 shows the discret, e velocity sets of the two models. It should be pointed out that there exist other

discrete velocity sets {Co.} that have the sufficient symmetry for hydrodynamics [6, 7]. A comparative study

of three three-dimensional LBE models including the fi[_een-velocity mode] (D3Q15), the nineteen-velocity

model (D3Q19), and the twenty-seven-velocity model (D3Q27), in terms of accuracy and computational

efficiency has been conducted by Meict aI. [111]. It was found that the nineteen-velocity model (D3Q]9)

offers a better combination of computational stability and accuracy. The D2Q9 and D3Q19 models will be

used in this study for force e_duation in two-dimensional and three-dimensional flows, respectively. Equat;ion

(1.3) is conveniently solved in two steps

1

collision: L(_, t) = fo(x_,t) .....:; ,_

streaming: f_(xi + Coit, t + 5t) = L(xi, l), (1.5b)

which is known as t;he LBGK scheme [1, 2]. The collision step is completely local and the streaming step

is uniform and requires little computational effort, which makes Eq. (1.5) ideal for parallel inlplementation.

The simplicity and compact nature of the LBGK scheme, however, necessitate the use of the square lattices

of constant spacing (5,_. = du) , and consequently lead to the unity of the local Courant-Friedrichs-Lewy

number, because 5t = _._:= 1.

1..2. 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 boundary, x_., denote
,,/! 3, 713

18 " "

..' " 9_ 15

7 .... t........ _-- --I -- --#

, ,/12

,/ 14 Y

F_G. ]. Discrete velocity _'et {e(_}. 0eft) Two-dimensional nine-velocity (D2Q9) 'model. (right) Three-dimensional

nineteen-vdocity (D3Q19) model

the intersect;ions of the boundary with various latl:.ice-to-lattice links. The fract.ion of an int:.ersected link in

the fluid region, A, is defined by

A .....Ilxf - x"'ll (1.6)

Ilxs - xbll

Obviously the horizontal or" vertical distance between x!, and x_,, is AS:_,on the square lattice, and 0 <_:_ < ].

In Eq. (1.5b), the value of ]_ (xi, t) needs to t)e constructed according to the location of the boundary and

the bounda,ry conditions if the grid point xi = :*:_ lies t)eyond the boundary, In the pa,st, the bounce.-.

b_ck bounda.ry condition has been vise to deal with a solid boundary in order to approxirnate the no-.slip

boundary condition at:. the solid boundary [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 2.-1]. However, it is

well understood that this bourlce-back boundary condition satisfies _;he no-slip boundary condition with a

second-order accuracy (for the Couette and Poiseuille flows) at the location one half lattice spacing (A = 1/2)

outside of' a boundary node where the boimce-back collision takes place; and this is only true with simple

boundaries of straight lines parallel to the lattice grid [19, 20, 21]. 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

tile flow field and the force on tile body at finite and higher Reynolds number. To circumvent this difficulty,

Mei and Shyy solved Eq. (1.2) in eurvilinear coordinates using a finite difIi?rence method to compute f_

[25]. He and Doolen used body--fitted curvilinear coordinates with interpolation throughout tile entire mesh,

except; at the boundaries where the bounce-back boundary condition is used [26]. In the recent works of

Filippova and Hgnel [27] and Mei e_ al. [28, 11], a second-order accurate boundary condition for curved

geometry was developed in conjunction with the use of Cartesian grids in order to retain the advantages

of the LBE method. An itlterpolation scheme is employed only at the boundaries to obtain ]'_(xi, t). The

detailed assessment on the impact of the boundary condition on the accuracy of the flow field has been given

in Ref. [28] for some two-dimensional flows and in Ref. [11] ['or some three-din, ensional flows.

Because the bounce-back type boundary conditions play an important role in lattice Boltzmann simula-.

tions, it is important for us to understand how the boundary conditions work. First of all, one must realize

that:, it:. is irnt)ossible 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, 30, 31] would be manifested in kinet.ie schemes [32, 19.2(., 21]. l_k)r example, in the Poiseuille and the Couette flows, the location

where hydrodynamic boundary conditions are satisfied are one-h_lf grid spacing away from the boundary

grids where the bounce-back boundary conditions are imposed [19, 20, 21]. For flows around an arbitrary

shaped body analytical solutions do not. exist. Nevertheless, substantial evidence shows that the bounce-
F_(',. 2. Layout of the r_gularly spaced lattices and curved wall boundary. The circles (o), discs (®), shaded discs (,*), and

diamo.nds (_) denote flv, id _odes, boundary locations (xw), solid nodes which are also boundary nodes (xb) inside solid, and

solid nodes, respectively.

back boundary conditions combined with interpolations, and including the one-half grid spacing correction

at boundaries, are in fact second-order accurate and capable of handling curved boundaries [26, 23, 24, 33].

This poin_ is also demonstrated in the present work.

1.3. Force evaluation and related works. In spite of numerous improvement for the LBE method

during the last several years, one important issue tha.t has not been systematically studied is the accurate

determination of the fluid dynamic force involving curved boundaries. Needless to say, accurate evaluation

of the force is crucial to "bhe study of fluid dynamics, especially in fluid-structure interaction. Several force

evaluation schemes, including momentmn exchange [14, 16] and integration of surface stress [26, 34], have

been used to evaluate the fluid dynamic force on a curved body in the context of the [,BE method.

He and Doolen [26] 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 velocity gradients. Even though a

body-fitted grid was used, an extrapolation was needed to obtain i;he stress in order to correct the halfgrid-cell spacing effect due to the bounce-back boundary condition. Filippova and Hfinel [27] developed

a second-order accura.te boundary condition for curved boundaries. However, the fluid dynamics force on

a circular cylinder asymmetrically placed in a two-dimensional channel was obtained by integrating the

pressure and deviat.oric stresses on the surt)ace of _he cylinder by extrapolating fi'om the nearby Cartesian

grids to t,he solid boundary [27, 34]. To gain insight into the method of surface stress integration, it is

instructive to examine the variation of the pressure on the surfa.ce of a circula.r cylinder at finite Reynolds

number obtained by using the LBE method for flow over a column of cylinders (see Ref. [28], and Sec. 3.2).

Figure 3 shows the pressure coefficient

Cp - p - p_'

1⁄2PU "_
1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5

0 o

, , , , , , , , I , , , , , , , , I , , , , , , , , I , , , , , , , ,

90 ° 180 ° 270 ° 360 c

9

FI('. 3. Flow past a colu'mn of 2D circut, ar cylinders. Distribution oj" the pressure coe_[ieient Cp o'n the surface of a 2D

circular cylinder of radius r -- 6.6, and center-to-center distance H/r -- 20. The stagnation point is located at 0 -- 180 ° . The

LBE *_sult denoted by s,;rnbols × is obtained with r -- O.(i and Re - 40. The solid line is the result obtained b_j usin9 a 3D

multi-block, body-fitted grid, and pressure-based Navier-Stokes solver with a much finer resolution.

on the surface obtained by using second-order extrapolation, where Po_ is the far upstream pressure. Only

those b()undary points, x_., intersected by the horizomal or vertical velocities, i.e., e_, e3, ers, and eT, are

considered in the result given by Fi& 3. If the boundary points intersected by the links in the diagonal

velocities, i.e., e2, e4_ e6, and es, are also considered, the variation of C1_ would be more noisy. The

components of the deviatoric stress tensor show a similar noisy pat, tern. It is not clear how the noise in

the pressure and stresses aitbct the accuracy of the fluid dynamic force in the stress-integration method.

While the programming in the extrapolation and integration is manageable in two--dimensional (3D) cases,

it is rather laborious in three-dimensional cases. In Fig. 3, the LBE result of Cp(0) (indicated by symbol

x) is compared with t;hat obtained by using a 3D multiblock, body-fitted coordina_ses, an(] pressure-based

Navier-Stokes solver [35, 36, 37] with a nmch finer resolution: 201 points around the cylinder and the smallest

grid size along the radial direction dr ..... 0.026 (relat;ive to _...... 1). Not; surprisingly, the result; obtained by

using the Navier-Stokes solver with body-fitted grid and high resolution is smoother than the LBE result

with a Cartesian grid of coarser resolution. Nevertheless, the LBE solution still essentially agrees with _.he

Navier-Stokes solution.

Instead of the st;ress-int;egration met;hod, Lad(] used the momentum-exchange method to compute _he

fluid force on a sphere in suspension flow [14]. 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 has an area of

unity for a cubic lattice. The fbrce on each link [halfway between two lattices at a:f and xb = (xf + e_jt)

in which xb resides in the solid region] results from the momentum exchange (per unit time) between two

opposing directions of the neighboring lattices

ec_.f a (x f

in which e_ -7 ---e_. Whereas the momentum-exchange method is very easy to implement computationally,

its applicability and accuracy for a curved boundary have not been systematically studied. To recapi_ulate,
there are two major problerns associated with the method of surface stress integration. First, the components

of stress tensor are often noisy on a curved surface due to limited resolution near the body and tile use of

Cartesian grids. The accuracy of such a method has not been addressed in the literature. Second, the

irnplementa.tion of the extrapolation for Cartesian eornponents of the stress tensor to the boundary surfa.ee

and t;he integration of t:he stresses on the surface of a three-dimensional geometry are very laborious in

comparison with the intrinsic sirnplicit:y of the lattice Boltzmann simulations for flow field. The problems

associated with the method of n_omentum exchange are as follows. (a) The scheme was proposed for the case

with A = 1/2 at every boundary intersection x_,. Whether this scheme can be applied to the cases where

A ¢ 1/2, when, for example, the boundary is not straight, needs to be investigated. (b) As in the ease of

stress-integration method; the resolution near a solid body is often lbnited and the near wall flow variables

can be noisy. If one uses the momentum--exchange method to compute the total force, it is not elea.r what

the adequate resolution is to obtain reliable fluid dynamic fbrce on a bluff body at a given (moderate) value

of Reynolds number, say, Re _ 0(102).

1.4. Scope of the present work. In what follows, two methods for t:he force evaluat;ion, i.e., the

st;ress-int;egration 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 be first, examined to assess 5he suitability of the

momentum-exchange method when A # 1/2 and analyze the errors incurred. The results on the drag three

for flow over a column of circular cylinders using these two methods will be subsequently assessed ['or the

consistency. The drag coefficient at Re = 100 are compared with the result of Fornberg [38 i obtained by

using a second--order a.eeurate finite difference scheme with sufficient grid resolution. For flow over a cylinder

asymmetrically placed in a channel a'b Re .... 100, the unst;eady drag anti lift coefficients are eompu_;ed and

compared with the results in _;he literaeure. The momentum-exchange method is further evaluated for threedimensional fully developed pipe flow and for a uniform flow over a two-dirnensional array of spheres at finite

Reynolds number. We found that; _;he simple momentum-exchange medlod for force evulua_;ion gives Nirly

reliable results for the two-dimensional and three-dimensional flows.

2. Methods for Force Evaluation in LBE Method.

2.1. Second-order accurate no-slip boundary condition for curved geometry. The analysis of

boundary conditions for a curved boundary in the lattice Boltzmann equation is accomplished by applying

the Chapman-Enskog expansion for t;he distribution function at the boundary. The following approximation

for the post-collision dist;ribution function on t;he right>hand side of Eq. (1.Sb) can lead t;o a second-order

aceura.te no-slip boundary condition [11, 27, .8 i

3

]; (xb, t) = (1 - X)/_:(xf, t) + :_f2(;ev, t) + 2w_.pTe_. _.,,, (2.1)

where

and

f_ I,"b, _) = wop(xf, t) 1 + . ,_vf) ÷ . -

c2 ' " zc 2c ]

3

= f2_q)(xs, _)+ ,,_'_,(x_., t)2;_-e,_. (_bs ....._s),

(2A -1)

UbS .....UH .....US(X s+e_`6t,t), X ..... (_---2) ' O<A<[,

1 3 (2A - 1) 1

"b_" .... }Z(2A - 3)W + ,T,A_,, _ .....(_-+ 1/'2) ' _ < _ < 1.

(2.'2)
The above treatment is applicable for both the two--dimensional and three--dimensional lattice Boltzmann

models.

By substitution of Eq. (2.2), Eq. (2.1) becomes

" _(eq),.x t'}]

3

+w_,p(xf.t)_e_ . (Ubf -- uj, -- 2uw). (2.4)

, , 7a , .

Thus, the above treatment of curved boundaries can be thought as a modification of the relaxation (the

viscous etIi_ct) near the wall (with the relaxation parameter ;Xi), in additional to a tbrcing term accounting

for the momentum-exchange effect due t.o the wall.

2.2. Force evahmtion based on stress integration. He and Doolen [26] evaluated the force by

integrating the total stresses on the boundary of the cylinder 0g_,

F = f d._,_. {----pl + p,,[(V :,*:) + (iV:'_t)r] } , (2.5)

J0f2

where fi is the unit out normal vector of the boundary 0_, I is the identity tensor of second rank, g : u

denotes the second rank tensor whose components are 0iu, and T is the tranpose operator. In Ref. [26], 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 a reliable velocity gradient. Oi'_l,j, Irl general, since u is not the primary

variable in the [,BE simulations and the evaluat.ion of u using _o: e_f_ based on fi_,'s suffers the loss of

accuracy due to the cancellation of two (:lose numbers in f_'s the evaluation of the derivative 0iu.i will

result in further degradation of the accuracy. I_'ilippova [34] used a similar integration scheme to obtain the

dynamic force on the body for the force on a circular cylinder [27] except that the deviatorie st.resses were

evaluated using the non-equilibrium part of the pz-;rticle dist:.ribution funct.ion [see Eq. (2.7) below]. However,

since a C.aitesian grid was used, the st.ress vectors on the surfaee of the body (with arbitrary A) have to

be computed through an extrapolation procedure based upon the information in the flow field. This leads

to further loss of accuracy tbr a finite lattice size d_ when the shear-layer near the wall is not sufficiently

resolved.

In Eq. (2.5), the pressure # can be easily evaluated using the equation of state p = c_o. For D2Q9 and

D3Q]9 models, c_ = 1/3 so that p = p/3. The deviatoric stress for two--dimensional incompressible flow

= (2.G)

(:an be evaluated using the non-equilibrium part of the dist.ribution function f_,,eq) .... [f<_ _ f}eq)]

( ) ( )

1 1 • eo/_ij (2.7)

O:

where e(_,i and e(_,j are ith and jth Cartesian component of t.he discrete velocity e_, respectively. For t.he

flow past a circular cylinder, a separate set of surmce points on the cylinder can be introduced in order

to carry out the numerical integration given by' Eq. (2.51). 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 secondorder extrapolation scheme based on the values of p and Tij at the neighboring fluid lattices. The forc_e

exerting on the boundary 0f_ is eomputed as

F' = f d..'t {pl + :u)+ , . (2.s)

•Jb_f2

It is worth commenting here that for the two-dimensional flow past, a cylinder, nearly half of the length of

the entire code is taken up by the above force evaluation procedure.
2.3. Method basedon the momentumexchange.In orderto employthemomentum--exchange

methodefficiently, twoscalararrays,w(i,j) and w_(i,j) are introduced. A value of 0 is assigned to w(i,j)

for the lattice site (i,j) tha.t are occupied by fluid; a value of I is assigned to w(i,j) for those lattice nodes

inside the solid t)ody. The array 'w5 (i, j) is set to zero everywhere except for those boundary nodes, Xb, where

a value of 1 is assigned. For a given nonzero velocity e,_, e(_ denotes "the velocity in "the opposite direction,

i.e., e_ - -e,_ (see Fig. 2). For a given boundary node xb inside the solid region with wv(i,j) ..... 1 and

w(i,j) = 1, the momentum exchange with all possible neighboring fluid nodes over a time step dt = 1 is

c_:/O

Simply summing the contribut;ion over all boundary nodes xv belonging to the body, the total force (acted

by the solid body on the fluid) is obtained as

all mb e_/0

In the momentum-exchange met;hod the force F is evaluated after the collision step is carried out; and

the value of j[_ at the boundary given by Eq. (2.1) has been evaluated. The monlentum exchange occurs

during the subsequent streaming step when f_,(x_,, t-t--St) and f_(x.f, t-+-5_) move to xf and xb, respectively.

As memioned in the introductory section, the efl>ct of the v_aria.ble A is not explicitly included, but it is

implicitly taken into account in the determination of f,,(xb, t + at). The applicability of Eq. (2.9) will be

examined and v_;lidated.

Clearly, the force is proportional to the number of boundary nodes xv in the above formula of F and

the number of the boundary nodes increase linearly with the size of the body in a two--dimensiona.1 flow.

However, since the force is normalized by pU2r in tile formula for CD in two--dimensions [see Eq. (3.9)], the

drug coefficient @ should be independent of r.

3. Results and Discussions. For straight walls, there is no doubt that Eq. (2.5) together with the

equation of state for pressure and Eq. (2.7) fbr r;j gives accurate results tbr the force provided that the f_'s

are accurately computed. To demonstrate the correctness of Eq. (2.9) based on the momentum exchange for

an arbitrary A, we first consider the pressure driven channel flow (see Fig. 4) for which exact solutions for

the velocity and stresses are known. The second case considered is the two--dimensional flow past a column of

circular cylinders at Reynolds number I_e = 100 and H/r = 20, where H is the distance between the centers

of two adjacent cylinders. The values of the drag computed using the two force evaluation methods are t;hen

compared with the result of ][:brnberg [38]. The dependence of the (]rag on the radius r in the momentumexchange method is examined to assess the reliability of this method. The third case is the two-dimensional

flow over a circular cylinder that is asymmetrically placed in a channel at Re = 100 (with vortex shedding).

The time dependence of the drag and lift coefficients is compared with results in the literat.ure.

We also consider two cases of three-dirnensional flow. The first case is the pressure driven flow in a

circular pipe for which the exact solutions for both the velocity profile and tile wa.ll shear stresses are known.

The assessment for the momentum-exchange method for three-dimensional flows will be made first in this

case. Firmlly, the momentum--excharlge method will be evaluated by considering the dra.g on a sphere due to

a uniform flow over a sphere in a finite domain. The details for t:he flow field computat;ion can be found in

Refs. [28, 11].

3.1. Two-dimensional pressure-driven channel flow. In the case of the channel flow, the force

on the top wall (y = H) at a given location x (i = N,_:/2--[--1, for example) can be evaluated using the
j N_j

j:=Ny .---]

X

j 2

j ===: ]

t_- L *t

\ ¢

P] ............................................................................................... _ / ) _2:t_......

-/ AS_

P2

FIG. 4. The channel flow co'_figu'rntion in the LBE simutations with an arbitrary A.

rnornentum--exchange method as follows. The wall is located between j = Ny and Nu - 1 (cf. Fig. 4). The

x and ;q cornt)onenls of the force oil the fluid at tile top wall near lhe ith node are

F_ = []6(i,J) + f_,(i - l,j - 1)] e6,.,: + [fs(i,J) + L(i + l,j - 1)] es# (3.1a)

/_), = []6(i,J) + f._,(i - l,j - 1)] e6,u + [fs(i,j) + A(i + l,j - 1)] cs,:_

+[]7(i,J) + L(i,j - 1)] or,y, (3.1b)

where c_,j denotes the jt.h Cartesian component of velocity e___.Since 5_. = 1, F,: and F,_ are, effectively, the

total shear and normal stresses, a_,j and u_y, which include the pressure and the deviatoric stresses, on the

fluid element at y = H.

Based on Eq. (2.7), the deviatoric component of the fluid shear stresses at j ....N_ - 1 (or y .....N_ -3 + A)

and Nv- 2 (or y .... N_ -4 + A) can be exactly evaluated based oil tile nonequilibriurn part of the dis_xibution

functions in the flow field if they are eorrectly given. A linear extrapolation of _uhe deviatoric shear stresses

to y = H = .N:_.....3 -+-2A yields

T(neq' %v(j ..... ,_7 _ i) + L_'_'k[%,u(j....N v - 1) - %y(j .... N - 2)] (3.2) _:g,W ...... • Y • Y

where the superscript "(neq)" denotes the value computed from f_,eq), and the subscript w refers to the

(neq)

value at the wall. The deviatoric normal stress, v,jy,_ , (!an be similarly computed. In a fully developed

channel flow, the normal component of the deviatoric stress V_y(y) is expected to be zero while the total

normal stress uw('y ) is equal to the negative of I;he pressure C-p). It needs "bo be pointed out that this

(neq) method of evaluat.ing r,_>.,, given by Eq. (3.2) for I;wo-dimensional channel flow is equivalent to "uhe method

of the surface stress integratk)n based on the ext;rapolat;ed pressure and the deviatoric stresses on I;he solid

wall except, that no numerical integration on the solid surface is needed.

After the velocity profile 'n_.(y) is obtained from f_, the shear stress %y on the wall can also be calculated

using the near wall velocity profile as

A

....pu.i-i--+----_i-['u_(j = I\\,, .....1) .....u;_,(j = N_ ....2)]. (3.3)

In the above, a linear extrapolation is employed to evaluate the velocity derivative (du_:/dy)lv=u at the wall.
'-£A_JLE3.1

Comparison of fluid stresses at y H in a two-dimensional pressure driven channel flow with dp/dx -1.0 × 10 -_ in

the lattice unit._, J\_y 35 and T 0.6 as a function of A. Column 2, 7-exact given by Eq. (3._); Cobumn 5', ---I_sj given ;cy _uJ '

by Eq. (3.13); Column _, (neq) -T£y,w given by b:q. (3.2); Column 5, -pz/(duj:/dy)iu=_t Eq. (3.3); Column 6, -b':t, given by

Eq. (3.1b); Cobumn _/, pressure p obtained in the simulation.

a ....._,_!,_× lo_' ....F,; × lo_ ...... ,-;'._ x lo" ....p,-_}?__ ,_ × lo" ....F, p

0.01 1.601 1.6333 1.6010 3.5294 O.3333 O.3333

0.02 1.602 1.6333 1.6020 2.5555 0.3333 0.3333

0.03 1.603 1.6333 1.6030 2.2309 0.3333 0.3333

0.04 1.604 1.6333 1.6040 2.0685 0.3333 0.3333

0.05 1.605 1.fi333 1.6050 1.9710 0.3333 0.3333

0.1 1.6:10 1.6333 1.6100 1.7760 0.3333 0.3333

0.2 1.620 1.6333 1.6200 1.6781 0.3333 0.3333

0.25 1.625 1.6333 1.6250 1.6583 0.3333 0.3333

0.3 1.630 1.6333 1.6300 1.6451 0.3333 0.3333

0.3333 1.633 1.6333 1,6330 1.6385 0.3333 0.3333

0.35 1.635 1.6333 1.6350 1.6357 0.3333 0.3333

0.4 1.640 1.6333 1.6400 1.6285 0.3333 0.3333

0.5 1.650 1.6333 1.6500 1.6184 0.3333 0.3333

0.6 1.660 1.6333 1.6600 1.6214 0.3333 0.3333

0.7 1.670 1.6333 1.6700 1.6244 0.3333 0.3333

0.8 1.680 1.6333 1.6800 1.6274 0.3333 0.3333

0.9 1.690 1.6333 1.6900 1.6305 0.3333 0.3333

0.95 1.695 1.6333 1.6950 1.6321 0.3333 0.3333

0.99 1.699 1.6333 1.6990 1.6335 0.3333 0.3333

Finally, the exact: solution for the fluid shear stress on the wall (y ..... I-/) is

._exact I eli, H. H = h_a .... 3 .+-2A (3.4)

based on the parabolic velocity profile or simple control volume analysis. This exact result call be used to

assess the accuracy of tile aforementioned methods for the force evaluation.

In the LBE simuladons_ the pressure gradient is enforced through t;he addition of an equivalent body

tbrce afl_er the collision step [26, 11]. While the velocity field given by the LBE solution can be unique,

the pressure field [thus the density field p(x,y)] can only be unique up to an arbitrary constant. In view

of Eq. (3.3), it is difficult to compare the stresses for different cases if p(i,j) converges to different values

in each case. To circumven_ this di_culty, the density field in the channel flow simulation is normalized by

p(i = 2,j = N;,ff2) at every time step. This normalization procedure results in p(x, y) = 1 throughout the

entire computational domain. It is also applied to tile three-dimensional flow in a circular pipe.

Table 3.1 compares t.he numerical values of the shear stress for a Wpical case (_"/Vy ..... 35, @/dx ...... 10 -_

_e×_,ct given by Eq. (3A), I;_ given by Eq. (3.13), r (neq} in the lattice units, and r .... 0.6) based on: 'xv,,, , <_,,,, given

by Eq. (3.2), and pt,(du_./dy)l:_:=H given by Eq. (3.3). Also listed is the comparison between I,} given by

Eq. (3.1b) and ----p. All computations are carried out with doub]e precision accuracy.

It is noted that. T (neq) is identical to ,_exact for all values of A. "\ closer examinat, ion of t,he shear

10
stress profile using Eq.. (2.7) across the channel reveals that; T_y,,_, (neq) __y), is also equal to the exact:, shear st;ress

profile _;_,:_ tYJ, which is linear, despit.e the errors in flle velocit;y profile 'a_(?j) for all values of A. A linear

extrapolation, Eq. (3.2), for _-_linear profile therefore gives t;he em-_et wall shear stress. Thus, _uhe exact;ness

oL z_,_, in the LBE simulation of cha,nnel flow indicates the reliability of the LBE solution tbr the stress

(neq] [ , _,_ _.[neq)

field v-,;j "kx,,9/ by using Eq. (2.7. However, as Fig. 3 indicates, the accuracy of integrating "7ij (a, y) to

obtain the fluid dynamic force in nontrivial geometries needs to be further investigated, as will be discussed

in the following sections.

For 0 < A < 1, the normal force _ given by Eq. (3.1b) based on the mornenturn-exehange method

agrees exactly with the pressure on the wall. This is a rather speeial quantity since deviatorie component of

the force is identically zero. Nevertheless, the method of _uhemomentum exchange does give a reliable value

for the normal stress.

For the shear (tangential) force, it is observed from Table 3.1 that for fixed @/dx, t_ does not etmnge

as A increases if'ore 0.01 to 0.99. On the other hand, the exac_ resuh _.ex_.c_,_y,_ = ?r(dp/dx)(A_. .....3 + 2A),

increases linearly with /A. Further computations were carried out over a range of Ny (= 35, 67, 99, and

131) and T (= 0.505, 0.51, 0.52, 0.6, 0.7, 0.8, 0.9, ].0, 1.2, 1.4, and 1.6). The result_s indicate that the

momentum-exchange method gives the shear stress on the top wall a.s

That is, F._ is independent of z and A. The error in E_ is zero when A = 1/3. The absolute error attains

the maximum when A = ], which gives the relative error of 4/3H for F;,_. Although the frequently used

momentum-exchange method is a natural choice for the three evaluation in conjunction with the bounce-back

boundary condition for A = 1/2, one must be aware that this method is not exaet arid the error in the force

evaluatk)n using the momentum-exchange method depends on A and the resolution.

The error in t<,. is due to the fact that "uhe derivatives of the velocity field are not; considered in the

boundary conditions. This can be understood by analyzing Eq. (3.1a). At the steady star:e, and with the

approximation that

f]_ _ f(eq).+ f_l) f(eq) ] 3 ...... ,_,,_, _(_. v)(,_,. _'} (3.6)

ct' " _ 9-" (_" J '

Equation (3.]a) at the top wall becomes

3

/%, _ 2u_,:_,.pTe2 - ('a_,s + u f - 2uw) , (3.7)

C

where the substitution of Eq. (2.4) for _,_ and J_ has been rnade. The only term ira the above equation which

has A dependence is ubf. When 0 <_:A 1/2, /_ is independent of A, and when 1/2 _< A < 1, F;_ weakly

depends on A because u_, = 0 in this ease [see Eqs. (2.3)]. In the case where F;,_is obtained by summing

over a set of symmetric lattice points, cancellations in the summation may further weaken the dependence

of F_: on A.

Table 3.1 also shows i:.hat for the shear stress based on the derivative of the velocity obt:ained by using

finite-difference, die loss of accuracy is q.ui'_e significant for small values of _ (< 0.05) when r .... 0.6. For

other" values of 2X (> 0.3), the accuracy is comparable wit:h t:hat of i;_. i-{owever, as shown in Fig. 5(a),

the aeeuracy of pp(du_:/dy)l _ n based on "_he near-wall velocity derivative deteriorates as the relaxation

time 7 increases (from 0.5] to 1.6). To see the cause of the increasing error in p,_(du_/d[q)l_,=rr , Fig. 5(b)

shows dirnensiordess wall velocity, u,,/u_, obtained by a three-point second-order Lagrangian extrapolation

] ]
0.006

0.004

0.002

(b)

h

0.2

1/a

0.5

0.7

FIG. 5. The LBE simulatioms of the ehannel flow, with A 0.2_ 1/3_ 0.5_ and 0.Z The pressure drop is O:_:p ....1.0 × 10 --(_

in lattice units. (a) l:_atio between the wall forcG puO_uxl_j=;_, evaluated by using Eq. (3.3), and the exact wdue ....%_y_wexact

-HO:ep/'2, 9iven by t':q. (3-I) as a function of T. (b) Normalized wall slip velocity u,w/uc as a function of "r.

of the near wall velocity u_(y) as a flmction of T. The increasing slip velocity u,_ on the wall with the

increasing relaxatk)n time v- was also observed in Ref. [15]. It is the result; of increasing partMe mean free

pat;h tha_ causes the deviation of _he kinetic solution from the hydrodynamk: solution. It is clear that _he

poor performance of pL,(d',_/dy)ly::::_ is associated with the increasing error in the near wall velocity profile

as r increases. Since the stress tensor ?;.i can be calculated directly from j'_ [see Eq. (2.7)i without the

need for directly completing velocity derivatives, the force evaluation method based on the evaluation of the

velocity gradient in _.he form of Eq. (2.6) is not recommended.

3.2. Steady uniform flow over a eohmm of cylinders. For a uniform flow over a column of circular

cylinders of radius r and center-to-center distance H (see the left part of Fig. 6 for illustration), symmetry

conditions for f,_'s are imposed at y ..... ±H/2. Most of the details of flow field simulation can be found in

Ref. [28]. The Reynolds number is defined by the diameter of the cylinder d as Re = Ud/lJ, where U is the

_miform velocity in the inlet. It must be noted that tbr a consistent, determination of the force, the upstream

boundary must be placed far upstream. A shorter distance between the cylinder and the boundary will result

in higher drag. In this study, it, is placed at about 20 radii to the left of the center of the cylinder. Re&_cing

the distance between the bo_mdary and the cylinder to 12.5 radii while keeping the rest of the computational

parameters fixed would increase the drag coefficient by about 1.8% at, Re - 100. The downstream boundary

is located about 25 30 radii behind the cylinder to allow sufficient wake development. The simulation is

terminated when the following criterion based on the relative L2-norm error in t;he fluid region ft is satisfied,

_: II'-(_, t + 1) - ..(x_, t)l?-

E2 = _-_-_--_) .......................................................................... < __. (3.s)

• _ ilu(x_,t + :1)11_ -

In this case, c ..... 10 -(_ was chosen h)r both Re .... 10 and 100.

Following Fornberg [38], the drag coefficient over a circular cylinder of radius r is defined as

i&[ (3.9)

CD ..... pU2r.

Figure 7(a) compares CD obtained from: momentum-exchange method, surface stress integration, and finite

difference result of Fornberg [38] using a vorticity-stream function formulation at Re = 100, H/r = 20, and

12
Z Uo

iii:: i .................................................................

Y

()

.................................... J

FIG. 6. Computational do'main for the uniform flow pa,_t a sphc'r_ o]' radius r, Th,e dashed lines indicate boundarie_s oj"

eo'mputatio'nal domain, (left) Unbounded domain i'n xg plane, and (right) bounded do'main in y_ plane.

radius r ranging from 2.8 to 13.2. t;br r > 8, both the momentum exchange and the stress integral;ion

methods give satisfactory result;s for @ in comparison with lhe value of 1.248 given in Ref. [38]. The

small differences in CD could l)e due to the fact that in Ref. [38], the computat;ional domain is much larger

in the downstream direction the downstream bolmdary condition is imposed at 300 radii behind the

cylinder in Ref. [38]_ as opposed to 25 30 radii here. This adds credence to the validity of Eq. (2.9) for

evaluating the total force on a body. The values of CD from the mornentum-exchange method have a. little

less variation than that from the stress integration. Accepting an error of less than 5%, reliable data for

CD can be ot)tained, using the momentum-exchange method, for r' > 5. That is, ten lattice spacings across

the diameter of the cylinder are necessary to obtain reliable values of the force. This is consistent with the

finding by Ladd [14]. In t:he range of 5 < r < 7, the stress-integration me_hod produces larger fluctuations

in the results than the momentum exchange method. ][:"or sm_,ller radius, i.e., coarser resolutions, while both

methods give poor results (due to insufficient resolution), the stress integration method yields much larger

errors.

Figllre ?(b) compares CD obtained from the methods of momentum exchange and the stress integration

for Re = :10. The momentum-exchange method seems to gives a converged result at larger r (> 8). Based on

the data for r > 8, an average value of CD _ 3.356 is obtained. In contrast, the stress--integration method has

a larger fluctuation than the large r result from the n_on_ent;um-exchange method even for r > 8. Averaging

over the resulis for r > 8, the stress integration gives CD _ 3.319. The difference between converged results

of two methods is about 1_. l%r r less than or around 5, the fluctuation in CD from _he stress-integration

method is much larger than thal in the momentum-exchange met;hod. The conclusions from the comparisons

in Fig. 7 are as follows: (i) both methods for force evaluation can give accurate results; (ii) the momentm_-

exchange method gives more consistent drag; and (iii) in the range of 10 < Re < 100, a resolution of ten

lattice spacings across the diameter of the cylinder are needed in order to obtain consistent and reliable

drag values. In other words, the lattice (grid) Reynolds number Re* (= U/_) should be less than 10 in the

calculations.

In the above results ]?resented in Figs. ?(a) and ?(b), the center of the cylinder is placed on a lat;tice

grid, thus t;he comput;ational mesh is symmetric with respect; t;o the geometry of the cylinder. To test the

effect of the mesh symmetry on the accuracy of t;he force evaluation, the calculation of the flow at Re ..... 10

is repeated with different vah_es of the cylinder center offset A_ in the z direction, or A_j in the ;q directioN.

The radius of the cylinder is deliberately chosen to be only 6.4 lattice grids. In order to preserve the mirror

13
O

1.8

1.7

1.6

1.5

1.4

1.3

1.2

1.1

2

........ R;=loo

X momentum exchange

[] [] stress inLegraLion

x

x D [] []

x

x D-xx 

[]

, , , i , , , i , , , i , , , i , , , i , , ,

4 6 8 10 12 14

T

0

3.51 (})) ..... Re=lO

3.41

3.31

3.21

X X momentum exehange

[] stress integration

[]

x

xN _ _-- _ --XD x x [] -_ []

..... =-_ .... [] [] [] []

[]

, , , i , , , i , , , i , , , i , , , i , , ,

4 6 8 10 12 14

F*C,. 7. The drag eoeJJicient for a uniform flow past a column of cylinders over a range of radius r. (a) [{e .... 100. The

dashed line indieates the value of CD 1.24: obtained i.n Ref. [38]. (b) Re 10. The da,shed lines indicate the values of CD

averaged over _ laryest radii.

symmet;ry of the flow in the y direction, we use different boundary conditions for upper and lower boundaries

(at y = ±H/2). For Ay = 0 while varying A_:, we use the symmetric boundary conditions, which maintain

the flow symmetry with respect to the center line in the x direction. For A_ = 0 while varying Av, we use

the periodic boundary conditions at y = ± H/2, which are equivalent to the symmetric boundary conditions

when A:_ = 0, but better reflect the flow symmetry when Ay ¢ 0. The results of the drag coefficient CD

are presented in Table 3.2. The variation of CD due to the change of the center of cylinder offset from a

grid poim is less than 19; when the cylinder diameter is only about 13 laltiee spacings. The outcome is

consistent with the expected truncation errors caused by mesh perturbation. We notice that the variation

in CD due to A;_, is about one order of magnitude smaller than that; due to Ay. This is precisely because

when Ay = 0 the mesh symmetry coincides with the flow symmetry in the y direction, and when Av¢ 0

the mesh symmet.ry is lost. This asymmetry due to A:_ _¢ 0 results in the change of the lift coefficient from

O(10 ---_4) to O(10--2), which is the same order of magnitude of _he variation in CD. It is our observation

that the accuracy of the force evaluation schemes used here is dictated by that of the boundary conditions

at the solid walls. The error due to symmetry of the computational mesh with respect to the geometry of

an object is well bounded. This is aJso observed in other independent studies [23, 33].

'raBL>: 3.2

The e_et of symmetry of the com.putational mesh on the force evaluation for the stead, y uniform flow over a column of

cylinders. The Reynolds number Re - 10 (:r - 0.6), the _ndius of the cylinder r - (i.4 (_n the lattice unit of &. - 1), and

H/r -- 20. The variation o]' (7o due to the change of the center of cylinder oJ'J_et j'rom a grid point is less than 17£

A, = 0, periodic boundary conditions at y = ±H/2

..................................................................................................................................................................

A:_ 0 0.2 0.4 I).6 0.8

cr, 3.3661 3.3637 3.3526 3.3526 3.3637

..................................................................................................................................................................

A v = 0, symmetric boundary conditions in y = ±H/2

..................................................................................................................................................................

Aa. 0 0.2 0.4 I).6 0.8

Co, 3.3661 3.3666 3.3646 3.3667 3.3692

..................................................................................................................................................................

It is worth noting t.hat the wall shea.r stress in the channel flow obtained by using the rnethod of

1,l
momentulnexchangehasarelativeerrorproportional totheresolutionacrossthechannel.Foraresolution

ofi10 20latticespacingsacrossthediameterconsideredhere,therelativeerrorinthedraga,ppears, however,

smallerthaninthechannel flowcase.At R.e= 100,withr > 10,theaveragevalueofthedragobtainedby

usingthemethodofmomentumexchangehasa1.7%relativeerrorcomparingwithFornberg'sdata;[38].If

theboundarylayerthicknessiestimatedroughlyto be3x 2r/_ _ 0, thereareonlyaboutsixlattice

spaeingsacrosstheboundarylayeroverwhichtheveloeidyprofilechanges substantially. Based on the insight

from the channel flow result, it is possible that the dev]atoric shear stresses on the surface of the cylinder

that are effectively incorporated in the method of momenturn exchange suf%r comparable levels of error as

in the channel flow. The eft%ctive error cancella:tion over the entire surface of the body may have contributed

to the good convergence behavior in the drag shown in Figs. 7(a) and 7(b).

3.3. Flow over an asymmetrically placed circular cylinder in channel with vortex shedding.

Sdl_%r and Turek [39] reported a set of benchmark results for a laminar flow over a circular cylinder of radius

r that is asymmetrically placed inside a channel. In tile present study, r = 12.8 is used and the center of

tlle cylinder coin(ides with a grid point. The distance fl'om the eenter of the cylinder to the upper wall

and lower wall is h+ .... 4.2r and h._ ----4.0r, respectively. This results in A+ ..... 0.76 for the upper wall and

A_ ----0.2 for the lower wall, respectively. The channel inlel has a parabolic profile and is placed at four

radii upstream of the cylinder center aecording to the specifieation of the benchmark test [39]. This results

in A = 0.2 tbr the inlet, boundary. A zeroth-order extrapolation tbr f_ is used at the exit boundary that is

located 40 radii downstream of the cylinder center. Thus there are a total of 564 × 105 square lattices in the

flow field. For Re = 2rU/v = 100 based on the average inlet velocity U, tile use of relaxation time T = 0.55

requires U _ 0.0fiS1.

At this Reynolds number, the flow becomes unsteady and periodic vortex shedding is observed. Figures 8(a), 8(b), and 8(c), respectively, show time-dependent behaviors of the li[% coemcient

CI: -----=,F_

t)U2r '

and the drag eoefl3cient CD [see Eq. (3.9)], and the pressure difference

A P = _'--t'---z-Pb

poU _ '

where p/and Pb are the pressures at the front and the back of the cylinder_ respectively, and P0 is the constant

density imposed at the entrance° The data of Co, CD, and AP are compared with the benchmark results in

Ref. [39]. We first note that the presei_t numerical value of Strouhal number St .... 2r/UT is about 0.3033,

where T is the period of the lift curve. This agrees very well with the range of St values (0.2950 0.3050)

given in Ref. [39]. We note that. the diffbrenee in Ct(t) between the mon_entum-exchange method and the

surface stress-integration method is indiscernible graphically. For _.he drag coefficient CD(t), it is interesting

to note that although there is about 0.25% difft, renee between the results given by the momentum--exchange

method a.nd the surfi_ce stress-integration method, both methods of foree evaluation give two peaks in the

Cr)(t) curves. Physically, these two peaks in the @ (t) curve correspond to the existence of a weaker vortex

and a stronger vortex alternately shed behind the cylinder. The difference in the strength of the vortices

results from the difference: h+/r ..... 4.2 and h.__./r ..... 4°0 in the passages between the cylinder and the channel

walls. There is no report on the occurrence of these two peaks in Ref. [39]. Instead, a range of the maxinmm

CD (from 3.22 to 3.24) by different researchers was given. The present value of the higher peak is well within

the range. It is interesting to note that, both peaks of CD(t) obtained by the momentum-exchange method

15
3.25

-1 3.17

0 0

i

(b) ...... t..... hange

---- stress integration

................................

.....................

I

,\J \7 , V , 7 -_

1000 2000 3000

t-t o

2.46 ..................................

, , , , , , , , , 1 , , , , , , , , , 1 , , , , , , , , ,

0 1000 2000 3000

t-t o

F_(:. 8. The 2D flow past a cylinder asymmetrically placed in a channel. The variations of the lij_ coe]ficient C'L, the

drag eoeJficient Up, and the pressure difference AP as functions of time _ (a]_er an initial run time to) are compared with

the benchmark results in t_ej'. /39]. At the time to, the lift. eoe]jieient 6'j,(t) attains its maximum value C)I _×. The dashed

horizontal lines indicate the upper and lower bounds in f_ef. [39]. The solid and dashed curves are the _esults obtained by using

momentum exchange and stress integratian, respectively. (a) The lij_ coej]icient CL (t). Note that the result._' obtained by using

the two methods are indistinguishable on the graph. (b) The drag coejficient 6:D(t). (e) The pressure difference AP(t). The

symbol >,', indicates the value of AP(_o + T/2) given i'n Table 3.3, where T (_ 1296.5) is the period of CL(t).

are also within the range, as shown in Fig. S(b). A fl_rther refined computation of the present problem using

a multiblock procedure [40] with r = 40 in the fine grid region yield nearly the same results for Cry(t) and

cL(t).

We compile in Table 3.3 the values of Strouhal number St, maximum and minimum drag coefficient

C_ _ and C_Dn'in , maximum and minimum lift coefficient (-'_'_ _, and C-m-Jin L , and the pressure diffi?rence AP

obtained by the LBE methods and other schemes of computational fluid dynamics given in Ref. [39]. The

value of AP is mea.sured at to + T/2, where to is the moment when CL(t) reaches its maximum value C_ 'a_,

and T is the periodicity of CL (t). For the LBE simulations, T is between ] 296 and 1297 (in the lattice unit

of 5t .....1). We use T ----1296.5 in the determination of the Strouhal number St. With a resolution much

coarser than those used in Ref. [39], the LBE results are well within the bounds given in Ref. [39]. This

clearly demonstrates the accuracy of the lattice Boltzm_nn method.

3.4. Pressure-driven flow in a circular pipe. The ste_-_dy-state flow field was obt.ained by using

D3Q]9 model with T ----0.52 [11]. Eq. (2.9) is used to evaluate the force on the boundary points along the

circumt_rence of the pipe over a distance of one lattice iu the axial direction. The resulting axial ibrce F_

16
TAmm 3.3

Vo.lue8 _max (?mln _m_x '-_mln of St_ C D , "D , C L , C L , amt AP ]'or the flow over a 2D cylinder asymmetrically placed in a channel

"Momentum" and "Stress" denote, respectivdy, the momentum-exchange method and the stress-integration method in the LBE

caleulations. The CFD results" are the bounds in [_e[., [39], which does not have data for 6D-_m_and ('L-'H_I_'"

method S 1; C}_ a'x C}_ in C}_,'ax C_ 'i_' AP

Mornentum 0.3033 3.2358 3. ] 7711 1.0045 -1.0347 2.49114

Stress 0.3033 3.2275 3.1708 1.0040 - 1.0340 2.4914

CFD 0.2950 0.3050 3.2200 3.2400 0.9900 1.0100 2.4600 2.5000

is, equivalently, the force given by _ 2rrrS_, where _,, is the wall shear stress and r is the pipe radius. For

a fully developed flow inside a circular pipe, the exact fluid shear stress at, the pipe wall is given by

I w , 57) = 7l-7 .2 diI: "

We examine t,he normalized axial force,

_J..... d_ (13.11)

.w_.2 d_

Figure 9 shows t_he normalized coefficient r_ over a range of r: 3.5 23.5. Except for r _ 5, _l is rather

close to 1. It was notieed in Ref. [11] that the accuracy of LBE solution for the pipe flow is not as good

as that fbr the two-dimensional channel flow due to the distribution of values of A around the pipe. The

accuracy of the drag is dict;ated by the accuracy of the flow field if the force evaluation method is exaet.

For the pipe flow, the error in t;'_, results from the inaeeuracy in the flow field and the errors in the force

evaluation scheme based on momentum exchange (as seen in the previous section for the two-dimensional

channel flow case). For r > 5, the largest error in F_. is about 3.5% and it occurs at r = 15.5. Again,

there is no systematic error in F;,:. Given the corr@exity of the boundary in this three-dimensional flow, _.he

results shown in Fig. 9 a.re satisfa.etory in the sense that it adds fllrther credence to the momentum-exchange

1.0

0.9

method ff)r ff}rce evaluation.

g-

[]

0 5 10 15 20 25

r

FIG. 9. The ratio '7 betwee*_ the tange'ntial foree [;_: on the pipe and its exact value (;vr '2Vp) over a ra'nge of pipe radius r.

17
3.5. Steadyuniform flow overa sphere.Tolimit thecomputational effort,a finitedomainof

-H/2 < y 51 H/2 and -H/2 _<_"z < H/2, with H/_* = 10 is used to compute the flow past a sphere of

radius r' (see Fig. 6). Two eases are considered: (a) the flow past a single sphere, and (b) the flow over a

two-dimensional array of spheres (all located at :r = 0) with the (:enter of the spheres forming square lattices.

In the former ease, the boundary eondi{;ions at jv ..... 1 (:{/ ..... ft/2 corresponds to j_ .... 2) for f,_'s are given by

the following linear extrapolation

The veloeit;y at jv ..... 2 is set as

(3._2)

u(j_, 2,j_) = u(j_., 3,j_). (3.13)

Similar treatment is applied at y = H/2 and z = ±H/2. In the latter ease, symmetry eonditions are posed

on fee's at Jv = 1 by using the values of f_'s at J:v = 3 (see Ref. [28] for the two--dimensional ease). At the

inlet, a uniform velocity profile is imposed at j_, ----1.5 (half way between the first and second lattices). The

upstream boundary is located at 7.5 radii to the left of the sphere center in all simulations.

D)r flow over a sphere, the drag eoefFmien{; is often expressed as

Fr 24 F,

C'D ........ f--77.2-------_= w--- _ _b ............................ (3.14)

_pt,-rc_ _-' t_e ' ' 6"JrrUpu '

where 0 accounts for the non-Stokesian effect of the drag. t,br two "hypes of t;he boundary conditions at

(_] .... ±H/2 and z ..... ±H/2), 0_ denot;es the non-Stokesian correction for t;he case where the symmetry

conditions are imposed at (?j ..... ±H/2 and z .... ±H/2) and .oo_ denotes the results for t:.he case where the

extrapolation for f_ is used at (y = :i:H/2 and z = ±H/2) in order to simulate the unbounded flow.

Figure 10(a) shows the non-Stokesian coefficient 0_ for r = 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 5.1, 5.2, 5.4, 5.6,

and 5.8, for H/r = l0 a.t Re = 10. The relaxation time is r = 0.7. With this range of r, the nurnber of

tile boundary nodes on the surface of the sphere increases roughly by a factor of (5.8/3) 2 _ 3.743 the actual

counts of the boundary nodes m!, gives a ratio 2370/546 = 4.3..5. The largest, diffbrenee is 1.98_ between r = 3.0

and r = 3.2 that have the least resohJtion in the eases investigated. For a uniform flow over an unbounded

sphere, an independent eomput;at;ion using a finite difference met;hod based on t;he vorticit:y-stream funct:ion

formulation wit;h high resolution gives a drag coefficient _5 _ 1.7986 at Re ..... 10. The largest difference

between this result and the LBE result;s is 1.36% at r ..... 3.2. If the LBE data for the drag is averaged over

the range of r, one obtains r/5 _-. 1.8086, which differs f'rorn 1.7986 by 0.54%. Hence, the LBE solutions with

3.0 < r < 5.8 yield very consistent values for the drag force. Figure 10(b) shows the non-Stokesian correction

factor (5_ tbr a uniform flow over a planar array of spheres for 3.1) < r <: 5.8 and H/r" = 10, at Re = 10. It

is important to note that with the improvement of the surface resolution by a factor of 4.3,5, there is little

systematic variation in c),(r). The largest deviation from the average value, _ _ 1.963, is 1.1_. at r = S.0. It

is (:lear that; the LBE solution gives reliable fluid dynamic forces on a sphere at r _ 3.5 for a moderate value

of lie. The set of data for cO, is inherently more consistent than that for _.<_ since the symmetry boundary

condition can be exactly specified at ._1..... 2t-f/2 and z ..... ±H/2, while the extrapolation conditions given by

Eqs. (3.12) and (3.13) do not guarantee the free s{;ream condition at; '.y .... ±H/2 and z ..... ±H/2. Yet, both

_ and _ exhibit remarkable self-consistency t?'om coarse to not-so-coarse resolutions.

4. Conclusions. Two methods for evaluating the fluid force in conjmletion with {;he method of la{;tiee

Boltzmann equation for solving fluid flows involving curve(] geometry have been examined. The momentun>

exchange method is very simple t.o implement. It. is shown in the channel flow simulation that momentumexchange method is not an exact method. The error in the wall shear stress is inversely proportional to

18
-o1.8g

1.82

1.81

1.80

1.79

1.78

(a)

[]

[]

[] []

i

Re=lO HJr=m

[]

[]

[]

[]

[]

[]

2.00

"_ 1.95

(b)n ' R LIo H/ =IO

[]

[]

[] []

[] []

[]

..... 0----

[]

[]

1.90 []

I,,,,I,,,,I,,,, I,,,,I,,,,I,,

3 4 5 3 4 5

T T

F_C. 10. Flow past sphere. Variation of the non-Stokesian correction ]_ctor _ -b_,/OTrrUpt, as a function of sphere

• (

radius r at Re 10. The dashed li'_.es are rabies o] O,r) averaged over" r. (a) The flow past a sin.gle sphere in an unbounded

field (If/r .--- oc). (b) The ]tow 'past a planar arr'ay of spheres (lt/r" ---. 10).

the resolution. In two- and three-dimensional flows over a bluff body, it; (:an give accurate drag values

when there are at least ten lattice spacings across ¢;he body at Re _ 100. The met;hod of integrating lhe

stresses on the surface of the body gives similar results when there is sufficient resolution but. it exhibits

much larger fluctuations {.ban that in the method of momentum exchange when the resolution is limit, ed. In

addition, the stress-integrat, ion method requires considerably more efforts in implementing the extrapolation

and int.egration on the body surface in comparison with the method of momentum exchange.

It is interest;ing to note that t:he momen_tum-exehange method is perhaps superior t;o the stress-integration

method because the former method is directly based on the distribut, ion f_nctions while the latter is derived

from f_rther processiug of t,he distribution functions. In addition, the momentmn-exchange method _lses

interpolations while the stress-integration method uses exr,rapolations. Often extrapolations a,re more noisy

and unstable than interpolations. Even with a coarse resolution that does not yield very aiccurate local

information, accurate force evaluation can be accomplished with the lattiee Boltzrnann method. Among the

two force evahmt;ion methods, the method of moment;urn exchange is recommended for force evaluation on

curved boundaries for it;s simplicity, accuracy, and robustness.

Acknowledgments. The authors are grat;eful to Prof. A.J.C. Ladd for his helpfid comments. This

work is partially supported by NASA Langley l:{esearch Center under the program of Innovative Algorit;hms

for Aerospace Engineering Analysis and Optimization. R. Mei also acknowledges partial support of lhe

Engineering Research Center (ERC) for Particle Science and Techno]ogy at the University of Florida, the

National Science Foundation (EEC-9402989), and Industrial partners of the ERC.

Appendix A. LBE Models in Two and Three Dimensions.

The nine-velocity (or 9-bit) LBE model on a _wo-dimensional square lattice, denoted as t;he D2Q9 model,

has been widely used for simulations of two-dimensional flows. For three-dimensional flows, there are several

cubic lattice models, such as _he fifteen-velocity (D3Q15), ninet;een-veloeity (D3Q19), and twent;y-sevenvelocity (D3Q27) models, which have been used in the lit;erature [10]. All these models have a rest particle

(with zero velocity) in the discretized velocity set {e_l-(_ = 0, 1, ..., (b ....._)}. For athermal fluids, the

equilibrium distributions for l.he D2Q9, D3Q15, D3Q19, and D3Q27 models are all of the following form

19
[6, 7]

f(fq)=vGp I + _(e_'u)+_ =c" j

where w_ is a weighting factor and eo is a discrete velocity, c - 5_/5, is the unit speed, a,nd 5:_, and at are

the lat;tice constant and the time step, respectively. The discrete velocities for t;he D2Q9 models are

{ (0, 0), ,_ ..... 0,

e<_= (±l, 0)c, (0, =:1)c, a .....1, a, 5, 7, (A.2)

(±1, ±1) c, (_ = 2, 4, 6, 8,

and the va,lues of the weighting factor wo are

g_ (t = O;

,,_ = }, _ ....1, a, a, 7, (f.a)

1

_, _ .... 2, 4_ 6, 8.

For the D3QI9 model, the discrete velocities are

(0, 0), c_ = 0,

_o, = (±i,o, o)< (o, ±1, o)(:, (0, 0, ±_)_, _ = :1< (A.4)

(±1, ±1, O) c (0 ±1, ±l)c, (±1, O, ±])G c_=7 18,

and the weighting factor w_ is given by [7]

wo = , (_ ..... 1 0, (A.5)

a_ (_..... 7 18.

The discrete velocity set;s {G_} for the D2Q9 and D3Q19 models are shown in Fig, 1.

The density and velocity ('.an be comput:ed from f,-,,

, .....Zs,, .....EI,7 '),

c_, (_

= }2 = }2 <,s2 <A.sb)

The speed of sound of _he above LBE models is

1

C s _ _C_

and the equation of sl;ate is that of an ideal gas such that

p .... c7 p. (A.7)

The viscosity of the fluid is

1/ ..... C s

for l;he discrete velocity model of Eq. (1.2). It should be noted thai; the equilibrium distribut;ion flmction

f:_eq) is in fact a Taylor series expansion of the Maxwellian f(0) [6, 7']. This approximation of .f(_q) in algebraic

form makes the LBE method valid only in the incompressible flow limit t_/c-----> 0.

2o
Equation(1.2)isoftendiscretizedinspacex andtimet into tile lattice Boltzrnann equation

1 (

fa:(xi + ea5_, t + at) - f,_(xi, t) = -_[f,:,a;i, t) - f},eq)(xi, t)] (A.8)

where r .... ),/at. Fbr this LBGK model [1, 2], t:he viscosity in "the Navier-Sl;okes equation derived from the

above lattice Boltzmann equation is

u= _...... c_O_. (i.9)

The -1/2 correction in the above formula for p comes from the second--order derivatives of f_, when f_(xi +

e_ at, t + at) in Eq. (A.8) is expanded in a Taylor series in _,. This correction in l/makes the lattice Boltzm an n

method formally a second-order method for solving incompressible flows [7]. Obviously, the physical and

computational stabilit;ies require that; r > 1/2.

REFEREN CES

[1] H. CttEN, S. CHEN, AND W.H. MATTttAEUS, Recovery of the Navier-Stohes equations using a latticegas Boltzmar_,n method, Phys. Rev. A, 4,5 (1992), pp. R5339 R5342.

[2] Y.H. Q_AN, D. D'HUMI_;a>;S, AND P. LALLEMAND, Lattice BGK models for Navier-Stokes W'_tation,

Europhys. Lett., 17 (1992), pp. 479 484.

[3] R. BFNZI, S. SUCCI, AN[) M. VERGASSOLA, The lattice t3ottzmann eq'_atior_,: Theory a_d applicatiorzs,

Phys. Rep., 222 (1992), pp. 145 197.

[4] S. CHEN AND G.D. DOOLEN, Lattice Boltz.ma_r_ method .[or.fl'_,id flows, Annu. t_ev. Fhfid Mech., 30

(1998), pp. 329 364.

[5] P.L. Bttz\TNA(-_AR, E.P. GROSS, AN13 M. KROOK, A model .for" collision processes iv, gases, l. Small

arwolit_Me processes i_ char:qed ar_,d ,_eutral o'ac-compor_ent system, Phys. Rev., 94 (1954), pp. 511

525.

[(3] X. HE AND L.--S. LUO, ./][pr'iori derivation of the lattice Boltzrn, a'nr_ eq,tation, Phys. R_ev. E, 55 (1997),

pp. R(3333 R633(3.

[7] X. HE AND L.--S. LUO, Theory of the lattice BoItzman_t method: From, the [3oltz'ma_rt equation to the

lattice Boltz'man,_ eq_,atior_., Phys. Rev. E, 5(3 (1997), pp. 6811 6817.

[8] T. ABE, Derivatior_. of the lattice Boltzmann method b,v meaus of the discrete ordinate m.ethod for the

Boltzmanv, equatior_., J. Comput. Phys., 131 (1997), pp. 241 24(3.

[9] X. SIIAN AND X. HE, Discretizatiov. of the velocity space i_ the soh, tion @ the Boltzro.ann cquatior_,

Phys. Rev. Lett., 80 (11998), pp. 65 (38.

[10] Y.H. QIAN, S. SuCCI, AND S.A. ORSZAC4, Recer_t advar_ces in lattice Bottz__,a_r_ co'm,p'_,tir_.9, in A_,rtv.al

Reviews of Corrtputatiortal Physics, Vol. Ill, edited by D. Staufli?r, (World Scientific, Singapore,

199,5), pp. 195 242.

[11] l:l. MEL W. StlYY, D. YU, AND L.-S. LUO, Lattice Boltzruamt method for 3-d flows with curved

bou_dar:q, J. Comput;. Phys., 161 (2000), pp. 680 699.

[12] D.P. Zm(_LER., Bo_u'udar"q co_Mitious for lattice Bolt._rnanv. simulations, J. S_;at. Phys., 71 (1993),

pp. 1171 1177.

[13] I. GINZBOURO AND P.M. ALDErt, Bov.,_dary flo_v conditior_, a_.al'qsis for the thrce-diruensiorzal lattice

Boltzrnanr_, model, J. Phys. II, 4 (1994), pp. 191 214.

21
[141A.J.C.LAJOD, Nu'merical simulation of particular suspensions via a discrctized BoItz'mann equation,

Part 2, N_merical results, J. Fluid Mech., 271 (1994), pp. 3111 339.

[:151 D.R. NOBLE, S. CHEN, J.G. GEO-H(_[ADIS, AN[) R.(). }:_UCKI[US, A consistent hydrodynamic boundary

condition for the lattice lloltzmann method, Phys. Fluid, 7 (1995), pp. 203 209.

[161 O. BEHREND, Solid boundaries in particle suspension simulations via lattice Bottzmann 'method, Phys.

Rev. E, 52 (1995), t)l;). 1164 1175.

[171 T. INAMUf{O, M. YOSH[NO, AND F. OG[NO, A non-slip boundary condition for lattice Boltzmann

simulations, Phys. Fluids, 7 (1995), pp. 2928 3000.

[181 S. C,HF, N, D. MAf{TfNF, Z, AND R. MEI, On boundary conditions in lattice Boltzmann method, Phys.

Fluids, 8 (1996), pp. 2527 2536.

[191 I. GINZBURG AND D. D'HuMI_;RI-:S, Local second-order boundary methods for lattice [)oltzmann m, odels,

J. Slat. Phys., 84 (1996), pp. 927 971.

[201 X. HE, Q. Zou, L.-S. Luo, AND M. DEMBO, Analytic solutions and analysis on non-slip boundary

condition for the lattice Boltzmann [dGK model, J. Star. Phys., 87 (1997), pp. 115 136.

[211 L.-S. Luo, Analytic solutions of linearized lattice t3oltzmann equation for simple flows, J. St;at. Phys.,

88 (1997), pp. 913 926.

[22] Q. zou AND X. HE, On pressure and velocity boundar'y conditions for the lattice Boltzm, ann BCK

model, Phys. Fluids, 9 (1997), pp. 1591 1598.

[23] M. BOUZTDI, M. FIaDAOUSS, AND P. LALLF, MAND, Momentum transfer of a Boltzmann-lattice fluid

with boundaries, Phys. Fluids, 13 (2001), pp. 4,52 3459.

[24] A.J.C. LADD AND R. VEaB,:aC, Lattice Boltzmann simulations of particle-fluid suspensions, J. Stat.

Phys., 104 (2001), pp. 1191 1251.

[251 1:I. ME[ AND W. SHYY, On the finite d_fference-based lattice Boltzmann method in curvilinear coordinates, J. Coil:lpuC. Phys., 143 (1998), pp. 426 448.

[261 X. HE AN[) G. DOOLEN, Lattice Boltzmann method on curviIinear coordinates syste'm,s: Flow around a

ci'_rular cylinder, J. Comput. Phys., 134 (1997), pp. 306 315.

[271 O. F_LWPOVA AND D. H;*iNEL, Grid refincro.ent.for lattice-BCg 'models, J. Comput. Phys., 147 (1998),

pp. 219 228.

[281 R. MEI, L..-.S. LUO, AND W. SHYY, An accurate curved boundary treatment in the lattice Boltzmann

method, J. Comput. Phys., 155 (1999), pp. 307 330.

[291 M.N. KOGa_'4, Rarefied Gas Dynamics (Plenum, New 5%rk, 1969).

[301 S. HAaRIS, An Introduction to the _lzeor'y of the Boltzrnann Equation (Holt, Rinehart and Winston,

New York, 1971).

[311 C. CEeC_CNA_q, The Boltzrnann Equation and its Applications (Springer-Verlag, New York, 1988).

[32] R. CORNUBERT, D. D'HuMIgaES, AND D. LF_VERMOaE, A Knudsen layer theory for lattice gases,

Physica. D, 47 (1991), pp. 241 259.

[331 D. D'HuM_gRES, M. BOUZIDI, AND P. LALL>;MANI), Thirteen--velocity 3D lattice Boltzrnann model,

Phys. Rev. E, 63 (2001), 066702.

[341 O. FI[LIPPOVA (private communication).

[351 W. SHYY, Computational Modeling for Fluid Flow and htte@tcial Transport, revised printing (Elsevier,

Amsterdam, 1994).

[361 W. SHYY, S.S. THAKUR., H. OUYANG, J. LIU, AND E. BLOSCU, Computational _I_chniqv.es for Complex Transport Phenomena (Cambridge University Press, Cambridge, 1997).

22
[37]S.THAKUH, W. SHYY_H.S.UDAYKUMAR_ ANDL. HILL,Mu[tibtock intc_faec t'reatmer_ts in a pressurebased flow solver, Numer. Heat Tr_msfer, Pa.rt B, 33 (]998), pp. 367 396.

[38] B. FORNBERG, Steady incompressible flow past a row of cir(mlar cytir_der, J. Fluid Mech., 225 (1991),

pp. 655 671.

[39] M. SCH;{FEa AND S. TUREK, Be.vchmark cornp,_tations of laminar.flow over a c:qlinder, in Ho'w Simulation with High-t-'erforrnanee Cornpvters H, edited by E.H. tfirschd, Notes in Numerical Fluid

Meehar_ies (Vieweg, Braunschweig, 1996), \%1. 52, pp. 547 566.

[40] D. Yu, R. MF. L aND W. SHYY, A multi-block lattice Bottzrna.nn .method fur fluid ]tows, to appear in

Int. J. Numer. Methods Fhfids (2002).

23
REPORT DOCUMENTATION PAGE Form Approved OMB No. 0704-0188

Publicreportingburdenforthiscollectionof informationisestimatedto average] hourperresponse, includhTgthetimeforreviewinginstructions,searchingexistingdatasources,

gatheringandmaintainlngthedataneeded,andcompletingandreviewingthecollectionof informatlon.Sendcommentsregardingthisburdenestimateoranyotheraspectof this

collectionofinformation,includingsuggestionsforreducingthisburden,to WashingtonHeadquartersServices,DirectorateforInformationOperationsandReports,1215Jefferson

DavisHighway,Suite1204,Arlington,VA222024302,andto theOfficeof Management andBudget,PaperworkReductionProject(07040188),Washington, DC20503.

1. AGENCY USE ONLY(Leaveblank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

July 2002 Contractor Report

4. TITLE AND SUBTITLE

FORCE EVALUATION IN THE LATTICE BOLTZMANN METHOD

INVOLVING CURVED GEOMETRY

6. AUTHOR(S)

Renwei Mei, Dazhi Yu, Wei Shyy, and Li-Shi Luo

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

ICASE

Mail Stop 132C

NASA Langley Research Center

Hampton, VA 23681-2199

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Langley Research Center

Hampton, VA 23681-2199

5. FUNDING NUMBERS

C NAS1-97046

WU 505-90-52-01

8. PERFORMING ORGANIZATION

REPORT NUMBER

ICASE Report No. 2002-22

10. SPONSORING/MONITORING

AGENCY REPORT NUMBER

NASA/CR-2002-211662

ICASE Report No. 2002-22

11. SUPPLEMENTARY NOTES

Langley Technical Monitor: Dennis M. Bushnell

Final Report

To appear in Physical Review E.

12a. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified Unlimited

Subject Category 34

Distribution: Nonstandard

Availability: NASA-CASI (301) 621-0390

12b. DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

The present work investigates two approaches for force evaluation in the lattice Boltzmann equation: the momentumexchange method and the stress-integration method on the surface of a body. The boundary condition for the particle

distribution functions on curved geometries is handled with second order accuracy based on our recent works. The

stress-integration method is computationally laborious for two-dimensional flows and in general difficult to implement

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 uniform 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.

14. SUBJECT TERMS

lattice Boltzmann method, force evaluation on fluid-solid interface, momentum-exchange

method, stress-integration method, boundary condition for curved geometries, accuracy,

3-D flows

15. NUMBER OF PAGES

28

16. PRICE CODE

A03

17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION

OF REPORT OF THIS PAGE OF ABSTRACT OF ABSTRACT

Unclassified Unclassified

NSN 7540-01-280-5500 StandardForm 298(Rev. 2-89)

PrescribedbyANSIStd. Z3918

298102