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

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

(杭州電子科技大學(xué)理學(xué)院，浙江 杭州 310018)

0 引 言

多相流是自然界和工程領(lǐng)域常見(jiàn)的物質(zhì)運(yùn)動(dòng)形態(tài)，計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)是探索多相運(yùn)動(dòng)規(guī)律的基本方法[1]。對(duì)于氣固兩相流數(shù)值模擬，傳統(tǒng)的CFD方法分為Euler-Euler，Euler-PDF和Euler-Lagrange。Euler-Euler方法采用連續(xù)性假設(shè)，將流體和顆粒都作為連續(xù)相處理，對(duì)模擬克努森數(shù)較低的稠密顆粒兩相流具有較好的效果，但模擬較高克努森數(shù)的稀疏顆粒兩相流時(shí)存在一定的局限性[2]；Euler-PDF方法通過(guò)描述顆粒概率密度函數(shù)(Probability Density Function，PDF)來(lái)刻畫顆粒運(yùn)動(dòng)的整體行為，理論上可以模擬任意克努森數(shù)的氣固兩相流動(dòng)，但顆粒PDF方程求解的困難限制了該方法的應(yīng)用[3]；Euler-Lagrange方法分別采用場(chǎng)和粒子的方法描述氣固兩相流動(dòng)，具有更清晰的物理圖像，在氣固兩相流模擬中應(yīng)用廣泛[4]。以上3種方法均采用Navier-Stokes方程描述連續(xù)相的運(yùn)動(dòng)，但Navier-Stokes方程受連續(xù)介質(zhì)假設(shè)和半經(jīng)驗(yàn)本構(gòu)關(guān)系的制約，所以亟需開(kāi)拓嶄新的模擬思路和計(jì)算方法。

近年來(lái)，基于氣體動(dòng)理學(xué)理論的介觀流體力學(xué)方法為氣固兩相流研究提供了新方案。其中，與格子Boltzmann方法結(jié)合的LBE-LGA方法和LBE-Lagrange方法較為流行。LBE-LGA方法中，氣相部分使用LBE方法求解，顆粒的運(yùn)動(dòng)通過(guò)格子氣自動(dòng)機(jī)方法求解[5]。LBE-Lagrange方法通過(guò)粒子分布函數(shù)得到流場(chǎng)速度和壓力分布，并根據(jù)流體當(dāng)前狀態(tài)計(jì)算氣固兩相的作用，與傳統(tǒng)Euler-Lagrange方法相比，該方法在處理顆粒相上沒(méi)有特別的優(yōu)勢(shì)[6]。Guo等[7]結(jié)合LBM和統(tǒng)一動(dòng)理學(xué)算法[8]的優(yōu)點(diǎn)提出離散統(tǒng)一氣體動(dòng)理論格式(Discrete Unified Gas-kinetic Scheme, DUGKS)模擬，在模擬連續(xù)和稀薄流領(lǐng)域獲得巨大的成功，并推廣到氣固兩相流問(wèn)題的研究[9]。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-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)

式中，ξ為粒子的速度，f為以速度ξ運(yùn)動(dòng)的速度分布函數(shù)，Ω為碰撞引起的變化，τ為粒子的弛豫時(shí)間，feq為依賴宏觀物理量的平衡函數(shù)：

(2)

式中，R為氣體常數(shù)，T為熱力學(xué)溫度，D為空間維度。宏觀物理量流體密度ρ和速度ug與分布函數(shù)f之間滿足：

(3)

顆粒在流體中的受力較為復(fù)雜，常見(jiàn)的有氣體阻力，重力和浮力，虛假質(zhì)量力和壓力梯度力等。對(duì)于球形重顆粒(ρp?ρ)，顆粒在流場(chǎng)中的受力以氣體阻力為主，于是顆粒在流體中的運(yùn)動(dòng)可以描述為：

(4)

(5)

式中：ρp為粒子的密度，dp為粒子的直徑，μ為動(dòng)力粘度，阻力系數(shù)CD采用Schiller和Nauman公式：

(6)

當(dāng)顆粒雷諾數(shù)Rep很小時(shí)，CD=24/Rep，這時(shí)稱為斯托克斯阻力系數(shù)[8]。

1.2 DUGKS-Lagrange耦合算法

對(duì)Boltzmann方程(1)及Lagrange方程(4)分別求解，即可獲得氣固兩相運(yùn)動(dòng)信息。二維問(wèn)題的計(jì)算網(wǎng)格如圖1(a)所示。控制體中心位置為Xc，控制體界面中點(diǎn)位置為Xb，顆粒所在位置記為Xp，初始時(shí)刻顆粒隨機(jī)分布在計(jì)算區(qū)域中。對(duì)于氣相Boltzmann-BGK方程(1)，本文采用DUGKS求解，通過(guò)反復(fù)求解以下2步算法實(shí)現(xiàn)：

(7)

(8)

構(gòu)造的函數(shù)可以由式(1)得到相對(duì)應(yīng)的關(guān)系：

(9)

(10)

顆粒Lagrange方程(4)是一組常微分方程，假設(shè)顆粒弛豫時(shí)間τp為常數(shù)，方程(4)具有如下形式的解析解[11]：

(11)

(12)

為了計(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)中的流體速度近似為：

(13)

式中，e,w,n和s為控制體界面，如圖1(b)所示。于是，顆粒位置與速度的解可以表示為：

圖1 DUGKS-Lagrange耦合算法網(wǎng)格系統(tǒng)

(14)

(15)

DUGKS-Lagrange耦合算法的核心是在使用DUGKS方法的過(guò)程中會(huì)獲得半個(gè)時(shí)間步長(zhǎng)的格點(diǎn)交界面處的分布函數(shù)，為顆粒相的運(yùn)動(dòng)提供基礎(chǔ)。算法的主要步驟如下：

(7)反復(fù)求解步驟2—步驟6，得到流場(chǎng)演化的速度和密度，直到滿足終止條件為止，終止條件為演化時(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ù)64×64的均勻網(wǎng)格，邊長(zhǎng)為10 cm，雷諾數(shù)為470，顆粒直徑為3.0 mm，與流體的密度比為1.21。開(kāi)始時(shí)將粒子靜止放置在方腔內(nèi)，之后以175 mm/s的速度向右拖動(dòng)上方托板，左右邊界和下邊界都為反彈邊界。

圖2(a)為文獻(xiàn)[12]使用Euler-Lagrange算法的模擬結(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]格子Boltzmann方法的結(jié)果更加連貫順滑，且與原始結(jié)果吻合得更好。

圖2 方腔流單顆粒運(yùn)動(dòng)軌跡模擬結(jié)果對(duì)比

2.2 顆粒群

使用DUGKS-Lagrange耦合算法對(duì)方腔內(nèi)的顆粒群進(jìn)行模擬。模擬使用的方腔與單顆粒模擬使用的方腔一致，網(wǎng)格為64×64的均勻網(wǎng)格，頂蓋的運(yùn)動(dòng)速度為1 m/s，方向向左。由于顆粒的弛豫時(shí)間影響因素較多，將影響因素用斯托克斯數(shù)St代替[13]：

(16)

式中，U為特征速度，L為特征長(zhǎng)度，模擬中選取St為0.027 78作為粒子的特征參數(shù)[12]，粒子的直徑為0.8 mm。

模擬實(shí)驗(yàn)中，先讓流體運(yùn)動(dòng)至穩(wěn)定狀態(tài)，再將粒子隨機(jī)均勻放入方腔的指定區(qū)域內(nèi)。為了避免粒子過(guò)早與壁面碰撞，將粒子在方腔中的位置設(shè)置為x×y∈[0.25,0.75]×[0.25,0.75]。判斷流體是否到達(dá)穩(wěn)定狀態(tài)的條件設(shè)為：

(17)

分別運(yùn)用DUGKS-Lagrange算法、Euler-Lagrange算法[12]、LBE-LGA算法[13]進(jìn)行模擬，選取與文獻(xiàn)[12]一致的無(wú)量綱時(shí)間t時(shí)刻顆粒群的運(yùn)動(dòng)位置，t的取值分別為2.50，3.75和5.00，得到3種算法的顆粒群運(yùn)動(dòng)圖分別如圖3—圖5所示。

圖3 t=2.50時(shí)，不同算法的顆粒群運(yùn)動(dòng)模擬圖

圖4 t=3.75時(shí)，不同算法的顆粒群運(yùn)動(dòng)模擬圖

圖5 t=5.00時(shí)，不同算法的顆粒群運(yùn)動(dòng)模擬圖

圖3中，t=2.50時(shí)刻，3種算法模擬的結(jié)果中放置于正方形區(qū)域的顆粒扭曲變形，顆粒群從四角處開(kāi)始擴(kuò)散。圖4中，t=3.75時(shí)刻，DUGKS-Lagrange算法與Euler-Lagrange算法的顆粒群四角旋起比LBE-LGA算法更為明顯。圖5中，t=5.00時(shí)刻，DUGKS-Lagrange算法模擬的結(jié)果中上方兩角處的顆粒被甩出，但仍能維持一個(gè)整體，與文獻(xiàn)[13]中的結(jié)果大不相同。整體上看，本文的DUGKS-Lagrange算法結(jié)果與文獻(xiàn)[12]的Euler-Lagrange算法結(jié)果基本吻合，優(yōu)于文獻(xiàn)[13]的LBE-LGA算法結(jié)果。

3 結(jié)束語(yǔ)

氣固兩相流中氣相與顆粒相的耦合十分關(guān)鍵，本文運(yùn)用DUGKS方法計(jì)算半個(gè)時(shí)間步長(zhǎng)的分布函數(shù)。在顆粒相耦合過(guò)程中，氣相對(duì)顆粒的作用力計(jì)算使用了半個(gè)時(shí)間步長(zhǎng)的演化信息，更好地模擬了顆粒在流場(chǎng)中的運(yùn)動(dòng)與擴(kuò)散特點(diǎn)。本文提出的DUGKS-Lagrange算法只模擬了氣相對(duì)固相的單相耦合作用，下一步將在氣固兩相流的模擬中加入固相對(duì)氣相的耦合作用。