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

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

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

?

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

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

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

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

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

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

1 模 擬

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

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

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

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

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

(1)

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

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

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

(2)

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

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

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

1.2 模擬方法及步驟

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

(3)

其中

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

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

(4)

其中

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

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

2 模擬結(jié)果分析

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

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

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

圖2 200 K下歸一化熱流自相關(guān)函數(shù)隨關(guān)聯(lián)時(shí)間的變化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

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

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

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

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

度。當(dāng)主客體分子間作用較弱時(shí),籠狀結(jié)構(gòu)的對(duì)稱性會(huì)約束甲烷分子在晶穴中的運(yùn)動(dòng),使其運(yùn)動(dòng)局域在晶穴中心,有著最小的勢(shì)能且相互間吸引力起主要作用，表明水分子被拉向晶穴的中心,導(dǎo)致晶格常數(shù)減小。隨著溫度和相互作用力強(qiáng)度的增加,客體分子的局域化特征減弱,更容易靠近晶穴附近的水分子,從而排斥水分子。甲烷分子在低溫和低頻區(qū)域與籠形水分子間的非簡(jiǎn)諧勢(shì)能更加顯著[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)表征不同頻率聲子振動(dòng)模式的強(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原子的振動(dòng)峰值分別在50～100 cm-1和300 cm-1附近,甲烷的振動(dòng)峰值與English等[8]的結(jié)果一致,分別對(duì)應(yīng)于水合物中大晶穴和小晶穴中甲烷分子的振動(dòng)。此外,甲烷分子的振動(dòng)在低頻區(qū)域局域化特征明顯,而在高頻區(qū)域局域化特征減弱。隨著甲烷分子與水分子間范德華作用強(qiáng)度的增加,受到主客體分子間共振散射的影響,甲烷分子的振動(dòng)峰值向高頻區(qū)域移動(dòng)。由圖5(b)可見(jiàn),不同主客體分子間范德華作用強(qiáng)度下O原子的VDOS變化不明顯,其兩個(gè)振動(dòng)峰值分別為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ν為每個(gè)聲子的能量,eV;1/(exp(hν/kBT)-1)為波色-愛(ài)因斯坦分布。

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

(2)高溫下的溫度相關(guān)性歸因于弛豫時(shí)間聲子的出現(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

國(guó)家自然科學(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)率分子動(dòng)力學(xué)模擬及分析[J]. 中國(guó)石油大學(xué)學(xué)報(bào)(自然科學(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.