太赫茲電磁波調控甲烷水合物分解*

朱智 閆韶健 段銅川 趙妍 孫庭鈺 李陽梅

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].