氣固兩相流DUGKS-Lagrange耦合算法

1、氣固兩相流DUGKS-Lagrange耦合算法摘要：離散統(tǒng)一氣體動理論算法(Discrete Unified Gas-Kinetic Scheme , DUGKS)比傳統(tǒng)動理學方法 擁有更豐富的流場時空演化信息，卻很少應(yīng)用于氣固兩相流領(lǐng)域&基于點源顆粒的介觀兩相流模 型，提出一種新的氣固兩相流耦合算法，并用于方腔中單顆粒運動及顆粒群擴散的數(shù)值模擬&采 用DUGKS-Lagrange方法求解氣固兩相流模型，相間作用的計算采用積分法則和差值函數(shù)進行時 間和空間的近似，在此基礎(chǔ)上提出新的氣固兩相流介觀耦合算法&氣固兩相流問題的分析與計算 結(jié)果表明，相較于經(jīng)典的格子Boltzmann方法，DUGKS提

2、供了更具潛力的氣固兩相流計算框架& 關(guān)鍵詞：氣固兩相流離散統(tǒng)一氣體動理學算法顆粒軌道模型氣固耦合數(shù)值模擬*引言多相流是自然界和工程領(lǐng)域常見的物質(zhì)運動形態(tài)，計算流體力學(Computational Fluid Dynamics， CFD)是探索多相運動規(guī)律的基本方法對于氣固兩相流數(shù)值模擬，傳統(tǒng)的CFD方法分為Euler- Euler，Euler-PDF和Euler-Lagrange& Euler-Euler方法采用連續(xù)性假設(shè)，將流體和顆粒都作為連續(xù)相 處理，對模擬克努森數(shù)較低的稠密顆粒兩相流具有較好的效果，但模擬較高克努森數(shù)的稀疏顆粒兩相流 時存在一定的局限性；Euler-PDF方法通過描述顆粒

3、概率密度函數(shù)(Probability Density Function， PDF)來刻畫顆粒運動的整體行為，理論上可以模擬任意克努森數(shù)的氣固兩相流動，但顆粒PDF方程求 解的困難限制了該方法的應(yīng)用(3; Euler-Lagrange方法分別采用場和粒子的方法描述氣固兩相流動，具 有更清晰的物理圖像，在氣固兩相流模擬中應(yīng)用廣泛也。以上3種方法均采用Navier-Stokes方程描述 連續(xù)相的運動，但Navier-Stokes方程受連續(xù)介質(zhì)假設(shè)和半經(jīng)驗本構(gòu)關(guān)系的制約，所以亟需開拓嶄新的 模擬思路和計算方法&近年來，基于氣體動理學理論的介觀流體力學方法為氣固兩相流研究提供了新方案&其中，與格子 Bo

4、ltzmann方法結(jié)合的LBE-LGA方法和LBE-Lagrange方法較為流行& LBE-LGA方法中，氣相部分 使用LBE方法求解，顆粒的運動通過格子氣自動機方法求解50 LBE-Lagrange方法通過粒子分布函 數(shù)得到流場速度和壓力分布，并根據(jù)流體當前狀態(tài)計算氣固兩相的作用，與傳統(tǒng)Euler-Lagrange方法 相比，該方法在處理顆粒相上沒有特別的優(yōu)勢(句。Guo等逍結(jié)合LBM和統(tǒng)一動理學算法的優(yōu)點提 出離散統(tǒng)一氣體動理論格式(Discrete Unified Gas-kinetic Scheme，DUGKS)模擬，在模擬連續(xù)和稀薄 流領(lǐng)域獲得巨大的成功，并推廣到氣固兩相流問題的研究

5、& DUGKS是一種求解Boltzmann方程的 二步格式有限體積方法，需要先求解半個時間步長控制體積邊界處的流體分布，與格子Boltzmann方法 相比，DUGKS具有更好的算法穩(wěn)定性和計算精度(10，同時該方法提供了更為豐富的流場信息&充分 挖掘DUGKS特有的優(yōu)勢，將有望發(fā)展更高精度的氣固兩相耦合方法&本文在DUGKS和Lagrange方法的基本框架下，采用積分法則和差值函數(shù)對相間作用項進行時間和空間的近似，以此為基礎(chǔ)，提出一種新的氣固兩相流耦合算法DUGKS-Lagrange,并用于方腔流動中 單顆粒與顆粒群的運動模擬，為氣固兩相流中顆粒的運動求解提供新的思路。1數(shù)值模擬方法DUGKS

6、-Lagrange耦合算法對氣固兩相流的模擬分為兩部分，一部分是氣相模型和算法，本文采 用Boltzmann方程描述氣相流場，并采用DUGKS進行求解。另外一部分是顆粒相的模型和算法，顆粒 相模型采用Lagrange方程，使用顆粒軌道方法追蹤顆粒軌跡。1.1 氣固兩相流模型氣固兩相流系統(tǒng)是氣體與顆粒相互作用的多場耦合系統(tǒng)，對于顆粒濃度較低的不可壓縮的氣相流 場，可以不考慮顆粒對流場的作用，在忽略外力作用的情況下，使用如下Boltzmann-BGK方程描述氣相 運動：(1)+( D 1 %/-/(1)式中，(為粒子的速度+為以速度(運動的速度分布函數(shù)(為碰撞引起的變化+為粒子的弛豫時間! +如為

7、依賴宏觀物理量的平衡函數(shù)：(2)+ = (,2nIiT)D/2 exP (2)式中，F(xiàn)為氣體常數(shù)T為熱力學溫度,D為空間維度。宏觀物理量流體密度4和速度ue與分布函數(shù)+ 之間滿足：(3)A = (4 )= j3(+d(，3()=(3)(4)顆粒在流體中的受力較為復(fù)雜，常見的有氣體阻力，重力和浮力，虛假質(zhì)量力和壓力梯度力等對 于球形重顆粒(.PpEp),顆粒在流場中的受力以氣體阻力為主，于是顆粒在流體中的運動可以描述為： 也= 匝=1(4)Atg,At+p; p)式中，心和p分別為顆粒位置和速度矢量，“，為顆粒所見流體速度矢量+p=ppdp/(18*+)為顆粒弛豫 時間，修正因子九為:(5)+

8、= # R-pRg = S? | g (5)24*式中：&為粒子的密度Sp為粒子的直徑,*為動力粘度，阻力系數(shù)#d采用Schiller和Nauman公式: d4 (1 + 0. 15Re* 687 ) Rep1 000(6)Cd (6)0. 44Rep.1 000當顆粒雷諾數(shù)R-p很小時,C = 24/Rep,這時稱為斯托克斯阻力系數(shù)。1. 2 DUGKSLagrange 耦合算法對Botzmann方程(1)及Lagrange方程(4)分別求解，即可獲得氣固兩相運動信息。二維問題的計 算網(wǎng)格如圖1(a)所示?？刂企w中心位置為瓦，控制體界面中點位置為X”顆粒所在位置記為X?,初始 時刻顆粒隨機分

9、布在計算區(qū)域中。對于氣相Botzmann-BGK方程(1),本文采用DUGKS求解，通過反 復(fù)求解以下2步算法實現(xiàn)：(7)(8)0+1/2 _I t (7)(8)=Jb+F1 / 2JF+產(chǎn)+1 _ . n+ 1/2 子0c 兒 Zro+1/2J C 2=J c + 2 (c$ W構(gòu)造的函數(shù)可以由式（1）得到相對應(yīng)的關(guān)系:（7, + = 2+ .t/2ln , 3./2 f eq+ _ 2+. +2+.+_（9）f：+_!+0 ,+-3f0式中，S為控制體界面中點處流體分布函數(shù)+0, +的梯度，F(xiàn)+1/2是流體分布函數(shù)微通量，可以表達成與 控制體積K和控制體積表面積S相關(guān)：po+1/2(10)

10、（：）f（! . + 1po+1/2(10)顆粒Lagrange方程（4）是一組常微分方程，假設(shè)顆粒弛豫時間弓為常數(shù)，方程（4）具有如下形式的 解析解11（:up (up (t) = up (0) exp ( . )+ exp ( . ) ) exp (u s(t )dif _p )_p _g 門 0_p /Xp (t) = !p (0) + _p 1 - exp ( ) (up (0) + ? 1 exp ( _ ) (u；(V)dt_G(11)(12) 1-expxp+1_xp + ( 1-expxp+1_xp + (up-u0+1/2+1/2( 13)(14)(15)（a）整體的網(wǎng)格劃分

11、 Xp5OXC為了計算式（11）與式（12）中的積分項，在傳統(tǒng)Euler-Lagrange及LBE-Lagrange方法中，積分處理 時，將積分中的氣相速度近似為當前時刻的當?shù)亓黧w速度）由數(shù)值積分公式可知，用當前時刻的流體速 度作近似只有一階精度，而中間時刻的中點法則近似具有二階精度）由式（7）可知，DUGKS可以提供 半個時間步長控制體積邊界處的流體分布，根據(jù)DUGKS的特點，本文將式（11）與式（12）積分項中的流 體速度近似為:0 Xp XeIITu2u0+1/2 = 1, 1 =0 XpXK1|T式中,-，T,0和；為控制體界面，如圖1（b）所示）于是，顆粒位置與速度的解可以表示為:圖

12、1 DUGKS-Lagrange耦合算法網(wǎng)格系統(tǒng)DUGKS-Lagrange耦合算法的核心是在使用DUGKS方法的過程中會獲得半個時間步長的格點 交界面處的分布函數(shù)，為顆粒相的運動提供基礎(chǔ)）算法的主要步驟如下:（1）初始化X處的分布函數(shù)f0+ *（2）使用插值法就可得到瓦處的分布函數(shù)代;（3）運用式（7）計算半個時間步長變化后的瓦處的分布函數(shù)f0+1/2*（4）運用式（3）恢復(fù)瓦處氣體的宏觀量u0+1/2,計算微通量;運用式(13)一式(15)計算顆粒的速度uf 和位移xf運用式(H)計算氣相在經(jīng)過一個時間步長后的分布函數(shù)f0+ ;反復(fù)求解步驟2步驟6,得到流場演化的速度和密度，直到滿足終止條

13、件為止，終止條件為演 化時間達到若干確定時刻&2模擬仿真與分析對頂蓋驅(qū)動的方腔氣固兩相流12(進行數(shù)值模擬，氣相為不可壓流體，離散相為剛性球形顆粒，先后 分析單顆粒運動與顆粒群的擴散兩種工況&2. 1 單顆粒為了觀察顆粒伴隨流體的運動，使用單個顆粒作為觀察對象&為了與文獻12(中的結(jié)果對比，工況 參數(shù)與文獻12( 一致&方腔為網(wǎng)格數(shù)64X 64的均勻網(wǎng)格，邊長為10 cm,雷諾數(shù)為470,顆粒直徑為 3.0 mm,與流體的密度比為1.21。開始時將粒子靜止放置在方腔內(nèi)，之后以175 mm/s的速度向右拖 動上方托板，左右邊界和下邊界都為反彈邊界&圖2(a)為文獻12(使用Euler-Lagra

14、nge算法的模擬結(jié)果，圖2(b)為文獻13(使用Boltzmann算法 的模擬結(jié)果，圖2(c)為本文的DUGKS-Lagrange算法模擬結(jié)果&由圖2(a)可知，在斯托克斯阻力的運 動下，顆粒隨流體運動，隨著時間的推移，顆粒向方腔的壁面移動，這與現(xiàn)有文獻中的模擬結(jié)果一致&運 用本文提出的DUGKS-Lagrange耦合算法得到的粒子運動軌跡比運用文獻13(格子Botzmann方法 的結(jié)果更加連貫順滑，且與原始結(jié)果吻合得更好&圖2(b)LBE-LGA圖2(b)LBE-LGA方腔流單顆粒運動軌跡模擬結(jié)果對比2.2顆粒群使用DUGKS-Lagrange耦合算法對方腔內(nèi)的顆粒群進行模擬&模擬使用的方腔

15、與單顆粒模擬使 用的方腔一致，網(wǎng)格為64X64的均勻網(wǎng)格，頂蓋的運動速度為1 m/s,方向向左&由于顆粒的弛豫時間 影響因素較多，將影響因素用斯托克斯數(shù)&代替13:(16)si。一(16)StL 1*L式中，9為特征速度,L為特征長度，模擬中選取&為0. 027 78作為粒子的特征參數(shù)12(，粒子的直徑為 0. 8 mm。模擬實驗中，先讓流體運動至穩(wěn)定狀態(tài)，再將粒子隨機均勻放入方腔的指定區(qū)域內(nèi)。為了避免粒子 過早與壁面碰撞，將粒子在方腔中的位置設(shè)置為dX d , 3 1. 0 X 107 , i = 1,2, + ,64 = 1,2, + ,64(17)0 1分別運用DUGKS-Lagrang

16、e算法、Euler-Lagrange算法、LBE-LGA算法13進行模擬，選取與文 獻12一致的無量綱時間，時刻顆粒群的運動位置*的取值分別為2. 50,3. 75和5. 00,得到3種算法 的顆粒群運動圖分別如圖3圖5所示&圖3 . = 2.50時，不同算法的顆粒群運動模擬圖圖4 . = 3.75時，不同算法的顆粒群運動模擬圖圖5 . = 5.00時，不同算法的顆粒群運動模擬圖圖3中，. = 2.50時刻，3種算法模擬的結(jié)果中放置于正方形區(qū)域的顆粒扭曲變形，顆粒群從四角處 開始擴散&圖4中，. = 3.75時刻，DUGKS-Lagrange算法與Euler-Lagrange算法的顆粒群四角旋起比 LBE-LGA算法更為明顯&圖5中，.= 5.00時刻，DUGK手Lagrange算法模擬的結(jié)果中上方兩角處的顆粒 被甩出，但仍能維持一個整體，與文獻13中