不同結構二氧化碳水合物熱物性的分子動力學模擬

李元超 焦麗君 尹曉霞 李 倩 華澤珍

（青島職業技術學院海爾學院 青島 266555）

0 引言

氣體水合物是在一定的溫度和壓力條件下形成、由水分子和客體分子組成的包絡狀固體化合物，水分子之間由氫鍵相連形成籠形結構，客體分子位于晶穴中心。常見的客體分子有甲烷、乙烷、丙烷和二氧化碳等。根據客體分子大小，通常分為結構I 型、II 型和H 型水合物。自然界中如天然氣水合物、二氧化碳水合物一般是I 型水合物，客體分子尺寸為4-5.5 ?；人工合成的水合物如環戊烷水合物、四氫呋喃水合物通常是II 型水合物，客體分子尺寸為6-7 ?；構成H 型水合物的大客體分子尺寸為8-9 ?，H 型水合物相對較少[1]。氣體水合物的結構不同，對應的熱力學性質，如相平衡特性有較大的區別，相比于二氧化碳水合物的相平衡條件（276.45-282.05K，1.81-3.83MPa），二氧化碳和環戊烷二元水合物的相平衡條件更加溫和（281.55-290.25K，0.15-1.92MPa）[2]，具體表現在相平衡溫度更高、壓力更低，能夠克服二氧化碳水合物合成壓力高、過冷度大的問題，使基于二氧化碳水合物的碳捕集、海水淡化以及空調蓄冷的低成本應用成為可能[3-7]。因此，研究不同結構水合物的熱物性對于水合物合成和分解的應用具有重要意義。

氣體水合物的熱物性測量有實驗和模擬兩種方法，目前研究較多的是甲烷水合物，其次是乙烷、丙烷和四氫呋喃水合物，二氧化碳水合物及二元水合物的熱物性報道較少[8-13]。實驗室人工合成的水合物通常具有晶格缺陷和空穴，Lee 等[14]用X 射線衍射方法測試了二氧化碳和環戊烷（CP）二元水合物的晶體結構和晶穴占有率，發現結構為II 型水合物，由于二氧化碳的存在，晶格常數相比于環戊烷水合物略有增加，其中環戊烷完全占據了大晶穴，二氧化碳占據了62%的小晶穴。實驗測量熱物性對于實驗設備的要求比較高，而且很難精確的合成不同結構和不同晶穴占有率的水合物，因此，本文采用分子動力學模擬的方法研究不同結構和晶穴占有率下氣體水合物的熱物性。

1 模型與方法

1.1 分子模型

利用Material Studio 建立不同結構的氣體水合物，結構類型和大小晶穴占有率如表1所示，初始構型如圖1所示，分別為I 型CO2水合物和II型CO2-CP 二元水合物，其中二氧化碳占據小晶穴，環戊烷占據大晶穴，圖1（b）是晶穴占有率為100%的情況。I 型水合物結構為2X·6Y·46H2O，II 型水合物結構為8X·16Y·136H2O，其中X 表示大晶穴，Y 表示小晶穴[1]。為了保證體系原子數大致相同，I 型水合物采用3×3×4 超晶胞結構，II 型水合物采用2×2×3 超晶胞結構，水分子數分別為1565和1632，晶穴數均為288。按照Takeuchi 等[15]提出的I 型和II 型水合物的晶胞常數及絕對坐標，分別放置氧原子和氫原子，并將客體分子放在晶穴中心。水分子、二氧化碳分子和環戊烷分子分別選擇SPC/E 模型、EMP2 模型和全原子力場模型，勢能函數見表2[16-18]。分子間相互作用力包括Lennard-Jones 范德華力和庫侖力，見公式（1），不同原子間的Lennard-Jones 勢能作用遵從Lorentz-Berthelot 作用規則，見公式（2）和（3）[19]。

表2 勢函數參數Table 2 Potential energy parameters

表1 模擬體系Table 1 Simulation system

圖1 初始構型（藍線-氫鍵，棒狀-共價鍵，球形-原子，氧原子-紅色，氫原子-白色，碳原子-黃色）Fig.1 Configuration of initial hydrate（The dotted blue line represents hydrogen bonding,the ball represents atom,the stick represents covalent bond,C-yellow,O-red,H-white）

1.2 模擬方法

采用開源軟件LAMMPS 進行分子動力學模擬。三個方向均采用周期性邊界條件，時間步長為1fs，分子間相互作用力的截斷半徑為12?，采用PPPM（particle-particle particle-mesh）方法計算長程靜電力，精度為10-6，采用velocity-verlet 積分算法求解運動方程，采用shake 算法固定水分子的鍵長和鍵角。首先進行能量最小化，然后將體系置于NVT系綜下運行500ps，采用Noose-Hoover 熱浴，阻尼系數設為0.1ps，使體系達到恒定溫度并保持平衡。接著將體系置于NPT 系綜下運行500ps，采用Noose-Hoover 熱浴和壓浴，溫度阻尼系數設為0.1ps，壓力阻尼系數設為1ps，使體系達到設定溫度和壓力。最后，將體系置于NVE 系綜下運行1ns。模擬結束后，采用平衡分子動力學模擬（EMD）方法，運用Green-Kubo 公式，編程計算熱導率，見公式（4）[20]。

其中，kB為波爾茲曼常數，J/K；V為模擬區域的體積，m3；T為模擬區域的溫度，K；＜J(t)·J(0)＞為熱流自相關函數（HCACF）。

聲子平均自由程是影響物體內部熱輸運的重要參數[21]，李佳等[22]通過對甲烷水合物的分子動力學模擬發現，分解過程的固-液界面厚度與甲烷水合物內部的聲子平均自由程十分接近，說明水合物的導熱特性對水合物分解過程影響較大。根據分子動力學理論，聲子平均自由程的計算方法如公式（5）所示[23]。

其中，k為熱導率，W·m-1·K-1；ρ為密度，kg/m3；βT為等溫壓縮系數，m·s2/kg；ν為聲學速度；，kg/s；cν為等容比熱，kJ/kg·K。以上參數均可以通過分子動力學模擬求解，密度取NPT 系綜下的統計平均值，等溫壓縮系數計算方法如公式（6）所示，在0.5ns 的NPT 系綜下計算系統體積的波動。

其中，V表示體系的體積，m3，T為體系溫度，K，kB為波爾茲曼常數，大小為1.38×10-23J/K。此外，等容比熱cν根據NVT 系綜下能量波動計算，見公式（7），其中E為體系的勢能，J/mol。

為了驗證模擬方法和計算方法的正確性，測試計算了I 型CO2水合物的熱物性，在相似的熱力學條件下，模擬值與文獻值相接近，具體對比將在結果與討論部分展開。

2 結果與討論

2.1 I 型CO2 水合物導熱特性

200K、20MPa 條件下，CO2水合物在模擬初始時刻和結束時刻的徑向分布函數如圖2所示，可以看出，水分子和二氧化碳分子之間均存在明顯的結構分布，水中氧原子之間的徑向分布函數峰值依次出現在2.73?，4.5?和6.46?，分別代表了水合物晶格結構五元環邊長和多面體籠形結構中水分子與對角水分子之間的距離，二氧化碳分子中碳原子之間的徑向分布函數峰值出現在6.76?，代表了兩個晶穴中心的距離，峰值與文獻報道一致[24]，證明了水合物結構的正確性。而且在模擬前后，水合物的結構峰值并未發生變化，只是略有下降，這是由于在分子間力的作用下，水分子和二氧化碳分子在初始位置上振動和旋轉導致的，水合物結構并未發生變化，證明了二氧化碳水合物在200K、20MPa條件下處于熱力學穩定區域，可以使用平衡分子動力學模擬方法（EMD）獲得熱物性。圖3 是I 型CO2水合物在200K、20MPa 條件下的熱流自相關函數（HCACF）和熱導率隨關聯時間變化曲線，可以看出，熱流在0.6ps 內快速收斂，熱導率也在1ps 內趨于穩定，平均熱導率為0.8498W·m-1·K-1。姚貴策[25]通過實驗測得的二氧化碳水合物的熱導率在0.45-0.53W·m-1·K-1（溫度范圍257.15K-273.15K），Jiang 等[26]通過分子動力學模擬測得的二氧化碳水合物的熱導率在0.75-0.95W·m-1·K-1（溫度范圍0-250K），熱導率與實驗和模擬結果均處于同一數量級，證明了模擬方法的可靠性。

圖2 200K，20MPa 條件下不同時刻徑向分布函數圖Fig.2 RDF plots of atoms at different time under 200K,20MPa

圖3 200K,20MPa 條件下（a）熱流自相關函數（b）熱導率隨關聯時間變化Fig.3 (a)HCACF and(b)thermal conductivity varies versus time under 200K,20MPa

圖4 是熱流自相關函數能譜，通過對熱流自相關函數作傅里葉變換得到，能夠反映出聲子能量的消散。可以看出不同頻率的聲子對水合物內部熱輸運的影響不同，高頻成分（＞15THz）大于低頻成分，說明在水合物導熱過程中高頻聲子（15-25THz）起到主導作用，而低頻聲子（＜12THz）幅值相對較小，表明在低頻區域主客體分子之間充分耦合，而高頻區域的聲子更容易消散，即聲學聲子的導熱對水合物導熱起主導作用。I 型二氧化碳水合物的聲子作用模式和甲烷水合物相類似[23]。

圖4 二氧化碳水合物熱流自相關函數能譜Fig.4 Power spectra of HCACF in CO2 hydrate

2.2 II 型CO2-CP 二元水合物導熱特性及對比分析

在200K、20MPa 條件下，計算II 型CO2-CP二元水合物在模擬前后的徑向分布函數，結果表明，模型正確而且水合物處于熱力學穩定區。熱流自相關函數和熱導率分別在1ps 和1.5ps 結束收斂，說明模擬方法對于II 型CO2-CP 二元水合物同樣具有適用性。由于篇幅限制，此處不再展示，主要分析熱導率以及影響水合物熱導率的聲子輸運模式。圖5 是不同條件不同結構水合物的熱導率，可以看出，對于I 型水合物，熱導率的模擬值均大于實驗值[25-27]，熱導率隨溫度變化不存在顯著的正相關或負相關趨勢，對于不同的水合物，表現為在一定的溫度范圍內熱導率隨溫度升高，而在另一范圍熱導率隨溫度降低。

圖5 二氧化碳水合物熱導率Fig.5 Thermal conductivity of carbon dioxide hydrate

為了分析不同結構二氧化碳水合物的導熱特性，圖6 繪制了20MPa、200K 條件下II 型CO2水合物（a）和II 型CO2-CP 水合物（晶穴占有率為50%（b）和100%（c））的熱流自相關函數能譜。可以看出，影響II 型水合物的聲子導熱模式除了高頻聲子外，低頻聲子區域（＜12THz）也存在明顯的峰值，而且當二氧化碳僅占據小晶穴時，低頻聲子區域的能譜峰值較低，隨著環戊烷分子占據大晶穴，能譜峰值提高，當二氧化碳和環戊烷完全占據大小晶穴時（晶穴占有率100%），低頻聲子區域的能譜峰值最高，幾乎接近于高頻聲子區域的能譜峰值。說明對于II 型水合物，隨著大晶穴的晶穴占有率的增加，低頻聲子導熱對水合物導熱起到的作用也變得顯著，高頻聲子和低頻聲子導熱均對水合物導熱起到重要作用。

圖6 二元水合物熱流自相關函數能譜Fig.6 Power spectra of HCACF in CO2-CP hydrate

2.3 不同結構二氧化碳水合物的熱物性

針對表1 列出的4 個不同體系，在50-200K溫度范圍內，計算了不同結構二氧化碳水合物的熱物性，包括密度、等溫壓縮系數、聲學速度、比熱、熱導率、熱擴散率和聲子平均自由程。表4 是不同結構二氧化碳水合物的密度，可以看出，當氣體水合物結構一定時，溫度越低，密度越大，這是由于溫度越低，分子的勢能越低，分子之間的相互作用力越大，水合物越穩定，此時分子熱運動緩慢，客體分子位于晶穴中心，主客體分子間吸引力較大，因此體積較小，密度較大。相同溫度下，不同結構氣體水合物的密度排序為：I 型CO2水合物，II 型CO2+CP 二元水合物（綜合晶穴占有率100%），II型CO2水合物（綜合晶穴占有率66.7%）和II 型CO2+CP 二元水合物（綜合晶穴占有率50%），可見晶穴占有率越高，客體分子對密度的貢獻越大，密度越大。Chialvo 等[28]使用水分子SPC/E 模型計算的I 型CO2水合物在270K、5MPa 下的密度為1.181g/cm3，與200K、20MPa 條件下的模擬結果（1.189g/cm3）相接近，而且滿足溫度越低，密度越大的規律。不同結構不同溫度二氧化碳水合物的等溫壓縮系數如表5所示，可以看出，相同結構的水合物，溫度越低，等溫壓縮系數越低。Ning 等[29]在20MPa、271.15K 條件下模擬了I 型CO2水合物的等溫壓縮系數，約為1.1×10-4MPa-1，與表5 中20MPa、200K 條件下的模擬結果0.8753×10-4MPa-1在同一個數量級，而且滿足溫度越高，等溫壓縮系數越大的規律。

表4 不同結構二氧化碳水合物的密度（g/cm3）Table 4 Densities of carbon dioxide hydrates with different structures

表5 不同結構二氧化碳水合物的等溫壓縮系數（×10-4 MPa-1）Table 5 Isothermal compressibility of carbon dioxide hydrates with different structures

聲學速度是影響固體熱輸運的重要參數，僅與密度和等溫壓縮系數有關，計算結果如表6所示。Jiang 等[26]在0.1MPa、150K 條件下通過分子動力學模擬得到的二氧化碳水合物的聲學速度是3300m/s，與表6 所列20MPa、150K 條件下的聲學速度3321.78m/s 相差不大。可以看出，溫度越低，聲學速度越高，而且I 型CO2水合物（晶穴占有率100%）的聲學速度大于II 型CO2+CP 二元水合物（晶穴占有率100%），這是由于II 型水合物主客體分子之間耦合增強的緣故，當晶穴占有率降低時，主客體之間的耦合作用減弱（見圖6），聲學速度增大[30]。定容比熱的計算結果如表7所示，可以看出，盡管水合物結構不同，但是隨溫度的變化趨勢是一致的，溫度越高，比容越大，這與Andersson 等[31]匯報的II 型四氫呋喃的比容以及黃文件等[32]匯報的天然氣水合物的比容隨溫度的變化規律一致，而且相比于I 型水合物，II 型水合物的比熱較小，而且隨著晶穴占有率的下降而降低。因此，通過實驗室合成的水合物或自然界中的水合物由于晶格結構缺陷，可能會存在比熱低于理論值的情況。模擬得到的熱擴散率如表8所示，對于I 型CO2水合物，大小介于2.83～6.40×10-7m2·s-1范圍內，與姚貴策[25]測量的實驗值（4～6×10-7m2·s-1）相差不大，而且可以看出，II 型水合物的熱擴散率相對較大，尤其是晶格有缺陷的情況，熱擴散率更大。聲子平均自由程是由上述熱物性參數根據公式求解得到的，計算結果如表9所示，可以看出，II 型二氧化碳水合物的聲子平均自由程高于I 型二氧化碳水合物，這可能會影響到二者的分解行為，仍需要進一步研究。

表6 不同結構二氧化碳水合物的聲學速度（m/s）Table 6 Acoustic velocity of carbon dioxide hydrates with different structures

表7 不同結構二氧化碳水合物的比熱容（kJ/kg·K）Table 7 Specific heat capacity of carbon dioxide hydrates with different structures

表8 不同結構二氧化碳水合物的熱擴散率（×10-7 m2·s-1）Table 8 Thermal diffusivity of carbon dioxide hydrate with different structures

表9 不同結構二氧化碳水合物的聲子平均自由程（?）Table 9 Phonon mean free path of carbon dioxide hydrates with different structures

3 結論與展望

采用平衡分子動力學模擬手段，在20MPa、50-200K 條件下，測試了I 型CO2水合物、II 型CO2水合物（小晶穴占有率100%）、II 型CO2+CP二元水合物（大小晶穴占有率均為50%）和II 型CO2+CP 二元水合物（大小晶穴占有率均為100%）的熱物性，包括密度、等溫壓縮系數、聲學速度、比熱、熱導率、熱擴散率和聲子平均自由程。得到以下結論：（1）在環戊烷占據大晶穴條件下，II 型二氧化碳水合物的聲子輸運模式與I 型二氧化碳水合物不同，對于II 型水合物，隨著大晶穴的晶穴占有率的增加，低頻聲子導熱對水合物導熱起到的作用也變得顯著，高頻聲子和低頻聲子導熱均對水合物導熱起到重要作用。（2）當氣體水合物結構一定時，溫度越低，密度越大，等溫壓縮系數越低，聲學速度越高，比熱容越小。晶穴占有率均為100%時，與II 型二元水合物相比，在相同的溫度和壓力條件下，I 型CO2水合物密度更大，等溫壓縮系數更小，聲學速度更大，比熱容更大。密度和比熱容均隨著晶穴占有率的降低而降低。

II 型CO2水合物及多組晶穴占有率下II 型CO2+CP 二元水合物熱物性的測試和導熱機理的分析，為熱力學促進劑環戊烷作用下二氧化碳水合物在空調蓄冷等領域的應用提供了基礎數據，對于蓄冷量等參數的計算提供了依據。但是，本文對導熱機理部分的研究尚不夠深入。接下來，需要從聲子輸運的角度進一步研究不同結構下水合物的導熱機理。另外，量子效應對分子動力學模擬的結果也需要進一步探討。