天然气水合物降压开采分解前缘移动数值研究

2020-12-31彭盈钰苏正刘丽华金光荣魏雪芹

海洋地质与第四纪地质 2020年6期

彭盈钰，苏正，刘丽华，金光荣，魏雪芹

1. 中国科学院天然气水合物重点实验室，中国科学院广州能源研究所，广州 510640

2. 中国科学院大学，北京 100049

天然气水合物是由水分子与气体分子组成的笼状结构化合物[1]。在自然界中，水合物中的气体分子主要是甲烷，因此是潜在的能源，同时也成为潜在的环境灾害因素。勘探结果表明，天然气水合物主要存在于海洋沉积物以及陆地冻土地区，其全球储量巨大，约为 2×1016m3[2]。

低温高压是天然气水合物稳定存在的必要条件，水合物的开采原理就是破坏水合物相平衡状态，使其分解为水和气体[3]。目前，常规的开采方法主要有3 种：降压法[4]、注热法[5]以及注入抑制剂[6]。现场试采、实验以及数值模拟结果表明，降压法是最具潜力的开采方法，而其他开采方法适合作为辅助手段来提高水合物的开采效率[7]。

开采过程会破坏天然气水合物的相平衡状态，使储层划分为不同的区域：分解区域、未分解区域，在两者之间存在一个特殊的过渡区域，即分解前缘。以水合物分解前缘为界的分解区域与未分解区域之间的流体性质、储层地质条件以及地层温度压力规律等都不相同[8]。在开采过程中水合物分解前缘移动规律与水合物开采动态有着紧密关联，分解前缘能直接反映了水合物开采特征。另外，由于水合物在储层中起到了一定的胶结支撑作用，随着分解前缘的移动，地层稳定性降低，可能引发地质灾害[9]。因此，在水合物开采过程中，研究水合物分解前缘移动规律具有重要指示意义。

目前，对于天然气水合物分解前缘移动问题的理论研究多基于Stefan 边界理论，即将水合物分解过程类比冰消融过程，是一个伴随相变的传热过程。Makogon 首次借鉴Stefan 问题，计算得到了降压分解水合物过程中压力分布的自相似解[10]。Verigin利用Stefan 移动边界问题，建立了一维半无限大水合物藏降压开采模型，模型考虑分解前缘两侧的气相流动以及分解前缘处气体质量守恒，根据模型与Stefan 问题的相似性对模型进行线性化自相似求解[11]。Ji 等在Verigin 模型基础上，建立水合物藏降压开采数学模型，假设水相静止，考虑了气相流动以及温度的变化，模型认为对流传热的作用比热传导强，分解区和水合物区的能量守恒方程中考虑了热对流以及节流和气体绝热效应，分解前缘处没有考虑水合物分解吸热的作用，将模型方程线性化处理后，自相似求解了分解前缘随时间的移动[12]。喻西崇等借鉴Ji 等提出的数学模型，利用自相似原理推导出分解前缘移动表达式与温度、压力分布表达式[13]。唐良广，李刚等将水合物分解过程看作移动边界问题，建立了水合物层温度分布的一维传热模型，模型考虑了分解区和水合物区的热传导以及分解前缘处的能量守恒，根据自相似求解得到不同时刻水合物藏温度分布以及分解前缘的位置[14]。张旭辉等建立二维热传导模型，研究水合物储层有热水管垂直穿过时水合物最大分解范围，采用分离变量法对模型进行求解，结果表明分解前缘的最大移动距离随温度的增大而增大[15]。刘乐乐建立水平一维降压-加热数学模型，将模型有限差分离散后数值计算得到分解相变阵面的位置与时间的平方根呈正比[16]。李明川等建立注热移动界面的三相一维传质模型，数值差分计算得到分解前缘移动速度前期较高，后逐渐降低[17]。另外，随着数值模拟器的成熟，Long 和 Tjok 利用 HydrateResSim 模拟水合物藏降压开采分解前缘移动，结果表明分解前缘的平均速度随绝对渗透率的增大而增大[18]。郑如意等数值模拟研究模型边界条件、渗透率、初始水合物饱和度、总热导率、井筒加热温度和井底压力等对水合物分解前缘移动的影响，发现分解前缘移动速度与渗透率、总热导率、井筒加热温度和边界供热成正比，相反，增大水合物初始饱和度和井底压力会降低分解前缘移动速度。此外，分解前缘移动规律也会随着参数的变化而变化[19]。

总的来说，前人采用数学建模以及数值模拟等多种方法对水合物分解前缘的移动规律进行了大量研究。与模型数值解相比，解析计算模型求解方便，计算过程也有助于深入了解某些物理变化的重要性。但大部分解析计算依据自相似原理，以分解前缘的移动与时间的平方根呈线性关系为前提假设，由此得到的分解前缘移动只在平均速度上存在差异，而分解前缘位置随时间变化规律已经确定，这将影响分解前缘移动规律探究。此外，利用现有数值模拟器探究水合物分解前缘移动规律时，要综合网格压力、相平衡压力以及水合物饱和度来判断分解前缘移动位置，这个过程相对复杂繁琐，并且单一储层模型的探究不具有普适性。

据此，本文建立了水合物降压分解一维三相数学模型，区别于传统自相似求解假设前提（即分解前缘的移动与时间的平方根呈线性关系），利用量级分析与偏微分方程无维化转化方法，解析计算探究了分解前缘随时间的移动规律，并由分解前缘移动计算相关产气量，对水合物开采动态进行了简单快速评估。

1 概念模型

模型为Class3 水合物藏。储层顶底板渗透率低，压力传递较慢，将顶部层和底部层认为是定温定压边界，允许发生热量和流水交换。

水合物分解是一个吸热的相变过程，会导致地层温度降低，此时，周围环境就会向水合物层传递热量。这种热补偿除了维持水合物分解所需热量外，还可以弥补吸热反应导致的地层温度下降。对于具有一定厚度的水合物层，可忽略顶底层围岩传入热量的影响，水合物分解相变前缘所吸收的热量主要来自于单位体积储层内能和未分解区热传导[20]。

随着时间的推移，水合物分解由近井区域向外扩散。假设水合物分解不是发生在整个储层内，而是发生在一定的狭窄区域内，可以将该区域视为一个表面，即所谓的水合物分解前缘[12]。它将储层分为两个部分，即水合物分解区与水合物未分解区。水合物分解区域自由气体在压力梯度的驱动下向井内流动，而分解前缘则向相反的方向移动。

根据前述假设，水合物分解过程可简化为分解前缘随时间向外移动的过程，而分解前缘就是天然气水合物发生分解的临界面，其厚度忽略不计，分解前缘处的地层压力即为该地层温度下的水合物相平衡压力[21]。在水合物层发生降压分解时，储层尺度内可采用水合物平衡分解模型，分解前缘界面（S（t）位置）把水合物层分为两个区域（图 1）：已分解气水区（r0＜x＜S（t））与未分解水合物区（S（t）＜x＜∞）。

根据概念模型的简化，本文进一步假设为：

（1）模型考虑三相（水合物、甲烷、水）两组分（甲烷、水），不考虑甲烷和水合物的溶解；

（2）由于气体与水之间的压力差不大，忽略毛细管压力的影响；

（3）在水合物开采中，扩散作用贡献小于对流作用，忽略气水扩散对水合物分解的影响；

（4）水合物分解所需能量主要包括分解前缘所在区域的单位体积储层内能以及未分解区域传导热。

图1 天然气水合物开采分解前缘迁移概念模型Fig.1 Conceptual model of migration of decomposition front of gas hydrate

2 数学模型

在数学模型中，将水合物储层看作x 方向的一维流体场。描述多孔介质中水合物分解过程的主要方程包括水合物分解区域和分解前缘处的质量守恒方程、能量守恒方程。其中，分解前缘移动相关方程参考Stefan 模型对冰水自由边界的描述[22]，其他方程则类似于Tsypkin[23-24]、Ahmadi[25]等提出的模型方程。

2.1 模型主要方程

水合物分解区域质量守恒方程：

当分解区与未分解区的压力梯度都比较小时，水合物分解区域与未分解区域的能量守恒方程简化为（体现水合物分解区域与未分解区域温度变化关系）：

式中：下标1，2 表示分解区与未分解区，λ 是导热系数，c是比热容，为拉普拉斯算子。

水合物分解前缘处质量守恒方程为：

水合物分解前缘处能量守恒方程为（体现水合物分解前缘上的温度变化）：

2.2 模型辅助方程

分解前缘压力为水合物分解为气水的相平衡压力[26]：

计算热传导系数方程：

2.3 模型初边值条件

模型初始条件与边界条件如表1 所示。

2.4 求解结果

我们将分解区气相质量守恒方程（1）与分解前缘气相质量守恒方程（3）作量级分析（具体步骤见附录A）。量级分析的思想是，如果方程是基于无量纲和归一化变量的表现形式，方程不同项的系数能用来度量这些项的重要程度[27]。无量纲化后的方程（1）中第二项系数远大于第一项系数1，可忽略方程（1）中第一项对方程的影响。另外，由于气水相运动黏度相差100 倍，在相同压力梯度下，水合物分解后，达到传输平衡，水相饱和度认为是不随时间变化。因此，在这个模型中水合物分解区流体流动简化为拟定常流动，得到分解区压力传导关系：

表1 初始条件与边界条件Table 1 Initial conditions and boundary conditions

根据分解区压力传导方程（7）与分解前缘质量守恒方程（3）得到分解前缘随时间移动规律：

由式（8）可知分解前缘的移动与流体相渗透率，水合物分解相平衡压力与井底压力之间的差值有关，与时间的平方根呈线性关系。

将分解区与未分解区能量守恒方程作无量纲转换后代入积分得到（具体步骤见附录A）：

将式（9），（8）代入式（4）并与式（5）联立得到一个用于求解分解前缘相平衡压的超越方程组：

单位横截面积生产井产气速率为：

由式（11）可以看到，产气速率与分解前移动距离成反比，说明随着天然气水合物的分解，产气量逐渐减少。

当分解前缘移动到S（t）时，单位横截面积生产井总产气体积为：

实际生产井总产气体积为：

式中：h为井射孔有效长度为井孔半径。

3 讨论

3.1 模型验证

Yousif 通过水合物砂岩样品降压分解实验，探究了分解前缘移动现象[28]。以Yousif 实验为依据，将实验中所设置的相关参数代入模型中，并通过对超越方程（10）以及式（5），（8）联立求解，得到分解前缘移动位置，实验结果与模型计算结果进行对比，如图2 所示，实验数据与模型计算结果近似，从而验证了模型结果的可靠性。

3.2 示例分析

以南海神狐海域天然气水合物藏[29-31]实际参数为例，模型计算所需的基本参数如表2 所示。不考虑水合物分解过程中冰的生成对模拟结果的影响，井底压力为 3 MPa。

图2 实验与模型下的分解前缘移动Fig.2 Model match of experimental results for location of hydrate dissociation front

表2 相关物性参数Table 2 The correlated parameters used for calculation

通过表2 所示模型参数，进行水合物储层降压开采模型计算，得到200 d 内分解前缘随时间变化规律见图3。在天然气水合物降压开采过程中，当地层压力低于水合物相平衡压力时，水合物开始分解，并出现水合物分解前缘。在开采60、120、200 d后，模型计算分解前缘随时间移动距离分别约为33.35、47.17、60.90 m。从图 3a 中，可以看到随着开采时间的推移，水合物分解前缘移动曲线斜率变小，说明分解前缘移动速度降低。

根据水合物分解前缘移动速率变化（图3b），看出水合物分解前缘移动速率在生产前期达到最大值，随开采时间推移，移动速率变慢最后将趋于平稳。出现这一现象的原因在于，当井底降压开始时，水合物平衡状态被打破，此时压降开始传递，井压与地层压力差作为水合物分解的主要驱动力。随着在地层中压力传递以及水合物分解过程能量消耗，导致地层温度下降，水合物相平衡压力与井压之间的压差减小，分解过程随之变慢；在开采后期，储层能量不足，水合物分解主要依靠储层热量的传导，分解前缘移动速率保持较低的平稳状态，在这种情况下分解主要受储层热物理性质的影响。此时，我们应该考虑储层注热等技术进一步促进水合物分解。

图3 水合物开采特征a. 分解前缘随时间移动规律，b. 分解前缘移动速率随时间变化，c. 产气速率随时间变化，d. 总产气量随时间变化Fig.3 Characteristics of hydrate productiona. The moving of decomposition front with time，b. the velocity of decomposition front varies with time，c. the rate of gas production varies with time， d. the volume of gas production varies with time

图3c 给出了甲烷产气速率（单位高度和单位宽度的产量）随时间变化，我们可以看到，生产井（单位高度和单位宽度）在200 d 内,最大产气速率约为250 m3/d，后期产气速率接近 20 m3/d。水合物分解气体体积变化与分解前缘移动有关，即井口产气速率受到水合物分解前缘移动的影响，因此，产气速率随时间变化趋势与分解前缘移动速率一致。同时，如图3d 所示，在200 d 内甲烷总产气量（单位高度和单位宽度的产量）为 18 000 m3。

图4 模型计算总产气量与试采结果对比Fig.4 The model calculated results vs the test result

为了进一步分析模型计算结果，将模型结果与2017 年南海神狐海域第一次水合物试采情况进行比对。根据中国地质调查局的报告，连续产气60 d后试采结束，累产气量3 ×105m3，平均日产气 5 000 m3[32]。另外，天然气水合物储层厚度为40 m，假设井射孔贯穿整个水合物层，井孔半径为0.15 m。我们将60 d试采数据进行曲线拟合，得到图4 所示模型计算总产气量与试采总产气量变化。从图4 可以看到，在开采10 d 内的两者产气量比较吻合，随着水合物开采过程的进一步推进，模型计算总产气值高于与实际开采结果，两者差异增大，但整体相对误差在可接受范围内。模型计算结果和实际试采结果之间存在较大差异的原因主要在于3 个方面，首先本文模型中认为井射孔贯穿整个水合物层，水合物降压分解出的甲烷气体能快速有效的向井口方向移动，减小了孔隙中气体积聚所导致的地层压力增大的影响，从而有利于分解甲烷气的产出、后续压降传递和水合物的进一步分解，同时，模型忽略了实际开采中气体的溶解与扩散等作用的消耗；其次本模型没有考虑水合物二次生成、冰的形成等对水合物分解过程以及甲烷气流动的影响，在实际开采中，当储层局部温度低于冰点以下时，会有冰的生成以及水合物的再次形成，固相物质的出现会降低地层孔隙度和地层绝对渗透率，从而影响压降的传递，水合物分解减慢，并阻碍流体向井口方向流动，导致产气量减小，最后，实际的天然气水合物藏是一个三维空间，水合物分解过程是发生在三维空间内的物理化学变化，本文仅是一维简单模型，不可避免与实际情况存在差距。因此，本文模型计算所得天然气水合物降压开采总产气量高于实际开采值，是对总产气量变化的乐观预测，能简单快速地为实际开采提供大方向的参考。

3.3 敏感性分析

为了分析水合物储层相关参数对水合物分解前缘移动距离的影响，本文采用单次单因子敏感性分析方法，除了变量参数外，其他参数均保持在前述的参考数值[33]。

3.3.1 储层初始温度

地层初始温度是影响天然气水合物开采的重要储层参数[34]。探究不同地层温度下单井降压开采天然气水合物分解前缘移动变化规律。为了不改变储层初始压力，将温度变化控制在满足水合物相平衡条件内。

在地层初始压力与生产井压不变的情况下，不同地层温度下（286.15～288.15 K）水合物分解前缘移动变化如图5 所示。可以看到，储层温度的改变显著影响水合物分解前缘的移动，当储层温度增大时，分解前缘移动距离增大，储层温度变化1 K 时，200 d 分解前缘移动距离在参考情况（287.15 K）基础上变化33%。出现这一现象的原因是，天然气水合物分解是一个吸热的过程，储层初始温度越高，储层所能提供水合物分解的能量越多，有利于水合物分解；此外，随着储层温度的升高，相应的水合物相平衡压力增大，水合物分解速度更快。因此，当水合物储层温度很低时，多考虑提高储层温度作为辅助手段来促进水合物分解，提高产气量。

图5 不同储层初始温度下水合物分解前缘移动距离Fig.5 The location of hydrate dissociation front at different initial reservoir temperatures

3.3.2 储层绝对渗透率

储层渗透率是反映流体运移能力的重要水力学参数[35]。最新的调查数据表明，南海储层类型主要为黏土质粉砂—低渗透粉砂，地层绝对渗透率平均为2～5 md[36]。据此，敏感性分析的地层渗透率范围为 2～5 md。

对于不同地层渗透率，图6 显示了分解前缘移动随时间变化。正如预期的那样，随着地层渗透率的减小，分解前缘移动速率变慢，移动距离减小。地层渗透率变化对水合物分解有明显的影响，较高的地层渗透率可以提高流体运移速率，有利于地层整体压降传递，为天然气水合物分解提供更大的驱动力，水合物分解前缘移动距离随之增大，产气量随之增加。当储层渗透率提高 1 md（以 2.9 md 情况为基准），分解前缘移动距离较参考情况增大17%。

3.3.3 储层孔隙度

另一个重要的地层参数是孔隙度，而渗透率又是与孔隙度有关函数，在研究储层孔隙度变化对水合物分解前缘移动影响时，根据Kozeny-Carman 方程[37]（孔隙度的三次方与渗透率之间成正比），在改变地层孔隙度的同时，地层渗透率也随之改变。

随着地层孔隙度的增大，分解前缘移动速率减小，移动距离也随之减小（图7）。在孔隙度为0.38时，此时水合物相平衡压力为3.05 MPa，井口和井口前端之间的压差很小，这意味着岩石的热量不足以用于水合物分解，此时只能通过水合物未分解区域的热量流入来获得额外的能量。在这种情况下，分解前缘移动由储层热物理参数决定。这是因为地层孔隙度很大时，含天然气水合物的沉积物单位体积潜热也较大，但由于天然气水合物含量高，比热低，热导率低，使得地层的热导率变小，因此,地层温度下降更快，从而不利于水合物分解。

图6 不同储层绝对渗透率下水合物分解前缘移动距离Fig.6 The location of hydrate dissociation front at different absolute permeability

图7 不同储层孔隙度下水合物分解前缘移动距离Fig.7 The location of hydrate dissociation front at different reservoir porosity

4 结论

（1）通过参数量级分析，将天然气水合物分解渗流场简化为拟定常场，积分求解得到分解前缘随时间移动规律，同时将求解温度场变化方程进行无维化转换，推导出温度分布方程，并根据Yousif 对含甲烷水合物砂岩样品降压分解实验进一步验证了模型结果的可靠性。

（2）通过南海神狐海域相关参数的实例分析，发现多孔介质中水合物分解前缘移动速率随时间减小，移动距离与时间的开平方呈线性关系；井口气体产量随时间减小，最后趋于稳定。另外，对比了南海水合物试采结果与模型计算的总产气量，发现模型计算值高于实际开采结果，且两者相对误差在25%范围内。

分析误差出现的原因，主要是建立的模型忽略了气体的溶解与扩散等作用，且没有考虑水合物二次生成、冰的形成等对水合物分解过程以及甲烷气流动的影响。因此，本文通过研究水合物分解前缘移动规律对产气速率、产气量进行了较为乐观的预测，为水合物开采潜力评价提供了一种新的简单快速的计算方法，同时对开采动态评估给出大方向的参考。

（3）通过对地层初始温度、绝对渗透率以及孔隙度敏感性分析发现，地层初始温度和渗透率与水合物分解前缘移动距离之间成正相关关系，初始地层温度对水合物分解过程影响显著，另外，地层孔隙度越大，分解前缘移动速率反而降低，移动距离减小，井口与分解前缘压差减小，此时分解前缘移动速率由储层热物理参数决定。

附录 A 模型求解

将分解区气相质量守恒方程（1）中基本物理量表示为特征值与无量纲数的形式：