環(huán)狀化合物-甲烷水合物穩(wěn)定性的分子模擬

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

（太原理工大學(xué) 化學(xué)化工學(xué)院，山西 太原 030024）

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

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

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

1 模型搭建和模擬方法

1.1 籠形單元組成

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

表1 客體分子范德華模型參數(shù)

圖1 客體分子占據(jù)籠形晶穴結(jié)構(gòu)

1.2 模型搭建

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

圖2 水合物模擬體系初始結(jié)構(gòu)

1.3 模擬方法

先對(duì)模擬體系進(jìn)行能量?jī)?yōu)化和幾何優(yōu)化，選用CVFF（一致性價(jià)力場(chǎng)）力場(chǎng)描述模擬體系，SPC（Simple point charge）模型[17]描述 H2O 分子，將其視作剛性分子，H-O-H平面夾角為104.52°，H-O鍵長(zhǎng)為固定值0.09570nm，水分子和客體分子、客體分子間非鍵作用使用Lennard-Jones作用勢(shì)描述，在MS中采用截?cái)嗲蚪拼?，截?cái)喟霃綖?.5nm，Ewald加和法處理長(zhǎng)程靜電作用力，精度設(shè)置為4.1868×103J/mol。

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

1.4 模擬最終結(jié)構(gòu)

將軟件輸出的最終結(jié)構(gòu)進(jìn)行比較分析，如圖3所示。

圖3 水合物模擬體系最終結(jié)構(gòu)

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

2 模擬結(jié)果分析

2.1 徑向分布函數(shù)（Radial Distribution Function，RDF）

穩(wěn)定的水合物主體框架水分子中的O原子RDF特征是：呈現(xiàn)多峰值曲線狀態(tài)，其峰高隨著O原子的增加逐漸降低?？梢来伺袛嗥錈o(wú)序狀態(tài)，也代表水合物區(qū)域密度和平均密度的比值；液態(tài)水分子中O原子go-o(r)在r值大于0.278nm時(shí)就趨向?yàn)?[18]，即區(qū)域密度等于平均密度表現(xiàn)為液態(tài)流動(dòng)，將A、B、C三種模型對(duì)應(yīng)壓力溫度值下的go-o(r)值對(duì)比分析，如圖4所示。

圖4 水合物O原子RDF分布

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

2.2 均方位移（MSD）

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

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

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

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

2.3 相互作用能

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

從能量角度分析，相互作用能越小，則體系最終結(jié)構(gòu)越穩(wěn)定，圖7為A、B、C對(duì)應(yīng)溫度壓力值下整個(gè)體系的相互作用能。

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

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

3 結(jié)論

在NPT系綜下分子動(dòng)力學(xué)模擬，本文通過(guò)最終結(jié)構(gòu)、RDF、MSD以及相互作用能對(duì)比分析，得出以下結(jié)論：

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

（2）CP模擬體系穩(wěn)定性優(yōu)于THF模擬體系，說(shuō)明環(huán)狀化合物與Ⅱ型水合物大晶穴的直徑比越接近1，水合物越穩(wěn)定；

（3）THF模擬體系穩(wěn)定性優(yōu)于THP模擬體系，說(shuō)明在高溫條件下，當(dāng)較大分子尺寸的環(huán)狀化合物分子含有O原子時(shí)，環(huán)狀化合物與水合物框架產(chǎn)生不平衡氫鍵作用，使水合物穩(wěn)定性降低。

符號(hào)說(shuō)明

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