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

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

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

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

表1 氣體水合物晶體特性

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

1 模擬方法

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

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

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

表2 各分子模型參數(shù)

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

(1)

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

(2)

(3)

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

2 結(jié)果與討論

2.1 徑向分布函數(shù)

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

(4)

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

V——模擬盒子的體積;

Nα和Nβ——兩種原子的數(shù)目;

r——兩個(gè)原子之間的距離;

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

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

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

2.2 均方位移

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

(5)

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

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

2.3 F3序參數(shù)

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

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

(6)

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

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

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

3 結(jié) 論

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

(1)通過均方位移和F3序參數(shù)，對比分析了不同結(jié)構(gòu)水合物的穩(wěn)定性，發(fā)現(xiàn)對于CO2客體分子，由于尺寸較小，在Ⅰ型水合物構(gòu)型下相對比較穩(wěn)定。

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

(3)本文從分子視角探索水合物分解特性和穩(wěn)定性的作用機(jī)理，為水合物的宏觀實(shí)驗(yàn)研究和工程應(yīng)用積累了理論基礎(chǔ)，同時(shí)在該工作的基礎(chǔ)上，將來可以繼續(xù)研究不同構(gòu)型和不同客體分子水合物的界面分解特性。