分子模型SPC/Fw-EPM2對CO2水合物動力學及熱力學特性表征

蘇保照 焦麗君,2 李裕斌 李 倩

（1.青島職業技術學院海爾學院 青島 266555；2.中國石油大學新能源學院 青島 266555）

0 引言

CO2水合物由水和二氧化碳分子組成，水分子之間通過氫鍵相連形成籠形結構，二氧化碳分子位于籠子中心，依靠范德華力相互作用[1]。CO2水合物能夠在0-10℃溫度條件下生成，與空調運行溫度相適宜，而且相變潛熱較大，約為500kJ/kg，遠大于冰的分解焓，是一種儲能性能良好的相變材料，未來極有可能取代冰蓄冷，解決空調儲能需要雙工況機組的技術難題[2,3]。而且與N2、H2和CH4等水合物相比，相平衡條件更加溫和，用于煙氣中二氧化碳捕集，能夠解決傳統碳捕集技術費用過高的問題[4-6]，為碳中和、碳達峰目標的實現提供技術支撐。

然而，由于CO2水合物合成過程存在誘導期長、過冷度大和生長緩慢等問題，限制了相關技術應用。為此，研究人員通過攪拌、鼓泡等動態強化技術和使用表面活性劑、納米粒子等靜態強化手段促進水合物合成，發現實驗條件對動力學促進效果影響較大，已有實驗研究中關于CO2水合物合成的動力學性質差別較大，而且作用機理尚不明確[7]。對于水合物分解動力學過程，研究人員發現傳熱方式、加熱溫度、降壓速度、添加劑種類和濃度等對水合物分解速率影響較大，但是分解過程動力學方程大多只能較好的擬合特定實驗，不具有普遍適用性[8-12]。由于實驗條件的限制，只能觀察實驗現象和宏觀動力學規律，若要洞悉水合物動力學行為背后的非平衡態傳熱傳質規律，還需要從分子水平進一步研究。

分子動力學模擬作為一種微觀研究手段，能夠從分子尺度表征微觀特性，但是對分子模型的依賴性較強[13]。水分子模型主要有剛性分子模型SPC/E、TIP4P/Ice 和TIP4P/2005，柔性分子模型SPC/Fw，二氧化碳分子模型包括TraPPE 和EPM2[14]。其中，柔性分子模型SPC/Fw 在鍵長伸縮和鍵角彎曲兩個方面具有自由度，相比于剛性分子模型，對溫度的描述更加準確，分子模型作用下的熱力學性質或許也更加準確。為此，使用分子模型SPC/Fw 對CO2水合物的微觀特性進行表征十分有必要。

1 模擬方法

1.1 模擬體系

為模擬CO2水合物熱力學及動力學特性，運用Material Studio 軟件建立水-水合物-二氧化碳三相體系，如圖1 所示。左側區域包括368 個水分子，采用棒狀模型，深色表示氧原子，淺色表示氫原子，氫鍵由虛線表示。右側區域包括192 個二氧化碳分子，采用球狀模型，深色表示氧原子，淺色表示碳原子。中間區域由2×2×6 超晶胞二氧化碳水合物組成，包括192 個二氧化碳分子和1104 個水分子，其中氧原子和氫原子按照Takeuchi 等報道的笛卡爾坐標放置，二氧化碳分子隨機放置在籠子中心[15]。模擬體系沿x、y 和z軸的尺寸為24.06?×24.06?×128.72?。分子模型選用SPC/Fw和EPM2，具體參數見表1 和表2[16,17]。分子間作用力包括范德瓦爾力、靜電力和鍵長伸縮、鍵角彎曲作用力，如公式（1）所示，不同原子之間的勢能參數采用Lorentz-Berthelot 規則計算，如公式（2）所示[18]。

圖1 模擬體系初始時刻構型Fig.1 Initial configuration of the simulated system

表1 勢能參數Table 1 Potential energy parameters

表2 鍵參數Table 2 Bond Parameters

式中，q為電荷，ε0為真空介電常數，r為原子之間的距離，ε和σ為勢能參數，k為彈力常數，r和θ分別為鍵長和鍵角。

1.2 模擬過程

分子動力學模擬運用開源軟件LAMMPS 進行，其中data 文件由模擬體系轉化成坐標數據得到，in 文件設置參數為：采用周期性邊界條件；分子間相互作用力采用長程靜電力，截斷半徑為12?，計算方法采用pppm（particle-particle particle-mesh），計算精度取10-6；時間步長為1fs，使用velocity-verlet 積分算法求解牛頓運動方程；溫度和壓力控制使用Nose-Hoover恒溫器和恒壓器處理，阻尼系數分別設為0.1ps 和1ps。模擬步驟為：首先使用共軛梯度算法對體系進行能量最小化，然后將體系置于NPT 系綜下弛豫100ps，使體系達到目標壓力和溫度，最后將體系置于NVE 系綜下運行約5ns。

2 結果與討論

2.1 微觀結構特性

徑向分布函數表示局部密度與體系密度的比值，可以表征體系的微觀結構，定義如公式（3）所示。

式中，N為分子總數，r為距離“參考分子”中心的距離，ΔN為介于r到r+δr之間分子數，ρ為系統密度，δr為設置的距離差，T為計算的總時間。水中氧原子間徑向分布函數用go-o(r)表示，二氧化碳中碳原子間徑向分布函數用gc-c(r)表示，結果如圖2 所示。可以看出，對于水中氧原子間徑向分布函數，當r=2.73?時，go-o(r)最高，即處于第一個峰值，與水合物中氫鍵形成的氧原子距離一致，當r=4.47?和6.45?時，處于第二和第三個峰值，與水合物中四面體網絡中氫鍵結構以及六元環中氧原子間距離相一致。而且與Kondori 等[13]報道的結果基本一致（三個峰值分別為2.75、4.5 和6.53?）。對于二氧化碳中碳原子間徑向分布函數，當r=6.63?時，gc-c(r)最高，表示籠子中二氧化碳分子之間的距離，由于分子振動或轉動，峰值在一定范圍內波動（6.45～6.87?），與文獻報道一致[19]。這表明使用分子模型SPC/Fw-EPM2 能夠準確描述CO2水合物的微觀結構特性。不止是徑向分布函數，由LAMMPS 軟件統計輸出沿z 軸的分子密度分布也能證明模型的可靠性，如圖3 所示，二氧化碳和水分子均有11 對高低峰，與圖1 所示體系構型一致，液態水區和二氧化碳集中區域也呈現出質量密度不變的特點。

圖2 初始時刻不同原子間徑向分布函數Fig.2 The radial distribution function between different atoms at the initial moment

圖3 初始時刻不同分子沿z 軸質量密度分布Fig.3 Mass density distribution of different molecules along the z-axis at the initial moment

2.2 分解動力學特性

本征動力學理論認為水合物在分解過程中存在活化能壘，可以作為定量描述水合物動力學特性的基本參數。計算方法如公式（4）和（5）所示[20]。

為擬合不同溫度條件下CO2水合物分解速率，運用SPC/Fw-EPM2 分子模型分別在3.5MPa、300-320K 初始熱力學條件下模擬，類水合物二氧化碳分子數目由Myshakin 等[22]報道的判別方法計算，數目隨時間變化如圖4 所示?？梢钥闯觯跏紲囟仍礁?，水合物分解速率越快，這與普遍的實驗結果一致。對不同溫度條件下，類水合物二氧化碳分子數目隨時間變化作線性擬合，得到分解速率，如表3 所示。將公式（4）和（5）合并，得到分解速率與溫度的依賴性表達式及活化能計算方法，如公式（6）所示。

表3 溫度與分解速率數據統計Table 3 Statistics of temperature and decomposition rate

圖4 不同熱力學條件下類水合物二氧化碳分子數目隨時間變化Fig.4 Number of hydrate-like carbon dioxide molecules change with time under different thermodynamic conditions

由此將分解速率自然對數與溫度倒數進行擬合，結果如圖5 所示，擬合優度為0.9790，擬合斜率為-6093.82±514.94，求得活化能ΔE為50.66±4.28kJ/mol，與實驗值在同一數量級（73kJ/mol）[23]，而且比前期使用SPC/E-EPM2 分子模型計算的活化能（18.48±0.45kJ/mol）更接近實驗值[14]，一方面表明分子模型SPC/Fw-EPM2 在描述CO2水合物的動力學特性方面的適用性和準確性，另一方面，也體現了模擬體系和研究方法的正確性。

圖5 分解速率對初始溫度的依賴性Fig.5 Dependence of decomposition rate on initial temperature

2.3 熱力學性質

盡管水合物的熱力學特性已經獲得了相對比較精確的研究，但是在模擬過程中由于熱力學條件的改變，相變動力學特性差別比較大，這一方面可能與分解界面的性質有關，另一方面也與分子模型的穩定性有關，即不同分子模型表征的水合物相平衡條件差別很大。Costandy 等使用四力點模型水分子模型TIP4P/Ice 和TIP4P/2005 及二氧化碳分子模型TraPPE，對CO2水合物相平衡特性進行了模擬和預測，盡管取得了較好的實驗一致性，但是由于TIP4P/Ice 和TIP4P/2005 均屬于剛性分子模型，在需要溫度梯度和顯熱熱流描述的模擬研究中并不適用。水合物的相平衡特性可以根據體系勢能隨時間變化的行為確定，對于圖1 所示三相模擬體系，若保持壓力不變，當溫度高于相平衡溫度時，水合物分解，分子間隨著距離的增大相互作用力減小，表現為勢能上升，反之，新相水合物合成放熱，勢能下降。但是由于合成過程成核所需誘導期長，受到計算資源限制，僅模擬5ns 時間，故僅針對體系的熱力學穩定性作相關說明。以3.5MPa、275-290K初始條件為例，模擬體系的勢能和溫度變化如圖6所示。可以看出，290K 初始溫度下，體系勢能顯著升高、溫度顯著下降，說明水合物分解。275K初始溫度下，體系勢能和溫度均不變，說明水合物處于穩定區。278 和280K 初始條件下，體系勢能在分解前期緩慢上升，隨著溫度的下降，當降到平衡溫度附近時，沒有足夠的溫差為分解過程提供動力，而且隨著二氧化碳氣體的釋放，體系壓力增大，水合物分解停止。從圖中看，平衡溫度約在275K左右，比相似壓力下測得的實驗值略低（3.395-3.667MPa、280.34-280.79K）[24]。

圖6 不同初始條件下（a）模擬體系勢能和（b）溫度隨時間變化Fig.6 (a)Potential energy and(b)temperature of the simulated system change with time under different initial conditions

Baghel 等[25]提出平衡溫度可以使用公式（7）計算[24]。

式中，T為瞬時溫度，Teq為平衡溫度，T0為初始溫度，k為擬合系數。以310K 初始溫度條件為例求解平衡溫度，擬合結果如圖7 所示，得到平衡溫度為279.85±0.33K，與實驗結果比較接近。正是由于柔性分子模型SPC/Fw 在鍵長伸縮和鍵角彎曲兩個方面具有自由度，相比于剛性分子模型，對溫度的描述更加準確，分子模型作用下的熱力學性質也更加準確。這一點在Li 等[26]對甲烷水合物的研究中也得到了證實。對于310K 初始溫度，水合物在4000ps 時已分解結束，體系壓力增大，因此平衡溫度也比275K 略大。

圖7 310K 初始溫度條件下模擬體系平衡溫度Fig.7 Equilibrium temperature of simulated system at the initial temperature of 310K

使用Clausius-Clapeyron 方程求解分解熱，如公式（8）所示，相平衡溫度和壓力關系取自文獻[24]，由公式（7）擬合后得到斜率為-10323。

式中，ΔHd為分解熱；R為通用氣體常數，大小為8.314J·mol-1·K-1；p為壓力，MPa；T為溫度，K；Z為壓縮因子，利用PR 方程計算[27]，大小為0.6954，分解熱為59.68kJ·mol-1，在文獻報道的數值范圍內（54.5-60.1kJ·mol-1）[28]。

此外，建立純水合物結構體系，運用平衡分子動力學模擬方法，借助Green-Kubo 公式，編程計算了水合物熱導率。計算公式如下：

式中，kB為波爾茲曼常數，大小為1.38×10-23J/K；V為模擬區域的體積，m3；T為模擬區域的溫度，K；為熱流自相關函數（HCACF）。270K、3.5MPa 條件下，熱導率隨時間變化如圖8 所示，后期熱導率趨于穩定，平均熱導率為0.84W·m-1·K-1，在Jiang 等[29]報道的二氧化碳水合物熱導率范圍內（0.75-0.95）W·m-1·K-1，表明分子模型SPC/Fw-EPM2 能夠較好的描述水合物的熱力學性質。

圖8 270K、3.5MPa 條件下熱導率隨時間變化Fig.8 Thermal conductivity changes with time at 270K,3.5Mpa

3 結論

考慮到柔性分子模型在溫度描述中的精確性，運用分子模型SPC/Fw-EPM2，在不同熱力學條件下，通過分子動力學模擬對CO2水合物相變過程動力學行為和熱力學特性進行了微觀表征，從微觀結構特性、分解活化能和熱力學穩定性三個方面表述了分子模型的適用性，為基于溫度梯度和顯熱熱流統計的漲落耗散分析提供了理論基礎，為水合物動力學行為的微觀機理探索提供了理論支撐。具體結論包括：

（1）分子模型SPC/Fw-EPM2 作用下，水中氧原子間徑向分布函數的第一、第二和第三峰值分別為2.73、4.47 和6.45?，二氧化碳中碳原子間徑向分布函數的峰值在6.45-6.87?之間波動。分子模型SPC/Fw-EPM2 能夠準確描述CO2水合物的微觀結構。

（2）3.5MPa、300-320K 初始熱力學條件下，CO2水合物分解速率對初始溫度的依賴性較強，活化能ΔE為50.66±4.28kJ/mol，分子模型SPC/Fw-EPM2 能夠準確描述CO2水合物的動力學特性。

（3）3.5MPa、270K 條件下CO2水合物處于熱力學穩定區域，3.5MPa、290K 條件下CO2水合物分解。310K 初始溫度條件下，模擬體系平衡溫度為279.85±0.33K，分解焓為59.68kJ·mol-1。3.5MPa、270K 條件下CO2水合物的平均熱導率為0.84W·m-1·K-1，均與實驗結果相接近，表明分子模型SPC/Fw-EPM2 能夠準確描述CO2水合物的熱力學性質。