賈雨彬 王樹青 朱玉穎

（1.東北石油大學(xué)石油工程學(xué)院；2.大慶油田有限責(zé)任公司采油四廠）

液固流化床反應(yīng)器（LSFB）以其良好的攪拌、傳熱和傳質(zhì)性能，在化工、能源、冶金、材料、醫(yī)藥、食品及環(huán)保等領(lǐng)域得到了廣泛應(yīng)用。 在許多工業(yè)過程中，LSFB中的顆粒往往具有不同的尺寸和密度，這直接影響到顆粒的沉降速度和顆粒的分離、 混合等流動(dòng)行為。 雙組分液固流化床（BLSFB）的分離與混合有均勻混合、非均勻混合、不完全分離和完全分離4種模式［1］，但是其顆?；旌虾头蛛x受到流體速度、顆粒密度、顆粒形狀和大小的影響［2～4］。

隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，采用數(shù)值模擬方法對(duì)二元固液流化床的流動(dòng)進(jìn)行研究已成為學(xué)術(shù)主流。 近年來，許多學(xué)者綜合利用計(jì)算流體動(dòng)力學(xué)（CFD）連續(xù)相模型和離散元模型（DEM）對(duì)雙組分顆粒流化床混合和分離進(jìn)行了數(shù)值模擬研究。Seibert K D和Burns M A研究了液固流化床內(nèi)不同大小與密度的顆粒的流動(dòng)現(xiàn)象［5］。 采用統(tǒng)計(jì)力學(xué)方法捕獲粒子的隨機(jī)運(yùn)動(dòng)，并結(jié)合顆粒體積分?jǐn)?shù)驗(yàn)證了液固流化床內(nèi)雙組分顆粒的彌散、分離和反混。 Renzo A D等利用DEM-CFD方法模擬再現(xiàn)了液固流化床層反演現(xiàn)象，研究了局部的粒子流場(chǎng)，顯示出固體在不斷形成和消失的混亂漩渦中的不規(guī)則運(yùn)動(dòng)是顆?；旌蠙C(jī)制；在強(qiáng)烈不均勻濃度分布的情況下， 流體-顆粒相互作用力沿著床層高度呈現(xiàn)恒定的趨勢(shì)［6］。

筆者采用MFIX-DEM方法對(duì)Galvin K P等的雙組分液固流化床流動(dòng)實(shí)驗(yàn)［7］進(jìn)行數(shù)值模擬。 首先，對(duì)雙組分顆粒在液固流化床內(nèi)的運(yùn)動(dòng)特性進(jìn)行分析， 得到顆粒分離與混合的基本運(yùn)動(dòng)規(guī)律；其次，對(duì)雙組分顆粒在混合與分離過程的受力情況進(jìn)行分析，得到顆粒實(shí)現(xiàn)混合和分離的基本機(jī)制；最后，引入顆粒擬溫度對(duì)顆粒隨機(jī)運(yùn)動(dòng)進(jìn)行分析，得到顆粒隨機(jī)運(yùn)動(dòng)的基本機(jī)制。

1 數(shù)學(xué)模型

1.1 連續(xù)相控制方程

連續(xù)性方程：

動(dòng)量方程：

其中，t是時(shí)間；αl是液相體積分?jǐn)?shù)；ρl是液相密度；Ul是控制方程中液相瞬時(shí)速度的系綜平均；pl是熱力學(xué)壓力；τl是液相應(yīng)力張量；M是固體顆粒相總相數(shù)；Il是液相與顆粒相的動(dòng)量交換相；g是重力加速度；下角i，j僅限于愛因斯坦求和約定。

1.2 顆粒運(yùn)動(dòng)方程

單個(gè)顆粒在流化床中的運(yùn)動(dòng)可以通過牛頓第二定律進(jìn)行描述：

其中，mp、up、Ip、ωp、Tp分別為顆粒p的質(zhì)量、速度、轉(zhuǎn)動(dòng)慣量、角速度、力矩；Fp，drag為顆粒運(yùn)動(dòng)過程中受到周圍流體給的曳力（包含顆粒所受到周圍流場(chǎng)的壓力梯度力）；Fp，coll為顆粒的碰撞力（軟球碰撞模型）［8］。 t時(shí)刻固相中的第m相中顆粒p在k-th網(wǎng)格內(nèi)所受到的曳力為：

其中，▽pl，k是在k-th網(wǎng)格中心的壓力梯度；Vp是顆粒p的體積；βlm是k-th網(wǎng)格內(nèi)局部的曳力系數(shù)；αsm是固相中m-th相顆粒的體積分?jǐn)?shù)；ul（Xk）是k-th網(wǎng)格內(nèi)液相速度；usm（Xk）是網(wǎng)格內(nèi)m-th相顆粒的局部平均速度。 這里的曳力系數(shù)采用適用于多組分的BVK曳力模型［9］，其表達(dá)式如下：

其中，F(xiàn)為無量綱力；Fstokes為stokes力；ym為mth相顆粒的體積分?jǐn)?shù)；m為平均雷諾數(shù)；αs為固相體積分?jǐn)?shù)；μl為液相動(dòng)力粘度系數(shù)；dpm為m-th相顆粒的直徑。

1.3 邊界條件及模擬參數(shù)

筆者對(duì)Galvin K P等的液固流化床內(nèi)雙組分顆粒流動(dòng)實(shí)驗(yàn)［7］進(jìn)行數(shù)值模擬，并將實(shí)驗(yàn)裝置簡(jiǎn)化成0.05m×1.2m的二維矩形。 初始條件下，顆粒在床層內(nèi)均勻堆積，其中輕顆粒在下層，重顆粒在上層。 采用速度入口和壓力出口邊界條件，并且假設(shè)入口處液相速度均勻分布，液相在壁面處采用無滑移邊界條件。 采用MFIX-DEM開源代碼對(duì)BLSFB內(nèi)的顆粒流動(dòng)進(jìn)行數(shù)值模擬， 其中在空間上采用二階精度的superbee離散方法， 在時(shí)間上采用隱式Euler方法。模擬基本參數(shù)見表1。模擬的網(wǎng)格數(shù)為10×240，采用自適應(yīng)時(shí)間步長(zhǎng)，其范圍是10-7～10-3，模擬時(shí)間為60s，如非特殊說明，取后10s的時(shí)均數(shù)據(jù)進(jìn)行分析。

表1 模擬基本參數(shù)

2 計(jì)算結(jié)果

2.1 模型驗(yàn)證

入口速度0.031m/s時(shí)顆粒體積分?jǐn)?shù)的軸向分布曲線與實(shí)驗(yàn)對(duì)比如圖1所示。 由圖1可知，輕、重顆粒體積分?jǐn)?shù)的模擬曲線與實(shí)驗(yàn)數(shù)據(jù)吻合較好，驗(yàn)證了模擬的有效性。

2.2 顆粒運(yùn)動(dòng)基本特性

圖2為入口速度0.031m/s時(shí)顆粒的瞬時(shí)運(yùn)動(dòng)圖，其中紅色圈代表輕顆粒，藍(lán)色圈代表重顆粒。 由圖2可知，顆粒的運(yùn)動(dòng)規(guī)律為：在初始時(shí)刻，輕顆粒均勻填充在底層，重顆粒填充在頂層；當(dāng)入口有流體注入時(shí)，兩種顆粒在流體的帶動(dòng)下一起向上運(yùn)動(dòng)，由于兩種顆粒與流體的密度差不同，而且在壁面處流體速度很小對(duì)顆粒的攜帶能力較弱，使得輕顆粒在床層中心向上運(yùn)動(dòng)，重顆粒沿著壁面向下運(yùn)動(dòng)；顆粒之間的碰撞與接觸使得兩種顆粒進(jìn)行混合；最終在顆粒與流體的密度差作用下使得兩種顆粒發(fā)生分離。

圖1 輕、重顆粒體積分?jǐn)?shù)的模擬曲線與實(shí)驗(yàn)數(shù)據(jù)

圖2 入口速度0.031m/s時(shí)顆粒的瞬時(shí)運(yùn)動(dòng)圖

2.3 顆粒受力分析

圖3為入口速度0.031m/s時(shí)顆粒的受力分布時(shí)均圖。 由圖3可知，除了碰撞力外，與軸向力相比顆粒所受到的徑向力更小，說明顆粒的軸向力對(duì)顆粒運(yùn)動(dòng)起主導(dǎo)作用。 通過顆粒碰撞力的分布可以發(fā)現(xiàn)， 不同方向上顆粒的碰撞力大小相仿，說明顆粒碰撞各向同性。 通過顆粒軸向力的分布可以發(fā)現(xiàn)，輕、重顆粒的曳力存在明顯差距，說明軸向上輕、重顆粒曳力的差異對(duì)顆粒分離起主導(dǎo)作用。

圖3 入口速度0.031m/s時(shí)顆粒的受力分布時(shí)均圖

圖4是入口速度0.031m/s時(shí)顆粒所受到的軸向力隨時(shí)間變化曲線。 由圖4可知，當(dāng)兩種顆粒由初始的分離狀態(tài)進(jìn)入混合狀態(tài)時(shí)，顆粒所受到的各種力均發(fā)生劇烈變化，說明當(dāng)兩種顆?；旌蠒r(shí)顆粒之間以及顆粒與流體之間發(fā)生劇烈的作用。在混合狀態(tài)時(shí)，碰撞力、曳力和壓力梯度力的數(shù)值關(guān)系依次是壓力梯度力大于顆粒間的碰撞力大于顆粒所受到的流體曳力，說明流體的紊亂運(yùn)動(dòng)帶動(dòng)顆粒的不規(guī)則運(yùn)動(dòng)從而引起顆粒之間的隨機(jī)碰撞，此時(shí)顆粒的隨機(jī)碰撞是顆粒在混合狀態(tài)下的主要運(yùn)動(dòng)形式。 在顆粒分離時(shí)，通過力的平均值和標(biāo)準(zhǔn)差可以發(fā)現(xiàn)，與混合狀態(tài)不同的是曳力大于顆粒間的碰撞力，并且輕、重顆粒所受到的曳力存在明顯差異，驗(yàn)證了輕、重顆粒之間的曳力差異是造成兩種顆粒分離的主要因素。 綜合分離和混合過程可以發(fā)現(xiàn)，壓力梯度力在總力中始終處于主導(dǎo)地位，說明流體的紊亂運(yùn)動(dòng)是顆粒隨機(jī)運(yùn)動(dòng)的主要因素。

圖4 入口速度0.031m/s時(shí)顆粒所受到的軸向力隨時(shí)間變化曲線

2.4 顆粒擬溫度

為了更好地分析局部區(qū)域內(nèi)顆粒的隨機(jī)脈動(dòng)，筆者引入顆粒動(dòng)理學(xué)中的顆粒擬溫度對(duì)顆粒的隨機(jī)脈動(dòng)做定量分析。 其中顆粒擬溫度θ的表達(dá)式如下［10］：

其中，N為計(jì)算網(wǎng)格內(nèi)的顆粒數(shù)；u′i為計(jì)算網(wǎng)格內(nèi)第i個(gè)顆粒的脈動(dòng)速度。

圖5為入口速度0.031m/s時(shí)顆粒擬溫度的散點(diǎn)分布圖。 由圖5可知， 重顆粒的顆粒擬溫度較大，兩種顆粒的擬溫度集中區(qū)域?qū)?yīng)的顆粒體積分?jǐn)?shù)與圖1中的顆粒體積分?jǐn)?shù)相當(dāng)， 說明在此工況下顆粒擬溫度主要受顆粒體積分?jǐn)?shù)的制約。 由圖1可知， 重顆粒的體積分?jǐn)?shù)在床層的底層相比輕顆粒更大，此時(shí)顆粒的接觸幾率更大。 另外，由受力分析可知重顆粒受到的力更大，顆粒的隨機(jī)運(yùn)動(dòng)更劇烈，所以重顆粒的顆粒擬溫度更大。

圖5 入口速度0.031m/s時(shí)顆粒擬溫度的散點(diǎn)分布圖

3 結(jié)束語

筆 者 采 用MFIX-DEM 方 法 對(duì)Galvin K P 等的雙組分液固流化床流動(dòng)實(shí)驗(yàn)進(jìn)行數(shù)值模擬。模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合良好， 并得到如下結(jié)果：

a. 在床層內(nèi)顆粒的分離與混合過程是在顆粒與流體的密度差帶動(dòng)下輕顆粒向上運(yùn)動(dòng)，重顆粒向下運(yùn)動(dòng)，最終兩種顆粒實(shí)現(xiàn)分離。

b. 壓力梯度力在顆粒的運(yùn)動(dòng)過程中占主導(dǎo)地位， 是顆粒隨機(jī)運(yùn)動(dòng)的主要因素； 碰撞力在顆粒的混合過程中起主導(dǎo)作用； 軸向曳力在顆粒的分離過程中起主導(dǎo)作用。

c. 顆粒擬溫度受顆粒的體積分?jǐn)?shù)制約，與輕顆粒相比重顆粒的顆粒擬溫度相對(duì)較大。