不可壓熱流體的離散統一氣體動理學簡化算法

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

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

0 引 言

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

1 數值計算方法

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

(1)

式中，fα為流體分布函數，當α=1,2時，fα分別表示流場分布函數和溫度分布函數，ξ為流場內流體粒子速度，τα為松弛時間，feq為Maxwellian平衡態分布函數：

(2)

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

(3)

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

圖1 二維網格示意圖

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

(4)

其中，

(5)

(6)

(7)

(8)

進一步使用如下離散格式求得下一時刻的分布函數：

(9)

其中，

(10)

在熱流問題的求解中，需要考慮傳熱與流動的耦合，根據Boussinesq假設，流體密度與溫度呈線性關系：

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

(11)

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

(12)

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

(13)

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

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

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

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

2 數值模擬結果與分析

2.1 仿真工況

圖2 自然對流方腔示意圖

2.2 仿真結果與分析

Ra分別為104，105，106時，采用本文算法以及DUGKS算法分別進行二維自然對流模擬，得到流場內溫度分布情況分別如圖3—圖5所示。從圖3—圖5可以看出，隨著Ra逐漸增大，溫度場內的中部等溫線由豎直逐步轉變為水平，熱端與冷端壁面附近的等溫線逐漸密集，即等溫壁面附近的溫度梯度變化越發明顯，方腔中部出現明顯的溫度分層現象。從圖5可以看出，在等溫壁面與絕熱壁面的交匯處，溫度呈現出先下降后升高的現象，從而導致等溫線出現一個較大的波動。

圖3 Ra=104下的方腔內溫度場對比

圖4 Ra=105下的方腔內溫度場對比

圖5 Ra=106下的方腔內溫度場對比

表1 不同Ra下，不同算法的模擬結果

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

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

(14)

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

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

(15)

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

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

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

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

(16)

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

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

3 結束語

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