Ⅰ型CH4/CO2水合物稳定性的分子动力学模拟

宋 军，徐 梁，王蕊蕊，赵伟龙，赵晨旭

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

在永久冻土区和深海海底蕴含着大量的天然气水合物。如我国南海地区，天然气水合物资源丰富。合理地开发利用天然气水合物资源，将成为解决我国未来能源短缺问题的重要途径[1]。水合物包含3种常见的结构：Ⅰ型体心立方结构、Ⅱ型面心立方结构和H型六方结构[2]。1 m3天然气水合物可能含有多达180 m3的气体（标准温度和压力），水合物中所蕴含的碳含量是全球化石燃料总碳量的两倍[3]。

采用CO2置换CH4水合物的方法既能收集CH4，又能储存CO2，近年来被广泛研究。余静薇等[4]研究发现Ar作为小分子气体添加剂促进了CO2置换CH4水合物，提高了CH4的置换率，在置换反应的初始阶段CH4的置换率增长较快，在置换反应中后期CH4的置换率变化平缓。孙建业等[5]通过实验发现置换反应包含快速反应和缓慢反应两个阶段，表层水合物分解过程影响快速反应阶段的CH4产气，孔隙水的扩散过程影响快速反应阶段的CH4产气。宋光春等[6]研究了温度、压力对置换的影响，发现温度对置换结果的影响大于压力对置换结果的影响。

分子动力学模拟（MD）基于分子动力学理论，其通过力场得到粒子运动轨迹、系统动态和热力学统计数据以及分子结构性质。分子模拟计算精度高，与实验数据往往高度吻合并且能获得实验无法获得的数据，被广泛应用于多种学科的研究。一些学者利用分子动力学模拟方法从微观角度对水合物结构特性进行了研究。张庆东等[7]发现温度对水合物结构稳定性的影响较大，其中温度对大孔穴结构稳定性的影响较大，而对小孔穴结构稳定性的影响较小。Kusalik等[8]通过MD方法研究了水合物生长速率，发现水合物最大生长速率高于冰生长速率的4倍，且Ⅰ型水合物界面支持Ⅱ型水合物的生长。Wu等[9]对CH4、四氢呋喃（THF）和CH4+THF水合物的成核过程进行了分子动力学模拟，结果发现CH4+THF混合客体体系的诱导时间明显短于纯CH4和纯THF体系。Sun等[10]利用MD方法，研究了水合物大小笼空缺对结构稳定性的影响，发现大笼缺失客体分子比小笼缺失客体分子对结构稳定性的影响更大。芦文浩等[11]利用MD方法，研究了273 K、283 K和293 K，压力2 MPa条件下环状化合物作为添加剂对CH4水合物的影响，发现同等条件下添加环戊烷体系最稳定。贾菊等[12]利用MD方法，研究了环戊烷、THF和四氢吡喃对煤层气水合物结构的影响，结果发现263 K、10 MPa时，环戊烷可以很大程度地提高水合物结构的稳定性。易莉芝等[13]利用MD方法，研究了NaCl溶液中CO2置换CH4水合物的过程，发现在NaCl溶液中被置换出的CH4分子易形成气泡，溶液相和水合物相形成CH4浓度梯度，水合物发生分解。

由于宏观实验无法实现CH4水合物和CO2水合物在不同晶穴占有率条件下的对比分析，同时模拟实验也较少研究晶穴占有率对置换过程中CO2水合物和CH4水合物结构稳定性的影响，为更深入地了解CO2置换CH4水合物的机理，本文应用MD方法，采用蒙特卡洛吸附模块分别建立CH4水合物模型和CO2水合物模型，在一定温度压力条件下，分析不同晶穴占有率对CH4水合物和CO2水合物结构稳定性的影响。通过结构变化、径向分布函数、均方位移和势能变化等方面，从分子角度分析两者的变化，为水合物置换技术应用提供参考。

1 模型的建立及模拟方法

1.1 水合物模型

Ⅰ型水合物模型属于体心立方结构，建立超晶胞2 × 2 × 2模型，晶格参数2.406 nm × 2.406 nm×2.406 nm，空笼结构由368个水分子组成。笼型结构水分子中氧原子的初始坐标来源于X射线单晶衍射实验[14]，利用Bernal-Fowler规则[15]对氢原子进行排列。水合物空笼的初始模型如图1所示，吸附客体分子采用蒙特卡洛方法吸附，100%占有率笼型吸附位点如图2中红色标记所示，客体分子主要是CH4和CO2。用此方法分别建立不同晶穴占有率（12.5%、25.0%、37.5%、50.0%、62.5%、75.0%、87.5%和100%）的CH4水合物与CO2水合物。

图1 水合物空笼初始构型

图2 笼中晶穴吸附位点

1.2 模拟方法

分子动力学计算应用牛顿运动定律，获取各时间下粒子运动状态的物理参量，通过位移、速度、加速度的变化情况来了解分子运动的微观细节[16]。模型选用CVFF力场进行计算，系综选择NPT，水分子模型选用SPC模型，描述长程静电作用力和范德瓦尔斯作用力采用Ewald方法，Lennard-Jones作用势描述分子间作用力，截断半径1.5 nm，加和法计算精度4.19 × 10-6kJ/mol。

模拟之前，先对模型进行能量优化和结构优化，而后将笼型结构中氧原子坐标进行固定，使用NVT系综，在253 K条件下进行50 ps的驰豫，控温方法使用Nose-Hoover法。为了减小计算过程中的误差，设置驰豫的时间步长为1 fs，每间隔1000步输出一张结构图。驰豫结束后解除氧原子固定，将晶穴占有率不同的模型在NPT系综中进行分子动力学模拟，依然设置温度253 K、压力2 MPa、模拟时间400 ps，温控方式和结构图输出间隔不变。

2 结果与讨论

2.1 构象分析

将模拟之后的结构图进行汇总，选取晶穴占有率为37.5%、50.0%和62.5%的CH4水合物和CO2水合物进行对比，如图3所示。初始的水合物笼型结构具有高度对称性，因此结构图的部分氧原子和氢键被完全遮挡，根据氧原子、碳原子的位移偏离程度可判断笼型结构的变形程度。由图3可知，晶穴占有率为50.0%和62.5%的CO2水合物的氢键和氧原子排序重叠度较高，客体分子仍处于笼中心位置。晶穴占有率为37.5%的CO2水合物氧原子和氢键位置发生轻微偏移，但整体结构相对稳定。晶穴占有率为62.5%的CH4水合物氢键和氧原子发生了大的位移偏离，笼型结构发生破裂。晶穴占有率为37.5%、50.0%的CH4水合物相比晶穴占有率为62.5%的CH4水合物，氧原子出现更大的位移偏离，氢键的排序混乱，水合物笼型结构完全解离。

图3 不同晶穴占有率的CO2水合物（(a)～(c)）和CH4水合物（(d)～(f)）最终构象

2.2 径向分布函数

径向分布函数（RDF）是分析结构无序性的重要方法。RDF中gα-β表示在粒子α相距r处粒子β出现的概率，其也能利用区域密度与平均密度的比值表示。对于液态水分子，其氧原子的径向分布逐渐趋近于1，表示水分子处于无序状态。如果水合物结构处于稳定状态，则水分子的氧原子径向分布曲线呈现多个峰值，并且这些峰的峰高随着r值的增加而逐渐减小。第二特征峰消失或者峰值降低，说明水合物处于不稳定状态。

图4是不同晶穴占有率下CH4水合物和CO2水合物的RDF曲线。由图4可知，RDF曲线在r=0.278 nm附近出现第一特征峰，代表笼型结构中两个氧原子之间的最近距离，在r=0.452 nm附近出现第二特征峰，代表氧原子之间相隔一个或者多个原子的距离，满足稳定结构的RDF。随着晶穴占有率增加，峰值升高、峰谷降低，说明水分子无序性降低，水合物稳定性增加。对比不同晶穴占有率的CO2水合物，RDF特征峰区别不明显，晶穴占有率为62.5%的体系结构稳定性优于37.5%和50.0%的体系。相同工况下，对比不同晶穴占有率的CH4水合物的径向分布，CH4分子晶穴占有率越小第一特征峰越低，第二特征峰也越低或消失，水合物结构无序程度增加。在相同工况下，晶穴占有率在37.5%与62.5%之间时，CO2水合物的RDF特征峰均明显优于CH4水合物。

图4 不同体系水合物氧原子RDF曲线

2.3 均方位移

在分子模拟计算过程中，原子在不停运动，其位置在不断变化，会相对于原始位置偏离一定距离。均方位移（MSD）用公式表示如下：

式中，N为系统中粒子总数；Ri(0)为i粒子0时刻的位置；Ri(t)为i粒子t时刻的位置。

MSD可以很清楚地反映粒子相对于初始位置的偏离程度。当体系为液态或者气态时，粒子所受的约束变小，此时的MSD曲线是一条随时间增长而上升的曲线，对应于水合物分解状态。稳定的水合物结构，粒子受到的约束大，MSD值会在0值附近较小波动，曲线整体接近一条水平线。

图5为CH4水合物不同晶穴占有率时氧原子的MSD曲线。由图5可知，晶穴占有率为75.0%、87.5%和100%的CH4水合物氧原子MSD值均很小，整条曲线接近水平，氧原子偏离程度很小，说明水合物笼型结构是稳定的。CH4水合物晶穴占有率低于75.0%后，占有率越低MSD值随时间增长的速率越快，水分子的无序程度越大，笼型结构解离，不能稳定存在。

图5 CH4水合物氧原子MSD曲线

图6为相同温压时晶穴占有率为37.5%、50.0%和62.5%的CH4水合物和CO2水合物氧原子的MSD曲线。由图6可知，相对于CO2水合物，CH4水合物在晶穴占有率为62.5%时的MSD值较小，但仍然高于同种占有率的CO2水合物，说明此时CH4水合物笼型结构无序化程度高于CO2水合物。占有率为37.5%、50%和62.5%时，CO2水合物的氧原子MSD值较为稳定，围绕0值附近波动，笼型结构未发生明显扭曲。因此，在一定温度、压力和晶穴占有率范围内，CO2可以形成稳定的CO2水合物，而CH4却不易形成稳定的CH4水合物，CH4水合物分解释放CH4而CO2仍然可以稳定存在于水合物中，从而达到CO2置换CH4的目的。

图6 不同体系水合物氧原子MSD曲线

2.4 势能变化

在分子模拟过程中，体系的势能有波动性，对于稳定的水合物结构，势能在一定范围平稳波动。图7为晶穴占有率为37.5%、50.0%和62.5%的CH4水合物和CO2水合物在模拟过程中势能随时间变化的曲线。由图7可知，晶穴占有率为37.5%的CH4水合物在140 ps附近势能从-15200 kJ/mol迅速升高到-13600 kJ/mol，水合物结构发生了快速分解，水分子的无序程度迅速升高。晶穴占有率为50.0%和62.5%的CH4水合物分别在180 ps附近和300 ps附近势能升高，结构同样也发生了解离。CH4水合物晶穴占有率越低，结构分解的时间越早，分解的速率也越快。而相对于相同占有率的CH4水合物，CO2水合物的势能在-15200 kJ/mol附近较小波动，表明体系整体结构较为稳定。

图7 不同体系水合物势能随时间的变化曲线

图8显示了CO2水合物势能随时间的变化。由图8可知，CO2水合物晶穴占有率低于37.5%时，势能在不同时间都升高，结构发生了解离；晶穴占有率在37.5%以上时势能波动小，结构较为稳定。

图8 CO2水合物势能随时间的变化曲线

3 结论

利用MD及蒙特卡洛吸附法建立了CH4水合物和CO2水合物模型，在NPT系综下研究了其稳定性，分析了最终构象、径向分布函数、均方位移以及势能的变化，得出以下结论：

（1）在253 K、2 MPa条件下，Ⅰ型CO2水合物晶穴占有率在37.5%以上时，水分子氢键组成的笼型结构整体能保持稳定，并且CO2的占有率越高水合物笼型结构越稳定。

（2）相同晶穴占有率时，CO2水合物比CH4水合物结构稳定性更强。晶穴占有率在37.5%与62.5%之间时，CH4水合物均发生了不同程度的分解并且晶穴占有率越低越先发生分解。此时，CO2水合物结构虽有不同程度的轻微扭曲，但整体结构依然保持稳定。

（3）在一定温度压力和晶穴占有率范围内，CO2水合物比CH4水合物稳定性好，解释了CO2能置换CH4水合物的原因，即CH4水合物发生分解释放CH4，而CO2稳定存在于水合物中，从而实现收集CH4并储存CO2的目的。