基于密度泛函理论的煤催化气化制甲烷过程研究

董脉帆，冯乐乐

（中国矿业大学 安全工程学院，江苏 徐州 221116）

近年来我国天然气消费量和需求量不断增加，对国外进口的依赖性持续加重，严重威胁我国能源安全[1-3]。而我国能源结构具有“富煤、贫油、少气”的特点，因此，煤制天然气是适合我国国情的煤炭利用技术。当前，煤制天然气技术主要分为两大类，即煤气化转化工艺和直接合成天然气工艺。前者的技术路线是先利用煤气化制备合成气，再通过气体转换反应、甲烷化反应制备天然气[4]；后者则是通过催化、加氢等方式直接一步合成煤基天然气[5]。

当前采煤技术的极限深度约为1500 m，而我国大约70%的煤炭资源分布在2000 m以深，大量煤炭难以通过传统的“建井-采煤-发电/化工”方式加以利用[6]。地下煤气化是指不经开采直接将煤层原位气化的技术，有望解决深部煤炭资源开采利用的难题，受到了研究者广泛关注[7-9]。地下煤气化是以煤层内的钻孔和燃烧形成的空腔作为反应场所，导致其气密性远不如地面上的成型反应器，因此无法安全地使用加氢气化工艺。此外，由于催化剂难以回收，无法采用价格昂贵的合成催化剂[10-11]。即便如此，研究者仍希望尽可能提高产物中的甲烷含量，从而使地下煤气化更好地满足煤制天然气的迫切需求。研究[12-13]表明，含有碱和碱土金属的生物质气化灰可以有效促进煤气化过程中H2等气体的生成；含有碱土金属的废弃蛋壳对煤气化也有显著的催化作用，可以增加H2和CO的产量。因此，含有碱和碱土金属的常见可弃催化剂成为众多学者研究煤催化气化制甲烷的重点关注对象。王会芳等[14-15]通过地下气化模型实验台，发现了含有Na+的废碱液对煤焦气化的碳转化率和甲烷化有明显的催化作用。张德成等[16]研究了CO2气氛下酒糟对煤焦气化的影响，发现酒糟中的Ca2+起到催化作用，提高了反应活性和甲烷产率。杜留娟等[17]研究发现，在煤气化过程中添加CaO可以大幅提高合成气中的甲烷含量。CHEN等[18]在地下煤气化模型的气化炉中使用较为廉价的Ca(OH)2作为催化剂，发现合成气中的甲烷含量增加。

为了掌握煤气化制甲烷的反应机理，贾建波等[19]利用密度泛函理论（DFT）探究了煤热解过程中的甲烷生成路径，发现反应温度对甲烷生成路径有显著影响。曾凡桂等[20]采用量子化学计算研究了霍林河褐煤热解产生甲烷的过程，总结了氧接脂碳类物质的热解、长链脂肪烃物质的二次热解、与芳香环体系直接或间接相连的甲基热解断裂和脂肪类物质的芳香化反应等4类反应途径。ZHANG等[21]研究了高挥发分烟煤热解产生甲烷的机理，发现甲烷生成过程中氢自由基辅助路径的活化能远低于其他反应路径。WANG等[22]利用量子化学计算探究了在低于热解温度下无烟煤产生甲烷的机理，发现在外界应力作用下，煤中含有的氢原子可以促进脂肪族中碳键的解离，有利于甲烷的形成。在煤的反应性研究方面，从头计算法、半经验法和DFT等方法经常被使用[23]。从头计算法存在计算效率低等问题[24]，半经验法存在计算繁琐、结果准确性较差等问题[25]，而DFT通过变量间的体系电子密度取代多电子波函数，极大地简化了计算过程，并且计算过程中考虑到了电子的相关效应，提高了计算结果的准确性[26]。

综上，研究人员证实了Na+和Ca2+对煤气化制甲烷的促进作用，但鲜有关注CO2气氛下煤基甲烷生成的催化反应机理。为此，本文基于DFT，选择乙基苯作为煤的模型化合物，研究在CO2氛围中的甲烷生成反应位置和路径，分析Na+和Ca2+对甲烷生成反应的催化作用，讨论不同温度下有无催化作用下的甲烷生成反应的活化能和速率，以期为煤气化制甲烷的催化剂选取提供参考。

1 模型与计算方法

1.1 煤的模型化合物

由于煤的分子结构较为复杂，为了能够从分子水平角度对煤催化气化制甲烷的反应过程进行探究，常采用具有部分煤结构特征的模型化合物来间接研究其反应机理。朱学栋等[27]发现甲烷的生成与煤结构中脂肪类结构有直接关联，且甲烷产量随着脂肪类结构含量的增加而提高。乙基苯是一种较为简单的化学结构模型（图1）。乙基苯既包含一个芳香环结构，也含有长脂肪烃侧链，非常具有代表性，是煤催化气化的较理想模型分子。因此，为了简化煤催化气化过程，本文选择了乙基苯作为煤的模型化合物，气化剂设置为CO2。

图1 乙基苯分子模型Fig.1 Ethyl benzene molecular model

1.2 计算方法

DFT是一种研究分子体系中电子结构的方法。其将能量与电子密度联系起来，并利用二者之间的交换-相关泛函表达式，使计算过程更加简便[28-30]。因此本文采用DFT对煤催化气化制甲烷的反应过程进行研究。相关计算在Gaussian软件中完成，分子结构在GaussView软件中构建。

本研究采用DFT方法中的B3LYP/6-31G(d,p)基组对反应物、过渡态物和生成物进行优化和参数表征。在计算过程中，基于优化后的过渡态结构，使用本征反应坐标（IRC）方法，寻找过渡态物所在的基元反应中的反应物和生成物。在检查优化结果过程中，需要判别反应物和生成物是否不存在虚拟频率以及过渡态物是否存在唯一虚拟频率且虚拟频率的振动方向是否对应基元反应中的反应物和生成物。基于自由体积理论，也就是在反应物和产物的分子数发生不对等的情况下，需要对所优化过程中计算出的能量值进行10.89 kJ/mol或-10.89 kJ/mol的修正。由于Gaussian软件的开发设置，其所计算背景被设置在298.15 K和101.33 kPa环境下，且无法进行修改设定。

1.3 热力学和动力学参数

单点能和吉布斯自由能是用来表征反应活性中的重要的特征参数。反应活化能可以定义为基元反应中的过渡态物的单点能减去反应物的单点能，同时，也表示反应发生所需要克服的最大能量障碍，计算公式如式(1)所示。

式中，Ea为反应活化能，kJ/mol；ETS为过渡态物的单点能，kJ/mol；EIS为反应物的单点能，kJ/mol。

单点能变（ΔE，kJ/mol）的大小可以反映出反应发生的难易程度，计算公式如式(2)所示。

式中，EIM为过渡态物或生成物的吉布斯自由能，kJ/mol。

吉布斯自由能变（ΔG，kJ/mol）用于表示反应发生的难易程度，计算公式如式(3)所示。对于任何反应，ΔG＜0时，可认为反应可以自发进行，而ΔG＞0时，反应不可自发进行。因此，本文基于反应活化能和吉布斯自由能变最小原则，来优选反应路径和判别反应发生的难易程度。

式中，GIM为过渡态物或生成物的吉布斯自由能，kJ/mol；GIS为反应物的吉布斯自由能，kJ/mol。

根据Eyring提出的化学反应的过渡态理论[31]，任何温度下对应的反应活化能（Ea,T，kJ/mol）的表达式如式(4)所示。

式中，ΔHθ为在298.15 K和101.33 kPa下基元过渡态物与基元反应物之间的焓变，kJ/mol；ΔGθ为在298.15 K和101.33 kPa下基元过渡态物与基元反应物之间的吉布斯活化自由能变，kJ/mol；ΔSθ为在298.15 K和101.33 kPa下的基元过渡态物与基元反应物之间的熵变，kJ/(mol·K)；n为基元反应物的分子个数；R为气体常数，8.31441 J/(mol·K)；T为温度，K。

由式(4)和式(5)变形可得出任何温度下的反应活化能与吉布斯自由能变以及熵变之间的关系式，如式(6)所示。

对于任何温度下的反应速率也可以根据过渡态理论进行计算[32]，如式(7)所示。

式中，k(T)为在任意温度下的反应速率，s-1；kB为玻尔兹曼常数，1.3806620×10-23J/K；h为普朗克常数，6.626176×10-34J·s；cθ为质量分数，取值100%。

2 结果与讨论

2.1 反应位置

为了辨别气化过程中反应发生的位置，需要找到反应物中的活性位点，从而为后续反应路径的寻找提供依据。活性位点是反应物最脆弱的反应位点，该处会最先发生反应。通过对优化后的反应物结构进行静电势分布分析，比较反应物中每个原子自带的电荷数，发现乙基苯与CO2反应物的活性位点，即可得到一个正电荷从无穷远处的位置移动到反应物上的某个位置所需要做的功的大小，正值（蓝色）代表做的功为正功，表明不易移动到该位置处；负值（红色）代表做的功为负功，表明容易移动到该位置处。由图2可知，CO2中的两端氧原子的周围以及煤的模型化合物中的苯环上的碳原子周围静电势能显示为负值，说明一个正电荷较为容易移动到该位置处。

图2 反应物的静电势分布Fig.2 Electrostatic potential distribution of reactants

对优化后的反应物结构进行了电荷分布分析，如表1所示。煤的模型化合物乙基苯中的氢原子以及CO2中的碳原子为正电荷，而煤的模型化合物中的大多数碳原子以及CO2中的氧原子为负电荷。由于C-19的正电荷的数值最大，在气化反应过程中，CO2中的碳原子优先与煤的模型化合物上的原子发生反应，结合图2可知，煤的模型化合物中的苯环结构周围的静电势能为负和侧链部分的静电势能数值略大于0，也就是说，这些位置可能被C-19结合。

表1 反应物中原子电荷分布Table 1 Atom charge distribution of reactants

2.2 反应路径

为了明确煤气化反应的具体过程，不仅需要明确反应位置，还需要探究反应路径。在上述分析的反应物的反应位置的基础上，对煤的模型化合物的可能会被结合的位置分别进行了尝试，列出了相应的反应路径，如图3所示。图3展示了CO2中的碳原子结合甲基、亚甲基、苯环上与侧链的相对位置为邻位、对位以及间位上的氢原子的各个反应变化过程，以及生成物生成情况。在不同的反应路径的基础上，分别对其进行了量子化学计算，如表2所示。

图3 5种可能的反应路径Fig.3 Five possible reaction pathways

由表2可知，各反应路径中对应步骤最大的反 应 活 化 能 依 次 为435.56 kJ/mol、373.23 kJ/mol、398.65 kJ/mol、403.03 kJ/mol和398.46 kJ/mol，最大的基元反应的过渡态物与反应物之间的吉布斯自由能变依次为449.21 kJ/mol、367.92 kJ/mol、417.96 kJ/mol、423.35 kJ/mol和424.57 kJ/mol。根据最小反应活化能和吉布斯自由能变原则，乙基苯与CO2反应过程中的反应路径优选为路径2，其反应过程为C6H5CH2CH3+CO2→C6H5CH(OCHO)CH3→C6H5CH(OH)CH3+CO→C6H5CHO+CO+CH4→C6H6+2CO+CH4。

表2 不同反应路径的量子化学计算结果Table 2 Quantum chemical calculation results of different reaction paths

2.3 无金属离子的反应特性

在优选出的反应路径2中，涉及到的甲烷生成的基元反应为1-苯乙醇裂解反应（C6H5CH(OH)CH3→C6H5CHO+CH4）。由图4可知，在无金属离子存在的情况下，1-苯乙醇裂解反应的活化能为325.96 kJ/mol，过渡态物与反应物的吉布斯自由能变为325.23 kJ/mol，反应物1-苯乙醇的单点能和吉布斯自由能分别为-1.01398×106kJ/mol和-1.01407×106kJ/mol。乙基苯在气化过程中形成的中间产物1-苯乙醇为甲烷生成反应的反应物，其反应过程是其侧链上的羟基结构上的氧氢键断裂，生成氢自由基，氢自由基与尾链的甲基结合，从而形成甲烷分子和苯甲醛。通过计算得出的甲烷生成反应的生成物与反应物之间的单点能变和吉布斯自由能变分别为-9.83 kJ/mol和-35.63 kJ/mol，而吉布斯自由能变小于0，说明在煤气化过程中，甲烷生成属于一种自发反应，但前提是必须克服转化为过渡态物的能量障碍。

图4 甲烷生成反应的分子结构和热力学参数Fig.4 Molecular structures and thermodynamic parameters of methane formation reaction

2.4 催化反应特性

为了促进甲烷生成反应发生，特别地在无金属离子存在情况下的反应物模型结构的基础上分别加入了Na+和Ca2+，讨论了金属离子对甲烷生成反应的催化作用。

2.4.1 Na+的影响

在上述煤气化模型基础上，加入Na+，探究了其对甲烷生成反应的影响。由图5可知，在Na+存在的情况下，甲烷生成的基元反应的活化能为310.84 kJ/mol，过渡态物与反应物的吉布斯自由能变为308.05 kJ/mol，反应物1-苯乙醇的单点能和吉布 斯 自由能 分 别为-1.43999×106kJ/mol和-1.44008×106kJ/mol。与无金属离子存在的情况相比，Na+的加入一方面使反应物结构的单点能和吉布斯自由能分别降低了42.02%和42.01%；另一方面，甲烷生成反应的反应活化能降低了4.64%，同时过渡态物与反应物的吉布斯自由能变降低了5.28%，说明Na+的加入有利于甲烷生成反应。Na+通过与苯环侧链上的氧原子成键，使反应过程中与氧原子相连的氢原子更加容易脱离，成为氢自由基，随后与甲基形成甲烷分子。

图5 Na+作用下甲烷生成反应的分子结构和热力学参数Fig.5 Molecular structures and thermodynamic parameters of methane formation reaction with Na+

2.4.2 Ca2+的影响

在无金属离子的煤气化模型上，加入Ca2+，探究了Ca2+对甲烷生成反应的影响。由图6可知，在Ca2+存在的情况下，甲烷生成的基元反应的活化能为296.25 kJ/mol，过渡态物与反应物的吉布斯自由能变为294.23 kJ/mol，反应物1-苯乙醇的单点能和吉布斯自由能分别为-2.79274×106kJ/mol和-2.79283×106kJ/mol。与无金属离子存在的情况相比，Ca2+的加入不仅使反应物结构的单点能和吉布斯自由能分别降低了175.42%和175.41%，且甲烷生成反应的反应活化能降低了9.12%，过渡态物与反应物的吉布斯自由能变降低了9.53%，说明Ca2+的加入也有利于甲烷生成反应。Ca2+同样与苯环侧链上的氧原子形成化学键，从而促进了与氧原子相连的氢原子脱离，向着苯环侧链末端处的甲基位置方向移动，最后形成甲烷分子。

图6 Ca2+作用下甲烷生成反应的分子结构和热力学参数Fig.6 Molecular structures and thermodynamic parameters of methane formation reaction with Ca2+

与Na+存在的情况相比，Ca2+的加入一方面降低了反应物结构的单点能93.94%和吉布斯自由能93.93%，另一方面降低了甲烷生成反应的反应活化能4.69%和过渡态物与反应物的吉布斯自由能变4.49%，这是由于Ca2+带有两个正电荷，比Na+自带的正电荷多，从而会影响到整体反应物结构上的静电势和原子电荷数分布情况，说明相比于Na+，Ca2+更加有利于促进甲烷生成反应的进行，这与GAO等[33]的实验结果趋势一致。

2.5 热力学特性

为了评估上述有无金属离子的甲烷生成反应在不同温度下的活化能变化，根据式(6)进行相关计算。在地下煤气化的现场试验中发现，核心区的温度高达2073.15 K[34]，因此，对温度在298.15～2073.15 K内进行活化能的计算。由图7可知，在所选定的温度区间内，无金属离子的甲烷生成反应的活化能始终大于Na+和Ca2+存在情况下的活化能，并且Ca2+存在的情况下的甲烷生成反应的活化能最低。依据无金属离子和Na+存在情况下活化能随温度变化的趋势来看，温度高于2073.15 K时，Na+将抑制甲烷生成反应的发生，而无金属离子和Ca2+存在的情况下活化能变化趋势相近，说明Ca2+对于甲烷生成反应的促进作用较为稳定，不易随着温度的升高而发生变化。

图7 甲烷生成反应的热力学参数Fig.7 Thermodynamic parameters of methane formation reaction

2.6 动力学特性

动力学参数分析可以反映反应过程中的反应速率，为评估反应进程快慢提供理论支撑。根据过渡态理论，基元反应中的活化能常常用于反应速率的计算，因此根据式(7)，仍选用上述相同的温度区间进行反应速率的计算。由图8可知，在温度不断升高的情况下，Ca2+存在情况下的反应速率始终最大，而无金属离子的存在情况下的反应速率始终较小。温度低于1266.27 K时，无金属离子存在情况下的反应速率较小；温度低于1197.22 K时，Na+存在情况下的反应速率较小；温度低于1137.93 K时，Ca2+存在情况下的反应速率较小。对于在无金属离子或金属离子（Na+、Ca2+）存在下，1-苯乙醇裂解生成甲烷的3个反应来说，在2073.15 K对应的反应速率大约是298.15 K的1044倍，说明温度对反应速率的快慢有极强的影响作用。在2.5节中所提及的在高于2073.15 K的温度下，Na+虽然会抑制甲烷生成反应的发生，但并不会影响到Na+和无金属离子存在情况下的反应速率的相对大小关系，说明较高温度可以使活化能发生变化，从而改变反应发生的难易程度，但对反应速率来说，基本没有受到影响。

图8 甲烷生成反应的动力学参数Fig.8 Kinetic parameters of methane formation reaction

3 结论

本文基于密度泛函理论，选择乙基苯作为煤的模型化合物，对煤催化气化过程进行了量子化学计算，并优选了乙基苯在CO2氛围下的反应路径。在无金属离子模型的基础上，探究了Na+和Ca2+对甲烷生成反应的影响，考察了温度对甲烷生成反应的活化能和速率的作用，得出了以下结论。

（1）以乙基苯和CO2作为反应物，优选出的反应路径中甲烷生成的基元反应为1-苯乙醇裂解反应（C6H5CH(OH)CH3→C6H5CHO+CH4），其 反 应活化能为325.96 kJ/mol。

（2）Na+和Ca2+均对甲烷生成反应起到了促进作用，其反应活化能与无催化剂相比，分别降低了4.64%和9.12%，且Ca2+的催化效果优于Na+。

（3）无论是否存在金属离子的催化，甲烷生成反应的活化能、速率分别随着温度的逐渐升高而增大、加快。