噴動床DEM模擬曳力模型評價

朱衛兵，朱潤孺，邢力超，孫巧群

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

噴動床最早由加拿大的Mathur和Gishler發明，最初用于谷物干燥，如今噴動床已經被廣泛的應用到干燥、造粒、粉碎、煤燃燒和氣化、鐵礦石還原、油頁巖熱解、焦炭活化、石油熱裂化等領域［1］.在典型的噴動床型中，矩形噴動床與柱錐型噴動床相比，不但具有較低的最小噴動速度和最大噴動壓降的優點，還具有氣固接觸效率好、傳熱效率高的優點，且幾何結構簡單，易于制造.

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

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

1 數學模型

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

1.1 氣相模型

對于氣相場采用標準k－ε模型，考慮氣體的湍流運動.連續性方程和動量方程如下:

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

湍流動能和湍流動量耗散方程如下:

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

表1 模擬參數Table 1 Parameters for particle and fluid simulations

1.2 顆粒運動模型

顆粒的運動遵循牛頓第二定律，考慮重力、接觸力、流體曳力和壓力梯度力的作用.平動和轉動運動方程為

式中:

式中:dp、e、g、I、k、mi、Mi、nij、v、sij、δn，ij、δt，ij、μ、μr、η和ωi分別為顆粒直徑、顆粒恢復系數、重力加速度、轉動慣量、彈性系數、單個顆粒小球質量、滾動摩擦力矩、顆粒小球i和j之間的法向單位矢量、顆粒小球速度、顆粒小球i和j之間的切向單位矢量、顆粒小球i和j之間的法向變形、顆粒小球i和j之間的切向變形、滑動摩擦系數、滾動摩擦系數、阻尼系數和顆粒小球角速度.μ和μ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

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

2 結果與討論

2.1 曳力系數分析

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

圖2 曳力系數比較Fig.2 Comparison of drag force coefficients

2.2 流動結構

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

圖4是不同曳力模型下的流動結構模擬結果，由圖可以看出 Arastoopour、Syamlal＆O＇Brien和Tsuji模型的模擬結果均呈現出一種周期性噴動現象，周期大約為0.16 s，并且可以看出顆粒在噴泉區的運動呈現出均勻且對稱性較好的顆粒群運動狀態，這兩點與圖5的實驗結果符合很好，而Di Felice模型的模擬結果沒有出現這2種現象，所以不再對其進行進一步的討論.

圖3 不同曳力模型1 s時的流動結構Fig.3 The instantaneous flow patterns of different drag models at 1.0 s

圖4 流動結構模擬結果Fig.4 Simulated result of flow patterns

圖5 流動結構實驗結果Fig.5 Experimental result of flow patterns

2.3 空隙率分布

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

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

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

式中:Ds和np為噴射區直徑和噴射區給定床高處顆粒數目.

2.4 軸線上的顆粒速度

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

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

由圖還可看出，在(1)區的低位3個模型和文獻［11］的模擬結果相同，而到了高位Tsuji模型模擬結果與實驗值符合更好;在(2)區Tsuji模型預測值整體與實驗值最為接近，且好于文獻［11］的模擬結果，Syamlal＆O＇Brien模型對最大值的預測最好;在(3)區，3個模型給出的變化趨勢和實驗值相符，數值上Arastoopour模型與實驗值相差較大.圖7 (b)給出的是氣相軸向速度分布，由圖可以看出:z＜60 mm區域，Tsuji模型預測的速度最大;60 mm＜ z＜90 mm區域，Syamlal＆O＇Brien模型預測值最大;而z＞90 mm區域，3個模型預測值幾乎相同.3種模型預測的氣相速度不同是造成3個模型預測顆粒速度差異的主要原因，因為氣相作用力是噴射區的最大加速力，氣相速度越大顆粒獲的加速力就越大，從而得到較大的速度.

2.5 噴射區顆粒速度

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

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

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

2.6 環隙區內顆粒速度

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

圖9 環隙區顆粒速度分布Fig.9 Particle velocities in annulus

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

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

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

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

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

3 結論

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

1)Arastoopour、Syamlal＆O＇Brien和Tsuji模型均可預測矩形噴動床內的流動結構.

2)Arastoopour、Syamlal＆O＇Brien和Tsuji模型對顆粒軸向速度的預測在趨勢上均與趙香龍等人實驗結果相符，其中Tsuji模型的預測值與實驗結果最為接近.

3)在環隙區內，顆粒碰撞頻率增加導致了vz的增大.滾動摩阻力矩的加入使得環隙區內顆粒軸向速度最大值的預測更加精確，而且模擬結果顯示該區內滾動摩阻力矩大于切向接觸力矩，因此在噴動床DEM模擬中加入滾動摩阻力矩很有必要.

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

［1］金涌，祝京旭，汪展文.流態化工程原理［M］.北京:清華大學出版社，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直接數值模擬的噴動流化床顆粒運動特性［J］.東南大學學報:自然科學版，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.