固定/鼓泡床中煤瀝青球升溫過程的CFD-DEM模擬

薛嘉哲 鐘文琪 邵應(yīng)娟 謝立宇 李開喜

(1東南大學(xué)能源熱轉(zhuǎn)換及其過程測控教育部重點實驗室， 南京 210096) (2中國科學(xué)院山西煤炭化學(xué)研究所， 太原 030001)

煤瀝青的氧化不熔化過程是瀝青基球狀活性炭制備的重要環(huán)節(jié)之一，其目的是提高煤瀝青球的軟化點以確保其在隨后的炭化活化過程中不發(fā)生熔融變形[1].在此過程中，熱量在顆粒與氣體、顆粒與顆粒之間相互傳遞，使得瀝青球顆粒在合適的條件(反應(yīng)器、升溫速率)下緩慢升溫，最終達(dá)到目標(biāo)溫度.目前，由于瀝青球的機(jī)械強(qiáng)度較低，其氧化不熔化過程大多采用固定床反應(yīng)器，雖保證了瀝青球的結(jié)構(gòu)完整，但存在氧化效率低、能耗大、易燒料等缺點[2].流化床反應(yīng)器作為另一種工業(yè)常見設(shè)備，由于其具有氣固間接觸充分、床層溫度均勻等優(yōu)點而被廣泛用于化工生產(chǎn)領(lǐng)域[3-4].王永邦[5]成功采用流化床反應(yīng)器實現(xiàn)了瀝青球的高效氧化，從而表明在流化床反應(yīng)器上實現(xiàn)瀝青球氧化不熔化過程的可行性.然而，在實現(xiàn)瀝青球流化床氧化不熔化的應(yīng)用前，尚缺乏對瀝青球顆粒在流化床內(nèi)升溫傳熱過程的研究，瀝青球顆粒在反應(yīng)器內(nèi)熱量傳遞的規(guī)律及反應(yīng)器結(jié)構(gòu)對傳熱過程的影響尚不明確.

計算流體力學(xué)-離散單元法(computational fluid dynamics-discrete element method, CFD-DEM)能夠跟蹤氣固系統(tǒng)內(nèi)的每一個顆粒，可用于諸如流化床、煉鐵高爐等稠密氣固體系的模擬，并獲得豐富的運動與傳熱信息.Hou等[6]使用2D-CFD-3D-DEM方法研究了Geldart A類顆粒的流態(tài)化傳熱規(guī)律；Zhou等[7]使用2D-CFD-3D-DEM方法分析了不同操作氣速、顆粒導(dǎo)熱系數(shù)、楊氏模量等條件下，流化床氣固傳熱過程中不同傳熱模式的占比；Norouzi等[8]使用2D-CFD-DEM方法結(jié)合實驗測量手段研究了不同高徑比下，床層顆粒的混合及循環(huán)特性.上述研究能夠從顆粒尺度探尋氣固運動及傳熱的機(jī)理，從而表明了該方法在模擬稠密氣固體系上的優(yōu)勢.然而，以往研究出于計算經(jīng)濟(jì)性考慮，床體往往較薄(1～5倍粒徑)，其前后壁面被設(shè)置為周期性邊界條件，且尚無針對瀝青球顆粒在流態(tài)化系統(tǒng)中運動及傳熱過程的研究.

為此，本文利用CFD-DEM方法，建立了固定床/流化床內(nèi)瀝青球顆粒運動與傳熱過程的三維數(shù)理模型，對床層內(nèi)的氣固傳熱過程進(jìn)行了研究.首先將固定床和鼓泡床中氣固傳熱過程的數(shù)值模擬結(jié)果與實驗結(jié)果進(jìn)行對比，以驗證本文所采用計算模型的準(zhǔn)確性.隨后對不同流型(固定床/鼓泡床)、不同高徑比、厚寬比下床層內(nèi)的氣固流動及傳熱行為進(jìn)行研究，獲得瞬時氣固流動及傳熱結(jié)構(gòu)、床層升溫速率及溫度均勻性等信息，并從顆粒尺度揭示了床型結(jié)構(gòu)對瀝青球升溫過程的影響.

1 數(shù)值模型

1.1 顆粒相模型

煤瀝青球的運動過程通過牛頓第二定律描述，包括平動與轉(zhuǎn)動2種運動形式，控制方程分別表示如下：

(2)

式中，mi和Ii分別為顆粒的質(zhì)量和轉(zhuǎn)動慣量;vi和wi為其運動的線速度與角速度;Fc,ij為周圍顆粒j碰撞所產(chǎn)生的彈性接觸力；Fd,ij為黏性阻尼力;nc為與其相碰撞的周圍顆粒總數(shù);Tt,ij和Tr,ij分別為由切向力和滾動摩擦力產(chǎn)生的轉(zhuǎn)矩[9]；Fd為周圍流體對顆粒施加的曳力.

對于氣固系統(tǒng)，顆粒所受氣固作用力包括曳力、虛擬質(zhì)量力、馬格努斯升力、剪切提升力等.在本文的模擬中，氣固兩相間的密度差較大，虛擬質(zhì)量力、馬格努斯升力和剪切提升力對顆粒運動的影響較小[10]，可忽略.本文考慮的顆粒i受力包括與周圍顆粒j碰撞所產(chǎn)生的彈性接觸力Fc,ij、黏性阻尼力Fd,ij、與周圍流體之間的曳力Fd和重力mig.顆粒間接觸力的計算公式見文獻(xiàn)[11].

對于系統(tǒng)內(nèi)的傳熱過程，本文考慮3種傳熱機(jī)制：顆粒-顆粒間的導(dǎo)熱、顆粒-流體間的對流換熱以及顆粒-環(huán)境間的輻射換熱[6-7, 12-13].此外，本文中顆粒的畢渥數(shù)Bi<0.1，表征內(nèi)部導(dǎo)熱熱阻相較于外部熱阻很小，可以忽略不計，因此集總參數(shù)假設(shè)成立[7].顆粒的能量平衡可以通過如下零維方程來描述：

(3)

式中，Cp,i、Ti分別為顆粒定壓比熱和溫度；Qij為顆粒-顆粒間的導(dǎo)熱量；Qig為顆粒-流體間的對流換熱量；Qi,rad為顆粒-環(huán)境間的輻射換熱量.

1.2 氣相模型

氣相場基于歐拉框架，考慮氣體的湍流運動，使用k-ε湍流模型將其模化.詳細(xì)的控制方程表示形式見文獻(xiàn)[14-15].

1.3 氣固耦合子模型

1.3.1 曳力模型

對于本文所研究的稠密氣固系統(tǒng)，氣固曳力使用Gidaspow模型進(jìn)行計算[16]：

(4)

(6)

(7)

式中，βi為動量交換系數(shù)；αg為流體的體積分?jǐn)?shù);ug、ρg和μg分別表示周圍流體的速度、密度和動力黏度；Vi、di、CD、Rei分別為固體顆粒的體積、直徑、曳力系數(shù)、雷諾數(shù).

1.3.2 傳熱模型

本文考慮的傳熱機(jī)理包括：顆粒-顆粒導(dǎo)熱、顆粒-流體對流換熱[17-19].對于本文所研究的問題，其最高溫度低于600 K，因此輻射換熱并不顯著，可以將其忽略[13,20].導(dǎo)熱采用Chaudhuri等[21]提出的公式進(jìn)行計算，表達(dá)式如下：

(8)

(9)

(10)

式中，ki、kj為顆粒i、j的導(dǎo)熱系數(shù)；FN為作用在顆粒i上的法向力；E*和r*為兩顆粒的有效楊氏模量和幾何平均直徑；Ti、Tj為顆粒i和j的溫度；E為顆粒的真實楊氏模量；v為顆粒i的泊松比；ri、rj為顆粒i和j的半徑.

對于本文研究中所涉及到的氣固對流換熱，采用Li等[22]提出的公式進(jìn)行計算，表達(dá)式為

Qig=higAi(Tg-Tis)

(11)

(12)

(13)

式中，hig為對流換熱系數(shù)；Tg為流體溫度；Ai為顆粒的表面積；kg為流體導(dǎo)熱系數(shù)；Nui為努賽爾數(shù)；Pr為流體普朗特數(shù)；n為模型經(jīng)驗常數(shù)，取值為3.5.

1.4 模擬對象及計算參數(shù)

如圖1所示，本文的計算對象為三維矩形床，床料為直徑1.2 mm、密度1 200 kg/m3的瀝青小球[23].在計算前，溫度為300 K的瀝青球顆粒在床體上方隨機(jī)生成，并在重力作用下自由下落.待所有顆粒的動能為零后，床層底部形成一個初始顆粒堆積層.隨后，底部均勻送入溫度為600 K的流化風(fēng)，操作氣速不同，床層將呈現(xiàn)出固定床/鼓泡床的流動形態(tài).本文考慮了溫度對氣相物性參數(shù)的影響[7]，計算參數(shù)設(shè)置如表1所示.顆粒相最小流化風(fēng)速計算式見文獻(xiàn)[24].

圖1 數(shù)值模擬對象的幾何結(jié)構(gòu)

計算時，先對三維氣相場進(jìn)行計算.待收斂后，根據(jù)顆粒所在位置，計算其當(dāng)?shù)鼐W(wǎng)格下的流場參數(shù)(如ug、Tg等)，由此計算氣固作用力、對流通量.隨后計算顆粒碰撞力及顆粒間接觸導(dǎo)熱量，求解固相運動方程和能量方程從而更新顆粒的位置、速度及溫度信息.待其計算完成后，再將結(jié)果傳遞給流場求解器，氣固兩相連續(xù)交替耦合求解.

2 計算結(jié)果與分析

2.1 模型驗證

本文計算所采用的冷態(tài)CFD-DEM模型的驗證見文獻(xiàn)[14,25].由于本文研究涉及到固定床及鼓泡床中的氣固傳熱過程，因此需對本文所采用計算模型中的傳熱部分進(jìn)行驗證，所參考的固定床和鼓泡床傳熱實驗工況設(shè)置見文獻(xiàn) [26-27].

圖2(a)為不同操作氣速Uf及不同床料質(zhì)量m條件下，鼓泡床中的床層平均顆粒溫度隨時間變化情況.可以看出，操作氣速為1.54、1.71 m/s時的溫降曲線與Patil等[26]的實驗測量結(jié)果吻合較好，

表1 模擬中使用的參數(shù)

(a) Patil鼓泡床傳熱實驗[26]

(b) Collier固定床熱球?qū)嶒瀃27]

氣速為1.20 m/s時，冷卻速率比實驗稍慢.這與Wang等[28]、Patil等[29]的計算結(jié)果相似.Patil等[29]指出，當(dāng)操作氣速較低時，顆粒的瞬時停滯效應(yīng)會導(dǎo)致顆粒間存在更長的接觸時間.傳熱模型中未考慮這一效應(yīng)導(dǎo)致了低氣速下產(chǎn)生細(xì)微的誤差.圖2(b)給出了固定床中熱球溫度隨時間變化，通過比較發(fā)現(xiàn)計算與實驗結(jié)果吻合較好.值得注意的是，固定床中布置的9個熱球的溫度變化過程僅有微小差異，這與Zhou等[7]、Wei等[30]的CFD-DEM模擬結(jié)果相一致.

2.2 固定/鼓泡床內(nèi)瀝青小球的升溫過程

2.2.1 固定床/鼓泡床傳熱過程對比

由于本文中固定床與鼓泡床的床層尺寸、床料質(zhì)量、操作氣速均不相同，為便于分析比較兩者間的差異，本文對時間進(jìn)行無量綱化處理，即以實際時間t與床層熱平衡所需時間t0的比值τ=t/t0作為無量綱時間.其中床層熱平衡所需時間t0通過Hou等[6]提出的公式進(jìn)行估算.

圖3給出了不同時刻固定床與鼓泡床中瀝青小球的溫度分布情況.在固定床內(nèi)，不同位置處的顆粒溫升并不同步，底部顆粒最先被加熱至目標(biāo)溫度(600 K)，并在隨后的時間內(nèi)保持這一高溫.王永邦[5]指出，這種床層內(nèi)顆粒溫度的不均勻性分布容易導(dǎo)致燒料情況的發(fā)生，對最終的成球效果造成一定影響.在鼓泡床內(nèi)，床層底部氣泡的裹挾作用和上部顆粒的返混作用共同使得床層內(nèi)顆粒能夠均勻地得到加熱，床內(nèi)的溫度分布更均勻.

(b) 鼓泡床(ug=0.2 m/s)

圖4定量表征了固定床與鼓泡床中，瀝青小球的傳熱特性及其隨時間演化規(guī)律.從圖4(a)、(b)中可以看出，2種床型中發(fā)生氣固對流傳熱的位置不同：①固定床中的氣固對流換熱最早發(fā)生在床層底部，待底部氣固傳熱達(dá)到熱平衡后，對流換熱發(fā)生的位置逐漸上移.由于氣固間溫差降低，峰值水平也隨之下降.②而鼓泡床內(nèi)的對流換熱則始終發(fā)生在床層底部Y/H0=0～0.3這一范圍(約20 mm)內(nèi)，這與文獻(xiàn)[31-32]的研究結(jié)果較為一致.

由圖4(c)、(d)中2種床型下顆粒的溫度分布可以看出，固定床中不同位置處的顆粒溫升具有明顯的不同步性，靠近底部的顆粒溫升更快；而鼓泡床內(nèi)的軸向溫度分布較之固定床明顯更為均勻，整體上呈現(xiàn)出有效傳熱區(qū)-均勻溫度區(qū)-拋灑區(qū)的分布規(guī)律.其中，有效傳熱區(qū)表征床層底部劇烈的氣固對流傳熱作用，拋灑區(qū)表征上升氣泡對底部高溫顆粒的裹挾作用，而均勻分布區(qū)則表征鼓泡床內(nèi)良好的顆粒返混效果.

圖4(e)、(f)為2種床型下的床層平均溫度及溫度標(biāo)準(zhǔn)差隨時間演化情況.從圖中可以看出，2種床型的溫升斜率均隨時間推移而平穩(wěn)下降，這是由于氣固間溫差降低，傳熱的驅(qū)動力減弱所導(dǎo)致.對于溫度均勻性而言，兩者均表現(xiàn)為先升后降的變化趨勢.值得注意的是，由于鼓泡床內(nèi)良好的氣泡摻混效果，其溫度標(biāo)準(zhǔn)差明顯低于固定床，表明其具有更好的溫度均勻性.

圖5給出了固定床及鼓泡床中，某一跟蹤顆粒的軸向坐標(biāo)、溫度及對流傳熱通量隨時間變化情況.由圖5(a)中可以看出，固定床內(nèi)3個高度(Y/H0=0.02、0.6、0.9)處跟蹤顆粒的對流傳熱通量均呈現(xiàn)先增后減的變化規(guī)律.到達(dá)峰值的時間隨顆粒所處高度的增加而延后，峰值水平也隨之逐漸降低.由圖5(b)可知，鼓泡床內(nèi)的跟蹤顆粒溫升主要發(fā)生在其進(jìn)入底部區(qū)域的時間段內(nèi)(τ=0.02～0.03、0.052～0.058、0.075～0.076)，且顆粒越靠近底部布風(fēng)板，其溫升越大.當(dāng)顆粒處于床層中上部時，對流傳熱并不強(qiáng)烈，這與圖4(b)的對流通量分布結(jié)果相吻合.

2.2.2 高徑比的影響

從圖6可以看出，氣泡的生長、合并現(xiàn)象隨床層膨脹高度的增加而逐漸明顯，床層整體溫度水平也因熱容量的增加而有所降低.對于圖6(a)的淺床流化床，床層內(nèi)同時存在2個細(xì)長高溫帶狀結(jié)構(gòu)，這是由于堆積高度較低(H0/b=0.38)時，氣泡趨向于從床層中心的兩側(cè)生成，固相的內(nèi)循環(huán)模式發(fā)生變化所致[8].

(a) 固定床對流傳熱通量沿軸向分布

(b) 鼓泡床對流傳熱通量沿軸向分布

(c) 固定床溫度沿軸向分布

(d) 鼓泡床溫度沿軸向分布

(e) 固定床平均溫度及標(biāo)準(zhǔn)差隨時間變化

(f) 鼓泡床平均溫度及標(biāo)準(zhǔn)差隨時間變化

(a) 固定床

(b) 鼓泡床

圖7給出了不同高徑比下升溫速率及溫度均勻性的定量比較結(jié)果.從圖7(a)中可以看出，高徑比提高，床層升溫速率隨之下降，但其下降的幅度與高徑比提升的幅度之間并非線性關(guān)系.進(jìn)一步觀察圖7(b)可以發(fā)現(xiàn)，當(dāng)高徑比提高至0.95左右時，升溫速率的降低幅度已經(jīng)很小.這說明有效氣固傳熱距離并不隨堆積高度發(fā)生明顯變化.此外，從圖7(c)的溫度標(biāo)準(zhǔn)差隨時間變化情況可以看出，床層溫度分布隨高徑比的降低而變得更加均勻，但當(dāng)高徑比較大時，均勻性差距也并不明顯.

為從顆粒尺度解釋高徑比與顆粒運動、傳熱之間的內(nèi)在關(guān)系，圖8給出了3種高徑比下，底部中心處跟蹤顆粒的運動及傳熱情況.由圖可以看出，高徑比的增加提升了跟蹤顆粒的平均上升距離，從而使得顆粒進(jìn)入底部區(qū)域的頻率降低.這直接導(dǎo)致流傳熱通量峰值出現(xiàn)的頻率降低，并最終體現(xiàn)為圖8(c)中顆粒升溫速率的降低.

(a) H0/b=0.38

(b) H0/b=0.57

(c) H0/b=0.76

(d) H0/b=0.95

(e) H0/b=1.14

(a) 床層平均溫度

(b) 不同高徑比間的相對變化趨勢

(c) 溫度標(biāo)準(zhǔn)差

(a) 軸向坐標(biāo)隨時間變化

(c) 溫度隨時間變化

此外，跟蹤顆粒在傳熱初始階段出現(xiàn)了一定的溫降，且其時長隨高徑比的增加而增加(Δt3>Δt2>Δt1).對此現(xiàn)象，本文認(rèn)為：①通過前文分析可知，對流換熱主要發(fā)生在床層底部(見圖4(b))，當(dāng)流體經(jīng)過床層底部后，其溫度已與床層中上部冷顆粒的溫度相當(dāng)；②跟蹤顆粒在床層底部被空氣加熱后，受氣泡裹挾作用而向上運動.在此過程中，其周圍的流體溫度已低于受到加熱的跟蹤顆粒本身的溫度，因此其受到低溫空氣的冷卻而出現(xiàn)溫降.隨高徑比增加，顆粒在床層中上部區(qū)域的停留時間越長，其與低溫空氣發(fā)生對流換熱的時間越長，從而出現(xiàn)了更大的溫降.

2.2.3 厚寬比的影響

Geldart[33]指出，由于壁面效應(yīng)的存在，準(zhǔn)二維流化床在氣泡尺寸、上升速度、破碎頻率上與三維床存在一定差異.為進(jìn)一步探究鼓泡床厚寬比對氣固流動及傳熱的影響，本節(jié)將前文使用的床層厚寬比由0.2提升至0.6，以探究兩者在傳熱特性上的差異.

圖9給出了厚寬比為0.6時氣泡形狀(定義為αg=0.75等值面氣泡邊界)、中心剖面(Z=0)處的固相速度矢量Vp及Y方向動能隨時間的變化情況.

(a) 氣泡形狀及固相速度矢量

(b) Y方向顆粒動能隨時間變化情況

可以觀察到此時氣泡形狀呈現(xiàn)出強(qiáng)烈的三維特征.通過圖9(b)的顆粒動能可以看出，δ/b=0.6時的平均動能較之δ/b=0.2時提升約46.3%.這是由于厚寬比增加使得壁面效應(yīng)減弱，氣泡上升速度提高[34-35].此外，不同厚度下的顆粒Y方向速度均隨時間稍有增加，這是由于溫度變化導(dǎo)致的氣相物性參數(shù)變化所致.

圖10給出了2種厚寬比下床層平均顆粒溫度及標(biāo)準(zhǔn)差的時間演化情況.利用Hou等[6]提出的平均升溫速率表達(dá)式對圖10(a)中的升溫速率進(jìn)行統(tǒng)計可以發(fā)現(xiàn)，三維床相較于準(zhǔn)二維床，其升溫速率提高了24.2%.對此現(xiàn)象本文認(rèn)為：氣泡上升速度的增加能夠更快地將高溫顆粒帶離底部區(qū)域，同時床層中上部的低溫顆粒將返混至底部區(qū)域得到加熱，從而提高了底部區(qū)域氣固間的溫差.當(dāng)忽略顆粒間導(dǎo)熱時，由顆粒升溫過程的熱平衡方程

可知，氣固溫差越大，傳熱的驅(qū)動力越大，顆粒溫升也越大.

從圖10(b)可以看出，在計算初始階段，底部空氣的加熱打破了初始床層的溫度均勻性，使得溫度標(biāo)準(zhǔn)差從零迅速提高.且床層厚度增加導(dǎo)致的升溫速率加快提高了標(biāo)準(zhǔn)差的峰值.在隨后的標(biāo)準(zhǔn)差下降過程中，三維床的下降幅度略高于準(zhǔn)二維床，這可能是由于三維床更為強(qiáng)烈的氣泡摻混作用所致.

(a) 床層平均溫度隨時間演化情況

(b) 溫度標(biāo)準(zhǔn)差隨時間演化情況

3 結(jié)論

1) 在固定床中，瀝青顆粒溫升從床層底部逐漸向上傳遞，顆粒沿軸向存在明顯的溫度梯度，床層底部顆粒長時間處于高溫狀態(tài)；而采用鼓泡床反應(yīng)器，底部被加熱的顆粒能夠及時被上升氣泡帶離底部區(qū)域，避免了局部高溫區(qū)域(燒料現(xiàn)象)的形成.

2) 當(dāng)高徑比提高至0.95以上時，其對鼓泡床溫升的影響較小；隨高徑比的提高，跟蹤顆粒在床層上方受到低溫空氣冷卻的時間增加，從而出現(xiàn)了更大的溫降.

3) 隨厚寬比的增加，壁面效應(yīng)的影響減弱，氣泡上升速度、顆粒動能均有所增加，且床層顆粒的內(nèi)循環(huán)速率加快，床層底部氣固間溫差有所增加，床層顆粒升溫速率提升約24.2%.