董京楠，班凡生，车阳，曲帅，王子健，刘智强

1 中国石油集团工程技术研究院有限公司，北京 100081

2 中国石油大学(北京)石油工程学院，北京 102249

0 引言

天然气水合物是以甲烷为主要成分，在高压力、低温度的条件下形成的一种笼型晶体，在全球范围内有着巨大的储量，也是公认的21世纪继天然气和石油能源之后最具开发潜力的清洁能源。在海底深水区和永久冻土深处广泛存在，有着丰富的储量[1-4]。在我国，除了东海、南海等海域均有大量的天然气水合物矿藏，青藏髙原等冻土区域也已证实了天然气水合物的存在，并成功获取了水合物天然岩心样品[5]。同时，据估计青藏高原冻土区天然气水合物资源具有1.2×1012～2.4×1014m3的巨大的储量[6]。

海洋天然气水合物的开采不仅受温度和压力的影响，同时开采过程的地质构造作用、海平面的升降以及海底工程都会使得开发十分复杂。基于天然气水合物在常温常压下极其不稳定性的特性，目前所提出的天然气水合物开釆方法基本都是利用降压、升温等手段打破沉积层水合物的相平衡条件，促使其分解，然后再收集溢出的天然气。具体的开发手段主要包括降压法、热激发法、化学试剂法和CO2置换开采法等[7-13]。水合物开采过程，温度压力的控制决定着开采的效率，而开采对储层的影响也不可忽视。由于海底表层沉积物大多属于黏土矿物，其强度较低，储层内的水合物晶体维持着基本骨架的稳定性。开采时温度和压力的变化会导致水合物分解成气体和水，水合物对骨架的胶结作用降低，从而降低冻土骨架的强度，进而导致井壁失稳。而无论是海上还是陆地，水合物的大规模开釆必然会导致水合物沉积层的强度降低，继而发生沉降，进而可能造成盖层的失稳，甚至会引发地震等危害[14-17]。

对于水合物分解引起的一系列问题，其中包含相变的传热问题是其中首先要面对的。在经典的传热问题中，主要是考虑材料的本身热导率、热容对传热过程的影响。同时，大多数研究都是在常温和高温环境下的传热过程。而对于水合物的开发，第一其环境温度较低，第二在开发过程中水合物会发生相变，同时吸收大量的热，会严重影响储层的温度场。这些问题相对于传统传热问题更加复杂，研究也相对较少，甚至是阻碍了水合物储层开发的关键理论基础之一。

本文主要研究在海洋区域天然气水合物开采时，井筒周围的地层温度场变化。通过理论分析建立了低温条件下含相变的岩石传热控制方程。在控制方程中，将比热和潜热合并成一个等效热容参数。基于Gibbs-Thomson方程[18-19]，提出了一种利用压汞数据得到孔隙内含冰量随温度变化的方法。采用所提出的模型，通过类比水合物开发过程水合物相变与冰水相变的本质，以水为介质来模拟水合物开采过程中井筒温度升高对周围地层的温度扰动，为进一步研究开采过程中储层的稳定性提供储层内相变传热的理论基础。

1 理论模型的建立

通常情况下岩石内的热量的传递由热导率决定，可以写出一般形式的热传导控制方程(即Poisson方程[20])：

其 中T为 温 度(℃)，t为 时 间(s)，ρ为 岩 石 密 度(kg·m-3)，c为 岩 石 比 热(J·kg-1·K-1)，k是 岩 石 热 导 率(W·m-1·K-1)，∇是梯度算子。而当孔隙中发生相变时，热量不仅会以比热的形式传递(体系温度发生变化，但组分相态不变)，还会以潜热形式发生转移(体系组分发生相变，但温度不变)。因此，含相变的控制方程可以写为：

其中L是水的潜热(3.33×105J·kg-1)，ρi是冰的密度(kg·m-3)，θi孔隙中冰的体积分数。当假设冻融过程都是瞬时的，将dθi/dt= (dθi/dT)(dT/dt)带入上式，则有

其中ρc-Lρidθi/dT这项定义为等效热容，从等式可以看出，孔隙含冰量随温度的变化使得传热过程变得复杂。

为了求解这个非线性方程，首先需要确定含冰量与温度间的关系(即dθi/dT)。冰在孔隙中发生相变，其相变温度不是一个常数，而是与孔隙的尺寸有关。岩石孔隙中的相变可以由Gibbs-Thomson方程表示[18]

其中ri是被冰晶所侵入的最小入孔半径，S是物质的摩尔熵(J·mol-1·K-1)，Vi表示冰的摩尔体积(m3·mol-1)，γiw是冰水间的界面应力(0.04 J·m-2)，α是冰水间的接触角(在这里取α=0)。Gibbs-Thomson方程揭示出在一定温度下，水只有在大于某一临界尺寸的岩石孔隙中才能够发生相变，孔隙越小，其间发生相变所需的温度越低。

由Gibbs-Thomson方程可以得到冰晶侵入孔径尺寸与温度的关系。但是需要注意的是，由于分离压的存在，在孔隙中，冰晶与壁面会存在一层较薄的液膜。这层液膜的厚度同样与温度有关[18]

其中ξ是表示分子间力的作用的特征长度(在这里取ξ～2.3Å)，D是 扩 散 系 数，并 且D=γsi-γsw-γiw(取D～ 0.33 J·m-2)。因此，发生相变的真实孔隙半径的大小是要略大于ri的(如图1)，并且有以下关系

图1 孔喉半径与温度压力间关系。r(T)：考虑液膜厚度的Gibbs-Thomson方程表示的温度和冰晶侵入孔径关系；r(P)：压汞法中Washburn方程表示的汞压与孔径关系Fig.1 Relation between the smallest pore access radius and temperature.r(T): relation between temperature and pore radius by Gibbs-Thomson equation considering the premelting film (G-T+e); r(P): relation between the pore radius and Hg pressure in mercury intrusion porosimetry (MIP).

所以，如果已知岩石的孔径分布，含冰量与温度间的关系就可以由方程(4)～(6)给出。在实验方法中 ，孔径分布数据一般是由压汞法取得，压汞法的原理就是Washburn方程

其中p是汞压(MPa)，σHg是汞的表面张力(0.48 mN·m-1)，αHg是汞与固体壁面的接触角(取cosαHg=0.765)。

结合方程(4)～(7)，可以得到一定孔径下相变温度与汞压间的关系，

由上式可以看出，含冰量随温度的变化可以用压汞法的数据来代替，即通过方程(8)将汞压转换成温度。然而，由于液膜的存在，压汞法求得的孔隙体积分数(θHg=VHg/V)要比同一孔隙中的含冰量(θi=Vi/V)要大。假设孔隙都是圆柱体孔隙，对应一定孔径的含冰量与真实孔隙体积之间的关系为:

应该注意的是，以上关于冰含量与温度间关系的推导是基于在每个孔隙中冻结和融化现象都独立存在这样一个假设。然而，由于结冰过程是在水中形成冰核，在此过程中产生冰水界面需要消耗表面能。因此，结冰过程往往是由外部孔隙向内逐步推进的，而不是在每个孔道中独立发生。换句话说，结冰过程依赖于孔隙的连接方式，一个孔隙的冻结依赖于它的入口孔隙半径大小。考虑一个入口孔径较小的大孔，温度降低时大孔不会先冻结，而需要等到入口孔喉结冰后，大孔才会结冰。相反，由于液膜的存在，融化过程不涉及新表面的生成，因此可以在每个孔隙中独立发生。因此，上述模型仅适用于不受孔隙连通性影响的融化过程。

2 储层的温度场模拟

在水合物开采过程中，井筒内流体温度相对较高，会对周围储层的温度场产生扰动。针对这一问题，建立如图2所示平面几何模型。模型中心的井眼直径设为200 mm，地层半径为10 m且为均质地层。注意到水合物开采过程中，井筒内的流体通过热对流方式将热量传递给井筒壁面，其后再以热传导方式在储层内传递，所以外边界设为自由边界，内边界设为对流边界。储层初始温度设为-4 °C，在使用隔热管的情况下，井内流体大部分热量被隔绝，那么设等效的流体温度为5 °C。

图2 有限元网格模型Fig.2 Finite element mesh model

利用实验室得到的一组砂岩的压汞实验数据求得其孔径分布曲线(如图3a所示，这里用岩石A来表示，峰值孔径约为10 μm)。同时，为了考虑其孔径分布对传热过程的影响，将岩石A的孔径分布曲线向右平移得到岩石B的孔径分布图(峰值孔径约为50 μm)。再通过方程(8)和(9)得到岩样A和B孔隙内含冰量与温度间的关系(如图3b所示)，岩石A和B的含冰量基本都在0 °C附近变化较快，推断在此区间内的相变反应剧烈，而岩样B含冰量大于岩石A则是由于岩样B的孔隙度较大。结合表1给出的各组分物质的物性参数以及传热过程的相关系数，得到等效热容值并代入，对地层的传热过程进行模拟。

表1 物性参数和传热系数[21]Table 1 Physical parameters and thermal coefficient in this work[21]

图3 a.岩样A和B的孔径分布曲线；b.孔隙内含冰量与温度间关系Fig.3 a.The pore distribution curve of sample A and B; b.ice content vs temperature in pores

在模拟结果中，基于模型轴对称的特点，以轴心为原点，选取任一半径为X轴。在1小时到30天之间的不同时刻，模拟得到地层温度场的演化云图(图4a)。取其沿X轴方向的地层岩石温度分布，得到图4b所示结果。在不同的t时刻，温度都呈梯状分布，沿X轴方向逐渐降低。同时，由于同一热源，初始相邻区域温差较大，温度的传播更快，其温度梯度更陡；随着热量向外的扩散，传播距离逐渐增大，同时温差逐渐减小，地层温度分布梯度逐渐变缓，可以预测一定地层的温度分布曲线是斜率逐渐变缓的过程。

图4 a.温度场演化；b.不同时刻地层温度分布Fig.4 a.the evolution of temperature field; b.temperature distribution of permafrost at different moments

对于岩样A和B，取井筒处同一点的温度变化曲线(如图5所示)。地层温度随时间逐渐升高，当达到零度附近时A，B曲线均出现比较平缓的一段，随后温度继续上升。B曲线在A曲线上方，表明其温度上升速率比A快。同时，假设温度随时间变化的速率小于0.01以视为温度曲线的平段区域(dT/dt＜0.01)，图中曲线B的平段要比曲线A长。

图5 岩石A和B的温度随时间的变化关系Fig.5 Temperature vs time for sample A and B

3 分析与讨论

在模拟结果中，相变潜热对地层的温度变化有极大的影响。在地层的传热过程中，孔隙内冰不断发生相变，平段区间的出现正是由于其孔隙中相变过程较为明显，大量的热量转化为物质的潜热，导致温度的上升极其缓慢。从控制方程(3)中也可以看出，在其他参数的一定的条件下，当含冰量随温度剧烈变化时，即潜热项的突然增大，等同于增大了整体的热容值，系统吸收相同的热量时其温度变化远没有之前明显。对于岩样A和B，在相似的孔径分布下，岩样B的孔径较大，其所对应的相变温度就较大，即在相同的热源下，岩样A比B更早的发生相变，即B曲线前期的升温速率要大于A。同时，在2种岩样的含冰量随温度的变化图中，零度附近B曲线的含冰量随温度的变化曲率大于A曲线，即在较小的温度区间内，岩样B孔隙中的含冰量变化等效增加的体系热容较大，所以升温更慢，温度曲线的平段更长。

由上可知，曲线的平段在模型中是由含冰量随温度的变化决定的，而含冰量随温度变化的关系是由压汞实验的孔径分布数据获得的，那么，通过测量材料在低温下的温度变化曲线，根据温度曲线平段的位置和长度就可以得到物质的孔径分布。值得注意的是虽然压汞法是实验室中一种较为常见的测量孔隙分布的方法，它的测量范围广，操作简单，但在理论模型中，利用压汞法来反求出岩石孔隙中的含冰量与温度间关系进行了两次变换，有液膜和温度压力间这两个过程转换误差的影响。因此，在条件允许的情况下，用更加精确的方法测量岩石孔径分布或者直接测出岩石孔隙中含冰量随温度的变化关系都可以对模型进行优化。

4 结论

针对水合物开发过程由于相变现象吸收大量的热，从而会造成储层温度场的扰动这一现象，类比冰水相变进行了储层温度场的扰动模型研究，并研究了不同孔径分布对储层温度随时间变化的规律的影响。主要得到了以下结论：

(1) 建立了含相变过程的传热控制方程。其中，将潜热和比热2种能量转化形式合并为等效热容，从而使控制方程简化为类似于不含相变过程的经典Poisson方程。另外，通过压汞数据可以预测含冰量与温度的关系。

(2) 当储层温度升高时，其地层温度变化曲线会在0 °C附近出现一段缓慢上升的阶段。这主要是因为在这一阶段相变剧烈，需要吸收大量热量，从而使温度变化并不明显。

(3) 储层的温度演化受岩石本身的孔隙结构影响巨大。当孔隙度较大时，由于骨架的缺失，岩石的整体热容会相应减少，从而使温度变化速率升高。同时，孔径分布也对温度平段的位置和长度产生。反之，通过测量储层的温度变化曲线，也可以反演得到岩石的孔径分布。