甲烷水合物聲子導熱及量子修正

劉明，徐哲

（1 勝利油田分公司石油工程技術研究院，山東省稠油開采技術重點實驗室，山東東營257000；2中國石油大學（華東）新能源學院，山東青島266580）

引 言

甲烷水合物是由甲烷和水形成的一種籠型結構的化合物[1-3]，其熱導率在甲烷水合物勘探、開采、儲運以及其他應用過程中起著重要的作用，因此研究甲烷水合物的導熱具有重要意義[4]。

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

綜上所述，為了合理描述甲烷水合物在整個溫度范圍內的導熱情況，本文采用分子動力學模擬計算得出甲烷水合物的熱導率，并在一定溫度范圍內對其進行量子修正。通過對甲烷水合物的導熱進行聲子模式分解探究了甲烷水合物的導熱機理，并研究了主客體分子相互作用對甲烷水合物導熱的影響。

1 計算方法

1.1 模擬方法

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

式中，S 為調整范德華力作用強弱的標度系數；εij和σij分別為勢阱深度和原子間平衡距離；qi和qj為原子所帶電荷量；rij為原子間距離；ε0為電常數。模擬中設置范德華力截斷半徑為10 ?，Bugel 等[18]研究認為截斷半徑為2.5σ 時，截斷產生的誤差就可以忽略。采用pppm 算法處理長程靜電力相互作用，Luty等[19]認為pppm 算法在保證精度的條件下可以提高計算效率。選用velocity-Verlet 積分算法求解運動方程。

熱導率的計算使用G-K關系式[20]

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

1.2 模擬方法

平衡分子動力學模擬是在2 個×2 個×2 個sⅠ型甲烷水合物晶胞內進行的，其結構如圖1 所示。晶胞的長度為11.877 ?。

圖1 甲烷水合物模擬結構Fig.1 Initial configuration of methane hydrate

MD 模擬分為三步：首先將整體放置于NVT 系綜弛豫1 ns，采用Nose-Hoover 熱浴法[21-22]使系統達到預設的目標溫度，時間步長設置為1 fs。然后將系統轉入NPT 系綜，設置目標壓力為10 MPa，在目標溫度下弛豫1 ns，采用Nose-Hoover熱浴及壓浴法使系統達到平衡，最后將系統轉入NVE 系綜運行2 ns，計算熱導率。

2 結果與討論

2.1 量子修正

當模擬溫度遠低于德拜溫度時，大量聲子模式沒有被完全激發出來，此時的量子效應不可忽略，經典的MD 方法無法求出精確的熱導率值，可以考慮引入量子修正解決這個問題。量子修正的主要思想是通過經典力學系統中的模擬溫度來計算量子系統中對應的修正溫度。Lukes 等[23]應用量子修正方法準確計算出了單壁碳納米管的熱導率。量子修正假設模擬溫度下的聲子總能量等于量子修正溫度下的聲子總能量，即

式中，Dtot(ω)為聲子態密度；TMD和Tq分別為模擬溫度和修正溫度；ω 為聲子頻率；1/2 項代表量子效應中的零點能量。量子修正后的熱導率表示為

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

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

圖3 量子修正系數Fig.3 Quantum correction coefficient

圖4 量子修正后的熱導率Fig.4 Thermal conductivity after quantum correction

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

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

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

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

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

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

不同溫度下甲烷水合物的熱流自相關函數如圖6 所示。可以明顯看到，HCACF 從1 開始迅速衰減，然后在零附近振蕩，約0.4 ps 后開始趨近于零，低溫時小幅振蕩持續時間更長。不同溫度下的熱流自相關函數在快速下降階段(聲學聲子主導)的變化基本吻合，而在之后的衰減振蕩階段(光學聲子主導)，開始出現顯著的差異，隨著溫度的升高，HCACF的振蕩幅度逐漸減小，衰減得更迅速，這表明隨著溫度的升高，聲子的熱導率受光學模式影響的比重逐漸上升。

圖5 弛豫時間擬合Fig.5 Relaxation time fit

圖6 熱流自相關函數Fig.6 HCACF at various temperatures

圖7 為甲烷水合物30～150 K 的熱流自相關函數的能譜，它是通過對HCACF進行傅里葉變換得到的。可以看到，高頻區域的幅值要遠遠大于低頻區域，表明低頻區域主客體分子之間的振動充分耦合，高頻區域的聲子更容易在水合物中消散。振蕩特性是由于局域化的光學模式引起的，即水分子的振動，能譜的峰值反映了孤立諧振子的局域化特性。隨著溫度的升高，能譜的峰值減小，表明水分子和甲烷分子之間振動的充分的耦合[32]。

圖7 熱流自相關函數能譜Fig.7 Power spectra of HCACF

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

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

2.3 作用力強度對導熱的影響

為了研究客體甲烷分子和主體水分子間范德華作用強度對熱導率的影響，分別計算了不同作用力強度下甲烷水合物的熱導率。為了分析主客體分子對導熱的影響，將總的熱導率拆分為水分子對導熱的貢獻kww，甲烷分子對導熱的貢獻kmm以及水分子和甲烷分子耦合作用對導熱的貢獻kwm。溫度為200 K 時不同范德華作用強度下甲烷水合物的熱導率見表2。

表2 不同作用力強度下甲烷水合物的熱導率Table 2 Heat conductivity of methane hydrate with various strengths

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

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

圖10 態密度重疊區域的能量Fig.10 Overlapped phonon energy of DOS

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

3 結 論

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

符 號 說 明

Dtot(ω)——聲子態密度

HCACF——熱流自相關函數

J——系統的微觀熱流，W/m2

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

kB——玻耳茲曼常數，J·K-1

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

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

r——原子間距離，?

S——標度系數，用來調整范德華作用力強弱

T——模擬系統的溫度，K

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

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

V——系統體積，?3

ε——勢阱深度，4.183J·mol-1

ε0——真空介電常數，F·m-1

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

τac,τopt——分別為聲學聲子和光學聲子弛豫時間，ps

ω——聲子頻率，THz

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

下角標

i,j——原子編號