清华大学学报(自然科学版)2000年第40卷第4期
CN11-2223/NJTsinghuaUniv(Sci&Tech),2000,Vol.40,No.415/34
5154
完全气体格子Boltzmann热模型*
孙成海, 王保国, 沈孟育
(清华大学工程力学系,北京100084)
文 摘:为建立一种具有任意比热比的完全气体多速度格子Boltzmann热模型,引入粒子的势能来调整压能与热力学能的关系;利用Chapman-Enskog方法从BGK型的格子Boltzmann方程推导出了Navier-Stokes方程和能量方程。对一维正弦波形式的能量衰减过程进行了模拟,测得的热扩散率与理论预测值相吻合。还模拟了绕加热平板的二维强制热对流问题,结果合理。
关键词:格子Boltzmann;完全气体;热模型中图分类号:O354
文献标识码:
A
度方向的个数;Njκ表示具有速度cjκ的粒子数密度。Fermi-Dirac型的平衡分布使格子气失去了Galilean不变性。这是格子气的缺陷。对于LB方法,平衡分布的选择有一定的自由度。如果选择适当可以避免这个问题[4]。平衡分布由下式给出N=deqjκ
κ221+cjκ v+(cjκ v)-v
cκ2cκD
2
,(1)
其中:D是空间维数,v是流体速度,dκ是由流体密度和能量决定的量。
BGK型格子Boltzmann方程为
Njκ(x+cjκΔt,t+Δt)-Njκ(x,t)=Kjκ,(2)其中:
Kjκ=-eq
(Njκ(x,t)-Njκ)f
(3)
文章编号:1000-0054(2000)04-0051-04
自从1986年出现了满足Navier-Stokes(N-S)方程的二维FHP格子气模型以来,已经对许多物理问题建立了相应的模型,其中包括质量扩散模型、热模型等[2,3]。格子气具有Fermi-Dirac型的平衡分布,从而导致了非Galilean不变性。为了克服这一缺陷,人们提出了格子Boltzmann(LB)方法[4,5]。文[6]提出了LB方法的质量扩散模型。由于能量方程是非独立的,所以不能模拟热输运问题。一些作者通过引入多重粒子速度使能量方程成为独立方程,从而建立了LB热模型[7,8]。这些模型一般只适用于比热比V=2的完全气体。作者通过引入粒子的势能来调整压能与内能的关系,从而建立了具有任意比热比的LB模型。
[1]
为碰撞项,x为网格节点,t为时间,f是松弛因子,稳定性要求f>2。
流体的密度、动量和能量定义为:
d=dv=de=
∑mN
κ,j
jκ
,(4)(5)(6)
∑mN
κ,j
jjκjκ
c,
∑mN
κ,j
2
cκ+ep,2
1 LB热模型
在规则的网格上,设cjκ是第j个网格方向的模为cκ的粒子速度;κ=1,2;j=1,…,b;b是粒子速
收稿日期:1998-06-27
作者简介:孙成海(1960-),男(汉),山东,副教授 *基金项目:国家自然科学基金项目(19672030,
19972037)和教育部留学回国人员经费项目
式中:m是粒子的质量,ep是粒子单位质量的势能。这里引入了粒子势能是为了能够增加一个自由
2
度,从而获得任意比热比。粒子的能量由动能mcκ
2
和势能mep两部分组成。下面确定ep和式(1)中的
dκ。平衡分布Njκ应该满足式(4～6)。其中式(5)自然满足,式(4)和(6)变为
d=de=
eq
∑mbd
κκ
κ
,(7)(8)
∑mbd
κ
c2+dep,
κ
52
清华大学学报(自然科学版)2000,40(4)
0jκ
另外,比热比为V的完全气体满足状态方程p=(V
-1)du,这里u=e-v2是比热力学能。而压力表
2
示如下
22
p=mbdκc-dv.(9)κ
D∑Dκ 在低Mach数条件下可忽略v,得到
2
(V-1)de=∑mbdκcκ.(10)
Dκ
当粒子有二级速度(κ=1,2)时,由式(7,8,10)解得:
d1=
dc,bm(c2-c1)
2
22
2
jκ0
,N cjκ+L FYΔ
=-c2+epcjκ
∑Lt2κκ,j
00jκ
Njκ+XT-Njκ cjκ+ F0
2LYΔ 0mepcjNjκ+XT-∑κ,jΔ
jκ
Njκ cjκ+ F0.
LYΔ
(20)
+
(21)
为了获得N-S方程,要求粒子速度4阶张量之
和是各向同性的,即
bcκ
cjκcjκcjκcjκ= ∑D(D+2)j
(WTVWUW+WTWWUV+WTUWVW)eTeUeVeW
其中:eT,T=1,…,D是空间的一组正交单位基向
量。在二维空间中,正六边形网格可以满足上式。方程(19～21)经过化简可得到如下的连续方程、N-S方程和能量方程:
v)=0,(22)Lt+div(d
+div(dvv)=-p+
LtΔ (_v)+[(_v)]T-div(_v),(23)ΔΔ
+dv e=XTf-(V-1) LtΔ
2222c1c2c1+c2
-d+de+
(V-1)2 2
1-(V-1)de-e(de)+dv ep,2ΔΔ
(24)其中:I是二阶单位张量,_是粘度,且
2
_=f-mbdκcκ.(25)∑D+κ
在式(24)中,右端最后一项d vep是由于粒子
Δ
势能的变化引起的。它破坏了能量方程的Galilean不变性。
4
D(V-1)de-dc1
d2=,
bm(c2-c1)
(V-
1)e.2
为了方便将宏观变量写成向量形式
ep=
1-Y=(d,dv,de),
T
(11)(12)(13)
采用Chapman-Enskog方法由方程(2)推导宏观守
eq
恒方程。将Njκ在Njκ附近进行渐近展开得
Njκ(x,t)=
eq
∑N
n=0∞
∞
njκ
(x,t)X,
n
(14)
其中:Njκ=Njκ,X为一小量。再将宏观变量对时间的导数进行渐近展开
Lt=
∑FX.
nn
n=0
(15)
T
2
令:Δt=XT,ηjκ=m,mcjκ,mc+ep。将
2κ
方程(2)进行Taylor展开,再将式(3,14,15)代入展开后的方程,考虑式(4～6),比较X的同次项可确定
出:
0
F=-∑ηjκcjκ Njκ,
κ,jΔ
1
F=-∑ηjκcjκ
κ,j
1jκ
(16)
0LNjκ0
N+2,(17)Njκ cjκ+ F
LYΔ
1LNjκNjκ=-Njκ cjκ+ F0.(18)
LYΔ
如果只取前2阶,将式(16～18)代入式(15)可得守恒方程
=-
Lt=-t2 算 例
2.1 用正弦波形测量热扩散率
这个算例是在流动速度很小的条件下进行的。
这要求压力梯度很小,即de很小。如果将其忽略,
Δ
则能量方程(24)可化简成:
+(V-1)v e=div(ae),(26)Lt2a=
孙成海,等: 完全气体格子Boltzmann热模型
53
12.(27)+(V
-1)1-2De2(V-1) 采用六边形网格,粒子的速度分别为1和
3
(见图1)。如果速度很小可以忽略对流项,且忽略热
扩散率a的变化,则式(
26)有解析解
e(x,t)=a′+b′exp(-at)sin(2πx/L).
(28)
22
2.2 对流扩散问题
在相对值u=0.5,d=2.5,e=2.0的均匀流场
中放入一加热板,观察其能量扩散情况。采用六边形网格,100×80个结点;上下边界采用周期性边界条件;在x=25,y=40为中心处,放一尺寸为30×4个节点、e=2.5的加热板;取V=1.4,f=1.0.在图4a和4b中分别给出了在相对值t=60和t=100时e的等值线。从中可以看出能量向下游扩散的情况。与u=0.5的流动速度相比,能量向下游扩散的速度相对实际情况要慢一些。这是由于能量方程的非Galilean不变性造成的。在方程(26)中对流项v
e的因数本应该为1,这里却是0.4。Δ
图1 粒子速度
如果在某一时刻t=t1(t1可取为振幅衰减一半所用的时间)测得e(x,t1),则可根据式(28)算出热扩散率a。这样就可以通过数值模拟测出模型所对应的扩散系数。取a′=2,b′=0.05,V=1.4;100×4个网格结点。图2给出了热扩散率随f的变化。实线是由式(27)计算出的理论值
,圆点是利用上述方法得到的测量值,两者吻合很好。图
3给出了t=0和t=400时相对值e的分布。与式(28)给出的解析解相一致。
(a)t=60
(b)t=100
图2 热扩散率随τ的变化
图4 e的等值线
3 总 结
建立了具有任意比热比的LB完全气体热扩散
模型,利用Chapman-Enskog渐进展开法得到了连续性方程、Navier-Stokes方程和能量方程。能量方程的形式比较复杂,经过简化得到了能量的对流扩散方程。通过数值模拟对热扩散率进行了检验,得到了与理论预测相一致的结果。热扩散率和粘度都是
3t0t=时e的分松弛因子的线性函数,可以灵活调整。对于平板加热
54
清华大学学报(自然科学版)
[8]
2000,40(4)
果。由于引进了非常数的粒子势能,能量方程的Galilean不变性受到了破坏。能量向下游扩散的速度比实际情况要慢一些。此外,本模型的压力还依赖于速度[见式(9)]。这两点不足之处作者在最近提出的局部自适应LB模型中得到了很好解决[9,10]。
ChenY,OhashiH,AkiyamaM.HeattransferinlatticeBGKmodeledfluid[J].81(1/2):7185.
JStatPhys,1995,
[9]SunChenghai.LatticeBoltzmannmodelsforhighspeedflows[J].PhysRevE,1998,58(6):7287.
7283
[10]SunChenghai,WangBaoguo,ShenMengyu.
[参考文献] References
[1]
FrischU,HasslacherB,PomeauY.Lett,1986,56:15051508.[2]
BernardinD,
Sero-GuillaumeO,
SunChenghai.
Multispecies2Dlatticegaswithenergylevels:diffusiveproperties[J].PhysicaD,1991,47:169188.[3]
SunChenghai.
Contributionà
L 'Etudedela
ThermodynamiquedesGazSurRéseaux[D].Nancy:L'InstitutNationalPolytechniquedeLorrain,1993.[4]
ChenH,ChenS,MatthaeusW.RecoveryoftheNavier-Stokesequationusingalattice-gasBoltzmannmethod[J].PhysRevA,1992,45(8):53395342.[5]
QianY,d'HumièresD,LallemandP.LatticeBGKmodelsforNavier-Stokesequation[J].Lett,1992,17:479484.
[6]孙成海.多组份流体质量扩散的格子Boltzmann方法
[J].力学学报,1998,30(1):2026.Sun
Chenghai.
Multispecies
lattice-boltzmannActaMechanica
modelsformassdiffusion[J].[7]
McNamaraG,
AlderB.
EurophysLatticegas
automatafortheNavier-Stokesequation[J].Phys
AdaptivelatticeBoltzmannmodelforcompressibleflows[J].TsinghuaScienceandTechnology,2000,5(1):4346.
Thermallattice-Boltzmannmodel
forperfectgas
SUNChenghai,WANGBaoguo,SHENMengyu
(DepartmentofEngineeringMechanics,TsinghuaUniversity,Beijing100084,China)
Abstract:Multi-speedthermallatticeBoltzmannmodelwasformulatedforaperfectgaswitharbitraryspecificheatratiothroughtheintroductionofaparticlepotentialenergywhichadjuststherelationbetweenthepressureenergyandthermodynamicenergy.TheNavier-StokesequationsandtheenergyequationwerederivedfromtheBGKtypelatticeBoltzmannequationbytheChapman-Enskogmethod.Thethermaldiffusivitypredictedusingaone-dimensional
simulationwithasinusoidalenergydistributioncomparedwellwiththetheoreticallypredictedvalue.Reasonableresultswerealsoobtainedfortwo-dimensionalforcedconvectionforuniformflowpassingoveraheatedplateinthecenterofachannel.
Keywords:latticeBoltzmann;perfectgas;thermalmodel
Sinica,1998,30(1):2026.(inChinese)
AnalysisoftheLattice
Boltzmanntreatmentofhydrodynamics[J].PhysicaA,1993,194:218228.
清华大学学报(自然科学版)2000年第40卷第4期
CN11-2223/NJTsinghuaUniv(Sci&Tech),2000,Vol.40,No.415/34
5154
完全气体格子Boltzmann热模型*
孙成海, 王保国, 沈孟育
(清华大学工程力学系,北京100084)
文 摘:为建立一种具有任意比热比的完全气体多速度格子Boltzmann热模型,引入粒子的势能来调整压能与热力学能的关系;利用Chapman-Enskog方法从BGK型的格子Boltzmann方程推导出了Navier-Stokes方程和能量方程。对一维正弦波形式的能量衰减过程进行了模拟,测得的热扩散率与理论预测值相吻合。还模拟了绕加热平板的二维强制热对流问题,结果合理。
关键词:格子Boltzmann;完全气体;热模型中图分类号:O354
文献标识码:
A
度方向的个数;Njκ表示具有速度cjκ的粒子数密度。Fermi-Dirac型的平衡分布使格子气失去了Galilean不变性。这是格子气的缺陷。对于LB方法,平衡分布的选择有一定的自由度。如果选择适当可以避免这个问题[4]。平衡分布由下式给出N=deqjκ
κ221+cjκ v+(cjκ v)-v
cκ2cκD
2
,(1)
其中:D是空间维数,v是流体速度,dκ是由流体密度和能量决定的量。
BGK型格子Boltzmann方程为
Njκ(x+cjκΔt,t+Δt)-Njκ(x,t)=Kjκ,(2)其中:
Kjκ=-eq
(Njκ(x,t)-Njκ)f
(3)
文章编号:1000-0054(2000)04-0051-04
自从1986年出现了满足Navier-Stokes(N-S)方程的二维FHP格子气模型以来,已经对许多物理问题建立了相应的模型,其中包括质量扩散模型、热模型等[2,3]。格子气具有Fermi-Dirac型的平衡分布,从而导致了非Galilean不变性。为了克服这一缺陷,人们提出了格子Boltzmann(LB)方法[4,5]。文[6]提出了LB方法的质量扩散模型。由于能量方程是非独立的,所以不能模拟热输运问题。一些作者通过引入多重粒子速度使能量方程成为独立方程,从而建立了LB热模型[7,8]。这些模型一般只适用于比热比V=2的完全气体。作者通过引入粒子的势能来调整压能与内能的关系,从而建立了具有任意比热比的LB模型。
[1]
为碰撞项,x为网格节点,t为时间,f是松弛因子,稳定性要求f>2。
流体的密度、动量和能量定义为:
d=dv=de=
∑mN
κ,j
jκ
,(4)(5)(6)
∑mN
κ,j
jjκjκ
c,
∑mN
κ,j
2
cκ+ep,2
1 LB热模型
在规则的网格上,设cjκ是第j个网格方向的模为cκ的粒子速度;κ=1,2;j=1,…,b;b是粒子速
收稿日期:1998-06-27
作者简介:孙成海(1960-),男(汉),山东,副教授 *基金项目:国家自然科学基金项目(19672030,
19972037)和教育部留学回国人员经费项目
式中:m是粒子的质量,ep是粒子单位质量的势能。这里引入了粒子势能是为了能够增加一个自由
2
度,从而获得任意比热比。粒子的能量由动能mcκ
2
和势能mep两部分组成。下面确定ep和式(1)中的
dκ。平衡分布Njκ应该满足式(4～6)。其中式(5)自然满足,式(4)和(6)变为
d=de=
eq
∑mbd
κκ
κ
,(7)(8)
∑mbd
κ
c2+dep,
κ
52
清华大学学报(自然科学版)2000,40(4)
0jκ
另外,比热比为V的完全气体满足状态方程p=(V
-1)du,这里u=e-v2是比热力学能。而压力表
2
示如下
22
p=mbdκc-dv.(9)κ
D∑Dκ 在低Mach数条件下可忽略v,得到
2
(V-1)de=∑mbdκcκ.(10)
Dκ
当粒子有二级速度(κ=1,2)时,由式(7,8,10)解得:
d1=
dc,bm(c2-c1)
2
22
2
jκ0
,N cjκ+L FYΔ
=-c2+epcjκ
∑Lt2κκ,j
00jκ
Njκ+XT-Njκ cjκ+ F0
2LYΔ 0mepcjNjκ+XT-∑κ,jΔ
jκ
Njκ cjκ+ F0.
LYΔ
(20)
+
(21)
为了获得N-S方程,要求粒子速度4阶张量之
和是各向同性的,即
bcκ
cjκcjκcjκcjκ= ∑D(D+2)j
(WTVWUW+WTWWUV+WTUWVW)eTeUeVeW
其中:eT,T=1,…,D是空间的一组正交单位基向
量。在二维空间中,正六边形网格可以满足上式。方程(19～21)经过化简可得到如下的连续方程、N-S方程和能量方程:
v)=0,(22)Lt+div(d
+div(dvv)=-p+
LtΔ (_v)+[(_v)]T-div(_v),(23)ΔΔ
+dv e=XTf-(V-1) LtΔ
2222c1c2c1+c2
-d+de+
(V-1)2 2
1-(V-1)de-e(de)+dv ep,2ΔΔ
(24)其中:I是二阶单位张量,_是粘度,且
2
_=f-mbdκcκ.(25)∑D+κ
在式(24)中,右端最后一项d vep是由于粒子
Δ
势能的变化引起的。它破坏了能量方程的Galilean不变性。
4
D(V-1)de-dc1
d2=,
bm(c2-c1)
(V-
1)e.2
为了方便将宏观变量写成向量形式
ep=
1-Y=(d,dv,de),
T
(11)(12)(13)
采用Chapman-Enskog方法由方程(2)推导宏观守
eq
恒方程。将Njκ在Njκ附近进行渐近展开得
Njκ(x,t)=
eq
∑N
n=0∞
∞
njκ
(x,t)X,
n
(14)
其中:Njκ=Njκ,X为一小量。再将宏观变量对时间的导数进行渐近展开
Lt=
∑FX.
nn
n=0
(15)
T
2
令:Δt=XT,ηjκ=m,mcjκ,mc+ep。将
2κ
方程(2)进行Taylor展开,再将式(3,14,15)代入展开后的方程,考虑式(4～6),比较X的同次项可确定
出:
0
F=-∑ηjκcjκ Njκ,
κ,jΔ
1
F=-∑ηjκcjκ
κ,j
1jκ
(16)
0LNjκ0
N+2,(17)Njκ cjκ+ F
LYΔ
1LNjκNjκ=-Njκ cjκ+ F0.(18)
LYΔ
如果只取前2阶,将式(16～18)代入式(15)可得守恒方程
=-
Lt=-t2 算 例
2.1 用正弦波形测量热扩散率
这个算例是在流动速度很小的条件下进行的。
这要求压力梯度很小,即de很小。如果将其忽略,
Δ
则能量方程(24)可化简成:
+(V-1)v e=div(ae),(26)Lt2a=
孙成海,等: 完全气体格子Boltzmann热模型
53
12.(27)+(V
-1)1-2De2(V-1) 采用六边形网格,粒子的速度分别为1和
3
(见图1)。如果速度很小可以忽略对流项,且忽略热
扩散率a的变化,则式(
26)有解析解
e(x,t)=a′+b′exp(-at)sin(2πx/L).
(28)
22
2.2 对流扩散问题
在相对值u=0.5,d=2.5,e=2.0的均匀流场
中放入一加热板,观察其能量扩散情况。采用六边形网格,100×80个结点;上下边界采用周期性边界条件;在x=25,y=40为中心处,放一尺寸为30×4个节点、e=2.5的加热板;取V=1.4,f=1.0.在图4a和4b中分别给出了在相对值t=60和t=100时e的等值线。从中可以看出能量向下游扩散的情况。与u=0.5的流动速度相比,能量向下游扩散的速度相对实际情况要慢一些。这是由于能量方程的非Galilean不变性造成的。在方程(26)中对流项v
e的因数本应该为1,这里却是0.4。Δ
图1 粒子速度
如果在某一时刻t=t1(t1可取为振幅衰减一半所用的时间)测得e(x,t1),则可根据式(28)算出热扩散率a。这样就可以通过数值模拟测出模型所对应的扩散系数。取a′=2,b′=0.05,V=1.4;100×4个网格结点。图2给出了热扩散率随f的变化。实线是由式(27)计算出的理论值
,圆点是利用上述方法得到的测量值,两者吻合很好。图
3给出了t=0和t=400时相对值e的分布。与式(28)给出的解析解相一致。
(a)t=60
(b)t=100
图2 热扩散率随τ的变化
图4 e的等值线
3 总 结
建立了具有任意比热比的LB完全气体热扩散
模型,利用Chapman-Enskog渐进展开法得到了连续性方程、Navier-Stokes方程和能量方程。能量方程的形式比较复杂,经过简化得到了能量的对流扩散方程。通过数值模拟对热扩散率进行了检验,得到了与理论预测相一致的结果。热扩散率和粘度都是
3t0t=时e的分松弛因子的线性函数,可以灵活调整。对于平板加热
54
清华大学学报(自然科学版)
[8]
2000,40(4)
果。由于引进了非常数的粒子势能,能量方程的Galilean不变性受到了破坏。能量向下游扩散的速度比实际情况要慢一些。此外,本模型的压力还依赖于速度[见式(9)]。这两点不足之处作者在最近提出的局部自适应LB模型中得到了很好解决[9,10]。
ChenY,OhashiH,AkiyamaM.HeattransferinlatticeBGKmodeledfluid[J].81(1/2):7185.
JStatPhys,1995,
[9]SunChenghai.LatticeBoltzmannmodelsforhighspeedflows[J].PhysRevE,1998,58(6):7287.
7283
[10]SunChenghai,WangBaoguo,ShenMengyu.
[参考文献] References
[1]
FrischU,HasslacherB,PomeauY.Lett,1986,56:15051508.[2]
BernardinD,
Sero-GuillaumeO,
SunChenghai.
Multispecies2Dlatticegaswithenergylevels:diffusiveproperties[J].PhysicaD,1991,47:169188.[3]
SunChenghai.
Contributionà
L 'Etudedela
ThermodynamiquedesGazSurRéseaux[D].Nancy:L'InstitutNationalPolytechniquedeLorrain,1993.[4]
ChenH,ChenS,MatthaeusW.RecoveryoftheNavier-Stokesequationusingalattice-gasBoltzmannmethod[J].PhysRevA,1992,45(8):53395342.[5]
QianY,d'HumièresD,LallemandP.LatticeBGKmodelsforNavier-Stokesequation[J].Lett,1992,17:479484.
[6]孙成海.多组份流体质量扩散的格子Boltzmann方法
[J].力学学报,1998,30(1):2026.Sun
Chenghai.
Multispecies
lattice-boltzmannActaMechanica
modelsformassdiffusion[J].[7]
McNamaraG,
AlderB.
EurophysLatticegas
automatafortheNavier-Stokesequation[J].Phys
AdaptivelatticeBoltzmannmodelforcompressibleflows[J].TsinghuaScienceandTechnology,2000,5(1):4346.
Thermallattice-Boltzmannmodel
forperfectgas
SUNChenghai,WANGBaoguo,SHENMengyu
(DepartmentofEngineeringMechanics,TsinghuaUniversity,Beijing100084,China)
Abstract:Multi-speedthermallatticeBoltzmannmodelwasformulatedforaperfectgaswitharbitraryspecificheatratiothroughtheintroductionofaparticlepotentialenergywhichadjuststherelationbetweenthepressureenergyandthermodynamicenergy.TheNavier-StokesequationsandtheenergyequationwerederivedfromtheBGKtypelatticeBoltzmannequationbytheChapman-Enskogmethod.Thethermaldiffusivitypredictedusingaone-dimensional
simulationwithasinusoidalenergydistributioncomparedwellwiththetheoreticallypredictedvalue.Reasonableresultswerealsoobtainedfortwo-dimensionalforcedconvectionforuniformflowpassingoveraheatedplateinthecenterofachannel.
Keywords:latticeBoltzmann;perfectgas;thermalmodel
Sinica,1998,30(1):2026.(inChinese)
AnalysisoftheLattice
Boltzmanntreatmentofhydrodynamics[J].PhysicaA,1993,194:218228.