CO2水合物導熱特性的分子動力學模擬

劉庭崧， 劉 妮， 洪春芳

( 1. 上海理工大學 能源與動力工程學院， 上海 200093; 2. 上海市動力工程多相流動與傳熱重點實驗室， 上海 200093)

1 引 言

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

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

2 構建模型與分子模擬

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

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

表1 缺陷結構

模擬使用Materials Studio軟件，模擬體系X、Y、Z方向都設置為周期性邊界. 選用CVFF力場來描述分子間相互作用力，H2O和CO2分子勢能參數如表2所示. 首先，依據最速下降法和共軛梯度法依次對水合物結構展開幾何優化，使體系的初始結構處于能量最小狀態. 模擬過程中，系統溫度調節采用Berendsen算法[13]，Ewald加和用于計算長程庫侖相互作用，范德華相互作用截斷半徑取為10 ?. 先使體系在NPT系綜下弛豫2 ps，達到設定平衡溫度，然后轉到NVT系綜持續計算100 ps. 所有的MD計算都使用Forcite軟件包進行，時間步長為1fs.

表2 相互作用參數

3 結果與討論

3.1 完整CO2水合物的導熱系數和功率圖譜分析

為驗證模型的準確性，首先模擬了268.15-278.15K范圍內完整CO2水合物的導熱系數，結果與萬麗華等[5]和本課題組[14]的實驗測試值較為符合，如圖2所示. 圖中顯示模擬結果與兩組實驗測量值均存在一定誤差，但與萬麗華等的結果更為接近. 兩組實驗數據相互存在差異，可能與晶體生成狀態、導熱測試方法等相關.

圖2 純質CO2水合物導熱性能Fig. 2 Thermal conductivity of pure CO2 hydrate

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

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

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

圖 4 200-280 K溫度區間內完整CO2水合物體系的運動功率圖譜. (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 結構缺陷CO2水合物的導熱系數

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

圖5 200-280 K溫度區間不同結構缺陷CO2水合物導熱模擬Fig. 5 Thermal conductivity simulations of CO2 hydrate with different defect structures in the temperature range of 200-280 K

3.3 結構缺陷CO2水合物聲子態密度對應頻率

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

圖6 200-280 K溫度區間內最大聲子態密度所對應頻率.(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 結 論

本文采用MD模擬方法對結構完整CO2水合物、不同結構缺陷CO2水合物在溫度為200-280 K的NVT系綜下進行了導熱模擬計算并總結其變化規律，主要結論如下.

1. 對于結構完整的CO2水合物，在溫度200-230 K時，體系導熱系數由0.4684 W·m-1·K-1變化到0.4836 W·m-1·K-1，溫度相關性較弱；而溫度230-280 K時，體系導熱系數由0.4836 W·m-1·K-1變化到0.7494 W·m-1·K-1，溫度相關性變強. 此外，CO2分子和H2O分子中O、 H原子的功率圖譜峰值對應頻率分別位于0-50 cm-1、170-250 cm-1和650-750 cm-1之間，主體分子對水合物體系的導熱貢獻更大.

2. 對于結構缺陷的CO2水合物，各類缺陷結構水合物導熱能力的大小順序基本表現為，λCO2> λDEL S1> λDEL L1≈ λDEL S2> λDEL L2> λEmpty> λDEL 3H2O> λDEL 3H2O，即晶穴占有率和籠形結構缺陷對體系導熱均有一定影響，空籠晶胞導熱系數約為完整晶胞導熱系數的86.67%，體系的導熱能力主要取決于主體結構的性質，客體分子能加強體系的導熱，但影響程度有限.