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

1、不可壓熱流體的離散統(tǒng)一氣體動(dòng)理學(xué)簡(jiǎn)化算法摘要：離散統(tǒng)一氣體動(dòng)理學(xué)算法是一種適用于連續(xù)流和稀薄流的統(tǒng)一數(shù)值算法，具有良好的數(shù)值 穩(wěn)定性和較高的數(shù)值精度，但算法結(jié)構(gòu)較為復(fù)雜，計(jì)算效率有待提升&針對(duì)不可壓縮熱流體流動(dòng), 在原始的離散統(tǒng)一氣體動(dòng)理學(xué)算法的基礎(chǔ)上，使用梯形法則代替中點(diǎn)法則求解對(duì)流項(xiàng)，提出一個(gè) 簡(jiǎn)化的離散氣體動(dòng)理學(xué)算法，用二維自然對(duì)流進(jìn)行數(shù)值模擬&理論分析和計(jì)算結(jié)果表明，對(duì)于二 維熱對(duì)流問(wèn)題的求解，和原始算法相比，改進(jìn)算法在保障穩(wěn)定性和數(shù)值精度下，計(jì)算效率提升了 30%左右&關(guān)鍵詞：離散氣體動(dòng)理學(xué)算法;介觀流體力學(xué)方法熱流動(dòng)數(shù)值模擬。引言流動(dòng)與熱交換廣泛存在于自然界及工程領(lǐng)域，表現(xiàn)形式多種

2、多樣&從自然界中的大氣循環(huán)、水循 環(huán)、風(fēng)霜雨雪的形成，到工業(yè)中微電子的散熱冷卻、航空航天飛行器的航行活動(dòng)，都與流動(dòng)和傳熱過(guò)程密 切相關(guān)&傳統(tǒng)的熱流體力學(xué)研究從宏觀角度出發(fā)，通過(guò)求解Navier-Stokes方程來(lái)研究熱流動(dòng)問(wèn)題，逐 漸發(fā)展成一套完善的解決熱流過(guò)程問(wèn)題的可靠方法1-2。但是，宏觀方法難以克服連續(xù)性假設(shè)的制約， 在稀薄流及微流動(dòng)領(lǐng)域的局限性日益突出&以格子Boltzmann方法(Lattice Boltzmann Method，LBM) 為代表的介觀流體力學(xué)方法從分子底層分析流動(dòng)與傳熱機(jī)理，為熱流體的數(shù)值研究提供了新的思 路傳統(tǒng)的LBM在可壓縮及稀薄流領(lǐng)域的研究也遇到難以突破的瓶頸

3、，拓展一種廣泛適用于連續(xù) 流和稀薄流的統(tǒng)一數(shù)值方法仍是一個(gè)尚未徹底解決的科學(xué)問(wèn)題& Guo等)*提出的離散統(tǒng)一氣體動(dòng)理學(xué) 方法！ Discrete Unified Gas Kinetic Scheme，DUGKS)D發(fā)展跨越微觀-宏觀尺度的統(tǒng)一計(jì)算流體力學(xué) 方法提供了新的可行性方案& DUGKS是一種基于有限體積的新方法，結(jié)合LBM和統(tǒng)一氣體動(dòng)理學(xué) (Gas Kinetic Scheme，GKS)方法的優(yōu)點(diǎn)，成功應(yīng)用于連續(xù)流和稀薄流的相關(guān)研究，并不斷被發(fā)展完 善句。Wang等)5*拓展了熱流耦合的DUGKS方法，通過(guò)動(dòng)力學(xué)邊界條件分別獨(dú)立求解流場(chǎng)內(nèi)流體的 速度和溫度，成功模擬熱流動(dòng)現(xiàn)象& Ya

4、ng等),*使用DUGKS求解器，將相場(chǎng)方法應(yīng)用于兩相流的模擬 中，有效地模擬了界面嚴(yán)重變形的幾種情況，并且精準(zhǔn)捕獲了許多微妙的細(xì)節(jié)& Wu等逆開(kāi)發(fā)了一個(gè)守 恒的DUGKS方法用于解決流場(chǎng)多尺度問(wèn)題，在非結(jié)構(gòu)化粒子速度空間中，宏觀量由介觀氣體分布函數(shù) 通量來(lái)更新，大大提高了效率&與LBM相比，DUGKS適用于任意努森數(shù)的流動(dòng)，網(wǎng)格適用靈活，其數(shù) 值精度和算法穩(wěn)定性均有所提升10-11，但在一個(gè)時(shí)間步長(zhǎng)內(nèi)需要演化多個(gè)分布函數(shù)，算法結(jié)構(gòu)和計(jì)算 效率需要進(jìn)一步優(yōu)化和提升&本文針對(duì)不可壓縮熱流體，通過(guò)優(yōu)化演化過(guò)程中分布函數(shù)微通量的計(jì)算 方法，提出一個(gè)簡(jiǎn)化的DUGKS，在不改變數(shù)值精度和穩(wěn)定性的前提下，

5、大幅提高了計(jì)算效率&1數(shù)值計(jì)算方法本文采用雙分布Boltzmann方程描述熱流體的輸運(yùn)過(guò)程，碰撞項(xiàng)采用Bhatnagar-Gross-Krook(BGK)碰撞模型:(1)竺1 ：fa-fe(1)式中，fo為流體分布函數(shù)，當(dāng)0=1,2時(shí)，:分別表示流場(chǎng)分布函數(shù)和溫度分布函數(shù)(為流場(chǎng)內(nèi)流體粒 子速度為松弛時(shí)間，feq為Maxwellian平衡態(tài)分布函數(shù):# p/ _ | ( | 2 1(2)一| (2(RT)D/2 exp(2RT ) (2)T/ _ ( |2=9(2(RT)D(2 exp(2RT ) 0 2式中,R為氣體常數(shù),D為空間維度,p為宏觀量密度,為流體速度,T為溫度)通過(guò)分布函數(shù)得到宏

6、觀 量密度、速度以及溫度:(3)p ifl* ,4 = 7：i ,T (3)方程(1)的求解采用類似于Richtmyer格式%12&的兩步算法，網(wǎng)格節(jié)點(diǎn)分別采用實(shí)心與空心點(diǎn)表示, f=,fc分別為網(wǎng)格格點(diǎn)和網(wǎng)格中心點(diǎn)上的流體分布函數(shù)，其結(jié)構(gòu)如圖1所示。格點(diǎn) 中心& 邊界好首先,求解半個(gè)時(shí)間步長(zhǎng)后格點(diǎn)q(.處流體分布函數(shù),碰撞項(xiàng)采用梯形法則，將方程(1)格點(diǎn) 中心& 邊界好首先,求解半個(gè)時(shí)間步長(zhǎng)后格點(diǎn)q(.處流體分布函數(shù),碰撞項(xiàng)采用梯形法則，將方程(1)離散為:丁)= + =()=.()n /J a n /? |【a /(4)其中=+2 rn+4 n+2J a ,ijJ a ,ij& #a ,i

7、j(=) = (7 :) / : = (xh)dS式(6)中，：為控制體表面的外法向量，S為控制體面積,微通量F *的計(jì)算也采用梯形法則，以通過(guò) 的控制體右側(cè)Jn的通量為例，其計(jì)算過(guò)程如下:M L ,n.=+2 打)M+ =( f a H (X+1 h+1 ) + f a , H (X+1 H ) /_J (:) j a, ij ( j ) os72(,+1 H+1+1 h然后，利用分布函數(shù)f=，2在格點(diǎn)處的值，采用梯形法則求出通過(guò)也的控制體在n + 1/2時(shí)刻的微通量:F：，1 (8)= J ()： (QdSSV*進(jìn)一步使用如下離散格式求得下一時(shí)刻的分布函數(shù)：(5)(6)(7)(8)$ n+

8、1/、_ $ + ,= /、廠 n+1 /J a ,j (c)= J a j (8) j V O a, j (8)(9)(10) TOC o 1-5 h z $=+1 _c=C 1 Zc=+1 $ + , = rnCn(10)：a, H J a, H2 Xa，H ! ： a，H J a, H2 H在熱流問(wèn)題的求解中，需要考慮傳熱與流動(dòng)的耦合，根據(jù)Boussmesq假設(shè),流體密度與溫度呈線性關(guān)系： pbpo*po3(T*To)(11)式中,4o,To分別為流場(chǎng)中的平均密度和平均溫度,3為熱擴(kuò)散系數(shù)，記重力加速度為g0。在Boussmesq 假設(shè)下，重力與外力加速度表示為：(12)#G=4go =

9、 &40*403(TTo)g。 a = go3(T*To(12)(13)在演化過(guò)程中，外力項(xiàng)可以合并到碰撞項(xiàng)中，于是將流場(chǎng)分布函數(shù)控制方程中的碰撞項(xiàng)修正為： 。三 *1 J*J1*Ji *1 &/-/eq +；(* 碰撞項(xiàng)的修正不改變演化過(guò)程，算法將繼續(xù)演化流場(chǎng)內(nèi)分布函數(shù)，(13)簡(jiǎn)化DUGKS算法的計(jì)算步驟如下：初始化控制中心處的流場(chǎng)宏觀量，計(jì)算平衡態(tài)Jq,并用其初始化控制中心處的分布函數(shù)J” (采用插值函數(shù)計(jì)算得到格點(diǎn)處的分布函數(shù)fa, H，根據(jù)分布函數(shù)之間對(duì)應(yīng)關(guān)系* + 一 2, a 2 $+32 reqJ a，H = 2,a + max +1) ej max + m ax e (jma

10、x + 1)13 max ej max )(14)本文采用梯形法則計(jì)算微通量，先后演化控制中心兀與格點(diǎn)如 不再需要計(jì)算邊界中心處的流場(chǎng) 及溫度分布函數(shù)。一個(gè)時(shí)間步長(zhǎng)的演化次數(shù)約為：(15)NNSDUGKS max X j max + max + 1 ) X (jmax + 1)12 max ej max) 理論上，本文算法的計(jì)算量相較于原始算法大約下降了 1/3(15)Ra = 106時(shí)，分別采用本文算法DUGKS算法進(jìn)行二維自然對(duì)流模擬，得到方腔豎直中線處水平 速度、方腔水平中線處的垂直速度以及溫度分布如圖6所示#從圖6可以看出，在流場(chǎng)和溫度場(chǎng)的模擬 中，本文算法與DUGKS算法取得的結(jié)果非

11、常吻合0.030.020.010.00-0.01-0.020.10、-K!Mffi8 6 4 2 0 2 4 .OOOOOOO 0.030.020.010.00-0.01-0.020.10、-K!Mffi8 6 4 2 0 2 4 .OOOOOOO O.O.O.1.0X(c)水平方向上的溫度T分布圖6 水平中線處速度分布和溫度分布比較3結(jié)束語(yǔ)f/min圖3結(jié)束語(yǔ)f/min圖7不同算法的演化逮度對(duì)比為了比較2種算法的演化速度，設(shè)置方腔內(nèi)水平速度殘差E為終止條件,(16)E=2 (16)(兀+1 000) | 2 ； II（C=+1(兀+1 000) | 2 ； II式中，誠(chéng)=（k，G =）為方腔內(nèi)的水平速度$在本文計(jì)算中將E取值為106 $本文算法%DUGKS算法達(dá)到相同終止條件 時(shí)的演化速度如圖7所示$從圖7可以看出，本 文算法的收斂速度明顯快于DUGKS算法，2種 算法達(dá)到終止條件所需時(shí)間分別為26. 80 min和 40. 20 min$由此可見(jiàn)，本文算法在不改變DUGKS 算法穩(wěn)定性的前提下，提高