液固流化床內雙組分顆粒流動數值模擬

賈雨彬 王樹青 朱玉穎

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

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

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

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

1 數學模型

1.1 連續相控制方程

連續性方程：

動量方程：

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

1.2 顆粒運動方程

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

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

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

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

1.3 邊界條件及模擬參數

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

表1 模擬基本參數

2 計算結果

2.1 模型驗證

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

2.2 顆粒運動基本特性

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

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

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

2.3 顆粒受力分析

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

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

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

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

2.4 顆粒擬溫度

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

其中，N為計算網格內的顆粒數；u′i為計算網格內第i個顆粒的脈動速度。

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

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

3 結束語

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

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

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

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