循環(huán)流化床內(nèi)顆粒混合特性的數(shù)值模擬

李 斌，宋小龍

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

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

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

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

1 數(shù)值模擬方法

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

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

1.2 數(shù)值模擬對象

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

1.3 數(shù)值模擬過程

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

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

2.1 混合過程序列圖

當(dāng)空氣進(jìn)口速度為29m／s時(shí)，模擬了床內(nèi)顆粒的混合過程.圖1和圖2分別為0～2s內(nèi)顆粒軸向混合和徑向混合的序列圖.由圖1可知，隨著流化氣體的加入，中間底部的顆粒在氣流卷吸作用下向上運(yùn)動(dòng)，同時(shí)兩側(cè)靠近底部的顆粒向下運(yùn)動(dòng)，并向被流化氣體帶走的顆粒的原始位置進(jìn)行補(bǔ)充，進(jìn)入空氣進(jìn)口的噴射區(qū)，當(dāng)向上運(yùn)動(dòng)的顆粒到達(dá)床層表面時(shí)，顆粒向兩側(cè)進(jìn)行擴(kuò)散.由圖2可知，到達(dá)床層表面的顆粒分別向?qū)?cè)進(jìn)行擴(kuò)散，顆粒較大范圍的混合首先在床層表面進(jìn)行.由于壁面附近的氣體速度較小，擴(kuò)散到壁面附近的顆粒重力大于氣體對顆粒的曳力，顆粒開始向下運(yùn)動(dòng)，在向下運(yùn)動(dòng)過程中不斷與壁面進(jìn)行碰撞，同時(shí)與周圍的顆粒進(jìn)行混合.當(dāng)顆粒運(yùn)動(dòng)到床層底部區(qū)域時(shí)，重新被卷吸到噴射區(qū)，在氣流的作用下向上運(yùn)動(dòng)，形成床內(nèi)顆粒的內(nèi)循環(huán)，這樣周而復(fù)始，實(shí)現(xiàn)了顆粒的整體混合.從圖2還可以看出，由顆粒直接徑向運(yùn)動(dòng)而導(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)主要受到氣體對顆粒的曳力、自身重力和碰撞力的作用，這就影響了顆粒在水平方向的運(yùn)動(dòng)，進(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可知，顆粒速度的運(yùn)動(dòng)趨勢與氣體基本一致.靠近床層底部的顆粒由于受到主氣流的卷吸作用以及對噴口中心區(qū)域進(jìn)行補(bǔ)充，噴口中心左側(cè)的顆粒向右運(yùn)動(dòng)，同時(shí)右側(cè)的顆粒向左運(yùn)動(dòng)，受以上雙重作用的影響，該區(qū)域顆粒水平方向速度大于其他區(qū)域.此外，在靠近床層表面，由于受中間氣體向兩側(cè)的推擠作用和受周圍顆粒的影響較小，靠近床層表面的顆粒水平方向速度也較大，床層表面的混合要好一些.在垂直方向上，噴口中心處的顆粒向上運(yùn)動(dòng)且速度較大，壁面處的顆粒速度為負(fù)，即向下運(yùn)動(dòng).同樣，顆粒水平方向速度明顯小于顆粒垂直方向速度，從而導(dǎo)致顆粒的徑向混合明顯弱于軸向混合.

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

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

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

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

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

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

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

3 顆粒混合程度的定量評價(jià)

隨著時(shí)間的推移，顆粒的混合程度不斷深入，從直觀上難以對床內(nèi)顆粒的混合程度進(jìn)行有效、準(zhǔn)確的評價(jià).為了更好地定量考察顆粒的混合程度，采用當(dāng)?shù)仡w粒特征濃度對顆粒的混合情況進(jìn)行分析［5］.將床體平均分成若干個(gè)參考單元，選定某一組分的顆粒為基準(zhǔn)，分別計(jì)算每個(gè)參考單元內(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時(shí)，表示參考單元內(nèi)只有1種組分顆粒；當(dāng)e＝0.5時(shí)，表示參考單元內(nèi)2種組分的顆粒各占一半，混合效果相對較好；當(dāng)0＜e＜0.5或0.5＜e＜1時(shí)，表示顆粒處于混合過程中.e越接近0.5，表示混合效果越好.

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

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

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

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

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

4 結(jié) 論

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

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

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

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

［1］田鳳國，章明川，顧明言，等.流化床微觀混合過程的定量評價(jià)［J］.動(dòng)力工程，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)局部顆粒混合特性的研究［J］.動(dòng)力工程，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é)報(bào)，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］張勇，金保升，鐘文琪.噴動(dòng)氣固流化床顆?；旌弦?guī)律的實(shí)驗(yàn)研究［J］.中國電機(jī)工程學(xué)報(bào)，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é)報(bào)，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)兵，邢力超，等.矩形噴動(dòng)床混合特性的三維數(shù)值研究［J］.中國電機(jī)工程學(xué)報(bào)，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é)報(bào)，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)的顆粒動(dòng)力學(xué)模擬提升管氣固兩相流動(dòng)特性［J］.高校化學(xué)工程學(xué)報(bào)，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)氣固兩相流動(dòng)特性的影響［J］.工程熱物理學(xué)報(bào)，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)流動(dòng)、傳熱與燃燒特性的DEM 數(shù)值模擬［D］.武漢：華中科技大學(xué)，2006.

［17］李斌，紀(jì)律.流化床爐內(nèi)顆?；旌系碾x散單元法數(shù)值模擬［J］.中國電機(jī)工程學(xué)報(bào)，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)流化床流動(dòng)及磨損特性的DEM 數(shù)值模擬［D］.保定：華北電力大學(xué)，2011.