一种新的水合物有效导热系数球形串并联模型

孙始财 谷林霖 杨震东 李艳敏

山东省土木工程防灾减灾重点实验室·山东科技大学土木工程与建筑学院

0 引言

天然气水合物（以下简称水合物）主要存在于海洋沉积物以及陆地冻土地区，因其全球储量巨大（约为2×1016m3）[1]，已成为国际认可的非常规战略能源之一。自2002年，水合物开发从勘查阶段进入勘查试采一体化阶段，近年正逐步向商业开采阶段迈进[2]。然而，由于水合物稳定域边界很敏感，在开发过程中温度发生变化时容易引起水合物分解或水合物二次生成，可能导致海底塌陷、滑坡和加剧温室效应等不可逆的环境问题[3-5]或管道堵塞、井喷和爆炸等安全问题[6]，掌握水合物的热物理学性质是安全、高效开发水合物资源的基础。导热系数是水合物重要的热物理学性质参数之一，水合物分解前缘主要受到传热传质的影响[7]，而传热的主导方式是导热，因而有必要聚焦水合物导热系数研究。

目前对水合物导热系数的研究主要集中在实验室测量和分子动力学模拟。例如，姚贵策等[8]实时测量在253～270 K、7 MPa时未经压实的甲烷水合物导热系数为0.44～0.46 W/(m·K)。黄犊子等[9-10]用十二烷基苯磺酸钠（SDS）溶液合成了甲烷水合物样品，并测得未经压实和通过2 MPa轴向压力压实的甲烷水合物导热系数分别约为0.33 W/(m·K)和0.57 W/(m·K)。上述研究[8-11]基本采用不压实或较低压力压实的水合物样品，其测得的甲烷水合物样品导热系数介于0.33～0.57 W/(m·K)。实际上，这些样品并不是纯甲烷水合物，其水合物固体颗粒之间存在孔隙，而孔隙中往往存在气体或者液体（水），从而影响导热系数测量精度。冰粉合成水合物（种冰法）的理论转化率可达到近100%[12]，因此作为水合物制备方法而被广泛使用，但其样品中仍不可避免地存在或多或少的孔隙[13]。研究者们为了得到更纯甲烷水合物的导热系数，对冰粉进行预先压实[14]或对已制备的样品进行高压压实[15-17]，以减小水合物颗粒间隙。然而，即使经过这样处理的样品也无法保证没有孔隙，即不能保证测量结果是纯水合物的导热系数。

为此，不少学者采用分子动力学模拟纯甲烷水合物的导热系数，其结果基本在0.56～0.75 W/(m·K)[18-21]，与温度成正相关性[18-19,22]；且在较高压力区间其与压力有显著的正相关[22]，而在相同温压条件下，晶穴占有率越大，导热系数越大[18,22-23]。虽然分子动力学模拟从理论上排除水合物孔隙的影响，并考虑了晶穴占有率、甲烷分子与水笼子的碰撞和范德华作用强度等因素，但其计算过程复杂，过于理论化且使用不方便。

综上，实验测量结果与其测量样品的品质关系甚大，而样品的品质由实验步骤、实验方法等条件综合决定。而分子动力学模型从理论上模拟水合物导热系数，其过程不依赖实验，但计算复杂，不利于工程实际应用。因此，笔者提出一种新的便于计算的水合物有效导热系数球形串并联模型，可以为水合物资源开发和衍生技术应用提供技术支撑。

1 球形串并联模型

对自然界水合物导热系数进行原位测量难度大，而实验室非原位测量需要保压保温取样保存技术[24]，其实施过程复杂，所以实验室制备与原位测量一体化的方案更具有可行性。但实验室制备的甲烷水合物样品往往不是纯甲烷水合物。因为甲烷水合物样品的固有特性和实验室样品制备方法的有限性，导致甲烷水合物样品中水合物颗粒间不可避免地存在孔隙。由于孔隙中残存着气体或水，其测得的导热系数不是纯甲烷水合物导热系数[25]。例如，水的导热系数约0.61 W/(m·K)[26]，与甲烷水合物导热系数接近，而气体的导热系数0.03 W/(m·K)[26]比甲烷水合物导热系数小一个数量级，故气体的存在导致测得的导热系数偏小。因此，笔者提出一种新的水合物导热系数球形串并联模型，以消除水合物颗粒间孔隙的影响。

1.1 模型假设

在样品制备初期，水合物颗粒间的孔隙比较大；随着水合物的生成，孔隙体积逐渐减小。当孔隙尺寸缩小到一定程度时，水合物逐渐停止生长[27-29]，此时水合物样品内部的孔隙结构确定。本文认为此时水合物颗粒孔隙中无流动的反应气和反应水，系统处于稳定。为了简化，将水合物颗粒间的孔隙假设为球形，且在样品中随机分布，球形孔隙外的区域皆为密实水合物，几何模型如图1-a所示。基于本文文献[27-29]，水合物样品中球形孔隙半径最小值为9.2 nm。假设半径为9.2 nm的球形孔隙是甲烷水合物不能继续生成的最小孔径，认为小于此孔径的孔隙中不会再有水合物生成。热阻网络图如图1-b所示。由于球形孔隙的对称性，几何模型仅绘制了球形孔隙的1/8，并以边长为L的立方体作为研究对象。通过此模型可以计算不同孔隙度下水合物样品热阻和导热系数。球形孔隙的球心位于坐标轴原点，球形孔隙半径为r1，热流量Q方向沿z轴正方向，立方体的热阻网络图如图1-b所示。根据热流路径可将模型分为两个区域分别计算。Ⅰ区热流依次经过孔隙流体相（气相或者水相）和水合物相，则区域Ⅰ的等效热阻RI为孔隙流体相热阻R1和水合物相热阻Rh的串联热阻（图1-b中左支路）；Ⅱ区仅含水合物相，即热流仅途径水合物相，故Ⅱ区等效热阻RII为水合物相热阻Rh（图1-b中右支路）。Ⅰ区和Ⅱ区等效热阻并联获得模型的总热阻。应注意，热阻与各区域（相）的导热系数、热流方向厚度及垂直热流方向的截面积有关。因此，水合物相热阻Rh和孔隙流体相热阻R1会因其所属区域（相）几何特点不同而不相等。各区域（相）的几何特点是由r1和L决定的。

图1 球形串并联模型示意图

1.2 模型建立

多相、多组分材料的导热能力通常用有效导热系数描述。本文中有效导热系数皆指包含水合物及孔隙流体（气体或水）双相的水合物样品导热系数。式（1）、（2）分别代表傅里叶定律和热阻原理：

式中λe表示水合物样品的有效导热系数，W/(m·K)；δ表示沿热流方向的单元体厚度，m；A表示垂直于热流方向的单元体截面积，m2；Rmol表示模型的总热阻，K/W；Rj表示j区域的等效热阻，K/W。

首先根据图1几何模型和热阻网络模型积分得到区域Ⅰ、Ⅱ的等效热阻方程为：

式中λh和λ1分别表示甲烷水合物和孔隙流体导热系数，W/(m·K)；μ1,h表示仅与水合物及孔隙流体自身导热系数相关的无因次中间变量；r1表示球形孔隙的孔径，m；L表示单元体边长，m；β表示孔隙度，指水合物样品中孔隙体积V1占样品总体积V0的百分数，即

Turner等[28]研究发现，当沉积物孔隙半径大于60 nm时，甲烷水合物相平衡不会发生变化，即影响水合物形成条件的最大孔隙半径为60 nm。因此当L为60 nm、r1为9.2 nm时，计算得到水合物样品的最小孔隙度β为0.001 9。将式（2）～（5）代入式（1），得到求解水合物样品有效导热系数的球形串并联模型方程：

1.3 参数取值

由于分子模拟的结果差异较大且过于理想化[18-23]；而实验测量结果依赖于样品品质和测量方法，这为参数选取制造了困难。根据式（6）可知，欲反算出纯水合物导热系数需要样品有效导热系数（λe）、样品孔隙度（β）及孔隙流体导热系数（λ1）这3个参数。Zhao等[14]采用预先压实的冰粉在温度273.35 K和压力10 MPa环境下合成甲烷水合物样品，并使系统温度在263.15 K和278.15 K之间周期振荡超过3次，最后对样品老化3天以上，从而确保冰粉完全转化为水合物。随后，在确保整个测试过程水合物样品不分解的前提下，采用双螺旋结构的Hot Disk探头测量不同温度压力下的甲烷水合物的导热系数。其预先压实的冰粉具有与冰相同的导热系数，很好地消除了冰粉颗粒间隙的影响[14]；其后续温度震荡、老化等步骤使可能残留的水（冰粉）得以完全转化为水合物。这种样品制备方法亦未采用SDS等表面活性剂，减少了测量干扰因素[9,15,17]。同时，Sun等[30-31]采取类似的方法制备水合物样品进行非原位测量，认为该方法制备的样品纯度高，且甲烷水合物样品的孔隙度介于19.8%～17.8%，故取平均孔隙度18.8%作为样品孔隙度（β）。将Zhao等[14]在263.15～279.15 K测得的水合物样品导热系数0.484 6 ～ 0.494 6 W/(m·K)作为样品有效导热系数（λe）。由于此样品在制备及测量过程中基本处于无液态水系统下，判断水合物颗粒间孔隙流体为甲烷气，故孔隙流体导热系数（λ1）取值 0.03 W/(m·K)[26]。用式（6）计算得到纯甲烷水合物导热系数随温度变化关系（图2），拟合这些计算值可以得到纯甲烷水合物导热系数的半经验方程为：

图2 纯甲烷水合物导热系数与温度变化图

式中T表示温度，K。

式（7）在4.0～9.5 MPa压力范围内计算纯甲烷水合物导热系数的误差为±0.007 W/(m·K)。Zhao等[14]认为在2～10 MPa范围内，气相压力对导热系数的影响很小；Rosenbaum等[20]亦发现在3.8～14.2 MPa范围内，压力与导热系数相关性很小。故式（7）忽略一定范围的压力对导热系数的影响。

将式（7）代入式（6），得到甲烷水合物样品有效导热系数球形串并联模型为：

式中B=λh/λ1-1 ；C1=1.948 89 ；C2=1.208 99 ；C3=1.240 70。根据孔隙流体的不同（水或者气体），λ1可取不同的参数，也可是关于温度的函数。式（8）适用于计算水合物颗粒间孔隙流体为单一相的水合物样品有效导热系数。当水合物处于生成或分解过程时，孔隙流体是气液两相共存[32]，可能发生微对流传热，不再是简单的导热问题，亦不是本文讨论的对象。

2 模型验证

由于标准的纯水合物样品及其导热系数尚未见公开报道，这给模型验证带来一定的困难。实验测量结果受到测量方法、样品制备方法和水合物颗粒间孔隙等多种因素影响，而纯理论的分子动力学模拟忽略了一些实际因素而不能完全反映真实情况。鉴于此，把实验测量结果和分子动力学模拟结果作为参考进行验证。首先将纯甲烷水合物导热系数半经验方程计算结果与实验测量结果[17,20,33]和分子动力学模拟结果[18,20,22,34]进行对比分析。从图3可以看出，Rosenbaum等[20]采用瞬态平面热源法（TPS）测得的261.5～277.4 K、4.6～6.3 MPa条件下压实的样品导热系数为0.68±0.01 W/(m·K)，比同温度下半经验方程[式（7）]计算结果低 0.08±0.01 W/(m·K)。除了测量误差的原因，也可能是在有限的实验温度压力范围内，逐渐缩小的孔隙阻碍了水合物继续形成，导致样品中存在少量的低导热系数气体而使实验结果偏小。Gupta等[33]和Li等[17]结合实验测量结果计算出零孔隙度纯甲烷水合物的导热系数分别是 0.70±0.04 W/(m·K)和大于 0.70 W/(m·K)。由于他们均未明确说明温度范围，猜测可能是实验温度的范围，即277.3～279.1 K和244.58～278.43 K。如图3所示，半经验方程计算结果比同温度范围内Li等[17]和Gupta等[33]的计算结果略高0.03～0.07 W/(m·K)。由于文献的计算模型没有考虑温度对导热系数的影响，本文模型考虑了此影响，而使模拟结果随温度增大而增大。但总体上来说，模拟结果比较接近文献结果。

图3 部分文献与本文得到的纯水合物导热系数图

除了文献实验测量结果，图3还给出了部分文献的分子动力学模拟结果[18,20,22,34]。从图3可以看出，分子模拟结果与实验测得的导热系数变化规律基本一致。但是分子模拟结果之间差异比较大，仅在265～280 K温度范围较为接近，因此采用该温度段分子模拟结果进行验证。从图3可以看出，半经验方程计算结果（图3中直线）比杨德伟等[18]的分子模拟结果约低0.01 W/(m·K)，比Rosenbaum等[20]分子模拟结果约高0.03 W/(m·K)，说明式（7）在263.15～279.15 K范围的计算结果比较准确。

另外，采用球形串并联模型［式（8）］计算本文文献[14]条件下甲烷水合物有效导热系数与文献实验测量结果的相对误差为0.05%～0.24%（表1），说明本文模型计算精度比较高。

表1 球形串并联模型模拟结果及误差表

式中εerror表示相对误差；λexper表示实验测得的甲烷水合物导热系数，W/(m·K)；λmodel表示球形串并联模型计算结果，W/(m·K)。

3 模型应用

本文提出的球形串并联模型既可以计算含不同孔隙流体及不同孔隙度的水合物样品有效导热系数，也可以根据实验测得的有效导热系数计算实验样品的孔隙度，籍以分析样品制备的纯度。

3.1 水合物样品有效导热系数

利用式（8）计算了孔隙仅含气或仅含水的水合物有效导热系数，λ1分别取值0.03 W/(m·K)[26]或0.57 W/(m·K)[35]，λh取值 0.765 4 W/(m·K)，结果如图4所示。由于水的导热系数与纯水合物的导热系数接近，所以孔隙流体仅为水时，水合物有效导热系数与纯水合物导热系数接近。而气体的导热系数远小于纯水合物导热系数，因而孔隙流体仅为气体时，水合物有效导热系数与纯水合物导热系数相比变化很大。相比于孔隙度参数，孔隙流体性质对水合物导热系数的影响更大。所以，样品制备过程如果水合物固体颗粒间隙比较大，必然影响测量精度。

图4 有效导热系数与孔隙度关系图

3.2 水合物样品的孔隙度分析

实验室中水合物样品制备方法直接决定样品品质，进而影响导热系数的测量精度。选取有代表性的几组水合物导热系数[8-9,14-17]，利用本文提出的球形串并联模型分别计算其孔隙度，结果如图5和表2所示，分析样品制备纯度。

表2 甲烷水合物样品的平均孔隙度与浮动值的计算结果表

图5 各文献样品孔隙度分析图

从图5及表2可知，姚贵策等[8]实验采用的样品未经压实，其有效导热系数最小，孔隙度介于20%～25%，平均孔隙度最大。而Waite等[16]持续用102 MPa压实4天所得样品的有效导热系数最大，计算得到的平均孔隙度明显小于其他样品。Li等[17]在 4 MPa、12 MPa和 16 MPa 压实压力下样品的平均孔隙度依次是17.23%、15.08%和14.05%，很明显地随压实压力的升高而减小。说明压实方法在一定程度上能够制备更纯的样品，但仍无法达到零空隙。溶解气法[9,15,17]制备的水合物样品孔隙度稍低（10%～18%），但是表面活性剂不参与水合物反应而残留在水合物颗粒间孔隙中可能会影响测量结果，进而使计算得到的孔隙度偏小。在没有很高的压实压力时，溶解气法制备的样品中残余气体较少。但考虑其样品可能受到活性剂的污染，种冰法[14]样品制备相对更好。从图5和表2可以看出，预先压实种冰法制备的样品[14]平均孔隙度为18.73%，浮动值为±0.002 3，说明每次制备的样品纯度基本一致，这样测量的结果更可靠。同样是采用种冰法制备样品，Waite等[16]实验采用的样品经过压实后平均孔隙度仅为7.34%，但是浮动值明显大于预先压实冰粉的样品[14]的浮动值，表明每次样品的纯度差异较大。种冰法可以使气体通过冰粉颗粒间隙相对均匀扩散，形成的水合物样品也相对均匀；而预先压实冰粉后冰粉颗粒更加密实，水合物样品中孔隙相对均匀和更小。在用SDS溶液制备的样品[9,15,17]中，黄犊子[9]实验采用的样品其孔隙度的浮动值低，即样品相对均一。Li等[17]压实压力和平均孔隙度都高于黄犊子[9]，但其浮动值接近，这可能是压实过程中造成更多裂隙，甚至可能造成水合物晶体破损。李栋梁[15]实验过程与Li等[17]类似，但浮动值差异较大；姚贵策等[8]实验采用的样品和Waite等[16]的制备方法差异很大，但浮动值接近。这些情况说明样品的孔隙度受多种因素影响。

基于上述分析，对水合物样品制备方法可以得到一些初步的认识。在每次制备的水合物样品差异程度上，预先压实冰粉优于一般压实操作；在水合物样品纯度上，压实样品优于未压实样品，种冰法优于溶解气法。

4 结论

1）本文提出的半经验方程在263.15～279.15 K的温度范围内可以比较准确地计算纯水合物导热系数，在4.0～9.5 MPa压力范围内的误差为±0.007 W/(m·K)。

2）利用球形串并联模型计算预先压实种冰法制备的水合物样品导热系数的值与实验测量结果误差为0.05%～0.24%，说明球形串并联模型可以用于预测已知孔隙度、孔隙流体类型及温度的样品有效导热系数。计算了含不同孔隙流体及不同孔隙度的水合物样品有效导热系数，发现孔隙流体类型对有效导热系数影响更为明显。

3）通过分析不同方法制备的水合物样品的孔隙度和纯度，初步获得了高纯度水合物样品制备方法的一些认识，即未经压实的水合物样品孔隙度介于20%～25%，普通压实后的样品孔隙度介于7.34%～17.23%，孔隙度总体上随压实压力的增大减小，从而可以提高样品纯度；样品中残余的表面活性剂使样品孔隙度减小，但是影响样品纯度；预先压实种冰法无需表面活性剂，可很好地制备高品质水合物样品。