曳力模型及鏡面系數(shù)對鼓泡床氣固流場的影響

白 勇，司 慧，趙思元，李 萍，賀 堯

（北京林業(yè)大學 工學院，北京 100083）

流化床內(nèi)氣固接觸比表面積大、反應(yīng)速率高、安全可靠，普遍應(yīng)用于燃燒、氣化、冶金及催化裂化等工業(yè)領(lǐng)域［1-3］，其中氣固兩相流的流動、輸送、傳熱傳質(zhì)和反應(yīng)過程等均具有一定規(guī)律。因此對鼓泡床內(nèi)氣固兩相流的運動特性、流場分布與氣固流動梯度場等研究具有重要意義［4-5］。在氣固流場中，顆粒受到最大的力是曳力，對其運動影響最大，其他作用力如浮力、重力等數(shù)值較小，經(jīng)常忽略不計［6-9］。張佳寶等［10］考察了4種不同機理的曳力模型對內(nèi)流反應(yīng)器中氣泡運動規(guī)律及其下降管中氣含率的影響，獲悉了DBS-Local曳力模型下下降管中的氣含率分布云圖。鄭曉野［11］利用Fluent軟件發(fā)現(xiàn)不同尺寸網(wǎng)格對鼓泡床流態(tài)化效果模型結(jié)果影響極大，利用所開發(fā)的曳力模型可較好地計算兩相流的曳力系數(shù)，能夠準確地預(yù)測床內(nèi)渦的分布。顆粒動力學理論常用來研究氣固兩相流之間的相互作用［1，12］。其中，在相同曳力模型條件下，顆粒與壁面摩擦動能耗散用無量綱數(shù)的鏡面系數(shù)（ψ）表征［13］，是氣固兩相流中邊界條件設(shè)置的關(guān)鍵量，ψ愈大，顆粒的動能耗散加劇，鼓泡床內(nèi)壁面處顆粒團聚現(xiàn)象愈嚴重［14-15］。綜上所述，對于流化床的曳力模型研究大多集中在固體停留時間和氣固兩相分布，而對于床層膨脹率、固含率、顆粒速度等較少，同時發(fā)現(xiàn)鮮有關(guān)于鼓泡床ψ的研究。原理上鼓泡床作為一種流化床，內(nèi)部氣固流場相對較易觀察。

本工作利用顆粒動力學理論和歐拉雙流體模型對鼓泡床內(nèi)氣固流場進行了數(shù)值模擬，考察了床層膨脹率、壓降以及顆粒速度特性的影響。

1 實驗裝置及數(shù)學模型

1.1 實驗裝置

為研究鼓泡床內(nèi)氣固流場特性，搭建了如圖1所示的冷態(tài)實驗裝置。

圖1 鼓泡流化床冷態(tài)實驗裝置Fig.1 Bubbling fluidized bed cold state experimental device.

鼓泡流化床內(nèi)徑0.1 m，高0.5 m，底部進氣，頂部排氣，材料選用有機玻璃，便于觀察內(nèi)部流化現(xiàn)象，在鼓泡床底部安裝有布風板，使用空氣壓縮機提供流化載氣，并通過氣體流量計檢測反應(yīng)器入口氣速，U型壓差計兩端設(shè)置測壓點，以測量床層壓降，利用測試儀記錄鼓泡床內(nèi)氣固流場與固含率分布。

1.2 曳力模型

采用De Felice，Gidaspow，Wen-Yu等3種曳力模型，考察不同曳力模型對鼓泡床內(nèi)的氣固流場的影響［10］。

1.3 幾何模型與網(wǎng)格劃分

利用Gambit軟件，構(gòu)造鼓泡流化床作為控制體，底部進氣，頂部排氣，并填入球形固體顆粒作為填料。為驗證網(wǎng)格獨立性并保證預(yù)測結(jié)果的可靠性，將控制體劃分3種總網(wǎng)格數(shù)量不同（10 800，64 000，120 000）的非結(jié)構(gòu)化六面體網(wǎng)格，端面為發(fā)散式四邊形。其中軸向網(wǎng)格劃分為等間距線段；氣體入口區(qū)域網(wǎng)格劃分采用三角形網(wǎng)格，劃分的網(wǎng)格數(shù)量較多；氣體出口區(qū)域網(wǎng)格劃分采用網(wǎng)格尺寸較大的四邊形網(wǎng)格。對控制體進行模擬，結(jié)果見圖2。由圖2可知，不同網(wǎng)格數(shù)量具有無關(guān)性。考慮到計算精度和時間，最后采用中等數(shù)量網(wǎng)格（64 000）。

圖2 網(wǎng)格獨立性檢驗Fig.2 Grid independence testing.

1.4 計算方法及邊界條件

選用Fluent軟件三維雙精度求解器，歐拉雙流體模型，算法選用相耦合SIMPLIC，適用于多相流。對接模式為非穩(wěn)態(tài)，求解最小時間步長為10-3s，允許的最大迭代殘差為10-4。模擬計算過程根據(jù)不同模型及參數(shù)設(shè)定分為不同工況，每種工況持續(xù)10 s。

在鼓泡流化床入口氣體采用均勻布氣，出口設(shè)置為出口流動邊界，其他面設(shè)置為壁面，相關(guān)的算例設(shè)置見表1。

2 計算結(jié)果與討論

結(jié)合控制體的殘差曲線和縱截面氣固動量互換計算曲線變化趨勢，當變化幅度小于10-4并在2 s后基本達到穩(wěn)態(tài)，圖3～圖9均采用控制體達到穩(wěn)態(tài)情況下（3～10 s）的平均值。

表1 模型參數(shù)設(shè)置Table 1 Model parameter settings

2.1 不同曳力模型的影響

圖3為應(yīng)用3種曳力模型的床層膨脹率隨表觀氣速變化的預(yù)測值及實驗值。由圖3可知，4條曲線變化趨勢基本一致，床層膨脹率與表觀氣速呈正相關(guān)，當表觀氣速為0.1，0.2，0.4 m/s時，De Felice模型比Gidaspow模型的床層膨脹率預(yù)測值準確度分別提高1%，4%，8.2%，表明De Felice模型對床層膨脹率的預(yù)測準確性隨表觀氣速增加而提高，Wen-Yu模型的預(yù)測結(jié)果偏低。在顆粒低濃度范圍，De Felice模型的曳力系數(shù)最小，氣體流動越順暢，容易聚集產(chǎn)生氣泡，且有變大趨勢，增強了卷吸夾帶作用，造成床層體積增大，床層膨脹率增加。反之，Wen-Yu模型的曳力系數(shù)較低，氣泡較少，尤其是大氣泡較難形成，預(yù)測的床層膨脹率明顯較低。

圖3 不同曳力模型下床層膨脹率隨表觀氣速的變化Fig.3 Variation of bed expansion ratio with apparent gas velocity under different drag models.

圖4為依次采用De Felice，Gidaspow，Wen-Yu曳力模型獲得的鼓泡床縱截面固體顆粒速度分布矢量。由圖4可知，De Felice和Gidaspow模型模擬計算所得床層膨脹量幾乎相同，大于Wen-Yu模型計算結(jié)果。De Felice模型計算得到顆粒分布較為均勻，氣泡處于中上部，下部較少；Gidaspow，Wen-Yu模型顆粒較多地聚集在頂部，阻礙氣體流出，氣泡基本處于中部及下部。采用De Felice，Gidaspow，Wen-Yu模型計算床層膨脹量結(jié)果分別為97.5%，92.6%，62.3%，而實驗測得結(jié)果為86.7%，可見Wen-Yu模型床層膨脹量偏小，與實驗結(jié)果不符。

圖4 縱剖面不同曳力模型的速度矢量Fig.4 Speed vector diagram of different drag models in longitudinal section.

床層壓降是評判曳力模型的重要標準。圖5為壓降隨表觀氣速的變化。由圖5a可知，3種曳力模型總壓降預(yù)測值基本保持不變，表觀氣速增加的對其影響并不明顯，而實驗值卻與表觀氣速呈正相關(guān)趨勢?，F(xiàn)有研究顯示氣固兩相流態(tài)化狀態(tài)達到穩(wěn)態(tài)時，總壓降不隨表觀氣速大小的改變而改變。但是壁面附近通常會有顆粒堆積現(xiàn)象，在高表觀氣速條件下測得實際壓降高于理論計算壓降。當表觀氣速增大時，鼓泡床內(nèi)氣泡更易破裂，部分顆粒會被破裂氣泡產(chǎn)生的氣流瞬間帶出鼓泡床，造成大氣速下壓力監(jiān)測線測量的壓力值小于模擬預(yù)測值的現(xiàn)象。由圖5b可知，3種曳力模型下鼓泡床內(nèi)的局部壓降預(yù)測值與實驗值，二者變化趨勢基本相同，均隨表觀氣速增加而降低。但預(yù)測值明顯與實驗值偏差較大，主要是Gidaspow和Wen-Yu模型，而De Felice模型的預(yù)測效果與實際比較接近。表觀氣速增大，導致床層膨脹變大，大部分顆粒被氣流帶到鼓泡床床層中上部，壓力檢測線測得壓力變小，表明床層局部壓降降低；De Felice模型預(yù)測的局部壓降小于其他兩種模型，因為它對床層膨脹的預(yù)測值最大。

圖5 壓降隨表觀氣速的變化Fig.5 Variation of pressure drop with apparent gas velocity.

圖6為不同曳力模型條件下固相速率隨表觀氣速的徑向變化曲線。

圖6 不同曳力模型條件下固相速率隨表觀氣速的徑向變化曲線Fig.6 Radial variation curves of solid phase velocity with apparent gas velocity under different drag model conditions.

由圖6可知，3種模型預(yù)測值較實驗值偏高，而De Felice模型的預(yù)測值與實驗值接近但略高于實驗值。實驗曲線分別在x/R=0，x/R=±0.5（x:顆粒距離鼓泡床軸線的水平距離；R:鼓泡床內(nèi)徑，0.5 m）附近存在一個波峰和兩個波谷，與Gidaspow模型模擬的曲線相似，變化趨勢基本一致，而Wen-Yu曲線與De Felice曲線相似，但是兩者皆與實驗有較大偏差。因此，De Felice模型的模擬結(jié)果明顯好于Wen-Yu模型，且更加接近實驗結(jié)果，而Gidaspow模型的預(yù)測結(jié)果與實際固含率曲線偏差最大。

2.2 不同鏡面系數(shù)的影響

采用De Felice曳力模型，獲得了不同ψ（0.01，0.03，0.05，0.10，0.20）下鼓泡床內(nèi)固相濃度梯度分布特性。圖7為在z=0.10 m截面處不同ψ下固含率大小對顆粒溫度影響關(guān)系曲線。

圖7 顆粒溫度隨固含率的變化Fig.7 Variation of particle temperature with solid content.Section:z=0.10 m.

由圖7可知，當固含率大于1%時，顆粒溫度變化隨著固含率的增大而變小，直至為零；而當固含率小于1%時，顆粒溫度變化幅度較大，大約在固含率為0.6%時有一極大值點；顆粒溫度與ψ的變化呈正相關(guān)。當固含率較小時，此時鼓泡床內(nèi)幾乎沒有顆粒團聚現(xiàn)象，單個顆粒相較于顆粒團與壁面碰撞或者顆粒間碰撞導致的動能耗散大小不同，沒有碰撞現(xiàn)象的顆粒高于有碰撞現(xiàn)象的顆粒溫度；當固含率逐漸變大時，顆粒更多的以顆粒團的形式運動，顆粒濃度的增加導致顆粒團與壁面和顆粒團之間的碰撞頻率增加，導致顆粒溫度下降。根據(jù)顆粒動力學理論，顆粒碰撞會產(chǎn)生黏度，與壁面碰撞后顆粒黏度會變小，而顆粒團聚現(xiàn)象的產(chǎn)生與顆粒濃度變大有關(guān)。

圖8為不同ψ下的縱截面固含率分布云圖。由圖8可知，當ψ=0.01時，縱截面固相很少附著在壁面附近，顆粒與壁面幾乎沒有碰撞與動量傳遞，產(chǎn)生了較大氣泡，總體濃度梯度變化不明顯，流化效果較好，ψ增大導致顆粒與壁面產(chǎn)生的碰撞逐漸加劇，當ψ＞0.05時顆粒的分布愈加不均勻，濃度梯度明顯，更多固體堆積壁面附近，氣泡外形也變得“尖銳”，與實驗現(xiàn)象相同，表明動能耗散加劇，擴散效果變差。

圖8 縱截面固含率分布云圖Fig.8 Longitudinal section solid content cloud map.t=8 s.

圖9為在z=0.05 m截面處顆粒徑向分布曲線。由圖9可知，除ψ=0.01曲線外，其他曲線走勢基本相同，中間區(qū)域變化趨勢明顯，在x/R=±0.5處存在2個波峰，而在x/R=±0處存在波谷。當ψ逐漸增大時，壁面附近固體團聚現(xiàn)象顯著，固含率較高，曲線波動范圍變大，計算結(jié)果略有偏大。總體來看ψ=0.10時預(yù)測所得結(jié)果與實驗結(jié)果基本吻合，而ψ=0.01的曲線較平緩，與實驗偏差最大。

綜上所述，鼓泡床內(nèi)氣固兩相流中的固體顆粒容易團聚，以顆粒團的形式存在，運動形式隨時處于變化。因此為使預(yù)測值更加貼近實驗值，ψ大小與顆粒團大小、聚集破裂的關(guān)系仍需進行研究。

圖9 顆粒徑向分布曲線Fig.9 Radial distribution curves of particles.

3 結(jié)論

1）當表觀氣速為0.1，0.2，0.4 m/s時，De Felice模型比Gidaspow模型的床層膨脹率預(yù)測值準確度分別提高1%，4%，8.2%，表明De Felice模型對床層膨脹率的預(yù)測準確性隨表觀氣速增加而提高。

2）Wen-Yu模型得到的床層膨脹量與實驗結(jié)果誤差較大，但流場分布與De Felice模型相似，Gidaspow，De Felice模型床層膨脹量接近，而氣固流場不同，考慮對總壓降和局部壓降的預(yù)測精度，De Felice模型預(yù)測結(jié)果與實驗值最接近。

3）不同ψ下顆粒徑向分布曲線（除ψ=0.01外）的波峰大致分別位于x/R=±0.5和x/R=±0處，隨著ψ增大，壁面附近固含率較高，曲線波動梯度增大。比較預(yù)測值與實驗值，當ψ=0.10時鼓泡床內(nèi)固含率分布結(jié)果與實驗值最接近，表明顆粒與壁面動量互換效果基本符合實際情況。

參 考 文 獻

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

［2］ Jia Dening，Bi Xiaotao，Jim L C，et al. Heat transfer in a pulsed fluidized bed of biomass particles［J］. Ind Chem Res，2017，56:3740-3756.

［3］ Wang Xiao，Si Hui.Conveying characteristics of dual pneumatic feeder used for biomass pyrolysis［J］. Bioresour J Citation Rep，2017，12（3）:5970-5983.

［4］ Mittal A，Mallick S S，Wypych P W. An investigation into flow mode transition and pressure fluctuations for fluidized dense-phase pneumatic conveying of fine powders［J］. Particuology，2014，7（5）:187-195.

［6］ Papari S，Hawboldt K. A review on the pyrolysis of woody biomass to bio-oil:Focus on kinetic models［J］. Renew Sustain Energy Rev，2015，19（2）:1580-1595.

［6］ 宋寧寧，唐曉津，侯栓弟. 鼓泡床流體力學特性研究［J］. 石油煉制與化工，2017，48（2）:89-92.

［7］ 劉學冰. 雙流化床鍋爐和雙流態(tài)高效潔凈燃燒工藝［J］. 潔凈煤技術(shù)，2016，22（3）:84-89.

［8］ 李競岌，楊欣華，楊海瑞，等. 鼓泡床焦炭型氮氧化物生成的試驗與模型研究［J］. 煤炭學報，2016，41（6）:1546-1553.

［9］ Sasithorn S，Benjapon C，Pornpote P. Arti ficial neural network model for the prediction of kinetic parametersof biomass pyrolysis from its constituents［J］. Fuel，2017，193（4）:142-158.

［10］ 張佳寶，崔麗杰，楊寧. 曳力模型和湍流模型對內(nèi)環(huán)流反應(yīng)器數(shù)值模擬的影響［J］. 化工學報，2018，69（1）:389-395.

［11］ 鄭曉野. 鼓泡流化床曳力模型及氣固流動特性數(shù)值研究［D］. 南京:南京航空航天大學，2016.

［12］ Johnson P C，Jackson R. Frictional-collisional constitutive relations for granular materials with application to plane slearing［J］. J Fluid Mech，2006，176:67-93.

［13］ 袁竹林，朱立平，耿凡，等. 氣固兩相流與數(shù)值模擬［M］.南京:東南大學出版社，2013:160-168.

［14］ Fronk B M，Garimella S. In-tube condensation of zeotropic fluid mixtures:A review［J］. Int J Ref，2013，36（2）:534-561.

［15］ Allen E E，Zhu C P，James S，et al. Multicomponent condensation reactions via ortho-quinone methides［J］. Org Lett，2017，19（7）:1878-1881.