布朗运动的粒子的尺寸一般都在微米量级以下,它受到四周流体分子的碰撞才不平衡，颗粒的运动才显著

来源: marketreflections 于 2010-11-18 12:36:13 [] [博客] [旧帖] [给我悄悄话] 阅读数 : (17605 bytes)
回答: 保守力场力 牛顿三定律 由 marketreflections 于 2010-03-18 08:37:24

第43 卷第8 期
2009 年8 月
浙江大学学报(工学版)
Journal of Zhejiang University (Engineering Science)
DOI: 10. 3785/j. iss 口. 1008-973X. 2009.08. 016
模拟颗粒布朗运动的格子Boltzmann 模型
聂德明1， 2 林建忠1， 2
(1.浙江大学力学系，浙江杭州310027; 2. 中国汁量学院，汁量测试工程学院，浙江杭州310018)
Vol. 43 No. 8
Aug. 2009
摘要:通过在格子Boltzmann 方法的迭代格式中附加描述分子热运动涨落的分布函数，建立了描述颗粒布朗运
动的涨落格子Boltzmann 模型，给出了分布函数满足的条件以及在D2Q9 格子模型下的具体表达形式.通过Chap
m旦旦-Enskog 展开推导，得到了考虑分子热运动涨落的宏观流体动力学方程.在此基础上，对单颗粒在流场中的布
朗运动进行了数值模拟，得到了颗粒运动的均方速度及速度、角速度时间自相关函数.结果表明，均方速度满足能
量均分定理，说明颗粒最终达到热平衡;颗粒速度、角速度时间自相关函数符合理论预测的t t- 2 衰减规律.数值
结果证明了所建立模型的正确性，为采用格子Boltzmann 方法模拟颗粒的布朗运动提供了有效的方法.
关键词:格子Boltzmann 方法;颗粒;布朗运动
中图分类号: 0359 文献标志码:^ 文章编号: 1008 - 973X(2009)08 -1438 - 05
Lattice Boltzmann model for particle Brownian motion
NIE De-ming1,2 , LIN Jian-zho 且g1， 2
(1. Depα rtment 01 Mechωúcs ， Zhejiαηg University , Hω19zhou 310027 , Chi ηα;
2. College 01 Metrology ωzd Measurement Eη gzneerz ηg ， Chi ηα] iliω19 University , Hω19zhou 310018 , Chi ηα)
Abstract: A fluctuating lattice Boltzmann model for particle Brownian motion was established by incorpora
ting a stochastic term into the lattice Boltzmann equation , which represents the thermally-induced fluctua
tions in the stress tensor. The conditions for the stochastic term were derived and the expressions of the
stochastic term for the D2Q9 lattice model were also presented. The fluctuating hydrodynamic equations
were derived from the lattice Boltzmann equation through Chapman-Enskog expansion. The Brownian mo
tion of a single circular particle was numerically investigated by the newly developed lattice Boltzmann
model. N umerical results including particle mean-square velocity , velocity autocorrelation function and an
gular velocity autocorrelation fu日ction were presented. The energy equipartition theorem was reproduced
by the results of mean-square velocity , which i日dicated that the particle was in thermal equilibrium. The
results showed that the velocity autocorrelation fu日ction and the angular velocity autocorrelation function
decayed as a power law of t- 1 and t- Z respectively , as theoretically stated. Numerical results showed the
accuracy and robustness of the present model , which was proved to be an effective numerical method for
the particle Brownian motion.
Key words: lattice Boltzmann method (LBM); particles; Brownian motion
布朗运动是一种分子热运动涨落现象， 11故布朗
运动的粒子的尺寸一般都在微米量级以下，因为只
有当颗粒很小时，它受到四周流体分子的碰撞才不
平衡，颗粒的运动才显著.布朗运动的研究在某些微
收稿日期: 2008 - 03 - 02. 浙江大学学报(工学版)网址www.journals.zju.edu.cn/eng
基金项目·国家自然科学基金重点资助项目C2005CCA06900)
作者简介.聂德明(1 979-) ，男，福建三明人，博士，从事多相流体力学研究. E-mail: nieinhz@gmail. com
通信联系人·林建忠，男，教授. E-mail: jzlin@sIp. 勾u. edu. cn
第8 期聂德明，等:模拟颗粒布朗运动的格子Boltzman日模型1439
尺度研究领域成为热点，例如生物分子传输、纳米流
体等.纳米流体是指以一定的方式和比例在液体中
添加纳米量级的金属或金属氧化物颗粒，形成一类
新的传热介质.实验结果显示，在液体中添加纳米粒
子会显著提高液体的导热系数和对流换热系数，这
使得纳米流体在强化传热领域具有广阔的应用背景
和潜在的巨大经济价值[lJ 纳米流体在流动与传热
过程中有着与传统两相流不一样的作用机制，这主
要是因为纳米颗粒在流体中受到的作用力与常规颗
粒不同.对于纳米颗粒，重力、相间阻力、浮力等不再
重要，布朗运动和范德华力是主要作用力.当颗粒做
布朗运动时，颗粒的无规则行走强化了热量的传递
过程，同时使得运动更加复杂和难以描述.
传统的求解布朗运动的数值方法有布朗动力学
和Stokes 动力学2 种方法，两者的共同点都是将布
朗颗粒运动的求解分成流体运动和颗粒运动2 部
分，分别利用N-S 方程和郎之万方程求解[2-5J 这2
种方法有几个不足之处:颗粒运动和流体运动之间
是单向捐合，没有考虑颗粒对流体的作用;由于颗粒
运动通过郎之万方程描述，颗粒的形状仅限于圆球
形，没有方向性，对于一些超常规粒子，无法直接描
述;在郎之万方程中，颗粒的受力由平均项和随机项
2 部分力组成，这本身就是一种假设，无法体现分子
热运动的本质.
涨落动力学( fluctuati吨hydrodynamics) 是求
解布朗运动的另一种方法.与传统的数值方法不同，
涨落动力学通过流体力学控制方程中附加的随机应
力项来描述流体中的热涨落[6J 在涨落动力学方法
中，颗粒受到的随机力不再像郎之万方程中一样直
接给出，而是通过求解控制方程得到流场，进而通过
颗粒周围流场来求解其受力，因此随机力实际上直
接由流体分子的热运动产生，从这一点来看，涨落动
力学方法较前面2 种方法更接近物理真实，同时采
确、合理地描述物理系统的微观动力学机制. 1i故布朗
运动的颗粒的尺寸一般小于微米量级，因此以颗粒
作为参考对象，需要重新考虑流体的连续性以及N
S 方程的适用性.从这个角度看，格子Boltzmann 方
法在颗粒的布朗运动数值模拟中具有独特的优势.
Ladd[刊]最早将涨落动力学和格子Boltzmann 方法
结合起来，在格子Boltzmann 演化方程中加上一项
描述分子热运动的分布函数，通过该分布函数来产
生随机应力，从而求解颗粒的布朗运动.但是， Ladd
并没有详细地阐述和推导两者之间的关系.
本文针对目前格子Boltzman日方法中最常用的
D2Q9 格子模型，从演化方程出发，通过多尺度展
开，推导得到了考虑分子热运动涨落的宏观流体动
力学方程，建立了描述颗粒布朗运动的涨落格子
Boltzmann 模型，给出了描述分子热运动涨落的分
布函数应该满足的条件以及表达形式，从而为采用
格子Boltzmann 方法模拟颗粒的布朗运动提供了有
效的方法.
1 数学模型
1. 1 D2Q9 格子模型
格子Boltzmann 方法的演化方程[9J 为
ji(X + c;!.泣， t+ L t) - ji(X ， 。二l [ji(X ， 。
T
j;O) (X ,t) ] ;二0 ，1，…， 8. (1)
式中:X 为空间坐标， Ci 为Z 方向的格子速度， ji 为粒
子分布函数， τ 为松弛时间. D2Q9 格子模型如图1
所示.在本文中，格子间距Lx 和时间步长Lt 均
取1，格子速度c 二Lx/ Lt 二l.式(1)中j;O) 为平衡
分布函数，可以表示为[9J
f;O) 二ωiP [1 + 3CiU + 4. 5 (CiU) 2 - l. 5u2]. (2)
式中:U 为流体速度， ωo 二4/9 ， ω1 二ω2 二ω3 二ω4 二
用这种方法更容易描述颗粒受到的随机力. Hauge 1/9 ， ω5 二ω6 二ω7 二ω8 二1/36.
等[7J 在忽略惯性项的假定下，通过傅里叶变换求解宏观量如密度、速度分别定义为
了圆球的布朗运动，并得出速度的时间白相关函数ρ 二三孔， ρu 二三fι (3)
为一条渐近曲线等重要结果. Sharma 等[8J 利用涨落
动力学方法求解了单颗粒子的布朗运动，并与理论
解进行对比，获得了很好的结果.
格子Boltzmann 方法是20 世纪80 年代中期提
出的流体计算和物理建模的方法.与传统数值方法
不同，该方法并非以连续介质理论为基础，而是基于
微观尺度上的统计力学Boltzmann 方程[9J 因此，与
有限差分、有限元等传统的计算流体力学方法相比，
格子Boltzmann 方法的最大优势在于能够更加准
V/J/
川、/l~
图1 D2Q9 格子模型
Fig.l D2Q9 lattice model
1440 浙江大学学报(工学版) 第43 卷
1. 2 涨落格子BoItzmann 模型
为了描述分子热运动的涨落，可以在演化方程
(1)中添加一项脉动项[叫:
!i (x +CiL泣， t+ L t) - !i (X ，t) 二l [ji(X ， t)
T
!i(O) (X ,t) J + l; i 二0 ，1，…， 8. (的
式中: !f 表示由于分子热运动引起的涨落，满足质
量守恒和动量守恒:
三!f 二0 ，三fici 二o. (5)
下面推导考虑分子热运动涨落的宏观流体动力
学方程.引入2 个宏观时间尺度t1 二εt 和tz 二êZt 以
及一个空间尺度X1 二α ，并对时间导数、空间导数
和分布函数进行多尺度展开:
d t 二ε +ε d α 二ε d1α ，
fi 二fi(O) + ε+ ê
Z !i(Z) + …, !f 二εf~.
将式(的左端泰勒展开，并代入多尺度展开式，
比较各阶系数得到
ε1 :dt1 !;O) + 二l!i(1) + 丘，
T
εZ dtZ!i(O) + (1- } ) (d t1 + d1aCia )!i(1) +
飞4牛乌τ /
f 仙+d川
τ
(6)
(7)
在推导式(7)时利用了式(6) .根据式(6) 求。阶、
1 阶以及2 阶速度矩，得到
d t1 三!i(O) +此三二o ,
d t1 三+d1β 三!i(O)川伊二0 ，
d t1 三!i(O) ωρ +d1Y三二
7 三刀1)川+三f~CiaC伊
从式(的、(9) 可以得到t1 尺度上的宏观方程为
d t1 p+d1α (pu α) 二0 ，
d t1 (严α )+d1卢(rr;g) )二o.
式中:
II~) 二三!i(川， Ciß 二严α川+户。aß
(8)
(9)
(10)
(11)
(1 2)
其中， ρ 为流体密度，户为流体压力，归为Kronecker
符号.式(11)、(1 2) 实际上就是欧拉方程.
由式(1 0) 可得[11J
7 三刀1)川+三f~CiaC伊二
C; [d1α (pu 卢)+d1卢(pu α) J+ο( 扩). (13)
式中:CS 为声速.根据式(1 0) 求。阶、1 阶速度矩并
利用式(1 3) 得到tz 尺度上的宏观方程为
dtZP 二0 ， (1 4)
dtZ (pu α)υd1β [d1α (pu β )+d1β (pua) J 二
内三!~C山(1 5)
式中:υ 二二(τ 1/2 )为流体运动站性系数.由式
(11) 、(1 2) 、(1 4) 、(1日可得宏观流体动力学方程为
dtp+ 孔(pu α) 二0 ， (1 6)
d t (pua)+ 马(puaU β) 二daP +
υJβ [d α(严β)+ 马(严α) J+d 月7 与(17)
式(17)与通常的N-S 方程相比，多了一项，即内二
τ三句，这就是涨落动力学方法中的随机应
力项[汀，是由分子热运动的涨落引起的.
根据文献[6J ，对于不可压流体，随机应力冲满
足以下条件:
〈σ与〉二0 ，
〈σ与(X1 ， t1) σ妇(Xz ， tZ)) 二
A( δJμ+δ川岛γ)δ (X1-XZ) δ (t1-tZ).
(18)
(1 9)
式中:A 二年kBT ， 0 表示取系综平均， μ 为流体动力
勃性系数， ι 为Boltzmann 常数， T 为流体温度.
式(5) 给出了卫应当满足的条件，由此可以假定
!ò 斗， !i 二!~二LL ， f; 二五二1σ
~T ~T
!~二!二+(σ[+σ;，， +σ川，
ζ圭T
!~二!~二+(σ[+σLσ[ )
ζ圭T
(20)
可以证明，式(20) 满足式(5) ，同时满足σ年二
τ三伊在计算过程中， σLJL 、凡通过式
(1 8) 、(1 9) 确定.从式(1 8) 、(1 9) 可以看出， σL 、凡、
σ[ 实际上是均值和方差已知的高斯随机数，在计算
过程中，它们通过一高斯随机数的生成函数未获得.
2 计算方法
计算以圆形单颗粒在无界流场中的布朗运动为
例，忽略颗粒重力、浮力，仅考虑流体作用于颗粒的
力.颗粒直径取d 二10 (格子单位) .
颗粒的运动可以通过牛顿定律求解:
d.Q
M'".~ 二F( 仆， I 丁一二T( t). (2 1)
dt dt
式中:.Q为粒子的转动角速度;M 为颗粒质量; I 为
转动惯量;F( t)， T( t) 分别为颗粒受力和力矩.
Ladd[10 J 以动量守恒为基础，提出了移动边界的
反弹壁面条件，假设壁面在边界格点和流体格点链
接的中点，则
第8 期聂德明，等:模拟颗粒布朗运动的格子Boltzman日模型
j i' (x ,t + 1)二ji(x ， t+)-2Bi(Ci' Ub). (22)
式中:下标u+"表示碰撞后川、J 分别表示流进和流
出粒子边界的方向; Bi 二3ρωi ;Ub 为边界的运动速
度， Ub 二Uo +QX 矶， U。为粒子质心的平动速度，
Xb 二x+cd2 - xo ， x。为粒子质心位置.在Xb 土流
体对粒子的力和力矩分别为
F(x+cd2 ， t) 二2Ci[ji(X ， t+)-Bi(Ci' Ub)J ,
T(x + cd2 ， t) 二Xb X F.
粒子受到的力和力矩分别为
I 二三T(x+ ci/2 ， 仆， F 二三F(x+ cd2 ,0
求和在所有边界点上进行. Ladd[10J 假设粒子内部充
满流体，反弹过程在壁面两边都发生.这个假设的缺
陷是粒子的密度不能比流体密度大. Aidu日等[12J 没
有使用这个假设，令粒子内部不再有流体，当粒子在
格子上移动时，如果某个格点由粒子内部到达粒子外
面，成为流体格子，根据动量定理，相当于流体对粒子
施加了一个力: Fr (x ， t) 二prU(X ,t) ， 相应的对粒子
的力矩为Tr (x ， t) 二(x-Xo)X Fr(x ， t) 反之，如果
某个格子由流体格子成为粒子格子，粒子所受的力和
力矩分别为Fp 仪， t) 二ρrU(X ， 。和Tp (x ， t) 二
(x-xo)XFp(x ,t) ， 总的力和力矩修正为
F 二三F(x+cd2 ,0 +三Fr(x ，t)+ 三Fp(x ， t) ，
I 二三T(x+cd2 ，t) + 三Tr (x ，t)+ 三Tp(x ， t)
然后通过牛顿运动定律(2 1)求解颗粒的运动.
1441
(23)
(24)
一一一(u(O) 叭的)
(v(O)ν(0))
4己t 4
辛辛
o 3
毛2
o
""
o
三O~
o v
;:s
500
当初始速度为零时颗粒均方速度与时间的关系
Particle mean-square velocity vs. time for zero
ini tial veloci ty
300 400
th:u
200 - AVA V
图2
Fig.2
x ， y 方向上均达到了热平衡状态，计算值与理论值
的最大误差不超过10% ，即计算结果与理论预测结
果吻合得较好.
根据布朗运动理论，当颗粒有一定初始速度做
布朗运动时，会在一段时间之后遗忘自己的初始速
度而达到热平衡状态.图3 给出了当u 2 (0) 二hT/
M、2走BT/M、3 走BT/M 时颗粒均方速度与时间的关
系，其中设定hT/M二1X10- 11 • 从图3 可以看到，
无论颗粒的初始速度大小，经过一段时间(t/TB=
200) 之后，颗粒最终会达到热平衡状态.
颗粒速度的时间白相关函数关系是布朗运动研
究的一个重要结果. Hauge 等[7J 在忽略惯性项的假
定下，从理论上求得圆球布朗运动的速度和角速度
的时间白相关函数，证明两者分别为t -3/2 和t -5/2 的
渐进曲线.对于圆形颗粒，两者退化为t- 1 和t- 2 的渐
进曲线.
图4 、5 分别给出了速度、角速度时间白相关函
数的计算结果，其中a ， b 为拟和参数.可以看出，计
算结果有一定的振荡，这和计算时间有关.因为文献
[7J 给出的结果是足够长时间的统计平均，对于有限
~ 10-10
2升
叫「
￥È [、、
s=- - .1 0,,--1"1r I.. 、、飞、ρ"、""、':r 蝇::-:..:._.'Z"古q
~ι丸'、'
。l 丁，..
;:s
初始速度不为零时颗粒均方速度与时间的关系
Particle mea 口square velocity vs. time for 口on
zero ini tial veloci ty
l ( ) 2
一一内的二kJIM
_._.- u'(0)=2kRTIM
u'(O)斗kRTIM
l ()
100
tlru
l () 2 l ()
l AV 2
图3
Fig.3
计算中设定松弛时间τ 二1 ，流体密度ρ 二1 ,
运动站性系数υ 二1/6. 图2 、3 给出了颗粒均方速
度与时间的关系，其中时间通过TB 二M/γ 无量纲
化， γ 为颗粒摩擦系数.布朗运动理论的基本前提是
把布朗颗粒视为流体中的大分子，参与流体分子的
热运动[13J 颗粒在流体中经过足够长的时间，会达
到同周围流体介质处于热平衡的状态.因此，根据能
量均分定理，有以下关系:
(u(O)u(O)) 二(v(O)v(O)) 二kBT/M. (25)
式中:U"V 分别为做布朗运动的颗粒沿x ， y 方向的
速度.
图2 给出了在颗粒初始速度为零条件下的均方
速度随时间变化的曲线.计算中分别设定A 二1X
10- 9 、2 X 10- 9 、3X10- 9 ( 格子单位) ，颗粒质量M 二
300(格子单位) .根据式(17)可得， kBT/M 二1X
10- 11 、2 X 10- 11 、3X10- 11 • 从图2 可以看出，做布朗
运动的颗粒经过一段时间。IrB=100) 振荡之后，在
3 结果及讨论
第43 卷
(u(O)u(t))
i>. (v(O)v(l))
-ar'
报一(工学1反)
宏观之间的介观层次上的数值计算方法，将该方法
与涨落动力学方法相结合，无疑发挥了格子Boltz
man日方法的最大优势，即编程容易、计算效率高和
边界条件实施简单.本文为采用格子Boltzmann 方
法模拟颗粒的布朗运动提供了有效的方法.
学
学
p工大
浙
1442
[lJ 宣益民，李强，姚正平.纳米流体的格子Boltzmann 模拟[J].
中国科学(E 辑:技术科学)， 2004 , 34(3): 280 - 287.
XU^N Yi-mi 口， LI Qiang , Y ^O Zheng-ping. Numerical
simulatio 丑。f 口旦旦o-fluids using lattice Boltzm旦旦旦method
[JJ. Science in China Seri创(E: Technological Sciences) ,
2004 , 34(3): 280 - 287.
[2J BR^DY J F , BOSSIS G. Stokesia 口dynamics [J]. Annual
Review of Fluid Mechanics , 1988 , 20: 111 - 157.
[3J BR^DY J F , BOSSIS G. Self-diffusion of Brownia 口par
ticles i 口concentrated suspensions u口der shear [J]. J ournal
of Chemical Physics , 1987 , 87(9): 5437 - 5448.
[4J FOSS D R , BR^DY J F. Structure, diffusion and rheol
ogy of Brownia 口suspensions by Stokesian dynamics
simulatio 口[J]. Journal of Fluid Mechanics , 2000 , 407:
167-200.
[5J ERM^K D L , MCC^MMON J A Brow口1a口dynamics
with hydrodynamic interactions [J]. Journal of Chemical
Physics , 1978 , 69 (4): 1352 - 1360.
[6J L^ND^U L D, LIFSHITZ E M. Fluid mechanics [MJ.
Lo 口do 口Pergamo 口Press ， 1959.
[7J H^UGE E H , M^RTIN-LOF A Fluctuating hydrody
namics and Brownia口mot lO口[J J. J ournal of Statistical
Physics , 1973 , 7(3): 259 - 28 1.
[8 J SH^RM^ N , P ^ T ^NK^R N A Direct numerical
simulatio 丑。f the Brownian motio 丑。f particles by using
fluctuating hydrodynamic equatio 口s [J J. J ournal of
Computational Physics , 2004 , 201(2): 466 - 486.
[9J CHEN S , DOOLEN G D. Lattice Boltzmann method for
fluid flows [J]. Annual Review of Fluid Mechanics ,
1998 , 30: 329 - 364.
[10J L^DD ^ J C. Numerical simulations of particulate sus
pe口S10 且s via a discretized Boltzmann equatio 口. Part 1:
Theoretical foundatio 口[JJ. Journal of Fluid Mechanics
, 1994 , 271: 285 - 309.
[11J 郭照立，郑楚光，李青，等.流体动力学的格子Boltz
mann 方法[MJ.武汉:湖北科学技术出版社， 2002:
49 - 50.
[12J ^IDUN C K , LU y , DING E J. Direct analysis of par
ticulate suspe口S10 且s with inertia using the discrete
Boltzmann equatio 口[JJ. Journal of Fluid Mechanics ,
1998 , 373: 287 - 31 1.
[13J 熊吟涛.统i十物理学[M J.北京:人民教育出版社，
1981.
参考文献(References) :
。υ
l
图4 速度时间自相关函数
Time-dependent particle velocity autocorrelation
。
但
主叶1.0
尘0.8
ζ乙0.6
;:::
号〉
~ 0.4
1己O.úll.'
;:::
号〉
;:::
l AV
图5 角速度时间自相关函数
Time-dependent particle rotational velocity auto
correlatio 口fu 口ctio 口
R
~
8
8
A
口(.0 (0) .0(1))
-br'
4 6
tlrR
国占口口
E
4
2
fu 口ctio 口
。υ
白V
。
0.8
4己
主叶0.6
冉
号。4
c: 0.2
o
q
、/
Fig.4
Fig.5
时间的计算，小范围的数值振荡是正常的;而且计算
结果的趋势非常接近2 条渐近曲线，因此可以推断，
随着计算时间的增加，速度、角速度时间白相关函数
必定呈t- 1 和t- Z 的渐进曲线.
语
本文通过在格子Boltzmann 方法迭代格式中附
加一项描述分子热运动涨落的分布函数，建立了颗
粒布朗运动的涨落格子Boltzmann 模型;通过多尺
度展开，推导得到了相应的宏观流体动力学方程，证
明了附加的分布函数实际上就是涨落动力学方法中
在宏观流体动力学N-S 方程里附加的随机应力张
量.本文对颗粒的布朗运动进行了数值模拟，给出了
颗粒运动的均方速度及速度相关函数，所得计算结
果与理论预测结果吻合较好.
涨落动力学方法是求解颗粒在流体中做布朗运
动的直接数值模拟方法之一，该方法不仅在物理上
更加接近颗粒运动本质，而且对于不同尺寸、形状的
超常多颗粒运动具有实施简单、直接的优势.然而应
用传统数值方法求解N-S 方程存在一定难度，编程
方面不易实施.格子Boltzmann 方法是介于微观和
4 生.:op士

请您先登陆，再发跟帖！