不可壓熱流體的離散統(tǒng)一氣體動理學(xué)簡化算法

蔡宏敏，徐江榮，梅魯浩，王 路

(杭州電子科技大學(xué)能源研究所，浙江 杭州 310018)

0 引 言

流動與熱交換廣泛存在于自然界及工程領(lǐng)域，表現(xiàn)形式多種多樣。從自然界中的大氣循環(huán)、水循環(huán)、風(fēng)霜雨雪的形成，到工業(yè)中微電子的散熱冷卻、航空航天飛行器的航行活動，都與流動和傳熱過程密切相關(guān)。傳統(tǒng)的熱流體力學(xué)研究從宏觀角度出發(fā)，通過求解Navier-Stokes方程來研究熱流動問題，逐漸發(fā)展成一套完善的解決熱流過程問題的可靠方法[1-2]。但是，宏觀方法難以克服連續(xù)性假設(shè)的制約，在稀薄流及微流動領(lǐng)域的局限性日益突出。以格子Boltzmann方法(Lattice Boltzmann Method，LBM)為代表的介觀流體力學(xué)方法從分子底層分析流動與傳熱機(jī)理，為熱流體的數(shù)值研究提供了新的思路[3-4]。傳統(tǒng)的LBM在可壓縮及稀薄流領(lǐng)域的研究也遇到難以突破的瓶頸，拓展一種廣泛適用于連續(xù)流和稀薄流的統(tǒng)一數(shù)值方法仍是一個尚未徹底解決的科學(xué)問題。Guo等[5]提出的離散統(tǒng)一氣體動理學(xué)方法(Discrete Unified Gas Kinetic Scheme，DUGKS)為發(fā)展跨越微觀-宏觀尺度的統(tǒng)一計算流體力學(xué)方法提供了新的可行性方案。DUGKS是一種基于有限體積的新方法，結(jié)合LBM和統(tǒng)一氣體動理學(xué)(Gas Kinetic Scheme，GKS)方法的優(yōu)點，成功應(yīng)用于連續(xù)流和稀薄流的相關(guān)研究，并不斷被發(fā)展完善[6]。Wang等[7]拓展了熱流耦合的DUGKS方法，通過動力學(xué)邊界條件分別獨立求解流場內(nèi)流體的速度和溫度，成功模擬熱流動現(xiàn)象。Yang等[8]使用DUGKS求解器，將相場方法應(yīng)用于兩相流的模擬中，有效地模擬了界面嚴(yán)重變形的幾種情況，并且精準(zhǔn)捕獲了許多微妙的細(xì)節(jié)。Wu等[9]開發(fā)了一個守恒的DUGKS方法用于解決流場多尺度問題，在非結(jié)構(gòu)化粒子速度空間中，宏觀量由介觀氣體分布函數(shù)通量來更新，大大提高了效率。與LBM相比，DUGKS適用于任意努森數(shù)的流動，網(wǎng)格適用靈活，其數(shù)值精度和算法穩(wěn)定性均有所提升[10-11]，但在一個時間步長內(nèi)需要演化多個分布函數(shù)，算法結(jié)構(gòu)和計算效率需要進(jìn)一步優(yōu)化和提升。本文針對不可壓縮熱流體，通過優(yōu)化演化過程中分布函數(shù)微通量的計算方法，提出一個簡化的DUGKS，在不改變數(shù)值精度和穩(wěn)定性的前提下，大幅提高了計算效率。

1 數(shù)值計算方法

本文采用雙分布Boltzmann方程描述熱流體的輸運過程，碰撞項采用Bhatnagar-Gross-Krook(BGK)碰撞模型：

(1)

式中，fα為流體分布函數(shù)，當(dāng)α=1,2時，fα分別表示流場分布函數(shù)和溫度分布函數(shù)，ξ為流場內(nèi)流體粒子速度，τα為松弛時間，feq為Maxwellian平衡態(tài)分布函數(shù)：

(2)

式中，R為氣體常數(shù)，D為空間維度，ρ為宏觀量密度，u為流體速度，T為溫度。通過分布函數(shù)得到宏觀量密度、速度以及溫度：

(3)

方程(1)的求解采用類似于Richtmyer格式[12]的兩步算法，網(wǎng)格節(jié)點分別采用實心與空心點表示，fn,fc分別為網(wǎng)格格點和網(wǎng)格中心點上的流體分布函數(shù)，其結(jié)構(gòu)如圖1所示。

圖1 二維網(wǎng)格示意圖

首先，求解半個時間步長后格點xn處流體分布函數(shù)，碰撞項采用梯形法則，將方程(1)離散為：

(4)

其中，

(5)

(6)

(7)

(8)

進(jìn)一步使用如下離散格式求得下一時刻的分布函數(shù)：

(9)

其中，

(10)

在熱流問題的求解中，需要考慮傳熱與流動的耦合，根據(jù)Boussinesq假設(shè)，流體密度與溫度呈線性關(guān)系：

ρ=ρ0-ρ0β(T-T0)

(11)

式中，ρ0，T0分別為流場中的平均密度和平均溫度，β為熱擴(kuò)散系數(shù)，記重力加速度為g0。在Boussinesq假設(shè)下，重力與外力加速度表示為：

(12)

在演化過程中，外力項可以合并到碰撞項中，于是將流場分布函數(shù)控制方程中的碰撞項修正為：

(13)

碰撞項的修正不改變演化過程，算法將繼續(xù)演化流場內(nèi)分布函數(shù)，直至滿足終止條件。

簡化DUGKS算法的計算步驟如下：

(7)重復(fù)步驟2—步驟6，直到達(dá)到終止條件；

(8)輸出流場內(nèi)的宏觀物理量。

2 數(shù)值模擬結(jié)果與分析

2.1 仿真工況

圖2 自然對流方腔示意圖

2.2 仿真結(jié)果與分析

Ra分別為104，105，106時，采用本文算法以及DUGKS算法分別進(jìn)行二維自然對流模擬，得到流場內(nèi)溫度分布情況分別如圖3—圖5所示。從圖3—圖5可以看出，隨著Ra逐漸增大，溫度場內(nèi)的中部等溫線由豎直逐步轉(zhuǎn)變?yōu)樗剑瑹岫伺c冷端壁面附近的等溫線逐漸密集，即等溫壁面附近的溫度梯度變化越發(fā)明顯，方腔中部出現(xiàn)明顯的溫度分層現(xiàn)象。從圖5可以看出，在等溫壁面與絕熱壁面的交匯處，溫度呈現(xiàn)出先下降后升高的現(xiàn)象，從而導(dǎo)致等溫線出現(xiàn)一個較大的波動。

圖3 Ra=104下的方腔內(nèi)溫度場對比

圖4 Ra=105下的方腔內(nèi)溫度場對比

圖5 Ra=106下的方腔內(nèi)溫度場對比

表1 不同Ra下，不同算法的模擬結(jié)果

在原始的DUGKS中，微通量的計算采用中心法則，需要對控制中心xc及邊界中心xb處的流場及溫度分布函數(shù)進(jìn)行演化，一個時間步長的演化次數(shù)約為：

NDUGKS=imax×jmax+(imax+1)×jmax+imax×(jmax+1)≈3(imax×jmax)

(14)

本文采用梯形法則計算微通量，先后演化控制中心xc與格點xn不再需要計算邊界中心處的流場及溫度分布函數(shù)。一個時間步長的演化次數(shù)約為：

NSDUGKS=imax×jmax+(imax+1)×(jmax+1)≈2(imax×jmax)

(15)

理論上，本文算法的計算量相較于原始算法大約下降了1/3左右。

Ra=106時，分別采用本文算法、DUGKS算法進(jìn)行二維自然對流模擬，得到方腔豎直中線處水平速度、方腔水平中線處的垂直速度以及溫度分布如圖6所示。從圖6可以看出，在流場和溫度場的模擬中，本文算法與DUGKS算法取得的結(jié)果非常吻合。

圖6 水平中線處速度分布和溫度分布比較

為了比較2種算法的演化速度，設(shè)置方腔內(nèi)水平速度殘差E為終止條件，

(16)

圖7 不同算法的演化速度對比

本文算法、DUGKS算法達(dá)到相同終止條件時的演化速度如圖7所示。從圖7可以看出，本文算法的收斂速度明顯快于DUGKS算法，2種算法達(dá)到終止條件所需時間分別為26.80 min和40.20 min。由此可見，本文算法在不改變DUGKS算法穩(wěn)定性的前提下，提高了算法的計算效率，提升約30%。

3 結(jié)束語

針對于不可壓熱流體，本文基于DUGKS提出了一種簡化算法。計算網(wǎng)格內(nèi)微通量時，采用梯形公式代替原始算法中的中點公式，在保持原始算法的數(shù)值穩(wěn)定性和精度下，有效解決了原始算法消耗較多計算資源的問題，提高了算法的計算效率。在低瑞利數(shù)情況下，本文算法的計算效率已取得不錯的成績，后續(xù)將在高瑞利數(shù)、復(fù)雜流體等方向展開進(jìn)一步研究。