周連勇， 趙永志

(浙江大學 能源工程學院 化工機械研究所，杭州 310027)

1 引 言

噴動床作為一種特殊的流化床，具有結(jié)構(gòu)簡單和接觸效率高等優(yōu)點[1,2]，已經(jīng)在農(nóng)作物干燥、生物制藥以及能源化工等領域得到廣泛應用[3-5]。噴動床的工作效率主要由其內(nèi)部復雜的氣固兩相運動決定，而常規(guī)的實驗方法難以深入探究氣固兩相運動的機理，無法系統(tǒng)地掌握流動結(jié)構(gòu)和動力學參數(shù)。近年來，隨著計算機性能的進步和數(shù)值模型的不斷發(fā)展，使用計算機數(shù)值模型已經(jīng)成為研究噴動床的重要手段，實驗研究中存留的部分問題也通過數(shù)值模擬得到了很好的解決。

目前，噴動床內(nèi)氣固兩相運動的數(shù)值模擬方法大致可以分類兩類,一類是以雙流體模型(two-fluid model)為代表的歐拉-歐拉(Euler-Euler)方法[6-8]，該方法把固體顆粒相和氣相都看作連續(xù)相；另一類是以CFD -DEM(computational fluid dynamics coupled with discrete element method)為代表的歐拉-拉格朗日(Euler-Lagrangian)方法[9-12]，該方法把固體顆粒相看做離散相，把氣相看作連續(xù)相。相較于雙流體模型，CFD -DEM可以追蹤單個顆粒運動，獲得顆粒尺度的微觀信息，在噴動床內(nèi)部氣固兩相運動機理的探究以及噴動床結(jié)構(gòu)的優(yōu)化設計中逐漸成為主流。根據(jù)顆粒周圍流場的解析程度，CFD -DEM又可以分為解析的(resolved)CFD -DEM[13-15]和未解析的(unresolved)CFD -DEM[16-18]，其中解析的CFD -DEM要求流體網(wǎng)格尺度是顆粒尺度的1/10以解析顆粒邊界[19]，其計算結(jié)果非常精確，但是計算量對于工程尺度的模擬而言過于巨大，而未解析的CFD -DEM框架下流體網(wǎng)格尺度可以是顆粒尺度的3倍以上[20,21]。因此，得益于計算效率的顯著優(yōu)勢，未解析的CFD -DEM廣泛應用于工業(yè)級噴動床的數(shù)值模擬當中，也是本文采用的研究方法。由于未解析的CFD -DEM流體網(wǎng)格分辨率低，氣固兩相間最重要的作用力，即曳力的計算需要引入額外的曳力模型，曳力模型的精度很大程度上決定了未解析CFD -DEM模擬的準確性，因此研究不同曳力模型對模擬結(jié)果的影響具有重要意義。目前，關于曳力模型的比較研究大多是在傳統(tǒng)流化床中進行的[22-24]，有些在噴動床中進行的其計算域的離散使用的也是結(jié)構(gòu)化網(wǎng)格[11]，在非結(jié)構(gòu)化網(wǎng)格噴動床中進行的比較研究還未見報道。

綜上，本文基于未解析CFD -DEM和非結(jié)構(gòu)化網(wǎng)格，使用7個曳力模型分別對錐底噴動床內(nèi)氣固兩相運動進行了數(shù)值模擬并比較不同曳力模型對模擬結(jié)果的影響，以期揭示不同曳力模型的作用效果，為曳力模型的選擇提供依據(jù)。

2 CFD -DEM模型

2.1 顆粒運動方程

DEM通過求解牛頓第二定律描述顆粒運動,

(1，2)

式中m和I分別為顆粒質(zhì)量和慣性張量，v和ω分別為線速度和角速度，F(xiàn)c為法向接觸力和切向接觸力的矢量和，F(xiàn)d為曳力，g為重力加速度，Tc為由切向接觸力引起的接觸力矩。顆粒運動方程詳見文獻[25-28]，本文不再贅述。

2.2 氣相控制方程

在未解析CFD -DEM框架下，氣相的運動由局部平均化Navier-Stokes方程描述，

(3)

(4)

(5)

式中Vp，i為顆粒體積，n為CFD單元的顆?？倲?shù)，ΔV為CFD單元的體積。

2.3 曳力模型

在多顆粒系統(tǒng)中，作用在單個顆粒上的曳力可以表示為

(6)

式中dp為顆粒直徑，αs為固含率，β為相間動量交換系數(shù)。

現(xiàn)有的曳力模型可以分為兩類，一類是通過床層壓降測量和顆粒沉降測試等實驗方法獲得的傳統(tǒng)曳力模型，如Wen-Yu模型[29]、Di Felice模型[30]、Gibilaro模型[31]、Gidaspow模型[32]、Huilin-Gidaspow模型[33](簡稱Huilin模型)和Syamlal-O’Brien(簡稱Syamlal模型)模型[34]；另一類是借助如格子玻爾茲曼方法(Lattice -Boltzmann method)和直接數(shù)值模擬(direct numerical simulation)等先進的數(shù)值計算方法得到的曳力模型，如BVK模型[35]。本文考慮了以上7個常用的曳力模型。

Wen-Yu模型[29]是最早出現(xiàn)的幾個曳力模型之一，其基于單顆粒曳力函數(shù)，引入顆粒體積分數(shù)修正因子對曳力模型進行修正，并且考慮了周圍粒子的影響。Wen-Yu模型[29]適用于稀相系統(tǒng)，其可表示為

(7)

(8)

(9)

Gidaspow模型[32]結(jié)合了Ergun方程[36]和Wen-Yu模型[29]，固含率較低時使用Wen-Yu模型[29]，固含率較高時使用Ergun方程[36]，其相間交換系數(shù)的計算如下,

(10)

式中CD的計算同Wen-Yu模型。

Huilin模型[33]也是Ergun方程[36]和Wen-Yu模型[29]的組合，利用一個平滑過渡函數(shù)Ψ解決了Gidaspow模型[32]在計算空隙率跨越 0.8 的稠密氣固兩相流時曳力在數(shù)值上不光滑的問題，其可表示為

(11)

(12)

式中CD的計算同Wen-Yu模型。

Di Felice[30]通過對顆粒堆積體內(nèi)滲流實驗數(shù)據(jù)的整理和擬合推導出Di Felice曳力模型[30]，該模型認為顆粒堆積體內(nèi)流體對顆粒的作用力可以表達為單顆粒在流體內(nèi)平穩(wěn)沉降時所受的曳力與孔隙函數(shù)的乘積，

(13)

γ=3.7-0.65exp [-(1.5-lg Rep,α)2/2]

(14)

(15)

式中 Rep,α的計算同式(9)。

Gibilaro模型[31]也適用于稀相系統(tǒng)的曳力模型，可表示為

(16)

式中 Rep,α的計算同式(9)。

Syamlal模型[34]是基于測量顆粒在流化床或沉降床中的終端速度得到的。

(17)

(18)

vr=0.5(A-0.06Rep)+0.5·[(0.06Rep)2+

0.12Rep(2B-A)+A2]0.5

(19)

(20,21)

(22)

BVK模型[35]是基于大規(guī)模數(shù)值計算推導得到的，可表示為

(23)

(24)

式中 Rep,α的計算同式(9)。

2.4 流體速度插值

根據(jù)2.3節(jié)的介紹可知求解曳力時需要獲得顆粒位置處流體的速度，傳統(tǒng)做法是直接使用顆粒所在網(wǎng)格儲存的平均速度作為顆粒位置處的流體速度，但是在未解析CFD -DEM框架下網(wǎng)格尺度要大于顆粒尺度的3倍以上，這導致網(wǎng)格尺度的平均速度無法準確反應顆粒局部的流體運動情況，因此將網(wǎng)格尺度的平均速度插值到顆粒位置處是非常重要的。本文采用基于梯度的插值方法[37,38]獲得顆粒位置處的流體速度，其中網(wǎng)格的速度梯度使用最小二乘法[39]獲得，具體插值表達式為

up=uc+u·rc p，rc p=xp-xc

(26，27)

式中up為顆粒位置處的流體速度，uc為網(wǎng)格儲存的平均流體速度，rc p為網(wǎng)格質(zhì)心指向顆粒位置的向量，xp為顆粒位置，xc為網(wǎng)格質(zhì)心位置。

3 模擬條件及求解

本文的研究對象為錐底噴動床，如圖1(a)所示。床體高度為890 mm，寬度為250 mm，厚度為10 mm，呈現(xiàn)出典型的準二維結(jié)構(gòu)。錐底角度為60°，噴口寬度為20 mm，氣室高度為100 mm，氣體從氣室底部進入，從床體頂部流出。如圖1(b)所示，整個計算域用非結(jié)構(gòu)化網(wǎng)格離散，其中最小網(wǎng)格的尺寸約為顆粒粒徑的3倍，最大網(wǎng)格的尺寸不超過顆粒5倍以保證計算精度和穩(wěn)定性。模擬中使用的顆粒為相同粒徑的球形玻璃顆粒，氣體為空氣，具體的物性信息及模擬中使用的參數(shù)列入表1。本文分別使用7個曳力模型進行數(shù)值模擬，每個算例除了曳力模型不同外其余設置均相同，入口邊界條件設置為速度入口，氣速為恒定的10 m/s，出口邊界條件設置為壓力出口，壁面設置為無滑移壁面。每個算例開始時顆粒均堆積在床層底部，然后在入口氣體的作用下形成噴動現(xiàn)象。

圖1 設備結(jié)構(gòu)尺寸及計算網(wǎng)格

表1 模擬中使用的參數(shù)

對于連續(xù)相，采用SIMPLE算法進行求解，時間項、對流項和擴散項分別使用Crank-Nicolson格式、QUICK格式和中心差分格式離散。對于顆粒相，采用自主開發(fā)的DEM代碼進行求解，其中顆粒運動的計算方法為顯示時間積分法。每個算例均計算15 s，其中前5 s的數(shù)據(jù)丟棄，只用后10 s的數(shù)據(jù)進行統(tǒng)計平均以排除初始效應。更多關于 CFD -DEM 耦合的細節(jié)可以參見文獻[25,40]，本文不再贅述。

4 結(jié)果與討論

4.1 床層壓降

在模擬中監(jiān)控氣室入口和床體頂部的壓力值，二值之差即為床層壓降。床層壓降及其相關分析可以反映床內(nèi)氣固兩相運動的劇烈程度，是衡量模型準確性的重要指標，因此本文也使用床層壓降來分析不同曳力模型對模擬結(jié)果產(chǎn)生的影響。

圖2為時均床層壓降及其標準差，可以看出，7種模型預測的時均壓降處于同一水平，均在1480 Pa左右，但相應的標準差卻呈現(xiàn)出明顯差異，其中Wen-Yu模型和Gibilaro模型預測的標準差較大，分別為260 Pa和272 Pa，說明在這兩種模型作用下兩相系統(tǒng)運動最劇烈。Di Felice模型預測標準差小于Wen-Yu模型和Gibilaro模型，約為195 Pa，Syamlal模型、Huilin模型以及Gidaspow模型預測的標準差相似，均為155 Pa左右，BVK模型預測的標準差最小，為132 Pa。更進一步地，可以通過離散傅里葉變換獲得床層壓降的頻域信息，如圖3所示。可以看出，除BVK模型外，其余模型預測的壓力信息都有明顯的主頻，大約為3 Hz，符合實驗室規(guī)模噴動床的主頻范圍。Wen-Yu模型和Gibilaro模型預測的主頻幅值最大，均超過了100 Pa，其次是Di Felice模型和Syamlal模型，介于50 Pa～100 Pa，最后是Gidaspow模型、Huilin模型和BVK模型，三者預測的主頻幅值均低于50 Pa。

圖2 時均床層壓降及標準差

通過對床層壓降的分析可知，Wen-Yu模型和Gibilaro模型預測的床層運動最劇烈，BVK模型預測的床層運動最平緩，其余四個模型的預測效果介于中間。此外，可以發(fā)現(xiàn)Gidaspow模型和Huilin模型的預測結(jié)果幾乎一致，這是因為本文研究的氣固兩相體系屬于密相體系，光滑過渡函數(shù)沒有產(chǎn)生效果。

4.2 噴動高度

對不同高度處顆粒軸向速度的時均處理可以得到每個曳力模型預測的噴動高度，如圖4所示?？梢钥闯?，Wen-Yu模型預測的噴動高度最大，達到了279 mm，其次是Di Felice模型和Gibilaro模型，分別為267 mm和261 mm,Gidaspow模型和Huilin模型預測的噴動高度相同，均為255 mm，再是Syamlal模型，預測噴動高度為249 mm，最后是BVK模型，預測噴動高度為最小的243 mm。

圖4 不同曳力模型預測的噴動高度

為了更直觀地展示不同曳力模型對噴動床噴動高度的影響，圖5給出了2 s后不同曳力模型模擬的瞬態(tài)快照?？梢钥闯鯳en-Yu模型和Gibilaro模型的算例中顆粒運動最劇烈，噴動高度最高，BVK模型計算得到的噴動高度最低，其余四種模型的計算結(jié)果介于中間，這與噴動高度定量統(tǒng)計的結(jié)果一致。此外，由于在計算顆粒位置處的流體速度時使用了梯度插值的方法，圖5展示的瞬態(tài)快照沒有呈現(xiàn)出明顯的網(wǎng)格效應，氣相和顆粒相的分界面比較光滑。

圖5 2 s后的顆粒運動快照

可以看出，對噴動高度的分析得到了與床層壓降分析幾乎一致的結(jié)論，Wen-Yu模型和Gibilaro模型由于稀相適用性導致其在本文研究的密相噴動床中預測得到最劇烈的氣固兩相運動，因此統(tǒng)計的壓降波動和噴動高度也最大，而通過LBM得到的BVK模型計算得到的氣固兩相運動最平緩，相應的壓降波動和噴動高度也是所有模型中最小的。

4.3 顆粒軸向速度

圖6展示了各個曳力模型的時均顆粒軸向速度的徑向分布(距離氣體入口190 mm處)。可以看出各曳力模型的預測結(jié)果均呈現(xiàn)出中心正，邊壁負的分布特征，說明顆粒在噴動床中心區(qū)域向上運動，沿徑向顆粒速度逐漸減小，直至近壁面區(qū)域顆粒向下運動，顆粒在噴動床錐部形成了典型的環(huán)-核流動結(jié)構(gòu)，這種結(jié)構(gòu)特征已得到廣泛報道[1,41-43]。對比各個曳力模型的模擬結(jié)果發(fā)現(xiàn)，Wen-Yu模型和Gibilaro模型預測的軸向速度最大，其次是Di Felice模型,Syamlal模型、Gidaspow模型和Huilin模型的預測結(jié)果非常接近，BVK模型預測的軸向速度最小。

圖6 190 mm高度處的顆粒時均軸向速度

5 結(jié) 論

本文基于未解析的CFD -DEM和非結(jié)構(gòu)化網(wǎng)格，在錐底噴動床中比較了Wen-Yu模型、Gidaspow模型、Huilin模型、Di Felice模型、Gibilaro模型、Syamlal模型和BVK模型七個曳力模型對模擬結(jié)果產(chǎn)生的影響。綜合床層壓降、噴動高度和顆粒速度特性三個方面，Wen-Yu模型和Gibilaro模型預測的氣固兩相運動最劇烈，其次是Di Felice模型、Syamlal模型、Gidaspow模型和Huilin模型，BVK模型預測的氣固兩相運動最平緩。此外，由于本文研究的氣固兩相體系屬于密相體系，Huilin模型的光滑過渡函數(shù)沒有產(chǎn)生效果，所以Gidaspow模型和Huilin模型在各個方面的預測結(jié)果基本一致。此外，如果想判定特定工況下最優(yōu)的曳力模型，還需要進一步的實驗加以對照。