循環(huán)流化床內(nèi)顆?；旌咸匦缘臄?shù)值模擬

李 斌，宋小龍

（華北電力大學(xué) 能源動力與機(jī)械工程學(xué)院，保定071003）

循環(huán)流化床以其高效、低污染、燃料適用廣等優(yōu)點在電力行業(yè)得到了較廣泛的應(yīng)用，但其床內(nèi)復(fù)雜的氣固兩相流還難以被人們充分了解.由于床內(nèi)顆粒的運動和混合情況會影響顆粒與氣體的傳熱、傳質(zhì)、床內(nèi)溫度分布、揮發(fā)分的析出和顆粒的燃燒［1］，因此深入研究床內(nèi)顆粒的運動和混合情況對循環(huán)流化床的優(yōu)化設(shè)計和運行具有重要意義.

目前，不少研究者對循環(huán)流化床內(nèi)的混合情況進(jìn)行了實驗研究［2-7］，但是由于床內(nèi)氣固兩相流流動的復(fù)雜性及實驗條件的制約，實驗方法在床內(nèi)氣固兩相流的研究中具有一定的局限性，且在宏觀現(xiàn)象的研究中難以反映過程的本質(zhì).數(shù)值模擬作為對理論分析與實驗研究互補(bǔ)的一種重要研究方法，能夠獲得顆粒豐富的信息.隨著計算機(jī)的發(fā)展，數(shù)值模擬在床內(nèi)氣固兩相流的研究中起到了十分重要的作用.因此，也有學(xué)者對循環(huán)流化床內(nèi)顆粒的混合過程進(jìn)行了數(shù)值模擬［8-12］，但是從氣體和顆粒運動速度的角度將軸向混合和徑向混合結(jié)合起來，并對顆粒整體混合過程進(jìn)行研究的文獻(xiàn)還比較少.

數(shù)值模擬的關(guān)鍵是建立正確數(shù)學(xué)模型，氣固兩相流數(shù)值模擬模型主要分為2類：一類是基于歐拉方法的顆粒相擬流體模型［13-14］；另一類是基于拉格朗日方法的顆粒軌道模型［15-16］.筆者將計算流體力學(xué)與屬于顆粒軌道模型范疇的離散單元法相結(jié)合，在已有工作［17］的基礎(chǔ)上應(yīng)用自行開發(fā)的程序?qū)Υ矁?nèi)顆粒的軸向混合和徑向混合過程進(jìn)行數(shù)值模擬，并采用顆粒特征濃度標(biāo)準(zhǔn)差對床內(nèi)顆粒的混合程度進(jìn)行定量評價.從氣體和顆粒運動速度的角度將軸向混合和徑向混合結(jié)合起來，對顆粒整體混合過程進(jìn)行分析，在微觀層次上揭示床內(nèi)顆?；旌蠙C(jī)理.

1 數(shù)值模擬方法

1.1 數(shù)學(xué)模型簡述

將計算流體力學(xué)與屬于歐拉-拉格朗日范疇的離散單元法結(jié)合起來，在歐拉坐標(biāo)系下考察氣相的運動，在拉格朗日坐標(biāo)系下考察顆粒的運動，然后再將兩者進(jìn)行耦合.即氣相采用考慮氣固兩相流耦合的Navier-Stocks方程，氣相湍流模型采用k-ε兩方程模型，并采用Simpler算法進(jìn)行求解；顆粒間碰撞采用軟球模型，根據(jù)牛頓第二定律建立每一個顆粒的運動方程，對每一個顆粒的運動情況進(jìn)行求解；氣相和顆粒相的耦合遵循牛頓第三定律，詳細(xì)求解過程參見文獻(xiàn)［17］和文獻(xiàn)［18］.

1.2 數(shù)值模擬對象

所模擬的對象為150mm×4mm×900mm 矩形截面準(zhǔn)三維流化床（床深度為顆粒直徑即4 mm），在床身底部中心位置設(shè)置一個空氣進(jìn)口，空氣進(jìn)口寬度為10mm，網(wǎng)格劃分為17×92.空氣密度為1.205 kg／m3，動力黏度為1.8×10-5kg／（m·s）.顆粒采用2 400個直徑均為4mm 的小球，其中小球的密度、恢復(fù)系數(shù)r、彈性系數(shù)k和摩擦因數(shù)分別為2 700kg／m3、0.9、800N／m 和0.3，黏性系數(shù)采用（m為等效質(zhì)量）計算.

1.3 數(shù)值模擬過程

首先在床內(nèi)隨機(jī)生成2 400個顆粒，將其在只考慮重力和顆粒間碰撞力的作用下進(jìn)行自由落體運動，顆粒最終的靜止位置作為模擬過程中顆粒的初始位置.為了清晰地模擬床內(nèi)顆粒的混合過程，在模擬過程中將初始位置的顆粒沿軸向和徑向平均分成2組，并標(biāo)記不同的顏色，在開始混合之前2種顆粒處于分離狀態(tài).

2 數(shù)值模擬結(jié)果

2.1 混合過程序列圖

當(dāng)空氣進(jìn)口速度為29m／s時，模擬了床內(nèi)顆粒的混合過程.圖1和圖2分別為0～2s內(nèi)顆粒軸向混合和徑向混合的序列圖.由圖1可知，隨著流化氣體的加入，中間底部的顆粒在氣流卷吸作用下向上運動，同時兩側(cè)靠近底部的顆粒向下運動，并向被流化氣體帶走的顆粒的原始位置進(jìn)行補(bǔ)充，進(jìn)入空氣進(jìn)口的噴射區(qū)，當(dāng)向上運動的顆粒到達(dá)床層表面時，顆粒向兩側(cè)進(jìn)行擴(kuò)散.由圖2可知，到達(dá)床層表面的顆粒分別向?qū)?cè)進(jìn)行擴(kuò)散，顆粒較大范圍的混合首先在床層表面進(jìn)行.由于壁面附近的氣體速度較小，擴(kuò)散到壁面附近的顆粒重力大于氣體對顆粒的曳力，顆粒開始向下運動，在向下運動過程中不斷與壁面進(jìn)行碰撞，同時與周圍的顆粒進(jìn)行混合.當(dāng)顆粒運動到床層底部區(qū)域時，重新被卷吸到噴射區(qū)，在氣流的作用下向上運動，形成床內(nèi)顆粒的內(nèi)循環(huán)，這樣周而復(fù)始，實現(xiàn)了顆粒的整體混合.從圖2還可以看出，由顆粒直接徑向運動而導(dǎo)致的徑向混合比較微弱，由噴射氣流卷吸顆粒形成床內(nèi)顆粒的內(nèi)循環(huán)成為影響顆?；旌系闹饕蛩?

2.2 氣體速度分布

圖1 0～2s內(nèi)顆粒軸向混合序列圖Fig.1 Sequence diagram of axial mixing within two seconds

圖2 0～2s內(nèi)顆粒徑向混合序列圖Fig.2 Sequence diagram of radial mixing within two seconds

圖3和圖4給出了氣體水平方向和垂直方向速度沿床層寬度（簡稱“床寬”）的分布，其中h為床層高度（簡稱“床高”）.從圖3可以看出，氣體水平方向速度在噴口中心兩側(cè)的方向相反，速度大小以噴口為中心基本呈對稱分布.在靠近床層底部區(qū)域，由于氣體受到主氣流的卷吸作用，噴口中心兩側(cè)的氣體向中心聚集，其他區(qū)域的氣體由中心向兩側(cè)擴(kuò)散，且隨著床高的增加，氣體水平方向速度減小.從圖4可以看出，氣體垂直方向速度呈現(xiàn)中心高、兩側(cè)低的趨勢，且隨著床高的增加而減小.對比圖3和圖4可以看出，氣體水平方向速度明顯小于垂直方向速度，顆粒在床內(nèi)主要受到氣體對顆粒的曳力、自身重力和碰撞力的作用，這就影響了顆粒在水平方向的運動，進(jìn)而影響了顆粒的徑向混合.

圖3 氣體水平方向速度分布Fig.3 Horizontal velocity distribution of gas

圖4 氣體垂直方向速度分布Fig.4 Vertical velocity distribution of gas

2.3 顆粒速度分布

圖5和圖6給出了0～2s內(nèi)床內(nèi)不同方向顆粒速度的平均值.由圖5和圖6可知，顆粒速度的運動趨勢與氣體基本一致.靠近床層底部的顆粒由于受到主氣流的卷吸作用以及對噴口中心區(qū)域進(jìn)行補(bǔ)充，噴口中心左側(cè)的顆粒向右運動，同時右側(cè)的顆粒向左運動，受以上雙重作用的影響，該區(qū)域顆粒水平方向速度大于其他區(qū)域.此外，在靠近床層表面，由于受中間氣體向兩側(cè)的推擠作用和受周圍顆粒的影響較小，靠近床層表面的顆粒水平方向速度也較大，床層表面的混合要好一些.在垂直方向上，噴口中心處的顆粒向上運動且速度較大，壁面處的顆粒速度為負(fù)，即向下運動.同樣，顆粒水平方向速度明顯小于顆粒垂直方向速度，從而導(dǎo)致顆粒的徑向混合明顯弱于軸向混合.

圖5 顆粒水平方向速度分布Fig.5 Horizontal velocity distribution of particles

圖6 顆粒垂直方向速度分布Fig.6 Vertical velocity distribution of particles

圖7為1.0s時顆粒的速度矢量圖.從圖7可以看出，顆粒在床內(nèi)的運動符合由顆?；旌闲蛄袌D及顆粒速度場得出的顆粒在床內(nèi)運動情況的結(jié)論，也證明了以上分析是正確的.

圖7 1.0s時顆粒的速度矢量Fig.7 Vector diagram of particle velocity at 1.0s

2.4 顆粒的體積分?jǐn)?shù)分布

圖8 1.0s時顆粒的體積分?jǐn)?shù)分布Fig.8 Distribution of particle volume fraction at 1.0s

圖8為1.0s時顆粒的體積分?jǐn)?shù)分布.由圖8可知，床內(nèi)底部區(qū)域顆粒的體積分?jǐn)?shù)大于床層上部，即下部為濃相區(qū)，上部為稀相區(qū)，且靠近壁面兩側(cè)的顆粒體積分?jǐn)?shù)大于中心區(qū)域，呈現(xiàn)中心稀、兩側(cè)濃的特點.

3 顆?；旌铣潭鹊亩吭u價

隨著時間的推移，顆粒的混合程度不斷深入，從直觀上難以對床內(nèi)顆粒的混合程度進(jìn)行有效、準(zhǔn)確的評價.為了更好地定量考察顆粒的混合程度，采用當(dāng)?shù)仡w粒特征濃度對顆粒的混合情況進(jìn)行分析［5］.將床體平均分成若干個參考單元，選定某一組分的顆粒為基準(zhǔn)，分別計算每個參考單元內(nèi)該組分顆粒的濃度（即該顆粒的特征濃度）.顆粒特征濃度e為

式中：Va、Vb分別為床內(nèi)2 種組分顆粒a和顆粒b在參考單元內(nèi)的體積，其中以顆粒a為基準(zhǔn).

筆者選定標(biāo)記為黑色的顆粒作為基準(zhǔn).e的大小反映了2 種組分顆粒在床內(nèi)局部區(qū)域的混合情況.當(dāng)e＝0或e＝1時，表示參考單元內(nèi)只有1種組分顆粒；當(dāng)e＝0.5時，表示參考單元內(nèi)2種組分的顆粒各占一半，混合效果相對較好；當(dāng)0＜e＜0.5或0.5＜e＜1時，表示顆粒處于混合過程中.e越接近0.5，表示混合效果越好.

為了考察床內(nèi)顆粒偏離完全混合的程度，將參考單元內(nèi)的e作為取樣樣本，以顆粒a的平均特征濃度作為參考，計算參考單元內(nèi)e與的偏差，記為e的標(biāo)準(zhǔn)差σ（e）.

式中：E為數(shù)學(xué)期望；由于只有2種組分的顆粒且數(shù)目和直徑均相等，因此選定為0.5.

將床體平均分成5×30個參考單元，計算顆粒所在區(qū)域的σ（e），對0～2s內(nèi)顆粒的混合過程進(jìn)行了定量分析，并討論了不同流化風(fēng)速對顆粒混合的影響.

圖9給出了空氣進(jìn)口速度v＝29m／s時σ（e）隨時間的變化趨勢.由圖9可知，隨著時間的推移，床內(nèi)顆粒軸向混合和徑向混合過程的σ（e）不斷減小，表明顆粒的混合過程不斷深入.在混合過程中，σ（e）會出現(xiàn)一定的波動，但隨著時間的推移，σ（e）的變化總趨勢是逐漸減小的，與文獻(xiàn)［8］中采用混合指數(shù)得到的結(jié)果基本一致.這是因為在顆?；旌线^程中總會出現(xiàn)混合與偏析2種過程，當(dāng)局部區(qū)域的混合惡化時，就會影響全場的σ（e），尤其是在顆?；旌陷^完全的時候.從圖9還可以看出，在同一時刻，顆粒軸向混合的σ（e）總是小于顆粒徑向混合的σ（e），表明床內(nèi)顆粒的徑向混合弱于軸向混合，與之前的分析相吻合.

圖9 空氣進(jìn)口速度為29m／s時顆粒特征濃度標(biāo)準(zhǔn)差隨時間的變化Fig.9 Change of standard deviation of particle characteristic concentration with time at inlet air velocity of 29m／s

圖10 不同空氣進(jìn)口速度下顆粒特征濃度標(biāo)準(zhǔn)差隨時間的變化Fig.10 Change of standard deviation of particle characteristic concentration with time at different inlet air velocities

圖10給出了空氣進(jìn)口速度分別為27 m／s和29m／s（對應(yīng)的流化風(fēng)速分別為1.8 m／s和1.93 m／s）時σ（e）隨時間的變化.由圖10可以看出，較大流化風(fēng)速下顆粒軸向混合和徑向混合的σ（e）均小于較小流化風(fēng)速下，表明適當(dāng)?shù)卦龃罅骰L(fēng)速可以增大床內(nèi)顆粒的混合速度，使床內(nèi)顆粒達(dá)到完全混合的時間縮短；且增大流化風(fēng)速對顆粒徑向混合的影響要比對軸向混合的影響明顯，分析認(rèn)為這是由于增大流化風(fēng)速后，雖然顆粒軸向運動速度有所增大，但是相應(yīng)的床層膨脹高度也會增加，內(nèi)循環(huán)尺度變大，導(dǎo)致顆粒軸向混合對流化風(fēng)速的敏感性低于徑向混合.

4 結(jié) 論

（1）通過數(shù)值模擬得到了床內(nèi)顆粒軸向混合和徑向混合過程的序列圖，形象地展示了顆粒在床內(nèi)的混合過程，混合過程首先發(fā)生在床層表面.

（2）顆粒與氣體的速度分布趨勢基本一致.整體上中心區(qū)域的顆粒向上運動，壁面處的顆粒向下運動，床層底部顆粒向中心處聚集，其他區(qū)域的顆粒向兩側(cè)運動.顆粒水平方向速度明顯小于垂直方向速度.顆粒的體積分?jǐn)?shù)分布表現(xiàn)為下部為濃相區(qū)、上部為稀相區(qū)，沿壁面兩側(cè)的顆粒體積分?jǐn)?shù)大于中心區(qū)域的顆粒體積分?jǐn)?shù).

（3）采用顆粒特征濃度標(biāo)準(zhǔn)差對床內(nèi)顆粒的混合過程進(jìn)行了定量評價，顆粒的徑向混合弱于軸向混合.

（4）增大流化風(fēng)速有助于加速顆粒的混合過程和提高混合質(zhì)量，顆粒軸向混合對流化風(fēng)速的敏感性低于徑向混合.

［1］田鳳國，章明川，顧明言，等.流化床微觀混合過程的定量評價［J］.動力工程，2007，27（1）：45-49. TIAN Fengguo，ZHANG Mingchuan，GU Mingyan，et al.Quantitative evaluation of micro-mixing processes in fluidized beds［J］.Journal of Power Engineering，2007，27（1）：45-49.

［2］黃群星，馬增益，池涌，等.循環(huán)流化床內(nèi)局部顆?；旌咸匦缘难芯浚跩］.動力工程，2004，24（1）：13-17. HUANG Qunxing，MA Zengyi，CHI Yong，etal.Study of local particles＇mixture property in CFB［J］.Power Engineering，2004，24（1）：13-17.

［3］嚴(yán)建華，朱建新，馬增益，等.圖像法用于流化床顆?；旌咸匦缘难芯浚跩］.高校化學(xué)工程學(xué)報，2006，20（5）：745-751. YAN Jianhua，ZHU Jianxin，MA Zengyi，etal.Research on mixing character of fluidized bed particles with image processing technique［J］.Journal of Chemical Engineering of Chinese Universities，2006，20（5）：745-751.

［4］羅婭，馬增益，朱建新，等.圖像追蹤算法用于流化床顆粒橫向混合特性的研究［J］.熱力發(fā)電，2005，34（11）：41-45. LUO Ya，MA Zengyi，ZHU Jianxin，etal.Study on behavior of particles＇cross mixing by using picture tracking algorithm［J］.Thermal Power Generation，2005，34（11）：41-45.

［5］張勇，金保升，鐘文琪.噴動氣固流化床顆?；旌弦?guī)律的實驗研究［J］.中國電機(jī)工程學(xué)報，2008，28（20）：8-14. ZHANG Yong，JIN Baosheng，ZHONG Wenqi.Experimental investigation on particle mixing in spoutfluid bed［J］.Proceedings of the CSEE，2008，28（20）：8-14.

［6］BOKKERS G A，ANNALAND V S M，KUIPERS J.Mixing and segregation in a bidisperse gas-solid fluidised bed：a numerical and experimental study［J］.Powder Technology，2004，140（3）：176-186.

［7］BUSCIGLIO A，VELLA G，MICALE G，etal.Analysis of the bubbling behaviour of 2Dgas solid fluidized beds：Part I.digital image analysis technique［J］.Chemical Engineering Journal，2008，140（1）：398-413.

［8］田鳳國，章明川，齊永鋒，等.流化床軸徑向混合特性的數(shù)值研究［J］.中國電機(jī)工程學(xué)報，2006，26（21）：119-124. TIAN Fengguo，ZHANG Mingchuan，QI Yongfeng，et al.A numerical investigation on axial／lateral mixing in fluidized beds［J］.Proceedings of the CSEE，2006，26（21）：119-124.

［9］朱潤孺，朱衛(wèi)兵，邢力超，等.矩形噴動床混合特性的三維數(shù)值研究［J］.中國電機(jī)工程學(xué)報，2010，30（17）：12-16. ZHU Runru，ZHU Weibing，XING Lichao，etal.A three-dimensional numerical investigation on particle mixing characteristics in rectangular spouted beds［J］.Proceedings of the CSEE，2010，30（17）：12-16.

［10］趙永志，鄭津洋.寬粒徑分布流化床的微觀尺度模擬與分析［J］.中國電機(jī)工程學(xué)報，2007，27（35）：55-61. ZHAO Yongzhi，ZHENG Jinyang.Micro-scale simulation and analysis of gas-solid fluidized bed with wide size distribution of particles［J］.Proceedings of the CSEE，2007，27（35）：55-61.

［11］RHODES M J，WANG X S，NGUYEN M，etal.Study of mixing in gas-fluidized beds using a DEM model［J］.Chemical Engineering Science，2001，56（8）：2859-2866.

［12］FENG Y Q，YU A B.Microdynamic modelling and analysis of the mixing and segregation of binary mixtures of particles in gas fluidization［J］.Chemical Engineering Science，2007，62（1）：256-268.

［13］郝振華，陳巨輝，白穎華，等.考慮顆粒旋轉(zhuǎn)的顆粒動力學(xué)模擬提升管氣固兩相流動特性［J］.高?；瘜W(xué)工程學(xué)報，2010，24（5）：776-782. HAO Zhenhua，CHEN Juhui，BAI Yinghua，etal.Numerical study on gas-solid flow in a riser using dynamic theory with consideration of particle rotation［J］.Journal of Chemical Engineering of Chinese Universities，2010，24（5）：776-782.

［14］王淑彥，趙云華，姜健，等.數(shù)值模擬顆粒旋轉(zhuǎn)對流化床內(nèi)氣固兩相流動特性的影響［J］.工程熱物理學(xué)報，2007，28（2）：262-264. WANG Shuyan，ZHAO Yunhua，JIANG Jian，etal.Numerical simulations of effect of particle rotation on gas and solid flow behavior in fluidized beds［J］.Journal of Engineering Thermophysics，2007，28（2）：262-264.

［15］XU B H，YU A B.Numerical simulation of the gassolid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics［J］.Chemical Engineering Science，1997，52（16）：2785-2809.

［16］任合斌.流化床內(nèi)流動、傳熱與燃燒特性的DEM 數(shù)值模擬［D］.武漢：華中科技大學(xué)，2006.

［17］李斌，紀(jì)律.流化床爐內(nèi)顆?；旌系碾x散單元法數(shù)值模擬［J］.中國電機(jī)工程學(xué)報，2012，32（20）：42-48. LI Bin，JI Lü.Numerical simulation of particle mixing in circulating fluidized bed with discrete element method［J］.Proceedings of the CSEE，2012，32（20）：42-48.

［18］紀(jì)律.循環(huán)流化床流動及磨損特性的DEM 數(shù)值模擬［D］.保定：華北電力大學(xué)，2011.