楊德偉， 劉雨文， 修 毓, 徐宏國， 苑昆鵬，徐 哲

(1.中國石油大學(xué)儲運(yùn)與建筑工程學(xué)院, 山東青島 266580; 2.中國石化勝利油田分公司濱南采油廠,山東濱州 256600；3.山東省青島市黃島區(qū)城市建設(shè)局工程建設(shè)管理中心,山東青島 266555)

?

甲烷水合物熱導(dǎo)率分子動力學(xué)模擬及分析

楊德偉1， 劉雨文2， 修 毓3, 徐宏國2， 苑昆鵬1，徐 哲1

(1.中國石油大學(xué)儲運(yùn)與建筑工程學(xué)院, 山東青島 266580; 2.中國石化勝利油田分公司濱南采油廠,山東濱州 256600；3.山東省青島市黃島區(qū)城市建設(shè)局工程建設(shè)管理中心,山東青島 266555)

采用平衡態(tài)分子動力學(xué)方法模擬分子甲烷水合物導(dǎo)熱性能,結(jié)合聲子態(tài)密度分析甲烷分子和水分子間的能量耦合過程;探究范德華相互作用對熱導(dǎo)率溫度相關(guān)性的影響。結(jié)果表明:熱導(dǎo)率隨著甲烷分子和水分子間范德華相互作用的增強(qiáng)而增大。相互作用的增強(qiáng)令甲烷分子的振動峰值向高頻區(qū)域移動,使得甲烷分子與水分子間的振動耦合作用增強(qiáng),VDOS匹配程度增加,進(jìn)而增大了甲烷水合物的熱導(dǎo)率。高溫下的溫度相關(guān)性歸因于弛豫時間聲子的出現(xiàn)導(dǎo)致的非彈性散射,低溫下主要受到光學(xué)聲子模式和低頻聲子的約束影響。模擬的熱導(dǎo)率的溫度依賴性與實(shí)驗(yàn)結(jié)果吻合較好。

甲烷水合物; 分子動力學(xué); 聲子; 熱導(dǎo)率

作為新能源的甲烷水合物是保證國家能源安全和保護(hù)環(huán)境而迫切需要解決的問題之一[1]。甲烷水合物不同于冰及其他分子晶體的是其熱導(dǎo)率在接近德拜溫度的高溫區(qū)域表現(xiàn)出類似非晶體的特性,而在低溫下呈現(xiàn)出晶體特性[2]。Jiang等[3-4]采用非平衡態(tài)分子動力學(xué)方法研究了溫度在30～260 K甲烷、二氧化碳和xenon水合物的熱導(dǎo)率。水合物熱導(dǎo)率的溫度相關(guān)性在溫度低于50 K時與實(shí)驗(yàn)結(jié)果一致。English等[5-9]采用Green-Kubo方法研究了不同類型的甲烷水合物的熱導(dǎo)率,并分析了長程靜電力的處理對結(jié)果的影響。研究發(fā)現(xiàn)sI型甲烷水合物熱導(dǎo)率變化的轉(zhuǎn)變溫度為90 K。甲烷水合物的熱導(dǎo)率不僅與籠形結(jié)構(gòu)的穩(wěn)定性有關(guān),還取決于甲烷分子和水分子間的相互作用。平衡態(tài)分子動力學(xué)方法(EMD)和非平衡態(tài)方法都是利用Green-Kubo 關(guān)系式積分熱流自相關(guān)函數(shù)計算導(dǎo)熱系數(shù),該方法能夠從平衡系統(tǒng)微小的溫度波動中將導(dǎo)熱系數(shù)提取出來。相比非平衡態(tài)方法,該方法由于使整個體系保持在同樣的溫度下,溫度波動帶來的誤差相對小很多。不同的是平衡態(tài)方法在NVE系綜(恒N、V、E系綜)完成導(dǎo)熱系數(shù)的計算,能夠模擬真實(shí)的物理情景。筆者采用平衡態(tài)分子動力學(xué)方法模擬分子甲烷水合物導(dǎo)熱性能,結(jié)合聲子態(tài)密度分析甲烷分子和水分子間的能量耦合過程;探究范德華相互作用對熱導(dǎo)率溫度相關(guān)性的影響。

1 模 擬

1.1 物理模型及勢函數(shù)

根據(jù)X射線衍射實(shí)驗(yàn)[10]結(jié)果構(gòu)建1×1×1的sI型甲烷水合物晶胞,然后向X、Y、Z三個方向拓展成2×2×2的sI型甲烷水合物晶胞,如圖1所示。模擬中3個方向上都采用周期性邊界條件。由于甲烷水合物的聲子平均自由程小于一個原胞的長度,因此大尺寸系統(tǒng)不會顯現(xiàn)出顯著的尺寸效應(yīng)。

采用TIP4P[11]勢能函數(shù)描述水分子之間的相互作用,甲烷采用OPLS[12]模型。目前常用的幾種描述甲烷分子間的勢能模型為:Tersoff勢能、Brenner勢能及Brenner勢能基礎(chǔ)上修正后的REBO勢能等。Tersoff勢能能夠準(zhǔn)確地模擬碳單鍵、雙鍵、三鍵勢能與鍵長的關(guān)系,但該勢能在描述一些特殊成鍵(如空位)時與實(shí)際物理情況不符。后面兩種勢能均在對Tersoff勢能修正的基礎(chǔ)上發(fā)展起來,且第三種勢能在第二種勢能的基礎(chǔ)上又做了進(jìn)一步的修

正,但是這兩種勢能均須做大量的插值計算,長時間計算會降低結(jié)果的精度。通過對這三種勢能的實(shí)際計算比較,發(fā)現(xiàn)在不考慮特殊成鍵的情況下,Tersoff勢能在長時間的NVE循環(huán)中能更好地保持溫度的恒定性,溫度基本上在平衡值上下波動。故采用Tersoff型勢函數(shù)模擬甲烷分子之間相互作用。模擬過程不考慮原子間長程庫侖力的作用。Tersoff勢能的表達(dá)式為

Vtersoff(rij)=fC(rij)[fR(rij)+bijfA(rij)].

(1)

式中,i,j為體系中的原子;rij為i,j原子成鍵鍵長;Vtersoff(rij)為i原子和其鄰近原子j之間的作用能;fR、fA分別為產(chǎn)生作用能的排斥項(xiàng)和吸引項(xiàng);fC為勢能截取函數(shù)為bij為鍵序函數(shù)。

圖1 2×2×2的sI型甲烷水合物晶胞Fig.1 Physical model of sI methane hydrate with size 2×2×2

Tersoff多體勢能只考慮了近程原子的作用力,即當(dāng)原子間距大于S=2.1 ?時,原子間的多體相互作用勢能被截斷。故在計算中,有時還要考慮遠(yuǎn)程力的作用,比如范德華力的作用。范德華作用勢能,即Lennard-Jones勢能的表達(dá)式為

(2)

式中,ε和σ分別為能量和距離參數(shù);χ為標(biāo)度系數(shù),用來調(diào)整范德華作用力強(qiáng)弱。

不同種類的分子之間的勢能參數(shù)可以采用Lorentz-Berthelot混合規(guī)則計算得出。本文中采用的分子間的L-J勢能參數(shù)見表1。采用精度為1×10-4的PPPM方法計算長程靜電力,截斷半徑為8.5 ?。采用Lennard-Jones(L-J)勢能模型描述甲烷分子和水分子之間的相互作用,截斷半徑為10.0 ?。由于Tersoff勢能在r為0.21 nm時已經(jīng)為零,而Lennard-Jones勢能在該處的值不為零,為了使勢能曲線在整個空間光滑,需要在其銜接處做樣條函數(shù)插值處理。

表1 分子間L-J勢能參數(shù)

1.2 模擬方法及步驟

采用平衡態(tài)分子動力學(xué)方法(EMD)計算甲烷水合物熱導(dǎo)率。EMD方法是利用基于線性響應(yīng)理論的Green-Kubo公式計算熱導(dǎo)率,其表達(dá)式為

(3)

其中

式中,λ為熱導(dǎo)率,W/(m·k);V為系統(tǒng)的體積,m3;T為模擬溫度,K;kB為Boltzmann常數(shù);Jq(τ)為τ時刻的總熱流,W·m;〈Jq(τ)·Jq(0)〉為熱流自相關(guān)函數(shù);εi為i原子的總能量,J。

采用Green-Kubo計算式計算導(dǎo)熱系數(shù),表示為

(4)

其中

式中，Vi為i原子的速度。

考慮到并行計算和采用的勢函數(shù),采用分子動力學(xué)程序LAMMPS進(jìn)行模擬[13]。首先采用共軛梯度算法對晶胞進(jìn)行能量最小化,然后采用Nosé-Hoover熱浴方法[14-15]使系統(tǒng)在NVT系綜下弛豫100萬步,調(diào)節(jié)系統(tǒng)到目標(biāo)溫度,之后轉(zhuǎn)到NPT系綜繼續(xù)弛豫100萬步,使系統(tǒng)達(dá)到平衡,然后再持續(xù)運(yùn)行400萬步統(tǒng)計熱流并進(jìn)行熱導(dǎo)率的計算,以確保熱流自相關(guān)函數(shù)收斂至零。模擬采用的時間步長為1 fs,運(yùn)動方程的積分選用Velocity-Verlet算法。

2 模擬結(jié)果分析

2.1 熱導(dǎo)率的溫度相關(guān)性

圖2給出了溫度為200 K時熱流自相關(guān)函數(shù)隨時間的收斂情況。從圖2可以看出,熱流自相關(guān)函數(shù)在開始時迅速衰減并在零值附近振蕩,在0.5 ps后逐漸收斂于零,表明在模擬時間內(nèi)水合物的熱導(dǎo)率已經(jīng)收斂。

圖3給出了客體甲烷分子100%填充和50%填充的水合物的熱導(dǎo)率,同時給出了實(shí)驗(yàn)上采用瞬態(tài)面熱源法的測量結(jié)果[5]。

圖2 200 K下歸一化熱流自相關(guān)函數(shù)隨關(guān)聯(lián)時間的變化Fig.2 Change of normalized heat current autocorrelation function with correlation time at temperature of 200 K

圖3 0.1 MPa不同溫度下甲烷水合物的熱導(dǎo)率Fig.3 Computed thermal conductivities for methane hydrate at different temperature for pressure of 0.1 MPa

計算得到的100%填充的水合物的熱導(dǎo)率在265 K時為0.75 W·m-1·K-1,比實(shí)驗(yàn)數(shù)據(jù)0.61 W·m-1·K-1高了10%,主要原因是合成的水合物試樣存在孔隙,且在模擬過程中存在晶格缺陷等。當(dāng)溫度從265 K減小到100 K時,甲烷水合物的熱導(dǎo)率隨著溫度減小,這是由于聲子非彈性散射的減小導(dǎo)致能量傳遞通道減少,進(jìn)而使得熱導(dǎo)率減小。一個顯著的特征是隨著溫度繼續(xù)減小熱導(dǎo)率開始增加。類似的行為也體現(xiàn)在實(shí)驗(yàn)結(jié)果中,但其幅值更小。然而,事實(shí)上很難直接和定量地比較分子動力學(xué)模擬結(jié)果和實(shí)驗(yàn)結(jié)果,因?yàn)閷?shí)驗(yàn)結(jié)果與樣品的特性和品質(zhì)有關(guān),例如孔隙度。此外,除了系統(tǒng)尺寸和靜電力之外,勢能函數(shù)的準(zhǔn)確性也給比較帶來了難度。在整個考慮的溫度范圍,100%填充率的甲烷水合物熱導(dǎo)率比50%填充率的水合物的熱導(dǎo)率大5%。晶格籠中包含甲烷分子導(dǎo)致了能量傳遞的增強(qiáng),表明客體分子和主體分子間發(fā)生了更多的振動耦合,盡管這種趨勢在模擬的誤差范圍內(nèi)。計算的甲烷水合物的熱導(dǎo)率在模擬的溫度區(qū)間內(nèi)對晶穴占有率不敏感。之前的分子動力學(xué)(MD)模擬也表明sI甲烷水合物的熱導(dǎo)率與晶穴占有率的相關(guān)性很微弱。100%占有率的甲烷水合物的熱導(dǎo)率行為主要?dú)w因于聲子輸運(yùn)中的短程弛豫時間,即更短的聲子平均自由程。這個發(fā)現(xiàn)意味著甲烷水合物和冰之間的晶格結(jié)構(gòu)的不同起著更重要的作用,這與之前對甲烷水合物低熱導(dǎo)率的解釋不同。換句話說,甲烷分子與水籠子的碰撞可能激發(fā)了載熱的聲學(xué)聲子模式,水分子的氫鍵網(wǎng)絡(luò)比空籠水合物更疏松,導(dǎo)致新的聲學(xué)模式的出現(xiàn)和能量的非局域化,導(dǎo)致熱導(dǎo)率的升高。這種特殊的熱行為可能是由于客體分子的相互作用引起的非簡諧聲學(xué)模式聲子的振動。所以,熱導(dǎo)率與空穴填充率的關(guān)系可以歸因于非簡諧聲子運(yùn)動的聲子非彈性散射。

2.2 主客體分子間作用力強(qiáng)度的影響

為了研究客體甲烷分子和主體水分子間范德華作用強(qiáng)度對熱導(dǎo)率的影響,分別計算了χ=1、2、3、4時甲烷水合物的熱導(dǎo)率。溫度為200 K、不同范德華作用強(qiáng)度下甲烷水合物的熱導(dǎo)率見圖4。

由圖4可知,甲烷水合物的熱導(dǎo)率隨著客體分子和主體分子間范德華作用強(qiáng)度的增強(qiáng)而增加,表明熱導(dǎo)率正比于甲烷分子和水分子間的相互作用強(qiáng)

度。當(dāng)主客體分子間作用較弱時,籠狀結(jié)構(gòu)的對稱性會約束甲烷分子在晶穴中的運(yùn)動,使其運(yùn)動局域在晶穴中心,有著最小的勢能且相互間吸引力起主要作用，表明水分子被拉向晶穴的中心,導(dǎo)致晶格常數(shù)減小。隨著溫度和相互作用力強(qiáng)度的增加,客體分子的局域化特征減弱,更容易靠近晶穴附近的水分子,從而排斥水分子。甲烷分子在低溫和低頻區(qū)域與籠形水分子間的非簡諧勢能更加顯著[16]。

圖4 200 K,0.1 MPa熱導(dǎo)率隨相互作用力強(qiáng)度χ的變化Fig.4 Computed thermal conductivities of methane hydrate for different values of χ at T=200 K,p=0.1 MPa

采用聲子態(tài)密度(VDOS)表征不同頻率聲子振動模式的強(qiáng)弱,可由速度自相關(guān)函數(shù)(VACF)的傅立葉變換得到。圖5為不同范德華作用強(qiáng)度下C原子和O原子的聲子態(tài)密度。

圖5 不同范德華作用強(qiáng)度下C原子和O原子的聲子態(tài)密度Fig.5 VDOS of carbon and oxygen atoms at different interaction strengths

由圖5可知, C原子的振動峰值分別在50～100 cm-1和300 cm-1附近,甲烷的振動峰值與English等[8]的結(jié)果一致,分別對應(yīng)于水合物中大晶穴和小晶穴中甲烷分子的振動。此外,甲烷分子的振動在低頻區(qū)域局域化特征明顯,而在高頻區(qū)域局域化特征減弱。隨著甲烷分子與水分子間范德華作用強(qiáng)度的增加,受到主客體分子間共振散射的影響,甲烷分子的振動峰值向高頻區(qū)域移動。由圖5(b)可見,不同主客體分子間范德華作用強(qiáng)度下O原子的VDOS變化不明顯,其兩個振動峰值分別為60 cm-1和300 cm-1。衡量C原子和O原子VDOS的重疊程度[17]的表達(dá)式為

(5)

式中,E為VDOS重疊區(qū)域的聲子能量,eV;h為普朗克常數(shù),eV·s;ν為聲子頻率,Hz;kB為玻爾茲曼常數(shù),eV/K;T為溫度,K;go(v)為原子重疊的VDOS;hν為每個聲子的能量,eV;1/(exp(hν/kBT)-1)為波色-愛因斯坦分布。

圖6為χ下C原子和O原子VDOS重疊區(qū)域的能量。由圖6可知,隨著主客體分子間范德華作用的增強(qiáng),C原子和O原子VDOS重疊區(qū)域的能量增加。這是因?yàn)橹骺腕w分子間范德華作用的增強(qiáng)令甲烷分子的振動峰值向高頻區(qū)域移動,使得甲烷分子與水分子間的振動耦合作用增強(qiáng),VDOS匹配程度增加,從而提高了甲烷水合物的熱導(dǎo)率。

圖6 不同χ下C原子和O原子VDOS重疊區(qū)域的能量Fig.6 Energy in VDOS overlapped region of carbon and oxygen atoms for different values of χ

3 結(jié) 論

(1)熱導(dǎo)率隨著甲烷分子和水分子間范德華相互作用的增強(qiáng)而增大。相互作用的增強(qiáng)令甲烷分子的振動峰值向高頻區(qū)域移動,使得甲烷分子與水分子間的振動耦合作用增強(qiáng),VDOS匹配程度增加,進(jìn)而增大了甲烷水合物的熱導(dǎo)率。

(2)高溫下的溫度相關(guān)性歸因于弛豫時間聲子的出現(xiàn)導(dǎo)致的非彈性散射,低溫下主要受到光學(xué)聲子模式和低頻聲子的約束影響。

[1] SLOAN E D. Fundamental principles and applications of natural gas hydrates[J]. Nature, 2003,426(6964):353-363.

[2] TSE J S, WHITE M A. Origin of glassy crystalline behavior in the thermal properties of clathrate hydrates: a thermal conductivity study of tetrahydrofuran hydrate[J]. The Journal of Physical Chemistry, 1988,92(17):5006-5011.

[3] JIANG H, MYSHAKIN E M, JORDAN K D, et al. Molecular dynamics simulations of the thermal conductivity of methane hydrate[J]. The Journal of Physical Chemistry B, 2008,112(33):10207-10216.

[4] JIANG H, JORDAN K D. Comparison of the properties of xenon, methane, and carbon dioxide hydrates from equilibrium and nonequilibrium molecular dynamics simulations[J]. The Journal of Physical Chemistry C, 2009,114(12):5555-5564.

[5] ROSENBAUM E J, ENGLISH N J, JOHNSON J K, et al. Thermal conductivity of methane hydrate from experiment and molecular simulation[J]. The Journal of Physical Chemistry B, 2007,111(46):13194-13205.

[6] ENGLISH N J. Effect of electrostatics techniques on the estimation of thermal conductivity via equilibrium molecular dynamics simulation: application to methane hydrate[J]. Molecular Physics, 2008,106(15):1887-1898.

[7] ENGLISH N J, JOHN S T. Mechanisms for thermal conduction in methane hydrate[J]. Physical Review Letters, 2009,103(1):015901.

[8] ENGLISH N J, JOHN S T, CAREY D J. Mechanisms for thermal conduction in various polymorphs of methane hydrate[J]. Physical Review B, 2009,80(13):134306.

[9] ENGLISH N J, JOHN S T. Guest and host contributions towards thermal conduction in various polymorphs of methane hydrate[J]. Computational Materials Science, 2010,49(4):S176-S180.

[10] KIRCHNER M T, BOESE R, BILLUPS W E, et al. Gas hydrate single-crystal structure analyses[J]. Journal of the American Chemical Society, 2004,126(30):9407-9412.

[11] ABASCAL J L F, VEGA C. A general purpose model for the condensed phases of water: TIP4P/2005[J]. The Journal of Chemical Physics, 2005,123(23):234505.

[12] JORGENSEN W L, MADURA J D, SWENSON C J. Optimized intermolecular potential functions for liquid hydrocarbons[J]. Journal of the American Chemical Society, 1984,106(22):6638-6646.

[13] PLIMPTON S. Fast parallel algorithms for short-range molecular dynamics[J]. Journal of Computational Physics, 1995,117(1):1-19.

[14] NOSE S. A unified formulation of the constant temperature molecular dynamics methods[J]. The Journal of Chemical Physics, 1984,81(1):511-519.

[15] HOOVER W G. Canonical dynamics: equilibrium phase-space distributions[J]. Physical Review A, 1985,31(3):1695.[16] SCHOBER H, ITOH H, KLAPPROTH A, et al. Guest-host coupling and anharmonicity in clathrate hydrates[J]. The European Physical Journal E: Soft Matter and Biological Physics, 2003,12(1):41-49.

[17] ZHANG X, HU M, TANG D. Thermal rectification at silicon/horizontally aligned carbon nanotube interfaces[J]. Journal of Applied Physics, 2013,113(19):194307.

(編輯 沈玉英)

Molecular dynamics simulation and analysis of thermal conductivity of methane hydrate

YANG Dewei1, LIU Yuwen2，XIU Yu3, XU Hongguo2，YUAN Kunpeng1，XU Zhe1

(1.CollegeofPipelineandCivilEngineeringinChinaUniversityofPetroleum，Qingdao266580,China； 2.BinnanOilProductionPlantofShengliOilfieldBranch,SINOPEC，Binzhou256600,China； 3.ProjectConstructionManagementCenterofCityConstructionBureauofHuangdaoDistrictinQingdao,ShandongProvince,Qingdao266555,China)

The heat conduction of methane hydrate was simulated using equilibrium molecular dynamics, and the thermal coupling between methane molecules and water lattices was studied by combining with the analysis of phonon density of states (VDOS). The influence of Van der Waals interaction strength on the temperature dependence of thermal conductivity was also investigated. The results show that the thermal conductivity increases proportionally with the enhancement of the Van der Waals interaction strength. With the increase of the interaction strength, the vibration peak of methane molecules shifts to a higher frequency region because of the stronger vibration coupling and the better matching of VDOS between methane molecules and water lattices, and then the thermal conductivity of methane hydrate is enhanced. The temperature dependence at high temperature may be attributed to inelastic scattering of the phonon caused by the appearance of phonons with the immediate relaxation time, while at low temperature it may be attributed to the confinement of the optic phonon modes and low frequency phonons. The calculated temperature dependence trend agrees well with the experimental results.

methane hydrate; molecular dynamics; phonon; thermal conductivity

2015-06-10

國家自然科學(xué)基金項(xiàng)目 (U1262112)

楊德偉(1964-),男,教授,博士,研究方向?yàn)闊崮芄こ碳夹g(shù)。E-mail:yangdw@upc.edu.cn。

1673-5005(2016)04-0141-05

10.3969/j.issn.1673-5005.2016.04.019

TK 124

A

楊德偉，劉雨文，修毓,等.甲烷水合物熱導(dǎo)率分子動力學(xué)模擬及分析[J]. 中國石油大學(xué)學(xué)報(自然科學(xué)版),2016,40(4):141-145.

YANG Dewei, LIU Yuwen, XIU Yu, et al. Molecular dynamics simulation and analysis of thermal conductivity of methane hydrate[J]. Journal of China University of Petroleum (Edition of Natural Science), 2016,40(4):141-145.