王蕊蕊，杨保亚，周 航，赵伟龙，宋 军，盛 伟

（河南理工大学 机械与动力工程学院，河南 焦作 454003）

甲烷水合物是甲烷（CH4）气体和液态水在特定的环境下，通过足够的压力驱动力和较高的过冷度生成的外观似冰包络型化合物。其中水分子利用氢键彼此结合，组成不同的空腔结构，CH4分子被包裹在这些空腔之中，主体空腔与气体分子之间受到范德华力的作用[1]。在标准状况下，1 m3水合物可储存160~180 m3CH4，储气密度极高，并且CH4在燃烧时热值高，污染少，因此甲烷水合物可作为未来的新兴能源之一。然而生成速率低、生成过程中诱导时间长和生成条件苛刻等问题，严重制约了甲烷水合物的工业化应用。

为了解决以上问题，已经有了相应的解决对策，如使用热力学促进剂四氢呋喃（THF）[2]、四丁基氯化铵（TBAB）[3]和四丁基溴化铵（TBAC）[4]等，这些添加剂分子通过占据水合物中较大的空腔，或者与水分子结合成半笼的方式，使其相平衡曲线朝着压力降低温度升高的方向移动，导致水合物在较缓和的工况下就能生成。动力学促进剂十二烷基硫酸钠（SDS）[5]和十二烷基苯磺酸钠（SDBS）[6]通过减小气液界面的表面张力，使更多的气体分子进入液相，增大了传质速率进而直接影响气体消耗量、气体吸收速率、诱导时间和水合物储存容量，大大提高了水合物生成速度和动力学特性。虽然传统的添加剂在促进水合物的生成方面已经取得了不错的效果，但会对环境造成一定危害。近些年来，高效且对环境无污染的纳米粒子已经成为水合物促进剂研究的新方向。氧化石墨烯（GO）作为一种新型的纳米粒子，在极性溶液中分散效果良好，并且保留了原本石墨烯中的炭网结构，导热效果优异[7]，可以快速的转移水合物生成过程中的热量，巨大的比表面积也为水合物提供了更多的成核位点。LIU等[8]研究了在静态条件下GO以及SDS和THF混合添加剂对二氧化碳水合物形成的影响。结果发现，加入一定量GO有效缩短了二氧化碳水合物形成诱导时间，相较于纯水缩短了74%。同时，在含有SDS/THF/GO混合添加剂的体系中，水合物发生二次成核，很大程度上增强了水合物的储气能力，最高可达到152%。此外，GO也被广泛用于CH4[9]、乙烯[10]和丙烷[11]水合物的生成中。

目前关于GO促进水合物的研究多局限于宏观的实验现象方面，对于从微观尺度分析其对水合物生成的研究较少。分子动力学（MD）模拟是从微观层面研究系统体系的一种方式，利用牛顿经典力学计算分子在不同时刻的相空间轨迹，模拟出系统随时间演变的过程。实验上达不到的条件可以通过模拟表示其物理现象和过程，并且可以从分子层面探究作用机理，已经有不少学者将MD模拟应用到水合物的研究中，LIAO等[12]模拟了蜡分子对甲烷水合物生长阶段的作用，发现蜡分子的加入给水合物提供了更多的成核位点，从而促进其生成。刘庭崧等[13]利用MD模拟的方法研究了不同因素（温度、压力和电解质溶液）对甲烷水合物生长速率的影响，发现在一定压力范围内，提高压力有助于水合物生长，而一定浓度的电解质溶液反而会抑制水合物生成。但是在已有水合物模拟的研究中，对添加剂分子的建模以及添加剂存在时整个模拟体系的搭建，尚存在不足。

本文采用MD模拟的方法，建立了在不同时期含有GO的模拟体系，从最终构象、均方位移和结合能的角度分析了GO对甲烷水合物的影响，为GO用做甲烷水合物生成的高效促进剂提供理论依据。

1 模型构建与参数设定

1.1 GO模型的构建

GO由氧化石墨剥离得到，主基面是单层的炭网结构，并负载了大量的环氧基团，在边缘处又引入了羧基以及羟基。由于制备过程以及实验条件等方面差异，迄今为止GO的精确结构还无法确定，根据Lerf-Klinowski结构（图1(a)）[14]与理论结构（图1(b)），使用MS 8.0 软件建立两种GO的模型，联合 Forcite模块，采用最速下降法和共轭梯度法进行模型的结构优化。优化后的模型不再是平整的单层结构，模型表面发生波浪形的卷曲，如图1(c)和图1(d)所示。

图1 GO结构和模型示意Fig. 1 Structures and models of GO

图中黄色虚线代表GO所负载的羟基羧基等基团，这些含氧基团能够与水分子结合形成O─H…O氢键，用于组成水合物笼，GO的加入能够与自由水结合提供更多的氢键，加速水合物笼的形成，目前在实验研究[8,15-16]中也已经证明GO能够大大缩短水合物的生成时间。由图1 可知，两种模型都包含相同的结构，可根据所需构建模拟体系的大小，选择合理的结构模型，便于简化计算。此外，两种GO模型都具有较大的比表面积，在实际实验过程中能够提供更多的成核位点[15]，能够促进水合物成核。

1.2 诱导生成时期模拟体系的构建

由于诱导生成时期的模型包含大量CH4和水分子，即构建的模拟体系尺寸较大，因此采用如图1(c)所示的GO模型建立了诱导生成时期CH4-GO-H2O模拟体系，如图2(a)所示。为了对照分析，另建立了未添加GO时的模型，如图2(b)所示。

图2 诱导生成时期添加GO和未添加GO的模拟体系Fig. 2 System models with GO and without GO in induction period

在模拟之前，需对整个体系进行优化，所采取的方法同上文1.1 节，最终得到体系最优的结构。整个模拟体系在COMPASS力场中进行计算，分子间的相互作用利用Lennard-Jones势函数表示。体系内分子存在的静电力采用Ewald求和法进行计算，Berendsen法控制体系压力[17]，Nose热溶法控制体系温度[18]。体系内所有分子间相互作用在1.5 nm的球形半径截止处被截断。完成以上工作后，为使整个体系达到平衡，选择NVT系综，在温度为278 K条件下先进行100 ps的弛豫计算。设置合理的时间步长，另规定运行1000 步输出1 次图片。之后在NVT系综下再进行600 ps的动力学计算，温度依旧设置为278 K ，模拟结果用于统计和机理分析。

1.3 快速生成时期模拟体系的构建

GO对甲烷水合物快速生成过程的影响，可以通过构建水合物切面与GO结合的模型来表示。使用MS软件中的Build模块，对结构优化后的I型甲烷水合物模型进行切割处理，目的是通过暴露不同的切面，来反应甲烷水合物的快速生成过程。由于水合物切面尺寸较小，为使模型构建的更加合理也为了节约时间减少运算量，使用图1(d)GO模型，借助Layer工具搭建模拟体系，如图3 所示。切面1 表示甲烷水合物已经开始生长但结构尚未完整；切面2表示水合物生长进行了一段时间，结构相对完整；切面3 表示水合物开始大量生成进入循环生成阶段。模拟步骤以及模拟参数的选择同上文1.1 节和1.2 节。

图3 甲烷水合物快速生长时期模拟体系Fig. 3 Simulation system of methane hydrate rapid growth period

2 结果与讨论

2.1 平衡态的判定

为减少计算误差，通常在正式MD模拟之前，要对弛豫后的模型进行平衡状态的判定，包括能量平衡和温度平衡。诱导生成时期弛豫过程中模拟体系的能量和温度随时间变化的曲线如图4 和图5所示。在初始阶段，整个体系的动能、势能、非键能以及总能量波动较大，模拟过程运行14 ps左右，所有能量曲线在某一固定值附近上下变化且变化较小，此时体系的能量已经达到平衡。体系的温度变化和能量变化趋势相似，并且最终的温度变化偏差低于10 K，表明体系已处于平衡状态。

图4 弛豫过程能量变化曲线Fig. 4 Energy curve of relaxation process

图5 弛豫过程温度变化曲线Fig. 5 Temperature curve of relaxation process

2.2 诱导生成时期构象及均方位移分析

水合物的形成包括气体溶解、诱导生成和快速生成阶段，目前本文模拟方法不能表示气体溶解的过程，主要通过添加剂对诱导生成时期和快速生成时期晶体构象的影响来分析GO对水合物生成的促进作用。诱导生成时期模拟体系的最终构象如图6所示，当模拟计算50 ps左右，气体分子朝着溶液层的方向扩散剧烈程度开始增大，直至计算结束大量的CH4分子进入到溶液中；另外通过对比发现水分子与GO表面基团结合，形成较强的氢键键能导致GO分子发生了一定程度的变形。

图6 诱导生成时期模拟体系最终构象Fig. 6 Final conformations of simulated system in induction period

同时，氢键会对溶液层的水分子产生“约束力”，影响水分子的扩散能力，使大量的自由水转化为结合水，然而从图6(a)并不能直观的看出这一过程。在动力学计算中可以通过分析均方位移（MSD，m2）[19]来表示水分子扩散的剧烈程度，其定义式如式(1)所示。

式中，Ri(0)和Ri(t)为分子i在初始时刻和不同时刻的位置；N为系统分子总数。

诱导生成时期甲烷水合物模拟体系的水分子MSD如图7 所示。

图7 不同甲烷水合物体系水分子的MSDFig. 7 MSD of water molecules in different methane hydrate systems

由图7 可知，相同模拟时长和工况条件下，未添加GO的模拟体系，水分子的MSD更大，约为添加GO时体系的2 倍，表明此时体系中自由水的含量较多，水分子运动也更加剧烈。由于水合物笼是依靠水分子之间氢键的彼此约束形成的，当水分子运动越剧烈，水合物笼越不易形成，进而影响整个水合物的生成速率，而GO通过表面基团与水分子结合成键的方式，在一定程度上约束了水分子运动，同时结合而成的大量氢键也更加有利于水合物笼的形成，可大大提高生成速率。

2.3 快速生成时期构象及结合能分析

快速生成时期模拟体系的最终构象如图8 所示。

图8 快速生成时期甲烷水合物不同切面模拟体系的最终模拟构象Fig. 8 Final conformations of methane hydrate in different sections during rapid growth

由图8 可知，3 个不同水合物切面下GO都发生了一定程度的变形，由于切面1 暴露了半笼型的结构，在氢键的吸引下，使得GO表面的边缘基团逐渐向表面贴近。而切面2 表示的是水合物快速生成时期已经进行了一段时间，结构相对完整，导致GO变形程度较小。对于切面3 来说，已经经历一个完整的水合物生成过程，并且即将进入循环生成阶段，此过程需要再次进行成核成笼。此时GO分子发生较大的形变，原因是GO能够为气体分子提供大量的成核位点，甲烷水合物倾向于在GO表面成核，此外水分子和含氧基团互相吸引形成氢键使得GO更加贴近于水合物切面。

从图8 所示的甲烷水合物晶体构象不能看出加入GO之后不同切面结合模型的结合能大小，需要通过进一步的计算来确定，其详细计算公式如式(2)和式(3)所示。

式中，Einteraction为相互作用能，kJ/mol；Etotal为模型总能量，kJ/mol；Ehdyrate为水合物切面的能量，kJ/mol；EGO为GO的能量，kJ/mol；Ebinding为结合能，kJ/mol。

不同甲烷水合物晶体切面模拟体系的相互作用能示意如图9 所示，作为对比，分析了未添加GO时不同切面的能量。

图9 不同甲烷水合物晶体切面模拟体系的相互作用能Fig. 9 Interaction energy of different methane hydrate crystal sections simulation systems

由图9 可知，添加GO大幅提升了不同甲烷水合物切面的相互作用能，意味着GO的加入增强了水合物之间所受力的作用，影响了水合物的结构进而影响了水合物生成速率。由于结合能和相互做用能互为相反数，从图中数据可以得出，不同切面模拟体系的结合能均为正值，说明GO和不同切面之间是引力作用，有利于GO参杂到水合物的结构中，进而导致甲烷水合物在GO表面进行异相成核，降低了在溶液中进行均相成核的能垒，加快了生长速率。

3 结论

本文利用MD模拟方法，模拟了GO对甲烷水合物诱导生成时期以及快速生成时期的影响，得到以下结论。

（1）GO表面以及边缘负载的含氧基团与水分子结合形成氢键，有利于水合物笼的形成，提高了生成效率。

（2）结合诱导生成时期模拟体系的最终构象以及水分子MSD的变化，发现当加入GO之后会束缚水分子的运动，此时大量的自由水转化为结合水，更易形成水合物的笼型结构。

（3）GO对甲烷水合物快速生成时期也有一定的促进作用，添加GO增强了与水合物切面之间所受力的作用，并且二者之间的引力，利于甲烷水合物在GO表面进行异相成核，达到快速生成的目的。