TBAB半籠型甲烷水合物穩(wěn)定性分子動力學(xué)模擬

張 強，梁海峰，吉夢霞，楊 琨

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

隨著常規(guī)化石能源的逐漸枯竭，人類對新能源的需求量越來越大，天然氣水合物作為一種非常規(guī)天然氣，資源總量非常龐大，估計為3.1 × 1015～7.6 ×1018m3，因此其開發(fā)利用得到了廣泛的關(guān)注。天然氣水合物是在低溫、高壓條件下烴類氣體分子和水形成的類冰狀晶體化合物，主體分子（水分子）通過氫鍵構(gòu)成籠型骨架，客體分子（烴類分子）填充在晶穴中，主客體分子間的作用力為范德華力。

目前發(fā)現(xiàn)的水合物晶體包括Ⅰ型、Ⅱ型和H型，自然界中存在的水合物以Ⅰ型和Ⅱ型為主，但大約只有33%的晶穴被客體分子占據(jù)[1]。 Saim等[2]指出甲烷水合物在純水體系中的客體分子晶穴占有率低于22.2%；Saruprias等[3]模擬了CO2水合物的分解過程，指出溫度和客體分子晶穴占有率是影響水合物穩(wěn)定性的重要因素；Sugahara等[4]指出了晶穴占有率對水合物穩(wěn)定性影響的規(guī)律；Yu等[5]研究表明水合物籠中CH4濃度越高， 水合物由亞穩(wěn)態(tài)結(jié)構(gòu)向穩(wěn)定結(jié)構(gòu)的轉(zhuǎn)變速度越快；蘆文浩等[6]研究了不同促進劑下甲烷水合物的穩(wěn)定性，表明溫度較低、客體分子與晶穴的直徑比接近1時，水合物穩(wěn)定性較好；耿春宇等[7]采用分子動力學(xué)方法模擬Ⅰ型甲烷水合物，指出由于客體分子占據(jù)水合物晶穴，晶穴間的相互作用能降低，水分子間的氫鍵作用加強，水合物的穩(wěn)定性提高， 并且在晶穴占有率高于87.5%的情況下，水合物表現(xiàn)出較好的穩(wěn)定性。 由此表明，較低的溫度和較高的晶穴占有率能提高水合物的穩(wěn)定性。

Fowler等[8]發(fā)現(xiàn)較大的客體分子破壞水合物籠型結(jié)構(gòu)后，占據(jù)多個半籠型晶穴，形成特殊的水合物，稱為半籠型水合物。 四丁基溴化銨（TBAB）水合物是最常見的一種半籠型水合物，由于其具有良好的可降解性和低毒性，能有效加快反應(yīng)速度和降低反應(yīng)條件，且不容易揮發(fā)，因此受到廣泛關(guān)注和研究。 Shimada等[9]發(fā)現(xiàn)TBAB水合物存在A型（TBAB·26H2O）和B型（TBAB·38H2O）兩種結(jié)構(gòu)，512晶穴能捕獲氣體小分子（CH4、N2、CO2等），形成氣體水合物。Arjmandi等[10]和Mohammadi等[11]測量了H2/N2/CH4/天然氣+ TBAB水溶液的相平衡數(shù)據(jù)。 Zhen等[12]研究表明， 由于TBAB半籠型水合物只有512晶穴能捕獲客體分子， 因此TBAB半籠型水合物的選擇性優(yōu)于Ⅰ型、Ⅱ型和H型水合物；Aladko等[13]指出，TBAB水合物中最穩(wěn)定的水合數(shù)為38H2O。 目前，TBAB半籠型水合物的相平衡實驗已得到充分研究，但由于晶體結(jié)構(gòu)的復(fù)雜性，其微觀結(jié)構(gòu)的動力學(xué)內(nèi)容尚未有研究，且水合物客體分子的晶穴占有率在實驗過程中難以人為控制， 因此本文基于分子動力學(xué)模擬，從微觀角度分析了B型TBAB水合物的穩(wěn)定性機理。

1 模型構(gòu)建和模擬方法

1.1 模型構(gòu)建

TBAB·38H2O水合物單晶胞屬于正交晶系，空間群為pmma， 水分子中氧原子和TBAB中原子的初始坐標(biāo)根據(jù)XRD結(jié)果確定[14]，水分子和溴離子構(gòu)成水合物籠，氫原子無序排列，且符合Bernal-Fowler規(guī)則[15]。模擬體系由2 × 2 × 2個TBAB半籠型水合物單晶胞構(gòu)成， 在x、y、z坐標(biāo)方向上采用周期性邊界條件，晶胞尺寸為4.2120 nm × 2.5286 nm × 2.4036 nm，體系包含608 個水分子，16 個四丁基氨陽離子（TBA+）和16個溴離子（Br-），48個512晶穴、32個51262晶穴和32個51263晶穴， 其中TBA+填充在由51262晶穴和51263晶穴形成的半籠晶穴中，512晶穴由CH4分子填充。 本文通過分子動力學(xué)模擬研究不同晶穴占有率（θ=100%、75%、50%、25%）在不同溫度下對TBAB半籠型水合物穩(wěn)定性的影響，模型初始結(jié)構(gòu)如圖1所示。

1.2 模擬方法

運用Materials Studio 8.0軟件，采用COMPASSⅡ力場 （Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies） 在正則系綜（NVT）下進行模擬， 使用Ewald加和法描述長程靜電作用力和范德瓦爾斯作用力，SPC（Simple Point Charge）模型描述水分子間相互作用，其中水分子視為剛性分子，鍵長和鍵角的值固定，Lennard-Jones勢能函數(shù)描述水分子和其它分子間非鍵相互作用。

在進行分子動力學(xué)模擬之前，先對模型采用最速下降法和共軛梯度法進行幾何優(yōu)化，使結(jié)構(gòu)處于能量最小狀態(tài)。 然后將溴離子和氧原子的坐標(biāo)固定，使用Nose-Hoover方法控制體系溫度分別在220 K、240 K、260 K、280 K和300 K下進行100 ps（1 ps=10-12s）的NVT馳豫，目的是完成結(jié)構(gòu)松弛，減小結(jié)構(gòu)誤差，時間步長為1 fs（1 fs = 10-15s），每1000步輸出一幀結(jié)構(gòu)圖；然后解除溴離子和氧原子的固定，在相同條件下，選擇NVT系綜進行動力學(xué)模擬，模擬時間為400 ps。圖2為分子動力學(xué)計算過程中體系能量變化圖，可以發(fā)現(xiàn)，在模擬的前100 ps內(nèi)，體系的能量穩(wěn)定在一個固定的值附近，說明此時體系已經(jīng)達到了穩(wěn)定狀態(tài)，后300 ps可以分析體系的相關(guān)性質(zhì)。對于每個體系， 分子動力學(xué)模擬過程均進行三次，保證模擬的準(zhǔn)確性。

2 結(jié)果與討論

2.1 構(gòu)象分析

模型的初始結(jié)構(gòu)是理想的晶體結(jié)構(gòu)，籠型結(jié)構(gòu)框架具有清晰的氫鍵網(wǎng)絡(luò)結(jié)構(gòu)和良好的對稱性（圖

1），在分析籠型結(jié)構(gòu)的穩(wěn)定性時，可以分析氫鍵的規(guī)整程度和氧原子的對稱程度來判斷水合物籠型骨架的變形程度，進而直觀的判斷水合物的穩(wěn)定性情況。

從圖3可以看出在220 K時，體系A(chǔ)和B的水合物氫鍵結(jié)構(gòu)變化很小，客體分子仍然位于晶穴的中心位置，水合物籠型完整，體系C和D的水合物籠型結(jié)構(gòu)發(fā)生扭曲，部分客體分子偏離中心位置，此時水合物的穩(wěn)定狀態(tài)相比A和B體系下降，穩(wěn)定性情況為A ＞B ＞C ＞D，由此得出晶穴占有率越高，水合物越穩(wěn)定。 圖4中在相同占有率情況下，隨著溫度的升高，水合物籠型結(jié)構(gòu)逐漸被破壞，水分子和TBAB逐漸偏離原始位置， 且溫度越高， 偏離的位移越大。240 K時CH4分子開始從籠中逸出，說明水合物開始分解；260 K時體系的籠型結(jié)構(gòu)開始坍塌，水合物兩側(cè)的CH4分子開始逸出， 中心位置的晶穴仍能保持籠型結(jié)構(gòu)，表明水合物的分解是由外側(cè)向內(nèi)側(cè)逐漸演變的；到300 K時，體系的水合物籠型完全坍塌，CH4分子發(fā)生聚集，水合物完全分解，并出現(xiàn)分相現(xiàn)象。 分析圖4可知，在相同占有率的情況下，隨著溫度升高，水合物穩(wěn)定性下降。

圖1初始構(gòu)象中水合物完整的籠型個數(shù)為48。分析不同溫度下ABCD四個體系的最終構(gòu)象圖，可以得出表1， 表中的數(shù)據(jù)為各構(gòu)象圖中完整的水合物籠的個數(shù)。 縱向?qū)Ρ缺碇械臄?shù)據(jù)，可以發(fā)現(xiàn)溫度一定時，晶穴占有率越高，完整的籠型個數(shù)越多，水合物越穩(wěn)定；同樣的，當(dāng)晶穴占有率相同時，隨著溫度的升高，完整籠型的個數(shù)減少，水合物穩(wěn)定性下降。 綜合表1的分析也會得出，晶穴占有率越高，溫度越低，水合物越穩(wěn)定。

表1 各體系在不同溫度下水合物穩(wěn)定性情況（籠型個數(shù)）

2.2 徑向分布函數(shù)

徑向分布函數(shù)（Radial Distribution Function，RDF）是分析水合物結(jié)構(gòu)特征最常用的評價參數(shù)，是晶體結(jié)構(gòu)無序化的一個衡量指標(biāo)，gαβ(r)表示在距離α粒子質(zhì)心r處出現(xiàn)β粒子的概率， 相當(dāng)于體系的區(qū)域密度與平均密度的比值。當(dāng)r值較小時，g(r)若出現(xiàn)幾個極大值， 則在這些r值的附近， 出現(xiàn)β粒子的幾率較大，此時結(jié)構(gòu)處于有序狀態(tài)；當(dāng)r值較大時，g(r)值接近1，說明此時是無規(guī)則狀態(tài)。 若水合物處于穩(wěn)定狀態(tài)， 則水分子中氧原子的排列具有周期性，RDF曲線有多個峰值，峰高隨著r的增大逐漸降低；若水合物處于失穩(wěn)狀態(tài)，第二特征峰及其以后的特征峰將降低或消失。 對于采用周期性邊界條件的體系，RDF的定義如公式（1）所示：式中：V為模擬體系的體積，10-30m3；Nα和Nβ分別為α粒子和β粒子的個數(shù)；niβ(r)表示以第i個α粒子為中心，在r～r+Δr（單位為10-10m）范圍內(nèi)β粒子的數(shù)目。

圖5為不同體系在不同溫度下的水分子中氧原子的徑向分布函數(shù)（RDF）。從圖中可以看出，RDF曲線在r= 2.73 × 10-10m處出現(xiàn)第一特征峰，表明最相近的水分子之間氧原子的距離為2.73 × 10-10m，第二特征峰出現(xiàn)在r= 4.52 × 10-10m，這與文獻[16]一致，驗證了模型的正確性。 隨著溫度的升高，峰值降低，峰谷升高，說明水分子無序性增大，籠型框架發(fā)生扭曲，水合物穩(wěn)定性下降。在A和B體系中，溫度升高到280 K時，第二特征峰消失，之后r值在1上下波動，與液態(tài)水的特征峰相似[17]，得出水合物開始分解。

對比在相同溫度下不同體系的徑向分布函數(shù)可以發(fā)現(xiàn)，CH4含量越大，RDF的特征峰越明顯，峰值相對較高，峰谷較低，表明水分子有序性較強，水合物穩(wěn)定性較好。 A和B體系差別不明顯，說明晶穴占有率在75%以上，水合物處于相對穩(wěn)定狀態(tài)，此時溫度對水合物的穩(wěn)定性占主導(dǎo)因素。 通過不同角度對比可以表明，溫度越低，晶穴占有率越高，水分子的有序性越強，水合物越穩(wěn)定。

2.3 均方位移

在分子動力學(xué)運動過程中，體系中的原子時刻在自由移動，每個瞬間各原子都會偏離原始位置到達某一位置。 原子位移平方的平均值就稱為均方位移（Mean Squared Displacement，MSD），如公式（2）所示。

式中：N為體系中的粒子總數(shù)；Ri(t)、Ri(0)分別表示i粒子在0、t時刻的位置。

均方位移反映了t時刻體系中原子的空間位置與初始位置的偏離程度。 穩(wěn)定的水合物，其分子處于相對固定的位置，MSD的值會在一個較低的值附近波動。 若水合物發(fā)生分解，其分子位置坐標(biāo)相對自由，MSD值會隨著時間的增加而逐漸遞增。

圖6為各體系在不同溫度下水分子中氧原子的均方位移圖，A和B體系中， 溫度在260 K以下時，MSD值很小，表明水分子的位置相對固定，水合物籠型結(jié)構(gòu)較完整，水合物處于穩(wěn)定狀態(tài)。 隨著溫度進一步升高，MSD迅速增大，此時水合物開始分解，水分子處于無序狀態(tài)，偏離初始位置較遠，水合物籠型結(jié)構(gòu)坍塌，且溫度越高斜率越大，由此得出溫度越高，水合物分解速率越快。

相對于A、B體系，C、D體系中只有220 K時MSD較小，說明低溫有利于水合物穩(wěn)定存在，且晶穴占有率越低，MSD越大，籠型結(jié)構(gòu)變形程度較大，水合物極易分解，原因是晶穴占有率越低，水合物中的空籠越多，由于缺少客體分子的支撐，空籠的穩(wěn)定性會下降，水合物的穩(wěn)定性降低。 對比四個圖發(fā)現(xiàn)溫度在260 K以下時，C、D體系水合物穩(wěn)定性遠遠低于A、B體系， 此時晶穴占有率對水合物穩(wěn)定性的影響大于溫度對水合物穩(wěn)定性的影響。

2.4 氫鍵的數(shù)目與鍵長

氫鍵是一些與電負性較大的原子成鍵的氫原子和附近電負性較大或者帶孤對電子的原子形成的較強的非鍵作用，其作用強度介于成鍵作用和非鍵作用之間。 體系中的氫鍵包括水分子之間的氫鍵和溴離子與水分子之間形成的氫鍵。 當(dāng)水合物處于穩(wěn)定狀態(tài)時，氫鍵數(shù)目與長度相對于初始結(jié)構(gòu)的氫鍵變化較小，若水合物發(fā)生分解，則氫鍵的數(shù)目與長度變化較大。本文所選的COMPASSⅡ力場對氫鍵的統(tǒng)計是基于對范德華力和靜電相互作用的綜合分析而實現(xiàn)的。

與RDF和MSD不同， 氫鍵的變化是描述水合物籠型結(jié)構(gòu)隨溫度和晶穴占有率變化而變化的直觀體現(xiàn)。 水合物初始結(jié)構(gòu)中，氫鍵的數(shù)目為1215，平均鍵長為1.82 × 10-10m。在分子動力學(xué)過程中，隨著水合物的分解，部分氫鍵斷裂，部分氫鍵被拉長，導(dǎo)致水合物籠型扭曲。 圖7為分子動力學(xué)之后水合物中氫鍵的數(shù)目和平均鍵長的變化， 從圖中可以看出，隨著溫度升高，氫鍵數(shù)目不斷減少，鍵長逐漸變長，表明水合物籠出現(xiàn)不同程度的扭曲變形；在同一溫度下，晶穴占有率越低，氫鍵數(shù)目越少，鍵長越長，水合物越不穩(wěn)定。 通過分析氫鍵的變化也可直觀地看出低溫和高占有率有利于水合物的穩(wěn)定存在。

3 結(jié)論

（1）溫度較低時，由于受到氫鍵的約束，水分子的位移較小，水合物的穩(wěn)定性較好。 260 K以下時，水合物處于相對穩(wěn)定狀態(tài)，280 K時水合物開始分解，到300 K時水合物已經(jīng)完全分解。

（2）晶穴中存在客體分子時，水合物籠受到客體分子的支撐，不容易坍塌，提高水合物的穩(wěn)定性。晶穴占有率在75%以上時，水合物的穩(wěn)定性較好，晶穴占有率低于50%時，水合物的穩(wěn)定性會明顯下降。

（3）溫度在260 K以下時，晶穴占有率對水合物穩(wěn)定性的影響大于溫度對水合物穩(wěn)定性的影響；晶穴占有率在75%以上時， 溫度對水合物的穩(wěn)定性影響占主導(dǎo)因素。

（4）水合物的分解過程是由外層向內(nèi)層逐層遞進的趨勢，水合物中部籠型的穩(wěn)定性最好。 水合物完全分解后會出現(xiàn)分相現(xiàn)象，甲烷分子會從水籠中逸出并發(fā)生聚集。