四氫呋喃－甲烷水合物穩(wěn)定性的分子動(dòng)力學(xué)模擬

張慶東，李玉星，王武昌

（中國(guó)石油大學(xué)（華東）儲(chǔ)運(yùn)與建筑工程學(xué)院，山東 青島266580）

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

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

1 模型搭建和模擬方法

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

1.1 模型搭建

Ⅱ型水合物晶胞是面心立方結(jié)構(gòu)，空間群為Fd3m，晶格參數(shù)a＝17.175（3）×10－10m［17］。氧原子初始坐標(biāo)來(lái)源于X射線單晶衍射實(shí)驗(yàn)［18］，氫原子隨意放入晶體中，并且滿足Bernal－Fowler規(guī)則［19］。

采用簡(jiǎn)單點(diǎn)電荷（simple point charge，SPC）模型［20］控制水的相互作用，控制水分子O—H鍵長(zhǎng) 為0.9570×10－10m，H—O—H鍵 角 為104.52°。甲烷分子采用剛性模型，C—H鍵長(zhǎng)為1.0895×10－10m。THF分子手動(dòng)搭建，并經(jīng)過(guò)幾何與能量最優(yōu)化處理，得到總能量最小、結(jié)構(gòu)最穩(wěn)定構(gòu)象。

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

1.2 模擬方法

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

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

式中，e為基本電荷量；rOO、rOH和rHH分別為不同水分子氧原子之間、氧原子和氫原子之間和氫原 子 之 間 的 距 離；UOO（rOO）、UOH（rOH）和UHH（rHH）則為相應(yīng)的相互作用勢(shì)；qO和qH分別為氧原子和氫原子的電荷量；ε和δ為L(zhǎng)ennard－Jones勢(shì)能參數(shù)。本文對(duì)各分子間的Lennard－Jones作用項(xiàng)采用截?cái)嗲蚪?，截?cái)喟霃絉cutoff＝15×10－10m，并采用Ewald求和法處理長(zhǎng)程靜電相互作用，Ewald精度為4.182×10－6kJ／mol，采用Lorentz－Berthelot混合規(guī)則處理交互作用參數(shù)。

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

2 模擬結(jié)果及分析

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

2.1 構(gòu)象分析

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

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

圖4為模擬溫度為283.15K時(shí)模型B的最終構(gòu)型，對(duì)比溫度在263.15K和273.15K時(shí)的最終構(gòu)型可以發(fā)現(xiàn)，模型B依然保持完好的籠型結(jié)構(gòu)和較小的形變，表明在客體種類(lèi)適宜和水合物晶穴占有率較高的情況下，水合物晶體在較高溫度下仍可保持穩(wěn)定，為水合物在較高溫度下儲(chǔ)存運(yùn)輸提供了指導(dǎo)。

圖1 模型初始構(gòu)型

圖2 273.15K時(shí)各模型最終構(gòu)型

圖3 263.15K時(shí)各模型最終構(gòu)型

2.2 徑向分布函數(shù)（radial distribution function，RDF）

徑向分布函數(shù)gα－β（r）表現(xiàn)的是距參考粒子α質(zhì)心r處出現(xiàn)粒子β的概率，可以解釋為系統(tǒng)的區(qū)域密度（local density）與平均密度（bulk density）的比，反映微觀結(jié)構(gòu)特征。參考分子的附近（r值?。﹨^(qū)域密度不同于系統(tǒng)的平均密度，但當(dāng)距參考分子遠(yuǎn)時(shí)區(qū)域密度應(yīng)該與平均密度相同，即當(dāng)r值大時(shí)RDF的值應(yīng)接近1。晶態(tài)、非晶態(tài)和液態(tài)都有自己的特征徑向分布函數(shù)。液態(tài)水中每個(gè)氧原子與周?chē)钹徑踉又g的距離為r＝2.78×10－10m，氧原子的RDF曲線在此處形成一個(gè)峰；當(dāng)r＞2.78×10－10m時(shí)，由于水中分子是均勻分布的，RDF曲線值應(yīng)該逐漸趨于1。穩(wěn)定的水合物晶體中水分子的排列具有周期性，氧原子RDF曲線有多個(gè)峰，各個(gè)峰高隨分子間距離的增加逐漸降低，且按晶體中水分子間隔距離規(guī)律分布。圖5和圖6為溫度為273.15K和263.15K時(shí)模擬結(jié)束后模型內(nèi)水分子中氧原子的徑向分布函數(shù)曲線，其中M和W分別為穩(wěn)定水合物晶體和液態(tài)水的氧原子徑向分布函數(shù)曲線。

圖4 283.15K時(shí)模型B的最終構(gòu)型

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

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

圖6 263.15K時(shí)氧原子徑向分布函數(shù)曲線

2.3 均方位移

分子動(dòng)力學(xué)系統(tǒng)中的原子由起始位置不停移動(dòng)，每一瞬間各個(gè)原子的位置皆不相同。粒子位移平方的平均值稱(chēng)為均方位移（mean squared displacement，MSD），見(jiàn)式（2）。

式中，括號(hào)表示平均值，Ri（t）和Ri（0）分別表示分子i在時(shí)間為0、t時(shí)的位置，N為體系中的分子總數(shù)。根據(jù)統(tǒng)計(jì)原理，只要分子數(shù)目夠多，計(jì)算時(shí)間夠長(zhǎng)，系統(tǒng)的任一瞬間均可當(dāng)作時(shí)間的零點(diǎn)，所計(jì)算的平均值應(yīng)該相同。

均方位移反映了系統(tǒng)中分子的空間位置與其初始位置的偏離程度，根據(jù)其大小可以判斷體系是處于固態(tài)還是液態(tài)。對(duì)于穩(wěn)定的水合物晶體，所有水分子均處于相對(duì)固定的晶格點(diǎn)上，氧原子的均方位移應(yīng)該在一個(gè)略大于0的值附近波動(dòng)，對(duì)于液態(tài)水，其均方位移隨模擬時(shí)間的延長(zhǎng)而增大。圖7和圖8給出了模擬溫度分別為273.15K和263.15K時(shí)液態(tài)水和各個(gè)模型水分子中氧原子的均方位移。

圖7 273.15K時(shí)氧原子均方位移曲線

圖8 263.15K時(shí)氧原子均方位移曲線

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

對(duì)比兩圖還可發(fā)現(xiàn)，模擬溫度由273.15K降低為263.15K時(shí)，各模型的氧原子均方位移均減小，表明溫度的降低抑制了水分子的運(yùn)動(dòng)，增強(qiáng)了籠型結(jié)構(gòu)的穩(wěn)固，有利于晶體結(jié)構(gòu)的穩(wěn)定性；模型B的氧原子均方位移變化很小，說(shuō)明對(duì)于已經(jīng)可以穩(wěn)定存在的晶體，溫度的降低對(duì)其穩(wěn)定性的影響不大。

2.4 相互作用能

圖9 相互作用能曲線

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

對(duì)比不同溫度下的Ebinding可發(fā)現(xiàn)，當(dāng)溫度降低時(shí)，各個(gè)模型的Ebinding均減小，說(shuō)明溫度的降低可以減小晶體內(nèi)主客體分子間的相互作用能，提高晶體的穩(wěn)定性的：模型A的Ebinding降低最大，說(shuō)明溫度的降低對(duì)其影響最大；模型B的Ebinding受溫度影響較小，進(jìn)一步證實(shí)了當(dāng)水合物處于穩(wěn)定狀態(tài)時(shí)，降低溫度對(duì)其穩(wěn)定性的影響并不大。

3 結(jié) 論

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

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

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

（3）對(duì)于不穩(wěn)定晶體，溫度的降低可以將水分子的均方位移降低到0附近，減小籠型結(jié)構(gòu)的變形，水合物晶體的分解速率降低，穩(wěn)定性提高；對(duì)于已經(jīng)穩(wěn)定的晶體，溫度的降低對(duì)已經(jīng)接近為0的均方位移影響不大，對(duì)水合物穩(wěn)定性的影響不明顯。

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

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

符號(hào)說(shuō)明

A——晶格參數(shù)，m

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

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

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

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

e——基本電荷，C

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

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

N——體系中的分子總數(shù)，個(gè)

qH——?dú)湓拥碾姾闪颗c基本電荷的比值，量綱為1

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

Rcutoff——截?cái)喟霃?，m

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

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

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

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

T——溫度，K

t——時(shí)間，s

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

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

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

ε——Lennard－Jones勢(shì)能參數(shù)，J／mol

δ——Lennard－Jones勢(shì)能參數(shù)，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］ 裘俊紅，郭天民.甲烷水合物在含抑制劑體系中的生成動(dòng)力學(xué)［J］.石油學(xué)報(bào)：石油化工，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］ 丁麗穎，耿春宇，趙月紅，等.Ⅰ型甲烷水合物晶體穩(wěn)定性的分子動(dòng)力學(xué)模擬［J］.計(jì)算機(jī)與應(yīng)用化學(xué)，2007，24（5）：570.

［13］ 萬(wàn)麗華，梁德青.客體分子數(shù)對(duì)甲烷水合物導(dǎo)熱性能影響的分子動(dòng)力學(xué)模擬［J］.化工學(xué)報(bào)，2012，63（2）：382－386.

［14］ 萬(wàn)麗華，顏克鳳，魯濤，等.分子動(dòng)力學(xué)模擬甲烷水合物結(jié)構(gòu)穩(wěn)定性［J］.武漢理工大學(xué)學(xué)報(bào)，2010，32（2）：25－27.

［15］ 馬應(yīng)海，茍?zhí)m濤，何曉霞，等.四氫呋喃水合物零度以上生成動(dòng)力學(xué)研究［J］.天然氣地球科學(xué)，2006（2）：244－248.

［16］ 涂運(yùn)中，蔣國(guó)盛，張凌，等.SDS和THF對(duì)甲烷水合物合成影響的實(shí)驗(yàn)研究［J］.現(xiàn)代地質(zhì)，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，F(xiàn)owler 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.