環狀化合物-甲烷水合物穩定性的分子模擬

蘆文浩，梁海峰，王 帥，賈 菊

（太原理工大學 化學化工學院，山西 太原 030024）

天然氣（NG）作為一種可提供高效熱值、燃燒無污染和儲量豐富的能源，被公認為未來重要的戰略資源。以水合物形式存儲天然氣（主要成分為甲烷）的固化天然氣技術 （Solidified natural gas，SNG），是最具有潛力、安全環保、高效解決能源供應問題的儲運方式。

Ⅱ型水合物小籠晶穴可以容納大量CH4，大晶穴填充熱力學促進劑，穩定性提高，其中以環狀化合物為主要的促進劑，可使水合物在溫和條件下穩定存在。許多學者就實驗探索已證明其可行性：胡亞飛等[1]構建環戊烷（CP）+CH4水合物實驗體系，證明CP可提升水合物穩定存在溫度；孫強等[2]添加環狀醚類四氫呋喃（THF）提純煤層氣中的CH4，實驗驗證生成水合物壓力明顯降低；梁德青等[3]測定H2+CH4+環狀醚類四氫吡喃（THP）+水（H2O）、H2+CH4+CP+H2O兩種體系的水合物相平和數據，結果表明加入環狀化合物，降低了生成壓力；實驗只能提供直觀過程，無法洞悉微觀結構改變和作用機理。水合物分子動力學模擬可為多組分體系中分析分子運動、作用和結構改變等微觀角度的探究提供有效途徑[4]，對此研究也取得一些成果：梅東海等[5]首次采用分子動力學計算H型水合物，客體分子填充CH4和環辛烷，得出H型水合物3種籠形晶穴都需要有客體分子占據，水合物才能穩定存在；Wu[6]利用GROMACS模擬軟件，計算CH4+THF混合客體分子水合物成核的分子動力學，結果說明體系形成的籠形簇團相對單組份更穩定；張慶東等[7-8]采用分子動力學模擬填充不同占有率的CH4和THF對水合物穩定性影響，結果表明水合物占有率決定水合物穩定性，隨著溫度降低，利于單組份CH4水合物穩定；殷開梁等[9]比較 Material studio（MS）軟件中各分子力場對于氫鍵作用描述的適用性，結果表明該軟件可以采用一般的非鍵作用近似反映氫鍵。耿春宇等[10]通過分子動力學模擬CH4對應Ⅰ型水合物，得出結論：因CH4占據晶穴，水合物晶穴的相互作用能降低，水分子間氫鍵作用加強，穩定性提高，且占據率高于87.5%水合物才能表現較好穩定性；丁麗穎等[11]模擬填充CH4的Ⅰ型水合物穩定性受晶穴占有率和溫度的影響，結果表明較低占有率（<0.625）和較高溫度（>320k）條件下，水合物基本瓦解。

綜上所述，有關環狀化合物-CH4水合物穩定性研究多數集中在實驗，少有利用分子動力學探尋CP和THP對Ⅱ型天然氣水合物的穩定作用，且高占有率利于水合物穩定性；本文采用Material Studio 8.0軟件，搭建占有率100%分別添加CP、THF和THP的二元體系CH4水合物，從分子角度比較其穩定性和解釋作用機理，為SNG技術實際應用提供參考。

1 模型搭建和模擬方法

1.1 籠形單元組成

水合物主體框架是以水分子氫鍵結合形成的一種空間點陣結構[12]，Ⅱ型水合物原始單晶胞包含136個水分子，構成16個小籠晶穴（512）和8個大籠晶穴(51264)。小籠晶穴由CH4分子占據，大籠晶穴分別由CP、THF和THP占據。客體分子手動搭建，并進行能量幾何優化，使其最接近真實狀態。表1為依據胡盛志等[13]優化方法計算得到的客體分子范德華模型相關參數，籠形晶穴占據結構如圖1所示。

表1 客體分子范德華模型參數

圖1 客體分子占據籠形晶穴結構

1.2 模型搭建

在周期邊界條件下，Ⅱ型水合物原始單晶胞在溫度T=123K下，晶格參數a=b=c=1.7175nm。晶胞中心包含對稱存在的四面體，與周邊6個十二面體形成Ⅱ型水合物主體框架，其氧原子坐標來自X射線衍射實驗[14]，氫原子分布滿足“冰規則”[15]。客體分子與主體框架間為范德華作用力（Van der Waals）[16]，模擬體系構建2×2×2的Ⅱ型Super cell晶胞，該體系中1088個H2O構建128個小籠晶穴（512）和64個大籠晶穴(51264)的主體框架。小籠晶穴填充128個CH4，大籠晶穴分別填充64個CP、THF和THP，本文3種模擬體系分別命名為A、B、C，便于觀察，將體系H原子隱藏，初始結構模型如圖2所示。

圖2 水合物模擬體系初始結構

1.3 模擬方法

先對模擬體系進行能量優化和幾何優化，選用CVFF（一致性價力場）力場描述模擬體系，SPC（Simple point charge）模型[17]描述 H2O 分子，將其視作剛性分子，H-O-H平面夾角為104.52°，H-O鍵長為固定值0.09570nm，水分子和客體分子、客體分子間非鍵作用使用Lennard-Jones作用勢描述，在MS中采用截斷球近似代替，截斷半徑為1.5nm，Ewald加和法處理長程靜電作用力，精度設置為4.1868×103J/mol。

對幾何優化結束的體系進行結構松弛 ，選水分子中的O原子固定，控溫方法使用Nose-Hoover，NVT（正則系綜）模擬設定T分別為273K、283K、293K，分子初始速依據Maxwell-Boltzmann隨機產生，模擬時間為 50ps（1ps=10-12s），步長為 1fs（1fs=10-15s），每5000步輸出一幀結構圖；解除O原子固定，在NPT（等溫恒壓）系綜下進行動力學模擬，控溫方式不變，采用Berendsen控壓方式調整體系壓力值，壓力設為2MPa，模擬時間為400ps，模擬結束后分析二元體系水合物結構。

1.4 模擬最終結構

將軟件輸出的最終結構進行比較分析，如圖3所示。

圖3 水合物模擬體系最終結構

由圖3水合物模擬體系最終結構看出：A、B可維系原有結構，氫鍵排序重疊度較高，客體分子仍能與籠形晶穴中心相對重合；壓力不變，隨著溫度升高，出現部分晶穴的扭曲加深，B較A在293K下水合物結構破壞嚴重；C模型體系在273K時完整度與相同條件下的A、B相比較差，部分客體分子因晶穴的扭曲加重位移較大，氫鍵分布出現小范圍的雜亂，隨著溫度的升高結構出現了不同程度的扭曲，283K溫度下水合物結構基本破壞，殘存的氫鍵不足以維系體系的穩定存在，客體分子位移較大，出現小范圍的聚集，293K條件下氫鍵雖然存在，但晶穴有序結構不存在，客體分子逃逸，出現聚集現象，預判此時從固體轉變為液體。

2 模擬結果分析

2.1 徑向分布函數（Radial Distribution Function，RDF）

穩定的水合物主體框架水分子中的O原子RDF特征是：呈現多峰值曲線狀態，其峰高隨著O原子的增加逐漸降低。可依此判斷其無序狀態，也代表水合物區域密度和平均密度的比值；液態水分子中O原子go-o(r)在r值大于0.278nm時就趨向為1[18]，即區域密度等于平均密度表現為液態流動，將A、B、C三種模型對應壓力溫度值下的go-o(r)值對比分析，如圖4所示。

圖4 水合物O原子RDF分布

由圖4可知O原子RDF分布在r=0.278nm時出現第一峰高，代表相鄰水分子最近的O原子間距離，第二特征峰高在r=0.452nm出現，滿足晶體結構RDF特征；T=273K時，A、B、C呈現出穩定水合物晶體的特征峰，但都出現了峰高降低和峰谷升高的現象，說明有部分O原子發生位移，即部分籠形晶穴扭曲，但整體框架結構能夠穩定存在，且第二峰高峰值ABC，則A穩定性要高于相同條件下的B，C模型峰谷升高，穩定性較差；T=283K時A同樣略優于B，C模型在r=0.278nm出現峰值后，峰谷升高，第二特征峰消失，此時與液體特征峰相似[15]；T=293K時，A和B模型RDF分布重合度相近，C模型RDF分布較T=283K對應分布：峰谷升高加劇，與液體特征峰分布更為接近，說明破壞程度更大，客體分子逃逸并發生聚集。

2.2 均方位移（MSD）

MSD表示指定分子位置與初始位置的偏離程度，依此可以判斷體系是固態還是液態，MSD分布圖像斜率K越大表示偏離初始程度越大。圖5為A、B、C對應壓力溫度值下的MSD圖像。

圖5為模擬體系所有條件下的O原子MSD分布，為更清晰比較，去除C模型283K、293K條件，呈現圖6。圖6可分析出C模型在283K、293K條件下，MSD隨著模擬時間增大，與液態水分子MSD斜率K相似，表明這兩種狀態下無序化程度嚴重，結構破壞變成了液態溶液；從圖6、7綜合比較，相同p值下，較穩定的模型體系MSD圍繞在接近零值的地方波動，結果 KC(293K)KB(293K)KA(293K)，KC(283K)KB(283K)KA(283K)，KC(273K)KB(273K)KA(273K)，說明對添加同種環狀化合物促進劑，隨著溫度降低，有利于水合物的穩定性，且體系在對應模擬條件下呈現KAKBKC；

圖5 水合物O原子均方位移

圖6 水合物O原子均方位移局部放大圖

2.3 相互作用能

分子動力學模擬結束后可得到模擬體系的總能量Etotal，主體框架能量Ecages以及填充的所有客體分子能量Eguests，從而得到主體框架和客體分子間的相互作用能Ebinding，滿足公式：

從能量角度分析，相互作用能越小，則體系最終結構越穩定，圖7為A、B、C對應溫度壓力值下整個體系的相互作用能。

從圖7可看出，A模型E(293K)E(283K)E(273K)，B、C模型呈現相同規律即低溫有利于水合物結構穩定；相互作用能越低，水合物結構越穩定，A模型整體穩定性較好，與MSD分析結果相一致；從分子角度分析3種促進劑對穩定性的影響，已知客體分子與晶穴直徑比值在0.76～1之間[19]，可維系穩定水合物結構，環狀化合物的填充，起到支撐水合物框架的作用。從表1可看出Ⅱ型大晶穴，直徑比值CP略高于THF，更接近1，穩定性得到提升，THP與晶穴比值優于CP和THF，理論上應相對更穩定，但模擬結果高溫區穩定性較差，原因分析：THP分子過大，高溫狀態下分子活躍程度較大，其分子組成中的O原子會與周邊水分子框架形成氫鍵作用，使得與周邊其他水分子間范德華作用力不均衡，客體分子牽動晶穴單元發生位移（氫鍵作用強于范德華作用力），從而使得模擬體系穩定性下降。

圖7 模擬體系不同溫度下相互作用能

3 結論

在NPT系綜下分子動力學模擬，本文通過最終結構、RDF、MSD以及相互作用能對比分析，得出以下結論：

（1）環狀化合物作為客體分子占據甲烷水合物大晶穴，能在溫和條件下維系水合物穩定存在。隨著溫度升高，模擬體系穩定性下降；

（2）CP模擬體系穩定性優于THF模擬體系，說明環狀化合物與Ⅱ型水合物大晶穴的直徑比越接近1，水合物越穩定；

（3）THF模擬體系穩定性優于THP模擬體系，說明在高溫條件下，當較大分子尺寸的環狀化合物分子含有O原子時，環狀化合物與水合物框架產生不平衡氫鍵作用，使水合物穩定性降低。

符號說明

a-晶格參數，nm；b-晶格參數，nm；c-晶格參數，nm；Ebinding-相互作用能，J/mol；Ecages-空晶穴能量，J/mol；Eguests-總客體分子能 J/mol；Etotal-模型總能量，J/mol；E（T）-對應溫度下相互作用能，J/mol；K-MSD分布斜率，m2/s；MSD—O 原子均方位移，m2；N-分子總數， 個；p-壓力，MPa；r-截斷半徑，nm；RDF-徑向分布函數，量綱為1；T-溫度，K。