四氫呋喃－甲烷水合物穩定性的分子動力學模擬

張慶東，李玉星，王武昌

（中國石油大學（華東）儲運與建筑工程學院，山東 青島266580）

天然氣水合物是在高壓、低溫條件下烴類氣體與水形成非化學計量的、外形似冰的包絡型化合物［1］，其中主體分子（水分子）之間通過氫鍵相互作用構成晶穴，客體分子（烴類分子）填充于晶穴中，主客體分子間的作用力為范德華力。自然界中的天然氣水合物一般以Ⅰ型［2（512）6（51262）］和Ⅱ型［16（512）8（51264）］為主［2］，且晶穴占有率不高（平均只有大約1／3的晶穴被氣體分子占據）［3］。Saito等［4］的研究結果表明甲烷水合物的儲氣量在純水體系中均達不到40%（晶穴占有率＜22.2%）。Stern等［5］通過實驗發現甲烷晶穴占有率達不到100%（最高為94.3%）。分子動力學（MD）模擬是一種處理分子體系的動態研究方法，可以模擬系統隨時間演化的微觀過程。它依據Netwon經典力學求解核骨架的行為，在計算每個瞬時的勢能上依據量子力學，在求算各種體系性質上依據統計力學。分子動力學模擬可以用來研究目前實驗手段難以考察的物理現象和過程，在分子水平上深入分析其反應機理。目前水合物的研究已廣泛采用MD模擬，主要有水合物穩定性機理的MD模擬研究［6－7］、天然氣水合物儲氫穩定性的MD模擬［8］、天然氣水合物傳熱系數的MD模擬計算［9］、動力學抑制劑對水合物的作用的微觀機理模擬［10］。Rodger等［11］應用分子動力學模擬深入研究了客體甲烷分子對Ⅰ型甲烷水合物穩定性的影響，結果表明沒有甲烷分子填充的水合物晶格是不穩定的。丁麗穎等［12］模擬了不同晶穴占有率和不同溫度時Ⅰ型甲烷水合物晶體穩定性。萬麗華等［13－14］研究了不同晶穴占有率對Ⅰ型甲烷水合物導熱性能和結構穩定性的影響。

四氫呋喃（THF）不僅可以和水在常壓、0℃條件下生成Ⅱ型水合物［15］，而且可以降低天然水合物的相平衡壓力［16］。為此，本文構造了8個Ⅱ型水合物晶體模型，分別采用四氫呋喃（THF）和甲烷組合作為客體分子填充晶胞孔穴，應用MD模擬方法，系統地研究客體分子種類、溫度、不同晶穴占有率（θ＝0，33.3%，66.7%，100%）對Ⅱ型水合物穩定性的影響。

1 模型搭建和模擬方法

構建的水合物模型為2×2×2個Ⅱ型水合物晶胞，每個單胞中包括136個水分子，由8個大孔穴（51264）和16個小孔穴（512）組成。THF分子尺寸比較大，只能占據大孔穴，CH4分子可以同時占據大孔穴和小孔穴。本模擬將THF和CH4按不同組合填充于大孔穴和小孔穴中，以研究THF和CH4對THF－CH4水合物結構穩定性的影響。

1.1 模型搭建

Ⅱ型水合物晶胞是面心立方結構，空間群為Fd3m，晶格參數a＝17.175（3）×10－10m［17］。氧原子初始坐標來源于X射線單晶衍射實驗［18］，氫原子隨意放入晶體中，并且滿足Bernal－Fowler規則［19］。

采用簡單點電荷（simple point charge，SPC）模型［20］控制水的相互作用，控制水分子O—H鍵長 為0.9570×10－10m，H—O—H鍵 角 為104.52°。甲烷分子采用剛性模型，C—H鍵長為1.0895×10－10m。THF分子手動搭建，并經過幾何與能量最優化處理，得到總能量最小、結構最穩定構象。

表1 客體分子在模型中的分配

1.2 模擬方法

本模擬采用較精確的全原子模擬，在正則系綜（NVT）條件下，水分子之間選用一致性價力場（consistent valence force field，CVFF），結構優化采用最速下降法（steepest descent）和共軛梯度法（conjugate gradient）。

采用Lennard－Jones勢能模型描述甲烷分子之間、THF分子之間以及它們與水分子之間的非鍵結相互作用，采用SPC力場［20］描述水分子之間的相互作用。SPC力場將水分子視為剛性分子，其鍵長與鍵角的值固定。此力場將二水分子的作用分為范德瓦爾斯與庫倫作用。僅O—O原子間有范德瓦爾斯作用。水分子的每個原子均帶有電荷，分子間原子對有庫倫靜電作用。SPC力場的形式見式（1）。

式中，e為基本電荷量；rOO、rOH和rHH分別為不同水分子氧原子之間、氧原子和氫原子之間和氫原 子 之 間 的 距 離；UOO（rOO）、UOH（rOH）和UHH（rHH）則為相應的相互作用勢；qO和qH分別為氧原子和氫原子的電荷量；ε和δ為Lennard－Jones勢能參數。本文對各分子間的Lennard－Jones作用項采用截斷球近似，截斷半徑Rcutoff＝15×10－10m，并采用Ewald求和法處理長程靜電相互作用，Ewald精度為4.182×10－6kJ／mol，采用Lorentz－Berthelot混合規則處理交互作用參數。

進行分子動力學模擬時，各分子起始速度由Maxwell－Boltzmann分布隨機產生，采用周期性邊界條件和時間平均等效于系綜平均等假設基礎上，運用velocity Verlet［22］算法求解牛頓運動方程。首先對體系內的分子進行最優化得到最優結構，然后對體系進行松弛處理，固定水分子中氧原子位置，使用Nose－Hoover熱容方法控制體系溫度T在263.15K、273.15K的條件下進行等溫等容模擬，模擬時間步長為1fs，每5000步輸出一幀，模擬時間長度為50ps。然后解除對氧原子的約束，在相同條件下對模擬體系進行等溫等容模擬，每200000步輸出一幀，模擬時間長度為4000ps，其中1000ps用于使體系達到平衡。

2 模擬結果及分析

經上述模擬過程得到晶體最終構型，比較在不同溫度下不同水合物籠型結構的變化，初步判斷晶體穩定性。通過計算水合物中水分子中氧原子的徑向分布函數、均方位移以及主客體分子間相互作用能，進一步分析溫度、晶穴占有率和客體分子對晶體穩定性的影響。

2.1 構象分析

圖1給出了模擬前各水合物模型的初始構型，對比8組模型結構可看出，各模型除晶胞內客體分子因方案設置不同外，其他結構均一致，說明各構型初始時籠子骨架主體水分子之間相互作用相同，排除了主體水分子對模擬結果的影響。

圖2、圖3分別給出了模擬溫度為273.15K、263.15K時模擬結束后各模型的最終構型。據圖可看出，在溫度為273.15K時，只有模型B的籠型結構完好并且形變很小，說明在晶穴占有率為100%、THF占據所有大孔穴、CH4占據所有小孔穴時最有利于晶體穩定，水合物可以穩定存在；其余模型均出現不同程度的晶格破裂，表明籠型結構坍塌，晶胞瓦解，水合物不能穩定存在。當溫度降低為263.15K時，各模型結構變形均減小，說明溫度越低，晶胞的穩定性越好，這與水合物在低溫下更穩定的性質相吻合；除模型B外，模型A亦能夠保持籠型結構完好且變形較小，表明此溫度下模型A可以穩定存在。

圖4為模擬溫度為283.15K時模型B的最終構型，對比溫度在263.15K和273.15K時的最終構型可以發現，模型B依然保持完好的籠型結構和較小的形變，表明在客體種類適宜和水合物晶穴占有率較高的情況下，水合物晶體在較高溫度下仍可保持穩定，為水合物在較高溫度下儲存運輸提供了指導。

圖1 模型初始構型

圖2 273.15K時各模型最終構型

圖3 263.15K時各模型最終構型

2.2 徑向分布函數（radial distribution function，RDF）

徑向分布函數gα－β（r）表現的是距參考粒子α質心r處出現粒子β的概率，可以解釋為系統的區域密度（local density）與平均密度（bulk density）的比，反映微觀結構特征。參考分子的附近（r值小）區域密度不同于系統的平均密度，但當距參考分子遠時區域密度應該與平均密度相同，即當r值大時RDF的值應接近1。晶態、非晶態和液態都有自己的特征徑向分布函數。液態水中每個氧原子與周圍最鄰近氧原子之間的距離為r＝2.78×10－10m，氧原子的RDF曲線在此處形成一個峰；當r＞2.78×10－10m時，由于水中分子是均勻分布的，RDF曲線值應該逐漸趨于1。穩定的水合物晶體中水分子的排列具有周期性，氧原子RDF曲線有多個峰，各個峰高隨分子間距離的增加逐漸降低，且按晶體中水分子間隔距離規律分布。圖5和圖6為溫度為273.15K和263.15K時模擬結束后模型內水分子中氧原子的徑向分布函數曲線，其中M和W分別為穩定水合物晶體和液態水的氧原子徑向分布函數曲線。

圖4 283.15K時模型B的最終構型

圖5 273.15K時水合物中氧原子徑向分布函數曲線

在溫度為273.15K和263.15K時，對比兩圖可發現：模型B與模型M的氧原子徑向分布函數函數曲線幾乎重合，表明在模擬結束后模型B的籠型結構骨架氧原子的位置并未發生變化，水晶體未發生坍塌或者變形，模型B可以在此模擬溫度下保持穩定；模型C、D、E、F、G和H第一個峰高小于M，第二、第三等峰高消失，且曲線與液態水W趨于一致，表明模擬結束后，模型中水分子的有序排列減小，氧原子分布無序程度增大，可判斷水合物籠型結構被破壞，水合物晶體逐漸分解；在溫度為273.15K時，模型A的第二、第三等峰消失，水合物晶體不能穩定存在，當溫度降低為263.15K時，模型A、M的氧原子徑向分布函數函數曲線幾乎重合，表明在溫度降低時，晶穴占有率為100%的模型最容易從不穩定狀態轉變為穩定狀態。

圖6 263.15K時氧原子徑向分布函數曲線

2.3 均方位移

分子動力學系統中的原子由起始位置不停移動，每一瞬間各個原子的位置皆不相同。粒子位移平方的平均值稱為均方位移（mean squared displacement，MSD），見式（2）。

式中，括號表示平均值，Ri（t）和Ri（0）分別表示分子i在時間為0、t時的位置，N為體系中的分子總數。根據統計原理，只要分子數目夠多，計算時間夠長，系統的任一瞬間均可當作時間的零點，所計算的平均值應該相同。

均方位移反映了系統中分子的空間位置與其初始位置的偏離程度，根據其大小可以判斷體系是處于固態還是液態。對于穩定的水合物晶體，所有水分子均處于相對固定的晶格點上，氧原子的均方位移應該在一個略大于0的值附近波動，對于液態水，其均方位移隨模擬時間的延長而增大。圖7和圖8給出了模擬溫度分別為273.15K和263.15K時液態水和各個模型水分子中氧原子的均方位移。

圖7 273.15K時氧原子均方位移曲線

圖8 263.15K時氧原子均方位移曲線

從圖7、圖8可以看出，當模擬溫度相同時，各模型氧原子均方位移大小關系相同，均為B＜A＜H＜E＜C＜F＜G＜D，表明客體分子的存在可以提高晶體的穩定性，晶穴占有率越高，晶體越穩定。進一步分析可發現以下幾點。①晶穴占有率為0時，模型D的均方位移均大于液態水的氧原子均方位移，表明氧原子移動劇烈程度大于液態水，水合物籠型結構迅速坍塌，說明晶穴占有率為0時的水合物晶體極易分解，水合物不能穩定存在。②晶穴占有率為33.3%時，均方位移C＜F＜G，對比3個模型的結構可發現：客體分子為THF的晶體比客體分子為CH4的晶體更容易保持穩定，客體分子種類對晶體穩定性有影響，CH4占據小孔穴的晶體比CH4占據大孔穴的晶體更容易保持晶體穩定，客體分子占據的孔穴尺寸對晶體穩定性有影響。③晶穴占有率為66.7%時，模型H的均方位移始終小于模型E，但當溫度降低時，兩者之間的差距減小，表明在此溫度范圍內，客體分子為CH4的晶體相對于客體分子為THF的晶體的穩定性更易受溫度的影響。④晶穴占有率為100%、溫度為273.15K時，模型B可以穩定存在而模型A難以穩定存在，溫度降為263.15K后兩模型均可穩定存在，由于模型A、B不同之處僅在于占據大孔穴的客體分子不同，因此可以判斷，溫度是通過影響模型A中大孔穴與CH4分子之間的相互作用進而影響模型A的穩定性的。

對比兩圖還可發現，模擬溫度由273.15K降低為263.15K時，各模型的氧原子均方位移均減小，表明溫度的降低抑制了水分子的運動，增強了籠型結構的穩固，有利于晶體結構的穩定性；模型B的氧原子均方位移變化很小，說明對于已經可以穩定存在的晶體，溫度的降低對其穩定性的影響不大。

2.4 相互作用能

圖9 相互作用能曲線

從圖9可以看出，相互作用能從小到大依次為B、A、H、C、E、G、F，進一步分析可以發現：模型B的相互作用能最小，說明模型B最穩定，這與前面分析所得結論是一致的；模型C、H的相互作用能小于模型E、F、G，說明THF比CH4更容易保持水合物晶體穩定性；模型E的相互作用能小于模型F、G及模型H的相互作用能小于模型C，說明對于CH4，客體分子晶穴占有率越高，相互作用能越小；模型A的相互作用能比模型C、E、F、G、H小，說明100%的晶穴占有率對相互作用能的降低有較大促進作用，該作用大于不同客體分子種類對相互作用能的影響。

對比不同溫度下的Ebinding可發現，當溫度降低時，各個模型的Ebinding均減小，說明溫度的降低可以減小晶體內主客體分子間的相互作用能，提高晶體的穩定性的：模型A的Ebinding降低最大，說明溫度的降低對其影響最大；模型B的Ebinding受溫度影響較小，進一步證實了當水合物處于穩定狀態時，降低溫度對其穩定性的影響并不大。

3 結 論

采用NVT系綜分子動力學模擬方法研究了不同客體、晶穴占有率以及溫度對Ⅱ型水合物晶體結構穩定性的影響，通過對最終構象、RDF、MSD和相互作用能的分析，得到以下結論。

（1）尺寸較大的THF與大晶穴的相互作用能小于尺寸較小的CH4與大晶穴的相互作用能，更有利于提高水合物晶體的穩定性。在晶穴占有率相同時，客體分子類型對水合物晶體穩定性有較大影響，THF比CH4更容易保持水合物晶體的穩定性。

（2）晶穴中客體分子的存在可以降低主客體分子間的相互作用能，減小水分子的均方位移，提高晶體穩定性。

（3）對于不穩定晶體，溫度的降低可以將水分子的均方位移降低到0附近，減小籠型結構的變形，水合物晶體的分解速率降低，穩定性提高；對于已經穩定的晶體，溫度的降低對已經接近為0的均方位移影響不大，對水合物穩定性的影響不明顯。

由于晶穴占有率為100%的CH4水合物模型的穩定性高于晶穴占有率較低的水合物模型，因此在實際生產中，應該盡可能地提高水合物晶穴占有率以增大水合物的穩定性及儲氣量，添加適量THF占據大晶穴以提高水合物的穩定性，并在保證水合物穩定性的基礎上，盡量提高水合物的儲存溫度以降低儲存成本。

在魚類分類學上，鱖魚屬鱸形目鮨科鱖亞科鱖屬，其中鱖屬有兩個種，一個是翹嘴鱖，該魚生長速度快，個體大，常見為2～2.5kg，最大個體重可達50kg，是人工養殖的主要品種；另一個是大眼鱖，該魚生長緩慢，個體較小，最大個體能長至重2kg。鱖魚最為突出的生物學特性是該魚屬典型的兇猛肉食性魚類，對餌料有較強的分辨能力，終生主要以活魚蝦為食，即使是剛開口的魚苗，也要攝食其他魚類的幼苗。長大后，除食活魚外，還兼食蝦類以及蝌蚪等。

符號說明

A——晶格參數，m

Ebinding——相互作用能，J／mol

Ecage——空胞腔能量，J／mol

Eguest——客體分子自身的能量，J／mol

Etotal——模型總能量，J／mol

e——基本電荷，C

gα－β（r）——徑向分布函數，量綱為1

MSD——氧原子的均方位移，m2

N——體系中的分子總數，個

qH——氫原子的電荷量與基本電荷的比值，量綱為1

qO——氧原子的電荷量與基本電荷的值，量綱為1

Rcutoff——截斷半徑，m

Ri（0）——分子在時間為0時的位置，m

Ri（t）——分子在時間為t時的位置，m

rHH——不同水分子中氫原子間的距離，m

rOH——不同水分子中氧原子與氫原子間的距離，mrOO——氧原子間的距離，m

T——溫度，K

t——時間，s

UHH（rHH）——不同水分子中氫原子間的相互作用勢能，J

UOH（rOH）——不同水分子中氧原子與氫原子間的相互作用勢能，J

UOO（rOO）——不同水分子中氧原子間的相互作用勢能，J

ε——Lennard－Jones勢能參數，J／mol

δ——Lennard－Jones勢能參數，m

［1］ Sloan E D.Clathrate Hydrate of Natural Gase［M］.2th ed.New York：Marvel Dekker Inc.，1990：27－64.

［2］ Subramanian S，Ballard A I，Kini R A.Structural transitions in methane＋ethane gas hydrate－PartⅠ：Upper transition point and applications［J］.ChemicalEngineering Science，2000，55（23）：5763－5771.

［3］ 裘俊紅，郭天民.甲烷水合物在含抑制劑體系中的生成動力學［J］.石油學報：石油化工，1998，14（1）：1－5.

［4］ Saim Y，Kawaski，Okui T.Methane storage in hydrate phase with water soluble guests［C］／／Second International Conference on Natural Gas Hydrate.France：Toulouse，1996：459－465.

［5］ Stern L A，Kirby S H，Durham W B.Peculiarities of methane clathrate hydrate formation and solid state deformation，including possible superheating of water ice［J］.Science，1996，273（27）：1843－1848.

［6］ Cheng Wei，Wu Hucai，Ye Xiaoqin，et al.Molecular dynamics study on the structureⅠhelium hydrate［J］.ProgressinNaturalScience，2004，14（11）：1015

［7］ Kenichior K，Hideki T，Koichiro N.On the stability of clathrate hydrates encaging polar guest molecules：Contrast in the hydrogen bonds of methylamine and methanol hydrates［J］.MolecularSimulation，1994，12（3－6）：241.

［8］ Strobel T A，Koh C A，Sloan E D.Water cavities of sh clathrate hydrate stabilized by molecular hydrogen［J］.JournalofPhysicalChemistry，2008，112（7）：1885－1887.

［9］ Inoue R，Tanaka H，Nakanishi K.Molecular dynamics simulation study of the anomalous thermal conductivity of clathrate hydrates［J］.JournalofPhysicalChemistry，1996，104：9569.

［10］ Bjorn K，Tatyana K，Kjetil A.Molecular dynamics simulations for selection of kinetic hydrate inhibitors［J］.JournalofMolecularGraphics＆Modelling，2005，23（6）：524－536.

［11］ Moon C，Taylor P C，Rodger P M.Molecular dynamics study of gas hydrate formation［J］.Journalofthe AmericanChemicalSociety，2003，81：451.

［12］ 丁麗穎，耿春宇，趙月紅，等.Ⅰ型甲烷水合物晶體穩定性的分子動力學模擬［J］.計算機與應用化學，2007，24（5）：570.

［13］ 萬麗華，梁德青.客體分子數對甲烷水合物導熱性能影響的分子動力學模擬［J］.化工學報，2012，63（2）：382－386.

［14］ 萬麗華，顏克鳳，魯濤，等.分子動力學模擬甲烷水合物結構穩定性［J］.武漢理工大學學報，2010，32（2）：25－27.

［15］ 馬應海，茍蘭濤，何曉霞，等.四氫呋喃水合物零度以上生成動力學研究［J］.天然氣地球科學，2006（2）：244－248.

［16］ 涂運中，蔣國盛，張凌，等.SDS和THF對甲烷水合物合成影響的實驗研究［J］.現代地質，2008（3）：485－488.

［17］ Kirchner M T，Boese R，Billups W E，et al.Gas hydrate single－crystal strucure analyses［J］.Journalofthe AmericanChemicalSociety，2004，26：9407－9412.

［18］ Hlooander F，Jeffrey G A.Neutron diffraction study of the crystal structure of ethylene oxide deuteron hydrate at 80K［J］.JournaloftheAmericanChemicalSociety，1977，66：4699－4705.

［19］ Bernal J D，Fowler R H.A theory of water and ionic solution，with particular preference to hydrogen and hydroxylions［J］.TheJournalofChemicalPhysics，1933，1：515－548.

［20］ Berendsen H J C，Grigera J R，Straastsma T P.The missing term in effective pair potentials［J］.TheJournal ofPhysicalChemistry，1987，91：6269－6271.

［21］ Allen M P.Tildeslay D J.Computer Simulation of Liquids［M］.2th ed.Oxford：Clarendon Press，1987：156－162.