基于双尺度耦合模拟的甲烷化催化剂多级孔结构探究

叶 银，张亚新

新疆大学化工学院，新疆 乌鲁木齐 830046

我国是天然气消费大国，天然气消耗量逐年增加，对进口天然气的依赖度偏高。二氧化碳加氢制甲烷工艺既能实现温室气体高效绿色资源化利用，增加天然气的供应，又能为弃风、弃光等不稳定新能源电力的转化存储提供一种有效方式[1]。所以从环保、能源战略安全以及能源产业结构优化的角度出发，CO2甲烷化技术具有极大的研究价值。

甲烷化过程是合成天然气工艺中的重要一环，具有反应迅速、热效应强的特点[2]，要让反应能够顺利进行，并且兼顾反应物的转化率和甲烷的产率，就必须考虑催化剂的作用。这些年来，从催化剂的制备和成型角度出发，对活性中心、助剂以及异形结构等方面已进行了大量的研究[3-5]，但对颗粒中的扩散限制关注较少。程永刚等[6]的研究表明，在反应过程中，甲烷化催化剂的主要反应区域为颗粒的外层，催化剂内部空间未能得到充分利用，通过一定的手段降低颗粒的内扩散影响将对反应过程产生积极影响。早期的研究[7-8]表明，催化剂颗粒最佳孔结构的设计非常重要。随着实验技术的发展，通过合理设计并合成具有理想孔隙结构、性能更佳的催化剂已成为可能。因此，多级孔结构的设计已成为研究热点[9]，利用大孔隙改善物质扩散，小孔隙提供更多的反应区域。

就多尺度研究而言，通过实验手段探究催化剂的孔隙结构对反应的影响是较为困难的，且脱离反应器去研究催化剂性能无法直观得到其工业价值，目前已有一些反应器-催化剂双尺度耦合模拟的研究。肖瑶等[10]采用几何映射法将甲醇制丙烯（MTP）固定床反应器与颗粒耦合，探究了催化剂粒径对反应过程的影响；Shi等[11]利用反应器-颗粒耦合模型，在反应器尺度下优化了重油加氢脱金属（HDM）催化剂的孔结构；曹军等[12]在考虑产物液膜传质阻力的情况下，构建了能反映反应器与催化剂颗粒双尺度上费托合成产物分布的多场耦合数学模型，讨论了产物液膜及催化剂粒径对产物分布的影响。

本研究运用双尺度耦合方法建立了甲烷化固定床反应器的床层-催化剂颗粒双尺度数学模型，同时考虑了床层及颗粒内的热量和物质传递以及化学反应过程；在一定的操作条件下，评价了催化剂颗粒的多级孔隙结构参数（大孔平均孔径、介孔平均孔径、大孔孔隙率和介孔孔隙率）对颗粒中反应及传质的影响，为催化剂的孔结构设计确定了优化区间。

1 数学模型

1.1 床层-催化剂颗粒耦合计算模型

本研究建立的甲烷化反应器模拟模型可以获得催化剂颗粒尺度和床层尺度内的物质传递状况、温度变化及化学反应过程，床层-催化剂颗粒耦合模型如图1 所示。该模型与传统均质模型的区别在于它考虑了颗粒相的存在，因此需要用不同的方程组去定义流体相和颗粒相，通过加载颗粒表面的边界条件将颗粒中的浓度、温度分布与流体相的相关梯度相耦合，颗粒模型与床层模型的控制方程及边界条件如下。

图1 物理模型示意Fig.1 Schematic diagram of physical

颗粒模型的控制方程：

针对球形催化剂颗粒，采用Wakao 和Smith 提出的宏观-微观孔隙模型[13-14]，该模型用大孔平均孔径（dM）、介孔平均孔径（dm）、大孔孔隙率（εM）和介孔孔隙率（εm）四个参数来描述催化剂颗粒的孔隙结构。颗粒的反应比表面积（Sg，单位质量催化剂的表面积）和密度（ρ）定义与孔隙结构直接相关[15]，见式（6）。在该模型中，同时考虑了克努森扩散和分子扩散，有效扩散系数为浓度、温度和孔径分布的函数[16]。采用扩散-反应二阶偏微分方程建立颗粒内的反应-传质模型，见式（1），各组分扩散仅沿颗粒半径方向，并视颗粒内部为恒定温度。Solsvik 等[17]研究发现，颗粒内部的黏性流和压降很小，可以忽略不计，因此在颗粒方程中不含有对流项和动量平衡。

床层模型为对流-扩散方程，模型主要考虑床层轴向的变化。保证物料守恒的同时，床层模型还计算了沿床层的温度变化，Qr为反应热，Qvd表示床层与壁面的换热量。由于物质穿过颗粒-流体界面进入颗粒的流量可能是由流体主体一侧对质量传递的阻力决定的，因此在该模型中提出了薄膜条件假设[18]，hD表示薄膜传质系数，可以根据流动中对流和扩散质量传递的相对贡献的舍伍德数来计算。

1.2 甲烷化反应动力学模型

模型的计算要求耦合化学反应，甲烷化反应过程较为复杂，且伴随多个副反应，本研究仅考虑如下三个主要的反应：

上述反应的动力学方程采用Xu 等[19]建立的模型，该模型是基于Ni-MgO-Al2O3催化剂的实验测量得到的本征动力学模型。

式中：下角标1，2，3 分别代表CO 甲烷化反应、CO2甲烷化反应以及水煤气变换反应。甲烷的产率和选择性定义为：

1.3 模型求解方法及验证

上述数学模型采用COMSOL 5.4 多物理场耦合仿真软件求解。通过COMSOL 软件中的组件耦合模块，使得床层尺度和颗粒尺度的实时数据在两个维度间进行数据投影映射，从而实现催化剂与床层模型的耦合。计算过程中，提取床层中某点的各项物理场数据，以此作为该处催化剂颗粒模型的边界条件，同时把床层模型得到的物质浓度数据传递给颗粒方程，在颗粒维度对催化剂颗粒中的反应速率、物质浓度进行求解，得到相应数据；同时，在薄膜条件的假设下，颗粒内部的物质浓度数据经过处理换算为颗粒表面的物质通量，从而将相应数据传递至床层尺度，以此确定床层方程中反应速率；各物理场数据通过耦合变量在两个尺度间相互换算，两尺度间的数据传递与模型求解交替进行，直至模型收敛求得反应器内各物理场分布。

本研究为瞬态模拟，模拟时长定为1 000 s。由于模型高度对称，只对1/8 反应器空间进行模拟计算。颗粒维度采用立方根序列进行离散，单元数定为6。三维反应器网格的划分采用扫掠划分，沿径向作等差分布，经网格无关性验证，最终确定网格数为23 860。

为检验本工作所选用动力学模型的可靠性，用上述动力学模型计算得到的CO 甲烷化反应和CO2甲烷化反应的模拟数据与文献数据[20-21]进行了比较，结果如图2 所示。由图2 可知，模拟数据与实验结果基本吻合，证明模拟准确可靠。

图2 模拟值与实验值对比Fig.2 Comparison of simulated values and experimental data

2 结果与讨论

本研究模拟的反应器尺寸参数参照Hwang 等[21]进行CO2甲烷化实验的固定床反应器的，催化剂粒径（dp）为5.4 mm 的球形颗粒，床层总体孔隙率（εbed）用经验公式（21）[22]进行计算，其余初始条件如表1 所示。

表1 初始条件Table 1 Initial conditions

2.1 参考模型的模拟结果

本研究需先确定一个参考模型，以此作为后期多级孔结构优化设计的对比基准。参考模型的颗粒总体孔隙率（εtotal）设定为0.5，颗粒大孔平均孔径（dM）为100 nm、大孔孔隙率（εM）为0.1；颗粒介孔平均孔径（dm）为5 nm、介孔孔隙率（εm）为0.4。计算所得反应器出口处CO2转化率如图3 所示。模拟至900 s 左右时反应器内各个物理场趋于稳态。提取1 000 s 时，反应各组分在床层尺度的浓度分布及温度分布如图4 所示。

图3 CO2 转化率随反应器运行时间变化Fig.3 Change of CO2 conversion with running time

图4 1 000 s 时床层中各物质浓度及温度分布Fig.4 The concentration and temperature distributions of various substances in the bed at 1 000 s

计算得到反应各组分在反应器不同位置处催化剂颗粒内部的浓度分布如图5 所示，由床层中心入口处（R=0 mm，L=15.9 mm）颗粒内组分浓度分布可以看出，CO2和H2进入颗粒后迅速反应，浓度快速下降；由于存在微量水汽变换反应和CO 甲烷化反应，CO 浓度先小幅升高，而后降低；在颗粒中心处，各组分浓度变化趋于平缓。床层中部（R=0 mm，L=30.9 mm）和出口端（R=0 mm，L=45.9 mm）的催化剂颗粒内组分浓度分布与床层入口端处（R=0 mm，L=15.9 mm）相比，浓度变化渐渐平缓，颗粒内的反应速率明显降低。

图5 反应器不同位置处颗粒内的物质浓度分布Fig.5 The concentration distributions of components in particles at different positions of reactor

图6 为床层入口处催化剂颗粒内各组分反应速率的径向分布。图7 为反应器中心处床层浓度及颗粒浓度的轴向分布。可以看到，颗粒内部与床层间存在较大的浓度梯度。反应可能发生在催化剂颗粒近外壁面区域，约为颗粒纵深20%的区域内，催化剂中心处反应速率较低，催化剂利用率较低，因而甲烷化反应是受内扩散严重影响的过程。从催化剂制备的角度来说，大孔孔隙将会降低内扩散阻力，使各物质更容易在颗粒内部进行输运，从而更好地利用颗粒内部空间，所以研究催化剂颗粒的多级孔结构是有必要的。

图6 颗粒内各物质反应速率的径向分布Fig.6 The reaction rate distribution of components in particles along radial direction of the reactor

图7 反应器内床层浓度及颗粒浓度的轴向分布Fig.7 Axial concentration distributions in the bed and in particles at the centerline of the reactor

2.2 多级孔结构探究

球形催化剂颗粒的宏观-微观孔隙模型[13-14]中，模型参数有大孔平均孔径、介孔平均孔径、大孔孔隙率和介孔孔隙率，本研究在设定颗粒总体孔隙率为0.5 的情况下，考察了颗粒结构参数对催化剂活性的影响，寻求更有效的催化剂颗粒孔结构。

保持大孔孔径为100 nm，介孔孔径为5 nm，模拟了床层入口处颗粒内CO2径向浓度分布随大孔孔隙率的变化情况，结果见图8。由图8 可知，大孔孔隙率的增大改善了颗粒内部的扩散阻力，提高了反应程度；但是，当大孔孔隙率大到一定程度时（εM=0.5），由于反应比表面积减少很多，反应程度反而会降低。不同大孔孔隙率时固定床出口的CO2转化率、CH4选择性和产率的模拟结果见图9。由图9 可知，反应器性能随着大孔孔隙率的增加先逐渐变好而后又急剧变差。当大孔孔隙率为0.4 时，甲烷的产率达到81.1%，相较于参考模型提高了约16%。这是由于大孔孔隙率控制着大孔和介孔的分布，影响反应的扩散和比表面积，如果颗粒孔隙率是由较高比例的大孔隙组成的，扩散速度会更快，但反应表面活性较低。因此如果大孔孔隙率过大，反应组分在催化剂颗粒内的扩散速度加快，导致反应速率会先增大，但在某个值上活性表面的损失太高，转化率将会降低。

图8 大孔孔隙率对颗粒内CO2 浓度分布的影响Fig.8 Effect of macroporous porosity on the distribution of CO2 concentration in particles

图9 大孔孔隙率对反应结果的影响Fig.9 Effect of the macroporous porosity on the reaction results

图10 为颗粒的比表面积（Sg）和CO2有效扩散系数（DCO2eff）随大孔孔隙率的变化。如图10 所示，随大孔孔隙率的增大，比表面积呈线性减小，扩散系数呈二次增大。因此，在活性成分没有改变的情况下，若在造粒过程中对颗粒孔隙结构进行优化，可制造出性能优异的催化剂。

图10 反应比表面积和CO2 有效扩散系数随εM 的变化情况Fig.10 The variations of reaction specific surface area and CO2 effective diffusion coefficient with εM

在孔隙率分布固定的情况下（εM=0.4，εm=0.1），模拟研究了大孔平均孔径和介孔平均孔径对反应的影响。CO2转化率、CH4产率、CH4选择性和热点温度随大孔孔径的变化如图11 所示。由图11 可知，随着dM的增大，转化率和选择性有所增加，当dM≥120 nm 时，产率增幅趋于平稳。因为克努森扩散系数与孔径为线性关系，随着大孔孔径的增大，有效扩散系数变大，从而导致产率有所增加，但反应比表面积随dM的增大而减小，则会抑制反应，而产率趋于平衡可能是由于在所探究的大孔孔径范围内，扩散的影响占主导，弥补了反应比表面积减小的影响。

图11 大孔孔径对反应结果的影响Fig.11 Effect of the macroporous pore diameter on the reaction results

介孔孔径对反应结果的影响如图12 所示。由图12 可见，CO2转化率以及CH4产率都随介孔孔径的增大而降低。虽然颗粒的介孔孔径越大，其比表面积越大[见式（6）]，但扩散系数随dm的增大而增大[见式（3）]，而扩散的影响起了主导作用。同时，考虑到热点温度不超过积炭温度（823 K）[23]，介孔直径为4.5 nm 左右比较合适。

图12 介孔孔径对反应结果的影响Fig.12 Effect of the mesoporous pore diameter on the reaction results

3 结论

a）利用多物理场耦合计算方法所建立的床层-颗粒双尺度耦合模型能直观地反映床层及颗粒内的传递行为以及在不同尺度上的化学反应过程。

b）床层-颗粒双尺度耦合模拟结果表明：在未经优化的球型颗粒（dp=5.4 mm）中，甲烷化反应受催化剂中内扩散阻力的影响，主要发生在催化剂颗粒近壁面区域，有效反应区域约为颗粒外层的20%区域，催化剂整体利用率较小；且颗粒内与床层间存在较大浓度梯度。

c）在给定的操作条件下，基于双尺度耦合模拟对球形甲烷化催化剂（dp=5.4 mm，εtotal=0.5）的孔结构探究发现：反应器甲烷产率随颗粒大孔孔隙率增大呈现先增大后减小的趋势，大孔孔隙约为0.4 时甲烷产率达到最大；甲烷产率随大孔平均孔径的变化波动较小，当dM≥120 nm 甲烷产率趋于平稳；介孔平均孔径越小，甲烷产率越大，当dm减小至4.5 nm 时产率变化变缓。εM和dm对反应结果的影响较大，采用多级孔结构催化剂有助于实现甲烷产率的提升；当εM为0.35～0.4，dM≥120 nm，dm为4～5 nm 时，甲烷产率为80%左右，较参考模型提高约16%。

符号说明