Finite-Difference Lattice Boltzmann Methods for Binary Fluids (Mathematical Aspects of Pattern Formation in Complex Fluids)

193

Finite-Difference Lattice Boltzmann Methods for

Binary

Fluids

Aiguo Xu

Department

_of

Physics, Yoshida-South Campus,

Kyoto University, SakyO-ku, Kyoto, 606-8501, Japan

In this proceedingwe summarize

our

recent studiesontwofluid lattice Boltzmann methods for

binary fluids. We first clearify Sirovich’s kinetic theory, then based on which three multispeed

discrete velocity modelsareformulated, whicharefor the Euler equations,isothermal Navier-Stokes

equations andthe completeNavier-Stokesequations, respectively. Eachformulateddiscretevelocity model, together with an appropriate finite-difference scheme, composes afinite-difference lattice

Boltzmann method. Thevalidity of the methods is verified by investigating (i)the Couette flow

and(ii) theuniform relaxation process of the two components.

PACSnumbers: $47.11.+\mathrm{j}$,$51.10.+\mathrm{y}$,$05.20.\mathrm{D}\mathrm{d}$

I. INTRODUCTION

Lattice Boltzmann Method (LBM) isanumericalscheme to simulate kinetic systems. The

LBM

recovers

the hydrodynamic descriptions in the small Knudsen number limit. It has become

aviableand promisingnumericalscheme for simulatingfluidflows. Thereareseveral options to discretizethe Boltzmann equation: (i) Standard LBM (SLBM)[I]; (ii)

Finite-Difference

LBM (LBM) [1-3]; (iii) Finite-Volume LBM$[1, 4]$;(iv)finite ElementLBM[I, 5]; etc. These kinds

of schemes

are

expectedtobe complementary in the LBMstudies.

Even though various LBMsformulticomponentfluids 8-20] have been proposedand

devel-oped, (i)mostexistingmethodsbelong totheSLBM[6,8-17, $\mathrm{a}\mathrm{n}\mathrm{d}/\mathrm{o}\mathrm{r}$basedonthe single-fluid

theory[8-15, 17, 18, 21]; (ii) inRefs. $[6, 7]$ twoSLBMsareproposed, but these two modelsare

not convenient (ifnotimpossible) tosimulatethermaland compressible systems,

even

isother-maland incompressible systems only if the two components have different particle

masses.

In this study

we

develop two fluid FDLBMs forthermal and compressible binary fluids.

II. FORMULATIONANDVERIFICATION OF THEFDLBMS

The formulation of aFDLBM consists of three steps: (i) select

or

design

an

appropriate

discretevelocity model(DVM), (ii) formulatethediscrete local equilibrium distribution

func-tion, (iii) choose afinite-difference scheme. The continuous Boltzmann equation has infinite

velocities,sotherotationalinvariance is automaticallysatisfied. Recoveringrotational invariant

macroscopic equationsfrom adiscrete finitevelocity microscopicdynamicsimposesconstraints

on

the isotropy of DVM used. In

our

studies, the proposed FDLBMs

are

based

on

the two DVMsdescribed below.

DVMI$:\mathrm{v}_{0}=0$, $\mathrm{v}_{\mathrm{k}1}=\mathrm{v}_{\mathrm{k}}[\cos(\frac{\mathrm{i}\pi}{6}),$$\sin(\frac{\mathrm{i}\pi}{6})],\mathrm{i}=1,2$,\cdots ,12, (1)

where$k$indicates the$k$-thgroup ofparticle velocities and$i$indicatesthe directionof theparticle

speed. It iseasyfindthat (i) its odd rank tensors

are

zero, and(ii) its initial four

even

rank tensors satisfy

$\sum_{0=1}^{12}v_{ki\alpha}v_{ki\beta}=6v_{k}^{2}\delta_{\alpha\beta}$, $\sum_{i=1}^{12}v_{ki}\alpha kv\dot{\iota}\beta v_{ki\gamma}v_{ki\delta}$ $= \frac{3}{2}v_{k}^{4}\Delta_{\alpha\beta\gamma\delta}$,

$\sum_{i=1}^{12}vk\dot{\iota}\alpha vki\beta v_{k:}v\cdot v_{ki\mu}v_{k\nu}\gamma b\delta|.=\not\supset^{v_{k}\Delta_{\alpha\beta\gamma\delta\mu\nu}}16$ , (2) $\sum_{\dot{l}=1}^{12}v_{k}|\alpha vk_{\dot{l}}\beta vki\gamma vk_{\dot{l}}\delta v_{k\dot{\cdot}\mu}v_{k:\nu}v_{k\lambda}|.v_{ki\pi}=\frac{1}{32}v^{8}k\Delta_{\alpha\beta\gamma\delta\mu\nu\lambda\pi}$,

where$\alpha$, $\beta$, $\cdots$ indicate$x$or$y$componentand

$\Delta_{\alpha\beta\gamma\delta}=\delta_{\alpha\beta}\delta_{\gamma\delta}+\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma}$ , (3) 数理解析研究所講究録 1413 巻 2005 年 193-200

$\Delta_{\alpha\beta\gamma\delta\mu\nu}=\delta_{\alpha\beta}\Delta_{\gamma\delta\mu\nu}+\delta_{\alpha\gamma}\Delta_{\beta\delta\mu\nu}+\delta_{\alpha\delta}\Delta_{\beta\gamma\mu\nu}+\delta_{\alpha\mu}\Delta_{\beta\gamma\delta\nu}+\delta_{\alpha\nu}\Delta_{\beta\gamma\delta\mu}$, (4)

$\Delta_{\alpha\beta\gamma\delta\mu\nu\lambda\pi}=\delta_{\alpha\beta}\Delta_{\gamma\delta\mu\nu\lambda\pi}+\delta_{\alpha\gamma}\Delta_{\beta \mathit{5}\mu\nu\lambda\pi}+\delta_{\alpha\delta}\Delta_{\beta\gamma\mu\nu\lambda\pi}+\delta_{\alpha\mu}\Delta_{\beta\gamma\delta\nu\lambda\pi}$

$+\delta_{\alpha\nu}\Delta_{\beta\gamma\delta\mu\lambda\pi}+\delta_{a\lambda}\Delta_{\beta\gamma\delta\mu\nu\pi}+\delta_{\alpha\pi}\Delta_{\beta\gamma\delta\mu\nu\lambda}$

.

(5)

It is clear that this DVM isisotropic up to,

at

least, its 9th rank tensor.

(6)

DVM2:$\mathrm{v}_{0}=0$,$\mathrm{v}_{\mathrm{k}1}=\mathrm{v}_{\mathrm{k}}[\mathrm{c}\mathrm{o}\mathrm{e}$ $( \frac{\mathrm{i}\pi}{4})$,$\sin(\frac{\mathrm{i}\pi}{4})],\mathrm{i}=1,2$,$\cdots,8$

.

Similarly, (i) its odd ranktensors

are

zero, and (ii) its initial

three even

ranktensorssatisfy

$\sum_{i=1}^{12}v_{k\iota\alpha}v_{k_{l}\beta}=4v_{k}^{2}\delta_{\alpha\beta}$, $\sum_{i=1}^{12}v_{ki\alpha}v_{ki\beta}v_{ki\gamma}v_{k\delta}|.=v_{k}^{4}\Delta_{\alpha\beta\gamma\delta}$,

(7)

$\sum_{\dot{l}=1}^{12}v_{k\dot{\iota}\alpha}v_{kt\beta}v_{ki\gamma}v_{ki\delta}v_{ki\mu k_{t}\nu}v$ $= \frac{1}{6}v_{k}^{6}\Delta_{\alpha\beta\gamma\delta\mu\nu}$

.

DVM 2 is isotropic up to its7thrank tensor.

We considerabinarymixturewithtwocomponents,$A$and$B$, where the

masses

and

temper-aturesof the two componentsare notsignificantlydifferent. The interparticlecollisions canbe divided into two kinds: collisions within the

same

species (self-collision) and collisionsamong differentspecies (cross-collision) [22]. Based

on

theDVM(1),the 2-dimensionalBGK[23]kinetic equationforspecies$A$reads,

$\partial_{t}f_{ki}^{A}+\mathrm{v}_{ki}^{A}\cdot\frac{\partial}{\partial \mathrm{r}}f_{ki}^{A}-\mathrm{a}^{A}\cdot\frac{(\mathrm{v}_{ki}^{A}-\mathrm{u}^{A})}{\Theta^{A}}f_{k}^{A(0)}|.=J_{k\iota}^{AA}+J_{k}^{AB}|$. (8)

where

$J_{k\dot{l}}^{AA}=-[f_{ki}^{A}-f_{ki}^{A(0)}]/\tau^{AA}$ , $J_{ki}^{AB}=-[f_{ki}^{A}-f_{ki}^{AB(0)}]/\tau^{AB}$ (9)

$f_{ki}^{A(0)}= \frac{n^{A}}{2\pi\Theta^{A}}\exp[-\frac{(\mathrm{v}_{kj}^{A}-\mathrm{u}^{A})^{2}}{2\Theta^{A}}]$, $f_{ki}^{AB(0)}= \frac{n^{A}}{2\pi\Theta^{AB}}\exp[-\frac{(\mathrm{v}_{k\dot{\mathrm{a}}}^{A}-\mathrm{u}^{AB})^{2}}{2\Theta^{AB}}]$ (10)

$\Theta^{A}=k_{B}T^{A}/m^{A}$, $\Theta^{AB}=k_{B}T^{AB}/m^{A}$ (11)

$f^{A(0)}$ and $f^{AB(0)}$

are

the correspondingMaxwellian distribution functions. $n^{A}$, $\mathrm{u}^{A}$, $T^{A}$

are

the

local

density, hydrodynamic velocity andtemperature ofspecies

A.

$\mathrm{u}^{AB}$, $T^{AB}$

are

the

hydrodynamic velocity and temperature of the mixture after equilibrationprocess. $\mathrm{a}^{A}$ is the

accelerationof species $A$due tothe effectiveexternalfield.

Forspecies$A$,we have

$n^{A}= \sum_{k_{\dot{1}}}$ $f_{ki}^{A}$ ,

$n^{A} \mathrm{u}^{A}=\sum_{ki}\mathrm{v}_{ki}^{A}f_{ki}^{A}$, $P^{A}(e_{\mathrm{i}\mathrm{n}\mathrm{t}}^{A}=n^{A}k_{B}T^{A})= \sum_{ki}\frac{1}{2}m^{A}(\mathrm{v}_{k}^{A}|.-\mathrm{u}^{A})^{2}f_{ki}^{A}$ (12)

where$P^{A}(e_{1\mathrm{n}\mathrm{t}}^{A})$isthelocal pressure(internal energy). For species$B$,

we

have similarrelations.

For themixture,we have

$\mathrm{u}^{AB}=(\rho^{A}\mathrm{u}^{A}+\rho^{B}\mathrm{u}^{B})/\rho$, _{$nk_{B}T^{AB}= \sum_{kt}\frac{1}{2}[(\mathrm{v}_{k\mathrm{s}}^{A}-\mathrm{u}^{AB})^{2}m^{A}f_{ki}^{A}+(\mathrm{v}_{kv}^{B}-\mathrm{u}^{BA})^{2}m^{B}f_{ki}^{B}]$}

(13) where$\rho^{A}=n^{A}m^{A}$, $n=n^{A}+n^{B}$_and $\rho=\rho^{A}+\rho^{B}$

.

Three sets of hydrodynamic quantities (forthetwo components$A$,$B$and forthe mixture)

are

involved, but only two sets of them

are

195

independent. So this is a tw0-fluid model. Without lossing generality, we focus on hydrody-namics of the two individualspecies. By expanding the local equilibrium distributionfunction $f^{AB(0)}$ around$f^{A(0)}$ to thefirstorderinflowvelocity and temperature, theBGK model(8-11)

becomes

$\partial_{t}f_{ki}^{A}+\mathrm{v}_{ki}^{A}\cdot\frac{\partial}{\partial \mathrm{r}}f_{ki}^{A}-\mathrm{a}^{A}\cdot\frac{(\mathrm{v}_{k\mathrm{n}}^{A}-\mathrm{u}^{A})}{\Theta^{A}}f_{kt}^{A(0)}=Q_{ki}^{AA}+Q_{ki}^{AB}$ (14)

$Q_{k\iota}^{AA}=-( \frac{1}{\tau^{AA}}+\frac{1}{\tau^{AB}})[f_{ki}^{A}-f_{ki}^{A(0)}]$ (15)

$Q_{k}^{AB}|.=- \frac{f_{ki}^{A(0)}}{\rho^{A}\Theta^{A}}\{\mu_{D}^{A}(\mathrm{v}_{ki}^{A}-\mathrm{u}^{A})\cdot(\mathrm{u}^{A}-\mathrm{u}^{B})$

$+ \mu_{T}^{A}[\frac{(\mathrm{v}_{kt}^{A}-\mathrm{u}^{A})^{2}}{2\Theta^{A}}-1](T^{A}-T^{B})-M^{A}[\frac{(\mathrm{v}_{k_{1}}^{A}\cdot-\mathrm{u}^{A})^{2}}{2\Theta^{A}}-1](\mathrm{u}^{A}-\mathrm{u}^{B})^{2}\}(16)$

where$\mu_{D}^{A}=\rho^{A}\rho^{B}/(\tau^{AB}\rho)$, $\mu_{T}^{A}=k_{B}n^{A}n^{B}/(\tau^{AB}n)$, $M^{A}=n^{A}\rho^{A}\rho^{B}/(2\tau^{AB}n\rho)$

.

Now,we go to the secondstep: formulate $f_{ki}^{A(0)}$. ThecontinuousMaxwellian$f^{A(0)}$ possesses

an infinite sequence of moment properties. The Chapman-Enskog analysis[24] shows that, requiring the discrete $f_{ki}^{A(0)}$ to follow the initial eight

ones

is sufficient to describe the

same

Navier-Stokesequations,

$\frac{\partial\rho^{A}}{\partial t}+\frac{\partial}{\partial r_{\alpha}}(\rho^{A}u_{\alpha}^{A})=0$, (17)

$\frac{\partial}{\partial t}(\rho^{A}u_{\alpha}^{A})+\frac{\partial}{\partial \mathrm{r}_{\beta}}(\rho^{A}u_{\alpha}^{A}u_{\beta}^{A})+\frac{\partial P^{A}}{\partial r_{\alpha}}-\rho^{A}a_{\alpha}^{A}-\frac{\partial}{\partial r_{\beta}}[\eta^{A}(\frac{\partial u^{A}}{\partial r_{\beta}}+\frac{\partial u_{\beta}^{A}}{\partial r_{a}}-\frac{\partial u}{\partial r}\mathrm{L}\delta_{\alpha\beta})\gamma A]$

$+ \frac{\rho^{A}\rho^{B}}{\tau^{AB}\rho}(u_{\alpha}^{A}-u_{\alpha}^{B})=0$, (18)

$- \partial \mathrm{e}_{\frac{A}{t}+\frac{\partial}{\partial r_{\alpha}}}\not\supset[(e^{A}+P^{A})u_{\alpha}^{A}]-\rho^{A}\mathrm{a}^{A}\cdot \mathrm{u}^{A}-\frac{\partial}{Tr_{\alpha}^{-}}[k^{A}\frac{\partial(k_{B}T^{A})}{\partial r_{\alpha}}+\eta^{A}u_{\beta}^{A}(^{\partial}*_{\beta}^{u^{A}}+\frac{\partial \mathrm{u}_{\beta}^{A}}{\partial r_{\alpha}}-\frac{\partial u_{\gamma}^{A}}{\partial r_{\gamma}}\delta_{\alpha\beta})]$

$+_{\tau}^{\epsilon^{A}}*_{\rho}^{B}[(u^{A})^{2}-\mathrm{u}^{A}\cdot \mathrm{u}^{B}]+_{\overline{\tau}n}n_{\mathrm{I}^{n}arrow k_{B}(T^{A}-T^{B})-n^{A}\frac{\rho^{A}\rho^{B}}{2\tau^{AB}n\rho}(\mathrm{u}^{A}-\mathrm{u}^{B})^{2}=0}^{AB}$, (19)

where

$e^{A}=e_{\mathrm{i}\mathrm{n}\mathrm{t}}^{A}+ \frac{1}{2}\rho^{A}(u^{A})^{2}$, $\eta^{A}=P^{A}\tau^{AA}\tau^{AB}/(\tau^{AA}+\tau^{AB})$, $k^{A}=2n^{A}\Theta^{A}\tau^{AA}\tau^{AB}/(\tau^{AA}+\tau^{AB})$

.

(20) Recall that $\mathrm{u}^{A}(\mathrm{u}^{B})$ is a small quantity. By usingEq. (17), $P^{A}=n^{A}k_{B}T^{A}$, and neglecting

thesecondand higherorderterms in $\mathrm{u}^{A}$,Eq. (18) showsthat thediffusionvelocity,$u_{\alpha}^{B}-u_{\alpha}^{A}$,

isrelatedto the gradients of$n^{A}$ and$T^{A}$

.

Thefirst threerequirements

on

$f_{ki}^{A(0)}$ are referredto Eq. (12) with$f_{k}^{A}|$.replacedby$f_{k\dot{l}}^{A(0)}$,and

theremainingfive

are

$\sum_{ki}m^{A}v_{ki\alpha}^{A}v_{kv\beta}^{A}f_{ki}^{A(0)}=P^{A}\delta_{\alpha\beta}+\rho^{A}u_{\alpha}^{A}u_{\beta}^{A}$ (21)

$\sum_{ki}\frac{1}{2}m^{A}(v_{ki}^{A})^{2}v_{ki\alpha}^{A}f_{ki}^{A(0)}=2n^{A}k_{B}T^{A}u_{\alpha}^{A}+\frac{1}{2}\rho^{A}(u^{A})^{2}u_{\alpha}^{A}$ (23)

$\sum_{kj}\frac{1}{2}m^{A}(v_{ki}^{A})^{2}v_{ki\alpha}^{A}v_{ki\beta}^{A}f_{ki}^{A(0)}=2P^{A}\Theta^{A}\delta_{\alpha\beta}+\frac{1}{2}P^{A}(u^{A})^{2}\delta_{\alpha\beta}$

$+3P^{A}u_{\alpha}^{A}u_{\beta}^{A}+ \frac{1}{2}\rho^{A}(u^{A})^{2}u_{\alpha}^{A}u_{\beta}^{A}$ (24)

$\sum_{k_{1}}\frac{1}{2}m^{A}(v_{ki}^{A})^{4}v_{ki\alpha}^{A}f_{ki}^{A(0)}=[12P^{A}\Theta^{A}+6P^{A}(u^{A})^{2}+\frac{1}{2}\rho^{A}(u^{A})^{4}]u_{\alpha}^{A}$ (25)

The requirement equation (25) contains the fifth order ofthe flow velocity $\mathrm{u}^{A}$

.

So

it is sufficienttoexpand$f_{k_{\dot{l}}}^{A(0)}$in polynomial uptothefifth orderof$\mathrm{u}^{A}$:

$f_{k\dot{\iota}}^{A(0)}=n^{A}F_{k}^{A} \{[1-\frac{(u^{A})^{2}}{2\Theta^{A}}+\frac{(u^{A})^{4}}{8(\Theta^{A})^{2}}]+\frac{v_{k_{2}\xi}^{A}u_{\xi}^{A}}{\Theta^{A}}[1-\frac{(u^{A})^{2}}{2\Theta^{A}}+\frac{(u^{A})^{4}}{8(\Theta^{A})^{2}}]$

$+ \frac{v_{ki\xi}^{A}v_{ki\pi}^{A}u_{\xi}^{A}u_{\pi}^{A}}{2(\Theta^{A})^{2}}[1-\frac{(u^{A})^{2}}{2\Theta^{A}}]+\frac{v_{k\not\in}^{A}v_{\mathrm{k}t\pi}^{A}v_{k\iota\eta}^{A}u_{\xi}^{A}u_{\pi}^{A}u_{\eta}^{A}}{6(9^{A})^{3}}[1-\frac{(u^{A})^{2}}{2\Theta^{A}}]$

$+ \frac{v_{ki\xi}^{A}v_{k\iota\pi}^{A}v_{ki\eta}^{A}v_{ki\lambda}^{A}u_{\zeta}^{A}u_{\pi}^{A}u_{\eta}^{A}u_{\lambda}^{A}}{24(\Theta^{A})^{4}}+\frac{v_{ki\epsilon^{v_{ki\pi}v_{ki\eta}v_{ki\lambda}v_{ki\delta}u_{\xi}u_{\pi}u_{\eta}u_{\lambda}u_{\delta}}}^{AAAAAAAAAA}}{120(\Theta^{A})^{5}}\}$

$+\cdots$ (26)

where

$F_{k}^{A}= \frac{1}{2\pi 8^{A}}\exp[-\frac{(v_{k}^{A})^{2}}{2\Theta^{A}}]$

.

(27)

Thetruncated equilibriumdistribution function$f_{ki}^{A(0)}(26)$ containsthefifthranktensorofthe

particlevelocity$\mathrm{v}^{A}$

and therequirement (22) contains its third rank tensor. Thus,a DVMbeing

isotropicup toits8th rank tensors is enoughtorecoverthephysical isotropyofthe continuous

Boltzmann equations to the Navier-Stokes level. So DVM (1) is an appropriate choice. To calculate the discrete $f_{ki}^{A(0)}$,

one

firstneeds calculate the factor $F_{k}^{A}$

.

$F_{k}^{A}$ isdeterminedbythe

eight requirements

on

$f_{ki}^{A(0)}$andthe isotropic propertiesof the DVM (1). We finally obtain

$\sum_{k\dot{\iota}}F_{k}^{A}=1$, $\sum_{k}F_{k}^{A}(v_{k}^{A})^{2}=\frac{\mathrm{e}^{A}}{6}$, $\sum_{k}F_{k}^{A}(v_{k}^{A})^{4}=\frac{2}{3}(\mathrm{e}^{A})^{2}$

$\sum_{k}F_{k}^{A}(v_{k}^{A})^{6}=4(\Theta^{A})^{3}$, $\sum_{k}F_{k}^{A}(v_{k}^{A})^{8}=32(\Theta^{A})^{4}$, $\sum_{k}F_{k}^{A}(v_{k}^{A})^{10}=320(\Theta^{A})^{5}$ (28)

Once

a

zerospeed, $v_{0}^{A}=0$, and other fivenonzeroones,$v_{k}^{A}(k=1,2,3,4,5)$ arechosen,$F_{k}^{A}$

$(k=0,1,2,3,4,5)$ will be fixed.

Wecome to the third step: finite-differenceimplementationof the discrete kinetic method. Thereare morethanonechoices$[2, 18]$ available. Onepossibility is shown below,

$f_{ki}^{A,(n+1)}=f_{ki}^{A,(n)}+[ \mathrm{a}^{A}\cdot\frac{(\mathrm{v}_{ki}^{A}-\mathrm{u}^{A})}{9^{A}}f_{ki}^{A(0)}+Q_{k\dot{\iota}}^{AA,(n)}+Q_{ki}^{AB,(n)}-\mathrm{v}_{ki}^{A}\cdot\frac{\partial f_{k\dot{l}}^{A,(n)}}{\partial \mathrm{r}}]\Delta t$ , (29)

wherethesecondsuperscripts$n$,$n+1$indicate the consecutive two iteration steps,$\Delta t$the time

step; thespatialderivatives

are

calculated$\mathrm{a}\epsilon$

$\frac{\partial f_{ki}^{A,(n)}}{\partial\alpha}=\{$

$(3f_{b,I}^{A,(n)}-4f_{ki,I-1}^{A,(n)}+f_{k\iota,t-2}^{A,(n)})/(2\Delta\alpha)$ if$v_{k\iota\alpha}^{A}\geq 0$

$(3f_{k_{\dot{l}},I}^{A,(n)}-4f_{ki,I+1}^{A,(n)}+f_{k\dot{\mathrm{r}},t+2}^{A,(n)})/(-2\Delta\alpha)$ if$v_{k\mathrm{i}\alpha}^{A}<0$

137

where$\alpha$$=x$,$y$, thethirdsubscripts $I-2$,$I-1$, $I$,$I+1$,$I+2$indicate consecutive mesh nodes

in the$\alpha$direction.

If the kinetic numerical scheme is required to recover the hydrodynamics only up to the isothermal Navier-Stokeslevel, Eqs. (17)-(18)orthe Euler level, Eqs.(17)-(19) with$\eta^{A}=k^{A}=$

$0$, followingthe sameprocedures, it is easy to find thatDVM 2 is enough. For the isothermal

Navier-Stokesequation, Eq. (28) is replaced by

$\sum_{ki}F_{k}^{A}=1$, $\sum_{k}F_{k}^{A}(v_{k}^{A})^{2}=\frac{\mathrm{e}^{A}}{4}$, $\sum_{k}F_{k}^{A}(v_{k}^{A})^{4}=(\mathrm{e}^{A})^{2}$

$\sum_{k}F_{k}^{A}(v_{k}^{A})^{6}=6(\mathrm{e}^{A})^{3}$ (31)

Forthe complete Euler equation, Eq. (28) is replaced by

$\sum_{ki}F_{k}^{A}=1$, $\sum_{k}F_{k}^{A}(v_{k}^{A})^{2}=\frac{\mathrm{e}^{A}}{4}$, $\sum_{k}F_{k}^{A}(v_{k}^{A})^{4}=(\mathrm{e}^{A})^{2}$

$\sum_{k}F_{k}^{A}(v_{k}^{A})^{6}=6(\Theta^{A})^{3}$. $\sum_{k}F_{k}^{A}(v_{k}^{A})^{8}=48(\Theta^{A})^{4}$ (32)

In summary, to

recover

the tw0-dimensional complete Navier-Stokes equations, a 2-dimensional 61 velocity $(\mathrm{D}2\mathrm{V}61)$ model is needed; a $\mathrm{D}2\mathrm{V}33$ model is sufficient to

recover

thetw0-dimensionalEuler equations; recovering the

tw0-dimensional

isothermal Navier-Stokes equationscan resort on asimpler$\mathrm{D}2\mathrm{V}25$model. In principle,

a

DVM withlower isotropy

can

be

replaced byonewith higherisotropy. But in practicalsimulations, onegenerallyneeds choose the simplest

one.

The validity of the formulated the FDLBMs is verified through two test examples. (The

Boltzmannconstant $k_{B}=1.$) The firstone is theisothermal and incompressibleCouetteflow

with a single component. In this case, $A=B$. The initial state of the fluid is static. The

distance between the two walls is $D$

.

At time $t=0$ they startto moveat velocities $U$, $-U$,

respectively. It is clear that all the three models $(\mathrm{D}2\mathrm{V}25, \mathrm{D}2\mathrm{V}33, \mathrm{D}2\mathrm{V}61)$ work for such a

system. The horizontalvelocity profilesof species $A$

or

$B$along

a

vertical lineagreewith the

following analytical solution,

u$= \gamma y-\sum_{j}(-1)^{j+1}\frac{\gamma D}{j\pi}\exp(-\frac{4j^{2}\pi^{2}\eta}{\rho D^{2}}t)\sin(\frac{2j\pi}{D}y)$, (33)

where$\gamma=2U/D$isthe imposed the shearrate,$j$isaninteger,the two walk locate at$y=\pm D/2$

.

(Forexample,seeFig. 1.)

The second

one

isthe uniform relaxation process, which is anideal processto indicate the equilibration behavior of the mixture. By neglecting the force terms and terms in spatial

derivatives, theNavier-Stokesequations (17)-(19) give

$\frac{\partial}{\partial t}\rho^{A}=0$, (34)

$\frac{\partial}{\partial t}(\mathrm{u}^{B}-\mathrm{u}^{A})=-\frac{1}{\rho}(\frac{\rho^{A}}{\tau^{BA}}+\frac{\rho^{B}}{\tau^{AB}})(\mathrm{u}^{B}-\mathrm{u}^{A})$ , (35)

$\frac{\partial(T^{B}-T^{A})}{\partial t}=-\frac{1}{n}(\frac{n^{A}}{\tau^{BA}}+\frac{n^{B}}{\tau^{AB}})(T^{B}-T^{A})+\frac{\rho^{A}\rho^{B}}{2k_{B}n\rho}(\frac{1}{\tau^{AB}}-\frac{1}{\tau^{BA}})(\mathrm{u}^{B}-\mathrm{u}^{A})^{2}$ (36)

(TheEulerequationsplay thesameroleasthe Navier-Stokesequationsinthis

case.

Both the

$\mathrm{D}2\mathrm{V}61$and$\mathrm{D}2\mathrm{V}33$works. ) Theflow velocities of the two componentsequilibrateexponentially

with time. (Forexample, seeFig. $2(\mathrm{a}).$) The equilibration offlowvelocities also affects that

of the temperatures. When the flow velocity difference is zero, the temperatures equilibrate exponentially with time. (For example, seeFig. $2(\mathrm{b}).$) The simulation resultsagreewell with

$\supset\geq$

$\mathrm{y}/\mathrm{D}$

FIG. 1: Horizontalvelocity profilesalongaverticalline for the twospecies, $A$and$B$, at time_$t=8$

.

Thesymbolsareforsimulationresults. Thesolid line correspondstothetheoreticalresult, Eq. (33).

Parameters used in the tw0-fluid FDLBM are $m^{A}=m^{B}=1$, $T=1$, $n^{A}=n^{B}=1$, $\gamma=0.001$,

$\tau^{AA}=\tau^{BB}=\tau^{AB}=\tau^{BA}=0.2$

.

_{Parameters used}$1\mathrm{n}|$Eq. (33) are$\eta=\eta^{A}=0.1$,_{$\rho=\rho^{A}=1$}

.

$\mathrm{t}$

$\mathrm{t}$

FIG. 2: Uniform relaxation processes, (a) Equilibration ofvelocities; (b) Equilibration of

tempera-tures. Thesymbolsare for simulation results. The solid lines possess the theoreticalslopes. Common

parametersfor the simulations in(a)and (b)are$n^{A}=10$_,$n^{B}=1$_,$m^{A}=1$_,$m^{B}=10,7^{-AA}=\tau^{BB}=1$,

$\tau^{AB}=10$_,$\tau^{BA}=1$

.

$\ln(\mathrm{a})$ the initial conditionsare $u_{x}^{A(0)}=-u_{x}^{B(0)}=-0.3$, $u_{y}^{A(0)}=u_{y}^{B(0)}=0$, and

$T^{A(0)}=1.3$_, $T^{B(0)}=0.7$

.

_{Theslopeofthe solid linein (b) is -11/20, which is consistent with Eq.}

(35). In (b) theinitial conditionsare$\mathrm{u}^{A(0)}=\mathrm{u}^{B(0)}=0$, and$T^{A(0)}=1.3$_,$T^{B(0)}=0.7$

.

_{The slopeof}

the solid line in (b) is -10.1/11, whichis consistentwiththe firstterm ofright-handsideofEq.(36).

The second superscript “(0)” denotesthe corresponding initial value. Thisfigure showsan example

where theparticlemassesof thetwo speciesaresignificantlydifferent.

III. CONCLUSIONS AND REMARKS

The Chapman-Enskog analysis shows what properties the discrete Maxwellian distribution function $f_{ki}^{A(0)}$ should follow. Those requirements tell the lowest order ofthe flow velocity

$\mathrm{u}^{A}$ in

the Taylor expansionof$f_{k\dot{v}}^{A(0)}$. The highest rank oftensors ofthe particlevelocity $\mathrm{v}^{A}$

intherequirements

on

the truncated $f_{k\dot{\iota}}^{A(0)}$ determines the needed isotropyof the DVM. The

incorporationof theforce temsmakesno additionalrequirement

on

the isotropy of the DVM.

199

interfacial tension is to modify the pressure tensors[14], which is implemented by changing the force terms[3]. The specific force terms or pressure tensors depend on the system under

consideration,whichareoutof the scope of thisLetter,but canbe resolved under thesame

tw0-dimensional61-velocitymodel$(\mathrm{D}2\mathrm{V}61)$

.

For binary fluids with disparate-mass components, say

$m^{A}\ll m^{B}$, only if the total

masses

and temperatures ofthe two species

are

not significantly

different, Sirovich’s kinetic theoryworks, so do the correspondingFDLBMs. (See Fig. 2 for

anexample.) When the masses$\mathrm{a}\mathrm{n}\mathrm{d}/\mathrm{o}\mathrm{r}$the temperatures ofthe two components are greatly

different, the tw0-fluid kinetic theory should bemodified. Inthose cases, the Navier-Stokes equationsand the FDLBMs

are

notsymmetric about the two components, butthe FDLBMs can still beresolved under the$\mathrm{D}2\mathrm{V}61$ model. The formulation procedure is straightforward.

A

more

detailed description is referred to $[25, 26]$

.

We finally emphasize that the numerical

errors

ffom the finite-difference schemes result in artificial viscosities in the simulation. The

comparisonofvarious finite-differenceschemes,discussion

on

numerical accuracy and stability are referred to Refe. [2, 3, 18].

Acknowledgments

Aiguo Xu acknowledgesProf. G. Gonnella for guiding him into the LBM field and thanks Profe. H. Hayakawa, V. Sofonea, S. Succi for helpful discussions. This work is partially sup-ported byGrant-in-Aids for Scientific Research(GrandNo. 15540393)and for the 21-th Century COE “Centerfor Diversity and Universality in Physics” ffom the MinstryofEducation,Culture and Sports,Scienceand Technology (MEXT) ofJapan.

[1] SucciS., The Lattice Boltzmann Equation, OxfordUniversity Press, New York, 2001; Benzi R.,

Succi S.Vergassola M., Phys. Rep. 222,3(1992);HigueraF.,Succi.S. BenziR.,Europhys.Lett.

9, (1989) 345.

[2] CaoN.,Chen S.,Jin S. MartinezD., Phys. Rev.$\mathrm{E}55$, (1997)R2l; SetaT.Takahashi R. , J. Stat.

Phys. 107, (2002) 557; Sofonea V. Sekerka R. F., J. Comp. Phys. 184, (2003) 422; Watari M.

TsutaharaM., Phys. Rev.$\mathrm{E}67$, (2003) 36306;KataokaT.TsutaharaM., Phys.Rev.$\mathrm{E}69$, (2004)

56702; KataokaT. TsutaharaM., Phys.Rev. $\mathrm{E}69$, (2004)R35701.

[3] SofoneaV.,Lamura A., GonnellaG.,CristeaA., Phys.Rev.$\mathrm{E}70$, (2004)046702.

[4] NannelliF.,andSucciS.,J. Stat.Phys.68 ,(1992) 401; Xi H., PengG.Chou S.H., Phys.Rev. $\mathrm{E}$

60,(1999) 3380;Ubertini S.,Bella G., SucciS.,Phys.Rev.$\mathrm{E}68$, (2003) 16701;Seealsoreferences in [1].

[5] LiY.,LeBoeufe E. BasuP.K., Phys. Rev. $\mathrm{E}69$,(2004) 065701(R).

[6] LuoL.S. Girimaji S. S., Phys.Rev. $\mathrm{E}66$,(2002) R35301; ibid, 67, (2003) 36302$[]$

[7] Malevanets A. and YeomansJ., Faraday Discuss.112, (1999) 237.

[8] GuoZ. Zhao T. S., Phys. Rev.$\mathrm{E}68$, (2003) R35302.

[9] Xu Aiguo, Gonnella G. Lamura A., Phys. Rev. $\mathrm{E}67$, (2003) 56105; Physica A331, (2004) 10;

PhysciaA 344, (2004) 750;c0nd-mat/0404205; Xu Aiguo,Commun.Theor. Phys.39, (2003)729.

[10] Gunstensen A. K. , Rothman D.H.,ZaleskiS.,and Zanetti G.,Phys.Rev. A43, (1991) 4320.

[11] FlekkoyE. G., Phys.Rev.$\mathrm{E}47$, (1993) 4247.

[12] GrunauD.,ChenS.,EggertK.,Phys. Fluids A 5, (1993) 2557.

[13] Shan X. Doolen G., J. Stat. Phys. 81, (1995)379; Phys.Rev.$\mathrm{E}54$, (1996)3614.

[14] Orlandini E., Osborn W. R., Yeomans J. M., Europhys. Ltt. 32, (1995)463; Osborn W. R.,

OrlandiniE.,Swift M.R.,Yeomans J.M.,Banavar J. $\mathrm{R}$,Phys.Rev.Lett., 75, (1995) 4031; Swift

M. R., Orlandini E., Osborn W. R., Yeomans J. M., Phys. Rev. $\mathrm{E}$, 54, (1996) 5041; Gonnella

G., OrlandiniE.,Yeomans J.M.,Phys.Rev. Lett. 78, (1997)1695, Phys.Rev. $\mathrm{E}58$, (1998) 480;

Lamura A.,GonnellaG., YeomansJ.M.,Europhys. Lett.45, (1999)314.

[15] Kendon V.M., Desplat J. C.,Bladon P., Cates M. E., Phys.Rev. Lett., 83, (1999) 576;Kendon

V. M.,CatesM.E.,PagonabarrageI.,DesplatJ.$\mathrm{C}$,Bladon P.,J. Fluid Mech. 440, (2001) 147. [16] YuH.,LuoL. S.,Girimaji S.S.,Int. J. Comp. Eng.Sci.3, (2002)73.

[17] Facin P.C., PhilippiP. $\mathrm{C}$,dosSantos L. O. E. inICCS2003, LNCS2657

editedbySloot P.M.

A. etal., 2003,pp.1007-1014.

[18] CristeaA.SofoneaV., CentralEuropeanJ. Phys. 2,(2004)382.

[19] ShmX.ChenH., Phys.Rev.$\mathrm{E}47$, (1993)1815; ibid,49, (1994)2941.

[20] Love P. J. Coveney P.V., Phys. Rev. $\mathrm{E}64$, (2001)21503; Love P. $\mathrm{J}.,\mathrm{M}\mathrm{a}\mathrm{i}\mathrm{l}\mathrm{l}\mathrm{e}\mathrm{t}$J. B., Coveney P. V., ibid, 64, (2001)61302; Chin J. Coveney P. V., ibid, 66, (2002) 16303; Gonzalez-Segredo $\mathrm{N}$,

NekoveeM., Coveney P.V., ibid,67, (2003)46304; Gonzalez-Segredo N. Coveney P. V., ibid, 69,

(2004)61501.

[21] Sofonea V.SekerkaR. F., PhysicaA299, (2001)494.

[22] SirovichL., Phys.Fluids5,(1962) 908; ibid,9, (1966)2323.

[23] BhatnagarP.L.,GrossE.P., KrookM., Phys. Rev., 94, (1954)511.

[24] Chapmann S. and Cowling T. G., The Mathematical Theory

_of

_Non-Unifom

Gases, 3rd ed.,

CambridgeUniversityPress, Cambridge1970.

[25] Xu Aiguo, Phys. Rev. $\mathrm{E}$(in press) [ $h\mathrm{f}\mathrm{f}\mathrm{l}p.\cdot//amiv$

.

0rg/abs/c0nd-mat/04060lB].