王佐 張家忠 王恒

(西安交通大學(xué)能源與動(dòng)力工程學(xué)院,西安 710049)

非正交多松弛系數(shù)軸對稱熱格子Boltzmann方法?

王佐 張家忠?王恒

(西安交通大學(xué)能源與動(dòng)力工程學(xué)院,西安 710049)

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

提出了一種模擬軸對稱熱流動(dòng)的非正交多松弛系數(shù)格子Boltzmann(MRT-LB)模型.通過采用非正交轉(zhuǎn)換矩陣,建立了基于D2Q9離散速度的求解流場的MRT-LB模型和基于D2Q5離散速度的求解溫度場的MRT-LB模型.Chapman-Enskog分析表明,該模型可以恢復(fù)對應(yīng)的柱坐標(biāo)下的宏觀連續(xù)方程、動(dòng)量方程和能量方程.與現(xiàn)有的軸對稱MRT-LB模型相比,本文采用的非正交轉(zhuǎn)換矩陣中含有更多的零元素,因而具有更高的計(jì)算效率.采用本文模型對幾種典型的軸對稱熱流動(dòng)問題,包括熱Womersley流動(dòng)、豎直圓柱體內(nèi)的Rayleigh-Bénard對流和環(huán)形柱體內(nèi)的自然對流進(jìn)行了數(shù)值模擬,通過等溫線圖和流線圖以及定量數(shù)據(jù)的對比,驗(yàn)證了本文模型的可行性和可靠性.并且數(shù)值模擬結(jié)果表明,相對現(xiàn)有模型,本文模型具有更好的數(shù)值穩(wěn)定性和計(jì)算效率.

軸對稱熱流動(dòng),格子Boltzmann方法

1 引 言

軸對稱流動(dòng)在換熱器、太陽能設(shè)備、電子器件冷卻等工業(yè)領(lǐng)域有著重要應(yīng)用.近年來有許多文獻(xiàn)對軸對稱流動(dòng)換熱現(xiàn)象進(jìn)行了數(shù)值研究[1,2].格子Boltzmann(LB)方法作為一種新興的數(shù)值方法,由于其分子動(dòng)力學(xué)背景,具有簡單高效、邊界處理簡單、并行性強(qiáng)等特點(diǎn),在許多領(lǐng)域獲得了重要應(yīng)用[3-5].其中單松弛系數(shù)模型格子Boltzmann(SRT-LB,也稱作LBGK)模型由于其實(shí)現(xiàn)簡單,具有最為廣泛的應(yīng)用.近年來,具有更加合理的理論基礎(chǔ)的多松弛系數(shù)格子Boltzmann(MRT-LB)模型也取得了飛速發(fā)展[6,7].

截止目前已有大量研究致力于將LB推廣到軸對稱流動(dòng)問題的研究.Halliday等[8]首先基于宏觀柱坐標(biāo)下的Navier-Stokes(NS)方程提出了一種適用于軸對稱流動(dòng)的LB模型.隨后,Lee等[9]發(fā)現(xiàn)Halliday等[8]的模型缺失了一些徑向速度相關(guān)項(xiàng),不能正確恢復(fù)宏觀的柱坐標(biāo)下的NS方程.為克服此缺點(diǎn),Lee等[9]提出了一種新的軸對稱LB模型.此后,Reis和Phillips[10,11]對Lee等[9]的模型進(jìn)行了部分改進(jìn).但這些模型中均存在梯度項(xiàng),對梯度項(xiàng)采用差分方法處理,降低了模型的數(shù)值準(zhǔn)確性和穩(wěn)定性,破壞了LB固有的計(jì)算局部化特性.與前面幾個(gè)基于宏觀柱坐標(biāo)下的NS方程推導(dǎo)的LB模型不同,Guo等[12]直接從描述柱坐標(biāo)流動(dòng)的Boltzmann方程出發(fā),推導(dǎo)出了一種柱坐標(biāo)下的LB模型.此模型不包含復(fù)雜的梯度項(xiàng),且其實(shí)施流程幾乎與標(biāo)準(zhǔn)LB完全相同.Li等[13]從二維笛卡爾坐標(biāo)下的標(biāo)準(zhǔn)LB模型出發(fā),將宏觀柱坐標(biāo)下的NS方程轉(zhuǎn)化為類似標(biāo)準(zhǔn)的笛卡爾坐標(biāo)下的NS方程加入與徑向坐標(biāo)r有關(guān)的多個(gè)源項(xiàng)的形式,并提出了針對此種形式宏觀方程的LB模型,詳細(xì)闡述了從標(biāo)準(zhǔn)LB模型出發(fā)構(gòu)造軸對稱流動(dòng)LB模型的過程.Zhou[14]的模型對Li等[13]的模型進(jìn)行了改進(jìn).

以上模型均針對等溫流場,為發(fā)展非等溫的軸對稱LB模型,Peng等[15]提出了一種采用LB模型求解速度場,而采用二階差分方法求解溫度場的混合模型.然而,Huang等[16]指出,隨著對流作用的增大,對流項(xiàng)采用二階差分方法容易導(dǎo)致發(fā)散,并對此混合模型進(jìn)行了改進(jìn).然而,采用混合模型在很大程度上喪失了LB簡單易行的優(yōu)勢.Li等[17]提出了一種求解溫度場的雙分布LB模型,其基本思想與文獻(xiàn)[13]類似.Zheng等[18]提出了另一種與Guo等[12]的方法類似的更加簡單的LB模型.

上述軸對稱模型均采用LBGK,即單松弛系數(shù)(SRT)碰撞模型.在該模型中,所有的物理量的弛豫時(shí)間均假設(shè)相同,這顯然與實(shí)際的物理過程不符.2000年,Lallemand和Luo[6]提出了一種多松弛系數(shù)LB模型(MRT-LB),在該模型中,通過引入轉(zhuǎn)換矩陣和碰撞矩陣,用不同的松弛系數(shù)來表示流體系統(tǒng)中不同物理量的弛豫.與SRT-LB相比,MRT-LB顯然具有更加合理的理論基礎(chǔ).事實(shí)上,大量的數(shù)值模擬結(jié)果表明,MRT-LB模型相對于SRT-LB模型具有更好的數(shù)值穩(wěn)定性和邊界處理準(zhǔn)確性[5].由于這些優(yōu)點(diǎn),MRT-LB模型先后被應(yīng)用到流動(dòng)傳熱[19]、多孔介質(zhì)流動(dòng)[20]、多相流[5]和微尺度流動(dòng)等[21]領(lǐng)域.這些MRT-LB模型均求解二維笛卡爾坐標(biāo)下流體相關(guān)問題.2010年,Li等[13]和Wang等[22]將MRT-LB模型推廣到軸對稱流動(dòng)領(lǐng)域,幾乎同時(shí)提出了針對柱坐標(biāo)下軸對稱流動(dòng)的不同的MRT-LB模型.

現(xiàn)有的所有MRT-LB模型均建立在正交化轉(zhuǎn)換矩陣的基礎(chǔ)上,該矩陣通過Gram-Schmidt正交化方法獲得,相應(yīng)的這些MRT-LB模型可稱作正交MRT-LB模型.然而,最近Premnath和Banerjee[23]以及Geier等[24]指出,Gram-Schmidt正交化方法并不是獲得轉(zhuǎn)換矩陣的惟一途徑,并且提出了一種適用于二維笛卡爾坐標(biāo)流場的非正交MRTLB模型.隨后,Liu等[25]根據(jù)這一思想,進(jìn)一步提出了一種適用于二維笛卡爾坐標(biāo)下的溫度場的非正交MRT-LB模型.與現(xiàn)有的正交MRT-LB模型相比,非正交的MRT-LB模型中轉(zhuǎn)換矩陣具有更多的零元素,因而計(jì)算效率獲得一定程度上的提高[25].

基于前人的工作,本文提出了一種針對軸對稱熱流動(dòng)的MRT-LB模型,該模型中流場和溫度場均采用非正交轉(zhuǎn)換矩陣,相較于現(xiàn)有的模型,本文模型中的轉(zhuǎn)換矩陣含有更多的零元素,這使得模型在保持?jǐn)?shù)值穩(wěn)定性良好的前提下其演化過程更加高效.

2 雙分布軸對稱MRT-LB模型

2.1 控制方程

柱坐標(biāo)下的流動(dòng)傳熱問題是一個(gè)三維問題,然而,當(dāng)周向速度為零時(shí),流動(dòng)可簡化為二維問題,其控制方程為[12,17]

其中,x和r為分別為軸向和徑向坐標(biāo);ux,ur分別為軸向和徑向速度;ax,ar分別為軸向和徑向方向的外力;ρ為密度,T為溫度,ν為黏性系數(shù),k為熱擴(kuò)散系數(shù).對于任意變量φ,

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

格子Boltzmann方法通過描述具有離散速度的流體粒子分布函數(shù)的變化過程來反映流體粒子的運(yùn)動(dòng)過程.基于非正交轉(zhuǎn)換矩陣,本文提出一種適用于軸對稱流場的多松弛系數(shù)格子Boltzmann方法,其演化方程為

其中,f(x,t)為節(jié)點(diǎn)x在t時(shí)刻的密度分布函數(shù);m和m(eq)分別為f和f(eq)對應(yīng)的矩空間分布函數(shù),其中f(eq)為密度平衡態(tài)分布函數(shù);e為離散速度方向矢量;S=diag(s0,s1,s2,s3,s4,s5,s6,s7,s8)為松弛系數(shù)矩陣;δt為時(shí)間步長;M為轉(zhuǎn)換矩陣;I為單位矩陣;F為外力項(xiàng).對于D2Q9離散速度模型,離散速度e為

轉(zhuǎn)換矩陣M為[23,25]

通過轉(zhuǎn)換矩陣可將粒子密度分布函數(shù)投影到速度矩空間.對于D2Q9模型,矩空間分布函數(shù)及其平衡態(tài)分布函數(shù)的計(jì)算如下：

與現(xiàn)有的正交MRT-LB模型中通過Gram-Schmidt正交化方法得到的轉(zhuǎn)換矩陣不同,本文模型中采用的轉(zhuǎn)換矩陣的獲得無需采用Gram-Schmidt正交化過程[23-25],相應(yīng)的MRT-LB模型可稱作非正交MRT-LB模型.顯然,與現(xiàn)有的轉(zhuǎn)換矩陣相比,本文模型中的矩陣具有更多的零元素,有望獲得更高的計(jì)算效率.

低Ma數(shù)條件下,平衡態(tài)分布函數(shù)可近似為

其中,w0=4/9,w1—4=1/9,w5—8=1/36,w5—8=1/36.為格子聲速,R和T分別為氣體常數(shù)和溫度.

外力項(xiàng)F可顯式的表示為

MRT-LB方法中分布函數(shù)的演化與LBGK模型類似,同樣可分為碰撞和遷移兩步,其中碰撞在矩空間進(jìn)行,而遷移在速度空間進(jìn)行：

碰撞

遷移

其中,碰撞后的分布函數(shù)f?=M-1m?,其中m?表示碰撞后的矩空間分布函數(shù).

獲得密度分布函數(shù)后,宏觀量密度ρ、速度u和黏性系數(shù)ν可通過密度分布函數(shù)計(jì)算：

Chapman-Enskog多尺度分析表明,此模型可以準(zhǔn)確恢復(fù)對應(yīng)的宏觀柱坐標(biāo)下的質(zhì)量和動(dòng)量守恒方程(見附錄A).

2.3 求解溫度場的MRT-LB模型

在LB方法中,二維的溫度場可以用D2Q4,D2Q5或者D2Q9離散速度求解.一般來說,采用D2Q4和D2Q5等計(jì)算量相對較小的離散速度模型足夠?qū)?biāo)量場進(jìn)行準(zhǔn)確模擬.本文采用D2Q5離散速度對溫度場進(jìn)行離散.

溫度分布函數(shù)的演化方程為

其中,h為溫度分布函數(shù),n和n(eq)分別為h和h(eq)對應(yīng)的矩空間分布函數(shù),其中h(eq)為溫度平衡態(tài)分布函數(shù),Ψ為源項(xiàng),N轉(zhuǎn)換矩陣,Θ=diag(σ0,σ1,σ2,σ3,σ4)-1為松弛系數(shù)矩陣.

D2Q5離散速度{ei|i=0,1,···,4}為

轉(zhuǎn)換矩陣N為[25]

平衡態(tài)分布函數(shù)為

其中,?∈(0,1)為D2Q5離散速度模型中的參數(shù).在本文模型中,?=1/2.

與標(biāo)準(zhǔn)笛卡爾坐標(biāo)下的能量方程相比,柱坐標(biāo)下的能量方程等號(hào)右邊多了一項(xiàng)這一項(xiàng)在本文模型中作為源項(xiàng)處理.為正確恢復(fù)宏觀柱坐標(biāo)下的能量方程,源項(xiàng)Ψ采取如下形式：

與求解速度場的MRT-LB模型類似,可通過轉(zhuǎn)換矩陣可將粒子溫度分布函數(shù)h(eq)投影到矩空間：

在實(shí)際應(yīng)用中,求解溫度場的MRT-LB其演化過程與求解速度場的MRT-LB完全類似,即碰撞在矩空間進(jìn)行,而遷移在速度空間進(jìn)行.

宏觀量計(jì)算如下：

通過Chapman-Enskog多尺度展開,可以證明本文模型可以恢復(fù)對應(yīng)的宏觀軸對稱能量方程(見附錄B).

3 數(shù)值模擬結(jié)果及討論

本節(jié)采用本文所提出的模型對幾種典型的軸對稱熱流動(dòng)進(jìn)行數(shù)值模擬以驗(yàn)證模型的可行性和準(zhǔn)確性,包括熱Womersley流動(dòng)、豎直柱體空腔中的Rayleigh-Bénard流動(dòng)和豎直環(huán)形柱體內(nèi)的自然對流.計(jì)算中各松弛系數(shù)選取如下：s0=s1=s2=1.0,s3=s6=1.1,s7=s8=1.2,s4和s5根據(jù)(17)式計(jì)算.σ0=1.0,σ3=σ4=1.5,σ1和σ2根據(jù)(25)式計(jì)算.數(shù)值模擬中溫度場和速度場的MRT-LB模型通過Prandtl數(shù)和松弛系數(shù)進(jìn)行耦合：

其中,Pr=ν/k為Prandtl數(shù).固體邊界采用非平衡外推方法[26]進(jìn)行處理.

3.1 熱Womersley流動(dòng)

本節(jié)首先采用本文模型對熱Womersley流動(dòng)問題進(jìn)行模擬.恒定壁溫下的Womersley流動(dòng)是工程中常見的一種流動(dòng)現(xiàn)象.該流動(dòng)是由周期性外力驅(qū)動(dòng)的平直管道中的流動(dòng).周期性外力ax=Gcos(ωt),其中G為外力幅值,ω=2π/Tp為角頻率,Tp為周期.Womersley流動(dòng)的特征數(shù)為Reynolds數(shù)和Womersley數(shù),分別定義如下：

其中,R為管道半徑,ν為流體黏性,u0=GR2/4ν為特征速度.

Womersley流動(dòng)的軸向速度的解析解為

其中,i為虛數(shù)單位,φ=α(i-1)/2,J0為第一類零階Bessel函數(shù).

溫度場的解析解為

計(jì)算中Re和α分別取為1200和8,Pr取為0.1,壁面溫度分別取為Tw=-1,1和3,流體初始溫度取為0.經(jīng)過網(wǎng)格無關(guān)性驗(yàn)證,計(jì)算網(wǎng)格取為Nx×Nr=82×42.流動(dòng)在經(jīng)過10個(gè)周期后達(dá)到相對的穩(wěn)定狀態(tài),此時(shí)軸向速度和溫度的分布與解析解應(yīng)當(dāng)相同.

圖1(a)為流動(dòng)達(dá)到相對穩(wěn)定狀態(tài)后一個(gè)周期內(nèi)t=Tp/4,Tp/2,3Tp/4和Tp時(shí)刻軸向速度剖面圖和相應(yīng)解析解的對比圖(為節(jié)省篇幅,由流動(dòng)的軸對稱特征,圖1(a)只展示了關(guān)于軸心r=0對稱的1/2區(qū)域,后面的圖示類似處理).由圖1(a)可見,本文模型計(jì)算結(jié)果與解析解符合得很好.圖1(b)顯示了壁面溫度分別為Tw=-1,1和3,流場初始溫度取為0,流動(dòng)達(dá)到相對穩(wěn)定狀態(tài)后溫度場剖面圖與解析解對比情況,其中圓形、正方形和菱形符號(hào)分別對應(yīng)Tw=-1,1和3時(shí)的數(shù)值解,實(shí)線為相應(yīng)的解析解.可以看出溫度場的計(jì)算結(jié)果也與解析解完全相符.

圖1 Womersley流動(dòng)軸向速度和溫度剖面圖 (a)不同時(shí)刻速度剖面圖;(b)不同壁面溫度時(shí)溫度剖面圖Fig.1.Axial velocity profiles of Womersley flow：(a)Velocity profiles at different times;(b)temperature profiles at different wall temperatures.

為驗(yàn)證本文模型相對LBGK模型的數(shù)值穩(wěn)定性,將兩模型進(jìn)行了對比.首先采用現(xiàn)有的軸對稱LBGK模型[23]求解流場和溫度場.參數(shù)選取以及計(jì)算網(wǎng)格均與本節(jié)前面所述相同.對熱Womersley流動(dòng)進(jìn)行模擬,通過逐步減小黏性系數(shù),發(fā)現(xiàn)在黏性系數(shù)ν=2.76×10-2時(shí)LBGK模型發(fā)散.而采用本文模型,當(dāng)黏性系數(shù)ν=2.76×10-2時(shí),模型依然穩(wěn)定并給出準(zhǔn)確結(jié)果(見圖2(a)).進(jìn)一步,黏性系數(shù)取更小的值時(shí),例如ν=2.0×10-2,本文模型依然穩(wěn)定(見圖2(b)).數(shù)值實(shí)驗(yàn)表明本文模型相對LBGK模型具有更好的數(shù)值穩(wěn)定性.值得指出的是,雖然熱Womersley流動(dòng)中速度場和溫度場是相互獨(dú)立的,但對于具體流動(dòng)介質(zhì),流體黏性和熱擴(kuò)散系數(shù)通過Pr數(shù)相互關(guān)聯(lián)(例如本例中Pr=ν/k=0.1),由于求解速度場的MRT-LB模型和求解溫度場的MRT-LB模型相互耦合,對于速度場數(shù)值穩(wěn)定性的驗(yàn)證也間接驗(yàn)證了溫度場模型的數(shù)值穩(wěn)定性.

圖2 低黏性系數(shù)時(shí)本文模型所得Womersley流動(dòng)軸向速度剖面圖 (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流動(dòng)

Rayleigh-Bénard流動(dòng)是常見的一種柱體內(nèi)的自然對流現(xiàn)象,許多文獻(xiàn)對此進(jìn)行過數(shù)值研究[18,27].自然對流的特征數(shù)為Ra數(shù)以及Pr數(shù)：

其中,g為重力加速度,β為熱膨脹系數(shù),L為特征尺寸.本文引入Boussinesq假設(shè),熱浮升力為ρa(bǔ)=-ρg(T-Tr),其中Tr為參考溫度.

在豎直圓柱腔體中的Rayleigh-Bénard自然對流問題中,充滿空氣的柱體腔體高為H,半徑為R,且H/R=1.柱體腔體上下面水平恒溫,溫度分別為Tc和Th,且Tc＜Th;其余壁面絕熱.計(jì)算中取Tc=0,Th=1,參考溫度Tr=(Ti+To)/2,Pr=0.7,Ra=5000,重力加速度g沿軸向負(fù)方向.對于自然對流,根據(jù)(17)式和(25)式以及Ra,Pr的定義,模型中未知的松弛系數(shù)可表示為

其中,L為特征尺寸,本例中取為H,Ma為Mach數(shù),為滿足不可壓流動(dòng)假設(shè),本文取Ma=0.1.

計(jì)算中收斂標(biāo)準(zhǔn)定義為

首先取流場的初始溫度為0,經(jīng)過網(wǎng)格無關(guān)性驗(yàn)證,計(jì)算網(wǎng)格取為200×100.圖3(a)和圖3(b)分別顯示了采用本文模型計(jì)算的等溫線和流線圖(為節(jié)省篇幅,根據(jù)軸對稱特性只顯示了1/2區(qū)域).在流場溫度初始溫度為0時(shí),腔體底部溫度高于周圍介質(zhì)溫度,在重力作用下,空氣由腔體底部沿腔體中軸線向上運(yùn)動(dòng),而在圓周壁面附近由上部向下部運(yùn)動(dòng),形成穩(wěn)定的渦.然后取流場的初始溫度為1,計(jì)算網(wǎng)格同樣取為200×100.圖4(a)和圖4(b)分別顯示了流動(dòng)達(dá)到穩(wěn)定狀態(tài)時(shí)的等溫線以及流線圖.可以看出,此時(shí)的流動(dòng)特征正好與初始溫度為0時(shí)相反,即空氣沿腔體中軸線由上向下運(yùn)動(dòng),而在圓周壁面附近由下向上運(yùn)動(dòng).這些流動(dòng)特征與文獻(xiàn)[18]以及文獻(xiàn)[27]完全相符.

為進(jìn)一步對算法進(jìn)行定量考核,表1給出了兩種情況下沿中軸線上的速度最大值,并與文獻(xiàn)[18]以及文獻(xiàn)[27]的計(jì)算結(jié)果進(jìn)行了對比,其中無量綱速度定義為由表1可以看出采用本文模型計(jì)算所得結(jié)果與現(xiàn)有文獻(xiàn)符合得很好.

表1 本文模型獲得的中軸線上最大速度與已有研究結(jié)果的對比Table 1.Comparisons of the maximum velocity in the present study with those of the previous studies.

圖3 初始溫度為0時(shí)的等溫線和流線 (a)等溫線;(b)流線Fig.3.Isotherms and streamlines as the initial temperature is 0：(a)Isotherms;(b)streamlines.

圖4 初始溫度為1時(shí)的等溫線和流線 (a)等溫線;(b)流線Fig.4.Isotherms and streamlines as the initial temperature is 1：(a)Isotherms;(b)streamlines.

3.3 豎直環(huán)形柱體的自然對流

本小節(jié)采用本文模型對豎直環(huán)形柱體內(nèi)的自然對流進(jìn)行數(shù)值模擬.環(huán)形柱體內(nèi)的自然對流是一個(gè)廣泛研究的問題[17,28-30].流動(dòng)區(qū)域由高度為H的內(nèi)外兩個(gè)豎直同心圓柱之間的封閉區(qū)域構(gòu)成,內(nèi)外圓柱的半徑分別為ri和ro,其中ro=2ri,且H/ri=2.內(nèi)外圓柱壁面溫度分別為Ti和To(Ti＞To),其他兩個(gè)水平邊界均為絕熱邊界.計(jì)算中分別取Ra=103,104和105,經(jīng)過網(wǎng)格無關(guān)性驗(yàn)證,對應(yīng)的計(jì)算網(wǎng)格為100×200,200×400以及300×600.其他參數(shù)的選擇如下：To=0,Ti=1,參考溫度Tr=(Ti+To)/2,Pr=0.712,重力加速度g沿軸向負(fù)方向.豎直環(huán)形柱體自然對流的流動(dòng)特征數(shù)、松弛系數(shù)的計(jì)算以及收斂標(biāo)準(zhǔn)與3.2節(jié)相同,不同的是特征尺寸L=ro-ri.

圖5給出了Ra=103,104和105時(shí)的等溫線圖.由圖5可見,當(dāng)Ra較小時(shí),腔體內(nèi)以熱傳導(dǎo)為主,表現(xiàn)為等溫線在腔體內(nèi)幾乎均勻分布且弱對流作用導(dǎo)致的等溫線變形不大.隨著Ra增大,對流開始起主要作用,導(dǎo)致腔體中央的等溫線扭曲程度變大.從圖6的流線也可以看出,隨著Ra增大,對流作用增強(qiáng),流動(dòng)呈現(xiàn)出更加復(fù)雜的特征.這些等溫線圖和流線圖均與文獻(xiàn)[17,28—30]符合得很好.

為與前人研究結(jié)果做定量比較,本文計(jì)算了內(nèi)圓柱壁面的平均Nu數(shù),其中Nu數(shù)定義為表2為采用本文模型獲得的內(nèi)圓柱壁面的Nu數(shù)與已有研究結(jié)果的對比,可以看出符合良好,證明了本文模型的準(zhǔn)確性.

圖5 不同Ra數(shù)時(shí)的等溫線圖 (a)Ra=103;(b)Ra=104;(c)Ra=105Fig.5.Isotherms at differentRanumbers：(a)Ra=103;(b)Ra=104;(c)Ra=105.

圖6 不同Ra數(shù)時(shí)的流線圖 (a)Ra=103;(b)Ra=104;(c)Ra=105Fig.6.Streamlines at differentRanumbers：(a)Ra=103;(b)Ra=104;(c)Ra=105.

表2 本文模型獲得的內(nèi)圓柱壁面的Nu數(shù)與已有研究結(jié)果的對比Table 2.Comparisons of the average Nusselt numbers in the present study with those of the previous studies.

由于本文模型采用了非正交轉(zhuǎn)換矩陣,前文提到,相對于現(xiàn)有的正交轉(zhuǎn)換矩陣,非正交轉(zhuǎn)換矩陣具有更多的零元素,理論上可以獲得更高的數(shù)值計(jì)算效率.本節(jié)采用豎直柱體環(huán)中的自然對流算例對本文模型和現(xiàn)有的基于正交矩陣的軸對稱MRTLB模型[22](其中,溫度場的求解采用文獻(xiàn)[18]模型的MRT格式)進(jìn)行對比.計(jì)算中Ra=103,104和105的網(wǎng)格尺度分別為100×200,200×400以及300×600.計(jì)算在同一臺(tái)3.40 GHz主頻的計(jì)算機(jī)上進(jìn)行.表3列出了現(xiàn)有模型和本文模型計(jì)算耗時(shí)的對比.由表3可見,對于網(wǎng)格數(shù)較少的情況,本文模型計(jì)算效率提高了3.5%.然而隨著計(jì)算網(wǎng)格量的增大,本文模型效率提高幅度變大,對于300×600的計(jì)算網(wǎng)格,相對現(xiàn)有模型,本文模型效率提高了4.6%.可以預(yù)期,對于對計(jì)算網(wǎng)格量要求更高的復(fù)雜流動(dòng),本文模型將獲得更高的效率,特別是對于計(jì)算周期長的數(shù)值模擬,計(jì)算時(shí)間的縮短將特別明顯.

表3 本文模型與現(xiàn)有模型的效率對比Table 3.Comparisons of the computation efficiency of the present model with that of the existing model.

4 結(jié) 論

本文提出了一種軸對稱熱流動(dòng)的MRT-LB模型,模型的速度場采用D2Q9離散速度求解,溫度場采用D2Q5離散速度求解.本文模型的轉(zhuǎn)換矩陣采用非正交轉(zhuǎn)換矩陣,相比現(xiàn)有MRT-LB模型的正交矩陣具有更多的零元素,因此本文模型相比現(xiàn)有的軸對稱MRT-LB模型具有更高的計(jì)算效率.Chapman-Enskog多尺度分析表明,本文模型可恢復(fù)對應(yīng)的宏觀柱坐標(biāo)下的連續(xù)方程、動(dòng)量方程以及能量方程.通過對熱Womersley流動(dòng)、豎直圓柱腔體中的Rayleigh-Bénard對流和豎直環(huán)形柱體內(nèi)的自然對流進(jìn)行數(shù)值模擬,證明了本文模型的可行性和準(zhǔn)確性.低黏性系數(shù)條件下的數(shù)值實(shí)驗(yàn)表明,相對現(xiàn)有軸對稱LBGK模型,本文模型具有更好的數(shù)值穩(wěn)定性.且與現(xiàn)有的軸對稱MRT-LB模型相比,本文模型具有更高的計(jì)算效率.

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

為證明本文模型可以正確恢復(fù)宏觀的柱坐標(biāo)下的Naiver-Stokes方程,采用Chapman-Enskog多尺度分析方法對求解速度場的MRT-LB模型進(jìn)行分析.首先對以下變量進(jìn)行如下形式的多尺度展開:

其中,ε為與Knudsen數(shù)同量級的展開系數(shù),對于宏觀尺度,0＜ε?1.

將(A1a)—(A1d)式代入方程(5)中,可得以下不同空間(時(shí)間)尺度上的方程:

由(A2b)式可以得到t1時(shí)間尺度上的方程:

由(A2c)式可以得到t2時(shí)間尺度上的方程:

根據(jù)(A4)式第一、二和三行方程,忽略其中O(Ma3)的高階項(xiàng),可得

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

將(A7a)—(A7c)式代入(A5b)式和(A5c)式,可得t2時(shí)間尺度上的方程:

聯(lián)合t1時(shí)間尺度上的方程(A4)與t2時(shí)間尺度上的方程(A8a)和(A8b),可得宏觀控制方程(1)—(3),其中動(dòng)力黏性系數(shù)以及體積黏性系數(shù)分別定義為

附錄B溫度場MRT-LB模型的Chapman-Enskog多尺度分析

同樣對溫度分布函數(shù)進(jìn)行與(A1a)—(A1d)式類似的多尺度展開,并將展開結(jié)果代入方程(18)中,可得以下不同空間(時(shí)間)尺度上的方程:

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

矩陣E=(Ex,Er)可顯式地表達(dá)為[25]

由(B1b)式可以得到t1時(shí)間尺度上的方程:

其中,Ψ01—Ψ41分別表示t1時(shí)間尺度上的源項(xiàng)Ψ1=(0,0,kTσ2,0,0)T的對應(yīng)分量.

由(B1c)式可得t2時(shí)間尺度上的方程:

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

對于不可壓流動(dòng),(B5a)和(B5b)式的時(shí)間相關(guān)項(xiàng)可以忽略[25],所以(B5a)和(B5b)式可進(jìn)一步簡化為

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

結(jié)合t1時(shí)間尺度上的方程(B3)的第一行方程與t2時(shí)間尺度上的方程(B7),并考慮到柱坐標(biāo)下的不可壓流動(dòng)的連續(xù)方程,可得軸對稱溫度場的宏觀控制方程(4),其中熱擴(kuò)散系數(shù)定義為

[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物理學(xué)報(bào)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

?國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(批準(zhǔn)號(hào)：2012CB026002)和國家重點(diǎn)科技支撐項(xiàng)目(批準(zhǔn)號(hào)：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