2021, Vol. 37
当今国际能源资源短缺，不仅要寻求节能环保措施，更要寻找可代替清洁能源。天然气水合物是全球第二大碳储库，广泛存在于冻土地带和海洋地层中，其中所蕴藏的天然气资源十分可观［1］。但若是对其盲目开发易破坏其地质结构，导致地层坍塌，甚至引发不可控制的洋灾害［2］，因此需对其进行深入研究。水合物作为蓄冷介质，其相变温度适宜，蓄冷密度大，储−释冷的传热效率高，流动性较好［3］，在蓄冷空调、油气开采、气体分离和消防灭火等领域被广泛利用。目前对于水合物研究主要是通过实验方法，但实验研究过程总存在种种限制。由于实验装置的限制，无法研究特殊条件下水合物的相关特性，严重影响其研究进程。此外，实验研究只能观测到水合物的宏观性质，无法从其微观结构进行分析。但经过多年的发展，分子动力学模拟（molecular dynamics，MD）具有强大的模拟计算体系的分子动力学模拟已经较为成熟，可以更为直观地提供很多实验过程中无法得到的微观细节。
利用各种模拟方法对模拟数据进行处理分析，为研究提供了大量有利结论，使其在许多科学和实用领域中得到广泛应用。目前对于水合物方面的分子动力学研究进展也有了长足的进步，近十年来基于模拟的水合物生成特性研究大幅增加。关于笼形水合物MD模拟计算作为一种新兴的研究方法，也逐渐成为其发展主流。
1 基本模拟方法 1.1 水合物生成微观原理水合物是由主体水分子借氢键结合形成笼形结晶，客体分子被包围在其晶格之中形成的稳定笼形晶体包络物。目前已发现的气体水合物的晶胞构型，按结构可分为Ⅰ型（SI）、Ⅱ（SII）型、H型（SH）三种［4］。水合物需要在高压低温的环境下生成，其生成过程也经历了复杂的多元变化，主要包括成核和生长两个阶段。水合物生成微观原理如图1所示。
|
图 1 水合物生成微观原理 Fig.1 Microscopic principle of hydrate formation |
首先，主体分子之间依靠氢键形成空的分子笼；之后，客体分子被输送到合适位置，并在合适的温度和压力条件下进入主体分子笼中，形成亚稳态分子簇；最后，该不稳定团簇不断发展达到临界尺寸，从而实现水合物的成核。当晶核满足临界状态条件后，水合物便会在其诱导下迅速生长。
1.2 水合物基本模拟方法MD模拟早在1957年就由Alder等［5］实现，研究了32～500个刚性小球分子系统的微观运动，并于1959年提出将其推广到更为复杂的分子体系，以模拟研究其结构和性质［6］。目前，MD模拟技术已经成熟，水合物分子动力学模拟是目前最广泛采用的计算庞大、复杂体系的计算方法，即通过合理的计算软件，选用合适的力场、势函数以及系综等搭建科学的计算模拟体系进行模拟计算。再通过可视化程序对其结果进行观测并分析，是目前水合物MD模拟方法的主流。具体水合物生成特性的分子动力学模拟研究流程如图2所示：
|
图 2 水合物生成特性的分子动力学模拟研究流程 Fig.2 Molecular dynamics simulation research process of hydrate formation |
目前研究使用的分子动力学模拟软件主要有Lammps、Gromacs、DEPOLY、Materials Studio、NAMD等。通过以上软件进行相关数据分析处理：
（1）氢键的形成。水合物成核初期首先是其原子间氢键的形成，其标准为：氧原子间距离小于0.35 nm，H−O···0角度小于30°。
（2）径向分布函数（radial distribution function，RDF）。RDF表示为指定粒子位置，可通过研究水合物生长过程中不同原子之间RDF第一峰的高度变化，探究分子的运动转移，判断水合物是否形成。
（3）均方位移（mean square displacement，MSD）。MSD体现的是体系内粒子空间位置偏离最初位置的程度。在水合物生长过程中，MSD会随着模拟的进程逐渐上升，在水合物生成完成后于稍高于0位置上下浮动，并保持稳定。
（4）Z轴原子密度。水合物形成后客体分子困于笼中，由于晶体结构完整，Z轴原子密度则会呈现周期分布：高低峰交替出现。
（5）扩散系数。在水合物的形成过程中，水分子之间形成的强大氢键网络将客体分子紧紧包围并挤压，其扩散系数远小于正常客体分子。
（6）能量分析。水合物生成是氢键形成，释放能量，整体势能减小，水合物形成后势能再次稳定达到平衡。
通过以上分析方法进行结构性质、扩散系数以及能量等分析，证明水合物中水分子之间所形成的氢键网络以及后续稳定笼形物的形成，从而证实水合物的成核以及生长机理。但目前对于水合物的分子动力学模拟研究中所涉及的大小和时间尺度仅为nm和ns级别，相对于水合物成核的实验研究来说，目前是不切实际的，仍需未来进行更多对比研究。
2 水合物成核分子动力学模拟 2.1 水合物成核机理水合物在成核过程需要提供较高驱动力［7］，是水合物生成的重要步骤。且由于水合物成核的随机性，其成核现象的时间尺度跨度较大，所形成笼形不断变化［8］，增加了其分子动力学模拟研究的难度。部分学者进行了天然气水合物的成核模拟［9-11］，发现其生长机理是可重复的：成核过程是一个从初始状态的无序、无定型到最后形成相对稳定的晶体状态的过程。水合物成核过程如图3所示。
|
图 3 水合物成核过程 Fig.3 Hydrate nucleation |
此外，Hawtin等［9］研究结果表明，其不稳定团簇达到临界尺寸（为至少包含三个笼型，即至少包含20个客体分子）时，才能诱导水合物的快速生长。
2.2 水合物成核的影响因素近年来的研究结果表明，水合物成核除了需要高驱动力外，还包括很多影响因素，包括客体分子溶解度、不同界面性质、第三组分介质的加入等。
2.2.1 溶解度对水合物成核的影响首先，客体分子的溶解度不同，其水合物的成核速率也会不同。Vatamanu等［10］研究了甲烷水合物的成核机理，发现局部富含甲烷的区域将更易成核，并且快速结晶，不需要晶体的“记忆效应”。Nada等［12］研究了稀溶液中甲烷水合物的成核，结果发现，稀溶液中生长出的笼胞会有20%的空笼，同时也证明了甲烷的浓度决定了水合物的生长速率。此外，Vatamanu等［13］研究结果表明：客体分子的初始溶解度越高，晶胞生长越迅速。但客体分子浓度过高时会产生超临界效应，气泡的产生会使系统趋于相分离，从而会在一定程度上抑制水合物的成核。
2.2.2 界面性质对水合物成核的影响首先，不同的系统界面形状会对成核产生影响。Walsh等［14］通过研究球形、薄片形以及圆柱形界面形状发现：界面的曲率越高，其气体溶解度也会更高，成核率也就更高。
此外，系统界面的大小也会对成核率产生影响。其主要原因是临界区域会对初始状态造成不稳定的影响，即大系统中的成核率更高。Sarupria等［15］则研究了无界面系统CH4水合物的成核，观察到在几百ns的时间尺度上，水合物便会形成亚临界水和甲烷分子簇，之后聚集产生临界水合物核，其成核速率得到明显增加。
Defever等［16］认为表面与水的亲疏性也会影响水合物的成核。其结果表明，水合物成核的表面效应主要来自于它们对系统中客体分子的影响，而客体分子更易在疏水表面聚集。
2.2.3 第三组分的加入对水合物成核的影响第三组分的加入改变了成核环境的结构，同时也会对水合物成核产生影响。Bai等［17］研究了二氧化硅表面对水合物成核的影响。结果发现，羟基化二氧化硅表面CO2水合物逐层出现，成核速率增高。六边形冰表面也可以促进CH4水合物的异质成核［18］。研究发现，伴随着甲烷分子的积累，由于冰表面结构的不相容，对有缺陷结构形成了诱导，使成核过程偏向于在某一特定方向产生，从而促进新的结晶相的异质成核。
除此之外，研究人员也对多孔介质中水合物的形成进行了研究。Liang等［19］进行了氢氧化硅介质中CH4水合物成核的分子动力学模拟。研究发现，尽管两个固相的结构不匹配，但由于局部结构的灵活性，使硅石表面连接性更好，气体水合物与硅石表面之间形成了半笼状结构。最终，甲烷水合物在羟化硅形成的孔隙内成核，从而从MD角度有效证明了多孔介质有助于促进笼形水合物的生成。
综合目前水合物分子动力学模拟研究表明，具有可重复性的水合物成核过程被大量研究，探究其确切成核机制也成为MD模拟的一大热点。此外，使水合物达到临界尺寸成核的影响因素很多，其界面性质是其成核过程的重要影响因素。通过不断优化界面结构，加入合适异质结构以促进水合物成核将会是未来研究的重要方向。
3 水合物生长分子动力学模拟水合物成核之后便进入生长阶段，也是水合物生成研究的重要模拟主题之一。而其生长的传热传质过程十分复杂，研究［7］表明其异质成核和生长可能是水合物形成的唯一机制。目前笼形水合物的生成模拟主要研究了其生长的微观机理，以及不同热力学条件和抑制剂加入等对其生长速率的影响。
3.1 水合物生长机理近年来部分学者分别进行了甲烷水合物 ［12, 20］、乙烷水合物［8］、二氧化碳水合物［21］等生长机理的模拟研究。本课题组进行了二氧化碳水合物生成的分子动力学模拟，在260 K、30 MPa下其不断生长过程如图4所示。
|
图 4 水合物成核后的不断生长过程（260 K，30 MPa） Fig.4 Hydrate growth after nucleation（260 K，30 MPa） |
Tung等［21］进行了I型甲烷水合物生成的分子动力学模拟，发现伴随系统压力增高，甲烷溶解度增大，其生长速率也随之增加。但伴随着甲烷的消耗，甲烷浓度的改变也会对甲烷水合物的动态生长产生重要影响。为此，Liang等［22］研究开发了一种新方法促进甲烷向增长界面的转移，这不仅限制了甲烷的大规模的无效转移，而且降低了甲烷到达超饱和状态的可能性，从而促进了甲烷水合物的生长。
Vatamanu等［23］研究了甲烷水合物的异质晶体生长。发现生长的晶体中含有许多缺陷，通常水合物中约20%的水合物笼未被占用，并且这不依赖于所施加的温度和压力。
鉴于水合物开采技术过程中出现的种种挑战，人们越来越关注水合物中二氧化碳和甲烷的分子交换过程。Yi等［24］探究了CH4−CO2混合水合物生长的微观机理。图5为CH4−CO2混合水合物的成长过程，图中：方形标注的为CO2分子；圆形标注的为CH4分子。伴随时间转移，先是CH4分子占据CO2分子原本所占小笼[图5（c）]，后是被排出的CO2分子占据CH4分子所占大笼[图5（f）]。通过该研究可发现：混合水合物生长过程中CO2优先进入水合物。这为甲烷水合物的置换技术提供了更为可行性的技术证明。
|
图 5 CH4−CO2混合水合物的生长过程 Fig.5 Growth process of CH4-CO2 mixed hydrate |
目前研究表明，水相中水溶性第三组分的存在也会影响水合物生长动力学。通过加入第三组分的抑制剂以解决水合物运输线路中造成的阻塞问题，是当前水合物生长特性研究的主要方向。
3.2.1 热力学抑制剂对水合物生长的影响目前研究认为，在水相中无机盐和有机分子（醇）可被用作热力学抑制剂以有效地防止水合物形成。Choudhary等［25］研究了在中等温度压力条件下，低浓度抑制剂（质量分数1%）氯化钠和甲醇水溶液中甲烷水合物形成动力学的影响。结果发现，在试验初期，与纯水系统相比，质量分数为1%NaCl和1%甲醇系统能够有效抑制水合物的生长。但随着实验进程不断发展，氯化钠系统中水合物形成的动力学超过纯水系统，而甲醇系统的水合物形成动力学仍低于纯水系统。
此外，氨基酸类热力学抑制剂因其天然可降解的绿色环保特性，逐步引起各方关注［26-27］，也有相关学者对其进行了分子动力学模拟研究［28-29］。其中，Maddah等［28］进行了添有不同浓度（质量分数小于1.5%）的甘氨酸、丙氨酸、脯氨酸和丝氨酸的水合物分子动力学生成模拟，结果发现，氨基酸抑制剂均起到了良好抑制作用，且其效果为：丝氨酸 > 甘氨酸 > 丙氨酸 = 脯氨酸。
3.2.2 动力学抑制剂对水合物生长的影响动力学抑制剂（KHIs）的存在会增加体系无序度，降低水的结晶度，从而抑制水合物的生成。研究表明，PVP作为主流动力学抑制剂对水合物生成产生良好的抑制作用［30-31］。Li等［30］研究结果如表1所示。结果表明，这三种抑制剂均起到了抑制作用，其中抑制剂PVP−A效果最优。这是因为在生成过程中，KHI与水分子之间相互结合，破坏水和水合物之间的氢键，此外，双键氧的强吸引力促使更多水分子与KHI的结合，从而阻碍水合物的生长。
|
表 1 不同抑制剂作用下系统平衡时间对比 Table 1 Comparison of system equilibrium time under different KHIs |
研究表明，二元外来分子的存在不仅可以促进成核，也有助于提高水合物结构的稳定性［32］。Liang等［33］发现羟基化二氧化硅表面存在促进了甲烷水合物晶体的生长，并且溶液中的晶体生长总在二氧化硅表面附近进行。这可能是由于二氧化硅表面上的硅烷醇基团与水分子相比具有较低的极性，因此二氧化硅表面为甲烷气泡形成提供了合理的位置。目前的模拟结果表明，矿物表面的存在可能在天然气水合物的晶体生长中起关键作用，为这种晶体生长过程的微观机制提供了新的见解。
由于水合物生长过程的复杂特性，目前对于水合物生长分子动力学的模拟研究相对较少，主要集中在通过加入抑制剂以增加其生成的平衡时间，从而抑制水合物的生长。目前对于抑制剂的MD研究主要集中在普通的热力学抑制剂、动力学抑制剂等，更多天然型的环保抑制剂，如天然生物质和糖类等仍待关注。并且更多生长过程的微观细节特性也仍待发现。
4 结 论本文基于水合物分子动力学模拟的研究现状，总结了水合物研究的基本模拟方法，并从微观角度探究了水合物MD模拟中成核和生长的具体生成特性。根据近年来大量的水合物分子动力学研究成果表明，其模拟结果无论是与实验结果相比较，还是对其进行物理化学分析，均可证明其结果的合理性，大大推进了水合物方向的研究进程。
目前对于水合物分子动力学模拟虽然有了大量研究，但距离实际应用仍有一定差距。笔者认为在以后的研究过程中，还有以下工作有待开展：
（1）探寻更为合适准确的模拟体系，配合以更为精准的力场、系综以及势能参数，通过强大的计算机技术进行高效水合物MD模拟。
（2）优化模拟生成环境，不仅考虑到不同热力学条件，界面性质，还可在模拟过程中加入合适的金属离子、纳米粒子、天然蛋白质、糖类物质、防聚剂等，通过模拟实验研究出合适的环保型材料以促进水合物生成特性研究；
（3）目前水合物分子动力学模拟方向较为单一，因此模拟过程中不仅要考虑到水合物结构与形成水合物的客体分子的多样性，更要了解水合物尚未明确的基本性质，如稳定性、导热特性、流动特性等。
（4）MD模拟研究为微观层面的探讨，而实验研究为宏观现象的分析。将模拟研究与实验研究相结合，通过模拟实验的初步筛选以减少实验台研究的重复度，通过对比实验研究以提升研究结果的准确性，将两者相结合以推动未来高效的水合物生成特性研究。
| [1] |
MAKOGON Y F, HOLDITCH S A, MAKOGON T Y. Natural Gas-hydrates—a potential energy source for the 21st Century[J]. Journal of Petroleum Science and Engineering, 2007, 56(1−3): 14-31. DOI:10.1016/j.petrol.2005.10.009 |
| [2] |
KVENVOLDEN K A. Gas hydrates-geological perspective and global change[J]. Reviews of Geophysics, 1993, 31(2): 173-187. DOI:10.1029/93RG00268 |
| [3] |
LIU N, ZHOU J L, GAO M, et al. An experimental study on flow and heat transfer characteristics of binary hydrate slurries in a horizontal tube[J]. International Communications in Heat and Mass Transfer, 2019, 102: 34-41. DOI:10.1016/j.icheatmasstransfer.2019.01.009 |
| [4] |
LICENCE P. Clathrate hydrates of natural gases[J]. Fuel, 2008, 87(13): 3158. |
| [5] |
ALDER B J, WAINWRIGHT T E. Phase transition for a hard sphere system[J]. The Journal of Chemical Physics, 1957, 27(5): 1208-1209. DOI:10.1063/1.1743957 |
| [6] |
ALDER B J, WAINWRIGHT T E. Studies in molecular dynamics. I. General method[J]. The Journal of Chemical Physics, 1959, 31(2): 459-466. DOI:10.1063/1.1730376 |
| [7] |
BARNES B C, SUM A K. Advances in molecular simulations of clathrate hydrates[J]. Current Opinion in Chemical Engineering, 2013, 2(2): 184-190. DOI:10.1016/j.coche.2012.12.002 |
| [8] |
WILSON D T, BARNES B C, WU D T, et al. Molecular dynamics simulations of the formation of ethane clathrate hydrates[J]. Fluid Phase Equilibria, 2016, 413: 229-234. DOI:10.1016/j.fluid.2015.12.001 |
| [9] |
HAWTIN R W, QUIGLEY D, RODGER P M. Gas hydrate nucleation and cage formation at a water/methane interface[J]. Physical Chemistry Chemical Physics, 2008, 10(32): 4853-4564. DOI:10.1039/b807455k |
| [10] |
VATAMANU J, KUSALIK P G. Observation of two-step nucleation in methanehydrates[J]. Physical Chemistry Chemical Physics, 2010, 12(45): 15065-15072. DOI:10.1039/c0cp00551g |
| [11] |
JACOBSON L C, MATSUMOTO M, MOLINERO V. Order parameters for the multistep crystallization of clathrate hydrates[J]. The Journal of Chemical Physics, 2011, 135(7): 74501. DOI:10.1063/1.3613667 |
| [12] |
NADA H. Growth mechanism of a gas clathrate hydrate from a dilute aqueous gas solution: a molecular dynamics simulation of a three-phase system[J]. The Journal of Physical Chemistry B, 2006, 110(33): 16526-16534. DOI:10.1021/jp062182a |
| [13] |
VATAMANU J, KUSALIK P G. Molecular Insights into the heterogeneous crystal growth of sI methane hydrate[J]. The Journal of Physical Chemistry B, 2006, 110(32): 15896-15904. DOI:10.1021/jp061684l |
| [14] |
WALSH M R, BECKHAM G T, KOH C A, et al. Methane hydrate nucleation rates from molecular dynamics simulations: effects of aqueous methane concentration, interfacial curvature, and system size[J]. The Journal of Physical Chemistry C, 2011, 115(43): 21241-21248. DOI:10.1021/jp206483q |
| [15] |
SARUPRIA S, DEBENEDETTI P G. Homogeneous nucleation of methane hydrate in microsecond molecular dynamics simulations[J]. The Journal of Physical Chemistry Letters, 2012, 3(20): 2942-2947. DOI:10.1021/jz3012113 |
| [16] |
DEFEVER R S, SARUPRIA S. Surface chemistry effects on heterogeneous clathrate hydrate nucleation: A molecular dynamics study[J]. The Journal of Chemical Thermodynamics, 2018, 117: 205-213. DOI:10.1016/j.jct.2017.08.021 |
| [17] |
BAI D S, CHEN G J, ZHANG X R, et al. Microsecond molecular dynamics simulations of the kinetic pathways of gas hydrate formation from solid surfaces[J]. Langmuir the ACS Journal of Surfaces and Colloids, 2011, 27(10): 5961-5967. DOI:10.1021/la105088b |
| [18] |
PIRZADEH P, KUSALIK P G. Molecular Insights into clathrate hydrate nucleation at an ice-solution interface[J]. Journal of the American Chemical Society, 2013, 135(19): 7278-7287. DOI:10.1021/ja400521e |
| [19] |
LIANG S, KUSALIK P G. The nucleation of gas hydrates near silica surfaces[J]. Canadian Journal of Chemistry, 2015, 93(8): 791-798. DOI:10.1139/cjc-2014-0443 |
| [20] |
TUNG Y T, CHEN L J, CHEN Y P, et al. The growth of structure i methane hydrate from molecular dynamics simulations[J]. The Journal of Physical Chemistry B, 2010, 114(33): 10804-10813. DOI:10.1021/jp102874s |
| [21] |
TUNG Y T, CHEN L J, CHEN Y P, et al. Growth of structure i carbon dioxide hydrate from molecular dynamics simulations[J]. The Journal of Physical Chemistry C, 2011, 115(15): 7504-7515. DOI:10.1021/jp112205x |
| [22] |
LIANG S, KUSALIK P G. Explorations of gas hydrate crystal growth by molecular simulations[J]. Chemical Physics Letters, 2010, 494(4-6): 123-133. DOI:10.1016/j.cplett.2010.05.088 |
| [23] |
VATAMANU J, KUSALIK P G. Heterogeneous crystal growth of methane hydrate on its sII [001] crystallographic face[J]. The Journal of Physical Chemistry B, 2008, 112(8): 2399-2404. DOI:10.1021/jp077583k |
| [24] |
YI L Z, LIANG D Q, ZHOU X B, et al. Molecular dynamics simulations for the growth of CH4-CO2 mixed hydrate
[J]. Journal of Energy Chemistry, 2014, 23(6): 747-754. DOI:10.1016/S2095-4956(14)60208-4 |
| [25] |
CHOUDHARY N, KUSHWAHA O S, BHATTACHARJEE G, et al.. Molecular dynamics simulation and experimental study on the growth of methane hydrate in presence of methanol and sodium chloride[J]. Energy Procedia, 2017, 105: 5026-5033. DOI:10.1016/j.egypro.2017.03.1008 |
| [26] |
SA J H, LEE B R, PARK D H, et al. Amino acids as natural inhibitors for hydrate formation in CO2 sequestration
[J]. Environmental Science & Technology, 2011, 45(13): 5885-5891. |
| [27] |
SA J H, KWAK G H, HAN K, et al. Gas hydrate inhibition by perturbation of liquid water structure[J]. Scientific Reports, 2015, 5: 11526. DOI:10.1038/srep11526 |
| [28] |
MADDAH M, MADDAH M, PEYVANDI K. Molecular dynamics simulation of methane hydrate formation in presence and absence of amino acid inhibitors[J]. Journal of Molecular Liquids, 2018, 269: 721-732. DOI:10.1016/j.molliq.2018.08.108 |
| [29] |
洪春芳. CO2水合物分解及导热特性的分子动力学模拟[D]. 上海: 上海理工大学, 2018. .
|
| [30] |
YI L Z, JIANG F, QIN H B, et al. Molecular dynamics method to simulate the process of hydrate growth in the presence/absence of KHIS[J]. Chemical Engineering Science, 2017, 164: 307-312. DOI:10.1016/j.ces.2017.02.029 |
| [31] |
孙慧翠, 王韧, 王建华, 等. 不同分子量聚乙烯吡咯烷酮对四氢呋喃水合物形成及生长的影响[J]. 天然气工业, 2018, 38(5): 125-132. DOI:10.3787/j.issn.1000-0976.2018.05.015 |
| [32] |
ERFAN-NIYA H, MODARRESS H. Molecular dynamics simulation of structure H clathrate-hydrates of binary guest molecules[J]. Journal of Natural Gas Chemistry, 2011, 20(6): 577-584. DOI:10.1016/S1003-9953(10)60242-3 |
| [33] |
LIANG S, ROZMANOV D, KUSALIK P G. Crystal growth simulations of methane hydrates in the presence of silica surfaces[J]. Physical Chemistry Chemical Physics, 2011, 13(44): 19856-19864. DOI:10.1039/c1cp21810g |