壓力變化對Ⅰ型甲烷水合物穩定性影響的分子動力學模擬

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

（1.中國石油大學（華東）化學工程學院，山東 青島 266580；2.中國石化 勝利石油管理局鉆井工藝研究院，山東 東營 257017）

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

近年來，國內外研究人員針對水合物開展了大量實驗和理論研究，水合物在微觀水平的成核機制、分解機制以及抑制劑對水合物生成及分解的影響機理是目前的研究熱點［2］。基于經驗勢的分子動力學模擬是研究水合物的有效方法［3-5］，Saruprias等［6］用分子動力學模擬CO2水合物的分解過程，指出客體分子晶穴占據率和溫度效應是影響水合物穩定性的重要因素；Roger［7］用分子動力學方法研究了客體分子對水合物穩定性的影響；丁麗穎等［8-9］對Ⅰ型甲烷水合物的穩定性進行了分子動力學研究，測試了不同溫度及客體分子占據率下甲烷水合物的穩定性，提出了判斷水合物穩定性的依據；Takeshi等［10］模擬了客體分子晶穴占據率對水合物穩定性的影響規律。目前，采用分子動力學方法研究水合物穩定性主要集中在溫度、客體分子及其晶穴占據率等方面，對壓力條件的影響研究較少。

本工作采用分子動力學方法研究了壓力變化對Ⅰ型甲烷水合物穩定性的影響規律，對體系構象、徑向分布函數、勢能變化及均方位移等參數進行了分析，從微觀角度揭示了水合物降壓分解機理。

1 計算模型

1.1 模擬體系

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

1.2 計算方法

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

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

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

2 結果與討論

2.1 體系構象分析

對Ⅰ型甲烷水合物超晶胞體系進行分子動力學模擬，得到不同壓力下體系的最終構象結果如圖2所示。

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

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

圖3展示了10 MPa、280 K條件下，由水分子組成的籠狀孔穴結構的解離過程。由圖3可見，當模擬進行至50 ps時，籠狀孔穴結構發生輕微變形，但基本保持穩定；隨著模擬的進行（50～120 ps），水分子擴散加劇，籠狀孔穴結構破裂開口，封存在籠狀孔穴結構內的甲烷分子逸出，成為自由分子。

圖3 籠狀孔穴結構的變化Fig.3 Structural change of the hydrate cage.p=10 MPa.

2.2 徑向分布函數

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

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

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

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

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

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

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

2.3 勢能變化

模擬過程中體系的勢能變化見圖7。由圖7可見，體系勢能在各個壓力下都能達到平衡，且壓力越低，體系越不穩定，勢能越大。當體系壓力較高（40，30，20 MPa）時，勢能在某一位置波動，表明在此條件下體系具有穩定的晶體結構，水分子和甲烷分子在各自的平衡位置振動和轉動，水合物未發生分解。當壓力降至10 MPa時，0～100 ps內隨著模擬的進行，勢能在緩慢增加，結合2.1節中對體系構象的分析，這是因為此階段籠狀孔穴結構正在解離，水分子排列由有序變為無序；固相水合物完全分解為氣液兩相后，勢能不再增加，最終穩定在某一位置波動。

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

2.4 均方位移

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

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

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

3 結論

1）體系壓力越低，穩定后的勢能越大，體系越不穩定。

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

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

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

［1］陳光進，孫長宇，馬慶蘭.氣體水合物科學與技術［M］.北京：化學工業出版社，2007：305-309.

［2］劉源.氣體水合物分子間相互作用及分解行為的理論研究［D］.大連：大連理工大學.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型甲烷水合物結構穩定性的分子動力學模擬［J］.化工學報，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］丁麗穎，耿春宇，趙月紅.Ⅰ型甲烷水合物晶體穩定性的分子動力學模擬［J］.計算機與應用化學，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，Ferdows M.Monte Carlo Approach to Structure and Thermodynamic Property of CO2Hydrate［J］.JSME Int J B，2005，48（4）：802-809.

［17］萬麗華，顏克鳳，李小森，等.熱力學抑制劑作用下甲烷水合物分解過程的分子動力學模擬［J］.物理化學學報，2009，25（3）：486-494.

［18］陳正隆，徐為人，湯立達.分子模擬的理論與實踐［M］.北京：化學工業出版社，2007：112-114.