氣固兩相流中顆粒速度分布函數(shù)統(tǒng)計(jì)分析

王 婧 ，王軍武

（1.中國科學(xué)院過程工程研究所，多相復(fù)雜系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 100190；2.中國科學(xué)院大學(xué)化學(xué)工程學(xué)院，北京 100049）

氣固兩相流在工業(yè)中有廣泛的應(yīng)用，因此，研究氣固兩相流的流動機(jī)制對于反應(yīng)器放大和設(shè)計(jì)具有重要的意義。近年來，計(jì)算機(jī)模擬已經(jīng)成為研究氣固兩相流的主要手段之一，其中連續(xù)介質(zhì)模型經(jīng)常應(yīng)用于工業(yè)級反應(yīng)器的模擬[1-2]，它通常采用顆粒動理論封閉固相應(yīng)力[2]。

在顆粒動理論中，顆粒速度分布函數(shù)是核心。而顆粒動理論關(guān)系到連續(xù)模型中本構(gòu)關(guān)系的封閉，因此確定顆粒速度分布函數(shù)的形式是建立正確本構(gòu)關(guān)系的關(guān)鍵。在氣固兩相流中，顆粒-流體間的黏性作用、顆粒-顆粒以及顆粒-固體壁面間的非彈性碰撞引起的能量耗散，使得系統(tǒng)處于非平衡狀態(tài)，導(dǎo)致顆粒速度分布有別于常規(guī)氣體中的麥克斯韋分布。研究表明，在氣固兩相流中，顆粒速度分布函數(shù)可以呈現(xiàn)麥克斯韋分布[3]、指數(shù)分布[4-5]、雙峰分布[6-8]、t分布[9]等形式。在循環(huán)流態(tài)床的頂端等顆粒濃度較低的地方，顆粒速度分布函數(shù)基本滿足麥克斯韋分布，而在環(huán)核結(jié)構(gòu)界面等非均勻結(jié)構(gòu)起主導(dǎo)的位置，顆粒速度函數(shù)一般為指數(shù)分布、t分布或雙峰分布。不同的速度分布會影響到顆粒相本構(gòu)關(guān)系的建立，但是從文獻(xiàn)調(diào)研中我們發(fā)現(xiàn)，顆粒速度分布在不同工況下具有不同的函數(shù)形式，為了完善更加全面的顆粒動理論，我們將分析這幾種速度分布函數(shù)在顆粒流體系統(tǒng)中的普適性。

利用基于連續(xù)介質(zhì)的計(jì)算流體力學(xué)方法（CFD）和離散單元法（DEM）耦合在一起，通過CFD-DEM耦合法模擬研究氣固流化床中粗顆粒的速度分布函數(shù)形式，并與以上各種分布函數(shù)進(jìn)行比較，進(jìn)而找到非平衡系統(tǒng)中相對合理的顆粒速度分布函數(shù)形式，為顆粒動理論的發(fā)展提供支撐。

1 數(shù)值實(shí)驗(yàn)

1.1 CFD-DEM方法

大量的研究表明，CFD-DEM方法可以模擬氣固流化床內(nèi)的流體力學(xué)行為，特別是密相粗顆粒流態(tài)化的行為[10-11]，因此，直接采用CFD-DEM方法，不對CFD-DEM方法本身的合理性進(jìn)行論證。

CFD-DEM方法[12-13]采用平均化的Navier-Stokes方程描述氣相的運(yùn)動，采用牛頓第二定律描述顆粒相的運(yùn)動。它在處理顆粒和顆粒之間的相互作用時(shí)，采用“軟球模型”[14-15]，這種處理方式考慮到顆粒碰撞有一定的接觸時(shí)間，通過彈性、阻尼和滑移的力學(xué)特征來描述顆粒碰撞過程的形變、緩沖和滑移的接觸過程。另外，在計(jì)算相間耦合時(shí)，采用Beestra關(guān)聯(lián)式[16]來計(jì)算氣固相間曳力。CFD-DEM模型的控制方程為式（1）—式（10）。

顆粒運(yùn)動受力方程：

式中，ma為顆粒質(zhì)量，ra為顆粒位置，Va為顆粒體積，p為壓力，β為相間曳力系數(shù)，εg為氣相濃度，ug為氣相宏觀速度，vs為單個(gè)顆粒的速度，g為重力加速度，F(xiàn)c，a為顆粒碰撞時(shí)受到的接觸力。

兩顆粒間碰撞受力：

式中，F(xiàn)ab，n為顆粒碰撞時(shí)法向接觸力，kn為法向剛性系數(shù)，δn為顆粒碰撞時(shí)發(fā)生的形變量，nab為法向單位向量，ηn為法向阻尼系數(shù)，vab，n為法向相對速度，Ra和 Rb分別表示顆粒a和顆粒b的顆粒半徑。

氣相的質(zhì)量守恒方程：

氣相的動量守恒方程：

式中，ρg為氣相密度，τg為氣相應(yīng)力。

氣固相間曳力Sp：

Vc為計(jì)算網(wǎng)格的體積。

相間曳力系數(shù)β：

理想氣體狀態(tài)方程：

式中，ρg為氣體密度，dp為顆粒直徑，Mg為氣體摩爾質(zhì)量，R為摩爾氣體常數(shù)，Tg為氣相溫度。

氣相應(yīng)力本構(gòu)關(guān)系：

式中，μg為氣相的剪切黏度，I表示單位矩陣。

1.2 模擬參數(shù)

采用CFD-DEM模型模擬D類顆粒的鼓泡流態(tài)化、湍動流態(tài)化、循環(huán)流態(tài)化這3種兩相流中常見的流型，模擬采用擬三維流化床，即氣相是二維流動，固相是三維流動，床層的厚度是顆粒直徑的6.1倍（采用周期性邊界）。我們將流化床的水平方向設(shè)置為無滑移的壁面，進(jìn)口使用速度邊界，出口使用壓力邊界。D類顆粒流態(tài)化的模擬參數(shù)見表1。

表1 D類顆粒流態(tài)化的模擬參數(shù)Tab.1 Simulation parameters for the fluidization of Geldart D particles

1.3 速度分布函數(shù)

如果顆粒流體系統(tǒng)是均勻流動或系統(tǒng)處在平衡狀態(tài)，顆粒速度分布函數(shù)通常滿足麥克斯韋分布[2]，

式中，c為顆粒微觀速度，f（c）表示顆粒速度分布函數(shù)，θ為顆粒溫度，u為顆粒的平均速度。

但是在氣固兩相流系統(tǒng)中，大多數(shù)情況下，體系處于非平衡態(tài)并且伴隨著熱傳導(dǎo)、擴(kuò)散、黏滯現(xiàn)象等傳遞過程。如果顆粒流體系統(tǒng)中，顆粒處于非平衡狀態(tài)，顆粒的速度分布就會偏離麥克斯韋分布呈現(xiàn)出長拖尾的趨勢，或者呈現(xiàn)多個(gè)麥克斯韋分布疊加的形狀等。

式（12）表示的是顆粒速度分布呈指數(shù)分布的形式：

式中k、α、β均為修正系數(shù)。

式（13）是呈雙峰分布的速度分布函數(shù)，它是基于能量最小多尺度模型（EMMS）思想得到的。

式中，f為密相體積分?jǐn)?shù)，nc為密相數(shù)密度，（1-f）為稀相體積分?jǐn)?shù)，nf為稀相數(shù)密度，且有（1-f）nf+fnc=n，n 是顆粒總的數(shù)密度，θf是稀相顆粒溫度，usf為稀相顆粒的平均速度，θc為密相顆粒的溫度，usc為密相顆粒的平均速度。

式（14）是呈t分布的速度分布函數(shù)，t分布與麥克斯韋分布相比具有長拖尾特征[17]。

這里 υ表示自由度，并且 υ＞2，（θ+Θ）是修正后的顆粒溫度。式（14）中指的是函數(shù)，它的定義如下：

2 結(jié)果與討論

2.1 水平方向的速度分布函數(shù)

通過CFD-DEM模型計(jì)算出顆粒x方向的速度，根據(jù)DEM數(shù)據(jù)統(tǒng)計(jì)得到x方向的顆粒速度分布函數(shù)，然后，將上述的幾種顆粒速度分布函數(shù)和數(shù)值實(shí)驗(yàn)得到的速度分布進(jìn)行對比。圖1是鼓泡流態(tài)化下x方向顆粒速度分布函數(shù)。從圖可以發(fā)現(xiàn)，t分布對顆粒速度函數(shù)的預(yù)測比較接近于模擬數(shù)據(jù)得到的速度分布。

圖1 鼓泡流態(tài)化下x方向顆粒速度分布函數(shù)Fig.1 PDF of x velocity in bubbling fluidization

圖2為湍動流態(tài)化下x方向的顆粒速度分布。而從圖2可以看出，湍動流態(tài)化下x方向的速度分布函數(shù)更接近指數(shù)分布或雙峰分布。與鼓泡流態(tài)化下的速度分布相比，隨著氣速的增加，可以發(fā)現(xiàn)速度分布的區(qū)間增大，在這種狀態(tài)下，指數(shù)分布或雙峰分布可以較為合理地回歸出模擬得到的速度分布函數(shù)。

圖2 湍動流態(tài)化下x方向的顆粒速度分布函數(shù)Fig.2 PDF of x velocity in turbulent fluidization

隨著氣速的繼續(xù)增大，速度分布的區(qū)間明顯增大，速度分布也出現(xiàn)了明顯的長拖尾趨勢。圖3是循環(huán)流態(tài)化下x方向的顆粒速度分布函數(shù)。從圖可以看出，t分布能夠預(yù)測出這種長拖尾特征，但是與模擬結(jié)果仍存在一定的偏差。我們可以看到在顆粒速度達(dá)到3 m/s以上由模擬數(shù)據(jù)得到的速度分布更離散，這是因?yàn)樵谶@個(gè)范圍的顆粒速度出現(xiàn)比較少，統(tǒng)計(jì)樣本數(shù)不足造成的。

圖3 循環(huán)流態(tài)化下x方向顆粒速度分布函數(shù)Fig.3 PDF of x velocity in circulating fluidization

從圖1—3可以發(fā)現(xiàn)，麥克斯韋分布整體上偏離模擬數(shù)據(jù)較多，這是因?yàn)辂溈怂鬼f分布僅僅適用于描述平衡系統(tǒng)中固體顆粒的速度分布函數(shù)，但是氣固流態(tài)化系統(tǒng)是一個(gè)典型的非平衡系統(tǒng)，因此需要用其他形式的分布函數(shù)來描述。

2.2 氣流方向的速度分布函數(shù)

采用2.1節(jié)同樣的方法，我們比較z方向的速度分布函數(shù)，見圖4—6。從圖4、5、6可知：CFD-DEM模擬數(shù)據(jù)統(tǒng)計(jì)出來的速度分布不再是對稱分布，在不光滑的地方較多，這可能與實(shí)驗(yàn)本身的模擬條件有關(guān)。我們模擬的是一種擬三維的微型流化床，這有可能使得模擬出的速度分布呈現(xiàn)明顯的邊界效應(yīng)。

圖4 鼓泡流態(tài)化下z方向顆粒速度分布函數(shù)Fig.4 PDF of z velocity in bubbling fluidization

圖5 湍動流態(tài)化下z方向顆粒速度分布函數(shù)Fig.5 PDF of z velocity in turbulent fluidization

圖6 循環(huán)流態(tài)化下z方向顆粒速度分布函數(shù)Fig.6 PDF of z velocity in circulating fluidization

z方向的速度分布函數(shù)顯示，在氣體流動的方向也就是z方向上氣相和固相的相互作用更加明顯，顆粒和顆粒之間的相互作用也比x方向更加劇烈，因此系統(tǒng)內(nèi)的流動非平衡特征較明顯，這有可能導(dǎo)致速度分布不對稱甚至出現(xiàn)一些波動。

根據(jù)圖中顯示的幾種速度分布函數(shù)對比可以發(fā)現(xiàn)，雙峰分布對于這種遠(yuǎn)離平衡態(tài)的速度分布有較好的預(yù)測，但是在分布的尾端，雙峰分布也不能很好地?cái)M合模擬得到的數(shù)據(jù)，表明基于EMMS模型的雙峰分布仍需要進(jìn)一步完善。

3 結(jié)論與展望

將DEM模型計(jì)算得到的顆粒速度分布和目前提出的各速度分布函數(shù)進(jìn)行比較，發(fā)現(xiàn)麥克斯韋分布不能準(zhǔn)確預(yù)測出系統(tǒng)內(nèi)顆粒速度分布的形貌，顆粒流體系統(tǒng)中顆粒的速度分布有“長拖尾”的特征。基于EMMS模型的雙峰分布，從物理意義上更貼近氣固兩相流的流動機(jī)制，適合于描述遠(yuǎn)離平衡態(tài)下的顆粒速度分布函數(shù)。

顆粒速度分布函數(shù)在不同方向上具有不同的分布特性，水平方向上氣固間作用較弱時(shí)顆粒速度滿足t分布，豎直方向上氣固耦合作用明顯，兩相流動控制機(jī)制協(xié)調(diào)競爭，顆粒速度滿足雙峰分布。本文中的模擬研究對今后顆粒動理論的發(fā)展提出了新的思路，要求顆粒動理論中可以建立不同方向上滿足不同速度分布特性的本構(gòu)關(guān)系，從而得到更加合理的本構(gòu)關(guān)系。