Ⅰ型CH4/CO2水合物穩(wěn)定性的分子動(dòng)力學(xué)模擬

宋 軍，徐 梁，王蕊蕊，趙偉龍，趙晨旭

（河南理工大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，河南 焦作 454003）

在永久凍土區(qū)和深海海底蘊(yùn)含著大量的天然氣水合物。如我國(guó)南海地區(qū)，天然氣水合物資源豐富。合理地開(kāi)發(fā)利用天然氣水合物資源，將成為解決我國(guó)未來(lái)能源短缺問(wèn)題的重要途徑[1]。水合物包含3種常見(jiàn)的結(jié)構(gòu)：Ⅰ型體心立方結(jié)構(gòu)、Ⅱ型面心立方結(jié)構(gòu)和H型六方結(jié)構(gòu)[2]。1 m3天然氣水合物可能含有多達(dá)180 m3的氣體（標(biāo)準(zhǔn)溫度和壓力），水合物中所蘊(yùn)含的碳含量是全球化石燃料總碳量的兩倍[3]。

采用CO2置換CH4水合物的方法既能收集CH4，又能儲(chǔ)存CO2，近年來(lái)被廣泛研究。余靜薇等[4]研究發(fā)現(xiàn)Ar作為小分子氣體添加劑促進(jìn)了CO2置換CH4水合物，提高了CH4的置換率，在置換反應(yīng)的初始階段CH4的置換率增長(zhǎng)較快，在置換反應(yīng)中后期CH4的置換率變化平緩。孫建業(yè)等[5]通過(guò)實(shí)驗(yàn)發(fā)現(xiàn)置換反應(yīng)包含快速反應(yīng)和緩慢反應(yīng)兩個(gè)階段，表層水合物分解過(guò)程影響快速反應(yīng)階段的CH4產(chǎn)氣，孔隙水的擴(kuò)散過(guò)程影響快速反應(yīng)階段的CH4產(chǎn)氣。宋光春等[6]研究了溫度、壓力對(duì)置換的影響，發(fā)現(xiàn)溫度對(duì)置換結(jié)果的影響大于壓力對(duì)置換結(jié)果的影響。

分子動(dòng)力學(xué)模擬（MD）基于分子動(dòng)力學(xué)理論，其通過(guò)力場(chǎng)得到粒子運(yùn)動(dòng)軌跡、系統(tǒng)動(dòng)態(tài)和熱力學(xué)統(tǒng)計(jì)數(shù)據(jù)以及分子結(jié)構(gòu)性質(zhì)。分子模擬計(jì)算精度高，與實(shí)驗(yàn)數(shù)據(jù)往往高度吻合并且能獲得實(shí)驗(yàn)無(wú)法獲得的數(shù)據(jù)，被廣泛應(yīng)用于多種學(xué)科的研究。一些學(xué)者利用分子動(dòng)力學(xué)模擬方法從微觀角度對(duì)水合物結(jié)構(gòu)特性進(jìn)行了研究。張慶東等[7]發(fā)現(xiàn)溫度對(duì)水合物結(jié)構(gòu)穩(wěn)定性的影響較大，其中溫度對(duì)大孔穴結(jié)構(gòu)穩(wěn)定性的影響較大，而對(duì)小孔穴結(jié)構(gòu)穩(wěn)定性的影響較小。Kusalik等[8]通過(guò)MD方法研究了水合物生長(zhǎng)速率，發(fā)現(xiàn)水合物最大生長(zhǎng)速率高于冰生長(zhǎng)速率的4倍，且Ⅰ型水合物界面支持Ⅱ型水合物的生長(zhǎng)。Wu等[9]對(duì)CH4、四氫呋喃（THF）和CH4+THF水合物的成核過(guò)程進(jìn)行了分子動(dòng)力學(xué)模擬，結(jié)果發(fā)現(xiàn)CH4+THF混合客體體系的誘導(dǎo)時(shí)間明顯短于純CH4和純THF體系。Sun等[10]利用MD方法，研究了水合物大小籠空缺對(duì)結(jié)構(gòu)穩(wěn)定性的影響，發(fā)現(xiàn)大籠缺失客體分子比小籠缺失客體分子對(duì)結(jié)構(gòu)穩(wěn)定性的影響更大。蘆文浩等[11]利用MD方法，研究了273 K、283 K和293 K，壓力2 MPa條件下環(huán)狀化合物作為添加劑對(duì)CH4水合物的影響，發(fā)現(xiàn)同等條件下添加環(huán)戊烷體系最穩(wěn)定。賈菊等[12]利用MD方法，研究了環(huán)戊烷、THF和四氫吡喃對(duì)煤層氣水合物結(jié)構(gòu)的影響，結(jié)果發(fā)現(xiàn)263 K、10 MPa時(shí)，環(huán)戊烷可以很大程度地提高水合物結(jié)構(gòu)的穩(wěn)定性。易莉芝等[13]利用MD方法，研究了NaCl溶液中CO2置換CH4水合物的過(guò)程，發(fā)現(xiàn)在NaCl溶液中被置換出的CH4分子易形成氣泡，溶液相和水合物相形成CH4濃度梯度，水合物發(fā)生分解。

由于宏觀實(shí)驗(yàn)無(wú)法實(shí)現(xiàn)CH4水合物和CO2水合物在不同晶穴占有率條件下的對(duì)比分析，同時(shí)模擬實(shí)驗(yàn)也較少研究晶穴占有率對(duì)置換過(guò)程中CO2水合物和CH4水合物結(jié)構(gòu)穩(wěn)定性的影響，為更深入地了解CO2置換CH4水合物的機(jī)理，本文應(yīng)用MD方法，采用蒙特卡洛吸附模塊分別建立CH4水合物模型和CO2水合物模型，在一定溫度壓力條件下，分析不同晶穴占有率對(duì)CH4水合物和CO2水合物結(jié)構(gòu)穩(wěn)定性的影響。通過(guò)結(jié)構(gòu)變化、徑向分布函數(shù)、均方位移和勢(shì)能變化等方面，從分子角度分析兩者的變化，為水合物置換技術(shù)應(yīng)用提供參考。

1 模型的建立及模擬方法

1.1 水合物模型

Ⅰ型水合物模型屬于體心立方結(jié)構(gòu)，建立超晶胞2 × 2 × 2模型，晶格參數(shù)2.406 nm × 2.406 nm×2.406 nm，空籠結(jié)構(gòu)由368個(gè)水分子組成?；\型結(jié)構(gòu)水分子中氧原子的初始坐標(biāo)來(lái)源于X射線單晶衍射實(shí)驗(yàn)[14]，利用Bernal-Fowler規(guī)則[15]對(duì)氫原子進(jìn)行排列。水合物空籠的初始模型如圖1所示，吸附客體分子采用蒙特卡洛方法吸附，100%占有率籠型吸附位點(diǎn)如圖2中紅色標(biāo)記所示，客體分子主要是CH4和CO2。用此方法分別建立不同晶穴占有率（12.5%、25.0%、37.5%、50.0%、62.5%、75.0%、87.5%和100%）的CH4水合物與CO2水合物。

圖1 水合物空籠初始構(gòu)型

圖2 籠中晶穴吸附位點(diǎn)

1.2 模擬方法

分子動(dòng)力學(xué)計(jì)算應(yīng)用牛頓運(yùn)動(dòng)定律，獲取各時(shí)間下粒子運(yùn)動(dòng)狀態(tài)的物理參量，通過(guò)位移、速度、加速度的變化情況來(lái)了解分子運(yùn)動(dòng)的微觀細(xì)節(jié)[16]。模型選用CVFF力場(chǎng)進(jìn)行計(jì)算，系綜選擇NPT，水分子模型選用SPC模型，描述長(zhǎng)程靜電作用力和范德瓦爾斯作用力采用Ewald方法，Lennard-Jones作用勢(shì)描述分子間作用力，截?cái)喟霃?.5 nm，加和法計(jì)算精度4.19 × 10-6kJ/mol。

模擬之前，先對(duì)模型進(jìn)行能量?jī)?yōu)化和結(jié)構(gòu)優(yōu)化，而后將籠型結(jié)構(gòu)中氧原子坐標(biāo)進(jìn)行固定，使用NVT系綜，在253 K條件下進(jìn)行50 ps的馳豫，控溫方法使用Nose-Hoover法。為了減小計(jì)算過(guò)程中的誤差，設(shè)置馳豫的時(shí)間步長(zhǎng)為1 fs，每間隔1000步輸出一張結(jié)構(gòu)圖。馳豫結(jié)束后解除氧原子固定，將晶穴占有率不同的模型在NPT系綜中進(jìn)行分子動(dòng)力學(xué)模擬，依然設(shè)置溫度253 K、壓力2 MPa、模擬時(shí)間400 ps，溫控方式和結(jié)構(gòu)圖輸出間隔不變。

2 結(jié)果與討論

2.1 構(gòu)象分析

將模擬之后的結(jié)構(gòu)圖進(jìn)行匯總，選取晶穴占有率為37.5%、50.0%和62.5%的CH4水合物和CO2水合物進(jìn)行對(duì)比，如圖3所示。初始的水合物籠型結(jié)構(gòu)具有高度對(duì)稱(chēng)性，因此結(jié)構(gòu)圖的部分氧原子和氫鍵被完全遮擋，根據(jù)氧原子、碳原子的位移偏離程度可判斷籠型結(jié)構(gòu)的變形程度。由圖3可知，晶穴占有率為50.0%和62.5%的CO2水合物的氫鍵和氧原子排序重疊度較高，客體分子仍處于籠中心位置。晶穴占有率為37.5%的CO2水合物氧原子和氫鍵位置發(fā)生輕微偏移，但整體結(jié)構(gòu)相對(duì)穩(wěn)定。晶穴占有率為62.5%的CH4水合物氫鍵和氧原子發(fā)生了大的位移偏離，籠型結(jié)構(gòu)發(fā)生破裂。晶穴占有率為37.5%、50.0%的CH4水合物相比晶穴占有率為62.5%的CH4水合物，氧原子出現(xiàn)更大的位移偏離，氫鍵的排序混亂，水合物籠型結(jié)構(gòu)完全解離。

圖3 不同晶穴占有率的CO2水合物（(a)～(c)）和CH4水合物（(d)～(f)）最終構(gòu)象

2.2 徑向分布函數(shù)

徑向分布函數(shù)（RDF）是分析結(jié)構(gòu)無(wú)序性的重要方法。RDF中g(shù)α-β表示在粒子α相距r處粒子β出現(xiàn)的概率，其也能利用區(qū)域密度與平均密度的比值表示。對(duì)于液態(tài)水分子，其氧原子的徑向分布逐漸趨近于1，表示水分子處于無(wú)序狀態(tài)。如果水合物結(jié)構(gòu)處于穩(wěn)定狀態(tài)，則水分子的氧原子徑向分布曲線呈現(xiàn)多個(gè)峰值，并且這些峰的峰高隨著r值的增加而逐漸減小。第二特征峰消失或者峰值降低，說(shuō)明水合物處于不穩(wěn)定狀態(tài)。

圖4是不同晶穴占有率下CH4水合物和CO2水合物的RDF曲線。由圖4可知，RDF曲線在r=0.278 nm附近出現(xiàn)第一特征峰，代表籠型結(jié)構(gòu)中兩個(gè)氧原子之間的最近距離，在r=0.452 nm附近出現(xiàn)第二特征峰，代表氧原子之間相隔一個(gè)或者多個(gè)原子的距離，滿(mǎn)足穩(wěn)定結(jié)構(gòu)的RDF。隨著晶穴占有率增加，峰值升高、峰谷降低，說(shuō)明水分子無(wú)序性降低，水合物穩(wěn)定性增加。對(duì)比不同晶穴占有率的CO2水合物，RDF特征峰區(qū)別不明顯，晶穴占有率為62.5%的體系結(jié)構(gòu)穩(wěn)定性?xún)?yōu)于37.5%和50.0%的體系。相同工況下，對(duì)比不同晶穴占有率的CH4水合物的徑向分布，CH4分子晶穴占有率越小第一特征峰越低，第二特征峰也越低或消失，水合物結(jié)構(gòu)無(wú)序程度增加。在相同工況下，晶穴占有率在37.5%與62.5%之間時(shí)，CO2水合物的RDF特征峰均明顯優(yōu)于CH4水合物。

圖4 不同體系水合物氧原子RDF曲線

2.3 均方位移

在分子模擬計(jì)算過(guò)程中，原子在不停運(yùn)動(dòng)，其位置在不斷變化，會(huì)相對(duì)于原始位置偏離一定距離。均方位移（MSD）用公式表示如下：

式中，N為系統(tǒng)中粒子總數(shù)；Ri(0)為i粒子0時(shí)刻的位置；Ri(t)為i粒子t時(shí)刻的位置。

MSD可以很清楚地反映粒子相對(duì)于初始位置的偏離程度。當(dāng)體系為液態(tài)或者氣態(tài)時(shí)，粒子所受的約束變小，此時(shí)的MSD曲線是一條隨時(shí)間增長(zhǎng)而上升的曲線，對(duì)應(yīng)于水合物分解狀態(tài)。穩(wěn)定的水合物結(jié)構(gòu)，粒子受到的約束大，MSD值會(huì)在0值附近較小波動(dòng)，曲線整體接近一條水平線。

圖5為CH4水合物不同晶穴占有率時(shí)氧原子的MSD曲線。由圖5可知，晶穴占有率為75.0%、87.5%和100%的CH4水合物氧原子MSD值均很小，整條曲線接近水平，氧原子偏離程度很小，說(shuō)明水合物籠型結(jié)構(gòu)是穩(wěn)定的。CH4水合物晶穴占有率低于75.0%后，占有率越低MSD值隨時(shí)間增長(zhǎng)的速率越快，水分子的無(wú)序程度越大，籠型結(jié)構(gòu)解離，不能穩(wěn)定存在。

圖5 CH4水合物氧原子MSD曲線

圖6為相同溫壓時(shí)晶穴占有率為37.5%、50.0%和62.5%的CH4水合物和CO2水合物氧原子的MSD曲線。由圖6可知，相對(duì)于CO2水合物，CH4水合物在晶穴占有率為62.5%時(shí)的MSD值較小，但仍然高于同種占有率的CO2水合物，說(shuō)明此時(shí)CH4水合物籠型結(jié)構(gòu)無(wú)序化程度高于CO2水合物。占有率為37.5%、50%和62.5%時(shí)，CO2水合物的氧原子MSD值較為穩(wěn)定，圍繞0值附近波動(dòng)，籠型結(jié)構(gòu)未發(fā)生明顯扭曲。因此，在一定溫度、壓力和晶穴占有率范圍內(nèi)，CO2可以形成穩(wěn)定的CO2水合物，而CH4卻不易形成穩(wěn)定的CH4水合物，CH4水合物分解釋放CH4而CO2仍然可以穩(wěn)定存在于水合物中，從而達(dá)到CO2置換CH4的目的。

圖6 不同體系水合物氧原子MSD曲線

2.4 勢(shì)能變化

在分子模擬過(guò)程中，體系的勢(shì)能有波動(dòng)性，對(duì)于穩(wěn)定的水合物結(jié)構(gòu)，勢(shì)能在一定范圍平穩(wěn)波動(dòng)。圖7為晶穴占有率為37.5%、50.0%和62.5%的CH4水合物和CO2水合物在模擬過(guò)程中勢(shì)能隨時(shí)間變化的曲線。由圖7可知，晶穴占有率為37.5%的CH4水合物在140 ps附近勢(shì)能從-15200 kJ/mol迅速升高到-13600 kJ/mol，水合物結(jié)構(gòu)發(fā)生了快速分解，水分子的無(wú)序程度迅速升高。晶穴占有率為50.0%和62.5%的CH4水合物分別在180 ps附近和300 ps附近勢(shì)能升高，結(jié)構(gòu)同樣也發(fā)生了解離。CH4水合物晶穴占有率越低，結(jié)構(gòu)分解的時(shí)間越早，分解的速率也越快。而相對(duì)于相同占有率的CH4水合物，CO2水合物的勢(shì)能在-15200 kJ/mol附近較小波動(dòng)，表明體系整體結(jié)構(gòu)較為穩(wěn)定。

圖7 不同體系水合物勢(shì)能隨時(shí)間的變化曲線

圖8顯示了CO2水合物勢(shì)能隨時(shí)間的變化。由圖8可知，CO2水合物晶穴占有率低于37.5%時(shí)，勢(shì)能在不同時(shí)間都升高，結(jié)構(gòu)發(fā)生了解離；晶穴占有率在37.5%以上時(shí)勢(shì)能波動(dòng)小，結(jié)構(gòu)較為穩(wěn)定。

圖8 CO2水合物勢(shì)能隨時(shí)間的變化曲線

3 結(jié)論

利用MD及蒙特卡洛吸附法建立了CH4水合物和CO2水合物模型，在NPT系綜下研究了其穩(wěn)定性，分析了最終構(gòu)象、徑向分布函數(shù)、均方位移以及勢(shì)能的變化，得出以下結(jié)論：

（1）在253 K、2 MPa條件下，Ⅰ型CO2水合物晶穴占有率在37.5%以上時(shí)，水分子氫鍵組成的籠型結(jié)構(gòu)整體能保持穩(wěn)定，并且CO2的占有率越高水合物籠型結(jié)構(gòu)越穩(wěn)定。

（2）相同晶穴占有率時(shí)，CO2水合物比CH4水合物結(jié)構(gòu)穩(wěn)定性更強(qiáng)。晶穴占有率在37.5%與62.5%之間時(shí)，CH4水合物均發(fā)生了不同程度的分解并且晶穴占有率越低越先發(fā)生分解。此時(shí)，CO2水合物結(jié)構(gòu)雖有不同程度的輕微扭曲，但整體結(jié)構(gòu)依然保持穩(wěn)定。

（3）在一定溫度壓力和晶穴占有率范圍內(nèi)，CO2水合物比CH4水合物穩(wěn)定性好，解釋了CO2能置換CH4水合物的原因，即CH4水合物發(fā)生分解釋放CH4，而CO2穩(wěn)定存在于水合物中，從而實(shí)現(xiàn)收集CH4并儲(chǔ)存CO2的目的。