異形催化劑床層中Sabatier反應的微-介尺度模擬

沈文豪，張亞新

（新疆大學化學化工學院，新疆烏魯木齊830046）

據國際能源署（IEA）報道，2018 年全球CO2排放量已增加至331 億噸，達歷史新高[1]。同時為解決可再生電能多效利用問題，有人提出通過電解水制H2，將電能以氫燃料形式儲存[2]。但存貯氫燃料危險性很大，CO2加氫開鍵還原產生甲烷（Sabatier 反應[3]）完美地解決了上述兩個問題，國際空間站也可利用Sabatier 反應回收CO2以產生維持人體生命特征的水[4]。

近年來，國內外對Sabatier 反應工藝的研究主要集中在催化劑的制備和反應機理上[5-7]，反應器的研究相對較少，針對此類強放熱反應，通常采用級聯絕熱固定床、流化床和三相漿態床等反應器。由于固定床具有結構簡單、成本低等優點，早已在工業生產中投入使用，例如合成氣甲烷化反應[8]。但容易造成飛溫，且催化劑的設計技術仍需不斷優化改進以提高轉化率，使固定床反應器仍具有較高的研究價值，段洪敏等[9]已經在裝載Ru/Al2O3催化劑的固定床反應器中，進行了寬范圍氫碳比的Sabatier反應拉偏實驗研究。

采用數值模擬方法對反應器進行研究被大家廣泛使用[10-12]，例如Ducamp 等[13]對裝載Ni/Al2O3催化劑的固定床反應器進行Sabatier 反應實驗和模擬研究，對比了CO2轉化率、CH4選擇性和溫升，誤差均小于10%，證明了模擬具有一定準確性。數值模擬的優勢在于可節省大量人工、實驗成本，并且能夠更直觀地展現實驗難以測試的催化劑內部反應過程[14]。應景濤等[15]對費托合成蛋殼型催化劑單顆粒內部反應體系進行數值模擬，通過研究顆粒內部反應、擴散和傳熱的基本規律，得到了目標催化劑較佳的活性組分厚度。床層建模困難和整床模擬計算量較大導致這些研究都是針對單一顆粒的研究，然而隨機堆積催化劑床層的研究對了解真實反應器內傳遞規律更為重要。國外學者Dixon 等[16-18]通過數值模擬方法對固定床床層的建模、空隙率分布、流體流動和傳熱規律，再到耦合反應動力學的多尺度模擬進行了探索，得到甲烷蒸汽重整等反應的基本規律。固定床的多尺度研究報道在國內很少，將隨機堆積固定床床層應用于Sabatier反應的研究更少，因此具有重要的研究意義。

本文首先將模擬結果和實驗結果進行對比，驗證模擬方法準確后，建立了與實驗同樣直徑D=31.7mm的圓筒，并使用DEM 方法分別在圓筒內建立隨機堆積167 個等體積V=13.5mm3的球形（dp=5.4mm）、圓柱形（dpc/h=1、1.3）催化劑顆粒幾何模型，再對床層和催化劑顆粒內的流體流動、物質傳遞、能量傳遞和Sabatier 反應動力學方程進行耦合計算，得到了不同形狀催化劑介觀床層和微觀顆粒內的三傳一反規律，并探討了改變催化劑形狀、降低壁溫和增加入口惰性球層、流速、溫度對熱點、CO2轉化率和顆粒內物質分布的影響，結論可為固定床反應器中Sabatier 反應的多尺度模擬和過程強化提供參考。

1 模型的建立

1.1 顆粒的隨機堆積

實際生產過程中，顆粒（催化劑）入料方式為重力作用下的自由落體運動，下落過程中顆粒與顆粒、顆粒與壁面之間的摩擦碰撞會使顆粒位置和角度產生隨機性。為近乎真實的模擬這一過程，采用DEM 方法，將顆粒材料設置為glass，壁面材料設置為steel，使顆粒在重力作用下，從床層頂部隨機下落堆積，設置停止標準為顆粒速度vp＜10-5m/s，再將生成的顆粒幾何文件導入CFD 軟件。我國某商業在售的甲烷化Ni 基Al2O3催化劑型號為J108-2Q 球形（Φ4.5～5.5mm）、W907 柱形（Φ(3～5)×(3～5)mm），結合Hwang 等[20]實驗所用催化劑尺寸，共建立了3 組分別堆積167 顆等體積球形（dp=5.4 mm）、柱形（dpc/h=1、1.3）顆粒的模型，其中球形、柱形（dp/h=1.3）顆粒的堆積模擬過程如圖1所示。

圖1 顆粒堆積模擬過程

采用Dixon[19]通過實驗數據擬合得到的隨機堆積顆粒床層空隙率計算公式進行驗證，經計算得到柱徑比dp/D=0.17 時球形顆粒床層空隙率為0.421，柱形顆粒床層空隙率為0.397。本文堆積的球形顆粒床層空隙率為0.432，柱形顆粒床層空隙率為0.411。驗證了使用DEM 方法建立的床層模型能夠代替真實催化劑床層，誤差主要來源于顆粒材料和床層高度取值的偏差。球、柱形顆粒床層空隙率計算表達式[19]分別見式(1)和式(2)。

1.2 床層模型

本文計算模型采用Hwang等[20]進行CO2甲烷化實驗的固定床反應器，床層示意圖如圖2(a)所示，半徑R=15.875mm，顆粒堆積高度為8R，CO2和H2從上端流入，生成的CH4和H2O與反應物混合后從下端流出。考慮到計算成本，本文除實驗驗證選用整床為多孔介質模型計算外，后續模擬過程都只建立了直徑D=2R、高H=2R的隨機堆積顆粒床層，進、出口分別設置高R、2R的區域以消除進出口效應的影響。將離散元軟件生成的顆粒幾何模型導入圓筒模型中，建立的球形顆粒床層模型如圖2(b)所示。

圖2 物理模型

1.3 數學模型

為研究固定床反應器中微-介尺度的三傳一反現象，已經建立了隨機堆積顆粒床層幾何模型，但顆粒內部微孔道很難實施建模，將其假設為多孔介質模型，在床層流體域和顆粒多孔介質域分別添加傳質、傳熱模型，假設床層流體域發生層流流動，求解動量守恒、能量守恒、質量守恒和反應動力學方程。

（1）動量守恒方程見式(3)和式(4)。

式中，I 為單位張量；ρ 為混合物密度，可由式(5)計算。

式中，M為混合物摩爾質量，kg/mol。

（2）能量守恒方程見式(6)。

式中，Cp為恒壓熱容，J/(kg·K)；λe為顆粒有效熱導率，W/(m·K)；wi為各物質的質量分數，%，表達式見式(7)。

式中，λm為流體熱導率，W/(m·K)；λp為固體顆粒熱導率，W/(m·K)。

（3）質量守恒方程見式(8)。

（4）反應動力學方程

在H/C=4 的情況下，幾乎可忽略副反應的影響，Sabatier反應占據主導地位，反應式見式(9)。

采用Rotaru 等[21]在Ni/Al2O3催化劑實驗中測定的CO2甲烷化本征反應動力學方程，表達式見式(10)。

式中，k 為速率常數；K 為吸附平衡常數；Keq為反應平衡常數，表達式見式(11)～式(14)。

CO2轉化率（XCO2）表達式見式(15)。

1.4 網格劃分與邊界條件

堆積顆粒之間的接觸點會增大網格數量，且接觸點附近的網格質量很差，會使計算不收斂。人們通常對顆粒進行縮小、切割、放大或橋接以解決接觸點問題[22]，但都會減小或增大顆粒體積，本文采用縮小的方法，將顆粒整體縮小1%，床層空隙率只下降1.2%，有研究表明對計算結果影響非常小[23-24]。

流體域和顆粒域都采用非結構化四面體網格，為使顆粒中參數計算更加精確，對顆粒域網格進行加密處理，在顆粒外壁面和內壁面各添加3層邊界層以詳細捕捉外壁面附近流體流動及內壁面附近反應、擴散過程。邊界層第一層厚度7.94×10-4m，拉伸因子1.2，厚度調節因子2，網格劃分如圖3所示。由于柱形顆粒線性幾何更為復雜，本文只對dpc/h=1的柱形顆粒床層進行網格無關性驗證，結果如表1所示。證明使用7236918單元網格計算已足夠精確，繼續細化對結果并無影響。采用同樣方法進行后續網格劃分，dpc/h=1.3 的柱形顆粒床層網格數量為6861283，球形顆粒床層網格數量為4679633。

圖3 網格劃分

表1 網格無關性驗證

采用速度進口、壓力出口邊界條件，入口為充分發展的流動，入口H/C=4，vin=0.0126m/s。顆粒內無對流，由于進出口壓差較小，假設顆粒內壓力恒為P=3atm（1atm=101325Pa），將催化劑顆粒設置為反應源和熱源，孔隙率εp=0.44，熱量通過壁面與外界進行交換，為保證催化劑保持較高活性，初始壁面溫度和入口溫度都為T=623.15K。考慮到多場耦合中非線性方程組求解計算量大，先單獨計算穩態流場，再將其他物理場耦合計算，采用直接迭代PARDISO 算法求解，瞬態求解500s，步長為10s，使用默認松弛因子值。

2 結果與討論

2.1 模型準確性驗證

建立Hwang 等[20]實驗所用反應器床層幾何模型，考慮計算成本原因，將完整床層假設為多孔介質模型，采用與實驗相同的基本參數和條件，進口溫度可調節，反應器壓力1atm（1atm=1.013×105Pa），H/C=5∶1，只改變進口流率（2.4～12.0mL/s），得到不同進口溫度（498.15K、523.15K、573.15K）下CO2轉化率與質量空速倒數（mcat/Fin）的關系。如圖4所示，采用上述建模方法得到的模擬結果與實驗值吻合較好，證明模擬準確可靠。為探究真實顆粒堆積床層及顆粒內部微觀反應特征規律，本文采用上述經驗證的模擬方法，只假設顆粒為多孔介質模型進行后續模擬分析。

圖4 模擬值與實驗值對比

2.2 異形顆粒床微-介尺度場分布

2.2.1 CO2摩爾分數分布

圖5 球、柱形顆粒床層CO2摩爾分數、流線分布對比

實際生產的Ni/Al2O3催化劑形狀多種多樣，其中球形、柱形顆粒比較典型，對比3種形狀顆粒床層內CO2摩爾分數和速度流線，如圖5所示。由圖5(a)可知，由于反應物接觸的顆粒越來越多，反應越來越充分，球形顆粒床層CO2摩爾分數沿軸向逐層減少；反應物沿顆粒半徑方向由外向內擴散，且擴散路徑中會在活性物質作用下不斷被消耗，導致顆粒內CO2摩爾分數呈現外層低內層高的環狀趨勢。由圖5(b)～(c)可知，由于反應和傳遞機理相似，dpc/h=1、1.3柱形顆粒CO2摩爾分數分布也類似，但部分顆粒內部CO2摩爾分數比球形低，這是因為雖然3 種顆粒當量直徑相同，反應物從球形顆粒外表面擴散到中心的距離都相等，而反應物從柱形顆粒外表面擴散到中心的距離可能長于球形顆粒，導致部分沿床層軸向擴散路徑較長的柱形顆粒內反應物停留時間更長，正向反應轉化更充分。圖5(c)表明dpc/h=1柱形顆粒床層CO2摩爾分數在壁面為20%左右，這可能是由于dpc/h=1柱形顆粒壁面空隙率較大，形成溝流，大部分反應物流經床層壁面時并未參與反應。

提取3 種床層同一位置顆粒內CO2摩爾分數分布情況如圖6所示，由于物質在顆粒內層停留時間長，反應更充分，球形和柱形顆粒內CO2摩爾分數都呈現外高內低的環狀趨勢，受流體流動方向的影響，環狀分布并不對稱。CO2摩爾分數減小趨勢在球形顆粒內一直呈圓環狀，而在柱形顆粒外層先呈類似磨平棱角的矩形環狀，再向顆粒中心逐漸過渡為圓環狀，這是由于棱角處活性物質少、反應速率小，CO2較集中，導致柱形顆粒棱角處會出現低轉化率區。

2.2.2 床層溫度分布

柱形顆粒和球形顆粒床層溫度分布規律幾乎一致，都會出現熱點區域，并隨反應進行向出口端移動，本文只提取球形顆粒場分布進行分析。當t=500s時，沿反應器軸向中心截面熱源分布如圖7所示，放熱量從床層進口端到出口端依次降低，顆粒邊界處放熱量高于顆粒中心，表明床層進口端和顆粒表面反應更為劇烈，反應速率更大。截面溫度場分布如圖8所示，床層進口端反應不穩定，溫度梯度大，由于壁面與外界換熱的存在，導致反應器靠近壁面的溫度低于床層中心溫度，床層出口區由于無熱量產生尤為明顯；流體總體流動方向至上而下，導致床層進口端產生的反應熱更趨于向床層下端傳遞堆積。結合上述原因，床層中心會出現一個約881K的球形熱點區域。

圖6 球、柱形顆粒內CO2摩爾分數分布

圖7 截面熱源分布

圖8 截面溫度場分布

提取床層中心軸線上溫度每隔100s 的變化，如圖9所示，在外界溫度、壁面換熱系數不變的條件下，持續反應產生的熱量會不斷堆積使床層整體溫度上升。當t=100s時，床層溫度延中心軸線先大幅升高，再緩慢降低，這是由于床層進口端反應更劇烈，此時，球形熱點區域在距床層進口端2R 的位置。隨反應的進行，入口反應逐漸趨于穩定，床層進口端溫升梯度逐漸減緩，熱量隨流體作用向出口端推移，t=500s 時熱點已移動至距床層進口端3R的位置，幾乎在床層末端。

2.3 操作條件影響分析

從床層介尺度角度出發，通過改變操作條件來提高反應器性能，能大大減少研發成本。基于上述球形顆粒床層模型，選擇了顆粒形狀、床層壁面溫度、入口惰性顆粒、入口流速和入口溫度作為變量，探討了操作條件對CO2分布和熱點溫度的影響。初始床層模型共5層帶有活性的催化劑，入口CO2摩爾分數80%、H2摩爾分數20%，入口流速vin=0.0126m/s，壁面溫度和入口溫度都為T=623.15K，操作壓力P=3atm。

圖9 軸線溫度隨時間的變化

2.3.1 顆粒形狀、壁溫的影響

出口CO2轉化率和平均溫度隨時間的變化如圖10 所示，CO2轉化率分3 個階段變化：第一階段為0～10s，轉化率幾乎呈斷崖式增大，這可能與床層內突然變化的溫度有關，由于數據存儲步長為10s，這一過程在本節不做詳細描述；第二階段為10～200s，由于溫度上升，不利于放熱反應正向進行，使CO2轉化率增大趨勢不斷減小，直至200s增至峰值26%左右；第三階段為200～500s，隨溫度的升高，CO2轉化率逐漸下降。反應不斷放熱，操作條件不改變，熱點溫度會持續上升，前200s 由于轉化率在不斷增大，熱點溫度上升越來越快；200s 時熱點溫度達到761K 左右，較入口溫度高138K；200s以后熱點溫度上升趨勢逐漸趨于平緩，這是由于200s 后床層溫度過高，不利于正向反應進行，導致反應放熱量減少。

圖10 球、柱形顆粒床層反應特性隨時間的變化

由圖10 可知，等體積情況下，當t=200s 時，球形顆粒比dpc/h=1 柱形顆粒床層CO2轉化率高0.6%，放熱更強，導致球形顆粒比dpc/h=1 的柱形顆粒床層熱點溫度高7K，這說明球形結構的催化劑顆粒有效利用率更高，與上一節得出的結論一致。dpc/h=1.3 柱形顆粒比dpc/h=1 柱形顆粒CO2轉化率更高，這是由于dpc/h=1.3柱形顆粒床層壁面孔隙率更小，流通過的反應物少。為了改善200s 以后轉化率降低的現象，將200s 以后的外界環境溫度從623.15K 降低至573.15K，反應至500s。由圖11可知，當外界環境溫度突然降低后，會導致床層內傳熱突然不穩定。由于降溫使平衡正向移動，產熱增大，熱點溫度小幅突增。此后由于換熱加強，熱點溫度較原來小，500s 時床層熱點溫度較之前減少10K，出口CO2轉化率增大2%。

圖11 降低壁溫的影響

2.3.2 入口條件的影響

由上一節可知，球、柱形顆粒床層整體特征規律一致，本節以球形顆粒床層為代表，探究入口條件對顆粒內物質分布的影響。圖12 為t=500s時，床層中心軸線上各物質摩爾分數分布情況，反應物和生成物的摩爾比都符合Sabatier 反應的計量系數比；顆粒中反應物呈現外層高中心低的規律，生成物恰好相反，使床層內物質整體分布規律呈波浪形式。如圖12 所示，H2整體分布波浪振幅比CO2小，顆粒內外H2摩爾分數差更值小，且每個顆粒內H2含量最低值向出口端偏移，而CO2含量最低值在顆粒幾何中心，這都是由于H2有效擴散系數較大，擴散速度更快的原因。CO2摩爾分數在進入床層前會有小幅增大，這也是由于H2有效擴散系數較CO2大，導致少量H2先進入顆粒內部參與反應，CO2便堆積在床層前端，圖5 也表現出此現象，床層前端CO2摩爾分數約為20.5%，但流體流過第一層顆粒后，富CO2現象就會消失，表明無規律流動會消除由于內擴散差異帶來的CO2聚集現象。

圖12 軸線各物質摩爾分數分布圖

圖13 增加惰性球層對CO2分布的影響

圖14 入口流速對CO2分布的影響

入口區增加兩層同樣大小的26 顆惰性球層，顆粒內CO2分布變化如圖13所示，由于改變了進入反應區的流體流動方向，導致第一層具有活性的催化劑前富CO2區明顯減小，有助于平衡正向移動；顆粒中心CO2摩爾分數減小，反應更充分，最終導致出口處CO2摩爾分數減小。

入口流速每增加1倍時CO2分布變化如圖14所示，增大流速加快了H2和CO2的流動，富CO2區也會減少，但是流速增大會使反應物在顆粒內停留時間減小，造成前兩層顆粒內CO2摩爾分數成倍增大；床層后端由于受入口流速增大的影響小，顆粒內部CO2摩爾分數波動不大，但出口總體CO2摩爾分數明顯增大，轉化率下降。

圖15 入口溫度對CO2分布的影響

入口溫度每增加20K 時CO2分布變化如圖15所示，入口溫度升高，使第一層催化劑活性增大，顆粒內CO2轉化更充分，CO2摩爾分數減小，第二層之后的區域由于反應放出大量熱，不利于反應正向進行，顆粒內CO2摩爾分數升高，轉化率降低，最終導致出口CO2摩爾分數升高，轉化率下降。

3 結論

分別建立167顆等體積的隨機堆積球形、柱形催化劑床層幾何模型，使用經實驗驗證準確的模擬方法，對催化劑微尺度和床層介尺度的Sabatier 反應和傳遞數學模型進行耦合模擬計算，主要得到以下結論。

（1）柱形催化劑中心CO2轉化率比球形更高，但催化劑棱角和反應器壁面會出現低轉化率區，導致球形催化劑床層出口CO2轉化率更高。將柱形催化劑外壁的絕對線形改為更符合催化劑內物質擴散規律的弧形，能夠使物質剛擴散進入催化劑時起到很好的均勻分布和過渡作用，減少棱角低轉化率區和不利流態。

（2）柱形催化劑床層內流體溝流、回流和滯留現象比球形更突出，床層中部都會出現近似球形的熱點區域，并隨反應的進行向床層出口緩慢移動，催化劑內反應物濃度呈現外層高內層低的環狀分布。由于擴散系數的不同，導致床層前端會產生富CO2區。

（3）CO2轉化率隨反應的進行會先增大后減小；降低壁面溫度50K，熱點溫度下降10K，CO2轉化率增大2%；增加惰性球層和入口流速，會阻止床層前端富CO2區的形成；增加惰性球層、減小入口流速和溫度，能提高催化劑內和床層出口CO2轉化率。

符號說明

c—— 濃度，mol/m3

D—— 床層直徑，mm

Dim—— 擴散系數，m2/s

dp—— 球形催化劑直徑，mm

dpc—— 柱形催化劑直徑，mm

h—— 柱形催化劑高度，mm

l—— 床層高度，mm

P—— 壓力，Pa

Qr—— 反應熱源，W/m3

R—— 床層半徑，mm

r—— 反應速率，mol/(m3·s)

T—— 溫度，K

V—— 催化劑體積，mm3

v—— 床層流體速度，m/s

vp—— 催化劑下落速度，m/s

x—— 摩爾分數，%

ε—— 床層空隙率

εp—— 催化劑空隙率

μ—— 動力黏度，Pa·s

ρ—— 混合物密度，kg/m3

下角標

ave—— 平均值

i—— CO2,H2,CH4,H2O

in—— 進口

out—— 出口