不同結構氣體水合物穩定性的分子動力學研究

尹曉霞,焦麗君,2,李元超,華澤珍,李裕斌

(1.青島職業技術學院海爾學院，山東 青島 266555;2.中國石油大學新能源學院，山東 青島 266580)

在能源短缺、環境污染的大背景下，天然氣水合物作為未來戰略能源的勘探和開采，受到越來越多的關注，可以利用二氧化碳水合物比甲烷水合物穩定的特點[1]，通過置換法開采實現天然氣水合物的開采，同時完成二氧化碳的封存。因此，二氧化碳水合物的穩定性和分解特性的研究尤為重要。此外，由于二氧化碳水合物相平衡壓力相對較高，在海水淡化、空調蓄冷、煙氣捕集等領域的應用受到限制，二氧化碳水合物和環戊烷等熱力學促進劑相結合，節能效果顯著[2-4]。

氣體水合物是一種由小分子氣體或者易揮發的液體與水相互作用，在低溫高壓條件下形成的類冰狀固體。主體水分子之間通過氫鍵形成籠形網絡，客體分子位于晶穴中心，與水分子通過微弱的范德華力連接，形成穩定的水合物[5]。常見的客體分子有：甲烷、乙烷、丙烷、二氧化碳、氫氣、氮氣、環戊烷、四氫呋喃、四丁基溴化銨等。根據客體分子的大小，水合物通常分為三種結構，分別是sI型、sII型和sH型水合物。晶體特性如表1所示[6]?？腕w分子不同，氣體水合物的成核特性、分解特性及穩定性差別較大。Zhang等[7]通過對七種不同客體分子在水溶液中成核過程發現，相同濃度下，客體分子尺寸越大，自擴散系數越低，越容易生成水合物。Kondori等[8]發現，加入熱力學抑制劑后，II型氣體水合物分解特性和穩定性均表現出明顯差異，抑制程度排序為甲醇>乙醇>甘油。張慶東[9]通過對比不同客體分子對氣體水合物穩定性的影響，發現對于含有甲烷的二元水合物體系，二氧化碳、乙烷、硫化氫和丙烷對水合物穩定性的促進作用逐漸增加。表明二元水合物可以增強甲烷水合物的穩定性。王浩瑄等[10]發現，客體分子所處晶穴結構不同，對甲烷水合物的穩定性影響較大，大晶穴內部的客體分子作用更加顯著。目前關于氣體水合物的穩定性研究大多集中在微觀模擬，而且以甲烷水合物為主，二氧化碳水合物研究較少。

表1 氣體水合物晶體特性

二氧化碳氣體分子直徑為0.33 nm，通常形成Ⅰ型水合物，可以同時占據大小晶穴[11]。環戊烷、四氫呋喃等分子直徑較大，通常占據Ⅱ型水合物的大晶穴[12]。當二氧化碳和環戊烷(CP)客體分子同時存在時，可以形成二元水合物，Lee和Kim通過X射線衍射實驗發現，環戊烷完全占據大晶穴，二氧化碳占據一半以上的小晶穴。他們指出，這種二元的水合物結構比單一二氧化碳水合物更加穩定[13]。目前的研究大多集中在宏觀實驗層面，缺乏微觀視角，本文從微觀角度出發，從分子水平上研究不同結構的水合物在相同熱力學條件下的穩定性。

1 模擬方法

為了研究不同結構氣體水合物的穩定性，本文建立了Ⅰ型和Ⅱ型二氧化碳水合物、Ⅱ型環戊烷水合物、以及二氧化碳-環戊烷二元水合物分子模型，如圖1所示。水合物分子模型通過Materials Studio軟件建立，根據X射線衍射結果確定水分子內氧原子和氫原子的空間坐標，然后將客體分子的中心放置在水分子構建的籠形結構中心[14]。為了保證晶穴數目相當，對于I型水合物，構建3X3X4超晶胞，其中包括1 565個水分子和288個晶穴；對于II型水合物，構建2X2X3超晶胞，其中包括1 632個水分子和288個晶穴。

圖1 水合物的初始構型 (圖中線表示氫鍵，球體表示氧原子和碳原子)

分子動力學模擬基于牛頓力學為理論基礎，通過分子勢能函數產生分子間相互作用，計算原子受力及加速度，在非常短的時間間隔內，得到各分子的速度及位置。因此，分子勢能函數的選擇對模擬結果的影響較大。本文中水分子模型選用三點模型SPC/E，該模型簡單且高效，能夠很好的預測液態水的自擴散系數，而且在水合物相平衡溫度、分解焓和活化能的預測等方面也能與實驗值較好的吻合[15-16]。二氧化碳分子模型采用EPM2模型，該模型能夠準確預測二氧化碳氣體在水中的溶解度以及氣液共存線[17-18]。環戊烷分子模型直接采用Materials Studio軟件庫中的CVFF力場模型。各分子模型的勢能參數如表2所示。

表2 各分子模型參數

分子間相互作用力包括Lennard-Jones勢能和庫侖力，見公式(1)

(1)

不同原子間的Lennard-Jones勢能作用遵從Lorentz-Berthelot作用規則，見公式(2)(3)

(2)

(3)

模擬使用LAMMPS軟件，時間步長設為1 fs，采用速度-Verlet 算法，應用周期性邊界條件消除界面影響。范德華力的作用截止半徑設為12?，采用PPPM(Particle-Particle Particle-Mesh)方法來計算長程靜電力，計算精度為10-6。采用Nose-Hoover控溫和控壓，溫度和壓力阻尼系數分別設為0.1 ps和1 ps，將水合物體系在控制溫度和壓力的條件下(NPT系綜)，模擬200 ps。

2 結果與討論

2.1 徑向分布函數

徑向分布函數(Radial Distribution Function(RDF)是指以一個參考原子α為中心，與其距離為r到r+Δr之間，β原子的數目與系統中平均數目的比值。從物理意義上說，可以理解為系統中區域密度與平均密度的比值，可以反映分子的微觀結構。見公式(4)

(4)

式中g(α-β)(r)——原子α和原子β之間的徑向分布函數;

V——模擬盒子的體積;

Nα和Nβ——兩種原子的數目;

r——兩個原子之間的距離;

niβ(r)——在距離原子α為r的條件下，β原子的數目。

圖2為模擬初期t=20 ps時不同水合物體系所對應的O-O原子及C-C原子的徑向分布函數。對于水中氧原子的徑向分布函數，所有體系第一個峰值出現在r=2.79附近，第一個峰值代表水合物晶穴中兩個相鄰氧原子之間的距離，因此所有體系的值均相同。這與文獻結果相差不大，證明了模型的正確性和模型方法的有效性[19-20]。第二和第三個峰值分別代表水合物中相隔一個和兩個水分子的距離，對于Ⅰ型水合物和Ⅱ型水合物有偏差，體現了水合物不同的網格結構，Ⅱ型水合物氧原子平均距離相對遠些。對于Ⅰ型二氧化碳分子中碳原子的徑向分布函數，在穩定狀態下兩個碳原子的最小距離是6.81，其值與文獻相一致[19]。對于Ⅱ型水合物，二氧化碳分子中碳原子在穩定狀態下的距離最小值為6.21，這種差異主要受到水合物結構的影響，Ⅱ型水合物中碳原子距離更近。而且，通過RDF的峰值越高，表明局部密度相對于平均密度越大，原子越密集。因此，從b圖可以看出，二元CO2+CP水合物的穩定性要強于Ⅱ型CO2水合物，這與環戊烷是常用的熱力學促進劑的說法相一致，環戊烷通過占據Ⅱ型水合物的大晶穴，可以有效的改善相平衡條件，使水合物變得更加穩定[21]。對于環戊烷中的碳原子分布，第一個峰值在r=3.93附近，指的是同一個分子中碳原子之間的距離，與水合物結構關系不大，因此研究意義不大。

圖2 270 K,5 MPa條件下t=20 ps時不同體系氧原子和碳原子徑向分布函數

2.2 均方位移

分子動力學模擬是基于牛頓力學原理的，分子受力產生加速度，結果是模擬過程的每一時刻分子的受力和位置都不相同。通常用在某時刻與初始時刻相比，分子位移平方的平均值來表示分子運動的程度，稱為均方位移(Mean Square Displacement MSD)。見公式(5)

(5)

圖3是在270 K 5 MPa熱力學條件下，不同體系水分子和客體分子的均方位移，從圖中可以明顯看出，Ⅱ型CO2水合物最不穩定，在模擬初期(t=30 ps)MSD即開始迅速上升，體系內分子運動劇烈，表明由固態的水合物體系逐漸變為液態，說明水合物發生分解。Ⅰ型水合物在模擬后期(t=150 ps)MSD開始上升，表明體系經歷了一段亞穩定狀態，最終開始分解。而對于Ⅱ型的CP水合物和CO2+CP二元水合物，在模擬期間，MSD值幾乎不變，尤其是水分子的MSD為0，表明水合物結構十分穩定，一直沒有分解。水分子的MSD為0，可以看出在模擬過程中，水分子之間通過氫鍵結構構建水合物晶格，水分子不能旋轉和移動。而客體分子的MSD大概在2左右，表明在水合物中客體分子可以自由的旋轉或振動，在晶穴內部發生小位移的偏移，但不能逃離水合物晶穴。

圖3 270 K,5 MPa條件下不同體系水分子和客體分子均方位移隨時間變化曲線

2.3 F3序參數

為了進一步分析水合物的分解情況和穩定性，定義了序參數[22]。見公式(6)

F3,i=〈(cosθjik|cosθjik|+cos2(109.472))2〉jk

(6)

式中θjik——以氧原子i為中心，與距離其3.5?內兩個氧原子j和k，三者形成的夾角。

借助此式來反映氫鍵網格偏離理想四面體的情況，從而監測水合物的穩定性。水合物狀態下，F3=0，液態水狀態下，F3=0.1。不同體系下F3序參數的隨時間變化趨勢如圖4所示。從圖中可以清晰的看出，結構Ⅱ型CO2水合物最不穩定，模擬初期即開始分解，分解時間大概50 ps，結構I型CO2水合物次之，在模擬中后期開始出現分解趨勢，分解時間大概100 ps，分解相對較慢。因此，結構Ⅰ型水合物相對結構Ⅱ型水合物更加穩定，這也可以解釋分子動力學模擬水合物合成過程中，甲烷和水作用主要生成Ⅰ型水合物，生成Ⅱ型水合物較少[23]。而結構Ⅱ型的CP水合物和二元水合物序參數均保持不變，水合物沒有發生分解，這與均方位移分析結果一致。

圖4 270 K,5 MPa條件下不同體系F3序參數隨時間變化曲線

3 結 論

本文使用分子動力學模擬方法，通過對徑向分布函數進行分析，充分證實了該模擬方法的正確性，同時從微觀層面進一步認識了水合物的基本構型，包括結構Ⅰ型、結構Ⅱ型水合物，以及不同客體分子的微觀構型。此外，碳原子密度分布、系統勢能、構像圖等多種參數均反映了上述性質，由于篇幅限制未完全列出。主要結論及展望概括如下：

(1)通過均方位移和F3序參數，對比分析了不同結構水合物的穩定性，發現對于CO2客體分子，由于尺寸較小，在Ⅰ型水合物構型下相對比較穩定。

(2)相對于CO2水合物，大分子環戊烷水合物更加穩定，而且將環戊烷和二氧化碳分別置于Ⅱ型水合物的大小晶穴后，發現對于二氧化碳水合物而言，加入熱力學促進劑的二元水合物的熱力學穩定性較單一二氧化碳客體分子水合物要好。

(3)本文從分子視角探索水合物分解特性和穩定性的作用機理，為水合物的宏觀實驗研究和工程應用積累了理論基礎，同時在該工作的基礎上，將來可以繼續研究不同構型和不同客體分子水合物的界面分解特性。