王 昊，駱仲泱，方夢(mèng)祥，王勤輝，岑建孟

（浙江大學(xué)能源清潔利用國(guó)家重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 310027）

與常規(guī)煤粉爐或?qū)尤紶t等相比，流化床爐氣化方式存在顯著差別。其根本原因是基本流動(dòng)、反應(yīng)原理上的差別。流化床是一種氣固兩相流態(tài)化原理的技術(shù)，爐內(nèi)存在大量的物料顆粒，這些顆粒不僅構(gòu)成十分復(fù)雜且粒徑分布十分廣泛，從幾微米到數(shù)十毫米不等，因此，形成了流化床獨(dú)特的空氣動(dòng)力特性。正是由于其不同的空氣動(dòng)力特性，使得流化床在傳熱、燃燒、污染物排放、灰渣特性等許多方面，具有不同的特點(diǎn)和規(guī)律[1-3]。由于試驗(yàn)手段本身的限制，靠實(shí)驗(yàn)手段難以完全掌握流化床內(nèi)的氣固流動(dòng)特性，尤其是流化床中顆粒的濃度以及運(yùn)動(dòng)分布規(guī)律（這些都是流化床的重要參數(shù)）。隨著計(jì)算技術(shù)的快速進(jìn)步以及人們對(duì)多相流顆粒模擬理論認(rèn)知的進(jìn)一步加深，使用計(jì)算流體力學(xué)（CFD）數(shù)值計(jì)算技術(shù)逐漸成為研究流化床內(nèi)氣固流動(dòng)特性的一種重要的輔助手段[4-11]。通過(guò)數(shù)值計(jì)算，可以快速獲得速度場(chǎng)、濃度場(chǎng)、溫度場(chǎng)的數(shù)據(jù)；對(duì)數(shù)據(jù)分析處理后可以研究流化床內(nèi)的氣固流動(dòng)規(guī)律。

圖1是由浙江大學(xué)提出的循環(huán)流化床多聯(lián)產(chǎn)系統(tǒng)[12-14]。該系統(tǒng)由蒸汽和再生煤氣鼓泡流化床熱解器及循環(huán)流化床氣化爐組成。其工藝流程為：煤在熱解爐內(nèi)首先與燃燒爐返料器過(guò)來(lái)的高溫循環(huán)熱載體灰混合發(fā)生裂解反應(yīng)，產(chǎn)生半焦和揮發(fā)分；析出的揮發(fā)分在旋風(fēng)分離器中與細(xì)灰顆粒分離；經(jīng)分離后的煤氣進(jìn)入煤氣凈化系統(tǒng)，依次經(jīng)過(guò)洗滌劑塔、電捕焦油器；收集的焦油可以出售或提取高附加值化工產(chǎn)品，或者經(jīng)催化加氫工藝制備高品位合成油；凈化后的煤氣部分通過(guò)再循環(huán)風(fēng)機(jī)加壓后送入熱解器底部作為流化風(fēng)，其余部分則進(jìn)入裝置進(jìn)一步凈化變成民用煤氣或者作為合成氣用于生產(chǎn)液體產(chǎn)品；熱解爐的半焦和低溫循環(huán)熱載體灰進(jìn)入燃燒爐；半焦燃燒的熱能一部分用于把循環(huán)灰加熱后再送回氣化爐以提供氣化爐所需熱量，另一部分熱能可以用于發(fā)電和供熱。

圖1 浙江大學(xué)1 MW循環(huán)流化床熱電氣多聯(lián)產(chǎn)實(shí)驗(yàn)裝置Fig.1 The 1 MW circulating fluidized bed heat-electricity-gas poly-generation test apparatus in Zhejiang University

MFiX軟件是由美國(guó)聯(lián)邦能源技術(shù)實(shí)驗(yàn)室（NETL）開(kāi)發(fā)的氣固兩相流開(kāi)源計(jì)算軟件，采用Fortran語(yǔ)言編寫(xiě)，用于流固系統(tǒng)的傳熱、化學(xué)反應(yīng)以及流體動(dòng)力學(xué)計(jì)算。該軟件包括雙流體模型TFM和DEM模型，主要解決流化床中的顆粒運(yùn)動(dòng)和燃燒問(wèn)題[15-23]。該軟件收錄了大量關(guān)于氣固兩相流研究的計(jì)算模型，同時(shí)又預(yù)留接口，方便研究人員對(duì)源程序進(jìn)行二次開(kāi)發(fā)。

Wang等人[24]研究了單固體顆粒相時(shí)二維、三維鼓泡流化床內(nèi)的氣固流動(dòng)特性，得到了流型、氣體速度、顆粒速度等數(shù)據(jù)，與實(shí)驗(yàn)結(jié)構(gòu)吻合較好。Wu等人[25]研究了不同曳力模型下的流化床燃燒過(guò)程，考察了不同的曳力模型對(duì)反應(yīng)速率的影響。Li等人[16]采用了開(kāi)源軟件MFiX對(duì)中試循環(huán)流化床系統(tǒng)進(jìn)行了模擬，得到數(shù)據(jù)與實(shí)驗(yàn)結(jié)果相吻合，但也發(fā)現(xiàn)了采用多固體顆粒相時(shí)會(huì)使計(jì)算量大大增加的情況。多固體顆粒相的多相流多顆粒相模擬求解困難是MFiX軟件存在的問(wèn)題。本文基于MFiX軟件，采用Fortran語(yǔ)言二次開(kāi)發(fā)將其拓展至多固體顆粒相系統(tǒng)模擬，再采用了OpenMPI并行技術(shù)，最終實(shí)現(xiàn)了流化床三維多顆粒相的熱解模擬。

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

1.1 連續(xù)性方程

式中：α為氣固相質(zhì)量分?jǐn)?shù)；ρ為相密度；u為速度；t為時(shí)間；下標(biāo)i、j為坐標(biāo)方向；下標(biāo)g和s分別表示氣相和顆粒相；下標(biāo)gs表示從氣相到顆粒相，下標(biāo)sg表示從顆粒相到氣相。

1.2 動(dòng)量方程

納維-斯托克斯（Navier-Stokes，N-S）方程是瞬時(shí)速度的動(dòng)量守恒方程，在數(shù)值計(jì)算的雙流體模型中，通過(guò)速度時(shí)均化，得到了N-S方程的時(shí)均動(dòng)量方程：

式中：Igs=-Isg，為氣體和固體顆粒相間的曳力；τ為應(yīng)力張量，主要由垂直于面的彈性力和平行于面的剪切力構(gòu)成。

1.3 能量方程

流動(dòng)過(guò)程中通常包含熱量的生成、傳遞和耗散，氣體和顆粒相的能量方程可用同一形式表達(dá)：

式中：m、n為第m或n相；Hm為物質(zhì)的焓，可以通過(guò)定壓比熱容計(jì)算化合物的焓值；λm為氣相和顆粒相的導(dǎo)熱系數(shù)，氣相的導(dǎo)熱系數(shù)可以通過(guò)不同組分導(dǎo)熱系統(tǒng)和含量推算而出；Qmn為氣體和不同顆粒相之間對(duì)流換的換熱量熱量；γ為輻射系數(shù)；Smn為化學(xué)反應(yīng)產(chǎn)生的源相。

1.4 組分輸運(yùn)方程

本文主要模擬的是煤的熱解和燃燒，它們都是多種化學(xué)反應(yīng)的集合。氣相內(nèi)部，氣相和顆粒相間都存在著組分的運(yùn)輸，均由多組分構(gòu)成。氣相包括CH4、H2、H2O、CO、H2O、CO2、N2和焦油氣8種，固相為固定碳、揮發(fā)分、水分和灰分4種。各組分滿足質(zhì)量守恒方程：

式中：X為組分的質(zhì)量分?jǐn)?shù)；R為組分輸運(yùn)方程的源相，表示均相或異相反應(yīng)的生成或消耗速度；J為組分?jǐn)U散流量。

1.5 顆粒溫度方程

式中：J為與碰撞有關(guān)的顆粒溫度的消耗；κs為顆粒溫度的擴(kuò)散系數(shù)；ΠΘ為氣固相之間的相互作用引起的顆粒相顆粒溫度的耗散。

1.6 化學(xué)反應(yīng)模型

化學(xué)反應(yīng)模型采用單元平均法來(lái)計(jì)算。具體方法為：通過(guò)使用網(wǎng)格單元中流體相的平均性質(zhì)來(lái)計(jì)算均相和非均相的反應(yīng)，同時(shí)忽略氣體擴(kuò)散速率對(duì)系統(tǒng)空間內(nèi)反應(yīng)的影響，即認(rèn)為反應(yīng)速率完全由反應(yīng)動(dòng)力學(xué)控制。

為簡(jiǎn)化計(jì)算，煤熱解歸納總結(jié)為3個(gè)異相反應(yīng)和5個(gè)均相反應(yīng)。異相反應(yīng)方程由煤的熱重曲線擬合得到，適用于500～1 200 ℃。氣體均相反應(yīng)采用文獻(xiàn)歸納方程[13]。

1）揮發(fā)分析出反應(yīng)：揮發(fā)分=0.235 8CO+0.109 5CO2+ 0.473 5CH4+ 0.806 9H2+1.229 3H2O+0.092 1Tar（Tar表示焦油）。

式中：R為顆粒相體積密度；Xs,V為顆粒相揮發(fā)分質(zhì)量分?jǐn)?shù)；VVM為中間變量。

2）半焦的氧化反應(yīng)：Cs+0.5O2=CO。

式中：kfm為膜層阻力系數(shù)；p(O2)為氧分壓；D(O2)為氧擴(kuò)散系數(shù)；Shm為舍伍德數(shù)；Sc為施密特?cái)?shù)；kam為灰層阻力系數(shù)；krm為表面反應(yīng)的阻力系數(shù)。

3）水分析出反應(yīng)：(H2O)s= H2O。

式中：Xs,m為顆粒相水分。

4）焦油裂解反應(yīng)：Tar=0.18CO+2.88CH4+7.74C。

式中：C為對(duì)應(yīng)物質(zhì)的摩爾濃度；εg為氣體體積分?jǐn)?shù)。

5）焦油氧化反應(yīng)：Tar+13.59O2=10.8CO2+5.76H2O。

7）H2燃燒反應(yīng)：H2+0.5O2=H2O。

2 物理模型

在三維笛卡爾坐標(biāo)系中對(duì)鼓泡流化床進(jìn)行模擬，網(wǎng)格采用結(jié)構(gòu)化網(wǎng)格，邊界網(wǎng)格使用實(shí)際幾何形狀切割結(jié)構(gòu)化網(wǎng)格，通過(guò)非正交項(xiàng)和內(nèi)部插值法計(jì)算出邊界處的速度節(jié)點(diǎn)位置。模擬網(wǎng)格如圖2所示。采用20倍煤顆粒直徑作為最大網(wǎng)格尺寸，逐漸減小進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證。網(wǎng)格尺寸為0.02 m×0.02 m×0.03 m時(shí)冷態(tài)軸向壓降曲線趨于一致，此時(shí)總網(wǎng)格數(shù)為35×35×250。同時(shí)，近壁面網(wǎng)格采取細(xì)化處理。

圖2 CFD模擬網(wǎng)格示意Fig.2 Schematic diagram of the CFD simulation mesh

恰當(dāng)?shù)倪吔鐥l件和初始條件的設(shè)定是獲得有效數(shù)值解的重要前提。本文模型中所有入口條件均設(shè)置為質(zhì)量入口條件，速度垂直于入口界面；出口均設(shè)置為壓力出口邊界條件，壓力為背景大氣壓力。氣體壁面邊界條件設(shè)定為無(wú)滑移邊界條件，即氣體壁面切向和法相的速度均為0；顆粒相壁面邊界條件設(shè)定為自由滑移邊界條件，即顆粒相壁面的法相速度為0，切向速度不受壁面力的影響，由動(dòng)量方程求得。模擬采用的關(guān)鍵模型見(jiàn)表1。

表1 模擬關(guān)鍵模型Tab.1 Key models of the simulation

鼓泡流化床內(nèi)初始床料高度為1.4 m，床料的體積分?jǐn)?shù)為0.5，最大堆積體積分?jǐn)?shù)為0.5，計(jì)算流體區(qū)域的氣相和顆粒相速度均為0。計(jì)算時(shí)間步長(zhǎng)為10-4s，每個(gè)步長(zhǎng)的迭代直到收斂為止，但不超過(guò)20次，不收斂時(shí)自動(dòng)減小時(shí)間步長(zhǎng)再次迭代直到收斂，再過(guò)渡到下一步計(jì)算時(shí)間。在模擬時(shí)間5 s后，對(duì)動(dòng)量和質(zhì)量偏微分方程采用二階差分方式計(jì)算，提高計(jì)算精度。模擬時(shí)間總共為30 s。

本文模擬了熱解爐實(shí)驗(yàn)裝置在500、600、700 ℃熱解溫度工況下的反應(yīng)氣固流動(dòng)特性，共模擬了3組工況，每組工況參數(shù)見(jiàn)表2。

表2 鼓泡流化床煤熱解模擬工況Tab.2 The simulation conditions for coal pyrolysis in bubbling fluidized bed

由于多聯(lián)產(chǎn)實(shí)驗(yàn)中熱解爐采用再循環(huán)熱解氣作為熱解爐的流化風(fēng)，本文中熱解爐的模擬采用簡(jiǎn)化方案。首先采用惰性氣體N2作為流化風(fēng)，待模擬至爐膛出口氣體組分分布穩(wěn)定，即熱解過(guò)程達(dá)到穩(wěn)態(tài)時(shí)，采用此時(shí)爐膛出口CH4、H2、CO、CO2（除N2、H2O、焦油）4種熱解氣成分進(jìn)行歸一化處理后作為新的流化風(fēng)成分。采用新的流化風(fēng)模擬至流動(dòng)穩(wěn)定的狀態(tài)來(lái)實(shí)現(xiàn)再循環(huán)熱解氣作為流化風(fēng)的模擬過(guò)程。煤的工業(yè)分析和元素分析見(jiàn)表3。

表3 工業(yè)分析和元素分析 w/%Tab.3 Industrial and elemental analysis results of the coal

3 模擬結(jié)果及討論

3.1 氣體組分分布

圖3—圖5給出了熱解溫度為500、600、700 ℃時(shí)熱解爐內(nèi)煤氣各組分分布。這是根據(jù)模擬計(jì)算所得到的氣體組分分布，排除焦油和水蒸氣組分后，對(duì)CO、H2、CO2、CH4采取體積歸一化處理的結(jié)果，主要用于和實(shí)驗(yàn)得到的數(shù)據(jù)做進(jìn)一步的對(duì)比。

圖3 500 ℃時(shí)均截面氣體組分分布Fig.3 The gas composition distribution on time averaged cross section at 500 ℃

圖4 600 ℃時(shí)均截面氣體組分分布Fig.4 The gas composition distribution on time averaged cross section at 600 ℃

圖5 700 ℃時(shí)均截面氣體組分分布Fig.5 The gas composition distribution on time averaged cross section at 700 ℃

從圖3—圖5可以看到，所有溫度段的熱解產(chǎn)物中的CO、H2、CO2在密相區(qū)的質(zhì)量分?jǐn)?shù)均高于稀相區(qū)，而CH4在密相區(qū)的質(zhì)量分?jǐn)?shù)低于稀相區(qū)。CH4不同于其他氣體組分分布特性，主要是模擬中考慮了煤焦油二次裂解所致。同時(shí)，也可以看出反應(yīng)器內(nèi)煤焦油的二次裂解主要發(fā)生在稀相區(qū)內(nèi)。煤焦油二次裂解的產(chǎn)物為CO、CH4、固定碳，其中CO的化學(xué)計(jì)量數(shù)較小，所以煤焦油二次裂解生成的CO較少，因而對(duì)稀相區(qū)CO的影響較少，所以CO和H2、CO2的組分分布規(guī)律具有一定的相似性。

從圖3—圖5還可以看到，CO、H2、CO2、CH4在密相區(qū)底部，尤其是靠近給料入口有一個(gè)明顯的高質(zhì)量分?jǐn)?shù)區(qū)域。在該區(qū)域內(nèi)，煤與床料混合均勻，反應(yīng)劇烈；同時(shí)在該區(qū)域內(nèi)也可以觀察到焦油和熱解水的生成，顯示出煤熱解過(guò)程的異相反應(yīng)均發(fā)生在環(huán)-隙結(jié)構(gòu)的底部。

隨著熱解溫度的升高，反應(yīng)速度加快，熱解爐內(nèi)氣體組分達(dá)到穩(wěn)態(tài)的速度隨之加快。500 ℃熱解時(shí)，氣體組分的穩(wěn)定區(qū)域在爐膛上部；600 ℃時(shí)，組分穩(wěn)定面上升到了爐膛中部的漸擴(kuò)段；而700 ℃時(shí)，所有的反應(yīng)基本在爐膛下部的密相區(qū)內(nèi)已經(jīng)完成，此時(shí)熱解產(chǎn)生的氣體可以部分充當(dāng)流化風(fēng)的作用。這也解釋了為何700 ℃時(shí)實(shí)驗(yàn)所需流化風(fēng)量有明顯的減少。

圖6—圖7為煤氣中CH4、H2體積分?jǐn)?shù)隨熱解溫度的變化。由圖6—圖7可以看到，本文建立的多聯(lián)產(chǎn)系統(tǒng)模型可以很好地預(yù)測(cè)煤氣中CH4的體積分?jǐn)?shù)且變化趨勢(shì)也與試驗(yàn)一致。在熱解溫度為600、700 ℃時(shí)，對(duì)于煤氣中H2體積分?jǐn)?shù)的模擬數(shù)據(jù)與試驗(yàn)數(shù)據(jù)相關(guān)性好；500 ℃時(shí)H2體積分?jǐn)?shù)的預(yù)測(cè)偏差較大，可能和揮發(fā)分反應(yīng)方程在低溫下的適用性有關(guān)。

圖6 熱解氣中CH4的體積分?jǐn)?shù)Fig.6 Volume fractions of CH4 in the pyrolysis gases

圖7 熱解氣中H2的體積分?jǐn)?shù)Fig.7 Volume fractions of H2 in the pyrolysis gases

將模擬獲得的熱解爐出口氣體組分除去水蒸氣后換算為標(biāo)準(zhǔn)狀況下的數(shù)據(jù)，得到熱解爐煤氣中焦油成分的質(zhì)量濃度。圖8為多聯(lián)產(chǎn)系統(tǒng)模擬熱解爐出口煤氣中焦油質(zhì)量濃度與多聯(lián)產(chǎn)試驗(yàn)的對(duì)比。由圖8可以看出，建立的熱解模型能較好地預(yù)測(cè)熱解氣中焦油的質(zhì)量濃度。

圖8 熱解氣中焦油的質(zhì)量濃度Fig.8 Mass cocentrations of tar in the pyrolysis gases

3.2 流化床截面顆粒體積分?jǐn)?shù)及速度分布

流化床不同截面（選定4個(gè)截面對(duì)應(yīng)高度分別為0.2、0.6、1.0、1.5 m）顆粒體積分?jǐn)?shù)及速度分布如圖9—圖10所示。

圖9 顆粒時(shí)均體積分?jǐn)?shù)Fig.9 The time averaged volume fraction of the particles

圖10 顆粒時(shí)均速度Fig.10 The time averaged speed of the particles

由圖9、圖10可以觀察到，熱解流化床爐不同高度處徑向截面上均出現(xiàn)了中心區(qū)域低顆粒體積分?jǐn)?shù)和近壁面區(qū)域高顆粒體積分?jǐn)?shù)的現(xiàn)象，即存在中心稀疏的環(huán)-核結(jié)構(gòu)。此外還可以發(fā)現(xiàn)，環(huán)-核結(jié)構(gòu)隨流化床高度的增加也有變化。在密相區(qū)0.2 m和0.4 m高度處，高體積分?jǐn)?shù)顆粒區(qū)域略有增大；0.6 m高度處，高體積分?jǐn)?shù)顆粒區(qū)域的面積大幅減小。而在漸擴(kuò)段上方1.5 m處，高體積分?jǐn)?shù)顆粒區(qū)域面積又有較大增加，直至0.2 m處顆粒相已經(jīng)極難觀測(cè)到。這主要是由于漸擴(kuò)段的存在，使一定質(zhì)量的床料長(zhǎng)時(shí)間地存在于漸擴(kuò)段及其上部區(qū)域形成內(nèi)部回流流動(dòng)，同時(shí)阻止了下部密相區(qū)沿邊界進(jìn)入漸擴(kuò)段的顆粒量。新的環(huán)-核結(jié)構(gòu)的形成也是由于相同的原因。因此，由于漸擴(kuò)段的存在流化床內(nèi)形成了以漸擴(kuò)段為界的相連的2個(gè)環(huán)-核結(jié)構(gòu)。

4 結(jié) 論

1）采用MFiX軟件，對(duì)浙江大學(xué)1 MW多聯(lián)產(chǎn)系統(tǒng)進(jìn)行簡(jiǎn)化數(shù)值模擬。使用熱解機(jī)理反應(yīng)模型和多固體顆粒相流體模擬方法建立了多聯(lián)產(chǎn)系統(tǒng)熱解爐熱解反應(yīng)三維模型，并對(duì)系統(tǒng)模擬參數(shù)進(jìn)行優(yōu)化。

2）對(duì)3個(gè)不同熱解溫度的工況進(jìn)行模擬，將計(jì)算的爐膛出口氣體組分的分布和焦油產(chǎn)率與多聯(lián)產(chǎn)實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。其中焦油產(chǎn)率和CH4的分布模擬結(jié)果與實(shí)驗(yàn)一致。H2的體積分?jǐn)?shù)的模擬數(shù)據(jù)和實(shí)驗(yàn)值在較低熱解溫度時(shí)有一定的誤差，在較高熱解溫度時(shí)趨于一致。

3）分析了不同熱解溫度下流化床內(nèi)的氣固流動(dòng)特性，以及固體顆粒的體積分?jǐn)?shù)和速度分布的特性，驗(yàn)證了流化床內(nèi)環(huán)-核結(jié)構(gòu)的存在。