基于CFD-DEM法的矩形噴動床與柱錐形噴動床顆粒噴動特性模擬

雷 琨，宋海萍，高亞平，徐 慶, ，吳 龍, ，王瑞芳,

(1. 天津科技大學(xué)機械工程學(xué)院，天津 300222；2. 低碳綠色過程裝備國際聯(lián)合研究中心，天津 300222)

微波噴動床是集微波加熱與噴動床優(yōu)勢于一體的微波聯(lián)合干燥裝置．其既具有微波加熱速度快的優(yōu)點，又能利用噴動床中物料的隨機運動改善微波加熱不均勻性的缺點．傳統(tǒng)噴動床為柱錐形結(jié)構(gòu)，但在微波噴動床結(jié)構(gòu)設(shè)計時，為了保證噴動床內(nèi)電磁場分布均勻性，矩形腔更適合于微波均勻加熱．本文作者已在前期基于電場分布均勻性對單饋口矩形微波噴動床結(jié)構(gòu)進行了設(shè)計[1]，采用多物理場軟件COMSOL Multiphysics 模擬研究了噴動床內(nèi)部的電場分布情況，得到電場強度大且電場分布均勻的單饋口矩形微波噴動床結(jié)構(gòu)，其長×寬×高為427 mm×427 mm×1 000 mm．針對該結(jié)構(gòu)，本文采用離散元氣固耦合法(CFD-DEM)對矩形噴動床內(nèi)小麥顆粒噴動特性進行模擬研究．

CFD-DEM 法可以同時模擬研究流場和顆粒的動力學(xué)行為，包括流場和顆粒本身的動力學(xué)特性以及流場與顆粒相互作用的動力學(xué)特性，其研究結(jié)果更具可靠性和真實性．近年來，CFD-DEM 法被越來越多地應(yīng)用于流態(tài)化技術(shù)方面研究，黎明等[2]通過CFDDEM 法模擬了流化床的工作機理；Kafui 等[3]采用CFD-DEM 法模擬了顆粒在流化床內(nèi)的流動特性；匡世波[4]開發(fā)了三維CFD-DEM 軟件，并利用該軟件預(yù)測了氣力輸送過程；田鳳國[5]通過數(shù)值模擬對內(nèi)循環(huán)流化床氣固兩相流進行了數(shù)值模擬，得出顆粒流動規(guī)律．Khawaja[6]通過CFD-DEM 法對流化床內(nèi)的聲波與顆粒流化的影響進行了研究，并得出系統(tǒng)中的阻尼受固相黏度的影響．Kerst 等[7]利用CFD-DEM 法對流化床內(nèi)結(jié)晶水的流化特性進行研究，證明了流化床內(nèi)結(jié)晶水流化特性簡化模型的可靠性．Fries 等[8]用CFD-DEM 法對噴霧造粒在流化床中物性變化進行研究，結(jié)果表明不同噴霧方式及不同材料對噴霧造粒過程都會產(chǎn)生很大的影響．Ebrahimi 等[9]通過CFDDEM 法研究大型噴動流化床的噴動特性，提出一種縮放方法，既可以減少計算時間又可以保證計算的精確性．Sutkar 等[10]用CFD-DEM 法研究液態(tài)噴射對噴動流化床內(nèi)顆粒與液滴之間的相互作用以及噴動流化特性的影響，并得到了物料水分的分布情況．Salikov 等[11]通過CFD-DEM 法研究了菱形噴動床內(nèi)氣流對噴動穩(wěn)定性的影響，在進口氣速和壓降等方面都與實驗保持了一致性．Liu 等[12]用CFD-DEM法對噴動床內(nèi)的高密度核燃料顆粒在不同密度下顆粒的床層表現(xiàn)進行研究，得出了顆粒在噴動床中的最小噴動速度、進氣速度、穩(wěn)定噴動速度、壓降以及燃料涂層隨顆粒密度的變化．

傳統(tǒng)的柱錐形噴動床結(jié)構(gòu)方面的設(shè)計理論與經(jīng)驗相對成熟，而矩形噴動床的相關(guān)研究則較少．本文通過將矩形噴動床與柱錐形噴動床在混合時間、單顆粒循環(huán)時間、顆粒運動及氣體壓力(表壓)等方面進行對比，探索矩形噴動床噴動特性的優(yōu)缺點，以便對其結(jié)構(gòu)進行改進．

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

在流固耦合系統(tǒng)中，連續(xù)的流體相與離散顆粒相之間相互影響，在建立模型時要采用雙向耦合的方法處理兩相間的相互作用．CFD-DEM 模型用歐拉方法描述流體相和離散顆粒相的運動．該模型中，流體運動用兩相耦合的Navier-Stokes 方程描述，即在單相Navier-Stokes 方程中加入空隙率以及動量交換源項；顆粒的運動則考慮顆粒間相互碰撞作用力以及流體對顆粒的作用力，通過直接跟蹤顆粒的運動軌跡進行描述．

1.1 氣固耦合流體動力學(xué)控制方程

流體在單獨流動的過程中要遵循三大守恒定律，即質(zhì)量守恒定律、動量守恒定律和能量守恒定律，其控制方程是三大守恒定律的不同形式的表達．CFDDEM 耦合時存在顆粒相，因此在流體動力學(xué)研究過程中要考慮空隙率的影響，其方程一般的通用形式為

式(1)中從左至右依次是瞬態(tài)項、對流項、擴散項、源項[13]，其展開形式為

式中：ε 為空隙率；φ 為通用變量；Γ 為廣義擴散系；S為廣義源項．

表1 給出了控制方程中φ、Γ 和S 與每種特定方程的對應(yīng)關(guān)系．

表1 通用控制方程中各符號的表達式 Tab. 1 Expressions of symbols in general governing equation

1.2 氣體與顆粒的作用力方程

氣體與顆粒之間存在浮力、壓力梯度力、曳力、Magnus 力、Basset 力、Saffman 力等多種類型相互作用力．但對于密相顆粒或者無黏性的干顆粒，且顆粒密度遠(yuǎn)大于氣體密度時，與曳力相比，其他作用力均較小，可忽略．因此在計算顆粒與氣體之間存在的雙向耦合力時，可以只考慮曳力作用．單顆粒上受氣體作用力[14]為

式中： FDi為當(dāng)前時步單顆粒i 所受的氣體作用力；為顆粒i 的等效直徑；ρg為氣體密度；ug為顆粒i 所在網(wǎng)格的氣體速度，可通過求解氣固耦合的氣相控制方程求出；up為當(dāng)前時步顆粒i 的運動速度，可由CFD-DEM 耦合計算求出；ε 為空隙率，即網(wǎng)格內(nèi)氣體所占體積的百分比；Cd為單顆粒曳力系數(shù)(顆粒阻力系數(shù))．

1.3 接觸模型

顆粒的接觸模型采用Hertz-Mindlin 無滑移接觸模型，該模型可以準(zhǔn)確地計算出兩個顆粒在彈性接觸時的接觸作用力．

設(shè)兩球形顆粒發(fā)生碰撞接觸，其法向疊合量為

式中：R1、R2分別為兩球形顆粒的半徑；r1、r2分別為兩球心的位置矢量．

合力和合力矩通過顆粒碰撞過程中切向接觸力和法向接觸力的疊加來計算．顆粒間的法向接觸力

其中

式中： E*為等效彈性模量；E1、E2分別為兩顆粒的彈性模量；R*為等效粒子半徑；分別為兩顆粒的泊松比．

顆粒間的切向接觸力

其中

式中：St為切向剛度；δ 為切向疊合量；G*=，為等效剪切模量，其中G1、G2分別為兩顆粒的剪切彈性模量．

2 幾何體模型和顆粒模型

矩形噴動床與柱錐形噴動床的幾何體模型見圖1，計算域為三維空間．小麥顆粒采用球形顆粒模型，由于小麥顆粒為不規(guī)則形狀，且接近橢圓形，因此在測量時采用了等球平均粒徑．幾何體與物料相關(guān)參數(shù)及模型參數(shù)見表2．

圖1 矩形噴動床與柱錐形噴動床模型 Fig. 1 Model of rectangular spouted bed and cylindrical spouted bed

表2 模型參數(shù) Tab. 2 Parameters of the model

3 結(jié)果與討論

3.1 模型驗證

為驗證CFD-DEM 耦合模型用于分析矩形噴動床噴動特性的正確性與可靠性，作者對前期設(shè)計的長、寬均為 300 mm、高為 800 mm，床層高度為120 mm 的矩形噴動床[17]進行了模擬，并將利用本模型模擬得到的最小噴動速度、循環(huán)時間和完全混合時間等參數(shù)與文獻[17-19]的數(shù)據(jù)進行對比分析．

3.1.1 最小噴動速度對比

最小噴動速度實驗結(jié)果[17]與本文模型模擬的最小噴動速度結(jié)果對比見圖2．最小噴動速度分別為0.17 m/s(文獻)和0.14 m/s(模擬)，模擬值與實驗值相差17.6%，吻合度較好．

圖2 最小噴動速度對比 Fig. 2 Comparison of the minimum spouting velocity

最小噴動速度通常也用Mathur-Gishler 提出的經(jīng)驗公式[20]確定，用此經(jīng)驗公式計算的結(jié)果為0.15 m/s，模擬值與經(jīng)驗公式計算值相差7.1%，具有較好的吻合度．

3.1.2 單顆粒循環(huán)時間對比

將顆粒從噴動床某一高度出發(fā)后經(jīng)過一個循環(huán)(即顆粒經(jīng)過環(huán)流區(qū)、噴動區(qū)和噴泉區(qū))再次回到該高度的時間定義為單顆粒循環(huán)時間．模擬得到的單顆粒循環(huán)時間為1.85 s，模擬結(jié)果與Luo 等[18]得到的矩形噴動床內(nèi)單顆粒循環(huán)時間0～2.5 s 吻合．

3.1.3 混合時間對比

顆粒在噴動床內(nèi)達到完全混合的時間也是表征噴動床噴動效果的因素．混合時間主要通過混合指數(shù)判斷．將顆粒標(biāo)定為兩種不同顏色，并等量分為上下兩層放置于矩形噴動床中．模擬結(jié)束后，取出一定量的顆粒，計算其中一種顏色顆粒的質(zhì)量占總顆粒質(zhì)量的質(zhì)量分?jǐn)?shù)，將其定義為混合指數(shù)M．當(dāng)混合指數(shù)為0.5 時認(rèn)為混合達到均勻．

模擬中，混合時間為5.2 s 時混合指數(shù)達到0.5，說明顆粒已經(jīng)混合均勻(圖3)．王喆[17]的研究表明，雖然有一定的波動性，但是當(dāng)混合時間為6 s 時，顆?；净旌暇鶆?，此結(jié)果與本文模擬結(jié)果吻合．朱潤儒[19]的研究發(fā)現(xiàn)，經(jīng)過7 s 后顆粒達到了充分混合，與本文模擬結(jié)果接近．

綜上，通過對矩形噴動床的最小噴動速度、單顆粒循環(huán)時間和混合時間的模擬及與文獻結(jié)果進行對比可知，本文采用CFD-DEM 耦合模塊對顆粒在噴動床中噴動特性的研究是可靠的．

圖3 混合指數(shù)模擬與文獻[17]對比 Fig. 3 Comparison of simulated mixed index and literature

3.2 矩形噴動床與柱錐形噴動床對比

物料的靜止床層高度為170 mm，對噴動速度分別為ums、1.1 ums、1.2 ums時的噴動情況進行觀察．顆粒的噴動程度隨著噴動速度的增加而增加．當(dāng)us=1.2 ums時，在噴動過程中會出現(xiàn)少量顆粒飛出噴動床的現(xiàn)象．因此，在矩形噴動床與柱錐形噴動床對比研究中，采用的噴動速度為1.1 ums．

3.2.1 噴動床內(nèi)氣體壓力

圖4 為us=1.1 ums、t=5 s 時，兩種噴動床內(nèi)氣體壓力的分布情況(XZ 平面，Y=0 mm)．兩種噴動床的氣體壓力分布均呈現(xiàn)出中心區(qū)域壓力較大，越接近壁面氣體壓力越小，且呈對稱分布．

圖4 矩形噴動床和柱錐形噴動床氣體壓力對比 Fig. 4 Gas pressure comparison of rectangular spouted bed and cylindrical cone spouted bed

圖5 所示為Z=369 mm(即柱錐交界面)、569 mm和769 mm 時，XZ 平面(Y=0 mm)上的床層壓力沿X軸的變化情況．可以看出：隨著高度的增加，床層壓力遞減，柱錐形噴動床的壓力值小于矩形噴動床內(nèi)的壓力值．

圖5 矩形噴動床和柱錐形噴動床在不同平面的床層壓力對比(us＝1.1 ums) Fig. 5 Pressure comparison of rectangular spouted bed and sylindrical spouted bed at different planes(us＝1.1 ums)

3.2.2 混合特性

圖6 為Z=120 mm 和Z=369 mm(即柱錐交界面)時，混合指數(shù)隨時間的變化過程．可以看出：在0～2.8 s，柱錐形噴動床的混合速率略高于矩形噴動床的混合速率，這是因為矩形噴動床內(nèi)的氣體壓力較大，氣流速度比柱錐形噴動床低，顆粒受氣流作用，其加速度要比柱錐形噴動床略低；在2.8 s 以后，矩形噴動床和柱錐形噴動床內(nèi)顆粒均達到了混合均勻狀態(tài)，即 M=0.5．一般加熱過程遠(yuǎn)遠(yuǎn)大于2.8 s，因此，前期的混合特性差異對加熱過程影響很小，可以忽略，即矩形噴動床與柱錐形噴動床在顆?；旌咸匦苑矫婊疽恢拢?/p>

圖6 us＝1.1 ums 顆?；旌咸匦詫Ρ?Fig. 6 Comparison of particle mixing characteristics at us＝1.1 ums

3.2.3 單顆粒循環(huán)時間

為了保證顆粒循環(huán)時間的準(zhǔn)確性，分別在Z=40、80、120、170 mm 的XY 平面內(nèi)隨機選取50 個顆粒，通過分析顆粒在運動過程中高度的變化，每相鄰兩次位于同一高度時為顆粒運動一個循環(huán)，記錄其循環(huán)時間．通過統(tǒng)計10 s 內(nèi)顆粒的循環(huán)次數(shù)和循環(huán)時間，計算出單顆粒平均循環(huán)時間．

不同高度顆粒的運動過程均呈現(xiàn)出良好的循環(huán)特性(圖7)，而且柱錐形噴動床與矩形噴動床內(nèi)保持了較好的一致性．

圖7 us＝1.1 ums 時的矩形噴動床與柱錐形噴動床單顆粒循環(huán)時間對比 Fig. 7 Comparison of cycle time between rectangular spouted bed and sylindrical spouted bed at us＝1.1 ums

表3 為不同噴動床的不同位置單顆粒循環(huán)時間比較．矩形噴動床的單顆粒平均循環(huán)時間為1.9 s，柱錐形噴動床的單顆粒平均循環(huán)時間為1.7 s，矩形噴動床的內(nèi)部壓力較大導(dǎo)致顆粒運動速度較小，其循環(huán)時間略長．但是，矩形噴動床在不同位置的單顆粒循環(huán)時間相差較小，說明矩形噴動床內(nèi)氣流速度比柱錐形噴動床更加平穩(wěn)．

表3 不同位置單顆粒循環(huán)時間比較 Tab. 3 Comparison of cycle time of single particle in different positions

3.2.4 顆粒運動情況

圖 8 表示進口氣速 us=1.1 ums，Z=369、569、769 mm 時，XZ 平面(Y=0 mm)上的顆粒軸向速度沿X 軸的分布規(guī)律．由此看出：越接近噴射區(qū)的顆粒軸向速度越大，越接近壁面處的顆粒軸向速度越??；隨著軸向位置的升高，因顆粒的動能減小使得顆粒的速度呈現(xiàn)下降趨勢．矩形噴動床與柱錐形噴動床的顆粒速度場分布基本一致．

圖8 矩形噴動床和柱錐形噴動床中顆粒軸向速度對比(us＝1.1 ums) Fig. 8 Axial velocity comparison of particles in different planes of rectangular spouted bed and cylindrical spouted bed(us＝1.1 ums)

考慮到矩形噴動床中流場分布受其結(jié)構(gòu)形狀的影響，有些位置的流速較小或者顆粒受力處于平衡狀態(tài)，可能會出現(xiàn)顆粒停滯現(xiàn)象，不能更好地參與到噴動循環(huán)中，物料的集聚會影響加熱的均勻性．本文在矩形噴動床內(nèi)選取噴動床底部(Z=10 mm)和中間(Z=80 mm)兩個平面中靠近壁面的邊緣位置分別隨機選取100 顆顆粒，觀察其位置變化情況，從而確定邊緣顆粒在噴動過程中是否出現(xiàn)停滯情況．如圖9 所示，模擬開始時，顆粒受入射氣體的影響移動到一個較高的位置，隨著穩(wěn)定噴動過程的形成，顆粒運動范圍也逐漸穩(wěn)定，但與圖7 相比，顆粒運動范圍明顯減?。ㄟ^分析矩形噴動床在XY 平面(Z=80 mm)的壓力分布，發(fā)現(xiàn)在矩形噴動床棱角處氣體壓力最大，如圖10 所示．因此，在矩形噴動床底部(Z=10 mm)和中間(Z=80 mm)的棱角處隨機取單顆粒分析其運動情況，如圖11 所示，發(fā)現(xiàn)在某一個時間段，顆粒的Z 坐標(biāo)基本不變或變化很小，說明顆粒出現(xiàn)停滯現(xiàn)象．因而，在矩形噴動床的結(jié)構(gòu)設(shè)計時要充分考慮棱角處顆粒運動問題．

圖9 靠近噴動床壁面位置顆粒的位置變化對比 Fig. 9 Comparison of position changes of particles near the pouted bed wall

圖10 矩形噴動床內(nèi)氣體壓力分布(Z=80 mm) Fig. 10 Distribution of gas pressure in rectangular spouted bed(Z=80 mm)

圖11 矩形噴動床棱角處顆粒位置隨時間的變化 Fig. 11 Change of the position of particles at the edges and corners of rectangular spouted bed with time

4 結(jié) 語

本文通過CFD-DEM 法對矩形噴動床與傳統(tǒng)柱錐形噴動床中小麥顆粒噴動特性進行對比研究，發(fā)現(xiàn)矩形噴動床內(nèi)的壓力較大．針對本文研究的兩種噴動床結(jié)構(gòu)和尺寸，顆?；旌暇鶆驎r間均為2.8 s．矩形噴動床的單顆粒平均循環(huán)時間為1.9 s，柱錐形噴動床的單顆粒平均循環(huán)時間為1.7 s，相差較?。畠煞N噴動床的顆粒速度場分布基本一致，顆粒的平均速度差異較小，但通過對靠近噴動床壁面位置顆粒運動研究，發(fā)現(xiàn)矩形噴動床在棱角處的顆粒運動隨機性較差，主要因為矩形噴動床棱角處氣體壓力較大．

研究說明，矩形噴動床和柱錐形噴動床的噴動特性基本一致，其結(jié)構(gòu)設(shè)計在保證微波場均勻分布的前提下可以利用柱錐形噴動床的設(shè)計理論．需要注意的是，對于棱角處顆粒隨機運動問題，需要在噴動床中內(nèi)設(shè)結(jié)構(gòu)方面進一步改進．