不同液固動(dòng)量交換系數(shù)模型對(duì)顆粒流化行為數(shù)值模擬的適應(yīng)性

孫銘陽(yáng)，于傳兵，呂東，韋魯濱

(1.中國(guó)恩菲工程技術(shù)有限公司，北京，100038；2.中國(guó)礦業(yè)大學(xué)(北京) 化學(xué)與環(huán)境工程學(xué)院，北京，100083)

由于其優(yōu)異的混合、傳質(zhì)和傳熱特性，液固流態(tài)化技術(shù)已被廣泛用于冶金、選礦、化工和環(huán)保行業(yè)[1-4]。國(guó)內(nèi)外學(xué)者針對(duì)液固流化床進(jìn)行了大量研究，以揭示流化床內(nèi)顆粒運(yùn)動(dòng)特點(diǎn)和床層膨脹等行為。胡新輝等[5]對(duì)液固流化床內(nèi)顆粒循環(huán)運(yùn)動(dòng)機(jī)理進(jìn)行了研究，提出了液固流化床內(nèi)流體和顆粒徑向速度分布的數(shù)學(xué)模型，并認(rèn)為顆粒循環(huán)運(yùn)動(dòng)由流體徑向速度分布不均導(dǎo)致；DOROODCHI等[6]借助CFD研究了傾斜板對(duì)顆粒在液固分選流化床內(nèi)流化行為的影響，并比較了CFD模型和理論模型的優(yōu)劣；張鍇等[7-8]改進(jìn)了雙流體模型，從而較全面地考慮了離散顆粒對(duì)液固兩相動(dòng)量的影響，并用該改進(jìn)的模型對(duì)入口速度突變時(shí)液固流化床內(nèi)固體顆粒時(shí)空分布特性進(jìn)行了詳細(xì)研究；沈志恒等[9]借助CFD對(duì)倒置液固流化床內(nèi)顆粒速度、濃度和床層孔隙率等分布特性進(jìn)行了模擬，發(fā)現(xiàn)倒置業(yè)液固流化床內(nèi)存在中間向下兩側(cè)向上的顆粒運(yùn)動(dòng)規(guī)律；張衛(wèi)義等[10]對(duì)液固流化床內(nèi)顆粒最小流化速度進(jìn)行了研究，提出了一個(gè)Ergun型最小流化速度模型，并與其他模型進(jìn)行了對(duì)比，在其所列舉的試驗(yàn)數(shù)據(jù)下，該模型的預(yù)測(cè)結(jié)果比Wen-Yu模型的優(yōu)。 液固動(dòng)量交換系數(shù)是影響顆粒流化過程數(shù)值模擬準(zhǔn)確性的重要因素。目前針對(duì)液固動(dòng)量交換系數(shù)對(duì)液固流化床數(shù)值模擬準(zhǔn)確性影響方面的研究較少。本文作者對(duì)不同粒度或密度的9種顆粒流化過程分別進(jìn)行試驗(yàn)研究，并結(jié)合試驗(yàn)結(jié)果對(duì)比分析不同液固動(dòng)量交換系數(shù)對(duì)顆粒流化行為數(shù)值模擬的準(zhǔn)確性，本文試驗(yàn)方法和數(shù)值模擬結(jié)論可為相關(guān)液固流化床數(shù)值模擬研究提供參考。

1 顆粒流化試驗(yàn)與結(jié)果

1.1 流化顆粒與流化試驗(yàn)裝置

流化顆粒分別選擇不同粒度、密度的石英砂和煤顆粒，共包括如表1所示的9種顆粒。

搭建了如圖1所示的顆粒流化試驗(yàn)系統(tǒng)。上升水流由平行管流體分布器射出后，沿柱體向上流動(dòng)并在徑向上逐漸均勻分布，經(jīng)第1層孔徑0.045 mm篩網(wǎng)后，流體軸向速度基本實(shí)現(xiàn)均勻分布。流化顆粒置于第2層孔徑0.045 mm篩網(wǎng)上，篩網(wǎng)孔徑遠(yuǎn)小于流化顆粒直徑，可以保證底層顆粒受到均勻向上的流體曳力，且不會(huì)沿近壁面區(qū)漏下。流化床溢流給到過濾箱，過濾箱內(nèi)亦布置1塊孔徑0.045 mm篩網(wǎng)，溢流經(jīng)篩網(wǎng)過濾后返回至清水池，從而實(shí)現(xiàn)帶出顆粒的回收和流化水的循環(huán)利用。

表1 流化顆粒物性Table 1 Properties of fluidized particles

圖1 液固流化試驗(yàn)系統(tǒng)Fig.1 Experimental system for liquid-solid fluidization

1.2 流化試驗(yàn)結(jié)果

令H為上升水流速度為uF時(shí)顆粒床層高度，則此時(shí)液固流化床層內(nèi)顆粒平均體積分?jǐn)?shù)可表示為

床層孔隙率為

床層有效密度則可按下式計(jì)算：

式中：ρl和ρs分別為液相和顆粒相密度。

圖2所示為不同顆粒所形成液固流化床層有效密度隨上升水流變化情況。從圖2可以看出：顆粒粒度越小，密度越小，所形成的液固流化床層密度對(duì)上升水流越敏感。當(dāng)2種顆粒密度相差較大時(shí)，在相同流化速度下，不同顆粒形成的床層密度流化曲線基本不能相交。通過回歸分析發(fā)現(xiàn)：床層密度ρb與上升水流速度u0呈現(xiàn)形式的自然指數(shù)函數(shù)關(guān)系，同時(shí)也可被形式的冪函數(shù)較好擬合。

圖2 液固流化床層密度與上升水流關(guān)系Fig.2 Bed effective density at various fluidized velocities

2 顆粒流化過程數(shù)值模擬方法

2.1 液固兩相流動(dòng)控制方程

采用 Euler-Euler法模擬液固流化床層內(nèi)固液兩相流動(dòng)。

液、固兩相的質(zhì)量守恒均滿足下式：

式中：αp，up和ρp分別為p相體積分?jǐn)?shù)、速度和密度；p和k分別表示固相或液相，當(dāng)p為固相時(shí)則k為液相，反之亦然；mpk和mkp分別為由p相到k相的轉(zhuǎn)移質(zhì)量和由k相到p相的轉(zhuǎn)移質(zhì)量；Sp為p相的質(zhì)量源相項(xiàng)。對(duì)于惰性顆粒流化過程，mpk，mkp和Sp可忽略，且液固流化床內(nèi)流動(dòng)可看作不可壓縮，因此，式(4)可進(jìn)一步簡(jiǎn)化為

液相與各顆粒相存在相互作用，液相動(dòng)量方程為

式中：αl和ul分別為液相體積分?jǐn)?shù)和速度；us為顆粒相的速度；τl為液相應(yīng)力張量；g為重力加速度；Ksl(ul-us)為液相與顆粒相的交換動(dòng)量，Ksl為液相與固相間的動(dòng)量交換系數(shù)；最后一項(xiàng)為液相因施加給離散固體顆粒相升力、壓力梯度力和湍流分散力等而受到的相應(yīng)反作用。

顆粒相s的動(dòng)量方程為

式中：αs為顆粒相體積分?jǐn)?shù)；τs為顆粒相應(yīng)力張量；Flift,s為顆粒相所受升力；Fvm,s為顆粒相所受附加質(zhì)量力；Ftd,s為湍流中顆粒相所受湍流分散力；Kls(ul-us)為顆粒相所受流體阻力。

式(7)中τs由應(yīng)力張量可求顆粒相微元發(fā)生形變時(shí)所受法向和切向應(yīng)力。通過顆粒相應(yīng)力張量得到由速度和壓力表示的顆粒相運(yùn)動(dòng)方程，進(jìn)而通過求解該方程得到顆粒相流場(chǎng)。顆粒相應(yīng)力張量具有以下形式：

式中：μs為顆粒相剪切黏度；I為單位張量；Ps為顆粒相壓力；λs為顆粒相體積黏度；這些變量為顆粒脈動(dòng)能或法向恢復(fù)系數(shù)的函數(shù)。

類比氣體分子動(dòng)理論，可以得到顆粒溫度的概念。顆粒溫度實(shí)際表征了顆粒相中固體顆粒隨機(jī)脈動(dòng)速度的強(qiáng)度，因此又稱顆粒脈動(dòng)能，其表達(dá)式如下：

式中：us,i為顆粒i方向的脈動(dòng)速度。

根據(jù)顆粒動(dòng)理學(xué)理論，顆粒脈動(dòng)能的輸運(yùn)方程為

2.2 計(jì)算域離散和求解方法

數(shù)值模擬中液固流化床幾何結(jié)構(gòu)如圖3所示，將液固流化床層計(jì)算域劃分為結(jié)構(gòu)網(wǎng)格。速度入口給入均勻穩(wěn)定的上升水流，上升水流速度根據(jù)試驗(yàn)值進(jìn)行設(shè)定。數(shù)值模擬選擇二維雙精度瞬態(tài)求解器。雖然表觀上升水流速度在0.1 m/s以下，但流化床內(nèi)的流動(dòng)是工程上典型湍流流動(dòng)，選擇RNGk-ε湍流模型來(lái)預(yù)測(cè)流化床內(nèi)液固兩相運(yùn)動(dòng)[11-12]。采用有限體積法(FVM)對(duì)控制方程進(jìn)行離散；選擇 QUICK插值方法，將控制單元節(jié)點(diǎn)處各物理量值表示成相鄰控制單元面上值；選擇 Phase coupled SIMPLE算法來(lái)進(jìn)行壓力與速度的耦合求解。當(dāng)液固兩相體積分?jǐn)?shù)殘差小于10-4，其余物理量的殘差小于10-3時(shí)，認(rèn)為該時(shí)間步內(nèi)的迭代達(dá)到收斂[13]。

圖3 數(shù)值模擬中液固流化床幾何模型Fig.3 Geometry and mesh of fluidized bed

3 不同液固動(dòng)量交換系數(shù)模型模擬結(jié)果

3.1 液固動(dòng)量交換系數(shù)模型

式(6)和式(7)中Kls=Ksl，為顆粒相與液相動(dòng)量交換系數(shù)。

式中：τs為顆粒弛豫時(shí)間；f為與顆粒相體積分?jǐn)?shù)和阻力系數(shù)等相關(guān)的函數(shù)。

由式(11)可知：動(dòng)量交換系數(shù)是曳力系數(shù)、顆粒相體積分?jǐn)?shù)和液固相對(duì)速度等物理量的函數(shù)。目前，液固兩相動(dòng)量交換系數(shù)的數(shù)學(xué)模型主要有Syamlal-O’Brien[14]、Wen-Yu[15]、Gidaspow[16]、Gibilaro[17]和Huilin-Gidaspow[18]。

Syamlal-O’Brien模型如下：

阻力系數(shù)CD采用DALLA VALLE表達(dá)式[19]求得：

Res為顆粒相雷諾數(shù)，有

式中：ds為s相顆粒直徑。

式(12)中ur,s為s相顆粒沉降末速修正，可按下式求得：

Wen-Yu模型數(shù)學(xué)形式如下：

Wen-Yu模型中曳力系數(shù)表達(dá)式為

Gidaspow模型為Wen-Yu模型和Ergun方程的組合，當(dāng)液相體積分?jǐn)?shù)α1＞0.8時(shí)，表達(dá)式為Wen-Yu模型；當(dāng)液相體積分?jǐn)?shù)α1≤0.8時(shí)，表達(dá)式為Ergun方程，具體形式如下：

Gibilaro模型具有以下形式：

Huilin-Gidaspow模型也是Wen-Yu模型和Ergun方程的組合，當(dāng)床層空隙率的變化跨越臨界點(diǎn)(0.8)時(shí)，利用一個(gè)平滑過渡函數(shù)來(lái)實(shí)現(xiàn)Wen-Yu模型和Ergun方程間的轉(zhuǎn)換，其形式如下：

式中：Ψ為平滑過渡函數(shù)，Ψ=0.5+arctant(262.5×(αs-0.2))/π。

3.2 不同液固動(dòng)量交換系數(shù)模擬床層高度對(duì)比

利用以上 5個(gè)液固動(dòng)量交換系數(shù)模型模擬0.25～0.35 mm石英砂在12.3 mm/s上升水流速度下流化過程。床層運(yùn)動(dòng)穩(wěn)定后，對(duì)比5個(gè)模型模擬所得床層內(nèi)固體顆粒體積分?jǐn)?shù)分布如圖4所示。從圖4可知：預(yù)測(cè)床層體積分?jǐn)?shù)最大的為Gibilaro模型，最小的為Syamlal-O’Brien模型；Gidaspow 和 Huilin-Gidaspow模型預(yù)測(cè)的床層內(nèi)顆粒體積分?jǐn)?shù)分布完全一致；Wen-Yu模型與Gidaspow 和Huilin-Gidaspow模型相比，預(yù)測(cè)的顆粒體積分?jǐn)?shù)略高。

圖4 12.3 mm/s流化速度下各模型預(yù)測(cè)0.25～0.35 mm石英砂顆粒體積分?jǐn)?shù)分布Fig.4 Distribution of 0.25-0.35 mm quartz particles at 12.3 mm/s fluidization rate

床層高度模擬值確定方法如圖5所示，用不同液固動(dòng)量交換系數(shù)模擬所得不同上升水流速度下床層高度與試驗(yàn)值對(duì)比如圖6所示。從圖6可知：Wen-Yu，Gidaspow 和Huilin-Gidaspow對(duì)床層高度的預(yù)測(cè)值與試驗(yàn)值最吻合，且Gidaspow 和Huilin-Gidaspow床層高度預(yù)測(cè)值完全相同；Wen-Yu，Gidaspow，Syamlal-O’Brien，Gibilaro和Huilin-Gidaspow這5個(gè)模型所得床層高度預(yù)測(cè)值與試驗(yàn)值間的均方根誤差分別為5.67，5.01，47.14，24.42和5.01，表明Gidaspow和Huilin-Gidaspow模型準(zhǔn)確度略高于Wen-Yu模型的準(zhǔn)確度。

圖5 預(yù)測(cè)的床層高度h確定方法示意圖Fig.5 Method of determining height of predicted fluidized bed

圖6 不同液固動(dòng)量交換系數(shù)模型模擬結(jié)果對(duì)比Fig.6 Predicted bed height by using different fluid-solid exchange coefficients

圖7所示為圖5中對(duì)應(yīng)的石英砂顆粒相速度矢量圖。從圖7可知：Gibilaro和Gidaspow模型得到的石英砂相只有在液固界面附近存在較大的速度變化，床層的其余位置處石英砂顆粒的速度很小，可認(rèn)為準(zhǔn)靜止地懸浮于上升水流中；Wen-Yu，Syamlal-O’Brien和 Huilin-Gidaspow模型得到的石英砂相存在由中心向上再流向壁面、沿壁面下降的循環(huán)流動(dòng)，這與試驗(yàn)中觀察到的床層流動(dòng)狀態(tài)吻合。因此，從床層高度預(yù)測(cè)值和床層內(nèi)液固兩相流態(tài) 2方面來(lái)講，Huilin-Gidaspow模型比其他4個(gè)模型更合適。

圖7 石英砂顆粒相速度矢量分布Fig.7 Velocity vector distribution of quartz particles phase

圖8 0.25～0.35 mm 石英砂在12.3 mm/s流化速度下流化過程Fig.8 Fluidization process of 0.25-0.35 mm quartz particles at 12.3 mm/s fluidization rate

圖8所示為0.25～0.35 mm石英砂在12.3 mm/s上升水流速度下流化過程。從圖8可知：隨著流動(dòng)時(shí)間的增長(zhǎng)，床層膨脹度逐漸增大；40 s時(shí)床層高度基本達(dá)到穩(wěn)定，但此時(shí)床層內(nèi)液固兩相流動(dòng)還未達(dá)到平衡，由石英砂相的流動(dòng)速度可以看出，0～100 s之間顆粒相運(yùn)動(dòng)比較混亂，沒有明顯規(guī)律；到達(dá)120 s后，床層內(nèi)液固兩相流動(dòng)基本實(shí)現(xiàn)穩(wěn)定，顆粒相表現(xiàn)出由中心向上，到達(dá)床層頂端后，向兩側(cè)壁面流動(dòng)，最后沿壁面低速區(qū)下降的循環(huán)流動(dòng)趨勢(shì)，與試驗(yàn)中觀察到的現(xiàn)象吻合。

圖9 4種顆粒流化床層顆粒體積分?jǐn)?shù)模擬值與試驗(yàn)值Fig.9 Comparison between simulated and experimental volume fraction of particles phase

進(jìn)一步選擇Hui-Gidaspow模型模擬0.15～0.20和0.50～0.63 mm石英砂顆粒以及 0.25～0.35和 0.90～1.00 mm精煤顆粒在不同上升水流速度下流化過程。圖9 所示為4種顆粒在不同上升水流速度下床層顆粒體積分?jǐn)?shù)試驗(yàn)值和預(yù)測(cè)值。從圖9可知：床層顆粒體積分?jǐn)?shù)預(yù)測(cè)相對(duì)誤差基本都在 5%以內(nèi)，說明在用Euler-Euler法模擬顆粒流化過程中，選擇 Hui-Gidaspow模型得到的模擬結(jié)果具有較高的準(zhǔn)確度。

4 結(jié)論

1) 5種液固動(dòng)量交換系數(shù)模型對(duì)粒度為 0.25～0.35 mm的石英砂顆粒流化床層高度預(yù)測(cè)準(zhǔn)確度由高到低依此為：Gidaspow模型和Huilin-Gidaspow模型，Wen-Yu模型，Gibilaro模型，Syamlal-O’Brien模型，其中Huilin- Gidaspow和Gidaspow模型均方根誤差都為5.01，Syamlal-O’Brien模型均方根誤差為47.14。

2) 利用Huilin-Gidaspow模型模擬得到的液固流化床層內(nèi)顆粒相存在由中心向上，到達(dá)床層頂端后，再向兩側(cè)壁面流動(dòng)，最后沿壁面低速區(qū)下降的循環(huán)流動(dòng)過程，與試驗(yàn)中觀察到的現(xiàn)象相符，而 Gidaspow模型模擬結(jié)果中顆粒相不存在循環(huán)流動(dòng)。從床層高度預(yù)測(cè)值和顆粒相運(yùn)動(dòng)狀態(tài) 2方面考慮，5個(gè)模型中Huilin-Gidaspow模型最適合液固流化床內(nèi)顆粒流化行為數(shù)值模擬。

3) 選擇 Huilin-Gidaspow模型對(duì) 0.15～0.20和0.50～0.63 mm石英砂顆粒以及 0.25～0.35和 0.90～1.00 mm精煤顆粒在不同上升水流速度下流化過程進(jìn)行了模擬，床層內(nèi)顆粒相體積分?jǐn)?shù)預(yù)測(cè)值相對(duì)誤差基本都處于5%以內(nèi)，說明利用Euler-Euler法，同時(shí)選擇 Hui-Gidaspow液固動(dòng)量交換系數(shù)模型模擬顆粒流化過程具有較高準(zhǔn)確度。