六回路循环流化床甲烷化过程数值模拟及优化

马继伟，张亚新,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 所示。回归模型预测值与残差值独立分布，且两者没有相互影响的趋势。表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)>入口温度>入口气速。回归模型优化操作参数为入口气速6.94 m/s，入口温度556.8 K，n(H2)/n(CO)为3.3，CH4产率为0.5880，在最优操作参数下CH4产率模拟值为0.5696，相对误差为3.13%，小于5.00%，表明该二次回归模型准确度较高。