关于I 型水合物晶体晶面对水合物生长动力学影响的实验设计与实践

赵 莉，方 彬，刘 冰，郑海霞，展凯云

（1. 中国石油大学（华东） 理学院，山东 青岛 266580；2. 中国地质大学（武汉） 工程学院，湖北 武汉 430074）

开设综合性研究型实验课程是一种新颖的教学模式，有助于激发学生的求知欲，培养创新能力、创新思维，扩宽知识面，提升综合运用各种交叉学科知识的能力[1-3]。分子动力学模拟是近年来飞速发展的一种微观尺度仿真方法，是研究凝聚态物理的有力工具，其在对水合物微观乃至宏观性质的研究中取得了重大突破，如对热力学和力学性质、形成与分解等的研究方面。这是一门结合物理、数学和化学的综合技术，它以经典力学、量子力学和统计力学为基础，利用数值求解经典力学运动方程得到体系的运动轨迹，并给出相应的结构信息与性质。本文所述的实验依托纵向科研项目，结合分子动力学的科研热点和专业发展前沿，以水合物生长动力学为研究载体，探索了水合物晶体不同晶面对水合物生长动力学的影响，生动展示了固体力学的基础理论，拓展了学生的科研视野。

1 实验内容设计依据

水合物为客体小分子与水分子在低温高压条件下形成的笼型化合物，客体分子位于由水分子构成的氢键网络（笼状结构）中，最常见的客体分子为甲烷，即形成甲烷水合物[4]。自然界中永久冻土和深海沉积物中蕴藏大量甲烷水合物资源，其中海洋沉积物中占据90%，保守估计为传统化石能源（煤、石油、天然气）的2倍[5-6]。由此可见，水合物资源具有巨大的能源潜力，是未来可替代传统能源的一种清洁能源。2017年，我国在南海海域水合物试采中取得了重大突破。目前，水合物已被我国列为第173个矿种，它的产业化是我国未来能源战略的重要组成部分。除此之外，水合物特殊的结构特性使其在二氧化碳封存[7]、气体混合物分离[8-9]、储氢[10-11]等领域具有广泛的应用。

现阶段，由于开采技术的限制以及复杂的沉积环境等不利因素，人们对水合物的研究依旧停留在实验室以及试开采阶段。作为潜在的可替代能源，天然气水合物资源的调查勘探以及未来的开发利用，都离不开对其形成条件和规律的研究。然而，水合物资源的开采研究是一项复杂的过程，涉及物理化学、工程地质学、环境学、机械制造等多种学科。本文主要针对的是水合物的生长研究，该过程对了解水合物在地层中的赋存模式具有至关重要的作用，是能源工业中天然气成藏研究的重要课题。Larsen等人[12]通过单晶生长实验观察到，生长速度慢的晶面发育较好，而生长速度快的晶面通常在生长过程中会消失，不同晶面生长的差异性直接影响到水合物的分布模式。但是，到目前为止该过程的微观机理尚不明确。

2 实验仪器

计算机、Material Studio软件（MS）、Gromacs开源软件包、VMD可视化软件，以及Origin绘图软件。

3 实验教学设计

3.1 实验预习

（1）实验开始之前，教师将固体力学相关书籍和学习资料、Gromacs和MS操作手册、VMD可视化软件以及Origin绘图软件等的使用说明发送给学生，引导学生自主学习。使学生充分了解固体力学的相关知识理论，掌握晶体类型以及晶面定义，掌握分子动力学模拟的理论基础和Gromacs输入文件的构建，能够对可视化软件及数据处理软件进行灵活运用。

（2）教师提前准备关于水合物生长研究的相关文献发给学生，让学生提前阅读，并在此基础上，引导学生自主查阅资料，理清在文献阅读中所遇到的问题，找出文献中的研究空白。

（3）学生通过阅读以上资料，以及自主查阅文献了解整个模拟流程和理论基础，并尝试进行实验设计。

3.2 实验操作

3.2.1 模型构建

模型构建是计算机模拟的关键过程，在本实验设计中，准确构建水合物单晶胞是研究的基础与关键。本文选取的是I型甲烷水合物，其模型搭建过程分两步：第一步从X射线晶体学研究结果中获取水合物晶体主体水分子中氧原子的初始位置；第二步在已经确定好位置的氧原子上插入相应的氢原子，并调整氢原子方向使其满足 Bernal-Fowler法则，同时保证单晶胞偶极矩和势能最小[13]。该晶胞包含8个甲烷分子与46个水分子（比例为1∶5.75）。在此单晶胞的基础上，通过 MS可视化软件切出我们所关注的三个晶体面（100）、（110）、（111），如图1所示，红色表示水分子，可通过氢键作用组成笼型结构，客体分子甲烷（灰色）则占据晶笼内部空间。相应的晶胞参数如表1所示。另外，对应单晶胞的三个晶面分别构建超晶胞：（100）晶面构建 3×3×2 超晶胞；（110）晶面构建 3×2×3超晶胞；（111）晶面构建1×2×4超晶胞。随后在Z方向将超晶胞体系延伸10 nm，留下10 nm空白空间，并在该区域随机放入甲烷与水分子，甲烷与水分子比例为1∶5.75，满足I型甲烷水合物比例（1∶5.75），三个晶面体系初始构型的参数与甲烷水分子数如表2所示。

表1 不同晶面所对应的晶体参数以及切除时的表面向量

表2 三个晶面体系初始构型的参数与甲烷水分子数

3.2.2 模型优化及计算模拟

在利用分子动力学模拟之前，需要对初始构型进行优化，使其能量最小，保证其初始构型稳定，随后进行100 ns的分子模拟用于数据采集与分析。采用带随机项的速度重新标度实现温度耦合，并利用Parrinello-Rahman扩展-集成压力控制对模拟盒子向量进行重新标度[14-15]。恒温以及恒压常数分别设置为0.2 ps和2.0 ps，并利用Particle-Mesh Ewald (PME)求和方法处理长程静电相互作用[16]。原子间的范德瓦尔斯相互作用用 Lennard-Jones势来描述，截断半径为12 Å，该半径应小于仿真体系向量的一半。采用跳蛙算法用于牛顿运动方程的积分，时间步长设置为2 fs[17]。在X、Y、Z三个方向均考虑周期性边界条件，水分子用 Tip4p/2005模型描述[18]，甲烷分子采用OPLS-UA势函数[19]，模拟采用NPT系综和周期性边界条件。温度压力分别维持在260 K和50 MPa，因为在该温压条件下水合物生长驱动力较大，有利于生长并缩短模拟时间。

为了便于结果讨论，我们将初始构型分为两个部分，如图2所示，左边为水合物晶核，右边为水气混合液相。同时，为了更好地表征水合物的生长过程，右边水气混合相被均匀分割为10等份，并分别标记为L1—L10。

图2 水合物生长动力学初始模型（以100表面体系为例）

为了确定在生长过程中哪些水分子和甲烷分子属于水合物相，必须要有判断标准。一般情况下，水氧原子可以在冰或水合物相中与相邻水的4个氧原子形成四面体单元的顶点，F3参数基于104.25° O—O—O角排列，这与液相水氧原子中是非常不同的[20]。因此，在水合物生长模拟过程中，可以用Baez和Clancy提出的F3阶参数表征水分子的局域状态[21]。该算法可以提供固态水分子氢键网络与标准四面体结构的偏离。

式中，θjik是3个水氧原子之间的角度，第i个水氧原子为中心，j和k水氧原子在i原子两边且位于i水氧原子第一水壳层以内（O—O原子径向分布函数的第一极小值，半径3.5 Å）。无序结构的F3值大于四面体结构的F3值，如笼形水合物。因此在我们的模拟中定义，当F3小于0.05时，水分子属于水合物相。此外，当晶笼中多个宿主水分子脱离分子间氢键约束时，甲烷分子很容易从不完全水合物氢键笼中逃逸。因此，假设甲烷分子（5.5 Å）在第一水合物壳层内的水合物相（F3＜0.05，定义如上）周围的水分子数目大于15，则认为甲烷分子仍处于水合物的多面体晶笼中。

4 结果与讨论

4.1 水合物生长动力学（以110体系为例）

水合物生长过程是一个相变过程，水气混合液相逐渐变为水气包裹且具有固定对称性的晶体结构，水合物生长过程是以水合物晶面为基底，以层方式逐步向液相生长的（见图3）。水合物/液相界面先由水分子形成半晶笼，然后吸附客体甲烷分子用以完成生长。可以看出，在开始阶段，由于甲烷分子随机分布于水气液相中，水合物生长较快，随后由于甲烷分子在水中溶解度较低，而此时甲烷分子又过饱和，于是甲烷分子快速聚集成纳米气泡，受周期性边界条件的影响，液相中甲烷气呈现圆柱体状，此时，水合物晶体生长驱动力较大，部分甲烷分子克服界面张力扩散至水合物晶体表面新形成的空腔中，该晶笼完成客体分子生长，但此过程需要时间较长，生长过程也较为缓慢。

4.2 不同晶面生长比较

图4 模拟时间结束水合物生长体系快照（Y-Z方向）

从图4可以看出，不同晶面水合物体系生长快慢不同。同为100 ns生长时间，100表面体系液相并未完全生长，而110及111表面体系生长均已完成。此外，100表面体系在模拟结束时，由于生成水合物量较少，大量甲烷气体汇聚并单独成为一相，整个体系有水合物相和液相，呈气相分离状态。110及 111两个晶面体系生长完成，所有液相均已转化为固相，但同时也有残余的少量气体分子聚集，说明生长过程中并不是所有的晶笼都被客体甲烷分子占据，存在部分未占据的水合物晶体。相比 110体系，111体系残余气体数更少。为更好地描述水合物生长的快慢，我们统计不同液相层L1—L5水合物相水气分子数，如图5所示。

图5所示水合物生长是从L1—L5依次生长的。相比之下，100表面体系只完成到L2层的生长，110、111表面体系已经完成了L1—L5层的生长。另外，110与111的晶面体系虽然都完成了生长，但其完成时间却有很大差异。110表面体系大概在 90 ns左右完成L5层的生长，而111表面体系在60 ns左右就完成L5层水合物生长。相比之下，111表面生长速率最大，110表面次之，而100表面最慢，这与实验观察所得结果基本一致，生长速率越慢，实际存在越多。

图5 L1—L5层中水合物相水分子以及甲烷分子数随时间的变化

为进一步阐明不同水合物晶面生长快慢的内在机制，统计了三种晶面在Z轴方向上的密度变化，如图6所示。而表3中统计了不同晶面上水分子以及甲烷分子的数量。可以看出，不同晶面的密度以及分子数是不同的。就密度来说，111面密度最大，110面次之，而 100面最小；从分子数上说，111面所包含的甲烷数最多。因此，晶体生长快慢与表面密度以及甲烷分子数有很大关系，表面甲烷数越多，生长越快。晶面结构与主客分子数共同决定了生长速率。

图6 不同晶面体系在Z轴方向上的密度变化

表3 不同晶面水合物晶体单位体积粒子数

4.3 实验内容拓展

本实验属于理论探究型实验，以实际问题为背景进行理论研究，并构建了水合物晶面生长的动力学模型。本实验采用的是分子尺度常用的仿真方法——分子动力学模拟。实验中，通过问题分析构建不同分子尺度模型，结合统计力学及涨落耗散等理论综合分析仿真结果，从纳米级尺度上模拟水合物晶体的生长过程，并解释了不同水合物晶面生长快慢的内在机制。学生在教师的指导下完成整个模拟实验过程，不仅了解了科研选题，进一步掌握了物理化学、化学工程、能源技术和环境技术等理论知识及研究的基本方法，同时还可激励学生的求知欲望和探索精神。此外，考虑课时量及学生情况的不同，还对实验内容进行了深层次的拓展，使模拟内容更加接近实际应用。扩展内容如下：

（1）由于 I型水合物晶体对称性为Pm3n，除文中提到的几种常见晶面外，还有多种晶面，学生可自行学习并使用软件构建初始构型，进行动力学模拟。

（2）由于真实环境下水合物地层温度压力会有所不同，学生可以研究不同温压条件下水合物的生长规律及不同晶面在该生长过程中的作用。

（3）在真实开采或流动保障过程中，会在流体中酌情添加促进剂或抑制剂，一般为有机或无机小分子（特殊情况下也有聚合物等），学生可自行查阅文献，自主设计实验，研究水合物在有添加剂情况下的生长情况及对晶面的影响。

（4）除传质过程外，水合物生长过程为放热过程，在考虑动力学过程的同时，学生也可通过更换系综来研究传热快慢对生长的影响。此外，水合物在相变过程中液固变化导致的导热系数改变，也可对水合物后续生长产生影响。

4.4 撰写实验报告

实验课后，要求学生对实验原理进行概括总结，对实验数据进行整理分析，并按照科研论文的格式撰写实验报告。此外，学生还需阐述对该课题的理解，并提出实验过程需改进的地方。对于优秀学生开展的拓展实验内容部分，由教师协助对报告进行修改，并以科技论文形式向有关专业期刊投稿。

5 结语

本综合性创新实验通过分子动力学模拟，研究了I型水合物晶体不同晶面对生长动力学的影响。研究发现，111面最快，110面次之，100面的生长最慢。该结论与实验结果一致，且对研究甲烷水合物的勘探开发利用具有重要意义。

该实验包含的知识体系复杂，涉及对物理化学、经典力学、分子化学、凝聚态物理、数学等多个交叉学科的综合运用，有助于培养学生的科研创新能力、问题解决能力、团结协作能力，以及知识体系综合运用能力。到目前为止，魏琦等6位本科生在专业期刊发表了相关的科研论文；刘晓旭、杨光瑞等4位学生通过自主查阅文献拟定了相关的大学生创新科研项目并获批。实践表明，通过设计性实验的综合训练，可有效锻炼学生的科研探索精神，并大大缩短其融入科研团队的时间。