張錫彥 ，梁海峰，邵偉強，張 華

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

天然氣相對于其他化石燃料具有污染排放量少和價格較低的優(yōu)勢[1]，因此如何大規(guī)模且安全有效地運輸天然氣成為亟待解決的問題[2]。天然氣水合物儲運方法具有安全經(jīng)濟、環(huán)境友好和儲氣容量高等優(yōu)勢[3]。I型、Ⅱ型和H型水合物的最大潛能分別為56體積、154體積和201體積[1,4-5]。比較發(fā)現(xiàn)，H型水合物在天然氣儲運方面比I型和Ⅱ型更具發(fā)展前景。

天然氣水合物的形成需要在高壓和低溫等特定的熱力學(xué)條件下進(jìn)行，使得氣體水合物物理和結(jié)構(gòu)特性的測量不準(zhǔn)確且成本高，給經(jīng)濟、安全和環(huán)境方面帶來各種挑戰(zhàn)[6]。分子動力學(xué)模擬可以在納米尺度上模擬天然氣水合物體系的成核、形成和分解，因此，近年來，該方法被廣泛應(yīng)用于天然氣水合物結(jié)構(gòu)以及穩(wěn)定性的研究[7]。丁麗穎等[8]模擬了不同晶穴占有率和溫度對I型甲烷（CH4）水合物晶體穩(wěn)定性的影響。萬麗華等[9-10]研究了客體分子數(shù)以及晶穴占有率對I型CH4水合物導(dǎo)熱性能和晶體結(jié)構(gòu)穩(wěn)定性的影響。張慶東等[11]研究了以四氫呋喃（THF）和CH4作為客體分子時，溫度和晶穴占有率對Ⅱ型水合物穩(wěn)定性的影響，得出了添加適量的THF占據(jù)大晶穴可以提高水合物穩(wěn)定性的結(jié)論。

H型水合物是RIPMEESTER等[12]于1987年發(fā)現(xiàn)的一種與I型和Ⅱ型結(jié)構(gòu)不同的新型氣體水合物。梅東海等[13]首次在分子水平上研究了H型水合物結(jié)構(gòu)及其穩(wěn)定性，提出了H型氣體水合物晶體的大、中、小3種晶穴均需被客體分子占據(jù)時才能保持穩(wěn)定。王璐琨等[14]通過分子動力學(xué)模擬為H型氣體水合物導(dǎo)熱系數(shù)的模擬計算提供了可靠的勢能模型和模擬算法。施軍錁等[15]考察了溫度對H型水合物晶體結(jié)構(gòu)穩(wěn)定性的影響，提出了H型水合物穩(wěn)定性隨溫度升高而降低并且籠形結(jié)構(gòu)有分解的趨勢。周生平等[16]采用分子動力學(xué)方法模擬了晶穴占有率和溫度變化對H型氣體水合物穩(wěn)定性的影響，表明晶穴填充物對H型水合物穩(wěn)定性的影響作用高于溫度。HAMID等[17]對CH4+大客體分子（新己烷、甲基環(huán)己烷（C6H11-CH3）和金剛烷等）二元H型水合物進(jìn)行了正則系綜（NVT）分子動力學(xué)模擬，研究了客體分子對H型水合物穩(wěn)定性的影響，得出CH4+大客體分子的存在有助于提高H型水合物的穩(wěn)定性。武文志等[1]在恒容條件下，通過實驗分析了不同過冷度對環(huán)庚烷（C7H14）、環(huán)庚酮（C7H12O）、環(huán)辛烷和甲基環(huán)己烷 4種大客體分子和CH4體系的H型水合物生成過程的影響，發(fā)現(xiàn)過冷度越大，水合物生成誘導(dǎo)時間越短，但是以上兩篇文獻(xiàn)均沒有對比不同大客體分子對H型CH4水合物的影響。中子衍射實驗從結(jié)構(gòu)角度表明了H型水合物中大晶穴理論上最多容納的CH4分子數(shù)為5個[18]。KUMAZAKI等[19]提出，在溫度t= 50 °C，壓力p= 1 GPa時，H型水合物大晶穴可容納2~3個CH4分子并且晶穴結(jié)構(gòu)穩(wěn)定，但此文獻(xiàn)沒有提出CH4分子數(shù)目變化對H型水合物穩(wěn)定性的影響。孫志高等[20]利用恒溫壓力搜索法測量了CH4和C6H11-CH3體系水合物形成的相平衡條件，提出了H型水合物相平衡壓力比I型低1.0 MPa以上。

在石油中常含有一些能形成H型水合物的烴類物質(zhì)，近年來學(xué)者已測定了許多大分子烴類水合物的形成條件[20]，但是對于大分子烴類水合物的研究多集中在單一促進(jìn)效果的探討，目前還鮮有不同大客體分子對H型甲烷水合物影響的報道。本文基于Materials Studio 8.0分子動力學(xué)平臺考察C7H12O、C6H11-CH3和C7H143種大客體分子對H型氣體水合物穩(wěn)定性的影響，同時研究在大籠中大客體分子數(shù)不變時，大籠中CH4分子數(shù)對水合物穩(wěn)定性的影響。

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

模擬系統(tǒng)由2 × 2 × 2的H型水合物單晶胞構(gòu)成，模型包含272個水（H2O）分子，24個512晶穴，16個435663晶穴和8個51268晶穴，將512和435663晶穴加起來研究，即40個小籠、8個大籠。本文研究的C7H12O、C6H11-CH3和C7H143種油狀大分子有機物只能占據(jù)51268晶穴，CH4由于分子結(jié)構(gòu)比較小既可以占據(jù)大晶穴也可以占據(jù)小晶穴。C7H12O、C6H11-CH3和C7H143種油狀大分子有機物以及H型水合物空籠的結(jié)構(gòu)如圖1所示。

圖1 大分子有機物和H型水合物空籠結(jié)構(gòu)示意Fig. 1 Macromolecular organics and empty cage structure of structure H hydrate

1.1 模型構(gòu)建

選用Materials Studio 8.0分子動力學(xué)模擬軟件進(jìn)行模型的構(gòu)建和求解，晶胞中氧原子的初始坐標(biāo)根據(jù)H型水合物晶體衍射試驗[12,21-24]確定。其中，水分子構(gòu)成水合物籠子，氫原子隨意放入晶體中并且呈無序排列，滿足Bernal-Fowler規(guī)則[25]。在x、y、z坐標(biāo)方向上采用周期性邊界條件。

首先采用蒙特卡洛吸附[26]在初始晶胞中的大晶穴分別填充8個C7H14、C6H11-CH3和C7H12O，每個小晶穴填充40個CH4，并將這3種模型命名為A、B、C。其次，占用6個大晶穴用來填充C7H14，小晶穴共填充40個CH4，剩余的兩個大晶穴分別放入0、1、2、3和4個CH4分子，并將上述5種模型分別命名為H-0、H-1、H-2、H-3和H-4（例如：H-1即為8個大晶穴中占用其中6個用來填充大客體分子C7H14，而在剩余的兩個大晶穴中，分別放入1個CH4分子）?？腕w分子在模型中的分配如表1所示。

表1 客體分子在模型中的分配Table 1 Distribution of guest molecules in model

1.2 模擬方法

本模擬在恒溫恒壓系綜（NPT）條件下，研究C7H12O、C6H11-CH3和C7H143種大分子有機物對H型水合物穩(wěn)定性的影響。模擬采用一致性價力場（Consistent valence force field, CVFF）來描述研究體系[27]；采用最速下降法和共軛梯度法[26]進(jìn)行結(jié)構(gòu)優(yōu)化；采用Ewald加和法描述長程靜電作用力[28]和范德瓦爾斯作用力，加和精度為4.18 × 10-4kJ/mol；采用Noose-Hoover熱浴法控制模型溫度；采用Berendsen方法控制壓力；采用簡單點電荷（Simple point charge，SPC）勢能模型[29]控制水分子之間的相互作用；采用Lennard-Jones勢能函數(shù)來描述水分子和其他分子間非鍵相互作用，截斷半徑為1.2 nm。

本模擬對表1所示的8種模型進(jìn)行松弛，將模型中的氧原子固定位置，模型溫度分別在290 K、300 K、310 K和320 K的條件下進(jìn)行NVT馳豫，時間步長為1 fs（1 fs = 10-15s），模擬時間為100 ps（1 ps = 10-12s），每5000步輸出一幀結(jié)構(gòu)圖，然后解除氧原子固定，在相同條件下，采用NPT系綜進(jìn)行動力學(xué)模擬，控制壓力為3 MPa、8 MPa[20]，模擬時間為400 ps，輸出步長為5000步。在以上溫度壓力條件下，研究了環(huán)庚酮、甲基環(huán)己烷和環(huán)庚烷3種油狀大分子有機物以及客體甲烷分子數(shù)對H型水合物穩(wěn)定性的影響。

2 結(jié)果與討論

經(jīng)過對上述8個模型進(jìn)行NVT馳豫以及NPT分子動力學(xué)模擬過程得到晶體最終構(gòu)型，比較不同溫度條件下A、B和C 3種水合物籠型結(jié)構(gòu)的變化，初步判定晶體結(jié)構(gòu)的穩(wěn)定性，通過計算水合物中主體分子氧原子的徑向分布函數(shù)、均方位移以及主客體分子的相對濃度（ Relative concentration，RC，代表特定法線方向距離r附近，A粒子數(shù)密度與系統(tǒng)內(nèi)該粒子總數(shù)密度的比值），來系統(tǒng)地分析溫度、壓力以及不同客體大分子和CH4分子數(shù)對H型水合物晶體穩(wěn)定性的影響。

2.1 最終構(gòu)型

在p= 8 MPa，T= 290 K、300 K、310 K和320 K時，A、B和C 3種模型的最終構(gòu)型分別如圖2、圖3和圖4所示。當(dāng)T= 290 K時，模型A、B和C都具有清晰的籠形結(jié)構(gòu)，氫鍵完整并且氧原子具有良好的對稱性，表明此時3種模型可以維持穩(wěn)定的水合物結(jié)構(gòu)。模型A和B的形變程度隨著溫度的升高而增大，通過對比發(fā)現(xiàn)：當(dāng)T≥ 300 K時，模型A開始分解；當(dāng)T≥ 310 K時，模型B開始分解。模型C在T= 300 K和310 K時依然可以穩(wěn)定存在。當(dāng)T= 320 K時，模型A、B和C都發(fā)生了分解，但模型C形變最小，水合物籠形結(jié)構(gòu)坍塌數(shù)目最少。

圖2 p = 8 MPa時模型A最終構(gòu)型Fig. 2 Final configuration of model A at p = 8 MPa

圖3 p = 8 MPa時模型B最終構(gòu)型Fig. 3 Final configuration of model B at p = 8 MPa

圖4 p = 8 MPa時模型C最終構(gòu)型Fig. 4 Final configuration of model C at p = 8 MPa

通過最終構(gòu)象可以發(fā)現(xiàn)溫度越低，水合物晶胞越穩(wěn)定，并且水合物結(jié)構(gòu)從邊緣向內(nèi)層逐漸分解，其坍塌導(dǎo)致客體分子釋放出來，逃逸出來的客體分子會聚集在一起形成團簇。模型C在較高溫度下仍可保持穩(wěn)定，表明此時含C7H14的H型水合物穩(wěn)定性最好。

2.2 徑向分布函數(shù)

徑向分布函數(shù)（Radial distribution function，RDF）表示在距離α粒子質(zhì)心r（單位：× 10-10m）處出現(xiàn)β粒子的概率，即系統(tǒng)區(qū)域密度與平均密度的比值。徑向分布函數(shù)（gαβ(r)）的表達(dá)式[7]如式(1)所示：

式中，V為系統(tǒng)體積，×10-30m3；Nα為α粒子的數(shù)目；Nβ為β粒子的數(shù)目；niβ(r)為以第i個粒子為中心，在r~(r+ Δr)范圍內(nèi)β粒子的數(shù)目。

對于水合物來說，當(dāng)水合物處于穩(wěn)定狀態(tài)時，由于氧原子處于初始位置附近，晶格排列有序，并且RDF曲線有多個特征峰，隨著r的增大，特征峰逐漸降低；當(dāng)水合物處于分解狀態(tài)或不穩(wěn)定狀態(tài)時，氧原子處于無序狀態(tài)，其RDF曲線的第二特征峰及其之后的特征峰將降低直至消失。

模擬時間400 ps，p= 8 MPa，T= 290 K、300 K、310 K和320 K時，A、B和C 3種模型主體水分子中氧原子的徑向分布函數(shù)如圖5所示。通過不同溫度下RDF曲線可以得出，當(dāng)模擬溫度不同時，模型C中氧原子的RDF曲線有多個特征峰，模型C水分子中氧原子的RDF曲線在r= 2.775 × 10-10m處出現(xiàn)第一特征峰，表明所有相鄰的氧原子都在所獲得的距離上，第二和第三特征峰分別出現(xiàn)在r= 4.575 × 10-10m和r= 6.525 × 10-10m處，這與文獻(xiàn)[17]一致，驗證了模型的準(zhǔn)確性。在T= 320 K時模型C中前3個特征峰依然存在，只是第二和第三特征峰相對比于低溫時峰高有所下降，峰谷有所升高，表明模型C在T= 320 K時才開始分解，并且分解程度不大；模型A在T= 300 K時，相較于T= 290 K，第二和第三特征峰的峰高明顯大幅下降，表明模型A在T= 300 K時水合物籠型部分結(jié)構(gòu)被破壞，水合物晶體逐漸分解；模型B在T= 310 K時，相較于T= 290 K和300 K時第二和第三特征峰的峰高大幅下降，表明模型B在T= 310 K時就開始分解。

圖5 不同溫度下模型A、B和C水分子中氧原子的徑向分布函數(shù)Fig. 5 RDF of oxygen atoms in water molecules of A, B and C models at different temperatures

T= 290 K、p= 3 MPa時H-0、H-1、H-2、H-3和H-4 5種模型水分子中氧原子的徑向分布函數(shù)如圖6所示。由圖6可知，H-3和H-4的RDF曲線有多個明顯的特征峰并且峰高最高，峰谷最低。

圖6 T = 290 K、 p = 3 MPa時模型H-0、H-1、H-2、H-3和H-4水分子中氧原子徑向分布函數(shù)Fig. 6 RDF of oxygen atoms in water molecules of H-0, H-1,H-2, H-3 and H-4 models at T = 290 K and p = 3 MPa

在恒定壓力的同一模型下，隨著溫度的升高，水分子中氧原子的RDF曲線的峰高降低，但是峰谷升高，由于氧原子的有序度小，這種行為表現(xiàn)為水分子和水合物籠中氧原子的低穩(wěn)定性，因此，當(dāng)壓力恒定，溫度升高時，水合物分解的可能性將增加。

2.3 均方位移

在分子動力學(xué)運動過程中，模型中粒子時刻不停移動，使得每個瞬間各個粒子的位置都不一樣。粒子位移平方的平均值稱為均方位移（Mean squared displacement,MSD，單位：× 10-20m2）[30]，其表達(dá)式如式(2)所示：

式中，N為模型中的粒子總數(shù)；Ri(t)為粒子i在時間為t時的位置；Ri(0)為粒子i在時間為0時的位置。

MSD反映了某一時刻模型中分子的空間位置與初始位置的偏離程度。由于分子不能在晶體（固體）中自由移動，所以穩(wěn)定的水合物晶體中氧原子的均方位移值接近一個恒定值，即MSD的斜率接近于0。而當(dāng)水合物處于分解或失穩(wěn)狀態(tài)時，由于分子的隨機運動，MSD的大小隨模擬時間線性增加，說明氧原子偏離初始位置，此時氧原子的位置坐標(biāo)相對自由，水合物處于分解或失穩(wěn)狀態(tài)。

模擬時間400 ps，p= 8 MPa，T= 290 K、300 K、310 K和320 K時，A、B和C 3種模型主體水分子中氧原子的MSD如圖7所示。

圖7 不同溫度下模型A、B和C水分子中氧原子的均方位移Fig. 7 MSD of oxygen atoms in water molecules of A, B and C models at different temperatures

由圖7可知，當(dāng)T不同時，A、B和C 3種模型中氧原子的均方位移都是模型C最小，并且在T≤ 310K時模型C的MSD接近一個恒定值。此時，水合物在初始位置上轉(zhuǎn)動或振動，偏離位置較小。由于分子不能在晶體（固體）中自由移動，可推測在溫度T≤ 310 K時模型C中水合物的籠形結(jié)構(gòu)完整，此時模型C是一個穩(wěn)定的水合物晶體。由于分子的隨機運動，MSD的大小隨模擬時間線性增加，說明氧原子偏離初始位置，此時氧原子的位置坐標(biāo)相對自由，水合物處于分解或失穩(wěn)狀態(tài)。根據(jù)MSD的線性變化可以得出模型A從T= 300 K開始分解，模型B從T= 310 K開始分解，模型C從T= 320 K才開始分解。

T= 290 K、p= 3 MPa時H-0、H-1、H-2、H-3和H-4 5種模型水分子中氧原子的均方位移如圖8所示。由圖8可知，H-3和H-4的MSD值較小并且斜率接近為零。

圖8 T = 290 K、p = 3 MPa時模型H-0、H-1、H-2、H-3和H-4水分子中氧原子的均方位移Fig. 8 MSD of oxygen atoms in water molecules of H-0, H-1,H-2, H-3 and H-4 models at T = 290 K and p = 3 MPa

通過對比不同溫度下同一模型的MSD，發(fā)現(xiàn)隨著溫度升高，MSD值增大，表明低溫有利于水合物晶體穩(wěn)定；通過對比同一溫度下不同油狀大分子有機物的MSD得出以C7H14（模型C）作為大客體分子的H型水合物穩(wěn)定性最好；通過控制溫度和壓力，減小大客體分子數(shù)，在剩余的大籠中分別填充0、1、2、3和4個CH4分子，發(fā)現(xiàn)水合物的穩(wěn)定性隨之升高。

2.4 相對濃度

相對濃度曲線可以用于分析系統(tǒng)內(nèi)各粒子的密度分布信息。

式中，ρr為在r處A粒子的相對濃度；ρi為r處附近區(qū)域的A粒子數(shù)密度，個/ × 10-30m3；ρtotal為A粒子在盒子內(nèi)的總數(shù)密度，個/ × 10-30m3。

通過前文分析得出C7H14為大客體分子時H型水合物穩(wěn)定性最好，條件最溫和，所以進(jìn)一步采用模型C來分析H型水合物晶體在穩(wěn)定和分解時模型內(nèi)主客體分子相對濃度的變化。

T= 290 K、p= 3 MPa（穩(wěn)定時）和T= 320 K、p= 3 MPa（分解時）模型內(nèi)主體水分子和大客體分子C7H14以及小客體分子CH4的相對濃度隨x軸的變化趨勢如圖9所示，其中縱軸代表分子的相對濃度，橫軸代表水合物盒子(1,0,0)方向的位置坐標(biāo)。由圖9(a)可知，水分子在T= 290 K時呈周期性對稱分布，這與初始建立的2 × 2 × 2周期性模型相對應(yīng)。穩(wěn)定時，水分子沿著x軸呈兩個峰高和一個峰谷周期性分布，但當(dāng)水合物分解時，只有r= 15 × 10-10m附近的峰高存在，外層的峰消失，說明水合物籠形結(jié)構(gòu)從外層開始分解。由圖9(b)和9(c)可知，穩(wěn)定時，C7H14和CH4分子沿著x軸對稱分布，有明顯的峰值。當(dāng)T= 320 K時，水合物籠形結(jié)構(gòu)外層開始坍塌，大籠內(nèi)的C7H14和CH4分子逃逸，在水合物模型內(nèi)自由移動，導(dǎo)致相對濃度曲線峰高降低、峰谷升高，相對濃度曲線對稱性與周期性分布消失。

在水合物晶體處于穩(wěn)定狀態(tài)時，模型內(nèi)H2O、C7H14以及CH4分子沿著(1,0,0)方向呈規(guī)則周期性分布。在水合物晶體處于失穩(wěn)狀態(tài)時，晶穴內(nèi)部的C7H14和CH4分子逃逸并在水合物模型內(nèi)自由移動，導(dǎo)致相對濃度曲線對稱性與周期性分布消失。

3 結(jié)論

在NPT系綜下，采用分子動力學(xué)模擬方法研究了C7H12O、C6H11-CH3和C7H143種油狀大分子有機物作為大客體分子對H型水合物穩(wěn)定性的影響，并在相同溫度壓力條件下，保持大客體分子數(shù)不變時，研究了客體CH4分子數(shù)對H型氣體水合物穩(wěn)定性的影響，得到以下主要結(jié)論。

（1）C7H12O、C6H11-CH3和C7H14作為客體分子可以占據(jù)H型CH4水合物大晶穴，并且油狀大分子有機物的存在有助于提高H型CH4水合物的穩(wěn)定性。

（2）H型水合物的穩(wěn)定性隨模擬溫度的升高而降低，不同油狀大分子有機物作用下，水合物穩(wěn)定性順序依次為C7H14> C6H11-CH3> C7H12O。

（3）水合物籠形結(jié)構(gòu)是由邊緣向中心逐層分解的。水合物籠形結(jié)構(gòu)坍塌而釋放出來的客體分子會聚集形成團簇。

（4）當(dāng)大籠中大客體分子數(shù)不變時，水合物的穩(wěn)定性隨著大籠中CH4分子數(shù)的增加而升高，大籠中最多可以穩(wěn)定存在4個CH4分子，且每個大籠中存在4個CH4分子的H型水合物穩(wěn)定性最高。