利用试井技术确定低渗透油藏CO2驱替前缘的方法

李友全，阎 燕，于伟杰

（中国石化胜利油田分公司勘探开发研究院，山东东营 257015）

CO2驱在中外油田中具有应用范围广、驱油效率高、成本低等特点［1-6］。CO2驱开采过程中，确定其前缘位置非常重要。但因CO2与原油间的相互作用，使CO2驱渗流规律较复杂［7-10］，确定CO2驱替前缘位置难度较大；数值模拟、试井分析等能在一定程度上解决该类问题［11-14］。目前，中外针对CO2驱试井分析的研究较少，且未有完善的试井模型。对于CO2驱试井理论研究成果主要集中在多区复合试井模型方面［15-21］，CO2驱低界面张力、降黏和原油膨胀等复杂的驱油机理未得到充分考虑，忽略了启动压力梯度的影响，难以真实地反映驱替过程中的压力动态。为此，笔者建立了低渗透油藏多组分CO2驱试井模型，采用有限体积法进行求解，分析了CO2驱试井曲线特征及曲线变化规律，实现了对CO2驱油藏动态参数的解释。

1 CO2驱数值试井模型

由于CO2驱油机理较复杂，使得地层中的流体性质及相态分布也变得复杂，进而对试井曲线产生较大影响。为准确认识试井曲线特征，采用状态方程、相平衡方程以及闪蒸计算描述气液两相间的传质现象及相态变化，基于组分模型建立CO2驱数值试井模型。

1.1 地质模型

基于胜利油田储层特征，建立单井及井网数值试井地质模型的主要步骤包括：①建立二维层模型。将沉积环境细分到微相单元，控制住砂体的几何形态；根据钻遇率确定砂体的分布面积；综合应用沉积、钻井、地震、测井和试井等多种方法确定砂体的厚度变化，从而确定砂体结构。②建立井模型。把各种储层信息转换成开发地质特征参数，在单井模型基础上得到该井的有效厚度、平均渗透率、平均孔隙度等参数，利用Lawson 算法实现Delaunay 三角剖分。③进行参数导入。根据网格剖分的坐标确定单井及储层渗透率、有效厚度和孔隙度等的相关参数，综合考虑计算速度和精度，如直井在模型的计算过程中优选出井筒周围的网格大小为0.5 m，而压裂直井在模型的计算过程中优选裂缝网格大小为1 m作为剖分标准。

1.2 数学模型

1.2.1 基本假设

油藏中存在油、气、水三相流体，nc个组分（包含CO2）；烃组分只分布在油相和气相中，水组分只存在于水相中；油藏处于恒温状态，流体流动处于热力学平衡状态；油、水两相渗流过程中考虑启动压力梯度的影响，符合非达西渗流特征，忽略气相中的启动压力梯度；岩石可压缩，且各向异性。

1.2.2 基本方程

通过物质质量守恒，可得每个组分所对应的质量守恒方程为：

对于油水两相，考虑启动压力梯度的运动方程为：

1.2.3 辅助方程

假设组分在相间的分布瞬间达到平衡，则可采用逸度守恒方程描述组分在相间的热力学平衡状态为：

液相和气相的逸度系数分别为：

（4）式也可写为：

为闭合方程系统引入约束方程式：

1.2.4 内边界条件

在数学模型求解过程中，对井筒的模拟采用井筒存储系数模型：

1.3 模型求解

1.3.1 基于有限体积方法的方程离散

在利用有限体积方法求解数学模型时，首先要对实际的连续介质进行离散化处理，即网格剖分，为了满足计算速度和精度的要求，对网格的形状、大小和数目等均进行优选。为考虑油藏及裂缝面的不同复杂边界，采用一维线网格描述裂缝，二维三角形网格描述油藏，即混合网格，采用二维非结构化网格对地质模型进行剖分。对模型求解区域进行网格剖分，离散得到若干子区域；把节点看作控制体积，控制体积上的物理量定义并存储在该节点上，并通过网格对数值模型进行离散化处理。

1.3.2 组分模型数值求解

组分模型的求解方法主要包括自然变量和摩尔变量2种。由于自然变量求解方法的未知量物理意义明确，编程实现方便，故基于自然变量对数值模型进行求解。

由组分模型可知，总共有2nc+4 个方程，变量个数为2nc+4 个，包括x1，x2，…，xnc，y1，y2，…，ync，p，Sg，So，Sw。为降低矩阵的维度，将方程分为主方程和辅助方程，将变量分为主变量和辅助变量，数值求解过程中，先通过求解主方程得到主变量，然后通过辅助方程获得辅助变量。

主方程为（1）式，主变量为y3,…,ync,p,Sg,So。辅助方程包括（4）式，（8）式，（9）式和（10）式，辅助变量为x1,…,xnc,y1,y2,Sw。

求解时基于非结构化网格对主方程进行离散，采用牛顿迭代法对非线性方程组进行全隐式求解。

1.3.3 闪蒸计算

组分数值模拟求解过程中，需要基于吉布斯能最小化原理判断单相是否稳定，若单相不稳定，则需要进行闪蒸计算，得到组分在油气相中的摩尔分数及油气饱和度。求解过程中，先通过二分法确定气液平衡常数的初始值，然后通过牛顿迭代法直接求解逸度守恒方程和摩尔分数的约束方程。求解过程中，逸度系数对xi和yi的导数需通过状态方程求取。实际闪蒸计算过程中，为了避免牛顿迭代法易造成局部收敛，难以得到精确解的问题，先通过逐次替换法逼近精确解，然后通过牛顿迭代法加快收敛速度。

2 压力响应规律分析

2.1 压力响应曲线特征

通过建立数值试井模型，求解得到了注入井压降试井压力动态曲线（图1）。从图1可以看出，在不同的渗流区域内曲线呈不同的延展规律。根据压力导数曲线的形态变化，将曲线分为5段：①井筒存储段。表征续流段影响，压力和压力导数曲线重合，斜率为1。②过渡段。井筒存储段结束后，受表皮效应影响，压力导数曲线出现峰值后下降，是压力波向CO2区传播的过渡段。③CO2区平面径向流阶段，压力导数曲线表现为水平。④CO2波及区流动段。压力导数曲线先上升后略有下降，造成压力导数曲线下降的原因是，关井导致井底压力在极短时间内变化较大；由于CO2黏度随着压力降低而降低，CO2波及区内的平均流度降低，最终引起压力导数曲线下降。⑤CO2未波及区流动段。由于该段流度小于CO2波及区，压力导数先上升，若油藏面积足够大，则CO2波及区和非波及区的流度相差较大，最终压力导数曲线将趋于平缓。

图1 注入井压降试井压力动态曲线Fig.1 Pressure behavior of drawdown well testing of injection well

2.2 CO2驱替前缘确定方法

为了更好地分析CO2含量在某个方向上的变化及其在压力特征曲线上的反映，采用矩形网格进行模拟，油藏尺寸为600 m×600 m，网格为30 m×30 m，见图2。

图2 CO2含量分布Fig.2 CO2content distribution

将网格点放大且标识每个网格，对应分析压力以及试井曲线的关系。试井曲线对应分析法的主要流程包括：①通过求解数值试井模型，可获得试井过程中每个网格的压力变化（图3a）；②将网格压力开始下降的时间作为该网格的岩石与流体性质（即流度）对试井曲线产生影响的时间；③通过分析试井过程中的岩石与流体性质分布，可以追踪油藏中岩石与流体变化对试井曲线的影响（图3b）。

图3 不同网格压力变化与双对数曲线的对应关系Fig.3 Relationship between pressure change and log-log plot of different grids

网格465 为CO2注入井所在网格，其CO2含量为100%；从网格466 开始，CO2浓度逐渐降低，地层流体黏度逐渐增大，流度逐渐减小，压力导数曲线先上升后略有下降；从网格469 开始，CO2含量进一步降低，流度大幅度减小，流度对压力导数曲线的影响占主导作用，使得压力导数曲线开始上翘；网格470为CO2驱替前缘。为了确定驱替前缘位置，在压力导数曲线后期阶段，于压力导数曲线近似直线处，作压力导数曲线的切线，切线与压力导数曲线开始相交处即为前缘位置。

2.3 参数敏感性分析

注入时间 由图4 可以看出，随着注入时间的增加，压力及压力导数曲线向右平移，CO2未波及区流动段的开始时间向后延迟；随着注入时间的增加，CO2波及区后期水平段、CO2未波及区水平段向下移动，说明注入时间越长，对油藏流体的降黏效果越明显，CO2波及区越大；另外，随着注入时间的线性增加，达到CO2未波及区的时间推迟，但推迟时间非线性增加，增加幅度越来越小。

图4 不同注入时间试井曲线对比Fig.4 Well testing curves at different injection time

渗透率 由图5可以看出，随着渗透率的增大，压力及压力导数曲线向下移动，说明地层压力降低；CO2未波及区流动段的开始时间提前，说明地层未波及区减小；随着渗透率的线性增加，压力下降幅度逐渐减小，压力导数曲线在CO2波及区的下降幅度减小。

图5 不同渗透率试井曲线对比Fig.5 Well testing curves with different permeability

3 实例分析

胜利油区樊142-7-斜4 井组于2013 年6 月开始实施CO2驱实验。该井组位于樊142-10 单元东北部，属于滩坝砂沉积、低孔特低渗透储层，该井组包括1 口注入井（樊142-7-斜4 井）和6 口生产井（樊142-6-2、樊142-6-3、樊141-1、樊142-8-斜4、樊142-8-3 和樊142-7-3 井）。利用注入井压降试井测试数据，通过模型拟合分析解释得到CO2驱替前缘位置即CO2波及区半径为79 m，表皮系数为0.1，井筒储集系数为1.4 m3/MPa，渗透率为0.15 mD，启动压力梯度为0.009 MPa/m，该解释结果反映了测试时刻CO2驱替前缘的推进状况。双对数曲线拟合情况见图6，拟合效果良好。

图6 注入井压降试井曲线拟合结果Fig.6 Log-log match result of drawdown test of injection well

根据数值试井模型、注水井及生产井的压力测试，对驱替前缘推进进行了跟踪预测。CO2注入初始阶段驱替前缘推进较缓慢，随着CO2注入量的增加，经过430 d 后，CO2驱替前缘推进速度加快，于630 d 后到达樊141-1 井，利用该预测结果设计了樊141-1 井的开井时间（图7）。该井于设计开井时间试产，稳定产油量由注气前2.0 t/d（泵抽）提高到5.2 t/d（自喷），CO2含量为48%，生产效果较好，进一步验证了试井解释结果的可靠性。

图7 樊141-1井CO2驱替前缘推进预测Fig.7 Well F141-1 CO2flooding displacement front predicting result

4 结论

基于组分模型建立了低渗透油藏CO2驱数值试井模型。通过数值求解绘制出CO2驱油藏试井典型曲线，曲线呈现井筒存储段、过渡段、CO2区平面径向流阶段、CO2波及区流动段和CO2未波及区流动段等5个主要流动段。分析了注入参数和渗透率等对曲线的影响规律。矿场应用表明，提出的CO2驱替前缘试井分析方法不仅能确定油藏渗透率，而且能有效确定驱替前缘位置，并可以进行驱替前缘推进情况预测。该方法对同类油田实施CO2驱动态监测有借鉴和指导作用。

符号解释

ϕ——孔隙度，%；np——相数量；i——序列号，其值为1,2,3,…,nc；xi,j——第i个组分在第j相中的摩尔分数；ρj——第j相的摩尔密度，g/mol；Sj——第j相的饱和度，%；t——时间，s；uj——第j相的速度，m/s；qj——源汇项；nc——组分数；vβ——流体的渗流速度，m/s；β——油相或水相；K——油藏渗透率，mD；Krβ——油相或水相相对渗透率，无量纲；μβ——流体黏度，mPa·s；G——低渗透油藏中流体的启动压力梯度，Pa/m；∇p——压力梯度，Pa/m；ρ——摩尔密度，g/mol；xi——油相摩尔分数；yi——气相摩尔分数；qi——第i组分的摩尔流量，mol/s；S——饱和度，%；o,w,g ——油、水、气相；——组分i呈液相的逸度，MPa；fiV——组分i呈气相的逸度，MPa；φiL——组分i的液相逸度系数，无量纲；φiV——组分i的气相逸度系数，无量纲；p——体系的平衡压力，MPa；qb——从地层流入井底的摩尔流量，mol/s；WI——井的生产指数，m3（/d·MPa）；Krj——油相、水相或气相相对渗透率，无量纲；μj——第j相的黏度，mPa·s；q——井筒到地面的质量流量，g/s；Mt——井筒混合物的摩尔质量，g/mol；ρt——井筒混合物的摩尔密度，g/mol；C——井筒存储系数，m3/MPa；pwf——井底流压，MPa。