_型甲烷水合物稳定性的分子动力学模拟

Jones 参数。
( ) u（rij） =
1 r13
Fra Baidu bibliotek
+
1 r14
+
1 r23
+
1 r24
+ 4q2 4πε0 r78
( ) － 2q2 4πε0
1 r18
+
1 r28
+
1 r37
+
1 r47
s + a1exp（ － b1r56）
（6）
+ a2［exp（ － b2r13） + exp（ － b2r14） + exp（ － b2r23） + exp（ － b2r24）］
动方程采用 Verlet leap － frog 算法，时间步长采用 0. 8fs，模拟时
间长度为 300ps，其中前 150ps 用于使系统达到平衡。
σij
=
1 2
( σi
+ σj )
（3）
参数 ε/K σ / A
εij = 槡εi εj
（4）
表 2 L － J 势能参数［15］
水分子 － 甲烷分子 346. 7 3. 2
关键词：甲烷水合物；分子动力学模拟；稳定性；晶穴占有率
The Impact of the Stability of DtructureⅠMethane Hydrate by Molecular Dynamics Simulation*
QI Ying － xia，TANG Cheng － wei，LIU Dao － ping，YANG Guang （College of Energy and Power Engineering，University of Shanghai for Science and Technology，Shanghai 200093，China）
法，截断半径为 0. 6nm，长程静电力则采用 Ewald 加和法计算，使 用洛仑茨 － 贝赛那（ Lorentz － Berthelot） 混合方法［14］确定相互交
叉作用参数，如公式 （3） 、（4） 所示，表 2 中给出了各分子间的
L － J势能参数。温度由单独的 Nose － Hoover 恒温水浴控制，运
9. 4，10. 9， 11. 9，14. 5
Gallaway，Ruska ［7］
6. 00 ± 0. 16
De Roo，Peters ［8］
6. 3，7. 4
Handa ［9］
6. 00 ± 0. 01
2 分子动力学模拟方法
2. 1 模拟系统和模拟方法的建立
本文采用经典的等温等容（ NVT） 系综进行 MD 模拟［10 ， － 12］
角度进行的光谱直接测量方法，它主要是通过 NMR 或者拉曼光
谱观测晶体的结 构 来 获 得 水 合 物 的 晶 穴 占 有 率 ，进 而 计 算 出 水
合数。众所周知，这 种 光 谱 直 接 测 量 法 目 前 还 存 在 着 一 定 的 技
术难题，因此获得的有效的实验数据较少，特别是甲烷水合物的
有效数据就更少。在下面的表 1 中我们列出了一些专家学者通
研究者
宏观直接测量
微观间接测量
Ripmeester，Ratcliffe ［5］
6. 05 ± 0. 06
Uchida，Hirano ［6］
6. 08，6. 13，6. 15， 6. 16，6. 2，6. 21，6. 22，
6. 28，6. 28，6. 3
6. 5，6. 9，7. 1， 7. 2，7. 6，8. 2，
相互作用进行描述，这主要是因为 MCY 势能模型在描述冰结构 的计算中证明是非常成功的［17］，而且冰的分子结构呈六方晶体
结构，与水合物的笼型结构极其相似。在这个模型中，水分子被
看作是刚性分子，其中 O － H 键长为 0. 09572 nm，H － O － H 键角
为 104. 52°，MCY 势能模型将分子间的势能表达为长程静电相互
作用和短程原子间的相互作用的贡献，如公式（6） 所示，表 3 中
给出了 MCY 的势能模型参数。
[ ( ) ( ) ] uL －J = 4ε
σ rij
12
－
σ6 rij
（5）
式中，rij 和 u（ rij ） 分别为分子 i 和 j 之间的距离和相互作用势
能，q 为水分子的电荷数，ε0 为真空电容率，ε 和 σ 为 Lennard －
·64·
广州化工
2010 年 38 卷第 2 期
nH
=
46 8
5.
75
（1）
而晶穴占有率 f 是指水合物中被客体分子占有的晶穴与所
有的晶穴之比，因此水合物的水合数也可以写成下面的关系式：
nH
= nwater nguest
= nwater ncoges × f
（2）
式中：ncages ———表示水合物晶体内的所有的晶穴数
Key words： methane hydrate； molecular dynamics simulation； stability； cage occupancies
目前，天然气水合 物 的 研 究 在 世 界 范 围 内 受 到 高 度 重 视 的 原因首先在于它被公认为是 21 世纪的重要后续能源，近二十年 来在海洋和冻土 带 发 现 的 天 然 气 水 合 物 资 源 特 别 巨 大 ，甲 烷 的 总资源约为 1. 81016 ～ 2. 1 × 1016 m3 ，有机碳储量相当于全球已探 明矿物燃料（ 煤、石油、天然气） 的两倍［1］。天然气也是一种洁净 能源，有利于环 保 和 国 民 经 济 的 可 持 续 发 展。天 然 气 水 合 物 是 自然界中天然气 存 在 的 一 种 特 殊 形 式，它 分 布 范 围 广、规 模 大、 能量密度高，研究、开发和利用水合物中的天然气资源已成为各 国政府在能源领域的当务之急。所谓天然气水合物是指以甲烷 为主的烃类气体与水在高压低温的条件下生成的非化学计量、 外观似冰的笼 型 水 合 物［2］，俗 称 可 燃 冰。近 些 年 来 人 们 对 气 体 水合物展开了广 泛 而 深 入 的 研 究，但 研 究 的 重 点 主 要 集 中 在 水 合物的实验和动 力 学 宏 观 方 面，但 这 些 实 验 研 究 还 不 足 以 解 释 水合物生成、分解的微观机理，而分子动力学模拟可以为研究水 合物的微观结构 提 供 强 有 力 的 工 具，而 计 算 化 学 在 模 拟 气 体 水 合物的生成和分解的过程研究方面具有不可比拟的优势［3 － 4］。
1. 2 实验测得的水合数
通过实验测量水 合 物 中 水 合 数 的 方 法 大 致 可 分 为 两 种 ，一
种是从宏观 角 度 进 行 的 间 接 热 力 学 测 量 方 法，它 主 要 是 利 用
PVT 参数，通过水合物在分解或生成过程中释放或消耗的气体
量和水量来确定水合物中的水合数。另一种测量方法是从微观
之间的势能利用 Lennard － Jones（ LJ） 势能函数来计算，如公式
（5） 所示。由于水分子为极性分子，考虑到水分子间的相互作用
对水合物的晶体结构具有至关重要的作用，所以我们采用 MCY （ Matsuoka，Clementi and Yoshimine） 势能模型［16］对水分子之间的
过实验测得的Ⅰ型甲烷水合物的水合数的结果。从表 1 中我们
可以看出，通过不同的实验获得的水合数会存在一定的偏差，利
用微观直接测量得到的水合数比利用宏观间接测量得到的水合
数的值要小一些。这可能是由于间接测量的水合物中含有一些
自由水，会导致最后计算得到的水合数有些偏高。
表 1 实验测得的甲烷水合物在 0 ℃ 时的水合数 nH
甲烷分子 － 甲烷分子 148. 6 3. 758
2. 2 势能模型
我们所研究的Ⅰ型甲烷水合物晶体的模拟体系中只包含两
种分子：水分子 和 客 体 甲 烷 分 子。为 了 减 少 计 算 机 工 作 量 和 运
行时间，我们把非极性甲烷分子看作是球形分子，把水分子看作
是带电荷的非线 性 的 极 性 分 子，甲 烷 分 子 之 间 和 甲 烷 与 水 分 子
Abstract： The possible hydration numbers of CH4 clathrate hydrates were investigated by computer simulations. Molecular dynamics simulations were performed for the various cage occupancies at the constant － NVT condition. The simulations results showed that even the partial occupation of the cages the clathrate host lattice remained stable. The system potential energy was negative and decreased as a function of cage occupancies. The pressure was positive and decreased roughly as a function of cage occupancies，except for the case when the host lattice was melting. The CH4 molecules exert on stabilization of the cages through the vander Waals intermolecular interactions. The simulation results were compared with experimental hydration numbers.
模拟体系只有单个晶胞组成，包含 46 个水分子和 0 ～ 8 个甲烷分
子。模拟体系盒子的原始尺寸为 1. 2 nm × 1. 2 nm × 1. 2nm ，采
用周期性边界条件（ 对于任意一个晶胞，其周围都有一个一模一 样的晶胞） ，利用 X 射线单晶衍射实验［13］得到分子的初始位置。
本次模拟针对各分子间的 Lennard － Jones 作用项采用硬球截断
1. 1 水合数
水合物中天然气 的 含 量 与 水 合 数 有 关，所 谓 水 合 数 就 是 指 每个晶胞中的水分子数与客体分子数之比。比如水分子数 nwater = 46，客体分子数 nguest = 8 时，甲烷水合物中的水合数为：
* 基金项目：上海市浦江人才计划（08PJ1408300） ，上海市重点学科建设项目（ S30503） 资助。 作者简介：祁影霞（1964 － ） ，女，博士，副研究员，主要从事能源与动力，气体水合物和分子动力学研究。E － mail：peggyqi@ yahoo. cn
1 水合物的晶体结构
天然气水合物中，水分子（ 主体分子） 形成一种空间点阵结 构，气体分子（ 客体分子） 则填充在水分子点阵间的晶穴中，气体 和水之间没有固定的化学计量关系。形成点阵的水分子之间靠 较强的氢键结合，而 气 体 分 子 和 水 分 子 之 间 的 作 用 力 为 范 德 华 力。到目前为止，已经发现的天然气水合物结构主要有三种，即 结构Ⅰ型，结构Ⅱ型和结构 H 型。其中结构 I 型水合物单晶是 体心立方结构，包含 46 个水分子，由 2 个小晶穴和 6 个大晶穴组 成。小晶穴为五边形十二面体（512 ） ，由 20 个水分子组成，形状 近似为球形；大晶穴是由 12 个五边形和 2 个六边形组成的十四 面体（51262 ） ，由 24 个水分子组成，为扁球形结构。
2010 年 38 卷第 2 期
广州化工
·63·
Ⅰ型甲烷水合物稳定性的分子动力学模拟*
祁影霞，汤成伟，刘道平，杨 光
（ 上海理工大学能源与动力工程学院，上海 200093）
摘 要：利用计算机模拟研究了 CH4 笼型水合物中的晶穴占有率对甲烷水合物稳定性的影响。模拟结果表明晶穴在部分被占
有的情况下，水合物的笼型结构依然能够保持稳定，水合物的势能和压力与晶穴占有率有关，势能随着晶穴占有率的增大而减小，且 表现为负值，压力则表现为正值，并且随着占有率的增大而减小，晶体结构被分解的情况除外，通过模拟结果与实验结果的比较，最后 推断出最有可能稳定存在的是晶穴被全部占有时的水合物。
+ a3［exp（ － b3r16） + exp（ － b3r26） + exp（ － b3r35） + exp（ － b3r45）］