劉明，徐哲

（1 勝利油田分公司石油工程技術(shù)研究院，山東省稠油開(kāi)采技術(shù)重點(diǎn)實(shí)驗(yàn)室，山東東營(yíng)257000；2中國(guó)石油大學(xué)（華東）新能源學(xué)院，山東青島266580）

引 言

甲烷水合物是由甲烷和水形成的一種籠型結(jié)構(gòu)的化合物[1-3]，其熱導(dǎo)率在甲烷水合物勘探、開(kāi)采、儲(chǔ)運(yùn)以及其他應(yīng)用過(guò)程中起著重要的作用，因此研究甲烷水合物的導(dǎo)熱具有重要意義[4]。

近年來(lái)，許多學(xué)者利用分子動(dòng)力學(xué)模擬方法(MD)研究了甲烷水合物的導(dǎo)熱特性[5-9]。Rosenbaum等[10]采用EMD 研究了水分子模型采用SPC/E、TIP4P-Ew、TIP4P-FQ 時(shí)甲烷水合物的熱導(dǎo)率，討論了水分子勢(shì)能函數(shù)的選擇對(duì)模擬結(jié)果的影響，他們認(rèn)為極化的水分子模型得到的結(jié)果與實(shí)驗(yàn)值更吻合。然而，當(dāng)選用極化的水分子模型時(shí)，存在計(jì)算量大、耗時(shí)多等缺點(diǎn)，因此結(jié)構(gòu)簡(jiǎn)單、計(jì)算效率高的非極化水分子模型被廣泛應(yīng)用于MD 模擬。Jiang等[11-12]采用非平衡分子動(dòng)力學(xué)模擬方法研究了當(dāng)溫度位于30～260 K 范圍內(nèi)時(shí)甲烷水合物的熱導(dǎo)率，模擬中水分子采用COS/G2和SPC/E 模型，研究發(fā)現(xiàn)當(dāng)溫度低于50 K 時(shí)體系的溫度變化趨勢(shì)與實(shí)驗(yàn)結(jié)果吻合，而高溫下模擬結(jié)果與實(shí)驗(yàn)值的變化趨勢(shì)相反。Krivchikov 等[13]把Xe 水合物和甲烷水合物熱導(dǎo)率的溫度特性分為了四個(gè)不同的區(qū)域，其中區(qū)域Ⅰ(2～6 K)內(nèi)熱導(dǎo)率與溫度的關(guān)系可以表示為k∝T1.2～1.4，區(qū)域Ⅱ(6～54 K)內(nèi)k∝T0.3，區(qū)域Ⅰ和區(qū)域Ⅱ的熱導(dǎo)率可用軟體勢(shì)模型解釋，在區(qū)域Ⅲ(54～94 K)，由于共振散射效應(yīng)，水合物的熱導(dǎo)率逐漸減小，區(qū)域Ⅳ為溫度大于94 K 的區(qū)域，該區(qū)域聲子的平均自由程達(dá)到最小值，熱導(dǎo)率與MD 模擬結(jié)果相吻合，因此，目前還沒(méi)有一種模型可以定量地描述水合物在整個(gè)溫度范圍內(nèi)的導(dǎo)熱情況。此外，在水合物的熱傳遞機(jī)理探究上，不同學(xué)者之間依然存在很大分歧。Tse等[14]針對(duì)甲烷水合物的熱輸運(yùn)提出了共振散射模型，認(rèn)為能量輸運(yùn)過(guò)程出現(xiàn)的熱耗散是由晶格聲子的振動(dòng)引起的能量轉(zhuǎn)移造成的，然而相關(guān)的實(shí)驗(yàn)研究表明[15]，方鈷晶體的熱輸運(yùn)機(jī)制與共振散射模型機(jī)制相互矛盾。Jiang 等[11-12]研究發(fā)現(xiàn)水合物的低導(dǎo)熱特性主要是由水合物中復(fù)雜的籠形結(jié)構(gòu)造成的。而English 等[16]采用平衡分子動(dòng)力學(xué)模擬方法研究了甲烷水合物的熱導(dǎo)率，他們認(rèn)為sI 甲烷水合物的低熱導(dǎo)率除了與水分子形成的晶格結(jié)構(gòu)有關(guān)外，還受主客體分子相互作用的制約，甲烷水合物在德拜溫度附近呈現(xiàn)的類玻璃的溫度依賴性主要受客體分子及主客體分子相互作用的控制。

綜上所述，為了合理描述甲烷水合物在整個(gè)溫度范圍內(nèi)的導(dǎo)熱情況，本文采用分子動(dòng)力學(xué)模擬計(jì)算得出甲烷水合物的熱導(dǎo)率，并在一定溫度范圍內(nèi)對(duì)其進(jìn)行量子修正。通過(guò)對(duì)甲烷水合物的導(dǎo)熱進(jìn)行聲子模式分解探究了甲烷水合物的導(dǎo)熱機(jī)理，并研究了主客體分子相互作用對(duì)甲烷水合物導(dǎo)熱的影響。

1 計(jì)算方法

1.1 模擬方法

MD 模擬計(jì)算熱導(dǎo)率有兩種方法：其中平衡方法依據(jù)Green-Kubo(G-K)線性響應(yīng)理論，通過(guò)對(duì)熱流自相關(guān)函數(shù)積分計(jì)算熱導(dǎo)率；非平衡方法通過(guò)在系統(tǒng)內(nèi)構(gòu)造熱流，根據(jù)傅里葉導(dǎo)熱定律來(lái)計(jì)算熱導(dǎo)率。本文借助LAMMPS 軟件[17]，采用平衡方法來(lái)進(jìn)行模擬。在計(jì)算中水分子采用TIP4P/2005 模型，H—O 鍵的鍵長(zhǎng)固定為0.9572 ?（1 ? =0.1 nm），H—O—H 的夾角固定為104.52°。甲烷分子采用OPLSUA 聯(lián)合原子模型，不同原子間的相互作用可以表示為

式中，S 為調(diào)整范德華力作用強(qiáng)弱的標(biāo)度系數(shù)；εij和σij分別為勢(shì)阱深度和原子間平衡距離；qi和qj為原子所帶電荷量；rij為原子間距離；ε0為電常數(shù)。模擬中設(shè)置范德華力截?cái)喟霃綖?0 ?，Bugel 等[18]研究認(rèn)為截?cái)喟霃綖?.5σ 時(shí)，截?cái)喈a(chǎn)生的誤差就可以忽略。采用pppm 算法處理長(zhǎng)程靜電力相互作用，Luty等[19]認(rèn)為pppm 算法在保證精度的條件下可以提高計(jì)算效率。選用velocity-Verlet 積分算法求解運(yùn)動(dòng)方程。

熱導(dǎo)率的計(jì)算使用G-K關(guān)系式[20]

式中，kB為玻耳茲曼常數(shù)；T 為模擬系統(tǒng)的溫度；V為系統(tǒng)體積；J為系統(tǒng)的微觀熱流，J(t)· J(0)稱為熱流自相關(guān)函數(shù)(heat current autocorrelation function，HCACF)。

1.2 模擬方法

平衡分子動(dòng)力學(xué)模擬是在2 個(gè)×2 個(gè)×2 個(gè)sⅠ型甲烷水合物晶胞內(nèi)進(jìn)行的，其結(jié)構(gòu)如圖1 所示。晶胞的長(zhǎng)度為11.877 ?。

圖1 甲烷水合物模擬結(jié)構(gòu)Fig.1 Initial configuration of methane hydrate

MD 模擬分為三步：首先將整體放置于NVT 系綜弛豫1 ns，采用Nose-Hoover 熱浴法[21-22]使系統(tǒng)達(dá)到預(yù)設(shè)的目標(biāo)溫度，時(shí)間步長(zhǎng)設(shè)置為1 fs。然后將系統(tǒng)轉(zhuǎn)入NPT 系綜，設(shè)置目標(biāo)壓力為10 MPa，在目標(biāo)溫度下弛豫1 ns，采用Nose-Hoover熱浴及壓浴法使系統(tǒng)達(dá)到平衡，最后將系統(tǒng)轉(zhuǎn)入NVE 系綜運(yùn)行2 ns，計(jì)算熱導(dǎo)率。

2 結(jié)果與討論

2.1 量子修正

當(dāng)模擬溫度遠(yuǎn)低于德拜溫度時(shí)，大量聲子模式?jīng)]有被完全激發(fā)出來(lái)，此時(shí)的量子效應(yīng)不可忽略，經(jīng)典的MD 方法無(wú)法求出精確的熱導(dǎo)率值，可以考慮引入量子修正解決這個(gè)問(wèn)題。量子修正的主要思想是通過(guò)經(jīng)典力學(xué)系統(tǒng)中的模擬溫度來(lái)計(jì)算量子系統(tǒng)中對(duì)應(yīng)的修正溫度。Lukes 等[23]應(yīng)用量子修正方法準(zhǔn)確計(jì)算出了單壁碳納米管的熱導(dǎo)率。量子修正假設(shè)模擬溫度下的聲子總能量等于量子修正溫度下的聲子總能量，即

式中，Dtot(ω)為聲子態(tài)密度；TMD和Tq分別為模擬溫度和修正溫度；ω 為聲子頻率；1/2 項(xiàng)代表量子效應(yīng)中的零點(diǎn)能量。量子修正后的熱導(dǎo)率表示為

本文計(jì)算得到的甲烷水合物的德拜溫度為165 K，此時(shí)的量子效應(yīng)不可忽略。圖2為計(jì)算得到的量子修正溫度和分子模擬溫度，可以看到在低溫時(shí)，修正后的溫度和模擬的溫度之間的差別很大，隨著溫度的升高兩者之間的差別逐漸減小。由于在量子力學(xué)中零點(diǎn)能量的存在，計(jì)算得到零點(diǎn)溫度為87 K，因此修正只限于MD 模擬中溫度大于87 K 的熱導(dǎo)率。圖3 為計(jì)算得到的修正系數(shù)，從零點(diǎn)溫度開(kāi)始，修正系數(shù)迅速增加；在德拜溫度附近，修正系數(shù)的增速趨于平緩；在高溫時(shí)，修正系數(shù)逐漸接近1，表明此時(shí)分子模擬可以得到準(zhǔn)確的熱導(dǎo)率。圖4為量子修正過(guò)后的熱導(dǎo)率，可以看到修正過(guò)后的熱導(dǎo)率與實(shí)驗(yàn)結(jié)果更加接近。

圖2 量子修正溫度和分子模擬溫度Fig.2 Quantum correction temperature and simulation temperature

圖3 量子修正系數(shù)Fig.3 Quantum correction coefficient

圖4 量子修正后的熱導(dǎo)率Fig.4 Thermal conductivity after quantum correction

2.2 甲烷水合物聲子模式分解

聲子是水合物導(dǎo)熱過(guò)程中的主要載體，是晶格振動(dòng)的量子化形式，研究聲子的性質(zhì)可以更深刻地了解材料的導(dǎo)熱過(guò)程，而聲子的弛豫時(shí)間屬于眾多聲子性質(zhì)中一個(gè)基本的物理量，只有知道了聲子的弛豫時(shí)間才可進(jìn)一步確定各種聲子模式對(duì)導(dǎo)熱的貢獻(xiàn)。對(duì)HCACF分解可以得到聲子的弛豫時(shí)間，如式(5)所示，聲子的弛豫時(shí)間一般包括兩部分，一部分為由聲學(xué)聲子引起的快速下降，一部分為由光學(xué)聲子引起的衰減[27-29]

通常情況下將聲子模式擬合為長(zhǎng)程聲子和短程聲子兩種模式，但當(dāng)溫度高于100 K 時(shí)，為了獲得更佳的擬合效果，在聲學(xué)聲子中引入了中程聲子項(xiàng)。而在光學(xué)聲子中引入常數(shù)項(xiàng)以表達(dá)聲學(xué)聲子及光學(xué)聲子之外的殘余振蕩[30]。

為了獲得聲子的弛豫時(shí)間，采用的擬合方法如下：如圖5(a)所示，為了計(jì)算光學(xué)聲子的弛豫時(shí)間，選擇HCACF 的峰值進(jìn)行分段擬合得到圖5(c)所示光學(xué)聲子的弛豫時(shí)間，從圖中可以看到，擬合可以分為短程、長(zhǎng)程光學(xué)部分以及常數(shù)項(xiàng)。將圖5(a)所示熱流自相關(guān)函數(shù)進(jìn)行傅里葉變換后，濾掉高頻部分，以忽略高頻光學(xué)聲子作用，再經(jīng)過(guò)傅里葉反變換得到低通濾波的熱流自相關(guān)函數(shù)[圖5(b)]。選擇相鄰波峰、波谷的時(shí)間中點(diǎn)進(jìn)行分段擬合得到聲學(xué)聲子的弛豫時(shí)間，如圖5(d)所示。從圖5(d)可以明顯看到，150 K時(shí)甲烷水合物的聲學(xué)聲子可以分為3段進(jìn)行擬合，即短程、中程和長(zhǎng)程聲學(xué)部分。

弛豫時(shí)間反映了相鄰原子之間的能量傳遞時(shí)間[31]。表1 給出了不同溫度下聲子的弛豫時(shí)間和光學(xué)模式峰值頻率，其中τsh,ac、τint,ac、τlg,ac分別代表短程聲學(xué)聲子、中程聲學(xué)聲子、長(zhǎng)程聲學(xué)聲子的弛豫時(shí)間；τsh,opt、τlg,opt分別代表短程光學(xué)聲子、長(zhǎng)程光學(xué)聲子的弛豫時(shí)間；ω 為光學(xué)模式的峰值頻率。從表1 可以看出，光學(xué)聲子和聲學(xué)聲子的弛豫時(shí)間均隨溫度升高而減小，這是因?yàn)閷?duì)于聲-聲U 過(guò)程散射，弛豫時(shí)間取決于頻率和溫度，三者之間的關(guān)系可以表示為τ-ω-2T-3。高溫下更多的聲子模式被激發(fā)，聲子之間的碰撞也會(huì)更加激烈，因此導(dǎo)致弛豫時(shí)間縮短。在150 K時(shí)，引入中程聲學(xué)聲子的弛豫時(shí)間τint,ac能夠使擬合效果更準(zhǔn)確，它的量級(jí)位于短程聲學(xué)聲子的弛豫時(shí)間τsh,ac和長(zhǎng)程聲學(xué)聲子的弛豫時(shí)間τlg,ac之間，是由高溫下主客體分子的耦合產(chǎn)生的。聲學(xué)聲子主導(dǎo)了甲烷水合物的熱輸運(yùn)過(guò)程，而光學(xué)聲子為聲學(xué)聲子提供了重要的散射通道。

表1 聲子弛豫時(shí)間和光學(xué)模式峰值頻率Table 1 Phonon relaxation time and peak frequency of optical modeontribution to thermal conductivity

不同溫度下甲烷水合物的熱流自相關(guān)函數(shù)如圖6 所示?？梢悦黠@看到，HCACF 從1 開(kāi)始迅速衰減，然后在零附近振蕩，約0.4 ps 后開(kāi)始趨近于零，低溫時(shí)小幅振蕩持續(xù)時(shí)間更長(zhǎng)。不同溫度下的熱流自相關(guān)函數(shù)在快速下降階段(聲學(xué)聲子主導(dǎo))的變化基本吻合，而在之后的衰減振蕩階段(光學(xué)聲子主導(dǎo))，開(kāi)始出現(xiàn)顯著的差異，隨著溫度的升高，HCACF的振蕩幅度逐漸減小，衰減得更迅速，這表明隨著溫度的升高，聲子的熱導(dǎo)率受光學(xué)模式影響的比重逐漸上升。

圖5 弛豫時(shí)間擬合Fig.5 Relaxation time fit

圖6 熱流自相關(guān)函數(shù)Fig.6 HCACF at various temperatures

圖7 為甲烷水合物30～150 K 的熱流自相關(guān)函數(shù)的能譜，它是通過(guò)對(duì)HCACF進(jìn)行傅里葉變換得到的?？梢钥吹剑哳l區(qū)域的幅值要遠(yuǎn)遠(yuǎn)大于低頻區(qū)域，表明低頻區(qū)域主客體分子之間的振動(dòng)充分耦合，高頻區(qū)域的聲子更容易在水合物中消散。振蕩特性是由于局域化的光學(xué)模式引起的，即水分子的振動(dòng)，能譜的峰值反映了孤立諧振子的局域化特性。隨著溫度的升高，能譜的峰值減小，表明水分子和甲烷分子之間振動(dòng)的充分的耦合[32]。

圖7 熱流自相關(guān)函數(shù)能譜Fig.7 Power spectra of HCACF

圖8 為甲烷水合物在30～150 K 的熱導(dǎo)率以及文獻(xiàn)中的結(jié)果[24]?？梢钥吹綇?0 K 到50 K，甲烷水合物的熱導(dǎo)率增大；隨著溫度的繼續(xù)增加，導(dǎo)熱率開(kāi)始下降，在100 K 達(dá)到最小值，然后隨溫度的增加繼續(xù)上升，這是由于高溫下聲子散射更加充分，聲子非彈性散射的增加導(dǎo)致能量傳遞通道增多，這與Krivchikov 等[24]實(shí)驗(yàn)結(jié)果變化趨勢(shì)一致。但是在數(shù)值上要比實(shí)驗(yàn)值偏大，這是由于本文模擬的水合物為滿晶穴，而實(shí)驗(yàn)樣品的晶穴占有率不可能達(dá)到100%造成的。與采用SPC/E 水分子模型的模擬結(jié)果相比，TIP4P/2005 模型得到的結(jié)果與實(shí)驗(yàn)值更接近，這與Jiang等[11]的結(jié)果一致。

圖8 量子修正后不同水分子模型下甲烷水合物的熱導(dǎo)率Fig.8 Thermal conductivity of methane hydrate under different water models after quantum correction

2.3 作用力強(qiáng)度對(duì)導(dǎo)熱的影響

為了研究客體甲烷分子和主體水分子間范德華作用強(qiáng)度對(duì)熱導(dǎo)率的影響，分別計(jì)算了不同作用力強(qiáng)度下甲烷水合物的熱導(dǎo)率。為了分析主客體分子對(duì)導(dǎo)熱的影響，將總的熱導(dǎo)率拆分為水分子對(duì)導(dǎo)熱的貢獻(xiàn)kww，甲烷分子對(duì)導(dǎo)熱的貢獻(xiàn)kmm以及水分子和甲烷分子耦合作用對(duì)導(dǎo)熱的貢獻(xiàn)kwm。溫度為200 K 時(shí)不同范德華作用強(qiáng)度下甲烷水合物的熱導(dǎo)率見(jiàn)表2。

表2 不同作用力強(qiáng)度下甲烷水合物的熱導(dǎo)率Table 2 Heat conductivity of methane hydrate with various strengths

由表2 可知，隨著主客體分子間范德華作用強(qiáng)度的增強(qiáng)，甲烷水合物總的熱導(dǎo)率、水-水作用、甲烷-甲烷作用、水-甲烷作用的貢獻(xiàn)均呈現(xiàn)增加的趨勢(shì)。同時(shí)，水-水作用的貢獻(xiàn)占總熱導(dǎo)率的比例最大，因此水分子形成的晶格主體對(duì)甲烷水合物的熱導(dǎo)率依然起著主導(dǎo)作用。此外，隨著作用力強(qiáng)度的增加，水-甲烷耦合作用在導(dǎo)熱中所占比例開(kāi)始增加。當(dāng)主客體分子相互作用較弱時(shí)，客體分子的局域化運(yùn)動(dòng)使得能量傳遞過(guò)程出現(xiàn)熱耗散，從而減少了熱量的傳遞，導(dǎo)致熱阻增加，熱導(dǎo)率減??；隨著主客體間相互作用的增強(qiáng)，客體分子運(yùn)動(dòng)的局域化特性不再明顯，客體分子的振動(dòng)與水分子晶格主體的聲學(xué)及光學(xué)聲子強(qiáng)烈耦合，從而導(dǎo)致熱導(dǎo)率的升高。

圖9 甲烷水合物的態(tài)密度圖Fig.9 DOS of methane hydrate with various strengths

圖10 態(tài)密度重疊區(qū)域的能量Fig.10 Overlapped phonon energy of DOS

圖9 為O 原子和C 原子的聲子態(tài)密度?？梢钥吹剑琌 原子和C 原子的振動(dòng)峰值分別在2 THz 附近。此外，甲烷分子的振動(dòng)在低頻區(qū)域局域化特征明顯，而在高頻區(qū)域局域化特征減弱。隨著甲烷分子與水分子間范德華作用強(qiáng)度的增加，受到主客體分子間共振散射的影響，甲烷分子的振動(dòng)峰值向高頻區(qū)域移動(dòng)，同時(shí)其峰值強(qiáng)度也減弱。氧原子的態(tài)密度隨著主客體分子間范德華作用強(qiáng)度增強(qiáng)，其峰值向高頻區(qū)域移動(dòng)。水分子和甲烷分子的態(tài)密度在低頻率范圍內(nèi)耦合程度較好，甲烷分子在低頻區(qū)域與籠形水分子間的非簡(jiǎn)諧勢(shì)能更加顯著。隨著相互作用力強(qiáng)度的增加，客體分子的局域化特征減弱，更容易靠近晶穴附近的水分子；甲烷分子的態(tài)密度強(qiáng)度下降，從而使主客體分子之間的耦合程度增強(qiáng)，熱導(dǎo)率增加。碳原子峰值頻率近似與作用力強(qiáng)度的平方根呈正比，在碳納米管共價(jià)鍵強(qiáng)度的研究中同樣發(fā)現(xiàn)這種變化規(guī)律。圖10 為不同作用力強(qiáng)度下甲烷水合物重疊區(qū)域的能量分布，隨著作用力強(qiáng)度的增加，C原子和O原子VDOS重疊區(qū)域的能量逐漸增加，與熱導(dǎo)率逐漸增加相吻合。

3 結(jié) 論

對(duì)分子動(dòng)力學(xué)模擬得到的甲烷水合物的熱導(dǎo)率進(jìn)行量子修正后可以更合理地描述水合物的導(dǎo)熱情況。甲烷水合物的導(dǎo)熱受聲學(xué)聲子及光學(xué)聲子共同作用的影響，其中聲學(xué)聲子在甲烷水合物的熱輸運(yùn)中起主導(dǎo)作用。隨著溫度的升高，聲學(xué)聲子及光學(xué)聲子的弛豫時(shí)間減少，更多的聲子模式被激發(fā)出來(lái)；同時(shí)光學(xué)模式的局域化特性減弱，主客體分子耦合更加充分。在甲烷水合物的導(dǎo)熱中，主體分子間的相互作用占主導(dǎo)地位，因此甲烷水合物導(dǎo)熱主要受水分子形成的晶格主體的影響，當(dāng)主客體分子相互作用增強(qiáng)時(shí)有助于甲烷水合物熱導(dǎo)率的提高。

符 號(hào) 說(shuō) 明

Dtot(ω)——聲子態(tài)密度

HCACF——熱流自相關(guān)函數(shù)

J——系統(tǒng)的微觀熱流，W/m2

k,kqc——分別為MD 模擬的熱導(dǎo)率和量子修正過(guò)后的熱導(dǎo)率，W·m-1·K-1

kB——玻耳茲曼常數(shù)，J·K-1

nac,nopt——分別為聲學(xué)聲子和光學(xué)聲子模式的種類

qi,qj——原子所帶電荷，e

r——原子間距離，?

S——標(biāo)度系數(shù)，用來(lái)調(diào)整范德華作用力強(qiáng)弱

T——模擬系統(tǒng)的溫度，K

TMD,Tq——分別為模擬溫度和修正溫度，K

U(rij)——原子i與j之間的相互作用勢(shì)，4.183J·mol-1

V——系統(tǒng)體積，?3

ε——?jiǎng)葳迳疃龋?.183J·mol-1

ε0——真空介電常數(shù)，F(xiàn)·m-1

σ——原子間平衡間距，?

τac,τopt——分別為聲學(xué)聲子和光學(xué)聲子弛豫時(shí)間，ps

ω——聲子頻率，THz

ω0,j——HCACF能量譜峰值的角頻率，rad·s-1

下角標(biāo)

i,j——原子編號(hào)