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

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

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

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

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

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

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

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

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

9、布在計(jì)算區(qū)域中。對(duì)于氣相Botzmann-BGK方程(1),本文采用DUGKS求解，通過反 復(fù)求解以下2步算法實(shí)現(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）得到相對(duì)應(yīng)的關(guān)系:（7, + = 2+ .t/2ln , 3./2 f eq+ _ 2+. +2+.+_（9）f：+_!+0 ,+-3f0式中，S為控制體界面中點(diǎn)處流體分布函數(shù)+0, +的梯度，F(xiàn)+1/2是流體分布函數(shù)微通量，可以表達(dá)成與 控制體積K和控制體積表面積S相關(guān)：po+1/2(10)

10、（：）f（! . + 1po+1/2(10)顆粒Lagrange方程（4）是一組常微分方程，假設(shè)顆粒弛豫時(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為了計(jì)算式（11）與式（12）中的積分項(xiàng)，在傳統(tǒng)Euler-Lagrange及LBE-Lagrange方法中，積分處理 時(shí)，將積分中的氣相速度近似為當(dāng)前時(shí)刻的當(dāng)?shù)亓黧w速度）由數(shù)值積分公式可知，用當(dāng)前時(shí)刻的流體速 度作近似只有一階精度，而中間時(shí)刻的中點(diǎn)法則近似具有二階精度）由式（7）可知，DUGKS可以提供 半個(gè)時(shí)間步長(zhǎng)控制體積邊界處的流體分布，根據(jù)DUGKS的特點(diǎn)，本文將式（11）與式（12）積分項(xiàng)中的流 體速度近似為: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方法的過程中會(huì)獲得半個(gè)時(shí)間步長(zhǎng)的格點(diǎn) 交界面處的分布函數(shù)，為顆粒相的運(yùn)動(dòng)提供基礎(chǔ)）算法的主要步驟如下:（1）初始化X處的分布函數(shù)f0+ *（2）使用插值法就可得到瓦處的分布函數(shù)代;（3）運(yùn)用式（7）計(jì)算半個(gè)時(shí)間步長(zhǎng)變化后的瓦處的分布函數(shù)f0+1/2*（4）運(yùn)用式（3）恢復(fù)瓦處氣體的宏觀量u0+1/2,計(jì)算微通量;運(yùn)用式(13)一式(15)計(jì)算顆粒的速度uf 和位移xf運(yùn)用式(H)計(jì)算氣相在經(jīng)過一個(gè)時(shí)間步長(zhǎng)后的分布函數(shù)f0+ ;反復(fù)求解步驟2步驟6,得到流場(chǎng)演化的速度和密度，直到滿足終止條

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

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

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

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