TBAB半籠型甲烷水合物穩定性分子動力學模擬

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

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

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

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

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

1 模型構建和模擬方法

1.1 模型構建

TBAB·38H2O水合物單晶胞屬于正交晶系，空間群為pmma， 水分子中氧原子和TBAB中原子的初始坐標根據XRD結果確定[14]，水分子和溴離子構成水合物籠，氫原子無序排列，且符合Bernal-Fowler規則[15]。模擬體系由2 × 2 × 2個TBAB半籠型水合物單晶胞構成， 在x、y、z坐標方向上采用周期性邊界條件，晶胞尺寸為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分子填充。 本文通過分子動力學模擬研究不同晶穴占有率（θ=100%、75%、50%、25%）在不同溫度下對TBAB半籠型水合物穩定性的影響，模型初始結構如圖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勢能函數描述水分子和其它分子間非鍵相互作用。

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

2 結果與討論

2.1 構象分析

模型的初始結構是理想的晶體結構，籠型結構框架具有清晰的氫鍵網絡結構和良好的對稱性（圖

1），在分析籠型結構的穩定性時，可以分析氫鍵的規整程度和氧原子的對稱程度來判斷水合物籠型骨架的變形程度，進而直觀的判斷水合物的穩定性情況。

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

圖1初始構象中水合物完整的籠型個數為48。分析不同溫度下ABCD四個體系的最終構象圖，可以得出表1， 表中的數據為各構象圖中完整的水合物籠的個數。 縱向對比表中的數據，可以發現溫度一定時，晶穴占有率越高，完整的籠型個數越多，水合物越穩定；同樣的，當晶穴占有率相同時，隨著溫度的升高，完整籠型的個數減少，水合物穩定性下降。 綜合表1的分析也會得出，晶穴占有率越高，溫度越低，水合物越穩定。

表1 各體系在不同溫度下水合物穩定性情況（籠型個數）

2.2 徑向分布函數

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

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

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

2.3 均方位移

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

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

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

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

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

2.4 氫鍵的數目與鍵長

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

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

3 結論

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

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

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

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