六回路循環流化床甲烷化過程數值模擬及優化

馬繼偉，張亞新,2

（1. 新疆大學 化工學院，新疆 烏魯木齊 830046；2. 新疆大學 省部共建碳基能源資源化學與利用國家重點實驗室，新疆 烏魯木齊 830046）

我國能源結構呈“富煤、貧油、少氣”的特點，因此煤炭在我國能源結構和國民經濟中的地位舉足輕重[1]。傳統煤炭在利用過程中易造成環境污染，煤制天然氣技術可實現煤炭的高效、清潔和安全利用，對實現2030 年碳達峰及2060 年碳中和目標具有積極推動意義。同時，煤制天然氣可有效彌補天然氣儲量小這一短板，對保障國家能源戰略安全及天然氣供需平衡具有重要意義[2]。

合成氣甲烷化為煤制天然氣重要工藝單元。甲烷化反應具有速率快和放熱強的特點[3]，如不能及時移熱，將導致床層溫度過高，使固體催化劑燒結，積炭嚴重，最終失去活性，并且局部高溫會對設備穩定運行構成威脅。甲烷化反應的強放熱特性及流化床反應器內復雜的氣固兩相流動，對甲烷化流化床反應器的設計及放大技術提出了更高的要求。流化床甲烷化技術一直落后于成熟的固定床技術，目前流化床甲烷化技術尚未實現工業化應用，但流化床反應器傳質、傳熱效率高且催化劑顆粒分布均勻，是強放熱甲烷化反應的理想反應器[4]，流化床甲烷化技術仍具有很好的發展前景。隨著計算機技術的發展，與實驗相比，數值模擬方法經濟靈活高效，可有效彌補實驗不足，已廣泛應用于流化床內反應過程研究中。目前，通過實驗及數值模擬方法針對實驗室規模甲烷化流化床反應器內反應特性進行了較多研究。LIU等[5]通過實驗探究了固定床與流化床甲烷化反應特性，結果表明，固定床存在明顯溫度梯度與熱點溫度，固定床出口轉化率明顯低于流化床。KOPYSCINSKI等[6-7]對Ni/Al2O3催化劑流化床甲烷化反應特性進行了實驗研究，依據實驗結果擬合得到了甲烷化反應動力學模型，并對流化床反應器進行了建模以研究流動及反應特性。SUN等[8]基于歐拉雙流體模型，模擬探究了不同工藝參數下實驗室規模流化床內氣固流動與甲烷化反應特性。張曉瑞等[9]針對小型流化床反應器進行了模擬研究，結果表明，甲烷化反應在入口段反應迅速，在10 mm處已基本反應完全。劉姣等[1]以小型甲烷化輸送床為研究對象，考察了操作參數對不同NiO含量（質量分數）的催化劑性能的影響。崔佃淼[10]針對實驗室小型循環流化床探究了操作參數對CO甲烷化反應性能的影響，結果表明，在高CO轉化率下床層溫度梯度也很小，高熱容的固體催化劑顆粒循環移除了大量熱量。但受實驗條件限制及反應器內復雜的氣固流動對反應器的設計及放大造成很大的困難，對于大型循環流化床甲烷化反應特性的研究較少，SUN等[11]針對底部直徑0.2 m，高1.2 m變徑循環流化床首次建立了全循環流化床甲烷化數值模型，分析了反應器內催化劑顆粒流化行為，并評價了操作參數對反應物轉化率及CH4產率的影響。車豪等[12]基于計算顆粒流體力學（CPFD）方法，對直徑0.3 m，高度8.5 m的循環流化床內甲烷化過程進行了模擬研究，獲得了反應器內部流場、溫度場以及組分濃度分布，并研究了不同操作條件對甲烷化反應特性的影響規律。

以上甲烷化技術只是針對單因素對甲烷化反應特性影響研究，未考慮操作參數間的交互作用，無法得到最優操作參數組合。此外，從以上甲烷化技術可以看出，循環流化床反應器內溫度梯度小，循環流化床甲烷化技術對提高CH4產率及反應器的穩定運行具有明顯優勢。對六回路循環流化床冷態研究中發現多回路的布置可有效提高氣固分離效率，可實現顆粒的快速流化讓更多的顆粒參加循環[13-14]，這對反應熱的移除是有利的，但目前尚未有多回路循環流化床甲烷化反應過程研究報道。

鑒于響應曲面法具有設計方法合理、實驗次數少、可對正交試驗結果進行統計與分析并通過構建多因素影響模型對操作參數進行預測優化的優點，本文基于CPFD方法，對六回路循環流化床甲烷化過程進行數值模擬，采用響應曲面法探究操作因素（入口氣速、入口溫度和n(H2)/n(CO)）對甲烷產率的影響，建立目標產物甲烷產率與操作參數的回歸模型，進一步得到最佳操作工藝參數。所開展工作可加強對六回路循環流化床內甲烷化反應特性的認識，將為甲烷技術的發展和反應器的設計提供參考。

1 模型及模擬

1.1 數學模型

1.1.1 控制方程

CPFD數值模擬方法通過在歐拉框架下基于連續介質假設求解流體運動，在拉格朗日框架下采用牛頓定律追蹤顆粒運動來反應離散顆粒場，并與連續性方程進行耦合實現兩相交互作用。主要數學方程如下[15-16]：

氣-固兩相流：

式中，εp為固體顆粒體積分數；εg為氣相體積分數。

氣相連續性方程：

式中，ρg為氣體密度，kg/m3；vg為氣體速度，m/s；Sg為氣體源項；t為時間，s。

氣相動量方程：

式中，p為氣體壓力，Pa；g為重力加速度，m/s2；F為氣相與固相動量交換率；τg為氣體應力張量；p為流體壓力梯度，Pa/m。

氣固動量交換率為：

式中，Dp為曳力系數；ρp為顆粒密度，kg/m3；mp為顆粒質量，kg；vp為顆粒速度，m/s；f為顆粒概率分布函數；Tp為顆粒溫度，K。

顆粒間的動量方程為：

式中，ps為壓力常數，Pa；ε為數量級為10-7的小量；εcp為顆粒緊密堆積體積分數；γ推薦值為2~5；τp為顆粒間作用力，Pa，采用顆粒正向作用力。

1.1.2 動力學方程

合成氣甲烷化過程主要反應為：

入口原料氣中CO2是微量的，為提高計算效率對模型進行簡化，僅對第一個反應進行模擬研究。

CO甲烷化反應動力學方程采用KOPYSCINSKI等[6]建立的模型，該模型是基于Ni/Al2O3催化劑通過實驗測量得到，CO甲烷化反應宏觀動力學方程為：

式中，Rm為CO甲烷化反應速率，mol/(kg·s)；K1為CO甲烷化反應速率常數，mol/((105Pa)0.5·kg·s)；KC和KOH為表面吸附有關的平衡常數，(105Pa)-0.5；pi為組分i的分壓，105Pa；T為溫度，K；R為氣體常數，取值8.314 J/(mol·K)。

甲烷產率（ ）定義為：

式中，wCH4,out為出口CH4質量分數；wCO,in為入口CO質量分數。

1.2 物理模型及參數設置

本文以六回路循環流化床反應器為模擬對象，反應器包括高度為5.80 m，長為0.92 m，寬為0.42 m的提升管和6 個H型軸對稱并聯布置的旋風分離器及返料器。因重點關注提升管內甲烷化反應特性，為節約計算成本，對模型進行簡化，采用軟件定義外部顆粒循環的方式代替旋風分離器與返料系統。即首先計算得到各出口顆粒質量流率，為保證反應器中催化劑物料量平衡，假設各分離器的分離效率為100%，將計算值通過對應的返料口由返料風全部送入提升管中，研究表明，上述對循環流化床中循環系統的簡化方法是可靠的[17-19]。初始時刻，催化劑顆粒堆積在爐膛底部，床層高度為309.00 mm，原料氣由提升管底部進入，原料氣中n(H2+ CO)/n(N2)為3.0，簡化的三維物理模型下初始催化劑顆粒堆積狀態及氣體進出口邊界如圖1 所示。模擬參數及邊界條件設置如表1 所示。催化劑顆粒粒徑分布如圖2 所示。

圖1 初始顆粒堆積狀態及氣體進出口邊界Fig. 1 Initial particle accumulation state and gas inlet and outlet boundary

圖2 催化劑顆粒粒徑分布Fig. 2 Catalyst particle size distribution

表1 模擬參數及邊界條件設置Table 1 Setting of simulation parameters and boundary conditions

1.3 網格無關性驗證

三維數值模型計算區域采用自適應笛卡爾網格進行劃分，網格質量在一定程度上影響數值計算精度，對簡化后的計算模型分別劃分網格數為77000、102000 和143000，不同網格數計算下CH4質量分數軸向分布如圖3 所示。由圖3 可知，102000與143000 網格數計算下，CH4質量分數沿軸向分布一致，說明當網格數為102000 時，網格數量對計算結果不會產生較大影響。在保證計算結果準確性下兼顧計算效率，最終選取102000 作為數值計算網格數。

圖3 網格無關性驗證Fig. 3 Grid independence verification

1.4 模擬方法有效性驗證

基于CPFD方法，本節以葉松壽[20]所建立CO甲烷化實驗裝置為研究對象，反應器直徑67.50 mm、高度1000.00 mm，實驗原料氣組成為n(H2+ CO)/n(N2)為3.0，其 中n(H2)/n(CO)為3.1，入 口 氣 體 溫 度593.0 K，初始催化劑溫度為623.0 K，出口壓力為101325 Pa。催化劑主要成分為Ni初始質量為125 g，其余氣固物性參數與實驗[20]保持一致。不同入口氣速下沿軸向CO轉化率模擬值與實驗值[20]對比如圖4 所示。由圖4 可知，模擬值與實驗值吻合較好，平均誤差均小于3.00%，驗證了數值模擬方法的準確性。雖然六回路循環流化床反應器與本節研究對象提升管形狀存在差異，前者為方箱形，后者為圓管形，但后續模擬研究中甲烷化動力學模型及模擬參數設置均與本節一致，因此本文針對六回路循環流化床甲烷化反應特性的研究是可靠的。

圖4 模擬值與實驗值對比Fig. 4 Comparison between simulated and experimental values

2 數值模擬結果及分析

甲烷化過程易受操作參數的影響。入口氣速直接影響流化床反應器內顆粒流化狀態，也是反應器處理能力的體現。強放熱的甲烷化反應是一個可逆反應，合適的入口溫度對保證反應平衡及催化劑活性至關重要。煤氣化產物中n(H2)/n(CO)通常為0.3~2.0，對實現CO的完全轉化這一比值較低[11]，為提高CO轉化率，對反應物中氣體組成有必要進行研究。為此，探究了操作參數（入口氣速、入口溫度和n(H2)/n(CO)）對反應器內甲烷化過程的影響，在對不同操作參數下CH4質量分數及CH4產率分布規律分析基礎上，通過響應曲面法對操作參數進行優化，可進一步認識循環流化床反應器內甲烷化反應特性。

2.1 入口氣速對甲烷化過程的影響

在n(H2)/n(CO)為2.5，入口溫度為540.0 K下，入口氣速從5.00 m/s增大至13.00 m/s時CH4質量分數沿反應器軸向分布如圖5 所示。由圖5 可知，CH4質量分數沿軸向逐漸增大，較低氣速下在一定高度處反應基本完全后CH4質量分數波動較小。隨著入口氣速的增大，需要更高的距離來完成反應，當入口氣速為5.00 m/s和9.00 m/s時，完成反應所需要的高度分別為0.70 m和3.25 m。相對高度為2.25 m，入口氣速為13.00 m/s時，CH4質量分數軸向分布明顯較小，這是因為隨著入口氣速的增大更多的催化劑顆粒被帶到提升管上部區域參加循環，提升管中催化劑顆粒體積分數的減小，減小了催化劑與氣體接觸面積且縮短了催化劑顆粒及反應物停留時間。在同一高度處，CH4質量分數隨著氣速的增大而減小，氣速過大不利于反應的進行。在底部由于受到催化劑顆粒返料影響，催化劑顆粒體積分數的波動使得提升管底部的CH4質量分數也出現較大波動。CH4產率隨入口氣速的變化情況如圖6 所示，可明顯看出，當入口氣速小于9.00 m/s時能獲得較高的CH4產率，且當入口氣速為7.50 m/s時，CH4產率最大為0.4566。當入口氣速大于9.00 m/s時，CH4產率出現明顯下降，因此，入口氣速取5.00~9.00 m/s。

圖5 不同入口氣速下CH4 質量分數軸向分布Fig. 5 Axial distribution of CH4 mass fraction at different inlet gas velocities

圖6 CH4 產率隨入口氣速的變化Fig. 6 Variation of CH4 yield with inlet gas velocity

入口溫度540.0 K，n(H2)/n(CO)為2.5，不同入口氣速下CH4和CO質量分數在1.59 m高度處徑向分布如圖7 所示。由圖7 可知，隨著氣速的增加，CH4質量分數逐漸減小，而CO質量分數逐漸增大。當氣速較低時，氣體與催化劑顆粒接觸充分，反應更加完全。CH4質量分數與CO質量分數徑向呈相反分布，氣速過大，顆粒及反應物停留時間縮短，不利于反應，隨著氣速的增加CO質量分數兩邊低中間高的趨勢更加明顯。

圖7 不同入口氣速下CH4 和CO質量分數徑向分布（軸向高度1.59 m）Fig. 7 Radial distribution of CH4 and CO mass fractions at different inlet gas velocities (axial height 1.59 m)

入口氣速9.00 m/s，入口溫度540.0 K，n(H2)/n(CO)為2.5，反應器不同高度截面處CH4質量分數與顆粒體積分數徑向分布如圖8 所示。由圖8 可知，隨著高度的增加CH4質量分數逐漸增大，且逐漸反應完全徑向分布更加均勻?？梢钥闯?，CH4質量分數與催化劑顆粒體積分數密切相關，具有相同的分布規律（中間低兩邊高），催化劑顆粒體積分數大的地方CH4質量分數大。

圖8 不同高度下CH4 質量分數和催化劑顆粒體積分數徑向分布Fig. 8 Radial distribution of CH4 mass fraction and catalyst particle volume fraction at different heights

2.2 入口溫度對甲烷化過程的影響

甲烷化反應是強放熱反應，溫度過高會導致固體催化劑高溫燒結失去活性并抑制甲烷化正向反應的進行，因此需要合適的溫度來保證反應的進行。入口氣速7.50 m/s，n(H2)/n(CO)為2.5，不同入口溫度下CH4質量分數軸向分布如圖9 所示。由圖9 可知，CH4質量分數沿軸向逐漸增加，在爐膛底部甲烷化反應劇烈，CH4質量分數迅速增加。當入口溫度為490.0 K時，CH4質量分數約為0.10；當入口溫度大于540.0 K時，出口區域CH4質量分數約為0.20 并趨于一致?？梢?，溫度增加50.0 K時，CH4質量分數增加1 倍。隨著入口溫度的增加，CH4質量分數增加明顯，CH4產率的變化如圖10 所示?？梢钥闯?，隨著入口溫度的增加，CH4產率先增加后保持一致；當入口溫度大于510.0 K時，CH4產率增速變緩，其中540.0 K時，CH4產率約為490.0 K的兩倍。入口溫度的增加有利于甲烷化反應的進行，但必須保證溫度低于失活溫度，因此入口溫度取510.0~570.0 K。

圖9 不同入口溫度下CH4 質量分數軸向分布Fig. 9 Axial distribution of CH4 mass fraction at different inlet temperatures

圖10 CH4 產率隨入口溫度的變化Fig. 10 Variation of CH4 yield with inlet temperature

2.3 n(H2)/n(CO)對甲烷化過程的影響

n(H2+ CO) /n(N2)為3.0，改 變n(H2)/n(CO)，探究了其對甲烷化過程的影響。入口溫度525.0 K，入口氣速為7.50 m/s，不同n(H2)/n(CO)下CH4質量分數軸向分布如圖11 所示。當反應穩定后，隨著n(H2)/n(CO)的增大，CH4質量分數先增大后減小。當n(H2)/n(CO)為0.5~3.0，且比值較低時，H2不足以轉化CO，生成的CH4量很少，但隨著比值增大，n(H2)增大，更多的CO參與反應提高了CH4質量分數。當n(H2)/n(CO)為3.0 時，CH4質量分數軸向分布最大。n(H2)/n(CO)在3.0~4.5 時，隨著n(H2)/n(CO)增大，n(CO)減小，導致CH4生成量逐漸減小。由圖12 可知，隨著n(H2)/n(CO)增大，CH4產率及CO轉換率先增加后基本一致；當n(H2)/n(CO)小于3.0，CH4產率及CO轉化率都較低，n(H2)/n(CO)為2.5 時，CH4產率已接近n(H2)/n(CO)為3.0 時的計算值，兩者相差0.1016；當n(H2)/n(CO)為3.0 時，CO轉化率為0.9653；當n(H2)/n(CO)為3.5 時，CO轉化率為0.9957，已基本完全反應。整體當n(H2)/n(CO)大于3.0 時，可獲得較高CH4產率及CO轉化率，且考慮到煤氣產物中n(H2)/n(CO)一般較低，因此n(H2)/n(CO)取2.5~3.5。

圖11 不同n(H2)/n(CO)下CH4 質量分數軸向分布Fig. 11 Axial distribution of CH4 mass fraction under different n(H2)/n(CO)

圖12 CH4 產率和CO轉化率隨n(H2)/n(CO)變化Fig. 12 Variation of CH4 yield and CO conversion with n(H2)/n(CO)

2.4 響應曲面優化實驗

2.4.1 響應面實驗方案及結果

通過單因素對甲烷化過程影響結果分析，根據Box-Behnken 中心組合試驗設計原理，對入口氣速（因素A）、入口溫度（因素B）和n(H2)/n(CO) （因素C）各選擇三個水平，采用三因素三水平響應面分析法，以CH4產率為響應值建立回歸模型對操作參數進行優化，各因素及水平選擇如表2 所示。不同操作參數下模擬值如表3 所示。

表2 響應面優化因素及水平選擇Table 2 Response surface optimization factors and level selection

表3 響應面實驗方案及結果Table 3 Response surface experimental schemes and results

2.4.2 回歸模型建立及方差分析

通過Design-Expert軟件對上述數據進行統計及結果擬合，得到CH4產率對入口風速、入口溫度、n(H2)/n(CO)的二次回歸模型，方程表達式如式(14)所示。

對二次回歸模型進行方差分析，方差分析結果如表4 所示。方差分析表中P值決定所建立二次回歸模型是否顯著，當P< 0.0500 時，代表所建立回歸模型顯著，且P值越小模型越顯著；當P< 0.0100 表示差異極顯著[21]。由結果分析可知，CH4產率與各因素所建立二次回歸模型P< 0.0001，說明所建立回歸模型擬合結果極顯著，回歸效果良好。回歸模型的擬合質量用R2表示，多元相關系數R2越大相關性越好，且當R2≥ 0.8000 時，代表回歸模型擬合質量良好，該回歸模型的相關系數R2= 0.9914，調整后的相關系數R2= 0.9803，這說明該回歸模型具有較高的可靠性及顯著性。

表4 響應面實驗方差分析Table 4 Analysis of variance of response surface experiment

CH4產率與各操作因素的回歸模型中，單因素方面，入口溫度和n(H2)/n(CO)的P值均< 0.0001，影響顯著，而入口氣速大于0.05，影響不顯著。單因素對CH4產率影響顯著性依次為n(H2)/n(CO) >入口溫度>入口氣速，交互項中三個P值均> 0.0500，說明三者交互作用對CH4產率影響不大；二次項中入口溫度和n(H2)/n(CO)的P值< 0.0100，說明入口溫度和n(H2)/n(CO)對CH4產率影響較大。

通過Design-Expert軟件將表3 中CH4產率模擬值與回歸模型預測值進行擬合，通過擬合效果對模型可靠度進行驗證，CH4產率模擬值與回歸模型預測值關系如圖13 所示，圖中的每個點均代表一個模擬值，斜線為回歸模型預測值，點越集中分布于斜線周圍則回歸模型的擬合程度越好。模擬值在斜線兩側緊密均勻分布，說明該回歸模型可靠，對結果可以進行有效預測。CH4產率回歸模型預測值與殘差值分布如圖14 所示?；貧w模型預測值與殘差值獨立分布，且兩者沒有相互影響的趨勢。表2 響應面實驗序號與殘差值分布如圖15 所示。殘差值在零點坐標軸兩側分布均勻，且均小于殘差極限值，表明該回歸模型預測的準確性高。

圖13 CH4 產率模擬值與預測值關系Fig. 13 Relationship between CH4 yield simulated value and predicted value

圖14 CH4 產率預測值與殘差值關系Fig. 14 Relationship between CH4 yield predicted value and residual value

圖15 響應面實驗序號與殘差值分布Fig. 15 Response surface experiment serial number and residual value distribution

2.4.3 響應面分析

在對回歸模型的方差分析基礎上，通過Design-Expert軟件建立入口氣速、入口溫度和n(H2)/n(CO)對CH4產率的三維響應曲面圖。自變量間相互影響程度可通過響應曲面曲率的大小來確定，響應曲面曲率越大兩自變量間交互作用越大，且響應曲面上某方向的坡度較陡，則該方向對應的自變量因素對響應結果的影響較大[22]。各因素對CH4產率的交互影響曲面圖如圖16 所示。從曲面各方向的坡度陡峭程度來看，各因素對CH4產率的影響權重為：n(H2)/n(CO) >入口溫度>入口氣速，這與方差分析結果一致。各因素對CH4產率影響的等高線圖如圖17 所示，橢圓形表示兩影響因素交互作用顯著，而圓形則與之相反。交互因素影響均不顯著，這與方差分析結果一致。但當交互因素中含有n(H2)/n(CO)時，圖像出現較大變化，說明n(H2)/n(CO)這一因素比其余操作參數對CH4產率有更大的影響。

圖16 各參數對CH4 產率的交互影響曲面圖Fig. 16 Interaction surface diagram between parameters and CH4 yield

圖17 各參數對CH4 產率影響的等高線圖Fig. 17 Contour map of effect of various parameters on CH4 yield

2.4.4 操作參數優化

利用Design-Expert軟件對CH4產率與各操作因素的回歸模型進行了優化求解，得到最優操作參數入口氣速6.94 m/s、入口溫度556.8 K、n(H2)/n(CO)為3.3 下CH4產率為0.5880，為了驗證回歸模型優化結果的準確性。采用上述最優操作參數進行模擬計算，得到CH4產率模擬值為0.5696，與回歸模型預測值相差0.0184，相對誤差為3.13%，小于5.00%，說明CH4產率模擬值與回歸模型預測值之間擬合效果非常好，表明CH4產率與各操作因素的回歸模型具有較高準確性。

3 結論

采用CPFD法對六回路循環流化床甲烷化過程進行數值模擬，通過操作參數：入口氣速、入口溫度和n(H2)/n(CO)對反應器內CH4質量分數分布及CH4產率影響分析基礎上，采用響應曲面法建立了操作參數與CH4產率的回歸模型并對操作參數進行了優化，得到以下結論。

（1） CH4與CO質量分數沿徑向呈相反分布。甲烷質量分數與催化劑顆粒體積分數密切相關，兩者沿徑向呈“中間小兩邊大”分布。入口氣速過大，減小了反應物與催化劑接觸面積并縮短反應物停留時間，不利于反應進行。當n(H2)/n(CO)大于3.0時，可獲得較高CH4產率及CO轉化率。

（2）影響CH4產率指標單因素權重為n(H2)/n(CO)>入口溫度>入口氣速?；貧w模型優化操作參數為入口氣速6.94 m/s，入口溫度556.8 K，n(H2)/n(CO)為3.3，CH4產率為0.5880，在最優操作參數下CH4產率模擬值為0.5696，相對誤差為3.13%，小于5.00%，表明該二次回歸模型準確度較高。