化学链燃烧中燃料气体在CaSO4(010)表面吸附扩散特性的分子动力学模拟

2020-03-04侯封校刘敦禹寇学森杨浩然王永贞

石油学报（石油加工） 2020年6期

侯封校，金晶，刘敦禹，寇学森，杨浩然，王永贞

(1.上海理工大学能源与动力工程学院上海市动力工程多相流与传热重点实验室，上海 200093；2.临沂大学土木工程与建筑学院，山东临沂 276000)

温室效应是全世界所面临的巨大挑战，化石燃料的广泛使用造成了CO2的大量排放，加剧了温室效应。为减缓全球气候变暖的趋势，碳减排势在必行。目前，碳减排的方式可以分为两类：一是采用清洁、低碳的供能方式；二是对产生的CO2进行捕集、封存或再利用。由于目前新能源等产能方式仍有一定的局限性[1-2]，且中国化石能源消费占比非常高，煤、石油、天然气的消费占比可达到58.25%、19.59%、7.43%[3]，因此，对传统化石能源燃烧所产生的CO2进行捕集对中国意义更大。碳捕集的方式可以分为燃烧前捕集、燃烧后捕集、富氧燃烧，但这些技术均会产生额外耗能，且运行成本较高[4-6]。而化学链燃烧通过载氧体的循环将传统的燃烧过程分为2个阶段，以此实现氮气和氧气的分离，且没有额外的分离耗能，还可实现能量的梯级利用，因此具有较好的发展前景[7-8]。

在化学链燃烧中，燃料与载氧体在燃料反应器中接触并进行反应，该反应器也被称为还原反应器。其内为还原性气氛，固体燃料进入还原反应器首先热解产生燃料气体，燃料气体与载氧体进行气-固异相反应；热解产生的焦炭一部分会被气化，另一部分仍会与载氧体接触，进行速率较低的固-固反应[9-11]。因此，载氧体在气-固和固-固反应中具有重要作用。

目前载氧体分为金属氧化物载氧体与非金属氧化物载氧体，其中非金属氧化物载氧体CaSO4具有载氧能力出色、价格低廉且易获得、环保性较好等优点，因此有较大的发展潜力。沈来宏等[12]从热力学平衡角度对比了CaSO4与金属氧化物载氧体Fe2O3、NiO，结果表明，Fe2O3载氧体仅在 Fe2O3-Fe3O4阶段有较好的反应性能，而CaSO4载氧体与NiO载氧体的反应性能接近，但CaSO4载氧体载氧量高、价格低廉，是较为理想的载氧体。Song等[13-14]分别采用固定床与流化床进一步探究CaSO4载氧体的性能，虽然热力学上CaSO4载氧体反应性能较好，但其反应速率较低；提高反应温度可以提高其反应性能，但当反应温度为1273 K时，会有6%左右的SO2析出，因此需通过增加反应时间来提高转化率。在反应温度为1223 K的实验中，CaSO4载氧体在反应初期表现了良好的反应能力，但是在反应末期其反应性能逐步降低，且积炭现象愈发严重，同时还有轻微的烧结现象。Deng等[15]采用H2为燃料进行CaSO4载氧体还原反应的数值模拟，结论一致。较低的反应温度与较大的载氧体粒径使鼓泡床中快速形成较大的气泡，从而使燃料与载氧体接触不充分，导致燃料转化效率较低。对于以煤为燃料的化学链燃烧反应，Zheng等[16]认为燃料反应器内的反应温度应低于1223 K，以防止CaSO4载氧体烧结或过多分解。根据CaSO4载氧体反应能力强但反应速率低的特点，掺杂Ni、Fe、Cu等对载氧体进行改性，可以在一定程度上增大载氧体的比表面积，同时降低硫释放的比率[17-19]。Zhang等[20]对比了CaO与CaCO3对CaSO4载氧体的硫释放特性的影响，表明CaO的固硫效果更好，且Ca/S摩尔比越高，固硫效果越好。而Zheng等[21]发现，在空气反应器中仍会有SO2释放，从反应平衡的角度，通过添加CaO可降低硫释放，在1173 K的反应温度下，CaO/CaS摩尔比为1∶1时可以将硫损失降到最低。Abad等[22]以煤作为燃料、CaSO4/CaO作为载氧体的实验结果表明：调控载氧体与燃料的比例可降低硫释放，该比例取决于固体燃料的含硫量；而对于固体燃料的在线气化化学链燃烧而言，CaSO4载氧体的氧传输性能是其他廉价载氧体(如钛铁矿)的4倍。综上，硫释放已不是制约钙基载氧体应用的主要因素，对于以CaSO4为载氧体的化学链燃烧系统，最为关键的问题是增大还原反应的反应速率。因此，揭示CaSO4载氧体反应速率较低的原因，对后续载氧体改性、操作条件优化都非常重要。但宏观实验只能从反应物、产物量变的角度对反应进行研究，研究的反应为总包反应，忽略了内部的吸附、扩散与反应等细节。Zhang等[23-24]对以CO与CH4为燃料的还原反应阶段进行了计算研究，说明了这2种燃料气体在CaSO4(010)表面的反应路径。但其主要针对反应阶段进行研究，而在反应前燃料气体分子在CaSO4载氧体上的吸附扩散行为未有报道。由于分子动力学在研究载氧体表面吸附扩散等特性时，可考虑温度、压力等宏观因素，因而被广泛应用[25-26]。

综上，为了研究燃料气体分子在CaSO4载氧体表面的吸附扩散特性，笔者采用CO作为典型燃料气体，使用分子动力学模拟了其在反应温度为1123 K、1173 K与1223 K时的动力学特性，分析了CO在CaSO4(010)表面的吸附行为与扩散特性，深入探讨了其作用机理；同时与CH4、H2进行了对比，研究了不同气体分子吸附扩散行为的机理。研究的结果可为CaSO4载氧体的改进与优化提供理论基础。

1 模型与方法

硫酸钙晶体在微观层面上最稳定的晶体表面为(010)面，采用3层CaSO4周期性晶体构成的CaSO4(010)表面作为气-固反应界面，其构型如图1所示。CaSO4(010)表面的大小为3.7428 nm×4.1946 nm，3个方向均采用周期性边界条件。为了表现CaSO4(010)表面的固体特性，在模拟中认为CaSO4(010)表面为刚性体系，保持其原子位置固定。此外，模拟分别采用30个CO、H2、CH4气体分子构成气相系统，气相系统z方向的高度(垂直于CaSO4(010)表面方向)采用理想气体方程进行确定。为了减少计算量，系统压力采用2.5331 MPa，由理想气体方程可得，在反应温度为1123 K、1173 K、1223 K时z方向高度分别为12.6761 nm、13.1987 nm、13.7213 nm。

图1 CaSO4(010)表面模型与燃料气体分子的几何构型Fig.1 The configuration of CaSO4(010) surface model and fuel-gas molecules(a) CaSO4(010)； (b) CO； (c) CH4； (d) H2

采用GROMACS软件[27]进行模拟，非键相互作用采用L-J(Lennard-Jones)12-6势描述，键合相互作用采用简谐势描述，其力场参数使用UFF(Universal Force Field)势[28]，CO与CaSO4(010)表面的相互作用能采用式(1)描述。

(1)

动力学模拟中，使用蛙跳算法对牛顿运动方程进行积分，时间步长为0.5 fs，模拟时间为10 ns，轨迹输出时间间隔为2 ps。范德华作用采用Cut-off方式进行计算，截断半径为1.5 nm；静电作用计算采用PME(Particle-Mesh-Ewald)方法，模拟中计算静电作用所需的原子电荷参数采用Mulliken电荷。在温度模拟中，采用钙基载氧体典型反应温度1123～1223 K[13-16]，通过Nose-Hoover热浴进行控温。根据能量输出结果，模拟中体系在2 ns左右已经达到平衡，因此对平衡后的8 ns进行采样分析。

2 结果与讨论

2.1 反应温度变化对CaSO4(010)表面吸附扩散特性的影响

选取CO作为典型燃料，考察反应温度变化时，CaSO4(010)表面对固体燃料热解产生的可燃性气体吸附扩散特性的影响。基于10 ns的NVT系综动力学模拟，选取体系平衡后的6 ns时间段内的轨迹进行能量分析。CO与CaSO4(010)表面之间的作用可以分为库仑作用与L-J非键作用，图2统计了在 6 ns 时间段内不同温度时两者作用能的概率分布。由图2可知，3个温度下，大部分能量分布均在能量小于0的范围内，表示CO在CaSO4(010)表面的吸附为放热反应。图2(a)表示库仑作用能概率分布情况，反应温度为1123 K时曲线峰值对应的作用能最低，1173 K时的次之，1223 K时的作用能最高；且各能量范围所对应的分布概率值也符合以上规律，1123 K时CO与CaSO4(010)表面的库仑作用能大概率存在于能量较低的位置，而1223 K时的高分布概率库仑作用能主要位于能量值接近于0的范围内。随着温度上升，高分布概率库仑作用能向横坐标正方向偏移，且值逐步趋于零，可知随着反应温度的升高库仑作用减弱，CO分子距离CaSO4(010)表面的距离增大。图2(b)为L-J非键作用能分布图，其规律与库仑作用能一致，高概率对应的作用能范围均随温度的上升向横坐标正方向偏移，表示其放热量逐渐减少。同时，随着反应温度的上升，L-J非键作用能峰值逐步下降，且峰值对应的作用能逐步向横坐标正向移动，说明随着温度升高，非键作用能减弱，表明该非键作用主要为范德华作用。而1173 K与1223 K时峰值所对应的概率较为相近。这是由于在模拟过程中，非键作用所包含的范德华作用随距离增大迅速衰减，因此在计算范德华作用时采用了 Cut-off 截断处理，超过截断距离的范德华作用为零。综上可知，CO在CaSO4(010)表面的吸附为物理吸附，随着反应温度的升高，其吸附能力逐步下降。

图2 不同反应温度下CO吸附于CaSO4(010)的能量概率(PE)分布图Fig.2 Energy probability (PE) distribution of CO adsorption on CaSO4(010) at different temperatures(a) Coulomb interaction； (b) Non-bond interaction

为了进一步分析CO在CaSO4(010)表面的吸附能力的变化，图3统计了CaSO4(010)表面法线方向CO分子的密度分布情况，统计时间范围为动力学平衡后的4 ns时间段内，z方向截断(Cut-off)距离为6 nm。图3中，z轴高度0～0.9389 nm范围为CaSO4(010)表面固体区域，大于0.9389 nm的范围为CO、H2、CH4气体区域。由图3可知，在CaSO4(010)表面附近富集的CO分子数远大于远离CaSO4(010)表面区域的CO分子数，该现象证明了CO分子在CaSO4(010)表面发生了明显的吸附现象。在z轴高度处于0.9389～1.6 nm范围内，通过与固体区域z值进行差减，得到在距离CaSO4(010)表面约0～0.6611 nm的区域，CO的密度有较大的峰值；在z轴高度大于1.6 nm的范围，即距离CaSO4(010)表面大于0.6611 nm区域，CO的密度保持稳定，基本为定值。这表明CaSO4(010)对CO的吸附作用产生了明显的分层现象，且在距离其表面约0.66 nm处为分层点。在距离CaSO4(010)表面小于0.66 nm的区域，由于其表面对CO分子的库仑作用与范德华作用较强，产生了明显的吸附作用，造成了CO分子的聚集并使CO分子运动性能减弱，从而形成高密度区域，即稳定吸附于CaSO4(010)表面，因此该区域为吸附层；而在距离CaSO4(010)表面大于0.66 nm的区域，CO分子的运动受CaSO4(010)表面吸附作用较小，其运动性能改变较小，从而密度分布较为均匀，因此该区域为体相层。随着温度的升高，吸附层中的CO密度减小，表明CaSO4(010)表面对CO吸附能力下降。该现象符合物理吸附随温度上升吸附能力减弱的特征，与相互作用能分析结果一致。而在体相层中，反应温度的变化对CO密度的影响较小，这是由于体相层密度较小，空间范围较大，吸附层分子逃逸至体相层后，对体相层密度的影响较为微弱。

图3 不同反应温度下CO吸附于CaSO4(010)垂直方向密度分布Fig.3 Vertical density distribution of CO adsorption on CaSO4(010) at different temperatures

在体系动力学模拟平衡后，对2 ns时间段内的轨迹进行采样，通过计算CO分子轨迹的均方位移，并结合Einstein公式计算以z方向为法线方向的xy平面二维扩散系数，其计算式如式(2)、(3)、(4)所示，燃料气体分子在xy平面方向3种温度下的扩散系数计算结果如表1所示。由表1可知，随着反应温度的升高，xy平面方向的扩散系数增大。该现象表明在反应温度为1123 K时，CO分子紧密吸附在CaSO4(010)表面，限制了其在xy平面方向的运动。当反应温度升高后，CaSO4(010)表面对CO分子的吸附能力减弱，CO分子因受到吸附力减少，其分子热运动加剧，更易迁移至新的反应位点。该结论与图3随反应温度升高CaSO4(010)表面吸附能力下降的结论一致。同时，该结论也表明CO在xy平面方向扩散符合分子碰撞理论，与CaSO4(010)表面对CO分子束缚能力(吸附能力)有较强的关联性，二者呈反比例关系。

表1 燃料气体分子在xy平面方向扩散系数(D)Table 1 The diffusion coefficient (D) of fuel-gasmolecules in xy plane

(2)

(3)

(4)

其中，D表示扩散系数，m2/s，下角标x、y表示扩散方向；t表示关联时间，s；t0表示初始时刻的时间，s；Δx、Δy分别表示x、y方向的位移，m；xi、yi分别表示第i个粒子的坐标；N表示粒子数。

由于CO分子在CaSO4(010)表面的吸附为物理吸附，因此CO分子与CaSO4(010)表面之间的弱相互作用为产生相互吸引的主要原因。基于平均非共价作用(average Noncovalent Interactions, aNCI)分析[29]，模拟平衡后统计1 ns时间段内体系的平均约化密度梯度(average Reduced Density Gradient, aRDG)散点图如图4所示，其中，纵坐标为平均约化密度梯度值，横坐标为sign(λ2)ρ，为电子密度函数Hessian矩阵本征值符号(sign(λ2))与电子密度(ρ)的乘积。该分析采用波函数分析软件包Multiwfn[30]进行。由图4知，在sign(λ2)ρ值为-0.04～0.04范围内出现了散点图的穗(Spike)，表明CO在CaSO4(010)表面的吸附是由弱相互作用主导的体系，进一步证明了该吸附是物理吸附。同时，在横坐标0～0.02的区域有明显的Spike，该现象表明体系中存在位阻作用，即表明该体系的弱相互作用不是由纯粹的范德华作用组成，静电作用所造成的弱相互作用更强。对比3个反应温度下的图像可知，随着反应温度的升高，Spike逐步减小，表明体系的弱相互作用逐步下降，CaSO4(010)表面对CO分子的吸引力变弱。横坐标-0.04～0区域的Spike对应范德华作用，该区域Spike随着反应温度升高迅速变小，表明范德华作用迅速衰减；靠近横坐标-0.03处的Spike对应较强的吸引作用，在该体系中即为静电作用，其随着反应温度升高也快速衰减。这2处Spike的变化规律与能量分布规律保持一致，说明反应温度升高使CO与CaSO4(010)表面相互作用减弱，导致CaSO4(010)表面对CO分子吸附作用降低。

图4 CO在CaSO4(010)表面在不同反应温度下aRDG散点图Fig.4 The aRDG scatter maps of CO adsorption on CaSO4(010) surface at different temperatures(a) 1123 K； (b) 1173 K； (c) 1223 K

2.2 不同热解燃料气体分子在CaSO4(010)表面的吸附扩散特性

为进一步研究化学链燃烧燃料反应器中固体燃料热解产生的不同燃料气体在CaSO4(010)表面的吸附扩散性质，选取了CO、CH4、H23种典型燃料气体进行模拟分析，其能量概率分布如图5所示。H2与CaSO4(010)表面的作用能均分布于横坐标轴的0点附近，表明H2吸附于CaSO4(010)表面所释放能量较小；而CH4的库仑作用能几乎分布在接近于横坐标轴0点处，非键相互作用能分布于横坐标轴-50～0处。相较于CO的能量分布，反应温度为1173 K时H2与CH4能量曲线峰面积小于CO，即表明H2、CH4与CaSO4(010)表面的作用能较弱，进而可知H2、CH4与CaSO4(010)表面的非键相互作用主要为范德华作用，二者在表面的吸附为物理吸附。根据图6的垂直方向气体密度分布可知：CaSO4(010)表面对CH4的吸附能力较弱，在吸附层吸附的CH4分子数量远低于吸附的CO分子数量，大部分CH4分子均分布于体相层；CaSO4(010)表面对H2吸附能力最弱，H2在吸附层中物质的量最小，在体相层中最多，因此相对于吸附层，体相层中H2的体积占比更高。3种气体分子所处的吸附体系中，吸附层与体相层的分层点仍为距离 CaSO4(010) 表面约0.66 nm处。

图5 不同燃料气体在CaSO4(010)表面的能量概率分布图Fig.5 Energy probability distribution of different fuel-gases adsorption on CaSO4(010)(a) Coulomb interaction; (b) Non-bond interaction

图6 不同燃料气体吸附于CaSO4(010)垂直方向密度分布Fig.6 Vertical density distribution of the different fuel-gases adsorption on CaSO4(010)(a) Mass; (b) Amount of substance

由表1可知，CH4与H2在xy平面方向的扩散系数均大于CO，表明二者在横向的扩散能力较强。其原因主要由于2个方面：一方面，被吸附的分子越多，其在扩散时受到的阻碍碰撞越多，使其扩散系数下降；另一方面，根据理想气体动力理论，气体分子的直径越小，质量越小，扩散系数越大。因此，三者的扩散能力由大到小依次为：H2、CH4、CO。

反应温度为1173 K时，CH4和H2在CaSO4(010)表面的aRDG散点图如图7所示。通过对比反应温度为1173 K时3种燃料气体的aRDG散点图(图7和图4)可知，CH4的静电作用要低于CO；而H2的位阻作用在3种燃料气体中最大，CH4次之，因此CO与CaSO4(010)表面结合较为紧密。