高压条件天然气水合物形成过程的分子动力学模拟

Molecular Dynamics Simulation of Gas Hydrate Formation Process under High Pressure Condition

摘要: 天然气水合物作为一种潜在的清洁能源已经引起广泛关注，为了更好地开发利用该资源必须对其形成机理进行充分了解和精准分析。本工作采用分子动力学方法旨在揭示1 GPa条件下天然气水合物的形成及结构演化。通过径向分布函数、均方根位移及环统计手段对研究体系的甲烷分子和水分子的行为及其周围环境进行分析。模拟结果表明，研究范围内温度对天然气水合物的结构影响不大，但对天然气水合物的生成数量具有明显的影响。环统计分析表明研究体系生成的五元环数量众多且相互联接而生成的六元环数量较少且孤立存在，两者数量比约为5:1。

Abstract: As a potential clean energy source, natural gas hydrate (NGH) has attracted wide attention. Sufficient understanding and accurate analyses for the complex formation mechanism in a natural gas hydrate system have to occur before the hydrate untilization. In this paper, molecular dynamics simulation is carried out to explore the formation and structure patterns under 1 GPa conditions. Radial distribution function, mean square displacement and ring statistics provide important insights into the behavior of the methane and water molecules present in the system. The simulations show that the structure of nature gas hydrate keeps unchanged at different temperatures, while the amount of gas hydrate is significantly influenced by temperature. Ring statistics were used to interpret that the five rings number is more and all connected, a small number of six rings for isolated form, the ratio is 5:1.

1. 引言

天然气水合物通常是指在一定温度和压强条件下甲烷分子与水形成的类冰状物质，其中，水分子作为主体分子通过氢键相互连接，形成大小不同的多面体笼子，这些多面体笼子通过点、线或面相互连接向三维空间扩展，形成水合物基本的笼形结构，甲烷分子则作为客体分子被包裹在笼中 [1]。每立方米的天然气水合物分解后可释放174立方米的甲烷气体，并留下0.8立方米的水，因此它是一种潜在的新型清洁能源 [2]。2017年5月18日，我国国土资源部宣布中国天然气水合物首次试采成功，标志着我国天然气水合物研究及应用已经达到新高度，这将改变我国的能源结构，改变世界能源格局 [3]。但是天然气水合物的基础研究依然薄弱，工业化应用技术也尚不成熟，达到商业应用及工业化使用必须解决天然气水合物的形成及分解机理问题 [4]。近几年，由于天然气水合物在能源、存储和分离气体等方面有很大的应用前景 [5] - [10]，关于天然气水合物的实验及理论研究因此出现了暴发性增长 [11]。虽然目前大部分文献针对天然气水合物的热力学形成条件进行了大量的研究，并且对于天然气水合物形成的热力学机理建立起了比较系统的基础理论数据库 [12]，但是，由于实验及计算条件的限制，人们对天然气水合物形成过程的动力学机理的研究一直以来前进缓慢，致使天然气水合物形成机理的研究十分有限 [13]。本工作利用分子动力学方法对天然气水合物形成过程在微观尺度上进行了研究，并利用环统计工具对天然气水合物演化过程中甲烷周围水分子的五元环和六元环进行统计，以此揭示天然气水合物在高压条件下的微观生成机理。

2. 计算方法及过程

分子动力学模拟采用的是基于牛顿力学的分子动力学软件LAMMPS [14]，在此之前LAMMPS已经广泛应用于天然气水合物的形成 [15] [16]、分解 [17] 及性质的研究 [15] [17] [18] [19]，并取得了可观的成果。环统计工具使用开源的R.I.N.G.S (Rigorous Investingation of Networks generated using simulations) [20] 代码，R.I.N.G.S工具已经广泛被用来分析晶体材料、玻璃材料等 [21] [22] [23] [24]，取得了不错的效果。

模拟过程中水分子的力场采用的是刚性力场TIP4P-2005 [25] 势，其键长和键角用SHAKE算法来固定，甲烷分子采用OPLS-UA [26] 力场，分子间的范德华力由LJ势指定，不同原子类型之间采用平滑粒子网格的Ewand求和，水分子及甲烷分子的势函数的详细参数列于表1中。构建的天然气水合物的模型为I型甲烷水合物的3 × 3 × 3的超胞，包含216个甲烷分子和1242个水分子，共计4806个原子，其结构如图1所示。整个模拟过程采用NPT系综，温度控制采用Nose-Hoover热浴，其热阻尼为0.1 ps，压力控制采用BAOROSTAT方法，其阻尼常数为1 ps。计算过程中使用周期性边界条件消除体系原子数量有限带来的误差，采用PPPM算法来处理长程库仑积分，时间步长为1 fs。

计算过程：1) 将系统温度设定为410 K，压强设为0.1 MPa，在此条件下平衡3 ns，使甲烷水合物分解；2) 保持压强不变，在2 ns的时间范围内，对系统降温，使温度达到310 K；3) 保持温度恒定，将系统压强在7 ns内逐渐升高并稳定在1 GPa；4) 保持压强恒定，将温度降至预定温度，分别为250、260、270、280、290、310 K，并平衡2 ns的时间；5) 最后，在恒温恒压条件下使系统弛豫4 ns，并选取最后2 ns的计算结果进行分析，水合物制备过程示意图如图2所示。

Table 1. Partial charge of CH4 and H2O molecules and atomic site-site Lennard-Jones potential parameters for CH4 Hydarte-H2O system

表1. 甲烷与水分子的部分电荷及甲烷水合物体系的原子位LJ势函数参数

红色小球——O原子；蓝色小球——水分子的H原子；黄色小球——C原子；白色小球——甲烷分子的H原子

Figure 1. Computational model of type I natural gas hydrate

图1. I型天然气水合物的计算模型

Figure 2. Illustration of the melting-quenching process for the preparation of hydrate

图2. 水合物制备过程示意图

3. 结果与讨论

3.1. 势能

为了表征分子动力学计算最后2 ns研究体系的结构是否为相对稳定的状态，将16~18 ns体系的势能进行分析。结果如图3所示，无论体系温度高低，16~18 ns模拟期间体系势能均十分稳定。对于不同的模拟时刻，温度越高体系具有的势能越高，且势能的起伏越大，符合基本物理规律。

Figure 3. Potential energy stabilizations for CH₄ hydartes

图3. 甲烷水合物势能的稳定性

3.2. 稳定结构的孤立面

为了更直观地反映不同温度对甲烷水合物形成的影响，将最后2 ns结构中水分子位移小于0.2 nm的认定为稳定的水分子，并提取该类水分子形成的孤立面，结果如图4所示。由图4可以明显看出，随着温度的升高，孤立面包络的体积明显变小，这表明在1 GPa高压条件下，温度决定了生成天然气水合物的数量，在温度较高时，甲烷可以与水分子形成混合物，但在1 GPa高压的条件下，必须同时有较低的温度条件才可以形成较稳定的甲烷水合物结构。

(a) 250 K (b) 260 K (c) 270 K (d) 280 K (e) 290 K (f) 310 K

Figure 4. Spatial distribution of the water molecules that are not “moving”

图4. 稳定水分子形成的孤立面

3.3. 径向分布函数

为了表征不同温度下天然气水合物形成过程中原子的位置信息而引入径向分布函数(对关联函数)进行表征，结果如图5所示。从图5(a)可以看出甲烷与甲烷分子之间的最近距离为0.5 nm，与天然气水合物晶体中甲烷之间的距离相吻合，且距离中心分子较远处的甲烷分子分布未趋于1，说明甲烷分子的位置已经相对固定并形成类似骨架结构；氢原子与氧原子的径向分布函数(图5(b))可以看出在0.1 nm附近有一个十分尖锐的峰，该峰代表水分子中的氢–氧键距离，距离中心氢原子1 nm以外的氧原子分布较均匀，形成明显的液态或无定型物质的特征径向分布函数，证明部分水分子在甲烷分子之间的空隙中以游离态存在，氧–氧径向分布函数可以给出相同的结论。

(a) C-C (b) H-O (c) O-O

Figure 5. Radial distribution function of different atoms

图5. 不同原子之间的径向分布函数

3.4. 均方根位移

为更清楚了解不同温度条件下形成的天然气水合物结构各原子的运动情况，分别研究甲烷和水分子的均方根位移(MSD)如图6。甲烷分子的均方根位移与水的均方根位移的数值差距为一个数量级，甲烷分子的均方根位移值均小于0.1 nm，可以认为无论温度高低研究体系中所有甲烷分子几乎都在其平衡位置附近振动，甲烷分子的均方根位移的斜率均趋近于零，由爱因斯坦关系可知，整个研究时间段内甲烷分子不会发生扩散；而水分子的均方根位移受温度影响变化较大且具有一定斜率。水在整个研究体系中无论体系的温度高低总有一部分水分子是可以随意移动的且扩散几率很高，随温度升高扩散系数明显增大，这也证明了图4中温度越高稳定结构孤立面包络体积减小相一致。

(a) 甲烷分子 (b) 水分子

Figure 6. MSD of (a) the methane molecules; and (b) the water molecules in the gas hydrate crystallization simulations

图6. 晶化后的气体水合物的(a) 甲烷分子的均方根位移；(b) 水分子的均方根位移

3.5. 环结构

众所周知，目前在自然界中发现的天然气水合物主要有三种晶体结构，分别为I型(5¹²6²)、II型(5¹²6⁴)和H型(5¹²6⁸) [27] [28]。I型和II型天然气水合物包含的五元环结构主要是图7(a)和图7(b)中所示的结构，未发现有图7(c)结构的五元环存在；I型的六元环结构主要是图7(e)和图7(f)所示的结构，未发现图7(d)所示的结构；而II型六元环结构主要为图7(d)和图7(e)结构。而高压下天然气水合物的结构中同时出现了图7所示的所有结构，证明在高压下形成天然气水合物结构为I型和II型的混合结构，与文献 [29] 压力大于1 Gpa出现II型结构的结论一致。

(a) 含0个完整H₂O (b) 含1个完整H₂O (c) 含2个完整H₂O (d) 含0个完整H₂O (e) 含1个完整H₂O (f) 含2个完整H₂O

Figure 7. Five-membered ring and six-membered ring structures in natural gas hydrate

图7. 天然气水合物中的五元环及六元环结构

3.6. 环数量(Rc)

为比较不同温度条件下研究体系生成天然气水合物中所包含环的数量，引入了归一化环数量(Rc) [30]，图8(a)可以清楚看出随温度升高五元环的Rc值不断下降，由此可知随温度上升，体系的五元环数量下降，且310 K时已经下降近50%；而由图8(b)可知，六元环的R_C值同样表现为随温度上升而下降；通过比较五元环Rc值与六元环Rc值可以获得两种不同的环的数量比约为5:1，介于I型的6:1和II型的3:1之间，可以推测，所形成的天然气水合物为I型和II型混合物，与环结构所得结论一致。

(a) 含10个原子的环 (b) 含12个原子的环

Figure 8. The rings per cell (Rc(n)) of nature gas hydrate

图8. 天然气水合物每单元所含环数量

4. 结论

采用分子动力学方法对天然气水合物在1 Gpa高压下的形成过程进行了模拟，利用均方根位移、对关联函数及环统计的方法对不同温度下形成天然气水合物的结构进行表征。结果表明：

1) 1 Gpa高压条件下所生成的天然气水合物由于受模拟时间限制，甲烷分子和部分水分子处于稳定振动态而部分水分子随温度升高扩散能力增加。

2) 对于稳定水分子和可移动水分子进行局部环境分析可知，稳定水分子周围形成稳定的多环结构，而可移动水分周围形成的环结构数量少、结构差。

3) 更详实的环统计结果表明，体系的五元环与六元环数量比值为5:1，介于I型与II型天然气水合物水合比之间，同时也发现体系存在II型天然气水合物的特征五元环。

4) 通过Rc可知体系的五元环数量多且相互联通，而六元环数量少联通性较差。

基金项目

海南省自然科学基金(519MS025)；海南省教改项目(Hnjg2019-25)。

NOTES

^*第一作者。

^#通讯作者。