不同結(jié)構(gòu)二氧化碳水合物熱物性的分子動(dòng)力學(xué)模擬

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

（青島職業(yè)技術(shù)學(xué)院海爾學(xué)院 青島 266555）

0 引言

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

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

1 模型與方法

1.1 分子模型

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

表2 勢(shì)函數(shù)參數(shù)Table 2 Potential energy parameters

表1 模擬體系Table 1 Simulation system

圖1 初始構(gòu)型（藍(lán)線-氫鍵，棒狀-共價(jià)鍵，球形-原子，氧原子-紅色，氫原子-白色，碳原子-黃色）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 模擬方法

采用開(kāi)源軟件LAMMPS 進(jìn)行分子動(dòng)力學(xué)模擬。三個(gè)方向均采用周期性邊界條件，時(shí)間步長(zhǎng)為1fs，分子間相互作用力的截?cái)喟霃綖?2?，采用PPPM（particle-particle particle-mesh）方法計(jì)算長(zhǎng)程靜電力，精度為10-6，采用velocity-verlet 積分算法求解運(yùn)動(dòng)方程，采用shake 算法固定水分子的鍵長(zhǎng)和鍵角。首先進(jìn)行能量最小化，然后將體系置于NVT系綜下運(yùn)行500ps，采用Noose-Hoover 熱浴，阻尼系數(shù)設(shè)為0.1ps，使體系達(dá)到恒定溫度并保持平衡。接著將體系置于NPT 系綜下運(yùn)行500ps，采用Noose-Hoover 熱浴和壓浴，溫度阻尼系數(shù)設(shè)為0.1ps，壓力阻尼系數(shù)設(shè)為1ps，使體系達(dá)到設(shè)定溫度和壓力。最后，將體系置于NVE 系綜下運(yùn)行1ns。模擬結(jié)束后，采用平衡分子動(dòng)力學(xué)模擬（EMD）方法，運(yùn)用Green-Kubo 公式，編程計(jì)算熱導(dǎo)率，見(jiàn)公式（4）[20]。

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

聲子平均自由程是影響物體內(nèi)部熱輸運(yùn)的重要參數(shù)[21]，李佳等[22]通過(guò)對(duì)甲烷水合物的分子動(dòng)力學(xué)模擬發(fā)現(xiàn)，分解過(guò)程的固-液界面厚度與甲烷水合物內(nèi)部的聲子平均自由程十分接近，說(shuō)明水合物的導(dǎo)熱特性對(duì)水合物分解過(guò)程影響較大。根據(jù)分子動(dòng)力學(xué)理論，聲子平均自由程的計(jì)算方法如公式（5）所示[23]。

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

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

為了驗(yàn)證模擬方法和計(jì)算方法的正確性，測(cè)試計(jì)算了I 型CO2水合物的熱物性，在相似的熱力學(xué)條件下，模擬值與文獻(xiàn)值相接近，具體對(duì)比將在結(jié)果與討論部分展開(kāi)。

2 結(jié)果與討論

2.1 I 型CO2 水合物導(dǎo)熱特性

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

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

圖3 200K,20MPa 條件下（a）熱流自相關(guān)函數(shù)（b）熱導(dǎo)率隨關(guān)聯(lián)時(shí)間變化Fig.3 (a)HCACF and(b)thermal conductivity varies versus time under 200K,20MPa

圖4 是熱流自相關(guān)函數(shù)能譜，通過(guò)對(duì)熱流自相關(guān)函數(shù)作傅里葉變換得到，能夠反映出聲子能量的消散?？梢钥闯霾煌l率的聲子對(duì)水合物內(nèi)部熱輸運(yùn)的影響不同，高頻成分（＞15THz）大于低頻成分，說(shuō)明在水合物導(dǎo)熱過(guò)程中高頻聲子（15-25THz）起到主導(dǎo)作用，而低頻聲子（＜12THz）幅值相對(duì)較小，表明在低頻區(qū)域主客體分子之間充分耦合，而高頻區(qū)域的聲子更容易消散，即聲學(xué)聲子的導(dǎo)熱對(duì)水合物導(dǎo)熱起主導(dǎo)作用。I 型二氧化碳水合物的聲子作用模式和甲烷水合物相類(lèi)似[23]。

圖4 二氧化碳水合物熱流自相關(guān)函數(shù)能譜Fig.4 Power spectra of HCACF in CO2 hydrate

2.2 II 型CO2-CP 二元水合物導(dǎo)熱特性及對(duì)比分析

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

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

為了分析不同結(jié)構(gòu)二氧化碳水合物的導(dǎo)熱特性，圖6 繪制了20MPa、200K 條件下II 型CO2水合物（a）和II 型CO2-CP 水合物（晶穴占有率為50%（b）和100%（c））的熱流自相關(guān)函數(shù)能譜?？梢钥闯?，影響II 型水合物的聲子導(dǎo)熱模式除了高頻聲子外，低頻聲子區(qū)域（＜12THz）也存在明顯的峰值，而且當(dāng)二氧化碳僅占據(jù)小晶穴時(shí)，低頻聲子區(qū)域的能譜峰值較低，隨著環(huán)戊烷分子占據(jù)大晶穴，能譜峰值提高，當(dāng)二氧化碳和環(huán)戊烷完全占據(jù)大小晶穴時(shí)（晶穴占有率100%），低頻聲子區(qū)域的能譜峰值最高，幾乎接近于高頻聲子區(qū)域的能譜峰值。說(shuō)明對(duì)于II 型水合物，隨著大晶穴的晶穴占有率的增加，低頻聲子導(dǎo)熱對(duì)水合物導(dǎo)熱起到的作用也變得顯著，高頻聲子和低頻聲子導(dǎo)熱均對(duì)水合物導(dǎo)熱起到重要作用。

圖6 二元水合物熱流自相關(guān)函數(shù)能譜Fig.6 Power spectra of HCACF in CO2-CP hydrate

2.3 不同結(jié)構(gòu)二氧化碳水合物的熱物性

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

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

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

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

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

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

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

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

3 結(jié)論與展望

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

II 型CO2水合物及多組晶穴占有率下II 型CO2+CP 二元水合物熱物性的測(cè)試和導(dǎo)熱機(jī)理的分析，為熱力學(xué)促進(jìn)劑環(huán)戊烷作用下二氧化碳水合物在空調(diào)蓄冷等領(lǐng)域的應(yīng)用提供了基礎(chǔ)數(shù)據(jù)，對(duì)于蓄冷量等參數(shù)的計(jì)算提供了依據(jù)。但是，本文對(duì)導(dǎo)熱機(jī)理部分的研究尚不夠深入。接下來(lái)，需要從聲子輸運(yùn)的角度進(jìn)一步研究不同結(jié)構(gòu)下水合物的導(dǎo)熱機(jī)理。另外，量子效應(yīng)對(duì)分子動(dòng)力學(xué)模擬的結(jié)果也需要進(jìn)一步探討。