甲烷水合物声子导热及量子修正

刘明，徐哲

（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——原子编号