太赫兹电磁波调控甲烷水合物分解*

朱智 闫韶健 段铜川 赵妍 孙庭钰 李阳梅

1) (上海理工大学光电信息与计算机工程学院，教育部医用光学仪器与设备重点实验室，上海 200093)

2) (国防科技创新研究院，太赫兹生物物理创新工作站，北京 100071)

甲烷水合物具有分布范围广、资源储量大、能量密度高、清洁无污染等特点，因此一经发现便得到人们的广泛关注.虽然如此，甲烷水合物的开采存在较多的困难，主要涉及笼型甲烷水合物的分解过程，因此本文通过分子动力学模拟，探索太赫兹电磁波对该水合物分解的特异性影响.通过分析甲烷水合物中水分子的振动谱，发现其区别于体相水在低温下的特异性振动峰.通过引入频率与该振动峰中心频率一致的太赫兹波发现，该电磁场会破坏甲烷水合物原有的氢键网络，降低甲烷周围水分子的配位数，从而使甲烷分子从水笼中脱离.进一步地，对序参数F4 的分析也证明，在太赫兹电磁刺激下，甲烷壳层水由原来的晶体水变成了液态水.同时证明了该频率具有相比于其他频率的绝对优势，因此具有其特异性.本文的发现有望通过非热效应来实现甲烷水合物的分解，高效、低能耗地提高开采效率，推动未来新能源的发展.

1 引言

甲烷(CH4)水合物(或称可燃冰)，是甲烷被水分子包围所形成的冰晶笼状物，具有污染小、贮藏量大、分布地域广、能量密度高等优点，被作为石油、天然气的新时代替代能源而备受期待，与它相关的科学研究受到了全世界的高度重视.它的开采、运输和贮存均存在很大的困难，重点涉及笼型水合物的形成与分解，同时也关乎着全球的气候变化，一但开采不当，将会产生巨大的负面影响[1]，因此，甲烷水合物的成核、生长和分解的分子动力学行为成为了研究的重点[2].

分子动力学模拟是一种研究分子行为的有效方法，科研人员基于该方法针对甲烷水合物的分解做了许多工作.Alavi 和Ripmeester[3]采用非平衡分子动力学模拟研究甲烷水合物的分解过程，结果表明，甲烷水合物的分解速率主要受笼状结构破坏、甲烷分子在水合物相与溶液相的交界面的逸出以及在水溶液中扩散的传热和传质的影响.颜克凤等[4]提出热激法分解甲烷水合物，并用分子动力学模拟方法系统地研究了高温液态水与甲烷水合物接触的分解机理.Yan 等[5]用分子动力学模拟基于“真空去除法”和常规方法分别研究了甲烷水合物的分解和减压机理，发现减压能够促进水合物的分解.Ding 等[6]通过分子动力学模拟总结得出甲烷水合物的分解过程，即首先水合物笼型结构扭曲变形直至破坏，接着甲烷从晶穴中逸出，最后无客体分子填充的空腔逐步瓦解崩塌.他们进一步研究表明，分解过程会形成一层液层，该液层阻碍了水分子以及甲烷分子的扩散，从而降低了水合物的分解速率[7].Yagasaki 等[8]研究了水合物抑制剂(甲醇和氯化钠)对甲烷水合物在水溶液中解离的动力学影响，结果表明，甲醇可以提高初始阶段的解离速率，而氯化钠降低了解离速率，这种差异源于两种热力学抑制剂对甲烷水合自由能的相反影响.Myshakin 等[9]利用分子动力学方法分别模拟了甲烷在甲烷水合物中的占有率对甲烷水合物分解过程的影响，发现甲烷的占有率越小，甲烷水合物的分解速率就越快.Bai 等[10]研究了CO2置换甲烷水合物的过程，得出客体分子的化学能、记忆效应及质量转移会影响置换效果，置换过程分两步进行，首先是甲烷水合物表面的分解，然后CO2在分解的表面上形成新的化合物.在Smirnov[11]的研究中，通过将甲烷水合物与沸石接触，改变氢键网络，进而破坏水合物的笼型结构，最终促进水合物分解.

以上的研究对甲烷水合物分解的基本物理过程进行了描述，并分析了外加的诸多因素如压强、化学物质等对分解速率的影响.由于笼状甲烷水合物中存在极性相互作用，近些年来研究人员开始关注电场对水合物分解的影响.Luis 等[12]研究了外部静电场对水合物的影响，发现当电场强度提高到1.5 V/nm 时，原有的氢键排序被打断，水分子变得有序，进而促进了水合物的分解.进一步地，Xu 等[13]采用分子动力学模拟方法在温度T=260 K 和压强P=100 bar (1 bar=105Pa)下，分别研究了静态与2.45 GHz—1.0 THz 振荡电场作用下甲烷水合物的分解过程.结果表明，静电场对甲烷水合物的影响存在强度阈值，超过阈值1.5 V/nm 的静电场才可以诱导水分子沿电场排列，从而推动水合物的分解.而对于振荡电场，随着电场频率的提高，其强度阈值也需要提高.虽然如此，以往的研究中均未考虑电场频率的特异性，由于有机分子的本征振荡频率以及分子间相互作用能所对应的频率均在太赫兹波段，因此特定的太赫兹场极有可能通过谐振作用改变分子的构象以及影响分子之间的弱相互作用，对此学者们进行了广泛的研究并已经在太赫兹场的生物效应方面取得了较大的进步[14-22].

本文通过分析甲烷水合物中水成分的振动吸收谱发现，甲烷水合物中水之间的氢键网络在10.3 THz 存在1 个不同于低温下体相水的吸收峰.进一步地，当对甲烷水合物引入10.3 THz 电磁刺激后，该电磁波会与氢键产生共振，极大地破坏水合物之间的氢键，水合物的水分子配位数随之降低.更进一步地，F4有序参数也能较好地反映出，随着氢键网络被破坏，水合物晶体结构发生较大的改变，甲烷分子从笼中脱离，水由晶态变为液态水.同时，证明了该频率相比其他频点对促进甲烷水合物的分解具备非常大的优势，因此可以看出，太赫兹电磁波对水合物分解具有特异的调控作用(如图1(a)).

2 模拟方法

首先运用Packmol[23，24]软件进行初始模型的构建，沿着模拟盒子的xyz三个方向得到8 × 4 × 4的半笼型甲烷水合物，其次通过Gromacs[25]软件对水分子添加虚原子，形成四位点水.将该甲烷水合物在NVT 系综下进行高温融化，温度为500 K，融化时间为5 ns，随后将融化后的构型与原构型放在一起，最终形成三相混合溶液体系，包含1024 个甲烷分子和5888 个水分子，如图1(b)上图所示.该体系所有的C 原子都用Lennard-Jones参数表征，范德瓦耳斯半径为1.0 nm，甲烷力场使用OPLS-AA 力场，水模型特别选用TIP4P/Ice[26]刚性水模型，因为该模型能够完美再现甲烷水合物的三相平衡温度.将初始构型能量最小化，先在NPT 系综下进行10 ns 的平衡，然后在NVT 系综下进行10 ns 的平衡，最后将得到的构型作为后面模拟的起始构型.在xyz三个方向应用周期性边界条件，为了保持恒定的温度，使用Nose-Hoover[27，28]恒温器来控制温度，截止半径为1.0 nm，积分步长为2 fs，以0.2 ps 的间隔时间记录原子轨迹，后续所有的模拟均使用软件Gromacs 2020 在恒定的NVT 条件下进行.

为了测试起始结构的稳定性，首先在温度260 K 下进行了200 ns 的NVT 系综模拟，结果如图1(b)中图所示.可以看出，甲烷水合物可以稳定存在，而且经过高温融化后聚集的甲烷气体在260 K 温度下逐渐散开.这表明，在温度为260 K且没有施加外电场的情况下，甲烷水合物逐渐稳定为成核的状态.

图1 (a) 太赫兹加速分解甲烷水合物的概念图，其中蓝色基底为甲烷水合物晶体，红色和白色小球为甲烷分子中的碳和氢原子，笼状物为包裹甲烷的水分子，其中黄色小球为水中的氧原子.(b) 上图为模拟的初始构型，其中绿线左侧为笼型甲烷水合物，右侧为高温融解后的甲烷和水，聚集的蓝色部分为甲烷气体，周围红色为水分子;中图为260 K 温度下NVT 平衡200 ns 后甲烷水合物的状态，其中交界部分已经出现成核的现象，同时聚集的甲烷分子散开，有继续成核的趋势;下图为260 K 温度下施加特定频率的太赫兹电磁刺激后，甲烷水合物的状态，大部分原有的甲烷水合物已经分解，并且有进一步分解的趋势Fig.1.(a) Conceptual graph describing terahertz wave accelerated decomposition of methane hydrate，where the blue substrate is methane hydrate crystals，the red and white balls are the carbon and hydrogen atoms，and the clathrate is the water molecules enveloping the methane.(b) Above:Initial simulated configuration.The left side of the green line is caged methane hydrate，while the right side is methane and water mixture after high temperature melting.The blue cluster therein is the methane gas，surrounded by water molecules (in red).Middle:State of methane hydrate after the NVT equilibration for 200 ns at a temperature of 260 K.The nucleation has occurred in the interface，and the initially gathered methane molecules have partly diffused and are expected to form more nucleation.Bottom:State of methane hydrate after a specific terahertz electromagnetic (THz-EM) stimulation at 260 K.Most of the original methane hydrate has been decomposed and developed into a methane cluster.

3 结果与讨论

为了研究特定频率THz 电磁波对甲烷水合物的影响，首先计算了甲烷水合物中水的振动吸收谱.如图2 所示，与体相水中水的吸收峰不同的是，甲烷水合物在频率为10.3 THz 时存在1 个明显的吸收峰，可见该特征频率的电磁波可以被甲烷水合物中的冰水吸收而几乎不能被体相水吸收.基于这一原理，对体系分别引入频率为10.3，5 和40 THz的电磁刺激，后两种频率作为非特征频率进行对比.我们知道，电磁波存在电场分量和磁场分量，由于磁分量对于热运动的水分子影响微乎其微，因此电场分量主导着电磁波与水分子的相互作用[14].所以在分子动力学模拟中，仅引入电场分量，强度均为2 V/nm，极化方向为x轴正方向，在温度260 K，NVT 系综下对体系进行研究.作为对比，在无电场引入时，分别在260 K 正常温度和330 K高温下进行两组模拟.

图2 水在260 K 温度下的振动吸收谱，其中红线为笼状结构甲烷水合物中水的振动吸收谱，黑线为体相水的振动吸收谱，可见甲烷水合物中存在10.3 THz 的吸收峰，而体相水对该频率下的太赫兹电磁刺激只有弱吸收，使用该频率的电磁刺激能够特异性影响甲烷水合物的氢键网络.内插图为甲烷水合物的笼状结构，中心蓝白色球棍结构代表甲烷分子，其外圈包围的为水分子Fig.2.Vibrational absorption spectra of water at a temperature of 260 K.The red line corresponds to the spectrum of water in the caged methane hydrate，while the black one denotes the spectrum of bulk water.There exists an absorption peak at 10.3 THz in methane hydrate but an absorption valley in bulk water.Hence，an EM stimulus at this specific frequency could alter the hydrogenbond network of methane hydrate.Inset:caged methane hydrate.The inner blue-white ball-stick structure denotes methane molecule，surrounded by water molecules in the outer.

首先考察水合物网络中氢键的变化.氢键是一些与电负性较大的原子(施体)成键的氢原子和附近电负性较大或者带孤对电子的原子(受体)形成的较强的非键作用，其作用强度介于成键作用和非成键作用之间.氢键是否形成主要是由氢原子-施体-受体所成角度θ与施体-受体之间距离d共同决定的，一般认为θ＜ 30°，d＜ 0.35 nm 时构成氢键.C 原子吸引电子的能力较弱，一般不形成氢键，所以本研究体系中的氢键主要来源于水分子之间形成的氢键.当水合物处于稳定状态时氢键的数目相对稳定，且数目远多于体相水中的氢键数目.体系氢键数目随时间变化的趋势如图3 所示，可以看出，在不加电磁波的情况下，水合物的氢键数目相对稳定;T=330 K 高温下，受水分子热运动影响，氢键断裂，水合物分解.而在温度260 K 下，当加入电磁波刺激后，频率为10.3 THz 的电磁波会与氢键产生共振，进而被吸收，导致水合物之间的氢键断裂.但是频率为5 和40 THz 的非特征频率电磁波由于不会完全被吸收，所以对于水合物的氢键网络影响较小.为了更加明显地看出水合物的状态，通过甲烷周围水的配位数来确定，由图4(a)可知，C—O 径向分布的第一峰谷位置出现在r=0.57 nm，因此以0.57 nm 为半径来确定每个甲烷分子第一壳层内水的配位数.Luis 等[12]提到，水合物为五边形十二面体时水的配位数为24，也就是五边形十二面体的顶点数，当配位数小于12[29]时，认为水合物已经完全分解.从图3 可知，水合物的配位数与氢键数变化相对应，高温和10.3 THz特征频率下的电磁波刺激均可以使水合物分解，而非特征频率的电磁波仅仅破坏了部分氢键，甲烷分子周围水分子的平均配位数仍然在20 附近，这表明甲烷分子仍然处于笼子中，并未完全分解.

图3 (a) 甲烷水合物在结晶态(右上插图)以及在太赫兹场刺激下分解后(右下插图)甲烷分子和水的空间分布形态;(b) 模拟体系总的氢键数目随时间的变化;(c) 引入电场强度为2 V/nm 情况下，体系氢键损失率与频率的关系;(d) 水合物中甲烷周围水分子的配位数，这个配位数是在0.57 nm 的壳层半径内计算的，该壳层半径对应于1 个稳定的包合物的C—O 分布函数中的第1 个最小值，内插图为单个甲烷分子被水包围的示意图;(e) 引入电场强度为2 V/nm 时，水分子配位数降低率与频率的关系Fig.3.(a) Spatial distribution of water and methane molecules in the caged methane hydrate (up-right inset) and decomposed mixture (down-right inset) after THz-EM stimulus.(b) Change of the total number of hydrogen bonds in simulated system with time.(c) Relationship between the hydrogen bond loss rate of system and the external electric field frequency at an intensity of 2 V/nm.(d) Coordination number of water molecules for the methane in hydrate.It is calculated within a shell radius of 0.57 nm，which corresponds to the first minimum value in the C—O distribution function of a stable clathrate.The inset describes a single methane molecule surrounded by water.(e) Relationship between the reduction rate of the coordination number and the introduced field frequency at an intensity of 2 V/nm.

径向分布函数[5，7]也是分析水合物结构特征常用的评价参数，是晶体结构无序化的1 个衡量指标，其公式如下:

式中，V为模拟 体系的体积;Nα和Nβ分别为α粒子和β粒子的个数;niβ(r)表示以第i个粒子为中心，在r—r+Δr范围内β粒子的数目.甲烷与水分子之间的径向分布函数第一主峰对应的距离为笼型空腔中甲烷与水笼的距离，峰值越高，说明甲烷与水之间的范德瓦耳斯相互作用越强，则甲烷所受的范德瓦耳斯力越大，被捕捉的概率越大，越有利于水合物的稳定.图4 给出了相同体系在不同条件下的径向分布函数，可以看出，当加入特征频率的电磁波时，第1 个峰值显著减少，表明甲烷分子周围的水分子急剧减少，而由图4(b)可知，水合物中C 原子之间分布的第一壳层距离为0.65 nm，但当加入频率为10.3 THz 的电磁波时，在r=0.37 nm处出现第一壳层，说明此时甲烷从水笼中溢出，甲烷出现聚集现象，而其他非特征频率的电磁波刺激下该现象不明显.

图4 (a) 260 K 温度下，不同频率的太赫兹电磁刺激对水分子中O 原子相对于甲烷分子中C 原子的径向分布函数的影响;(b) 260 K 温度下，不同频率的太赫兹电磁刺激对甲烷分子中C 原子相对周围甲烷分子中的C 原子的径向分布函数的影响Fig.4.(a) Effect of THz-EM stimulation at different frequencies on the radial distribution function (RDF) of O atoms in water molecules w.r.t.the C atom in a methane molecule at 260 K;(b) effect of the stimulations on the RDF of C atoms in surrounding methane molecules w.r.t.the C atom in a methane molecule at 260 K.

在水合物的模拟研究中，F4有序参数曾被多次提到，它可以量化水合物从晶相到液相的转化.随着水合物的熔化，水分子的结晶度以及甲烷-水分子的配位逐渐发生变化，F4参数可以探测水分子的类似水合物的结晶度.F4统计的是水体系中0.35 nm[30]内两个相邻水分子最外侧氢原子以及水分子中氧原子组成的二面角.其公式如下:

其中Φi是由H—O···O—H 构成的二面角(如图5(a)).已有的结果表明，水合物、液态水、冰的F4有序参数分别为0.7，—0.04，—0.4[31-36].本研究发现，未加入电磁波的状态下，水合物保持稳定状态，因此F4值维持在0.7 附近;当加入频率为5 或40 THz的电磁波时，F4值从0.7 下降到0.5，说明此时水合物出现拉伸变形状态，但仍保持着晶体结构;而当加入频率为10.3 THz 的电磁波时，F4值从0.7下降到0 附近，表明此时水合物的晶体结构遭到破坏(如图5(b)).另外，不同条件下体系F4参数数值也符合水合物有序参数的取值范围.结合前面对氢键断裂的时间以及水合物的径向分布函数的研究可以看出，随着氢键网络的破坏，水合物晶体结构发生改变，甲烷分子从笼中脱离，水由晶体状态变为液态.随后，对比了不同强度下10.3 THz电磁波对水合物分解的影响，可以看出，随着电场强度的增加，水合物的分解速度加快，分子分解的数量增多，相应地，F4值也不断降低(如图5(c)和图5(d)).特别地，当电磁波的强度超过1.6 V/nm时，甲烷水合物可以在200 ns 内完全分解.

图5 (a) F4 值参数的示意图，红色虚线为氢键，字母H，O 表示两个水分子的氢和氧原子位置，水分子中氧原子间的距离在0.35 nm 内，多面体内两端的H—O···O—H 为距离最远的一对H，Φi 为扭转角;(b) 不同条件下F4 值随时间的变化量;(c) 频率为10.3 THz 的电磁波，不同强度下F4 值随时间的变化量;(d) 频率为10.3 THz，强度不同的电磁波作用下的F4 值Fig.5.(a) Schematic diagram of the F4 value parameter.The red dashed line denotes the hydrogen bond，and the H and O letters locates the hydrogen and oxygen atoms in two water molecules.The distance between two oxygen atoms is within 0.35 nm.The H pair at both ends in the H—O···O—H polyhedron accounts for the largest distance. Φi is the torsion angle.(b) Change of the F4 value with time under different conditions.(c) Variations of the F4 value with time under different filed intensities but the same frequency of 10.3 THz.(d) Relation between the F4 value and the EM field intensity at the same frequency of 10.3 THz.

4 结论

本文提出了一种利用太赫兹波特异性地且高效地分解可燃冰的机制，理论证明了在10.3 THz，1.6 V/nm 电磁场刺激下，可以在200 ns 内完全分解可燃冰的主要成分—甲烷水合物.物理机理为，甲烷水合物存在1 个与体相水不同的特征水吸收峰(中心频率为10.3 THz)，因此可以对该频率下的电磁场进行有效谐振吸收，进而水合物中的氢键网络遭到破坏，促进甲烷分子从水笼中脱离.本文分别利用氢键数目、甲烷周围水分子的配位数、有序参数F4等在引入太赫兹场前后的变化充分验证了该作用机制以及有效场频率的特异性.本文的研究有望提高未来对可燃冰的开采效率，降低天然气在运输贮存中的成本，推动新能源的发展[37].