考虑动态相对渗透率曲线的油藏数值模拟方法*

付 强 杜志敏 王硕亮

（1.油气藏地质与开发工程国家重点实验室 四川成都 610500； 2.西南石油大学石油与天然气工程学院 四川成都 610500；3.中海油研究总院有限责任公司 北京 100028； 4.中国地质大学（北京） 北京 100083）

相对渗透率曲线是油藏工程和数值模拟的基础数据。目前油藏工程和数值模拟方法都认为相对渗透率曲线和毛管力曲线仅是润湿相（非润湿相）饱和度的函数，并不受润湿相（非润湿相）流速的影响。众多研究学者经实验验证已得到共识，稳态法测试得到的两相相对渗透率（非混相）曲线形态，不受流体流动速度的影响[1-3]。Ekwere J Peters[4]指出：由于测试方法存在差异，在油藏工程计算和数值模拟计算过程中，如果实际储层中的驱替过程与稳态法测试过程相同，则在进行油藏工程计算和数值模拟计算时应该采用稳态法相对渗透率曲线；如果实际储层中的驱替过程与非稳态法测试过程相同，则在进行油藏工程计算和数值模拟计算时应该应用非稳态法相对渗透率曲线。显然，非稳态法测试得到的相对渗透率曲线应用范围更为广泛。但是非稳态法相对渗透率曲线会受到流体流动速度的影响。Peters、Flock[5]指出，实际油藏的稳定系数往往会大于10 000，远远超出了稳定驱替的临界流速[6]，实际油藏数值模拟时，应该针对不同流速赋予不同的相对渗透率曲线。Leverett[7]提出，在室内实验过程中，测量得到的相对渗透率曲线与流速是有关系的。Rapoport[8]通过大量室内试验，证实了Leverett的观点，认为流体流动速度越大，润湿相和非润湿相的相对渗透率数值越大。众多研究学者（Avraam[9-10]、Avraam and Payatakes[11-12]、

Tsakiroglou[13]、Gutierez[14]、Tallakstad[15]、Erpelding[16]）利用相对渗透率曲线测试实验装置对混合润湿岩心进行了测试，测试得到了驱替过程和吸入过程的相对渗透率曲线，认为流速变化对驱替过程相对渗透率曲线和吸入过程相对渗透率曲线形态的影响规律是不同的，驱替过程中，相对渗透率曲线主要受到毛管指进现象的影响，吸入过程中，相对渗透率曲线主要受到黏性指进的影响。随着先进的无损检测技术的发展，核磁共振技术和CT扫描技术在相对渗透率测试领域的应用逐渐受到学者的重视。Krause[17-18]通过高精度的CT扫描，得到了岩心内部的孔喉结构分布，并利用图像分析方法，得到了不同驱替速度的相对渗透率曲线，Krause指出相同渗透率的岩心，驱替速度不同得到相对渗透率曲线不同。Konstantin Romanenko[19]利用核磁共振技术，测试得到了不同流速的相对渗透率曲线，明确指出流速会对相对渗透率曲线的形态产生影响，并且根据大量实验结果，提出了相对渗透率测试的合理流速范围。

实际油田地下储藏的流体流动速度的分布范围较广，水相流速最小值接近0，而水相流速最大值约为25 m／d，最大水相流速与最小水相流速相差很大。当油井高速开采后，流体的流动速度在空间上的分布变得更加不均匀。所以，常规方法测试得到的相对渗透率曲线（测试流速单一），已不能满足油田开发的实际要求。尤其针对海上油田，单井产能高，流体流动速度对相对渗透率曲线的影响已不可忽略，必须考虑低速、高速条件下相对渗透率曲线的获取及应用。

1 不同流速条件下相对渗透率测试实验

1.1 岩心流体和测试设备

实验设备包括注入控制系统、温度压力控制系统、岩心加持系统、出口回压控制系统和数据采集系统等5个模块以及Lab View实验控制软件。实验装置示意图如图1所示。

实验采用非稳态相对渗透率曲线测试方法。为了避免测试过程中流体流动对岩心孔喉结构的伤害，本次实验用岩心为5块真实岩心，分别来自中国海上某油田4口提液井（A井、B井、C井、D井），其中，A井取心2块（1、2号岩心渗透率分别为958、1 573 mD），B井取心1块（岩心渗透率为479 mD），C井取心1块（岩心渗透率为1 198 mD），D井取心1块（岩心渗透率为2 265 mD）。岩心润湿性都为中性润湿；岩心直径为2.5 cm，岩心长度为10～15 cm，岩心孔隙度为0.25～0.32，具体参数见表1。

图1 动态相对渗透率曲线测试实验装置图Fig.1 Dynamic relative permeability curve test device

表1 动态相对渗透率曲线测试实验参数Table1 Parameters of dynamic relative permeability curves

1.2 测试步骤

相对渗透率曲线测试实验步骤如下：

1）岩心抽真空，饱和地层水；

2）将岩心放入恒温箱中，排空阀门及管线，并联接好流程，在20℃下恒温10 h；

3）饱和原油，记录出口端产水量；

4）老化24 h后，将装满实验用水的容器抽空15 min；

5）将地层水注入，采用非稳态恒压法进行实验；

6）记录出口端产油量和产水量，实验进行到出口端含水率大于99.95%或者累积注水50 PV时结束；

7）绘制相对渗透率曲线；

8）变化驱替速度，重复上述过程。

实验用油是原油和柴油复配的混合液体，黏度为10 mPa·s；实验用水矿化度为4 000 mg／L，黏度为1.0 mPa·s；实验用油和水的黏度分别为3、0.59 mPa·s。

1.3 实验方案

采用不同驱替速度进行相对渗透率曲线测试，为了保证实验结果能够表征储层中的流体流动状态，采用相似准则和数值模拟2种方法，对相对渗透率曲线的测试流速进行了论证（表2）。选用中国海油某油藏实际温度60℃进行测试。

表2 动态相对渗透率曲线测试流速论证表Table2 Velocity calculation results of dynamic relative permeability curve test

1.4 不同驱替速度相对渗透率测试结果

通过对5块岩心不同驱替速度的非稳态相对渗透率测试，得到了5块岩心不同驱替速度条件下的相对渗透率曲线并计算岩心的驱油效率。

从5块岩心的测试结果（图2）可以看出，由于相对渗透率曲线测试的驱替速度不同，相对渗透率曲线形态差异明显；并且有3组岩心相对渗透率曲线形态存在一个共同的规律，当流速小于1 mL／min时，随驱替速度的增加，水相、油相相对渗透率均增加，残余油饱和度降低；当流速大于1 m L／min时，随驱替速度的增加，水相、油相相对渗透率均降低。

利用相对渗透率曲线计算得到岩心的驱油效率曲线，从驱油效率曲线（图2）可以看出，岩心绝对渗透率越大，驱油效率越高。岩心渗透率相同时，随着驱替速度的增加，驱油效率先增加后降低。以渗透率等于1 573 mD的5组岩心驱油效率结果（图2b）为例，可以看到，当驱替速度等于1 m L／min时，驱油效率达到极大值。其他渗透率的岩心测试结果也存在类似规律，说明在制定提液措施时，提液幅度并不是越大越好，而是存在一个最优值，最佳的提液幅度应该结合储层物性、储层非均质性和剩余油分布状态进行综合分析。关于驱替速度对驱油效率的影响规律，目前很多学者也发现了相似的结论[20-22]。

图2 5块岩心动态相对渗透率曲线与驱油效率测试结果Fig.2 Test results of dynamic relative permeability curve and oil displacement efficiency of 5 cores

将渗透率相近的5块岩心的所有相对渗透率曲线测试结果进行汇总，将油相相对渗透率或者水相相对渗透率作为因变量，将驱替速度（V）和含水饱和度（Sw）作为自变量，采用多元非线性回归方法，分别对实验数据进行拟合，得到了表征动态相对渗透率新型函数。

当渗透率小于500 mD时：

当渗透率大于500 mD小于1 000 mD时：

当渗透率大于1 000 m D时：

2 考虑动态相对渗透率曲线的油藏数值模拟方法

目前常用的数值模拟方法中，相对渗透率曲线是“静态”的，油相和水相的相渗透率只是含水饱和度的函数。静态相对渗透率曲线不能很好地反映变流速条件下的油水两相运移规律。根据上述研究成果建立动态相对渗透率计算公式，并代入到传统的数值模拟器中，对传统数值模拟器进行了修正。

以油水两相黑油模型为基础，将水相、油相运动方程中的静态相对渗透率替换为动态相对渗透率，对建立的新型油相、水相运动方程开展差分求解。考虑动态对渗透率曲线的数值模拟模型的前提假设条件为：

1）油藏内共有油、水两相；

2）油水之间没有质量交换；

3）流体流动规律符合达西流动；

4）油藏内温度恒定不变；

5）瞬时相平衡。

油、水两相模型渗流控制方程如下，其中

油相：

水相：

式（4～5）中：γo=ρog；γw=ρwg。

则油相控制方程的差分离散方程式可简写为

水相控制方程的差分离散方程式可简写为

3 数值模拟计算结果验证

3.1 数值模拟基本参数

根据上述理论方程的计算过程编写代码，建立考虑动态毛管力和动态相对渗透率的数值模拟计算程序。利用理想的一注一采模型，分析常规数值模拟和考虑动态毛管力、动态相对渗透率数值模拟的差异。数值模拟模型的基础取值如表3所示。

表3 数值模拟基础参数取值Table3 Basic parameters of numerical simulation

3.2 均质模型计算结果

为了分析动态相对渗透率对数值模拟计算结果带来的影响，利用本文所述的数值模拟器考虑动态相对渗透率和不考虑动态相对渗透率2种模式，计算得到油井含水95%时的含水饱和度分布结果（图3）。

图3 油井含水率95%时传统数值模拟方法与本文新方法含水饱和度分布对比Fig.3 Comparison of water saturation distribution of traditional numerical simulation method and new method when water cut of oil well is 95%

沿着注入井和采出井的连线方向提取了含水饱和度数值，得到了油井高含水时（油井含水率达到95%）的含水饱和度沿着注采井方向的数值（图4）。

图4 油井含水率95%时注采井连线含水饱和度数值Fig.4 Water saturation value on the line of injection well and production well when the water cut of oil well is 95%

从图4可以看出，油井含水率95%时，由于油井附近存在含水饱和度上翘的现象，考虑动态相对渗透率与不考虑动态相对渗透率的计算结果存在较大差异。水井端附近，含水饱和度和含水饱和度梯度都较高；油井端附近，含水饱和度和含水饱和度梯度也较高；油水井间的含水饱和度低于生产井和注入井附近的含水饱和度，含水饱和度剖面呈现出一个与压力剖面类似的“漏斗”形态。

均质模型油井的含水率曲线见图5。对比考虑动态相对渗透率数值模拟器的含水率曲线和不考虑动态相对渗透率的含水率曲线，可以看到，2种方法计算得到的见水时间相同。2条含水率曲线的差异主要体现在含水快速上升段，考虑动态相对渗透率数值模拟方法计算得到的含水率略高于不考虑动态相对渗透率数值模拟方法计算值。在高含水阶段，考虑动态相对渗透率数值模拟方法计算得到的含水率上升速度缓慢，与传统数值模拟方法计算得到含水率逐渐接近，但始终高于传统数值模拟方法计算得到的含水率。

图5 均质模型油井含水率曲线对比图Fig.5 Comparison chart of water cut curve of homogeneous model oil well

3.3 数值模拟计算结果的验证

为了进一步验证考虑动态相对渗透率数值模拟方法的准确性，建立了实际井的数值模拟模型。选择中国海上某油田A6H和B1井作为实例验证的单井，分别利用传统数值模拟方法和本文方法对A6H和B1井数值模拟模型进行数值模拟计算，得到的单井含水率曲线如图6、7所示。

通过对比A6H井实际含水率曲线和数值模拟计算得到的含水率曲线可以看到，本文方法计算得到的含水率曲线与实际生产曲线更为接近，说明本文方法具有较高的准确性和可靠性。

图6 中国海上某油田A6H单井含水率曲线Fig.6 Water cut curve of Well A6H in a China offshore oilfield

图7 中国海上某油田B1井提液后含水率曲线Fig.7 Water cut curve of Well B1 in a China offshore oilfield after enhanced liquid measurement

从图7可以看出，B1井在2011年10月8日和2012年1月2日分别进行了2次提液措施，B1井提液后，含水率曲线出现了2次先下降后上升的“漏斗”形态，采用传统数值模拟方法，无法精细描述提液后含水下降的现象；采用本文提出方法，可以相对准确的拟合出“漏斗”形态。通过实际B1井的模型，说明本文方法可以有效描述油井提液后的含水率变化。

4 结论

1）以中国海上某油田储层实际岩心为研究对象，开展了5个驱替速度条件下的非稳态相对渗透率曲线测试实验。分析实验结果可知，驱替速度越大，水相、油相相对渗透率数值均先变大后变小，存在临界流速，当驱替速度等于临界流速时，岩心的产液指数数值最大。

2）将水相、油相控制方程中的静态相对渗透率计算方法更改为动态相对渗透率计算方法，对考虑动态相对渗透率的基本控制方程进行了差分求解，在每个时间步计算相对渗透率时，综合考虑了该时间步流体流动速度和含水饱和度对相对渗透率的影响。建立了新型的能够考虑动态相对渗透率的数值模拟器。

3）利用平面均质模型，分析了本文方法与静态相对渗透率数值模拟方法的计算结果差异，发现考虑动态相对渗透率会增加油井附近和水井附近的含水饱和度，考虑与不考虑动态相对渗透率数值模拟方法计算得到的见水时间相似，考虑动态相对渗透率数值模拟方法计算得到的含水率高于静态相对渗透率数值模拟方法计算得到的含水率。

4）经过与物理模拟方法和实际生产数据的对比，说明本数值模拟器准确可靠。

符号注释

Sw—含水饱和度，f；

So—含油饱和度，f；

T—时间，s；

K—渗透率，mD；

Kro—油相相对渗透率，无因次；

Krw—水相相对渗透率，无因次；

Krod—动态油相相对渗透率，无因次；

Krwd—动态水相相对渗透率，无因次；

V—驱替速度，m／s；

ρ—密度，kg／m3；

μ—黏度，Pa·s；

p—压力，Pa；

γ—重度，N／m3；

D—网格所在高度，m；

q—产量，m3／d；

φ—孔隙度，f；

RK—渗透率下降系数；

i，j，k—表示x轴、y轴、z轴方向上的网格序数；

n—时间步长；

T—传导率，kg／（mPa·s）；

Φn+1o—第n+1个时间步的油相的势，m2／s；

Φn+1w—第n+1个时间步的水相的势，m2／s。