环状化合物-甲烷水合物稳定性的分子模拟

芦文浩，梁海峰，王 帅，贾 菊

（太原理工大学 化学化工学院，山西 太原 030024）

天然气（NG）作为一种可提供高效热值、燃烧无污染和储量丰富的能源，被公认为未来重要的战略资源。以水合物形式存储天然气（主要成分为甲烷）的固化天然气技术 （Solidified natural gas，SNG），是最具有潜力、安全环保、高效解决能源供应问题的储运方式。

Ⅱ型水合物小笼晶穴可以容纳大量CH4，大晶穴填充热力学促进剂，稳定性提高，其中以环状化合物为主要的促进剂，可使水合物在温和条件下稳定存在。许多学者就实验探索已证明其可行性：胡亚飞等[1]构建环戊烷（CP）+CH4水合物实验体系，证明CP可提升水合物稳定存在温度；孙强等[2]添加环状醚类四氢呋喃（THF）提纯煤层气中的CH4，实验验证生成水合物压力明显降低；梁德青等[3]测定H2+CH4+环状醚类四氢吡喃（THP）+水（H2O）、H2+CH4+CP+H2O两种体系的水合物相平和数据，结果表明加入环状化合物，降低了生成压力；实验只能提供直观过程，无法洞悉微观结构改变和作用机理。水合物分子动力学模拟可为多组分体系中分析分子运动、作用和结构改变等微观角度的探究提供有效途径[4]，对此研究也取得一些成果：梅东海等[5]首次采用分子动力学计算H型水合物，客体分子填充CH4和环辛烷，得出H型水合物3种笼形晶穴都需要有客体分子占据，水合物才能稳定存在；Wu[6]利用GROMACS模拟软件，计算CH4+THF混合客体分子水合物成核的分子动力学，结果说明体系形成的笼形簇团相对单组份更稳定；张庆东等[7-8]采用分子动力学模拟填充不同占有率的CH4和THF对水合物稳定性影响，结果表明水合物占有率决定水合物稳定性，随着温度降低，利于单组份CH4水合物稳定；殷开梁等[9]比较 Material studio（MS）软件中各分子力场对于氢键作用描述的适用性，结果表明该软件可以采用一般的非键作用近似反映氢键。耿春宇等[10]通过分子动力学模拟CH4对应Ⅰ型水合物，得出结论：因CH4占据晶穴，水合物晶穴的相互作用能降低，水分子间氢键作用加强，稳定性提高，且占据率高于87.5%水合物才能表现较好稳定性；丁丽颖等[11]模拟填充CH4的Ⅰ型水合物稳定性受晶穴占有率和温度的影响，结果表明较低占有率（<0.625）和较高温度（>320k）条件下，水合物基本瓦解。

综上所述，有关环状化合物-CH4水合物稳定性研究多数集中在实验，少有利用分子动力学探寻CP和THP对Ⅱ型天然气水合物的稳定作用，且高占有率利于水合物稳定性；本文采用Material Studio 8.0软件，搭建占有率100%分别添加CP、THF和THP的二元体系CH4水合物，从分子角度比较其稳定性和解释作用机理，为SNG技术实际应用提供参考。

1 模型搭建和模拟方法

1.1 笼形单元组成

水合物主体框架是以水分子氢键结合形成的一种空间点阵结构[12]，Ⅱ型水合物原始单晶胞包含136个水分子，构成16个小笼晶穴（512）和8个大笼晶穴(51264)。小笼晶穴由CH4分子占据，大笼晶穴分别由CP、THF和THP占据。客体分子手动搭建，并进行能量几何优化，使其最接近真实状态。表1为依据胡盛志等[13]优化方法计算得到的客体分子范德华模型相关参数，笼形晶穴占据结构如图1所示。

表1 客体分子范德华模型参数

图1 客体分子占据笼形晶穴结构

1.2 模型搭建

在周期边界条件下，Ⅱ型水合物原始单晶胞在温度T=123K下，晶格参数a=b=c=1.7175nm。晶胞中心包含对称存在的四面体，与周边6个十二面体形成Ⅱ型水合物主体框架，其氧原子坐标来自X射线衍射实验[14]，氢原子分布满足“冰规则”[15]。客体分子与主体框架间为范德华作用力（Van der Waals）[16]，模拟体系构建2×2×2的Ⅱ型Super cell晶胞，该体系中1088个H2O构建128个小笼晶穴（512）和64个大笼晶穴(51264)的主体框架。小笼晶穴填充128个CH4，大笼晶穴分别填充64个CP、THF和THP，本文3种模拟体系分别命名为A、B、C，便于观察，将体系H原子隐藏，初始结构模型如图2所示。

图2 水合物模拟体系初始结构

1.3 模拟方法

先对模拟体系进行能量优化和几何优化，选用CVFF（一致性价力场）力场描述模拟体系，SPC（Simple point charge）模型[17]描述 H2O 分子，将其视作刚性分子，H-O-H平面夹角为104.52°，H-O键长为固定值0.09570nm，水分子和客体分子、客体分子间非键作用使用Lennard-Jones作用势描述，在MS中采用截断球近似代替，截断半径为1.5nm，Ewald加和法处理长程静电作用力，精度设置为4.1868×103J/mol。

对几何优化结束的体系进行结构松弛 ，选水分子中的O原子固定，控温方法使用Nose-Hoover，NVT（正则系综）模拟设定T分别为273K、283K、293K，分子初始速依据Maxwell-Boltzmann随机产生，模拟时间为 50ps（1ps=10-12s），步长为 1fs（1fs=10-15s），每5000步输出一帧结构图；解除O原子固定，在NPT（等温恒压）系综下进行动力学模拟，控温方式不变，采用Berendsen控压方式调整体系压力值，压力设为2MPa，模拟时间为400ps，模拟结束后分析二元体系水合物结构。

1.4 模拟最终结构

将软件输出的最终结构进行比较分析，如图3所示。

图3 水合物模拟体系最终结构

由图3水合物模拟体系最终结构看出：A、B可维系原有结构，氢键排序重叠度较高，客体分子仍能与笼形晶穴中心相对重合；压力不变，随着温度升高，出现部分晶穴的扭曲加深，B较A在293K下水合物结构破坏严重；C模型体系在273K时完整度与相同条件下的A、B相比较差，部分客体分子因晶穴的扭曲加重位移较大，氢键分布出现小范围的杂乱，随着温度的升高结构出现了不同程度的扭曲，283K温度下水合物结构基本破坏，残存的氢键不足以维系体系的稳定存在，客体分子位移较大，出现小范围的聚集，293K条件下氢键虽然存在，但晶穴有序结构不存在，客体分子逃逸，出现聚集现象，预判此时从固体转变为液体。

2 模拟结果分析

2.1 径向分布函数（Radial Distribution Function，RDF）

稳定的水合物主体框架水分子中的O原子RDF特征是：呈现多峰值曲线状态，其峰高随着O原子的增加逐渐降低。可依此判断其无序状态，也代表水合物区域密度和平均密度的比值；液态水分子中O原子go-o(r)在r值大于0.278nm时就趋向为1[18]，即区域密度等于平均密度表现为液态流动，将A、B、C三种模型对应压力温度值下的go-o(r)值对比分析，如图4所示。

图4 水合物O原子RDF分布

由图4可知O原子RDF分布在r=0.278nm时出现第一峰高，代表相邻水分子最近的O原子间距离，第二特征峰高在r=0.452nm出现，满足晶体结构RDF特征；T=273K时，A、B、C呈现出稳定水合物晶体的特征峰，但都出现了峰高降低和峰谷升高的现象，说明有部分O原子发生位移，即部分笼形晶穴扭曲，但整体框架结构能够稳定存在，且第二峰高峰值ABC，则A稳定性要高于相同条件下的B，C模型峰谷升高，稳定性较差；T=283K时A同样略优于B，C模型在r=0.278nm出现峰值后，峰谷升高，第二特征峰消失，此时与液体特征峰相似[15]；T=293K时，A和B模型RDF分布重合度相近，C模型RDF分布较T=283K对应分布：峰谷升高加剧，与液体特征峰分布更为接近，说明破坏程度更大，客体分子逃逸并发生聚集。

2.2 均方位移（MSD）

MSD表示指定分子位置与初始位置的偏离程度，依此可以判断体系是固态还是液态，MSD分布图像斜率K越大表示偏离初始程度越大。图5为A、B、C对应压力温度值下的MSD图像。

图5为模拟体系所有条件下的O原子MSD分布，为更清晰比较，去除C模型283K、293K条件，呈现图6。图6可分析出C模型在283K、293K条件下，MSD随着模拟时间增大，与液态水分子MSD斜率K相似，表明这两种状态下无序化程度严重，结构破坏变成了液态溶液；从图6、7综合比较，相同p值下，较稳定的模型体系MSD围绕在接近零值的地方波动，结果 KC(293K)KB(293K)KA(293K)，KC(283K)KB(283K)KA(283K)，KC(273K)KB(273K)KA(273K)，说明对添加同种环状化合物促进剂，随着温度降低，有利于水合物的稳定性，且体系在对应模拟条件下呈现KAKBKC；

图5 水合物O原子均方位移

图6 水合物O原子均方位移局部放大图

2.3 相互作用能

分子动力学模拟结束后可得到模拟体系的总能量Etotal，主体框架能量Ecages以及填充的所有客体分子能量Eguests，从而得到主体框架和客体分子间的相互作用能Ebinding，满足公式：

从能量角度分析，相互作用能越小，则体系最终结构越稳定，图7为A、B、C对应温度压力值下整个体系的相互作用能。

从图7可看出，A模型E(293K)E(283K)E(273K)，B、C模型呈现相同规律即低温有利于水合物结构稳定；相互作用能越低，水合物结构越稳定，A模型整体稳定性较好，与MSD分析结果相一致；从分子角度分析3种促进剂对稳定性的影响，已知客体分子与晶穴直径比值在0.76～1之间[19]，可维系稳定水合物结构，环状化合物的填充，起到支撑水合物框架的作用。从表1可看出Ⅱ型大晶穴，直径比值CP略高于THF，更接近1，稳定性得到提升，THP与晶穴比值优于CP和THF，理论上应相对更稳定，但模拟结果高温区稳定性较差，原因分析：THP分子过大，高温状态下分子活跃程度较大，其分子组成中的O原子会与周边水分子框架形成氢键作用，使得与周边其他水分子间范德华作用力不均衡，客体分子牵动晶穴单元发生位移（氢键作用强于范德华作用力），从而使得模拟体系稳定性下降。

图7 模拟体系不同温度下相互作用能

3 结论

在NPT系综下分子动力学模拟，本文通过最终结构、RDF、MSD以及相互作用能对比分析，得出以下结论：

（1）环状化合物作为客体分子占据甲烷水合物大晶穴，能在温和条件下维系水合物稳定存在。随着温度升高，模拟体系稳定性下降；

（2）CP模拟体系稳定性优于THF模拟体系，说明环状化合物与Ⅱ型水合物大晶穴的直径比越接近1，水合物越稳定；

（3）THF模拟体系稳定性优于THP模拟体系，说明在高温条件下，当较大分子尺寸的环状化合物分子含有O原子时，环状化合物与水合物框架产生不平衡氢键作用，使水合物稳定性降低。

符号说明

a-晶格参数，nm；b-晶格参数，nm；c-晶格参数，nm；Ebinding-相互作用能，J/mol；Ecages-空晶穴能量，J/mol；Eguests-总客体分子能 J/mol；Etotal-模型总能量，J/mol；E（T）-对应温度下相互作用能，J/mol；K-MSD分布斜率，m2/s；MSD—O 原子均方位移，m2；N-分子总数， 个；p-压力，MPa；r-截断半径，nm；RDF-径向分布函数，量纲为1；T-温度，K。