液固流化床內(nèi)雙組分顆粒流動數(shù)值模擬

賈雨彬 王樹青 朱玉穎

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

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

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

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

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

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

連續(xù)性方程：

動量方程：

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

1.2 顆粒運動方程

單個顆粒在流化床中的運動可以通過牛頓第二定律進行描述：

其中，mp、up、Ip、ωp、Tp分別為顆粒p的質(zhì)量、速度、轉(zhuǎn)動慣量、角速度、力矩；Fp，drag為顆粒運動過程中受到周圍流體給的曳力（包含顆粒所受到周圍流場的壓力梯度力）；Fp，coll為顆粒的碰撞力（軟球碰撞模型）［8］。 t時刻固相中的第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相顆粒的體積分數(shù)；ul（Xk）是k-th網(wǎng)格內(nèi)液相速度；usm（Xk）是網(wǎng)格內(nèi)m-th相顆粒的局部平均速度。 這里的曳力系數(shù)采用適用于多組分的BVK曳力模型［9］，其表達式如下：

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

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

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

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

2 計算結(jié)果

2.1 模型驗證

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

2.2 顆粒運動基本特性

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

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

圖2 入口速度0.031m/s時顆粒的瞬時運動圖

2.3 顆粒受力分析

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

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

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

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

2.4 顆粒擬溫度

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

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

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

圖5 入口速度0.031m/s時顆粒擬溫度的散點分布圖

3 結(jié)束語

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

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

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

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