周連勇， 趙永志

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

1 引 言

噴動床作為一種特殊的流化床，具有結(jié)構(gòu)簡單和接觸效率高等優(yōu)點[1,2]，已經(jīng)在農(nóng)作物干燥、生物制藥以及能源化工等領(lǐng)域得到廣泛應(yīng)用[3-5]。噴動床的工作效率主要由其內(nèi)部復(fù)雜的氣固兩相運動決定，而常規(guī)的實驗方法難以深入探究氣固兩相運動的機理，無法系統(tǒng)地掌握流動結(jié)構(gòu)和動力學(xué)參數(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)化設(shè)計中逐漸成為主流。根據(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īng)用于工業(yè)級噴動床的數(shù)值模擬當中，也是本文采用的研究方法。由于未解析的CFD -DEM流體網(wǎng)格分辨率低，氣固兩相間最重要的作用力，即曳力的計算需要引入額外的曳力模型，曳力模型的精度很大程度上決定了未解析CFD -DEM模擬的準確性，因此研究不同曳力模型對模擬結(jié)果的影響具有重要意義。目前，關(guān)于曳力模型的比較研究大多是在傳統(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ù)的整理和擬合推導(dǎo)出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ù)值計算推導(dǎo)得到的，可表示為

(23)

(24)

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

2.4 流體速度插值

根據(jù)2.3節(jié)的介紹可知求解曳力時需要獲得顆粒位置處流體的速度，傳統(tǒng)做法是直接使用顆粒所在網(wǎng)格儲存的平均速度作為顆粒位置處的流體速度，但是在未解析CFD -DEM框架下網(wǎng)格尺度要大于顆粒尺度的3倍以上，這導(dǎo)致網(wǎng)格尺度的平均速度無法準確反應(yī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ù)值模擬，每個算例除了曳力模型不同外其余設(shè)置均相同，入口邊界條件設(shè)置為速度入口，氣速為恒定的10 m/s，出口邊界條件設(shè)置為壓力出口，壁面設(shè)置為無滑移壁面。每個算例開始時顆粒均堆積在床層底部，然后在入口氣體的作用下形成噴動現(xiàn)象。

圖1 設(shè)備結(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)計平均以排除初始效應(yīng)。更多關(guān)于 CFD -DEM 耦合的細節(jié)可以參見文獻[25,40]，本文不再贅述。

4 結(jié)果與討論

4.1 床層壓降

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

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

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

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

4.2 噴動高度

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

圖4 不同曳力模型預(yù)測的噴動高度

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

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

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

4.3 顆粒軸向速度

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

圖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模型預(yù)測的氣固兩相運動最劇烈，其次是Di Felice模型、Syamlal模型、Gidaspow模型和Huilin模型，BVK模型預(yù)測的氣固兩相運動最平緩。此外，由于本文研究的氣固兩相體系屬于密相體系，Huilin模型的光滑過渡函數(shù)沒有產(chǎn)生效果，所以Gidaspow模型和Huilin模型在各個方面的預(yù)測結(jié)果基本一致。此外，如果想判定特定工況下最優(yōu)的曳力模型，還需要進一步的實驗加以對照。