不同结构气体水合物稳定性的分子动力学研究

尹晓霞,焦丽君,2,李元超,华泽珍,李裕斌

(1.青岛职业技术学院海尔学院，山东 青岛 266555;2.中国石油大学新能源学院，山东 青岛 266580)

在能源短缺、环境污染的大背景下，天然气水合物作为未来战略能源的勘探和开采，受到越来越多的关注，可以利用二氧化碳水合物比甲烷水合物稳定的特点[1]，通过置换法开采实现天然气水合物的开采，同时完成二氧化碳的封存。因此，二氧化碳水合物的稳定性和分解特性的研究尤为重要。此外，由于二氧化碳水合物相平衡压力相对较高，在海水淡化、空调蓄冷、烟气捕集等领域的应用受到限制，二氧化碳水合物和环戊烷等热力学促进剂相结合，节能效果显著[2-4]。

气体水合物是一种由小分子气体或者易挥发的液体与水相互作用，在低温高压条件下形成的类冰状固体。主体水分子之间通过氢键形成笼形网络，客体分子位于晶穴中心，与水分子通过微弱的范德华力连接，形成稳定的水合物[5]。常见的客体分子有：甲烷、乙烷、丙烷、二氧化碳、氢气、氮气、环戊烷、四氢呋喃、四丁基溴化铵等。根据客体分子的大小，水合物通常分为三种结构，分别是sI型、sII型和sH型水合物。晶体特性如表1所示[6]。客体分子不同，气体水合物的成核特性、分解特性及稳定性差别较大。Zhang等[7]通过对七种不同客体分子在水溶液中成核过程发现，相同浓度下，客体分子尺寸越大，自扩散系数越低，越容易生成水合物。Kondori等[8]发现，加入热力学抑制剂后，II型气体水合物分解特性和稳定性均表现出明显差异，抑制程度排序为甲醇>乙醇>甘油。张庆东[9]通过对比不同客体分子对气体水合物稳定性的影响，发现对于含有甲烷的二元水合物体系，二氧化碳、乙烷、硫化氢和丙烷对水合物稳定性的促进作用逐渐增加。表明二元水合物可以增强甲烷水合物的稳定性。王浩瑄等[10]发现，客体分子所处晶穴结构不同，对甲烷水合物的稳定性影响较大，大晶穴内部的客体分子作用更加显著。目前关于气体水合物的稳定性研究大多集中在微观模拟，而且以甲烷水合物为主，二氧化碳水合物研究较少。

表1 气体水合物晶体特性

二氧化碳气体分子直径为0.33 nm，通常形成Ⅰ型水合物，可以同时占据大小晶穴[11]。环戊烷、四氢呋喃等分子直径较大，通常占据Ⅱ型水合物的大晶穴[12]。当二氧化碳和环戊烷(CP)客体分子同时存在时，可以形成二元水合物，Lee和Kim通过X射线衍射实验发现，环戊烷完全占据大晶穴，二氧化碳占据一半以上的小晶穴。他们指出，这种二元的水合物结构比单一二氧化碳水合物更加稳定[13]。目前的研究大多集中在宏观实验层面，缺乏微观视角，本文从微观角度出发，从分子水平上研究不同结构的水合物在相同热力学条件下的稳定性。

1 模拟方法

为了研究不同结构气体水合物的稳定性，本文建立了Ⅰ型和Ⅱ型二氧化碳水合物、Ⅱ型环戊烷水合物、以及二氧化碳-环戊烷二元水合物分子模型，如图1所示。水合物分子模型通过Materials Studio软件建立，根据X射线衍射结果确定水分子内氧原子和氢原子的空间坐标，然后将客体分子的中心放置在水分子构建的笼形结构中心[14]。为了保证晶穴数目相当，对于I型水合物，构建3X3X4超晶胞，其中包括1 565个水分子和288个晶穴；对于II型水合物，构建2X2X3超晶胞，其中包括1 632个水分子和288个晶穴。

图1 水合物的初始构型 (图中线表示氢键，球体表示氧原子和碳原子)

分子动力学模拟基于牛顿力学为理论基础，通过分子势能函数产生分子间相互作用，计算原子受力及加速度，在非常短的时间间隔内，得到各分子的速度及位置。因此，分子势能函数的选择对模拟结果的影响较大。本文中水分子模型选用三点模型SPC/E，该模型简单且高效，能够很好的预测液态水的自扩散系数，而且在水合物相平衡温度、分解焓和活化能的预测等方面也能与实验值较好的吻合[15-16]。二氧化碳分子模型采用EPM2模型，该模型能够准确预测二氧化碳气体在水中的溶解度以及气液共存线[17-18]。环戊烷分子模型直接采用Materials Studio软件库中的CVFF力场模型。各分子模型的势能参数如表2所示。

表2 各分子模型参数

分子间相互作用力包括Lennard-Jones势能和库仑力，见公式(1)

(1)

不同原子间的Lennard-Jones势能作用遵从Lorentz-Berthelot作用规则，见公式(2)(3)

(2)

(3)

模拟使用LAMMPS软件，时间步长设为1 fs，采用速度-Verlet 算法，应用周期性边界条件消除界面影响。范德华力的作用截止半径设为12Å，采用PPPM(Particle-Particle Particle-Mesh)方法来计算长程静电力，计算精度为10-6。采用Nose-Hoover控温和控压，温度和压力阻尼系数分别设为0.1 ps和1 ps，将水合物体系在控制温度和压力的条件下(NPT系综)，模拟200 ps。

2 结果与讨论

2.1 径向分布函数

径向分布函数(Radial Distribution Function(RDF)是指以一个参考原子α为中心，与其距离为r到r+Δr之间，β原子的数目与系统中平均数目的比值。从物理意义上说，可以理解为系统中区域密度与平均密度的比值，可以反映分子的微观结构。见公式(4)

(4)

式中g(α-β)(r)——原子α和原子β之间的径向分布函数;

V——模拟盒子的体积;

Nα和Nβ——两种原子的数目;

r——两个原子之间的距离;

niβ(r)——在距离原子α为r的条件下，β原子的数目。

图2为模拟初期t=20 ps时不同水合物体系所对应的O-O原子及C-C原子的径向分布函数。对于水中氧原子的径向分布函数，所有体系第一个峰值出现在r=2.79附近，第一个峰值代表水合物晶穴中两个相邻氧原子之间的距离，因此所有体系的值均相同。这与文献结果相差不大，证明了模型的正确性和模型方法的有效性[19-20]。第二和第三个峰值分别代表水合物中相隔一个和两个水分子的距离，对于Ⅰ型水合物和Ⅱ型水合物有偏差，体现了水合物不同的网格结构，Ⅱ型水合物氧原子平均距离相对远些。对于Ⅰ型二氧化碳分子中碳原子的径向分布函数，在稳定状态下两个碳原子的最小距离是6.81，其值与文献相一致[19]。对于Ⅱ型水合物，二氧化碳分子中碳原子在稳定状态下的距离最小值为6.21，这种差异主要受到水合物结构的影响，Ⅱ型水合物中碳原子距离更近。而且，通过RDF的峰值越高，表明局部密度相对于平均密度越大，原子越密集。因此，从b图可以看出，二元CO2+CP水合物的稳定性要强于Ⅱ型CO2水合物，这与环戊烷是常用的热力学促进剂的说法相一致，环戊烷通过占据Ⅱ型水合物的大晶穴，可以有效的改善相平衡条件，使水合物变得更加稳定[21]。对于环戊烷中的碳原子分布，第一个峰值在r=3.93附近，指的是同一个分子中碳原子之间的距离，与水合物结构关系不大，因此研究意义不大。

图2 270 K,5 MPa条件下t=20 ps时不同体系氧原子和碳原子径向分布函数

2.2 均方位移

分子动力学模拟是基于牛顿力学原理的，分子受力产生加速度，结果是模拟过程的每一时刻分子的受力和位置都不相同。通常用在某时刻与初始时刻相比，分子位移平方的平均值来表示分子运动的程度，称为均方位移(Mean Square Displacement MSD)。见公式(5)

(5)

图3是在270 K 5 MPa热力学条件下，不同体系水分子和客体分子的均方位移，从图中可以明显看出，Ⅱ型CO2水合物最不稳定，在模拟初期(t=30 ps)MSD即开始迅速上升，体系内分子运动剧烈，表明由固态的水合物体系逐渐变为液态，说明水合物发生分解。Ⅰ型水合物在模拟后期(t=150 ps)MSD开始上升，表明体系经历了一段亚稳定状态，最终开始分解。而对于Ⅱ型的CP水合物和CO2+CP二元水合物，在模拟期间，MSD值几乎不变，尤其是水分子的MSD为0，表明水合物结构十分稳定，一直没有分解。水分子的MSD为0，可以看出在模拟过程中，水分子之间通过氢键结构构建水合物晶格，水分子不能旋转和移动。而客体分子的MSD大概在2左右，表明在水合物中客体分子可以自由的旋转或振动，在晶穴内部发生小位移的偏移，但不能逃离水合物晶穴。

图3 270 K,5 MPa条件下不同体系水分子和客体分子均方位移随时间变化曲线

2.3 F3序参数

为了进一步分析水合物的分解情况和稳定性，定义了序参数[22]。见公式(6)

F3,i=〈(cosθjik|cosθjik|+cos2(109.472))2〉jk

(6)

式中θjik——以氧原子i为中心，与距离其3.5Å内两个氧原子j和k，三者形成的夹角。

借助此式来反映氢键网格偏离理想四面体的情况，从而监测水合物的稳定性。水合物状态下，F3=0，液态水状态下，F3=0.1。不同体系下F3序参数的随时间变化趋势如图4所示。从图中可以清晰的看出，结构Ⅱ型CO2水合物最不稳定，模拟初期即开始分解，分解时间大概50 ps，结构I型CO2水合物次之，在模拟中后期开始出现分解趋势，分解时间大概100 ps，分解相对较慢。因此，结构Ⅰ型水合物相对结构Ⅱ型水合物更加稳定，这也可以解释分子动力学模拟水合物合成过程中，甲烷和水作用主要生成Ⅰ型水合物，生成Ⅱ型水合物较少[23]。而结构Ⅱ型的CP水合物和二元水合物序参数均保持不变，水合物没有发生分解，这与均方位移分析结果一致。

图4 270 K,5 MPa条件下不同体系F3序参数随时间变化曲线

3 结 论

本文使用分子动力学模拟方法，通过对径向分布函数进行分析，充分证实了该模拟方法的正确性，同时从微观层面进一步认识了水合物的基本构型，包括结构Ⅰ型、结构Ⅱ型水合物，以及不同客体分子的微观构型。此外，碳原子密度分布、系统势能、构像图等多种参数均反映了上述性质，由于篇幅限制未完全列出。主要结论及展望概括如下：

(1)通过均方位移和F3序参数，对比分析了不同结构水合物的稳定性，发现对于CO2客体分子，由于尺寸较小，在Ⅰ型水合物构型下相对比较稳定。

(2)相对于CO2水合物，大分子环戊烷水合物更加稳定，而且将环戊烷和二氧化碳分别置于Ⅱ型水合物的大小晶穴后，发现对于二氧化碳水合物而言，加入热力学促进剂的二元水合物的热力学稳定性较单一二氧化碳客体分子水合物要好。

(3)本文从分子视角探索水合物分解特性和稳定性的作用机理，为水合物的宏观实验研究和工程应用积累了理论基础，同时在该工作的基础上，将来可以继续研究不同构型和不同客体分子水合物的界面分解特性。