朱衛(wèi)兵，朱潤(rùn)孺，邢力超，孫巧群

(哈爾濱工程大學(xué) 航天與建筑工程學(xué)院，黑龍江 哈爾濱 150001)

噴動(dòng)床最早由加拿大的Mathur和Gishler發(fā)明，最初用于谷物干燥，如今噴動(dòng)床已經(jīng)被廣泛的應(yīng)用到干燥、造粒、粉碎、煤燃燒和氣化、鐵礦石還原、油頁巖熱解、焦炭活化、石油熱裂化等領(lǐng)域［1］.在典型的噴動(dòng)床型中，矩形噴動(dòng)床與柱錐型噴動(dòng)床相比，不但具有較低的最小噴動(dòng)速度和最大噴動(dòng)壓降的優(yōu)點(diǎn)，還具有氣固接觸效率好、傳熱效率高的優(yōu)點(diǎn)，且?guī)缀谓Y(jié)構(gòu)簡(jiǎn)單，易于制造.

顆粒在噴動(dòng)床內(nèi)運(yùn)動(dòng)，受到流體曳力、自身重力和顆粒間接觸力等作用力，其中流體曳力是對(duì)噴動(dòng)床進(jìn)行數(shù)值研究時(shí)，唯一的相間相互作用力，因此曳力模型的選擇是數(shù)值結(jié)果精確與否的關(guān)鍵.Du Wei［2］等人采用雙流體模型(two fluid model，TFM)對(duì)柱錐型噴動(dòng)床的曳力模型選擇進(jìn)行了研究，Li Jie＆Kuipers J A M［3］對(duì)Hill等模型在流化床離散元法(discrete element method，DEM)模擬中的應(yīng)用做了比較，而對(duì)于矩形噴動(dòng)床相同問題的研究還未見報(bào)道，由此可見對(duì)矩形噴動(dòng)床內(nèi)顆粒運(yùn)動(dòng)現(xiàn)象的研究遠(yuǎn)不如柱錐床深入.噴動(dòng)床DEM模擬中，Tsuji［4］模型和Ergun(ε＞0.8)+Wen－Yu(ε＜0.8)聯(lián)合模型［5］是應(yīng)用最廣泛的曳力模型，而其他曳力模型，如Arastoopour［6］，Di Felice［7］和Syamlal＆O＇Brien［8］等模型還未見有研究報(bào)道采用.噴動(dòng)床中，滾動(dòng)的顆粒會(huì)受到滾動(dòng)摩擦的作用，但在以往的噴動(dòng)床DEM模擬中，滾動(dòng)摩擦很少被考慮，只考慮滑動(dòng)摩擦對(duì)顆粒運(yùn)動(dòng)的影響，這顯然會(huì)使模擬與真實(shí)物理現(xiàn)象間存在差異.

針對(duì)噴動(dòng)床DEM模擬中的曳力模型選擇問題，以DEM數(shù)值模擬為手段，首次將Arastoopour、Di Felice、Syamlal＆O＇Brien曳力模型和Beer＆Johnson［9］滾動(dòng)摩阻力矩表達(dá)式應(yīng)用到矩形噴動(dòng)床DEM(pseudo－three－dimensional)［10］模擬中，對(duì)矩形噴動(dòng)床內(nèi)顆粒運(yùn)動(dòng)進(jìn)行研究，并將模擬結(jié)果與趙香龍［11］等人的實(shí)驗(yàn)結(jié)果比較，具體分析了不同曳力模型下的床內(nèi)流動(dòng)結(jié)構(gòu)和空隙率、顆粒速度、力矩分布，以期為矩形噴動(dòng)床數(shù)值研究中曳力模型的選擇提供有益參考.

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

CFD－DEM模型中，氣相場(chǎng)采用歐拉法，固相場(chǎng)采用拉格朗日法.對(duì)于相間耦合，文中分別選擇了Arastoopour、Di Felice、Syamlal＆O＇Brien和Tsuji等曳力模型.

1.1 氣相模型

對(duì)于氣相場(chǎng)采用標(biāo)準(zhǔn)k－ε模型，考慮氣體的湍流運(yùn)動(dòng).連續(xù)性方程和動(dòng)量方程如下:

式中:ε、u、p、ρg和μ分別為孔隙率、氣相速度、壓力、密度和粘性.fd=β(u－v)，曳力系數(shù)β由式(11)～(14)確定.

湍流動(dòng)能和湍流動(dòng)量耗散方程如下:

式中:k、εt分別為湍流動(dòng)能和湍流動(dòng)量耗散率.μt=cμρgk2/εt.方程中 c1、c2、cμ、σε和 σk的值見表1.

表1 模擬參數(shù)Table 1 Parameters for particle and fluid simulations

1.2 顆粒運(yùn)動(dòng)模型

顆粒的運(yùn)動(dòng)遵循牛頓第二定律，考慮重力、接觸力、流體曳力和壓力梯度力的作用.平動(dòng)和轉(zhuǎn)動(dòng)運(yùn)動(dòng)方程為

式中:

式中:dp、e、g、I、k、mi、Mi、nij、v、sij、δn，ij、δt，ij、μ、μr、η和ωi分別為顆粒直徑、顆粒恢復(fù)系數(shù)、重力加速度、轉(zhuǎn)動(dòng)慣量、彈性系數(shù)、單個(gè)顆粒小球質(zhì)量、滾動(dòng)摩擦力矩、顆粒小球i和j之間的法向單位矢量、顆粒小球速度、顆粒小球i和j之間的切向單位矢量、顆粒小球i和j之間的法向變形、顆粒小球i和j之間的切向變形、滑動(dòng)摩擦系數(shù)、滾動(dòng)摩擦系數(shù)、阻尼系數(shù)和顆粒小球角速度.μ和μr的值見表1.

1.3 曳力模型

Arastoopour模型:

Di Felice模型:

式中:X=3.7－0.65exp［－0.5(1.5－lg Re)2］，

Syamlal＆O＇Brien模型:

式中:vrm=0.5{A－0.06Re+［(0.06Re)2+;式(11)～(13)中，Re= (dp|u－v|ρg)/μ.

Tsuji模型:

1.4 模型及求解

圖1 幾何模型Fig.1 Geometry of a vessel

模擬的對(duì)象為矩形噴動(dòng)床，床體尺寸見圖1，其中床厚15 mm，噴口寬9 mm.模擬所用網(wǎng)格尺寸為3 mm×2 mm×5.2 mm.計(jì)算顆粒直徑2 mm，顆粒密度2 380 kg/m3，顆粒填充高度100 mm，氣相參數(shù)選擇空氣在20℃，1.01×105Pa時(shí)的參數(shù)，表觀氣速Ug=1.58 m/s，固相顆粒計(jì)算時(shí)間步長(zhǎng)為10－6s.采用交錯(cuò)網(wǎng)格技術(shù)，參照SIMPLE算法，對(duì)氣相場(chǎng)和固相場(chǎng)連續(xù)交替耦合求解.對(duì)于模擬時(shí)最小噴動(dòng)氣速的確定，采用逐漸降低表觀氣速的方法.表觀氣速由1.6 m/s降到0.8 m/s后發(fā)現(xiàn)，4種模型的最小噴動(dòng)速度均小于模擬工況.

2 結(jié)果與討論

2.1 曳力系數(shù)分析

圖2為不同曳力模型曳力系數(shù)隨空隙率變化曲線，氣固兩相滑移速度選擇15m/s.由圖可見:1)空隙率大于0.8區(qū)域，不同曳力模型計(jì)算出的曳力系數(shù)值差異較小;2)空隙率在0.66～0.77內(nèi)Tsuji模型曳力系數(shù)最大;3)空隙率大于0.47區(qū)域，Arastoopour、Syamlal＆O＇Brien模型分別和Tsuji模型有2個(gè)交匯點(diǎn)，在交匯點(diǎn)區(qū)間內(nèi)，Arastoopour、Syamlal＆O＇Brien模型曳力系數(shù)均小于Tsuji模型曳力系數(shù); 4)空隙率小于0.8區(qū)域，Di Flice模型曳力系數(shù)隨著空隙率的減小而迅速增加，Arastoopour模型也顯現(xiàn)出類似的趨勢(shì)，但在數(shù)值上較小;5)空隙率小于0.47區(qū)域，Tsuji模型曳力系數(shù)最小.由以上分析可知空隙率變化造成了不同模型曳力系數(shù)的差異，而氣固曳力是噴動(dòng)床DEM模擬相間耦合的唯一作用力，并且噴動(dòng)床內(nèi)空隙率在不同區(qū)域相差較大，所以選擇不同的曳力模型必將造成模擬的結(jié)果不同.

圖2 曳力系數(shù)比較Fig.2 Comparison of drag force coefficients

2.2 流動(dòng)結(jié)構(gòu)

圖3為不同曳力模型在1 s時(shí)顆粒分布的模擬結(jié)果，由圖可以看出在同一時(shí)刻不同模型噴泉高度由高到低排列為:Di Felice、Arastoopour、Tsuji和Syamlal＆O＇Brien模型，Di Felice模型的噴射區(qū)直徑最大，相應(yīng)環(huán)隙區(qū)范圍要小于其他3個(gè)模型，這是由于其在空隙率小于0.55區(qū)域曳力系數(shù)較大所致.

圖4是不同曳力模型下的流動(dòng)結(jié)構(gòu)模擬結(jié)果，由圖可以看出 Arastoopour、Syamlal＆O＇Brien和Tsuji模型的模擬結(jié)果均呈現(xiàn)出一種周期性噴動(dòng)現(xiàn)象，周期大約為0.16 s，并且可以看出顆粒在噴泉區(qū)的運(yùn)動(dòng)呈現(xiàn)出均勻且對(duì)稱性較好的顆粒群運(yùn)動(dòng)狀態(tài)，這兩點(diǎn)與圖5的實(shí)驗(yàn)結(jié)果符合很好，而Di Felice模型的模擬結(jié)果沒有出現(xiàn)這2種現(xiàn)象，所以不再對(duì)其進(jìn)行進(jìn)一步的討論.

圖3 不同曳力模型1 s時(shí)的流動(dòng)結(jié)構(gòu)Fig.3 The instantaneous flow patterns of different drag models at 1.0 s

圖4 流動(dòng)結(jié)構(gòu)模擬結(jié)果Fig.4 Simulated result of flow patterns

圖5 流動(dòng)結(jié)構(gòu)實(shí)驗(yàn)結(jié)果Fig.5 Experimental result of flow patterns

2.3 空隙率分布

圖6給出了不同床層高度下沿徑向的空隙率分布.

圖6 空隙率分布Fig.6 Voidage profiles

由圖可知3個(gè)模型的模擬結(jié)果均有以下現(xiàn)象: 1)隨著床層高度的增加，噴射區(qū)的空隙率逐漸降低，且最小空隙率出現(xiàn)點(diǎn)外移，這與 Lim＆Mathur［12］給出的式(16)預(yù)測(cè)趨勢(shì)相符;2)沿床體徑向，空隙率成下降趨勢(shì)，在噴射區(qū)和環(huán)隙區(qū)交界處呈現(xiàn)出很大的梯度，而在環(huán)隙區(qū)內(nèi)變化較平緩.這2種現(xiàn)象和He［13－14］等人實(shí)驗(yàn)觀察得到的結(jié)論相同.

式中:Ds和np為噴射區(qū)直徑和噴射區(qū)給定床高處顆粒數(shù)目.

2.4 軸線上的顆粒速度

圖7(a)為顆粒軸向速度分布，實(shí)驗(yàn)結(jié)果表明沿軸線顆粒軸向速度由3個(gè)區(qū)組成:1)加速區(qū)，2)平緩區(qū)，3)減速區(qū).這種現(xiàn)象有別于以往的實(shí)驗(yàn)結(jié)果［13－15］.采用標(biāo)準(zhǔn)k－ε模型和Beer＆Johnson表達(dá)式加入湍流和滾動(dòng)摩擦2種因素的影響，區(qū)別于文獻(xiàn)［11］采用Jones＆Launder低雷諾數(shù)k－ε模型和Brilliantov＆Poeschel滾動(dòng)摩阻力矩表達(dá)式，由圖可見2種方法均能很好地預(yù)測(cè)這種分布.

圖7 軸線上軸向速度分布Fig.7 Velocities along spout axis

由圖還可看出，在(1)區(qū)的低位3個(gè)模型和文獻(xiàn)［11］的模擬結(jié)果相同，而到了高位Tsuji模型模擬結(jié)果與實(shí)驗(yàn)值符合更好;在(2)區(qū)Tsuji模型預(yù)測(cè)值整體與實(shí)驗(yàn)值最為接近，且好于文獻(xiàn)［11］的模擬結(jié)果，Syamlal＆O＇Brien模型對(duì)最大值的預(yù)測(cè)最好;在(3)區(qū)，3個(gè)模型給出的變化趨勢(shì)和實(shí)驗(yàn)值相符，數(shù)值上Arastoopour模型與實(shí)驗(yàn)值相差較大.圖7 (b)給出的是氣相軸向速度分布，由圖可以看出:z＜60 mm區(qū)域，Tsuji模型預(yù)測(cè)的速度最大;60 mm＜ z＜90 mm區(qū)域，Syamlal＆O＇Brien模型預(yù)測(cè)值最大;而z＞90 mm區(qū)域，3個(gè)模型預(yù)測(cè)值幾乎相同.3種模型預(yù)測(cè)的氣相速度不同是造成3個(gè)模型預(yù)測(cè)顆粒速度差異的主要原因，因?yàn)闅庀嘧饔昧κ菄娚鋮^(qū)的最大加速力，氣相速度越大顆粒獲的加速力就越大，從而得到較大的速度.

2.5 噴射區(qū)顆粒速度

圖8為噴射區(qū)不同床高處的顆粒軸向速度分布.

圖8 噴射區(qū)顆粒速度分布Fig.8 Particle velocities in spout

由圖可見實(shí)驗(yàn)結(jié)果呈現(xiàn)出類似立方函數(shù)的分布，vz隨床層高度的增加而增大，Tsuji模型的模擬結(jié)果和實(shí)驗(yàn)結(jié)果趨勢(shì)上符合較好，而其他2個(gè)曳力模型均有類似于San Jose［16］和趙香龍［17］提到的拋物線形式分布.由圖還可看出31.4 mm和45.6 mm兩個(gè)床層上Tusji模型結(jié)果與實(shí)驗(yàn)值符合的最好，60.8 mm床層上Syamlal＆O＇Brien模型結(jié)果與實(shí)驗(yàn)值符合最好，而在91.2 mm床層上模擬值均與實(shí)驗(yàn)結(jié)果存在較大差異.

2.6 環(huán)隙區(qū)內(nèi)顆粒速度

圖9為環(huán)隙區(qū)顆粒軸向速度分布，取豎直向下為正方向.由圖可見模擬與實(shí)驗(yàn)結(jié)果均顯示，在噴射區(qū)與環(huán)隙區(qū)交界處vz開始迅速增大，達(dá)到最大值后隨r的增大而減小.圖10為環(huán)隙區(qū)速度矢量分布，由圖9、10可見徑向上顆粒速度呈現(xiàn)出加速和減速兩個(gè)區(qū)域，出現(xiàn)這種兩區(qū)分布的原因是:在環(huán)隙區(qū)內(nèi)顆粒受到的主要作用力為接觸力，張勇等人［18］的研究結(jié)果表明顆粒碰撞頻率在環(huán)隙區(qū)存在一峰值，與圖9的顆粒速度分布相同，顆粒碰撞頻率反應(yīng)了顆粒間接觸次數(shù)的多少，接觸次數(shù)多，則顆粒所受接觸力可能就大，從而使其速度增大，反之亦然，由此可知顆粒碰撞頻率增加導(dǎo)致了vz增大.

圖9 環(huán)隙區(qū)顆粒速度分布Fig.9 Particle velocities in annulus

圖10 噴動(dòng)床局部顆粒矢量分布Fig.10 Vector fields of particle velocity in the local area of spouted bed

對(duì)于vz最大值的預(yù)測(cè)，Arastoopour模型的預(yù)測(cè)值與實(shí)驗(yàn)值相差最大，由91.2 mm至30.4 mm床高分別為實(shí)驗(yàn)值的1.47倍、1.4倍、1.36倍和1.15倍，與實(shí)驗(yàn)值符合最好的是Tsuji模型分別為1.12倍、1.14倍、1.12倍和1.14倍.雖然對(duì)于vz最大值的預(yù)測(cè)存在誤差，但3個(gè)模型的模擬結(jié)果均好于TFM［2］和無滾動(dòng)摩擦DEM［19］模擬結(jié)果的10倍和5倍誤差，這說明滾動(dòng)摩阻力矩的加入使速度最大值的預(yù)測(cè)更加精確.

圖11 91.2 mm時(shí)床高顆粒所受力矩比例分布Fig.11 Proportion of torque adding to particles at vessel height of 91.2 mm

圖12 91.2 mm時(shí)床高顆粒所受滾動(dòng)摩阻力矩分布Fig.12 The profile of torque at vessel height of 91.2 mm

圖11給出的是環(huán)隙區(qū)91.2 mm床高上顆粒所受滾動(dòng)摩阻力矩和切向接觸力矩比例分布，由圖可以看出，在環(huán)隙區(qū)滾動(dòng)摩阻力矩要大于切向接觸力矩，由此更能說明，在矩形噴動(dòng)床DEM模擬中加入滾動(dòng)摩阻力矩的必要性.圖12是91.2 mm床高近壁面處滾動(dòng)摩阻力矩分布，由圖可見隨r(r＞25 mm區(qū)域)的增大，顆粒所受滾動(dòng)摩阻力矩增大，Arastoopour模型的預(yù)測(cè)值最大，Syamlal＆O＇Brien模型預(yù)測(cè)值最小，Tsuji模型介于兩者之間，這與圖9中減速區(qū)顆粒速度梯度變化趨勢(shì)相同，即Arastoopour模型梯度最大，Tsuji模型次之，Syamlal＆O＇Brien模型最小，滾動(dòng)摩阻力矩與速度梯度兩者之間的關(guān)系類似于牛頓流體中粘性應(yīng)力和變形速率之間的關(guān)系，這說明在顆粒間相互作用中，滾動(dòng)摩阻力矩既可以是的驅(qū)動(dòng)力也可以是阻力.由圖9和圖12還可看出:vz峰值出現(xiàn)先于Mi峰值，這說明在相互影響關(guān)系上，速度梯度起主導(dǎo)作用.

3 結(jié)論

文中分別采用 Arastoopour、Di Felice、Syamlal＆O＇Brien和Tsuji模型結(jié)合DEM(pseudo－three－dimensional)方法對(duì)矩形噴動(dòng)床內(nèi)顆粒運(yùn)動(dòng)進(jìn)行了模擬，并與趙香龍等人的實(shí)驗(yàn)結(jié)果比較分析得出如下結(jié)論:

1)Arastoopour、Syamlal＆O＇Brien和Tsuji模型均可預(yù)測(cè)矩形噴動(dòng)床內(nèi)的流動(dòng)結(jié)構(gòu).

2)Arastoopour、Syamlal＆O＇Brien和Tsuji模型對(duì)顆粒軸向速度的預(yù)測(cè)在趨勢(shì)上均與趙香龍等人實(shí)驗(yàn)結(jié)果相符，其中Tsuji模型的預(yù)測(cè)值與實(shí)驗(yàn)結(jié)果最為接近.

3)在環(huán)隙區(qū)內(nèi)，顆粒碰撞頻率增加導(dǎo)致了vz的增大.滾動(dòng)摩阻力矩的加入使得環(huán)隙區(qū)內(nèi)顆粒軸向速度最大值的預(yù)測(cè)更加精確，而且模擬結(jié)果顯示該區(qū)內(nèi)滾動(dòng)摩阻力矩大于切向接觸力矩，因此在噴動(dòng)床DEM模擬中加入滾動(dòng)摩阻力矩很有必要.

4)滾動(dòng)摩阻力矩與顆粒速度梯度兩者之間是相互作用的，速度梯度起主導(dǎo)作用.

［1］金涌，祝京旭，汪展文.流態(tài)化工程原理［M］.北京:清華大學(xué)出版社，2001:360－362.

［2］DU Wei，BAO Xiaojun.Computational fluid dynamics (CFD)modeling of spouted bed:assessment of drag coefficient correlations［J］.Chemical Engineering Science，2006，61(5):1401－1420.

［3］LI Jie，KUIPERS J A M.Gas－particle interactions in dense gas－fluidized beds［J］.Chemical Engineering Science，2003，58(3－6):711－718.

［4］TSUJI Y，KAWAGUCHI T，TANAKA T.Discrete particle simulation of two－dimensional fluidized bed［J］.Powder Technology，1993，77:79－87.

［5］ZHONG W Q，XIONG Y Q.DEM simulation of gas－solid flow behaviors in spout－fluid bed［J］.Chemical Engineering Science，2006，61:1571－1584.

［6］ARASTOOPOUR H，PAKDEL P，ADEWUMI M.Hydrodynamic analysis of dilute gas－solids flow in a vertical pipe［J］.Powder Technology，1990:62(2):163－170.

［7］Di FELICE R.The voidage functions for fluid－particle inter－action system［J］.International Journal of Multiphase Flow，1994，20(1):153－159.

［8］SYAMLAL M，O＇BRIEN T J.Simulation of granular layer inversion in liquid fluidized beds［J］.International Journal of Multiphase Flow，1988，14(4):473－481.

［9］ZHOU Y C，WRIGHT B D，YANG R Y.Rolling friction in the dynamic simulation of sandpile formation［J］.Physica A—Statistical Mechanics and its Applications，1999，269: 536－553.

［10］XU B H，YU A B.Numerical simulation of the gas－solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics［J］.Chemical Engineering Science，1997，52:2785－2809.

［11］ZHAO X L，LI S Q，LIU G Q.DEM simulation of the particle dynamics in two－dimensional spouted beds［J］.Powder Technology，2008，184:205－213.

［12］LIM C J，MATHUR K B.Modeling of particle movement in spouted beds［M］.Cambridge:Cambridge University Press，1978:104－109.

［13］HE Y L，LIM C J，GRACE J R，ZHU J X，QZN S Z.Measurements of voidage profiles in spouted beds［J］.Canadian Journal of Chemical Engineering，1994，72(4):229－234.

［14］HE Y L，QIN S Z，LIM C J，GRACE J R.Particle velocity profiles and solid flow patterns in spouted beds［J］.Canadian Journal of Chemical Engineering，1994，72(8):561－568.

［15］MATHUR K B，EPSTEIN N.Spouted Beds［M］.New York:Academic Press，1974:265－270.

［16］SAN JOSé M J，MARTIN O，ALVAREZ S，IZQUIERDO M A，BILBAO J.Solid cross－flow into the spout and particle trajectories in conical spouted beds［J］.Chemical Engineering Science，1998，53(20):3561－3570.

［17］ZHAO X L，YAO Q，LI S Q.Effects of draft tubes on particle velocity profiles in spouted beds［J］.Chemical Engineering and Technology，2006，29(7):875－881.

［18］張勇，金保升，鐘文琪.基于顆粒尺度DEM直接數(shù)值模擬的噴動(dòng)流化床顆粒運(yùn)動(dòng)特性［J］.東南大學(xué)學(xué)報(bào):自然科學(xué)版，2008，38(1):110－115.

ZHANG Yong，JIN Baosheng，ZHONG Wenqi.Particlescale based DEM simulation of particle motion spout－fluid bed［J］.Journal of Southeast University:Natural Science E－dition，2008，38(1):110－115.

［19］KAWAGUCHI T，SAKAMOTO M，TANAKA T.Quasithree－dimensional numerical simulation of spouted beds in cylinder［J］.Powder Technology，2000，109(1):3－12.