劉庭崧， 劉 妮， 洪春芳

( 1. 上海理工大學(xué) 能源與動(dòng)力工程學(xué)院， 上海 200093; 2. 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室， 上海 200093)

1 引 言

氣體水合物是由客體分子和主體水分子形成的冰狀結(jié)晶化合物，其通常在低溫和高壓條件下形成[1, 2]. 目前，關(guān)于氣體水合物導(dǎo)熱性能的實(shí)驗(yàn)研究取得一系列成果. Rosenbaum等[3]采用單面瞬態(tài)平面熱源法,測(cè)得溫度為276K時(shí)CH4水合物導(dǎo)熱系數(shù)約為0.68 W·m-1·K-1. Huang等[4]采用瞬變平面熱源法測(cè)得CH4水合物在263.15-278.15K范圍內(nèi)的導(dǎo)熱系數(shù)約為0.56 W·m-1·K-1. 萬(wàn)麗華等[5]測(cè)得晶體形態(tài)和自保護(hù)作用下CO2水合物的導(dǎo)熱系數(shù)在溫度區(qū)間264.68-274.49 K增幅不大，而在275.47-282.04 K區(qū)間內(nèi)，導(dǎo)熱系數(shù)增長(zhǎng)明顯. 當(dāng)前關(guān)于水合物導(dǎo)熱系數(shù)的工作主要集中在CH4水合物，而對(duì)于CO2水合物導(dǎo)熱系數(shù)的分析相對(duì)較少. 為使CO2水合物更好應(yīng)用于工程實(shí)際，掌握其導(dǎo)熱系數(shù)隨各因素的變化規(guī)律非常重要. 分子動(dòng)力學(xué)(Molecular Dynamics，MD)為研究分子水平上的微觀機(jī)制提供了一套強(qiáng)有力的計(jì)算模擬方法. 基于此，MD方法成為研究水合物體系熱物理性質(zhì)的一種有效途徑. Jiang等[6]則通過(guò)NEMD方法比較了COS/G2和SPC/E兩種H2O分子模型對(duì)模擬結(jié)果的影響，模擬顯示根據(jù)COS/G2模型得出的水合物導(dǎo)熱系數(shù)更接近測(cè)試值. Wan等[7]指出晶體的低導(dǎo)熱特性主要?dú)w因于主體結(jié)構(gòu)，H2O分子對(duì)體系導(dǎo)熱的貢獻(xiàn)更大，主要因?yàn)槠溥\(yùn)動(dòng)引起的熱流更大. 籠形結(jié)構(gòu)中的客體分子也會(huì)使體系導(dǎo)熱進(jìn)一步增強(qiáng). 李期斌等[8]對(duì)高壓下CH4水合物的導(dǎo)熱機(jī)理做了相關(guān)分析，結(jié)果顯示，各種結(jié)構(gòu)水合物中H2O分子排布相似，但導(dǎo)熱系數(shù)存在一定差異. 對(duì)比功率圖譜對(duì)應(yīng)的最大聲子態(tài)密度得知，較高的壓力強(qiáng)化了CH4分子的擾動(dòng)，促進(jìn)水合物中的聲子傳熱.

MD模擬作為觀察微觀結(jié)構(gòu)的有力工具，從分子水平角度分析CO2水合物. 本文利用MD模擬在NVT系綜下進(jìn)行水合物的導(dǎo)熱模擬. 本次研究將采用改進(jìn)的Miiler-Plathe[9]模型，改進(jìn)算法模型結(jié)合了MD模擬熱導(dǎo)率需要的多種優(yōu)點(diǎn)：避免了出現(xiàn)人為的“溫度壁面”對(duì)系統(tǒng)構(gòu)型造成影響；X、Y、Z方向可進(jìn)行周期性擴(kuò)展，也適用于連續(xù)介質(zhì)；在整個(gè)模擬結(jié)構(gòu)上引起的溫度梯度相當(dāng)小，能同時(shí)保持能量和動(dòng)量守恒，保證所得導(dǎo)熱系數(shù)的合理性；局部區(qū)間收斂更快，很大程度上節(jié)約了模擬時(shí)間成本.

2 構(gòu)建模型與分子模擬

首先依據(jù)X射線單晶衍射數(shù)值確定水合物結(jié)構(gòu)H2O中O原子初始坐標(biāo)[10]，晶格中H原子排列符合Bernal-Fowler規(guī)則[11, 12]，CO2位于每個(gè)籠形中間. 然后在X、Y、Z三個(gè)方向擴(kuò)展成2×2×2超晶胞結(jié)構(gòu)，其晶胞參數(shù)為23.754 ? ×23.754 ? ×23.754 ?，如圖1所示. 此外，自然界中或?qū)嶒?yàn)室內(nèi)生成的水合物往往會(huì)存在結(jié)構(gòu)缺陷，還研究了不同種類結(jié)構(gòu)缺陷對(duì)水合物導(dǎo)熱特性的影響，基于完整的單晶胞CO2水合物結(jié)構(gòu)做出如表1所示的改變，形成結(jié)構(gòu)缺陷的CO2水合物，進(jìn)而以2×2×2的超晶胞結(jié)構(gòu)為模擬對(duì)象. 做出表中改變的原因在于研究晶穴占有率和CO2分子占據(jù)大籠或小籠對(duì)晶胞導(dǎo)熱特性的影響以及研究體系導(dǎo)熱的決定性因素在于主體還是客體.

圖1 完整CO2水合物超晶胞Fig. 1 Complete CO2 hydrate supercell

表1 缺陷結(jié)構(gòu)

模擬使用Materials Studio軟件，模擬體系X、Y、Z方向都設(shè)置為周期性邊界. 選用CVFF力場(chǎng)來(lái)描述分子間相互作用力，H2O和CO2分子勢(shì)能參數(shù)如表2所示. 首先，依據(jù)最速下降法和共軛梯度法依次對(duì)水合物結(jié)構(gòu)展開幾何優(yōu)化，使體系的初始結(jié)構(gòu)處于能量最小狀態(tài). 模擬過(guò)程中，系統(tǒng)溫度調(diào)節(jié)采用Berendsen算法[13]，Ewald加和用于計(jì)算長(zhǎng)程庫(kù)侖相互作用，范德華相互作用截?cái)喟霃饺?0 ?. 先使體系在NPT系綜下弛豫2 ps，達(dá)到設(shè)定平衡溫度，然后轉(zhuǎn)到NVT系綜持續(xù)計(jì)算100 ps. 所有的MD計(jì)算都使用Forcite軟件包進(jìn)行，時(shí)間步長(zhǎng)為1fs.

表2 相互作用參數(shù)

3 結(jié)果與討論

3.1 完整CO2水合物的導(dǎo)熱系數(shù)和功率圖譜分析

為驗(yàn)證模型的準(zhǔn)確性，首先模擬了268.15-278.15K范圍內(nèi)完整CO2水合物的導(dǎo)熱系數(shù)，結(jié)果與萬(wàn)麗華等[5]和本課題組[14]的實(shí)驗(yàn)測(cè)試值較為符合，如圖2所示. 圖中顯示模擬結(jié)果與兩組實(shí)驗(yàn)測(cè)量值均存在一定誤差，但與萬(wàn)麗華等的結(jié)果更為接近. 兩組實(shí)驗(yàn)數(shù)據(jù)相互存在差異，可能與晶體生成狀態(tài)、導(dǎo)熱測(cè)試方法等相關(guān).

圖2 純質(zhì)CO2水合物導(dǎo)熱性能Fig. 2 Thermal conductivity of pure CO2 hydrate

接著在溫度為200-280 K的NVT系綜下，對(duì)完整結(jié)構(gòu)的CO2水合物進(jìn)行100 ps的導(dǎo)熱計(jì)算和動(dòng)力學(xué)計(jì)算. 圖3為完整CO2水合物導(dǎo)熱能力與溫度的依賴關(guān)系，發(fā)現(xiàn)在200-230 K范圍內(nèi)，體系導(dǎo)熱系數(shù)隨溫度上升呈現(xiàn)弱線性增長(zhǎng)的趨勢(shì)，其值從0.4684 W·m-1·K-1變化到0.4836 W·m-1·K-1，平均增長(zhǎng)率為1.07%；溫度230-280 K范圍內(nèi)，其值從0.4836 W·m-1·K-1變化到0.7494 W·m-1·K-1，平均增長(zhǎng)率達(dá)9.25%，漲幅明顯變大. 上述現(xiàn)象表明在溫度較低時(shí)體系導(dǎo)熱的溫度相關(guān)性較弱，而溫度進(jìn)一步上升，此相關(guān)性變強(qiáng). 其主要原因?yàn)?，籠形結(jié)構(gòu)中包含的CO2分子隨溫度升高發(fā)生劇烈運(yùn)動(dòng)，使體系能量傳輸加快，主體H2O與客體CO2之間振動(dòng)耦合加強(qiáng)，為能量輸送提供了新的渠道，導(dǎo)致水合物導(dǎo)熱系數(shù)升高. 盡管水合物的導(dǎo)熱能力隨溫度上升會(huì)存在一定程度地提升，但其導(dǎo)熱系數(shù)始終處于一個(gè)較低的水平. 此現(xiàn)象可以解釋為，在水合物晶體結(jié)構(gòu)中，分子的平移運(yùn)動(dòng)和旋轉(zhuǎn)運(yùn)動(dòng)在很大程度上被限制，幾乎處在一個(gè)固定位置上振動(dòng)來(lái)傳遞能量. 此種傳遞能量的方式是非簡(jiǎn)諧的，使體系導(dǎo)熱系數(shù)整體偏低.

圖3 200-280 K溫度區(qū)間完整CO2水合物體系導(dǎo)熱模擬Fig. 3 Thermal conductivity simulations of complete CO2 hydrate system in the temperature range of 200-280 K

采用聲子態(tài)密度可表征不同頻率振動(dòng)模式的強(qiáng)弱特性，進(jìn)而了解晶體的傳熱過(guò)程. 聲子態(tài)密度可由粒子速度自相關(guān)函數(shù)經(jīng)傅里葉變換求出功率圖譜得到[15]，表示為功率圖譜曲線上的縱坐標(biāo). 在CO2水合物中客體CO2分子的功率圖譜表示為整個(gè)模擬系統(tǒng)聲子熱運(yùn)動(dòng)的平動(dòng)所做貢獻(xiàn)，主體H2O分子中O和H原子的功率圖譜則分別表示H2O的轉(zhuǎn)動(dòng)和平動(dòng)為聲子熱運(yùn)動(dòng)所做貢獻(xiàn)[8]. 圖4a為模擬溫度200-280 K下CO2分子的功率圖譜，橫坐標(biāo)表示經(jīng)傅里葉變化后的頻率(1cm-1=2.998×1010Hz)，縱坐標(biāo)代表歸一化的聲子態(tài)密度. 最大聲子態(tài)密度體現(xiàn)為功率圖譜上的峰值，對(duì)比該值對(duì)應(yīng)的橫坐標(biāo)頻率即可比較各參數(shù)對(duì)聲子運(yùn)動(dòng)的作用大小，其值對(duì)應(yīng)頻率越高，表明對(duì)體系導(dǎo)熱貢獻(xiàn)越大. 由圖4a可知，CO2分子平動(dòng)引起的最大聲子態(tài)密度對(duì)應(yīng)頻率主要位于頻率較低區(qū)域0-50 cm-1，說(shuō)明CO2分子對(duì)體系導(dǎo)熱特性的貢獻(xiàn)程度有限. 同時(shí)發(fā)現(xiàn)溫度上升時(shí)，最大聲子態(tài)密度所對(duì)應(yīng)的橫坐標(biāo)頻率略有增長(zhǎng)，表明溫度的上升有利于體系導(dǎo)熱性能的加強(qiáng)，但體現(xiàn)在對(duì)應(yīng)橫坐標(biāo)頻率上的增長(zhǎng)幅度有限. 圖4b, c表示模擬溫度200-280 K下H2O分子中O原子和H原子的功率圖譜，與CO2分子的功率圖譜相比，具有更多的高峰，圖形波動(dòng)更大. O原子和H原子的最大聲子態(tài)密度對(duì)應(yīng)頻率分別集中在170-250 cm-1和650-750 cm-1之間，其值遠(yuǎn)大于CO2平動(dòng)引起的最大聲子態(tài)密度所對(duì)應(yīng)的橫坐標(biāo)頻率，說(shuō)明水合物中主體結(jié)構(gòu)對(duì)體系導(dǎo)熱特性的貢獻(xiàn)較大. 顯然，在溫度升高時(shí)，O原子和H原子運(yùn)動(dòng)引起最大聲子態(tài)密度對(duì)應(yīng)頻率增長(zhǎng)幅度較CO2分子變化顯著，說(shuō)明溫度對(duì)主體H2O分子導(dǎo)熱性能作用較大，綜合表現(xiàn)為模擬溫度區(qū)間水合物的導(dǎo)熱能力隨溫度升高有所加強(qiáng).

圖 4 200-280 K溫度區(qū)間內(nèi)完整CO2水合物體系的運(yùn)動(dòng)功率圖譜. (a) CO2分子 (b) H2O分子中O原子 (c) H2O分子中H原子Fig. 4 Power spectra of complete CO2 hydrate system in the temperature range of 200-280 K. (a) CO2 molecules (b) O atom in H2O molecules (c) H atom in H2O molecules

3.2 結(jié)構(gòu)缺陷CO2水合物的導(dǎo)熱系數(shù)

模擬參數(shù)設(shè)置同完整CO2水合物一致，在溫度為200-280 K的NVT系綜下，對(duì)所有結(jié)構(gòu)缺陷的水合物進(jìn)行導(dǎo)熱計(jì)算和動(dòng)力學(xué)計(jì)算. 圖5是不同缺陷結(jié)構(gòu)水合物在溫度為200-280 K范圍內(nèi)導(dǎo)熱系數(shù)的變化曲線. 所有結(jié)構(gòu)中，完整晶胞結(jié)構(gòu)體系的導(dǎo)熱性能最好，各種水合物結(jié)構(gòu)的缺陷均使體系導(dǎo)熱性能有所削弱. 不同結(jié)構(gòu)晶胞導(dǎo)熱系數(shù)的大小順序基本表現(xiàn)為：λCO2> λDEL S1> λDEL L1≈ λDEL S2> λDEL L2> λEmpty> λDEL 3H2O> λDEL 6H2O，表明水合物晶穴占有率和籠形結(jié)構(gòu)缺陷對(duì)體系導(dǎo)熱系數(shù)產(chǎn)生某種程度的影響，并且后者的影響更大. 客體分子數(shù)量越多，一方面使體系的宏觀密度增加，從而導(dǎo)熱加強(qiáng)；另一方面，主客體分子相互作用加強(qiáng)，CO2分子與籠形結(jié)構(gòu)中H2O分子的碰撞激發(fā)了載熱的聲學(xué)聲子模式，其運(yùn)動(dòng)的不和諧性使有客體分子存在的籠形結(jié)構(gòu)比空籠結(jié)構(gòu)更疏松，出現(xiàn)能量的非局域化，兩者的綜合作用促進(jìn)水合物結(jié)構(gòu)傳熱過(guò)程的加強(qiáng). 但晶穴占有率的增加對(duì)體系導(dǎo)熱的影響是有限的，CO2分子數(shù)量增加，并不會(huì)使體系的導(dǎo)熱系數(shù)大幅度增加. 主體分子的籠形結(jié)構(gòu)對(duì)體系導(dǎo)熱能力影響顯著，對(duì)比空籠結(jié)構(gòu)和完整晶胞的導(dǎo)熱值，前者約為后者的86.67%，表明水合物的低導(dǎo)熱特性主要由籠形主體決定；客體分子的存在只能很小程度上增強(qiáng)體系的導(dǎo)熱能力，相較同溫度下完整CO2晶胞結(jié)構(gòu)，晶胞結(jié)構(gòu)SDEL S1和SDELS2導(dǎo)熱系數(shù)分別平均下降2.08%和6.22%，晶胞結(jié)構(gòu)SDEL L1和SDELL2導(dǎo)熱系數(shù)分別平均下降6.40%和7.50%，說(shuō)明體系導(dǎo)熱能力受整體晶穴占有率和客體分子占據(jù)籠形形態(tài)影響；當(dāng)籠形結(jié)構(gòu)被破壞時(shí)，體系導(dǎo)熱能力會(huì)被明顯削弱，晶胞結(jié)構(gòu)SDEL 3H2O和SDEL 6H2O導(dǎo)熱系數(shù)分別平均下降20.13%和24.03%，結(jié)構(gòu)破壞越嚴(yán)重，導(dǎo)熱系數(shù)越小，這可能是H2O分子的缺陷導(dǎo)致了導(dǎo)熱過(guò)程中聲子的大量散射，其產(chǎn)生的影響等效于聲子平均自由行程的減小，使水合物導(dǎo)熱減弱.

圖5 200-280 K溫度區(qū)間不同結(jié)構(gòu)缺陷CO2水合物導(dǎo)熱模擬Fig. 5 Thermal conductivity simulations of CO2 hydrate with different defect structures in the temperature range of 200-280 K

3.3 結(jié)構(gòu)缺陷CO2水合物聲子態(tài)密度對(duì)應(yīng)頻率

結(jié)構(gòu)缺陷CO2水合物中CO2分子、H2O分子中O原子和H原子功率圖譜與圖4a-c類似，由于功率圖譜中存在多個(gè)較高波峰，通過(guò)圖形難以辨認(rèn)最大波峰所對(duì)應(yīng)的橫坐標(biāo)頻率，因此將模擬溫度200-280 K下各結(jié)構(gòu)CO2分子、H2O分子中O原子和H原子的運(yùn)動(dòng)功率圖譜最大峰值所對(duì)應(yīng)的橫坐標(biāo)頻率進(jìn)行整合,如圖6a-c所示. 各結(jié)構(gòu)中CO2分子、H2O分子中O原子和H原子運(yùn)動(dòng)功率圖譜最大峰值所對(duì)應(yīng)的橫坐標(biāo)頻率大致滿足：SCO2> SDEL S1> SDEL L1> SDEL S2> SDEL L2> SEmpty> SDEL 3H2O> SDEL 6H2O，基本與各類型水合物導(dǎo)熱系數(shù)相對(duì)應(yīng). 當(dāng)水合物中存在CO2分子缺陷失時(shí)，最大聲子態(tài)密度對(duì)應(yīng)頻率略小于完整CO2水合物最大聲子態(tài)密度對(duì)應(yīng)頻率；而當(dāng)水合物中存在H2O分子缺失時(shí)，與完整CO2水合物相比對(duì)應(yīng)頻率則相差較大，表明主體H2O分子對(duì)體系導(dǎo)熱的貢獻(xiàn)較大. H2O分子缺失時(shí)運(yùn)動(dòng)功率圖譜峰值對(duì)應(yīng)的頻率較小，說(shuō)明H2O分子缺失將導(dǎo)致聲子的大量散射，使聲子的平均自由行程有所減小，弱化整個(gè)系統(tǒng)的導(dǎo)熱. 同樣，各結(jié)構(gòu)缺陷水合物的最大聲子態(tài)密度所對(duì)應(yīng)的頻率隨溫度變化趨勢(shì)與完整水合物的變化趨勢(shì)保持一致，都表現(xiàn)為隨溫度的上升而不斷增大，說(shuō)明水合物結(jié)構(gòu)的缺陷不影響體系導(dǎo)熱的溫度相關(guān)性.

圖6 200-280 K溫度區(qū)間內(nèi)最大聲子態(tài)密度所對(duì)應(yīng)頻率.(a) H2O中O原子 (b) H2O中H原子 (c) CO2分子Fig. 6 The frequencies corresponding to the maximum phonon state density in the temperature range of 200-280 K. (a) O atoms in H2O (b) H atoms in H2O (c) CO2 molecules

4 結(jié) 論

本文采用MD模擬方法對(duì)結(jié)構(gòu)完整CO2水合物、不同結(jié)構(gòu)缺陷CO2水合物在溫度為200-280 K的NVT系綜下進(jìn)行了導(dǎo)熱模擬計(jì)算并總結(jié)其變化規(guī)律，主要結(jié)論如下.

1. 對(duì)于結(jié)構(gòu)完整的CO2水合物，在溫度200-230 K時(shí)，體系導(dǎo)熱系數(shù)由0.4684 W·m-1·K-1變化到0.4836 W·m-1·K-1，溫度相關(guān)性較弱；而溫度230-280 K時(shí)，體系導(dǎo)熱系數(shù)由0.4836 W·m-1·K-1變化到0.7494 W·m-1·K-1，溫度相關(guān)性變強(qiáng). 此外，CO2分子和H2O分子中O、 H原子的功率圖譜峰值對(duì)應(yīng)頻率分別位于0-50 cm-1、170-250 cm-1和650-750 cm-1之間，主體分子對(duì)水合物體系的導(dǎo)熱貢獻(xiàn)更大.

2. 對(duì)于結(jié)構(gòu)缺陷的CO2水合物，各類缺陷結(jié)構(gòu)水合物導(dǎo)熱能力的大小順序基本表現(xiàn)為，λCO2> λDEL S1> λDEL L1≈ λDEL S2> λDEL L2> λEmpty> λDEL 3H2O> λDEL 3H2O，即晶穴占有率和籠形結(jié)構(gòu)缺陷對(duì)體系導(dǎo)熱均有一定影響，空籠晶胞導(dǎo)熱系數(shù)約為完整晶胞導(dǎo)熱系數(shù)的86.67%，體系的導(dǎo)熱能力主要取決于主體結(jié)構(gòu)的性質(zhì)，客體分子能加強(qiáng)體系的導(dǎo)熱，但影響程度有限.