基于不同曳力模型的鼓泡流化床CFD-DEM數(shù)值模擬與試驗研究

白玲，韓晨，徐云峰,周嶺*，張鈴杰，施衛(wèi)東,3

(1. 江蘇大學國家水泵及系統(tǒng)工程技術研究中心，江蘇 鎮(zhèn)江 212013； 2. 昌吉學院物理系，新疆 昌吉 831100； 3. 南通大學機械學院，江蘇 南通 226000)

鼓泡流化床廣泛應用于煤燃料脫硫脫硝、顆粒燃燒以及機械零部件的涂裝等國民生產領域[1].多年來，國內外相關學者開展了大量的試驗測量對流化床最低流化速度、床層高度、顆粒局部聚集等內部流動機理進行研究.但是由于氣固兩相流動的復雜性以及不確定性，試驗研究難以獲得微觀尺度的顆粒運動信息[2].

隨著計算機技術以及并行計算的發(fā)展，數(shù)值模擬方法成為研究流化床內部流動特性的主要手段.目前，對氣固兩相流動的數(shù)值計算主要有2種方法：第一種是歐拉-歐拉法，該方法將顆粒相看作擬連續(xù)相，基于傳統(tǒng)流體力學理論，分別解兩相的連續(xù)性方程；第二種方法是歐拉-拉格朗日法[3]， DEM屬于歐拉-拉格朗日法中的一種，該模型通過計算牛頓運輸方程以描述顆粒的移動和旋轉運動，同時能夠兼顧顆粒-顆粒、顆粒-壁面之間的碰撞，它突破了傳統(tǒng)的DPM方法僅能計算固相體積分數(shù)小于10%的稀相流動的瓶頸，因而該模型在稠密兩相流中得到廣泛應用.顆粒動力學理論認為曳力是氣固兩相流系統(tǒng)中最主要的力，當流體施加在顆粒上的曳力平衡顆粒自身重力時，顆粒懸浮[4].不同學者根據(jù)大量的試驗結果總結出多種不同的曳力模型，目前在DEM模型架構下，應用于流化床中氣固兩相流的曳力模型主要有Ergun，Wen & Yu，Syamal- O′Brien，Gidaspow，Di Felice，Huilin & Gidaspow這6種[5-8].由于試驗環(huán)境以及流化床內部流動的復雜性，在計算顆粒兩相流場時對曳力模型的選擇一直是困擾眾多學者的難題.鄭曉野等[9]針對低顆粒濃度范圍內曳力下降的問題通過光滑函數(shù)改進了曳力模型，并將模擬結果和Gidaspow，Syamal-O′Brien模型對比，發(fā)現(xiàn)改進的曳力模型能更好地描述兩相流動； VISURI等[10]依據(jù)歐拉雙流體模型對二維鼓泡流化床進行數(shù)值仿真，指出要依據(jù)系統(tǒng)性能選擇曳力模型以提高數(shù)值模擬精度.

由于缺乏對機理的理解，到目前為止，尚未建立令人滿意的通用方法來解釋這些機理并預測兩相流流動過程.文中基于DEM數(shù)值模擬以及準二維流化床高速攝影試驗，分別對比床層高度、空泡直徑，顆粒速度分布以及壓力波動，研究6種常用的曳力模型數(shù)值計算結果，以期為大尺度流化床數(shù)值計算和優(yōu)化設計提供理論參考.

1 數(shù)學模型

DEM模型屬于軟球模型，該模型在考慮顆粒相互干涉的同時，考慮顆粒相、連續(xù)相在計算微元中的體積分數(shù).假設εg為氣相體積分數(shù)，流體相的質量、動量守恒方程分別為

(1)

(2)

對于顆粒相的移動和轉動的描述主要通過解經典牛頓力學方程，即

(3)

(4)

文中采用彈性-阻尼器模型對碰撞過程中顆粒所受的短程作用力進行計算，該模型包含胡克彈性模型以及牛頓阻尼器模型.

2 試驗和數(shù)值模擬方法

為了直觀地觀察流化床內部的氣泡及顆粒流動特性，因而考慮將三維流化床簡化為準二維流化床，減少數(shù)值計算時間.在江蘇大學流體中心實驗室搭建了一個流化床高速攝影試驗臺，如圖1所示.

圖1 流化床高速攝影試驗Fig.1 Test rig for fluidized bed high speed photogra-phy experiment

空氣由空氣壓縮機加壓，經冷凍式干燥機干燥后通過高精度質量流量控制器進行流量調節(jié)，再進入流化床中.在長、寬、高分別為15,2,100 cm的準三維流化床內布置36 500個平均直徑為2.5 mm的玻璃球，初始堆積高度為17 cm，在流化床側邊2,22,40 cm高處分別布置3個壓力監(jiān)測點.壓力傳感器的型號是CY200,采樣頻率為1 kHz,流化床底部中心開一個孔，通入初始速度為350 L/min的空氣.質量流量為0.007 kg/s,為了防止顆粒被吹出，在流化床頂部布置了一張濾網(wǎng)，由于流化床較高，可以忽略濾網(wǎng)對整個系統(tǒng)壓力降的影響，通過正面放置的高速攝影機對流化床內鼓泡的運動以及顆粒的流動進行拍攝.高速攝影的拍攝頻率是1 000 s-1.

為了使得流體相的網(wǎng)格尺寸大于顆粒相的，在劃分流化床網(wǎng)格模型的時候，通過網(wǎng)格無關性分析最終確定網(wǎng)格尺寸為顆粒直徑的2.5倍，網(wǎng)格總數(shù)為22 880個的六面體網(wǎng)格模型，同時在模擬過程中流化床入口為中間位置采用了和9個圓形面積相等的矩形，矩形的長和寬分別為5.6 mm和5.0 mm，模擬過程中邊界條件設置為質量流量進口和壓力出口，為了描述湍流對流動的影響，采用標準k-ε湍流模型描述湍流對流動的影響， 進口處的湍流強度設定為5%.

3 討 論

3.1 鼓泡形態(tài)

圖2為在0～450 ms時間段內6種曳力模型的數(shù)值模擬結果與試驗結果的對比.在t=0 ms時刻，試驗中流化床顆粒床層初始堆積高度為17 cm,根據(jù)鼓泡的形態(tài)，在0

圖2 進口速度350 L/min下0～450 ms時間段內數(shù)值模擬結果(紅色)和試驗結果(黑色)對比Fig.2 Comparison of experimental(black) and numerical results (red ) from 0 to 450 ms with inlet speed of 350 L/min

3.2 床層高度

圖3為流化床床層高度的不同曳力模型與試驗對比結果.從圖中可以看出流化床床層高度的變化趨勢是先上升后下降，這主要是因為在0～300 ms階段，流化床中多余的氣體在流化床內部形成鼓泡，由于持續(xù)不斷的通氣，使得鼓泡不斷向上運動，推動了鼓泡頂部的顆粒層也向上運動，在300 ms時刻鼓泡發(fā)生了破裂，此時氣體逸出了顆粒床層，夾帶了部分顆粒彈射出流化床頂部自由表面，此時床層高度繼續(xù)向上運動，在350 ms時刻床層高度達到最大值，此時流化床顆粒床層內部的壓力與大氣壓一致，大部分的氣體都通過流化床頂部較大的顆粒間隙逸出流化床，流化床的高度開始下降，根據(jù)不同曳力模型的數(shù)值計算結果，Gidaspow曳力模型與試驗結果吻合度最高，Huilin & Gidaspow，Syamal-O′Brien和Ergun曳力模型都較高地預測了床層高度的變化趨勢，Di Felice曳力模型低估了床層高度的變化，而Wen & Yu曳力模型在0～350 ms階段高估了床層高度的變化，但是在鼓泡破裂后床層高度的回落過程中又低估了試驗值，因此綜合整個變化趨勢，Gidaspow曳力模型的表現(xiàn)最好，而在鼓泡形態(tài)方面表現(xiàn)突出的Huilin & Gidaspow曳力模型則高估了試驗值Hbed的變化.

圖3 床層高度的不同曳力模型數(shù)值模擬結果和試驗對比Fig.3 Comparison between experiment and simulation results of different drag models for bed height

3.3 鼓泡特性

床層高度和鼓泡面積是衡量鼓泡流化床內部流動性能的重要標尺.為了量化鼓泡的尺度，采用等效氣泡直徑公式計算，即

(5)

式中：Abub為鼓泡在高速攝影平面的投影面積.

圖4為流化床床層鼓泡面積的不同曳力模型與試驗對比結果，可以看到試驗中鼓泡面積基本上在300 ms時刻達到最大值，Gidaspow模型在200 ms時刻之前與試驗值比較吻合且小于試驗值，但在200 ms時刻之后預測值高于模擬值，但總的趨勢是一致的，同時可以看出基本上所有的曳力模型預測值在200 ms時刻之后都高于試驗值，其中Wen & Yu曳力模型與Di Felice曳力模型與試驗結果比較接近，Huilin & Gidaspow，Syamal-O′Brien和Ergun曳力模型都高估了鼓泡面積直徑的變化.

圖4 鼓泡直徑的不同曳力模型數(shù)值模擬結果和試驗對比Fig.4 Comparison of bubble diameter between expe-riment and simulation results of different drag models

從圖2中可以看到試驗中鼓泡發(fā)展到后期出現(xiàn)了第二級小氣泡，因此為了分析曳力模型對于第二級小氣泡的預測情況，直接采用小氣泡的面積來進行對比，主要是因為出現(xiàn)的第二級小氣泡在橫向方向尺度較大而在縱向方向尺度較小，其縱橫比較小，如果采用等效直徑的方法不能夠體現(xiàn)縱橫比小這一特性.從圖5可以看出Huilin & Gidaspow與Gidaspow曳力模型與試驗過程中第二級小氣泡的面積變化趨勢一致，其他曳力模型沒能夠符合小氣泡面積持續(xù)增大這一趨勢，特別是Di Felice曳力模型沒有預測出第二級小氣泡的現(xiàn)象，綜合發(fā)現(xiàn)Huilin & Gidaspow與Gidaspow曳力模型表現(xiàn)較好.

圖5 第二級小氣泡面積的不同曳力模型數(shù)值模擬結果和試驗對比Fig.5 Comparison of second-stage small bubble area between experiment and simulation results of different drag models

3.4 顆粒速度分布

圖6為6種曳力模型的軸向時均速度va對比，通過對比可以看出，軸向時均速度的分布為中心速度高，兩側速度低，而且時均速度大小的變化趨勢在橫向方向是從邊壁區(qū)到中心區(qū)先增加然后減小然后再增加，而且邊壁區(qū)的速度方向是向Z軸負向即顆粒的運動方向是向下運動的，這說明了顆粒的返混現(xiàn)象，由于鼓泡的運動以軸向運動為主，因此中心軸向速度較大.根據(jù)數(shù)值模擬結果Ergun模型的中心速度比其他曳力模型高許多，Gidaspow的中心速度比其他曳力模型稍低，其他4種模型的速度變化分布基本一致.

圖6 z=10 cm處z向速度的6種曳力模型的數(shù)值模擬結果對比Fig.6 Comparison of simulation results of the six drag models with z-direction velocity at z=10 cm

4 結 論

文中在自行搭建的流化床高速攝影試驗臺的基礎上，基于CFD-DEM方法，針對目前流化床中稠密氣固兩相流中Ergun，Wen & Yu，Syamal- O′Brien，Gidaspow，Di Felice以及Huilin & Gidaspow共6種曳力模型進行了數(shù)值模擬，并將模擬結果與試驗結果進行了對比.

總體上，在鼓泡形態(tài)的預測方面Huilin & Gidaspow與Gidaspow模型能夠對鼓泡的形態(tài)以及鼓泡破碎后出現(xiàn)的二級小氣泡進行比較準測的預測，但是在流態(tài)化特性的對比中，Gidaspow模型能夠在床層高度以及鼓泡的當量直徑方面具有明顯的優(yōu)勢，而Huilin & Gidaspow，Syamal-O′Brien模型以及Di Felice模型較高的估計了床層高度變化，因此，在基于CFD-DEM的流化床數(shù)值模擬中，Gidaspow模型能夠對稠密氣固兩相流動進行較為準確的預測，在空泡形態(tài)和壓力波動與試驗吻合較好，為工業(yè)中流化床的設計與優(yōu)化提供一定的理論參考.