基于格子Boltzmann方法的方腔內(nèi)自然對(duì)流與換熱的數(shù)值模擬

姚壽廣, 王公利, 程清芳, 周長(zhǎng)江

(江蘇科技大學(xué) 能源與動(dòng)力工程學(xué)院,江蘇 鎮(zhèn)江 212003)

格子Boltzmann方法是近二十年來(lái)發(fā)展起來(lái)的基于介觀層次的數(shù)值模擬方法.它與傳統(tǒng)的數(shù)值模擬計(jì)算方法不同,由于其具有獨(dú)特的介觀粒子背景,所以與其他數(shù)值計(jì)算方法相比,具有許多獨(dú)特的優(yōu)點(diǎn).在數(shù)值模擬計(jì)算中格子Boltzmann計(jì)算方法的演化過程簡(jiǎn)單清晰,程序結(jié)構(gòu)也非常簡(jiǎn)潔,基于這些優(yōu)點(diǎn),格子Boltzmann方法受到了國(guó)內(nèi)外許多專家學(xué)者的廣泛關(guān)注.格子Boltzmann方法不僅可以處理等溫流動(dòng)模型,同時(shí)也可以很好的處理有溫度變化的物理問題.通常將處理對(duì)流換熱問題的格子Boltzmann 模型稱之為熱格子Boltzmann模型[1],包括3種模型:雙分布函數(shù)模型(double distribution function, DDF)、多速模型(multi-speed)和混合模型(hybrid)[2],其中雙分布函數(shù)模型因其良好的數(shù)值穩(wěn)定性和模擬精度,在應(yīng)用中具有一定的優(yōu)勢(shì). 雙分布函數(shù)模型的基本思想是采用密度分布函數(shù)模擬速度場(chǎng),溫度分布函數(shù)則是用來(lái)模擬溫度場(chǎng).文中基于該雙分布函數(shù)模型模擬了方腔內(nèi)自然對(duì)流與換熱問題,其中對(duì)不同瑞利數(shù)Ra的溫度和速度分布進(jìn)行分析,并與傳統(tǒng)的數(shù)值模擬計(jì)算結(jié)果進(jìn)行了比較.

1 熱格子Boltzmann方法

格子Boltzmann模型通常由3個(gè)部分組成[3]:離散速度模型、平衡態(tài)分布函數(shù)以及分布函數(shù)的演化方程.三者的選擇相互獨(dú)立,但又相互影響.構(gòu)造格子Boltzmann模型時(shí),需要選擇合適的平衡態(tài)分布函數(shù),而平衡態(tài)分布函數(shù)的具體形式又與離散速度模型的構(gòu)造有關(guān).離散速度的對(duì)稱性決定了相應(yīng)的格子Boltzmann模型能否還原到所要求解的宏觀方程.因此離散速度模型的構(gòu)造非常重要,離散速度的個(gè)數(shù)太少可能會(huì)導(dǎo)致某些物理量不滿足守恒定律,個(gè)數(shù)太多又可能造成計(jì)算資源的浪費(fèi).

1.1 速度格子Boltzmann方程

采用的九速正方形離散速度模型,即D2Q9模型[4](圖1).在連續(xù)Boltzmann方程的基礎(chǔ)上,通過對(duì)時(shí)間和空間的離散得到模型的平衡態(tài)分布函數(shù).得到的格子Boltzmann演化方程:

圖1 D2Q9模型Fig.1 D2Q9 model

fα(r+eαδt,t+δt)-fα(r,t)=

根據(jù)離散速度模型確定平衡態(tài)分布函數(shù)：[5]

d0，d1,d2為參數(shù),滿足如下關(guān)系式:

1.2 溫度格子Boltzmann方程

在進(jìn)行溫度場(chǎng)的模擬過程中,建立的格子Boltzmann方程采用與速度分布函數(shù)同樣的九速正方形模型,則模型的演化方程為:

Tα(r+eaδt,t+δt)-Ta(r,t)=

式中：T為宏觀溫度.

1.3 模型的耦合

通常在討論流動(dòng)問題時(shí),會(huì)作如下假設(shè)[5]:①忽略流動(dòng)過程中的黏性熱耗散損失;②密度外的其他物性參數(shù)為常數(shù);③對(duì)密度變化僅考慮動(dòng)量方程中與體積力有關(guān)的項(xiàng),其余各項(xiàng)中的密度都作為常數(shù)考慮.

通常將以上假設(shè)稱之為Boussinesq假設(shè),其流體運(yùn)動(dòng)方程組為:

式中：ν,χ分別為運(yùn)動(dòng)粘度系數(shù)和熱擴(kuò)散系數(shù).通過Boussinesq假設(shè)可以將速度演化方程和溫度演化方程耦合為一個(gè)復(fù)合模型,而速度場(chǎng)和溫度場(chǎng)的耦合則是在速度演化方程的右端增加一個(gè)如下的外力項(xiàng)來(lái)實(shí)現(xiàn)的.

根據(jù)平衡態(tài)分布函數(shù)計(jì)算宏觀速度、溫度:

運(yùn)動(dòng)粘度系數(shù)ν和熱擴(kuò)散系數(shù)χ分別為:

1.4 邊界條件的處理

在使用格子Boltzmann方法進(jìn)行模擬計(jì)算時(shí),需要根據(jù)已知的邊界宏觀條件,來(lái)確定邊界節(jié)點(diǎn)上的分布函數(shù).這時(shí)候就需要合適的邊界處理格式.不同的邊界處理格式對(duì)數(shù)值模擬的精度、穩(wěn)定性以及計(jì)算效率都有很大的影響.文中采用非平衡態(tài)外推格式,其基本思路為:將邊界節(jié)點(diǎn)上的分布函數(shù)分解為平衡態(tài)和非平衡態(tài)分布函數(shù)兩部分,其中,平衡態(tài)分布函數(shù)部分根據(jù)邊界條件的定義獲得,而非平衡態(tài)部分則根據(jù)非平衡態(tài)外推來(lái)進(jìn)行確定.

2 方腔內(nèi)自然對(duì)流

圖2 物理模型Fig.2 Physical models

普朗特?cái)?shù)Pr和瑞利數(shù)Ra作為描述自然對(duì)流問題的2個(gè)最基本的無(wú)量綱 參數(shù),在此處的定義為

3 模擬結(jié)果

圖3～5為耦合分布函數(shù)模型的模擬結(jié)果畫出的不同的瑞利數(shù)下自然對(duì)流流線圖和等溫線圖.

a) 自然對(duì)流流線

b) 等溫線圖3 Ra=103自然對(duì)流流線圖和等溫線圖Fig.3 Ra=1 000 velocity vector diagram and the temperature distribution

a) 自然對(duì)流流線

b) 等溫線圖4 Ra=104自然對(duì)流流線圖和等溫線圖Fig.4 Ra=10 000 velocity vector diagram and the temperature distribution

a) 自然對(duì)流流線

b) 等溫線圖5 Ra=105自然對(duì)流流線圖和等溫線圖Fig.5 Ra=100 000 velocity vector diagram and the temperature distribution

從以上數(shù)值結(jié)果可以發(fā)現(xiàn):

當(dāng)Ra較小時(shí),在方腔中間出現(xiàn)近似圓形的渦流[8].熱量的傳遞主要是由于冷熱壁面之間的熱傳導(dǎo)所致的,等溫線接近于垂直.隨著Ra的增大,渦流的形狀從近似圓形逐漸轉(zhuǎn)變?yōu)闄E圓形.熱量傳遞方式也逐漸從熱傳導(dǎo)向?qū)α鲹Q熱轉(zhuǎn)變,對(duì)流換熱作用開始占統(tǒng)治地位.只有靠近冷熱壁面的等溫線保持垂直,方腔中心的等溫線逐漸趨向水平.當(dāng)Ra為105時(shí),在方腔中產(chǎn)生的漩渦偏離方腔中心,且分裂為兩個(gè)渦.這主要因?yàn)殡S著Ra的增加,浮升力也大幅度增加,導(dǎo)致渦流的卷起作用更為強(qiáng)烈.而此時(shí)方腔的等溫線更加接近于水平線.這些數(shù)值模擬的結(jié)果與物理真實(shí)較為一致.

為了更好地印證數(shù)值解的正確性,將文獻(xiàn)[7]所得到冷壁面的平均Nussdt數(shù)Nuave、最大Nusselt數(shù)Numax及其位置yNu、x方向速度分量ux在通過方腔中央的豎直線上的最大值umax及其位置ymax、y方向速度分量uy在通過方腔中央的水平線上的最大值vmax及其位置xmax的結(jié)果作為參考值,與數(shù)值解進(jìn)行了定量對(duì)比,發(fā)現(xiàn)數(shù)值解與參考解基本保持一致,誤差較小(表1).

表1 數(shù)值解與參考值的對(duì)比Table 1 Comparison of numerical solution and a reference value

4 結(jié)論

基于雙分布格子Boltzmann模型,建立了適合于流體流動(dòng)和換熱的熱格子Boltzmann模型,在建立溫度分布函數(shù)時(shí),采用了D2Q9離散速度模型.并以此熱格子Boltzmann模型模擬了方腔內(nèi)自然對(duì)流的形成及其演化.通過模擬計(jì)算發(fā)現(xiàn):

1)由于采用的是不可壓雙分布函數(shù)模型,所以在一定程度上降低了由可壓縮效應(yīng)帶來(lái)的數(shù)據(jù)誤差,提高了模擬的精確程度.

2)新的熱格子Boltzmann模型在離散速度模型的選擇上具有更大的靈活性和適應(yīng)性.

3)數(shù)值模擬的計(jì)算結(jié)果與傳統(tǒng)的計(jì)算結(jié)果具有很好的一致性.通過以上分析可以發(fā)現(xiàn),熱格子Boltzmann模型在處理流體流動(dòng)與傳熱方面存在著獨(dú)特的優(yōu)點(diǎn),同時(shí)也可以證明本文數(shù)值模擬計(jì)算方法和程序是切實(shí)有效的.

參考文獻(xiàn)(References)

[1] Shen X, Chen H. Lattice Boltzmann model for simulation flows with multiple phases and components[J].PhysicalReviewE,1993,47:1815-1819.

[2] Lallemand P, Luo L S. Theory of the lattice Boltzmann method: acoustic and thermal properties in two and three dimensions[J].PhysicalReviewE,2003(68):036706.

[3] Wolf-Gladrow D A. Lattice-gas cellular automata and lattice Boltzmann models:an introduction[C]. New York:Springer,2000.

[4] Qian Y H,Lallemand P.Lattice BGK models for Navier-Stokes equation[J].EurophysicsLett,1992(17):479-484.

[5] Guo Z, Shi B, Zheng C. A coupled lattice BGK model for the Boussinesq equations[J].InternationalJournalforNumericalMethodsinFluids,2002,39:325-342.

[6] 郭照立,鄭楚光,李青,等. 流體動(dòng)力學(xué)的格子Boltzmann方法[M]. 武漢:湖北科學(xué)技術(shù)出版社,2002:105-110.

[7] Markatos N C, Pericleous K A. Lamminar and turbulent natural convection flow in an enclosed cavity[J].InternationalJournalofHeatandMassTransfer,1984，27(5)：755-772.

[8] 陶文銓. 數(shù)值傳熱學(xué)[M].2版.陜西西安:西安交通大學(xué)出版社,2001:95-97.