甲烷化反应器催化剂积炭过程的模拟研究

2020-08-07沈文豪张亚新

高校化学工程学报 2020年3期

沈文豪, 张亚新, 赵玲

(1. 新疆大学化学化工学院, 新疆乌鲁木齐 830046;2. 华东理工大学化学工程联合国家重点实验室, 上海 200237)

1 前言

我国煤炭资源丰富，国家多项政策指明了煤炭行业供给侧改革的发展方向是清洁高效利用[1]。大规模煤制天然气项目已在新疆、内蒙古等富煤地区投产，该工艺中甲烷化反应段对甲烷产率影响最大。副反应中甲烷吸热会产生丝状炭附着在催化剂内孔道表面，长时间的积累会导致催化剂内孔道严重堵塞[2-3]。张沛等[4]采用扫描电子显微镜(SEM)等表征方法，对商用选择性催化还原(SCR)脱硝催化剂的失活和再生进行研究，发现砷和硫酸盐严重堵塞了催化剂微孔道，酸洗再生能有效脱除砷和硫酸盐，提高SCR 脱硝催化剂活性。KOO 等[5]采用不同Mg/Al 混合比控制的载体，提高了甲烷重整制合成气催化剂的抗结焦能力。DEMICHELI 等[6-7]在实验中测定了甲烷在Ni/Al2O3催化剂上形成积炭的反应动力学方程，使数值模拟方法研究甲烷化反应积炭过程成为可能。

采用 CFD 方法对固定床反应器内“三传一反”(动量传递、质量传递、热量传递、化学反应过程[8])规律进行研究能够节省大量的实验和人工成本，并获得反应器内部更详细的特征参数[9-10]。CHEIN 等[11]采用数值模拟方法研究了固定床管式反应器中甲烷化反应的影响因素，发现CO 甲烷化反应对温度最敏感，提高反应压力和降低空速能提高甲烷产率。SON 等[12]采用数值模拟方法研究了操作条件对积炭效应的影响，但模拟也仅停留在反应器宏观尺度和催化剂-流体域介观尺度。对固定床反应器的模拟研究通常将催化剂床层简化为多孔介质模型或者规则排列的小球[13-14]，耦合催化剂内和床层流体域的多尺度研究也几乎只在国外文献中报道，DIXON 等[15-17]在催化剂内部反应过程的研究中做了很大贡献，但大都针对于蒸汽重整等反应。掌握降低积炭效应的技术，必需建立更趋于真实的随机堆积固定床床层，并将介观流体域和微观催化剂域的多物理场耦合，研究多尺度下甲烷化反应催化剂积炭形成过程。

本文首先使用DEM 软件生成颗粒(催化剂)直径dp=0.012 m、床层与颗粒直径比N=4、高径比H/D=5的随机堆积球形颗粒床层，并提取颗粒位置信息，通过编程在CFD 软件中建立床层几何模型，将计算域分为床层流体域和催化剂颗粒域，建立“三传一反”多物理场耦合数学模型，通过数值模拟分析了微观颗粒域的积炭形成规律。最后讨论入口温度400、450、500 K，催化剂温度770、810、850 K，入口流速0.005、0.010、0.015 m·s-1和入口氢气占比0、20%、40%分别对甲烷化反应催化剂积炭过程的影响。

2 模型的建立与验证

2.1 几何模型

在随机堆积颗粒的固定床反应器中，较大管径比将导致床层颗粒数倍增，给数值求解带来巨大的计算量，而小管径比固定床反应器在强放热、吸热反应和实验室研究阶段应用也比较广泛，两者径向流动均匀性存在差异，但最终整体特征规律仍具有一定相似性。

本文建立了小管径比固定床反应器床层，结构如图1 所示，甲烷从上端进入，下端流出，床层总高为12 dp，直径为4 dp，将床层分为入口、催化剂堆积和出口3 个区域，入口区域高2 dp，催化剂堆积区域高5 dp，出口区域高5 dp。使用DEM 软件在催化剂堆积区域生成直径为dp=12 mm 的球型颗粒，并提取颗粒位置信息，建立几何模型如图2 所示，DEM 参数设置如表1。

图1 床层结构示意图Fig.1 Schematic diagram of the bed structure

图2 床层堆积模型Fig.2 Schematic diagram of the bed packing model

2.2 数学模型

2.2.1 控制方程

甲烷化反应在颗粒内放出的热量成为积炭反应热源，颗粒内由无数微孔道组成，现有建模手段和计算能力很难实现颗粒内成千上万微孔道的建模，积炭又覆盖在内孔道表面，改变了催化剂的孔隙率。因此，做出以下合理假设：

(1) 流体为理想气体；

(2) 流体为不可压缩、稳态的层流；

(3) 催化剂颗粒为多孔介质模型；

(4) 反应所需要的热量来自颗粒；

(5) 积炭对催化剂活性的影响非常小可忽略，主要影响孔隙率的变化。

基于以上假设，建立多物理场耦合数学模型，控制方程如下：

(1) 动量守恒方程：

表1 DEM 参数设置Table 1 Parameters for DEM

积炭会降低催化剂孔隙率，进而影响渗透率，初始渗透率采用多孔介质Kozeny-Carman 渗透率模型，描述催化剂孔隙率和渗透率变化的表达式[18]如下：

式中：MC为积炭的摩尔质量，0.012 kg·mol-1；ρdepos为积炭的密度，2 000 kg·m-3。

2.2.3 边界条件与计算方法

入口初始甲烷摩尔分数为1，温度为500 K，平均流速为0.015 m·s-1，出口压力为101 325 Pa，将反应源添加在催化剂颗粒域，壁面设置为无滑移和热绝缘边界条件，催化剂区域设置为热源，温度为850 K。催化剂初始孔隙率采用Dixon 实验测得的Ni/Al2O3球型催化剂孔隙率0.44[19]。

考虑到多物理场耦合中非线性方程组求解的难度，首先将物理场解耦，单独计算稳态流场，最后将稳态流场与其他物理场耦合计算，采用直接迭代PARDISO 算法求解，瞬态求解8 000 s，步长为60 s，使用默认松弛因子值。

2.3 验证

2.3.1 几何模型准确性验证

早期学者提出球型固定床床层整体空隙率[20]与径向空隙率分布[21]的经验公式，并与实验数据吻合较好，得到学者们的广泛认可。整体空隙率经验公式如下：

经计算，本文中DEM 堆积的床层整体空隙率为0.47，拟合公式计算空隙率为0.44，相差3%，误差主要来源于堆积的随机性以及计算空隙率时床层高度的取值。

径向空隙率经验公式如下：

图3 径向空隙率截取示意图Fig.3 Schematic diagram of the radial void fraction

图4 径向空隙率与文献对比Fig.4 Comparison of the radial void fraction results with literature

2.3.2 网格无关性验证

由于反应主要发生在催化剂堆积区，为减少计算量将进出口区域划分为六面体网格，催化剂堆积区均划分为四面体网格，筒体和颗粒边界划分3 层边界层。反应、扩散等行为在颗粒域更为复杂，对其网格进行加密处理，由于在颗粒中心施加温度边界，为更精确计算边界处参数，将颗粒中心网格再次细化，网格划分如图5 所示。直接划分网格会导致壁面接触点附近网格质量急剧降低，将所有颗粒直径减小1%，床层空隙率只下降1.2%，对结果影响较小。采用同样的网格划分方法划分1 040 000、1 570 000 和2 060 000网格，比较3 种网格床层压降变化，如图6 所示，结果表明，继续细化网格对结果并无较大影响，为降低计算成本，将网格数量确定为1 570 000。

图5 网格划分Fig.5 Mesh generation

图6 网格数量对压力变化的影响Fig.6 Effects of grid number on pressure variation

3 结果与讨论

3.1 模拟有效性验证

KUVSHINOV 等[23]在填装Ni/Al2O3催化剂，外部持续加热维持等温状态的管式反应器中进行甲烷积碳反应实验，测得H2占比为0、0.15、0.30、0.40 时，温度从760 K 增加到860 K 时，反应速率的变化规律。本文采用CFD 方法同样进行等温反应的模拟，将模拟值与文献值对比，如图7 所示，在氢气质量分数为0 时，模拟值与试验值有很小的误差，其余参数下二者结果基本吻合，证明模拟方法准确可靠。

图7 模拟值与实验值对比Fig.7 Comparison of simulation with experimental values

图8 多物理场分布Fig.8 Multiple physical field distribution

3.2 多物理场分布

将流体域与颗粒域多物理场耦合计算，得到如图8 所示的多物理场分布示意图。从速度流线图可以看出，流体流动在出入口区很稳定，在床层堆积区会沿颗粒表面无规则流动，局部速度会增大到入口流速的10 倍，颗粒尾部还会出现尾流和回流现象，流体流经颗粒多孔介质区域，在刚接触颗粒、离开颗粒和颗粒内的速度几乎为0。反应器壁面(右)是压力分布，压降发生在床层堆积区，进出口压降为0.104 Pa。反应器壁面(左)是温度分布，入口甲烷温度为500 K，随着距入口距离的增加，床层区温度逐渐上升，第1 层颗粒温度从500 K 迅速上升至850 K，到第2 层颗粒时温度基本达到850 K。颗粒表面是反应刚开始时甲烷浓度分布，甲烷在颗粒上半部分浓度高，下半部分浓度低，在入口区域浓度为14.4 mol·m-3，底部最小浓度能达到 14 mol·m-3。

3.3 场分布随时间的变化

3.3.1 热态场分布

分别提取了不同时刻床层中心截面温度、甲烷摩尔分数、反应速率和颗粒中心截面孔隙率的分布。图9 为不同时刻床层中心截面温度分布，当t=0 s 时，入口温度为500 K，颗粒中心温度为850 K，其他区域温度均为293.15 K。当t=60 s 时，各区域温度基本达到稳定，床层堆积区和流体区温度迅速上升至850 K，但由于入口区不断有500 K 的甲烷流入，使入口区和床层堆积区产生500～850 K 的温度梯度。

图9 温度随时间的变化Fig.9 Temperature profiles as a function of time

图10 甲烷摩尔分数随时间的变化Fig.10 Profiles of methane mole fraction as a function of time

图10 为甲烷摩尔分数随时间的变化，当t=60 s时，颗粒内少部分甲烷已经开始反应。当t=2 000 s时，由于甲烷向催化剂内层扩散，扩散阻力增大，部分甲烷反应生成了碳和氢气，甲烷摩尔分数开始减小，最小为0.93；气流方向至上而下，单个颗粒上半部分甲烷供给更快，导致甲烷在颗粒上半部分多、下半部分少。当t=4 000 s 时，颗粒内层甲烷摩尔分数最小为0.90，此时外层孔道已经严重堵塞，颗粒内部传质受外部气流影响减小，导致甲烷在颗粒内近似呈圆环状减小。当t=8 000 s 时，颗粒几乎完全堵塞，甲烷在颗粒内呈圆环状减小的趋势更明显，不同的是颗粒内层甲烷摩尔分数最小已下降至0.85。

图11 为反应速率随时间的变化，当t=60 s 时，最小反应速率为 0.006 1 mol·m-3·s-1；当 t=2 000 s时，颗粒内甲烷被消耗减少，最小反应速率减小至0.005 3 mol·m-3·s-1；当 t=4 000 、8 000 s 时，最小反应速率分别为 0.004 8、0.004 0 mol·m-3·s-1，分布规律同甲烷摩尔分数分布规律相同，反应速率在颗粒内也近乎呈圆环状减小。

图11 反应速率随时间的变化Fig.11 Variation of reaction rate as a function of time

3.3.2 孔隙率分布

如图12 所示，反应从60 s 进行到8 000 s，颗粒内平均孔隙率值从0.44 下降到0.04左右。随反应的进行，颗粒外层比内层积炭更严重，当t=4 000 s 时，颗粒内、外层孔隙率值最大相差0.05，对比图12 和11，发现反应速率越快的区域，孔隙率越小。提取了相同时间间隔(1 000 s)床层中心轴线上颗粒内的孔隙率分布，如图 13 所示，表明反应越充分，颗粒内孔隙率分布越不均匀，呈现出外层低、内层高的规律。提取前 8 000 s每1 000 s 颗粒内平均孔隙率下降值，如图14 所示，表明孔隙率下降程度随时间的增加在不断减缓，4 000 s 前，催化剂未完全堵死，反应剧烈，平均每1 000 s 颗粒内平均孔隙率值下降0.08；4 000 s 之后，颗粒内平均孔隙率值下降到0.1 以下，积炭反应减慢，平均每1 000 s 颗粒内平均孔隙率值下降0.025。

图12 中心截面孔隙率随时间的变化Fig.12 Variation of central section porosity as a function of time

图13 中心截线孔隙率随时间的变化Fig.13 Variation of central line porosity as a function of time

图14 相同时间间隔孔隙率下降趋势Fig.14 Decrease of porosity under same time intervals

3.4 稳态时场分布规律

根据3.3 节描述，发现时间进行到8 000 s 时反应程度已经很小，为增加结论的准确性，将计算时间增加到12 000 s，提取稳态时甲烷、氢气摩尔分数和颗粒内孔隙率分布云图。如图15 所示，当t=12 000 s 时，流体区甲烷摩尔分数高达 99%，氢气含量很少，从颗粒外层到内层，甲烷含量越来越低，氢气含量越来越高。如图16 所示，当t=12 000 s 时颗粒内平均孔隙率值为0.015，此时颗粒内、外层孔隙率值相差0.02，颗粒外层积炭仍然更严重。这种现象是颗粒内扩散和反应同时造成的，颗粒外层反应速率大，生成的氢气和碳也越多，同时大量甲烷也会稀释氢气，使颗粒外层的氢气减少；颗粒内层甲烷内扩散的同时也不断地反应产生氢气，最后集中在颗粒中心，内层反应程度反而越来越小，积炭程度也越来越小。

图15 甲烷和氢气摩尔分数对比（t = 12 000 s）Fig.15 Comparison of mole fraction of methane and hydrogen (t = 12 000 s)

图16 孔隙率分布(t = 12 000 s)Fig.16 Porosity distribution (t = 12 000 s)

3.5 不同操作参数对积炭的影响

实际生产中，最容易控制的是操作参数，且通过改变操作参数来提高生产效率经济成本最低。由图16 可知，每一个颗粒内孔隙率分布都表现出对称性，表明颗粒中心截线上孔隙率的变化能够反映颗粒内平均孔隙率的变化规律。因此，提取当t=12 000 s，稳态时，床层中心轴线上颗粒内孔隙率的变化，研究了不同入口氢气占比、温度、流速和催化剂温度对积炭效应的影响，如图17～20 所示。

如图17 所示，入口组分无氢气时，颗粒内平均孔隙率值为0.02，积碳非常严重，加入氢气，由于增加了生成物使促使反应逆向进行，颗粒内平均孔隙率值明显提高；当氢气占比 40%时，颗粒内平均孔隙率值增大到0.29，同时，加入氢气会使颗粒内、外层甲烷浓度梯度减小，积炭分布越均匀。如图18 所示，颗粒内温度为770 K 时，颗粒内平均孔隙率值为0.23；温度从770 K 上升到810 K，颗粒内平均孔隙率下降约0.15；温度从810 K 升高到850 K，颗粒内平均孔隙率值下降约0.06，表明催化剂升温会加速积炭的形成，继续升温对积炭的影响越来越小。如图19 所示，入口速度从0.005 m·s-1增加到0.010、0.015 m·s-1对颗粒内平均孔隙率几乎无影响，但颗粒中心处孔隙率值会下降0.001 左右，这是由于增大流速会增大反应器轴向压差，在更大压力驱动下相同时间甲烷内扩散越多，生成的积炭也越多。如图20 所示，入口温度的升高只导致第一层颗粒中心的孔隙率减小，因为入口温度越高，进入床层的气体能够更快达到反应区温度，更快生成积炭。这与温度场分布规律一致，由图10 可知，第1 层颗粒温度变化梯度较大，第2 层温度基本达到稳定。

图17 入口氢气含量对孔隙率的影响Fig.17 Effects of inlet hydrogen contents on porosity

图18 催化剂温度对孔隙率的影响Fig.18 Effects of catalyst temperature on porosity

图19 入口速度对孔隙率的影响Fig.19 Effects of inlet velocities on porosity

图20 入口温度对孔隙率的影响Fig.20 Effects of inlet temperature on porosity

图21 更直观的反应了4 种操作参数对颗粒内平均孔隙率的影响，入口氢气从0 增加到40%，颗粒内平均孔隙率值从0.02 上升至0.29，增加了0.27。催化剂温度从770 K 增加到850 K，颗粒内平均孔隙率值从0.24变为0.02，减小了0.22。当入口流速从0.005 m·s-1增加到0.015 m·s-1时，颗粒中心孔隙率减小 0.01，平均孔隙率几乎无变化。当入口温度从400 K 增加到500 K，第1 层颗粒中心孔隙率减小0.005，其他颗粒中心、平均孔隙率均无变化。