非正交多松弛系数轴对称热格子Boltzmann方法∗

王佐 张家忠 王恒

(西安交通大学能源与动力工程学院,西安 710049)

非正交多松弛系数轴对称热格子Boltzmann方法∗

王佐 张家忠†王恒

(西安交通大学能源与动力工程学院,西安 710049)

(2016年9月3日收到;2016年11月22日收到修改稿)

提出了一种模拟轴对称热流动的非正交多松弛系数格子Boltzmann(MRT-LB)模型.通过采用非正交转换矩阵,建立了基于D2Q9离散速度的求解流场的MRT-LB模型和基于D2Q5离散速度的求解温度场的MRT-LB模型.Chapman-Enskog分析表明,该模型可以恢复对应的柱坐标下的宏观连续方程、动量方程和能量方程.与现有的轴对称MRT-LB模型相比,本文采用的非正交转换矩阵中含有更多的零元素,因而具有更高的计算效率.采用本文模型对几种典型的轴对称热流动问题,包括热Womersley流动、竖直圆柱体内的Rayleigh-Bénard对流和环形柱体内的自然对流进行了数值模拟,通过等温线图和流线图以及定量数据的对比,验证了本文模型的可行性和可靠性.并且数值模拟结果表明,相对现有模型,本文模型具有更好的数值稳定性和计算效率.

轴对称热流动,格子Boltzmann方法

1 引 言

轴对称流动在换热器、太阳能设备、电子器件冷却等工业领域有着重要应用.近年来有许多文献对轴对称流动换热现象进行了数值研究[1,2].格子Boltzmann(LB)方法作为一种新兴的数值方法,由于其分子动力学背景,具有简单高效、边界处理简单、并行性强等特点,在许多领域获得了重要应用[3-5].其中单松弛系数模型格子Boltzmann(SRT-LB,也称作LBGK)模型由于其实现简单,具有最为广泛的应用.近年来,具有更加合理的理论基础的多松弛系数格子Boltzmann(MRT-LB)模型也取得了飞速发展[6,7].

截止目前已有大量研究致力于将LB推广到轴对称流动问题的研究.Halliday等[8]首先基于宏观柱坐标下的Navier-Stokes(NS)方程提出了一种适用于轴对称流动的LB模型.随后,Lee等[9]发现Halliday等[8]的模型缺失了一些径向速度相关项,不能正确恢复宏观的柱坐标下的NS方程.为克服此缺点,Lee等[9]提出了一种新的轴对称LB模型.此后,Reis和Phillips[10,11]对Lee等[9]的模型进行了部分改进.但这些模型中均存在梯度项,对梯度项采用差分方法处理,降低了模型的数值准确性和稳定性,破坏了LB固有的计算局部化特性.与前面几个基于宏观柱坐标下的NS方程推导的LB模型不同,Guo等[12]直接从描述柱坐标流动的Boltzmann方程出发,推导出了一种柱坐标下的LB模型.此模型不包含复杂的梯度项,且其实施流程几乎与标准LB完全相同.Li等[13]从二维笛卡尔坐标下的标准LB模型出发,将宏观柱坐标下的NS方程转化为类似标准的笛卡尔坐标下的NS方程加入与径向坐标r有关的多个源项的形式,并提出了针对此种形式宏观方程的LB模型,详细阐述了从标准LB模型出发构造轴对称流动LB模型的过程.Zhou[14]的模型对Li等[13]的模型进行了改进.

以上模型均针对等温流场,为发展非等温的轴对称LB模型,Peng等[15]提出了一种采用LB模型求解速度场,而采用二阶差分方法求解温度场的混合模型.然而,Huang等[16]指出,随着对流作用的增大,对流项采用二阶差分方法容易导致发散,并对此混合模型进行了改进.然而,采用混合模型在很大程度上丧失了LB简单易行的优势.Li等[17]提出了一种求解温度场的双分布LB模型,其基本思想与文献[13]类似.Zheng等[18]提出了另一种与Guo等[12]的方法类似的更加简单的LB模型.

上述轴对称模型均采用LBGK,即单松弛系数(SRT)碰撞模型.在该模型中,所有的物理量的弛豫时间均假设相同,这显然与实际的物理过程不符.2000年,Lallemand和Luo[6]提出了一种多松弛系数LB模型(MRT-LB),在该模型中,通过引入转换矩阵和碰撞矩阵,用不同的松弛系数来表示流体系统中不同物理量的弛豫.与SRT-LB相比,MRT-LB显然具有更加合理的理论基础.事实上,大量的数值模拟结果表明,MRT-LB模型相对于SRT-LB模型具有更好的数值稳定性和边界处理准确性[5].由于这些优点,MRT-LB模型先后被应用到流动传热[19]、多孔介质流动[20]、多相流[5]和微尺度流动等[21]领域.这些MRT-LB模型均求解二维笛卡尔坐标下流体相关问题.2010年,Li等[13]和Wang等[22]将MRT-LB模型推广到轴对称流动领域,几乎同时提出了针对柱坐标下轴对称流动的不同的MRT-LB模型.

现有的所有MRT-LB模型均建立在正交化转换矩阵的基础上,该矩阵通过Gram-Schmidt正交化方法获得,相应的这些MRT-LB模型可称作正交MRT-LB模型.然而,最近Premnath和Banerjee[23]以及Geier等[24]指出,Gram-Schmidt正交化方法并不是获得转换矩阵的惟一途径,并且提出了一种适用于二维笛卡尔坐标流场的非正交MRTLB模型.随后,Liu等[25]根据这一思想,进一步提出了一种适用于二维笛卡尔坐标下的温度场的非正交MRT-LB模型.与现有的正交MRT-LB模型相比,非正交的MRT-LB模型中转换矩阵具有更多的零元素,因而计算效率获得一定程度上的提高[25].

基于前人的工作,本文提出了一种针对轴对称热流动的MRT-LB模型,该模型中流场和温度场均采用非正交转换矩阵,相较于现有的模型,本文模型中的转换矩阵含有更多的零元素,这使得模型在保持数值稳定性良好的前提下其演化过程更加高效.

2 双分布轴对称MRT-LB模型

2.1 控制方程

柱坐标下的流动传热问题是一个三维问题,然而,当周向速度为零时,流动可简化为二维问题,其控制方程为[12,17]

其中,x和r为分别为轴向和径向坐标;ux,ur分别为轴向和径向速度;ax,ar分别为轴向和径向方向的外力;ρ为密度,T为温度,ν为黏性系数,k为热扩散系数.对于任意变量φ,

2.2 求解速度场的MRT-LB模型

格子Boltzmann方法通过描述具有离散速度的流体粒子分布函数的变化过程来反映流体粒子的运动过程.基于非正交转换矩阵,本文提出一种适用于轴对称流场的多松弛系数格子Boltzmann方法,其演化方程为

其中,f(x,t)为节点x在t时刻的密度分布函数;m和m(eq)分别为f和f(eq)对应的矩空间分布函数,其中f(eq)为密度平衡态分布函数;e为离散速度方向矢量;S=diag(s0,s1,s2,s3,s4,s5,s6,s7,s8)为松弛系数矩阵;δt为时间步长;M为转换矩阵;I为单位矩阵;F为外力项.对于D2Q9离散速度模型,离散速度e为

转换矩阵M为[23,25]

通过转换矩阵可将粒子密度分布函数投影到速度矩空间.对于D2Q9模型,矩空间分布函数及其平衡态分布函数的计算如下：

与现有的正交MRT-LB模型中通过Gram-Schmidt正交化方法得到的转换矩阵不同,本文模型中采用的转换矩阵的获得无需采用Gram-Schmidt正交化过程[23-25],相应的MRT-LB模型可称作非正交MRT-LB模型.显然,与现有的转换矩阵相比,本文模型中的矩阵具有更多的零元素,有望获得更高的计算效率.

低Ma数条件下,平衡态分布函数可近似为

其中,w0=4/9,w1—4=1/9,w5—8=1/36,w5—8=1/36.为格子声速,R和T分别为气体常数和温度.

外力项F可显式的表示为

MRT-LB方法中分布函数的演化与LBGK模型类似,同样可分为碰撞和迁移两步,其中碰撞在矩空间进行,而迁移在速度空间进行：

碰撞

迁移

其中,碰撞后的分布函数f∗=M-1m∗,其中m∗表示碰撞后的矩空间分布函数.

获得密度分布函数后,宏观量密度ρ、速度u和黏性系数ν可通过密度分布函数计算：

Chapman-Enskog多尺度分析表明,此模型可以准确恢复对应的宏观柱坐标下的质量和动量守恒方程(见附录A).

2.3 求解温度场的MRT-LB模型

在LB方法中,二维的温度场可以用D2Q4,D2Q5或者D2Q9离散速度求解.一般来说,采用D2Q4和D2Q5等计算量相对较小的离散速度模型足够对标量场进行准确模拟.本文采用D2Q5离散速度对温度场进行离散.

温度分布函数的演化方程为

其中,h为温度分布函数,n和n(eq)分别为h和h(eq)对应的矩空间分布函数,其中h(eq)为温度平衡态分布函数,Ψ为源项,N转换矩阵,Θ=diag(σ0,σ1,σ2,σ3,σ4)-1为松弛系数矩阵.

D2Q5离散速度{ei|i=0,1,···,4}为

转换矩阵N为[25]

平衡态分布函数为

其中,ϖ∈(0,1)为D2Q5离散速度模型中的参数.在本文模型中,ϖ=1/2.

与标准笛卡尔坐标下的能量方程相比,柱坐标下的能量方程等号右边多了一项这一项在本文模型中作为源项处理.为正确恢复宏观柱坐标下的能量方程,源项Ψ采取如下形式：

与求解速度场的MRT-LB模型类似,可通过转换矩阵可将粒子温度分布函数h(eq)投影到矩空间：

在实际应用中,求解温度场的MRT-LB其演化过程与求解速度场的MRT-LB完全类似,即碰撞在矩空间进行,而迁移在速度空间进行.

宏观量计算如下：

通过Chapman-Enskog多尺度展开,可以证明本文模型可以恢复对应的宏观轴对称能量方程(见附录B).

3 数值模拟结果及讨论

本节采用本文所提出的模型对几种典型的轴对称热流动进行数值模拟以验证模型的可行性和准确性,包括热Womersley流动、竖直柱体空腔中的Rayleigh-Bénard流动和竖直环形柱体内的自然对流.计算中各松弛系数选取如下：s0=s1=s2=1.0,s3=s6=1.1,s7=s8=1.2,s4和s5根据(17)式计算.σ0=1.0,σ3=σ4=1.5,σ1和σ2根据(25)式计算.数值模拟中温度场和速度场的MRT-LB模型通过Prandtl数和松弛系数进行耦合：

其中,Pr=ν/k为Prandtl数.固体边界采用非平衡外推方法[26]进行处理.

3.1 热Womersley流动

本节首先采用本文模型对热Womersley流动问题进行模拟.恒定壁温下的Womersley流动是工程中常见的一种流动现象.该流动是由周期性外力驱动的平直管道中的流动.周期性外力ax=Gcos(ωt),其中G为外力幅值,ω=2π/Tp为角频率,Tp为周期.Womersley流动的特征数为Reynolds数和Womersley数,分别定义如下：

其中,R为管道半径,ν为流体黏性,u0=GR2/4ν为特征速度.

Womersley流动的轴向速度的解析解为

其中,i为虚数单位,φ=α(i-1)/2,J0为第一类零阶Bessel函数.

温度场的解析解为

计算中Re和α分别取为1200和8,Pr取为0.1,壁面温度分别取为Tw=-1,1和3,流体初始温度取为0.经过网格无关性验证,计算网格取为Nx×Nr=82×42.流动在经过10个周期后达到相对的稳定状态,此时轴向速度和温度的分布与解析解应当相同.

图1(a)为流动达到相对稳定状态后一个周期内t=Tp/4,Tp/2,3Tp/4和Tp时刻轴向速度剖面图和相应解析解的对比图(为节省篇幅,由流动的轴对称特征,图1(a)只展示了关于轴心r=0对称的1/2区域,后面的图示类似处理).由图1(a)可见,本文模型计算结果与解析解符合得很好.图1(b)显示了壁面温度分别为Tw=-1,1和3,流场初始温度取为0,流动达到相对稳定状态后温度场剖面图与解析解对比情况,其中圆形、正方形和菱形符号分别对应Tw=-1,1和3时的数值解,实线为相应的解析解.可以看出温度场的计算结果也与解析解完全相符.

图1 Womersley流动轴向速度和温度剖面图 (a)不同时刻速度剖面图;(b)不同壁面温度时温度剖面图Fig.1.Axial velocity profiles of Womersley flow：(a)Velocity profiles at different times;(b)temperature profiles at different wall temperatures.

为验证本文模型相对LBGK模型的数值稳定性,将两模型进行了对比.首先采用现有的轴对称LBGK模型[23]求解流场和温度场.参数选取以及计算网格均与本节前面所述相同.对热Womersley流动进行模拟,通过逐步减小黏性系数,发现在黏性系数ν=2.76×10-2时LBGK模型发散.而采用本文模型,当黏性系数ν=2.76×10-2时,模型依然稳定并给出准确结果(见图2(a)).进一步,黏性系数取更小的值时,例如ν=2.0×10-2,本文模型依然稳定(见图2(b)).数值实验表明本文模型相对LBGK模型具有更好的数值稳定性.值得指出的是,虽然热Womersley流动中速度场和温度场是相互独立的,但对于具体流动介质,流体黏性和热扩散系数通过Pr数相互关联(例如本例中Pr=ν/k=0.1),由于求解速度场的MRT-LB模型和求解温度场的MRT-LB模型相互耦合,对于速度场数值稳定性的验证也间接验证了温度场模型的数值稳定性.

图2 低黏性系数时本文模型所得Womersley流动轴向速度剖面图 (a)ν=2.76×10-2;(b)ν=2.0×10-2Fig.2.Axial velocity profiles of Womersley flow at low viscosities：(a)ν=2.76×10-2;(b)ν=2.0×10-2.

3.2 竖直圆柱腔体中的Rayleigh-Bénard流动

Rayleigh-Bénard流动是常见的一种柱体内的自然对流现象,许多文献对此进行过数值研究[18,27].自然对流的特征数为Ra数以及Pr数：

其中,g为重力加速度,β为热膨胀系数,L为特征尺寸.本文引入Boussinesq假设,热浮升力为ρa=-ρg(T-Tr),其中Tr为参考温度.

在竖直圆柱腔体中的Rayleigh-Bénard自然对流问题中,充满空气的柱体腔体高为H,半径为R,且H/R=1.柱体腔体上下面水平恒温,温度分别为Tc和Th,且Tc＜Th;其余壁面绝热.计算中取Tc=0,Th=1,参考温度Tr=(Ti+To)/2,Pr=0.7,Ra=5000,重力加速度g沿轴向负方向.对于自然对流,根据(17)式和(25)式以及Ra,Pr的定义,模型中未知的松弛系数可表示为

其中,L为特征尺寸,本例中取为H,Ma为Mach数,为满足不可压流动假设,本文取Ma=0.1.

计算中收敛标准定义为

首先取流场的初始温度为0,经过网格无关性验证,计算网格取为200×100.图3(a)和图3(b)分别显示了采用本文模型计算的等温线和流线图(为节省篇幅,根据轴对称特性只显示了1/2区域).在流场温度初始温度为0时,腔体底部温度高于周围介质温度,在重力作用下,空气由腔体底部沿腔体中轴线向上运动,而在圆周壁面附近由上部向下部运动,形成稳定的涡.然后取流场的初始温度为1,计算网格同样取为200×100.图4(a)和图4(b)分别显示了流动达到稳定状态时的等温线以及流线图.可以看出,此时的流动特征正好与初始温度为0时相反,即空气沿腔体中轴线由上向下运动,而在圆周壁面附近由下向上运动.这些流动特征与文献[18]以及文献[27]完全相符.

为进一步对算法进行定量考核,表1给出了两种情况下沿中轴线上的速度最大值,并与文献[18]以及文献[27]的计算结果进行了对比,其中无量纲速度定义为由表1可以看出采用本文模型计算所得结果与现有文献符合得很好.

表1 本文模型获得的中轴线上最大速度与已有研究结果的对比Table 1.Comparisons of the maximum velocity in the present study with those of the previous studies.

图3 初始温度为0时的等温线和流线 (a)等温线;(b)流线Fig.3.Isotherms and streamlines as the initial temperature is 0：(a)Isotherms;(b)streamlines.

图4 初始温度为1时的等温线和流线 (a)等温线;(b)流线Fig.4.Isotherms and streamlines as the initial temperature is 1：(a)Isotherms;(b)streamlines.

3.3 竖直环形柱体的自然对流

本小节采用本文模型对竖直环形柱体内的自然对流进行数值模拟.环形柱体内的自然对流是一个广泛研究的问题[17,28-30].流动区域由高度为H的内外两个竖直同心圆柱之间的封闭区域构成,内外圆柱的半径分别为ri和ro,其中ro=2ri,且H/ri=2.内外圆柱壁面温度分别为Ti和To(Ti＞To),其他两个水平边界均为绝热边界.计算中分别取Ra=103,104和105,经过网格无关性验证,对应的计算网格为100×200,200×400以及300×600.其他参数的选择如下：To=0,Ti=1,参考温度Tr=(Ti+To)/2,Pr=0.712,重力加速度g沿轴向负方向.竖直环形柱体自然对流的流动特征数、松弛系数的计算以及收敛标准与3.2节相同,不同的是特征尺寸L=ro-ri.

图5给出了Ra=103,104和105时的等温线图.由图5可见,当Ra较小时,腔体内以热传导为主,表现为等温线在腔体内几乎均匀分布且弱对流作用导致的等温线变形不大.随着Ra增大,对流开始起主要作用,导致腔体中央的等温线扭曲程度变大.从图6的流线也可以看出,随着Ra增大,对流作用增强,流动呈现出更加复杂的特征.这些等温线图和流线图均与文献[17,28—30]符合得很好.

为与前人研究结果做定量比较,本文计算了内圆柱壁面的平均Nu数,其中Nu数定义为表2为采用本文模型获得的内圆柱壁面的Nu数与已有研究结果的对比,可以看出符合良好,证明了本文模型的准确性.

图5 不同Ra数时的等温线图 (a)Ra=103;(b)Ra=104;(c)Ra=105Fig.5.Isotherms at differentRanumbers：(a)Ra=103;(b)Ra=104;(c)Ra=105.

图6 不同Ra数时的流线图 (a)Ra=103;(b)Ra=104;(c)Ra=105Fig.6.Streamlines at differentRanumbers：(a)Ra=103;(b)Ra=104;(c)Ra=105.

表2 本文模型获得的内圆柱壁面的Nu数与已有研究结果的对比Table 2.Comparisons of the average Nusselt numbers in the present study with those of the previous studies.

由于本文模型采用了非正交转换矩阵,前文提到,相对于现有的正交转换矩阵,非正交转换矩阵具有更多的零元素,理论上可以获得更高的数值计算效率.本节采用竖直柱体环中的自然对流算例对本文模型和现有的基于正交矩阵的轴对称MRTLB模型[22](其中,温度场的求解采用文献[18]模型的MRT格式)进行对比.计算中Ra=103,104和105的网格尺度分别为100×200,200×400以及300×600.计算在同一台3.40 GHz主频的计算机上进行.表3列出了现有模型和本文模型计算耗时的对比.由表3可见,对于网格数较少的情况,本文模型计算效率提高了3.5%.然而随着计算网格量的增大,本文模型效率提高幅度变大,对于300×600的计算网格,相对现有模型,本文模型效率提高了4.6%.可以预期,对于对计算网格量要求更高的复杂流动,本文模型将获得更高的效率,特别是对于计算周期长的数值模拟,计算时间的缩短将特别明显.

表3 本文模型与现有模型的效率对比Table 3.Comparisons of the computation efficiency of the present model with that of the existing model.

4 结 论

本文提出了一种轴对称热流动的MRT-LB模型,模型的速度场采用D2Q9离散速度求解,温度场采用D2Q5离散速度求解.本文模型的转换矩阵采用非正交转换矩阵,相比现有MRT-LB模型的正交矩阵具有更多的零元素,因此本文模型相比现有的轴对称MRT-LB模型具有更高的计算效率.Chapman-Enskog多尺度分析表明,本文模型可恢复对应的宏观柱坐标下的连续方程、动量方程以及能量方程.通过对热Womersley流动、竖直圆柱腔体中的Rayleigh-Bénard对流和竖直环形柱体内的自然对流进行数值模拟,证明了本文模型的可行性和准确性.低黏性系数条件下的数值实验表明,相对现有轴对称LBGK模型,本文模型具有更好的数值稳定性.且与现有的轴对称MRT-LB模型相比,本文模型具有更高的计算效率.

附录A速度场MRT-LB模型的Chapman-Enskog多尺度分析

为证明本文模型可以正确恢复宏观的柱坐标下的Naiver-Stokes方程,采用Chapman-Enskog多尺度分析方法对求解速度场的MRT-LB模型进行分析.首先对以下变量进行如下形式的多尺度展开:

其中,ε为与Knudsen数同量级的展开系数,对于宏观尺度,0＜ε≪1.

将(A1a)—(A1d)式代入方程(5)中,可得以下不同空间(时间)尺度上的方程:

由(A2b)式可以得到t1时间尺度上的方程:

由(A2c)式可以得到t2时间尺度上的方程:

根据(A4)式第一、二和三行方程,忽略其中O(Ma3)的高阶项,可得

将(A6a)—(A6d)式代入(A4)式第四至第六行方程,可求得:

将(A7a)—(A7c)式代入(A5b)式和(A5c)式,可得t2时间尺度上的方程:

联合t1时间尺度上的方程(A4)与t2时间尺度上的方程(A8a)和(A8b),可得宏观控制方程(1)—(3),其中动力黏性系数以及体积黏性系数分别定义为

附录B温度场MRT-LB模型的Chapman-Enskog多尺度分析

同样对温度分布函数进行与(A1a)—(A1d)式类似的多尺度展开,并将展开结果代入方程(18)中,可得以下不同空间(时间)尺度上的方程:

其中,D1=∂t1I+Eα∂α1,Eα=N(eiαI)N-1,Θ′=

矩阵E=(Ex,Er)可显式地表达为[25]

由(B1b)式可以得到t1时间尺度上的方程:

其中,Ψ01—Ψ41分别表示t1时间尺度上的源项Ψ1=(0,0,kTσ2,0,0)T的对应分量.

由(B1c)式可得t2时间尺度上的方程:

由(B3)的第二行和第三行可得

对于不可压流动,(B5a)和(B5b)式的时间相关项可以忽略[25],所以(B5a)和(B5b)式可进一步简化为

将(B6a)和(B6b)式代入(B4)式,并考虑到(Ψ11,Ψ12)T=(0,kTσ2)T,可得:

结合t1时间尺度上的方程(B3)的第一行方程与t2时间尺度上的方程(B7),并考虑到柱坐标下的不可压流动的连续方程,可得轴对称温度场的宏观控制方程(4),其中热扩散系数定义为

[1]Vynnycky M,Maeno N 2012Int.J.Heat Mass Transfer55 7297

[2]Grosan T,Pop I 2011Int.J.Heat Mass Transfer54 3139

[3]Huang H,Hong N,Liang H,Shi B C,Chai Z H 2016Acta Phys.Sin.65 084702(in Chinese)[黄虎,洪宁,梁宏,施保昌,柴振华2016物理学报65 084702]

[4]Aidun C K,Clausen J R 2009Annu.Rev.Fluid Mech.42 439

[5]Li Q,Luo K H,Kang Q J,He Y L,Chen Q,Liu Q 2015Prog.Energy Combust.Sci.52 62

[6]Lallemand P,Luo L S 2000Phys.Rev.E61 6546

[7]d’Humières D,Ginzburg I,Krafczyk M,Lallemand P,Luo L S 2002Philos.Trans.R.Soc.London A360 437

[8]Halliady I,Hammond L A,Care C M,Good K,Stevens A 2001Phys.Rev.E64 011208

[9]Lee T S,Huang H,Shu C 2006Int.J.Mod.Phys.C17 645

[10]Reis T,Phillips T N 2007Phys.Rev.E75 056703

[11]Reis T,Phillips T N 2008Phys.Rev.E2008 77 026703

[12]Guo Z L,Han H F,Shi B C,Zheng C G 2009Phys.Rev.E79 046708

[13]Li Q,He Y L,Tang G H,Tao W Q 2010Phys.Rev.E81 056707

[14]Zhou J G 2011Phys.Rev.E84 036704

[15]Peng Y,Shu C,Chew Y T,Qiu J 2003J.Comput.Phys.186 295

[16]Huang H,Lee T S,Shu C 2007Int.J.Numer.Methods Fluids53 1707

[17]Li Q,He Y L,Tang G H,Tao W Q 2009Phys.Rev.E80 037702

[18]Zheng L,Shi B C,Guo Z L,Zheng C G 2010Comput.Fluids39 945

[19]Meng X,Guo Z L 2015Phys.Rev.E92 043305

[20]Liu Q,He Y L 2015Physica A429 215

[21]Li Q,He Y L,Tang G H,Tao W Q 2011Microfluid.Nanofluid.10 607

[22]Wang L,Guo Z L,Zheng C G 2010Comput.Fluids39 1542

[23]Premnath K N,Banerjee S 2009Phys.Rev.E80 036702

[24]Geier M,Schönherr M,Pasquali A,Krafczyk M 2015Comput.Math.Appl.70 507

[25]Liu Q,He Y L,Li D,Li Q 2016Int.J.Heat Mass Transfer102 1334

[26]Guo Z L,Shi B C,Zheng C G 2002Chin.Phys.B11 366

[27]Lemembre A,Petit J P 1998Int.J.Heat Mass Transfer41 2437

[28]Li L K,Mei R W,Klausner J F 2013Int.J.Heat Mass Transfer67 338

[29]Kumar R,Kalam M A 1991Int.J.Heat Mass Transfer34 513

[30]Venkatachalappa M,Sankar M,Natarajan A A 2001Acta Mech.147 173

PACS:47.11.—j,02.60.Cb DOI:10.7498/aps.66.044701

Non-orthogonal multiple-relaxation-time lattice Boltzmann method for axisymmetric thermal flows∗

Wang Zuo Zhang Jia-Zhong†Wang Heng

(School of Energy and Power Engineering,Xi’an Jiaotong University,Xi’an 710049,China)

3 September 2016;revised manuscript

22 November 2016)

Axisymmetric thermal flows in cylindrical systems are widely encountered in engineering practices.Typically,axisymmetric thermal flows belong in three-dimensional(3D)problems.However,taking advantage of the axisymmetric condition,the 3D axisymmetric flows can be reduced to quasi two-dimensional(2D)problems in the meridian plane,which significantly reduces the computational requirements and avoids treating the curved boundary.In recent years,various 2D lattice Boltzmann(LB)models,including single relaxation time LB(SRT-LB,or LBGK)and multiple relaxation time LB(MRT-LB)models,for axisymmetric thermal flows have been proposed.In the LB community,it is well accepted that the MRT-LB is superior to the LBGK in terms of numerical stability.The existing MRT-LB model for axisymmetric thermal flows are developed based on orthogonal basis vectors obtained from the combination of the lattice velocity components,i.e.,the transform matrix in the existing MRT-LB is an orthogonal one.Unlike the existing MRT-LB model,in this paper,a non-orthogonal multiple-relaxation-time lattice Boltzmann(MRT-LB)method of simulating axisymmetric thermal flows is proposed.In the proposed MRT-LB method,the velocity field is solved by a D2Q9 discrete velocity set while the temperature by a D2Q5 discrete velocity set.The main advantage of the present MRT-LB model is that the transform matrix of the model is a non-orthogonal one,which is comprised of some proper non-orthogonal basis vectors obtained from the combination of the lattice velocity components.The non-orthogonal transform matrix of the present MRT-LB model contains more zero elements than the classical orthogonal transform matrix,and thus the present MRT-LB model is expected to be more efficient than the existing orthogonal-based MRT-LB model.The equilibrium velocity and temperature moments of the present MRT-LB model are expressed by mapping the equilibrium distribution functions onto their moment spaces through using the non-orthogonal transformation matrix.Also the vectors in the forcing term are modified according to the matrix mapping.Through the Chapman-Enskog analysis,it is demonstrated that the macroscopic governing equations in the cylindrical coordinate can be recovered from the present MRT-LB model.Then several numerical tests,including thermal Womersley flow,Rayleigh-Bénard convection in a vertical cylinder and natural convection in a vertical annulus,are conducted to validate the present model.It is found that the present numerical results are in good agreement with the analytical solutions and/or other numerical results reported in the literature.Numerical stability is also tested,and the results suggest that the present MRT model shows better numerical stability than its LBGK counterpart.Moreover,the numerical results also indicate that the present MRT-LB model is more computationally efficient than the existing MRT-LB model for axisymmetric thermal flow.These findings indicate that the present MRT-LB model can serve as a powerful method of computing the axisymmetric thermal flows.

axisymmetric thermal flow,lattice Boltzmann method

:47.11.—j,02.60.Cb

10.7498/aps.66.044701

∗国家重点基础研究发展计划(批准号：2012CB026002)和国家重点科技支撑项目(批准号：2013BAF01B02)资助的课题.

†通信作者.E-mail：jzzhang@mail.xjtu.edu.cn

*Project supported by the National Fundamental Research Program of China(Grant No.2012CB026002)and the National Key Technology Research and Development Program of China(Grant No.2013BAF01B02).

†Corresponding author.E-mail：jzzhang@mail.xjtu.edu.cn