段小龍，任 紅，高光海，許俊良，仇性啟

（1.中國(guó)石油大學(xué)（華東）化學(xué)工程學(xué)院，山東 青島 266580；2.中國(guó)石化 勝利石油管理局鉆井工藝研究院，山東 東營(yíng) 257017）

氣體水合物是水與甲烷、乙烷、CO2和H2S等小分子氣體在高壓、低溫環(huán)境中形成的非化學(xué)計(jì)量性的籠狀晶體物質(zhì)。空的水合物晶格就像一個(gè)高效的分子水平的氣體存儲(chǔ)器，每立方米水合物可儲(chǔ)存160～180 m3天然氣［1］。隨著世界油氣資源的日益枯竭，天然氣水合物作為人類(lèi)的潛在能源在世界范圍內(nèi)受到高度重視。

近年來(lái)，國(guó)內(nèi)外研究人員針對(duì)水合物開(kāi)展了大量實(shí)驗(yàn)和理論研究，水合物在微觀水平的成核機(jī)制、分解機(jī)制以及抑制劑對(duì)水合物生成及分解的影響機(jī)理是目前的研究熱點(diǎn)［2］?；诮?jīng)驗(yàn)勢(shì)的分子動(dòng)力學(xué)模擬是研究水合物的有效方法［3-5］，Saruprias等［6］用分子動(dòng)力學(xué)模擬CO2水合物的分解過(guò)程，指出客體分子晶穴占據(jù)率和溫度效應(yīng)是影響水合物穩(wěn)定性的重要因素；Roger［7］用分子動(dòng)力學(xué)方法研究了客體分子對(duì)水合物穩(wěn)定性的影響；丁麗穎等［8-9］對(duì)Ⅰ型甲烷水合物的穩(wěn)定性進(jìn)行了分子動(dòng)力學(xué)研究，測(cè)試了不同溫度及客體分子占據(jù)率下甲烷水合物的穩(wěn)定性，提出了判斷水合物穩(wěn)定性的依據(jù)；Takeshi等［10］模擬了客體分子晶穴占據(jù)率對(duì)水合物穩(wěn)定性的影響規(guī)律。目前，采用分子動(dòng)力學(xué)方法研究水合物穩(wěn)定性主要集中在溫度、客體分子及其晶穴占據(jù)率等方面，對(duì)壓力條件的影響研究較少。

本工作采用分子動(dòng)力學(xué)方法研究了壓力變化對(duì)Ⅰ型甲烷水合物穩(wěn)定性的影響規(guī)律，對(duì)體系構(gòu)象、徑向分布函數(shù)、勢(shì)能變化及均方位移等參數(shù)進(jìn)行了分析，從微觀角度揭示了水合物降壓分解機(jī)理。

1 計(jì)算模型

1.1 模擬體系

采用Ⅰ型甲烷水合物的2×2×2超晶胞作為模擬體系（如圖1所示），在x，y，z 3個(gè)坐標(biāo)方向上使用周期性邊界條件，控制晶胞尺寸為2.452 nm×2.452 nm×2.452 nm，體系包含368個(gè)水分子和64個(gè)甲烷分子。采用文獻(xiàn)［11］報(bào)道的方法搭建Ⅰ型水合物籠狀孔穴結(jié)構(gòu)，氧原子初始位置根據(jù)X單晶衍射實(shí)驗(yàn)［12-13］結(jié)果確定，晶格中氫原子按無(wú)序性排列，使其滿足Bernal-Fowler規(guī)則［14］。

1.2 計(jì)算方法

模擬在Materials Studio平臺(tái)下進(jìn)行，采用NPT系綜，選用相容化合價(jià)力場(chǎng)［15］，使用Ewald加和法計(jì)算范德瓦爾斯相互作用和長(zhǎng)程靜電相互作用，并分別采用Nose熱溶方法和Berendsen方法控制體系的溫度和壓力。本工作使用Lennard-Jones勢(shì)能函數(shù)描述甲烷分子與水分子之間以及甲烷分子之間的相互作用；采用簡(jiǎn)單點(diǎn)電荷勢(shì)能模型控制水分子之間的相互作用，將水分子視為剛體分子，其鍵長(zhǎng)與鍵角固定。在進(jìn)行分子動(dòng)力學(xué)模擬前，先用最速下降法和共軛梯度法對(duì)體系結(jié)構(gòu)進(jìn)行幾何優(yōu)化，使體系的初始結(jié)構(gòu)處于能量最小狀態(tài)。模擬過(guò)程的時(shí)間步長(zhǎng)取1 fs，每個(gè)溫壓點(diǎn)的模擬時(shí)長(zhǎng)取1 ns。

圖1 Ⅰ型甲烷水合物的初始構(gòu)象Fig.1 Initial structure of SI methane hydrate.

本工作主要對(duì)溫度為280 K、晶格占據(jù)率為1、壓力分別為10，20，30，40 MPa下的Ⅰ型甲烷水合物體系進(jìn)行分子動(dòng)力學(xué)模擬，分析不同壓力條件下體系構(gòu)象、徑向分布函數(shù)、體系勢(shì)能和均方位移等參數(shù)的變化，并總結(jié)出相關(guān)規(guī)律。

2 結(jié)果與討論

2.1 體系構(gòu)象分析

對(duì)Ⅰ型甲烷水合物超晶胞體系進(jìn)行分子動(dòng)力學(xué)模擬，得到不同壓力下體系的最終構(gòu)象結(jié)果如圖2所示。

圖2 不同壓力下體系的最終構(gòu)象Fig.2 Final structures of the mathane hydrate under different pressure.

由圖2可看出，壓力分別為40 MPa和30 MPa時(shí)，水分子基本保持有序性排列，由水分子構(gòu)成的籠狀孔穴結(jié)構(gòu)僅出現(xiàn)輕微的變形；壓力為20 MPa時(shí)，水合物的籠狀孔穴結(jié)構(gòu)變形程度較嚴(yán)重，同時(shí)少量的甲烷分子從籠狀孔穴中逸出并出現(xiàn)聚集傾向；當(dāng)壓力降至10 MPa時(shí)，由水分子構(gòu)成的有序晶體結(jié)構(gòu)已經(jīng)不存在，甲烷分子聚集在一起，體系呈氣液兩相分離的狀態(tài)。

圖3展示了10 MPa、280 K條件下，由水分子組成的籠狀孔穴結(jié)構(gòu)的解離過(guò)程。由圖3可見(jiàn)，當(dāng)模擬進(jìn)行至50 ps時(shí)，籠狀孔穴結(jié)構(gòu)發(fā)生輕微變形，但基本保持穩(wěn)定；隨著模擬的進(jìn)行（50～120 ps），水分子擴(kuò)散加劇，籠狀孔穴結(jié)構(gòu)破裂開(kāi)口，封存在籠狀孔穴結(jié)構(gòu)內(nèi)的甲烷分子逸出，成為自由分子。

圖3 籠狀孔穴結(jié)構(gòu)的變化Fig.3 Structural change of the hydrate cage.p=10 MPa.

2.2 徑向分布函數(shù)

徑向分布函數(shù)是衡量凝聚態(tài)物質(zhì)有序性的重要參數(shù)。水分子中的氧原子和甲烷分子中的碳原子的徑向分布函數(shù)分別見(jiàn)圖4和圖5。由圖4可看出，氧原子的徑向分布函數(shù)分別在0.275 nm和0.450 nm處出現(xiàn)兩個(gè)特征峰：第一特征峰代表以氫鍵相連的兩個(gè)氧原子之間的距離，即籠狀結(jié)構(gòu)中五邊形的邊長(zhǎng)；第二特征峰代表水合物中相隔一個(gè)或多個(gè)水分子的氧原子的距離。這與文獻(xiàn)［16-17］報(bào)道的結(jié)果一致。

圖4 氧原子的徑向分布函數(shù)Fig.4 Radial distribution functions of oxygen atom(gO—O(r))under different pressure.

圖5 碳原子的徑向分布函數(shù)Fig.5 Radial distribution functions of carbon atom(gC—C(r))under different pressure.

隨著壓力的下降，氧原子徑向分布函數(shù)特征峰展寬變寬，峰值減小，說(shuō)明甲烷水合物中水分子排列的有序性降低，晶體結(jié)構(gòu)受到破壞；當(dāng)體系壓力為40，30，20 MPa時(shí)，氧原子徑向分布函數(shù)的特征峰值較接近，這是由于體系壓力較高，水合物中的籠狀孔穴結(jié)構(gòu)尚未解離，水分子依然呈有序性排列；當(dāng)體系壓力降至10 MPa時(shí)，徑向分布函數(shù)特征峰高度下降程度最為明顯，說(shuō)明水合物籠狀孔穴結(jié)構(gòu)被破壞，水合物晶體結(jié)構(gòu)逐漸分解。

由圖5可以看出，體系壓力為40，30，20 MPa時(shí)，碳原子的徑向分布函數(shù)在0.690 nm處出現(xiàn)第一特征峰值，代表了兩個(gè)籠狀孔穴中碳原子之間的距離［17］，說(shuō)明此時(shí)水合物籠狀孔穴結(jié)構(gòu)穩(wěn)定，甲烷分子分別處于各自籠狀的孔穴中；而當(dāng)壓力降至10 MPa時(shí)，碳原子徑向分布函數(shù)第一特征峰值移至0.410 nm處，代表甲烷分子從籠狀孔穴中逸出并聚集后的碳原子之間的距離。

圖6為模擬過(guò)程中不同時(shí)段的碳原子徑向分布函數(shù)。由圖6可看出，碳原子的徑向分布函數(shù)在最初階段（0～50 ps）只在0.680 nm處有一特征峰值，此時(shí)水合物保持穩(wěn)定狀態(tài)，甲烷分子處于籠狀孔穴中；在模擬的中間階段（50～120 ps），在0.410 nm處開(kāi)始出現(xiàn)另一特征峰值，兩個(gè)特征峰值并存，且隨著模擬的進(jìn)行，0.410 nm處的特征峰值逐漸升高，0.680 nm處的特征峰值逐漸下降，說(shuō)明籠狀孔穴結(jié)構(gòu)正在逐步解離，甲烷分子逸出孔穴并發(fā)生聚集；在模擬的最后階段（120～1 000 ps），0.680 nm處的特征峰值消失，只在0.410 nm處存在一個(gè)特征峰值，此時(shí)水合物已徹底分解，所有甲烷分子聚集在一起。

圖6 不同模擬時(shí)段碳原子的徑向分布函數(shù)Fig.6 Radial distribution functions of carbon atom in different time interval(t).

2.3 勢(shì)能變化

模擬過(guò)程中體系的勢(shì)能變化見(jiàn)圖7。由圖7可見(jiàn)，體系勢(shì)能在各個(gè)壓力下都能達(dá)到平衡，且壓力越低，體系越不穩(wěn)定，勢(shì)能越大。當(dāng)體系壓力較高（40，30，20 MPa）時(shí)，勢(shì)能在某一位置波動(dòng)，表明在此條件下體系具有穩(wěn)定的晶體結(jié)構(gòu)，水分子和甲烷分子在各自的平衡位置振動(dòng)和轉(zhuǎn)動(dòng)，水合物未發(fā)生分解。當(dāng)壓力降至10 MPa時(shí)，0～100 ps內(nèi)隨著模擬的進(jìn)行，勢(shì)能在緩慢增加，結(jié)合2.1節(jié)中對(duì)體系構(gòu)象的分析，這是因?yàn)榇穗A段籠狀孔穴結(jié)構(gòu)正在解離，水分子排列由有序變?yōu)闊o(wú)序；固相水合物完全分解為氣液兩相后，勢(shì)能不再增加，最終穩(wěn)定在某一位置波動(dòng)。

圖7 體系的勢(shì)能變化Fig.7 Potential energy change of the system.

2.4 均方位移

均方位移反映了體系中的粒子在分子動(dòng)力學(xué)模擬過(guò)程中空間位置與其初始位置的偏離程度。圖8和圖9分別為不同壓力下水分子中氧原子和甲烷分子中碳原子的均方位移隨時(shí)間的變化曲線。從圖8和圖9可看出，隨壓力的下降，氧原子和碳原子的均方位移均增大；當(dāng)壓力為40，30，20 MPa時(shí)，氧原子和碳原子的均方位移均上升至略大于零的數(shù)值后趨于穩(wěn)定，這表明水分子在固定的位置發(fā)生微小振動(dòng)，體系具有穩(wěn)定的晶格結(jié)構(gòu)；而壓力降至10 MPa后，隨著模擬的進(jìn)行，氧原子和碳原子的均方位移均逐漸增加，且均方位移增長(zhǎng)的變化率隨時(shí)間的延長(zhǎng)而逐漸增大，在最后階段（700～1 000 ps）均方位移曲線均趨于一條直線。氣體分子和液體分子在擴(kuò)散過(guò)程中其均方位移與時(shí)間呈線性關(guān)系［18］，這表明甲烷水合物經(jīng)歷了籠狀孔穴結(jié)構(gòu)破裂、甲烷分子逸出并最終分解為氣液兩相的過(guò)程。

圖8 不同壓力下氧原子的均方位移Fig.8 MSD of oxygen atom under different pressure.

圖9 不同壓力下碳原子的均方位移Fig.9 Mean square displacement(MSD) of carbon atom under different pressure.

3 結(jié)論

1）體系壓力越低，穩(wěn)定后的勢(shì)能越大，體系越不穩(wěn)定。

2）隨壓力的降低，氧原子徑向分布函數(shù)特征峰值減小，展寬變寬，水分子排列有序性下降。

3）體系壓力為40，30，20 MPa時(shí)，氧原子和碳原子的均方位移穩(wěn)定在略大于零的某一數(shù)值，體系較穩(wěn)定；當(dāng)壓力降至10 MPa時(shí)，均方位移急劇增大，在最后階段（700～1 000 ps）曲線趨于一條直線。

4）通過(guò)對(duì)體系構(gòu)象、勢(shì)能變化、徑向分布函數(shù)及均方位移等參數(shù)進(jìn)行分析，從微觀角度揭示了Ⅰ型甲烷水合物降壓分解機(jī)理：隨體系壓力的降低，由水分子以氫鍵連接形成的籠狀孔穴結(jié)構(gòu)變形并逐漸解離，甲烷分子從籠狀孔穴結(jié)構(gòu)中逸出并聚集，固相水合物最終分解為氣液兩相。

［1］陳光進(jìn)，孫長(zhǎng)宇，馬慶蘭.氣體水合物科學(xué)與技術(shù)［M］.北京：化學(xué)工業(yè)出版社，2007：305-309.

［2］劉源.氣體水合物分子間相互作用及分解行為的理論研究［D］.大連：大連理工大學(xué).2012.

［3］Brodholt J，Wood B.Molecular-Dynamics Simulations of the Properties of CO2-H2O Mixtures at High-Pressure and Temperatures［J］.Am Mineral，1993，78（5/6）：558-564.

［4］Yan Kefeng，Mi Jianguo，Zhong Chongli.Molecular Dynamics Simulation of Inhibiting Effects on Natural Gas Hydrate［J］.Acta Chim Sin，2006，64（3）：223-228.

［5］梅東海，李以圭，陸九芳，等.H型甲烷水合物結(jié)構(gòu)穩(wěn)定性的分子動(dòng)力學(xué)模擬［J］.化工學(xué)報(bào)，1998，49（8）：662-670.

［6］Saruprias S，Debenedetti P G.Molecular Dynamics Study of Carbon Dioxide Hydrate Dissociation［J］.J Phys Chem A，2009，113（10）：1913-1921.

［7］Roger P M.Stability of Gas Hydrate［J］.J Phys Chem，1990，94（15）： 6080-6089.

［8］丁麗穎，耿春宇，趙月紅.Ⅰ型甲烷水合物晶體穩(wěn)定性的分子動(dòng)力學(xué)模擬［J］.計(jì)算機(jī)與應(yīng)用化學(xué)，2007，24（5）：569-574.

［9］Ding Liying，Geng Chunyu，Zhao Yuehong，et al.Molecular Dynamics Simulation on the Dissociation Process of Methane Hydrates［J］.Mol Simulat，2007，33（12）：1005-1016.

［10］Takeshi Sugahara，Hiroki Mori，Jun Sakamoto，et al.Cage Occupancy of Hydrogen in Carbon Dioxide，Ethane，Cyclopropane，and Propane Hydrates［J］.Open Therm J，2008（2）：1-6.

［11］Storr M T，Taylor P C，Monfort J P，et al.Kinetic Inhibitor of Hydrate Crystallization［J］.J Am Chem Soc，2004，126（5）：1569-1576.

［12］Rawn C J，Rondinone A J，Chakoumakos B C，et al.Neutron Powder Diffraction Studies as a Function of Temperature of Structure Ⅱ Hydrate Formed from Propane［J］.Can J Phys，2003，81（1/2）：431-438.

［13］Gutt C，Asnussen B，Press W，et al.The Structure of Deuterated Methane-Hydrate［J］.J Chem Phys，2000，113（11）： 4713-4721.

［14］Greathouse J A，Cygan R T，Simmons B A.Vibrational Spectra of Methane Clathrate Hydrates from Molecular Dynamics Simulation［J］.J Phys Chem B，2006，110（13）：6428-6431.

［15］Geng Cunyu，Wen Hao，Zhou Han.Molecular Simulation of the Potential of Methane Reoccupation During the Replacement of Methane Hydrate by CO2［J］.J Phys Chem A，2009，113（18）：5463-5469.

［16］Ota M，F(xiàn)erdows M.Monte Carlo Approach to Structure and Thermodynamic Property of CO2Hydrate［J］.JSME Int J B，2005，48（4）：802-809.

［17］萬(wàn)麗華，顏克鳳，李小森，等.熱力學(xué)抑制劑作用下甲烷水合物分解過(guò)程的分子動(dòng)力學(xué)模擬［J］.物理化學(xué)學(xué)報(bào)，2009，25（3）：486-494.

［18］陳正隆，徐為人，湯立達(dá).分子模擬的理論與實(shí)踐［M］.北京：化學(xué)工業(yè)出版社，2007：112-114.