异形催化剂床层中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—— 出口