低渗复合气藏斜井非达西渗流试井分析

伍锐东， 张春光， 杨志兴， 石美雪， 姜瑞忠

( 1. 中海石油(中国)有限公司 上海分公司，上海 200335； 2. 中国石油大学(华东) 石油工程学院，山东 青岛 266580 )

0 引言

低渗气藏常采用酸化等措施改善近井区物性条件，加剧储层的平面非均质性，呈现复合气藏的渗流特征。研究低渗复合气藏的渗流特征需考虑低渗区气体的非达西渗流、储层的应力敏感性和各向异性等因素[1-4]，采用传统渗流模型描述此类气藏存在局限性。

低渗气藏非达西渗流曲线由非线性段和拟线性段组成，非线性段对开发的影响不可忽略[5]。冯文光[6]分析流固表面作用、孔喉结构等因素对天然气非达西渗流的影响。李允等[7]基于实验研究启动压力梯度对气藏开发的影响，并建立拟启动压力梯度渗流模型。考虑启动压力梯度和应力敏感性，冯青[8]建立气水两相数值试井模型，并对相关参数进行敏感性分析。刘广峰等[9]应用幂律方程及拟线性方程的分段模型，描述流体的非达西渗流。欧阳伟平等[10]基于三参数连续非线性渗流方程建立致密气藏数值试井模型，表征气体的非线性渗流特征。低渗非达西渗流模型可归结为3类：拟启动压力梯度模型，忽略渗流的非线性段，放大渗流阻力；分段模型，涉及临界点判断，应用复杂；连续非线性渗流模型，多基于实验拟合，表征形式多样。为研究复合气藏的渗流特征，许峰等[11]与孟凡坤等[12]分别建立双重和三重介质复合气藏水平井不稳定渗流模型，谢飞[13]与陈军等[14]进一步研究复合油气藏水平井、垂直裂缝井的渗流模型。孙高飞[15]与黄雨等[16]分别考虑储层的应力敏感性和启动压力梯度进行低渗复合气藏试井分析。

目前，常规复合气藏的研究较为成熟，涉及低渗复合气藏斜井非达西渗流的研究较为鲜见，常采用简化的拟启动压力梯度模型进行表征，未综合考虑储层的应力敏感性和各向异性，不能充分体现低渗复合气藏的渗流特征。基于低渗气藏非达西渗流和应力敏感性微观机理，笔者考虑渗透率的各向异性，完善渗流运动方程，建立低渗复合气藏斜井非达西渗流数学模型，采用有限元方法进行Matlab编程求解，绘制井底压力动态曲线，对比不同模型渗流规律并对相关参数进行敏感性分析，应用模型进行实际斜井动态压力拟合并解释相关储层参数，验证模型实用性。

1 试井模型建立与求解

1.1 物理模型

在低渗气藏开发过程中，受工艺措施或气藏特点影响，易形成物性差异较大的复合气藏。

建立低渗复合气藏斜井物理模型(见图1)。假设条件：(1)气藏顶底封闭不渗透，气藏半径为re，地层等厚且厚度为h；(2)储层划分为内、外区，内区半径为r1，外区半径为r2，地层初始拟压力为mi；(3)斜井贯穿整个气藏，井斜角为θ，井长度为L，井中心的纵向深度为zw，以恒定地面产量qsc生产；(4)内、外区交界面不存在附加压力降，等温渗流，流体为单相微可压缩，忽略重力和毛管力；(5)考虑井筒储集系数、表皮因数及渗透率的各向异性，气藏水平方向初始渗透率为Knhi，纵向初始渗透率为Knvi(n=1,2，1代表内区，2代表外区)；(6)内区流体遵循达西渗流，外区流体考虑非达西渗流及储层应力敏感性。

图1 低渗复合气藏斜井物理模型Fig.1 Physical model of inclined well in the low-permeability composite gas reservoir

1.2 运动方程

低渗气藏渗透率低、束缚水饱和度高，气相吸附作用增强，孔道壁出现边界层，气源压力较低时气体出现明显的非达西渗流。低渗气藏非线性渗流模式见图2，其中δ为表边界层厚度，r0为毛细管半径，ν为通过岩心的流速，ζ*为最小启动压力梯度，ζ为拟启动压力梯度，p为压力。低渗储层流体渗流速度与压力梯度的关系曲线见图2(b)，由非线性段ζ*α与拟线性段αβ组成，α为临界点，η为临界点对应的压力梯度，β为渗流曲线的末端点。

图2 低渗气藏非线性渗流模式Fig.2 The nonlinear flow model of low-permeability gas reservoir

从微观角度解释非达西渗流的产生原因，考虑剪切应力及边界层厚度的影响，修正Hagen-Poisseuille方程[17]：

(1)

式中：K为岩心渗透率；μ为流体黏度；τ0为流体屈服应力；δ/r0表征边界层的影响，8τ0/3r0表征流体屈服应力的影响。

通过微观实验拟合，式(1)变形简化，得到宏观上的连续非线性渗流运动方程：

(2)

式中：c1、c2为通过实验拟合得到的参数。

为增强连续非线性渗流运动方程调整的灵活性，令a=-c2/c1，b=1/c1，则式(2)简化为

(3)

式中：a为非线性因数，反映流体屈服应力及边界层对渗流的影响；b为启动压力梯度的倒数。

当0

低渗储层具有较强的应力敏感性，随气藏开采的进行，岩石的有效应力不断增加，渗透率不断减小。为描述渗透率随压力的变化关系，采用PEDROSA O A定义的渗透率模量[18]：

(4)

式(4)积分可得

K=Kie-γ(pi-p),

(5)

式中：Ki为储层初始渗透率；γ为渗透率模量；pi为原始地层压力。

综上，考虑低渗气藏应力敏感性和非达西渗流的运动方程为

(6)

1.3 数学模型

1.3.1 模型建立

基于物理模型假设，联立非达西运动方程、质量守恒方程及状态方程，并引入拟压力函数可得低渗复合气藏斜井非达西渗流无因次微分方程：

(7)

(8)

式(7-8)中：m1D、m2D分别为量纲一的内区和外区压力；xD、yD、zD为量纲一的方向坐标；tD为量纲一的时间；rD、r1D分别为量纲一的渗流半径和内区半径；γmD为按拟压力定义的量纲一的渗透率模量；η1,2为内外区导压系数比；λmlD为量纲一的非线性渗流中间变量，l=x,y,z，分别代表x、y、z方向变量。

初始条件为

m1D|tD=0=m2D|tD=0。

(9)

定产内边界条件为

(10)

式中：zwD为量纲一的斜井中心纵向深度；εD为量纲一的微变量。

顶底封闭外边界条件为

(11)

水平封闭外边界条件为

(12)

式中：reD为量纲一的最大边界半径。

水平定压外边界条件为

m2D|rD=reD=0。

(13)

界面连接条件为

m1D(r1D,tD)=m2D(r1D,tD),

(14)

(15)

式中：M1,2为内外区流度比。

1.3.2 模型求解

非线性渗流数学模型解析求解较复杂，采用有限元法进行数值求解[19]。应用Galerkin法得到不存在源汇项时外区的有限元方程为

(16)

式中：Ni为形函数的分量。

通过格林公式分部积分，得到内部单元和具有封闭条件外边界单元的有限元方程：

(17)

将有限元方程变换成矩阵形式：

(18)

对有限元方程的矩阵形式进一步简化得

(19)

(20)

(21)

式中：Ke为非线性系数矩阵；Fe为单元载荷向量。

式(19)为该模型外区的有限元单元平衡方程，采用相同方法可得到内区的有限元单元平衡方程。

无限导流模型更符合斜井的实测压力分布，但求解困难，通常选取均匀线源模型等价压力点表征斜井压力分布[20-21]。采用CINCO L H等[20]提出的方法描述井底压力，井底等效压力点取值为

xD=0.3LDsinθw,yD=1,zD=-0.2LDcosθw,

(22)

应用Delta函数可得到与井位置有关的线源/汇公式，把斜井所在节点考虑成单元内源汇项，对井节点所在单元积分并无因次化，得到单元源汇项的有限元方程为

(23)

式中：d为线源划分的节点数。

通过Laplace变换并根据杜哈美原理，得到考虑井筒储集效应与表皮因数影响的井底压力解为

(24)

式中：S为表皮因数；s为Laplace变量；CD为井筒储集系数。

采用Stehfest数值反演，得到考虑井筒储集效应与表皮因数的斜井井底压力解为

(25)

(26)

2 渗流规律分析

2.1 流态划分

利用Matlab编程进行模型求解，绘制模型拟压力动态曲线，并与达西渗流模型和拟启动压力梯度模型进行对比[22-23](见图3)。模型基本参数：CD=100，S=1，hD=200，θ=60°，r1D=15，η1,2=0.2，γmD=0。非线性渗流模型考虑非线性特征和启动压力梯度，取a=0.5，bmlD=100；拟启动压力梯度模型仅考虑启动压力梯度，取a=0，bmlD=100。两个模型的量纲一的拟启动压力梯度为0.01。

根据拟压力动态曲线特征划分7个流动阶段：①井筒储存，拟压力及其导数曲线呈斜率为1的重合直线；②表皮效应过渡，拟压力导数曲线呈驼峰状；③井斜角控制，随井斜角的增大，逐渐呈水平井的渗流特征，井斜角大于60° 后逐渐出现早期垂向径向流和中期线性流，导数曲线分别呈水平直线和斜率为0.5的直线；④内区拟径向流，拟压力导数曲线呈值为0.5的水平直线；⑤内外区拟径向流过渡，拟压力导数曲线斜率与内外区导压系数比有关；⑥外区拟径向流，不考虑非达西渗流与应力敏感性等因素，拟压力导数曲线呈水平直线；⑦边界影响，受封闭边界影响，拟压力及其导数曲线剧烈上翘(见图3)。

图3 不同模型拟压力动态曲线Fig.3 Dynamic pseudo-pressure curves of different models

仅考虑外区的低渗特性，拟启动压力梯度模型和非线性渗流模型的拟压力及其导数曲线在外区拟径向流阶段发生上翘，非线性渗流模型的上翘幅度明显低于拟启动压力梯度模型的。

2.2 敏感性分析

分别对非线性因数a、井斜角θ、内区半径r1D及内外区导压系数比η1,2进行敏感性分析[24-26]。

2.2.1 非线性因数

非线性因数对拟压力动态曲线的影响见图4。由图4可知，受非线性渗流和拟启动压力梯度的影响，拟压力动态曲线在外区拟径向流阶段发生上翘。随非线性因数的增大，非线性段曲线的斜率减小，最小启动压力梯度逐渐趋于0，非线性渗流带来的附加渗流阻力减小，曲线上表现为外区拟径向流阶段上翘时间延迟，上翘幅度减小。非线性因数对低渗储层生产动态的影响明显，忽略非线性渗流特征的达西渗流模型或拟启动压力梯度模型，难以客观描述储层压力的动态特征。

图4 非线性因数对拟压力动态曲线的影响Fig.4 Effect of the nonlinear parameter on dynamic pseudo-pressure curve

2.2.2 井斜角

井斜角对拟压力动态曲线的影响见图5。由图5可知，在井斜角控制阶段，压力波尚未传远，井斜角对拟压力动态曲线的影响较为明显。井斜角较小时，拟压力动态曲线呈近似直井的特征；随井斜角增大，拟压力动态曲线逐渐呈近似水平井的特征。当井斜角大于60°后逐渐出现早期垂向径向流和中期线性流，其拟压力导数曲线分别呈一条水平线和斜率为0.5的直线。随井斜角的增大，井长度逐渐增加，早期径向流持续时间变长，定产生产时产量均匀分布，井底压力降低，拟压力及其导数曲线下移。不考虑储层应力敏感性和非达西渗流时，拟压力导数曲线在内区和外区的拟径向流阶段汇聚成一条直线。

图5 井斜角对拟压力动态曲线的影响Fig.5 Effect of the well angle on dynamic pseudo-pressure curve

2.2.3 内区半径

内区半径对拟压力动态曲线的影响见图6。由图6可知，内区半径影响内区拟径向流的持续时间和内外区拟径向流过渡阶段的开始时间：内区半径越小，内区拟径向流的持续时间越短，过渡阶段出现的时间越早，当内区半径小到一定程度时，压力波迅速传播到内边界，曲线上不能体现内区拟径向流的特征。2.2.4 内外区导压系数比

图6 内区半径对拟压力动态曲线的影响Fig.6 Effect of the inner radius on dynamic pseudo-pressure curve

内外区导压系数比对拟压力动态曲线的影响见图7。由图7可知，导压系数比主要影响内外区拟径向流过渡阶段曲线的斜率：内外区导压系数比小于1.0时，即外区渗透性好于内区的，内外区拟径向流过渡阶段的拟压力导数曲线呈下凹的特征；内外区导压系数比大于1.0时，即内区渗透性好于外区的，过渡阶段拟压力导数曲线呈上凸的特征。随内外区导压系数比的增大，外区的渗流能力相对逐渐变差，渗流阻力相对增大，内外区拟径向流过渡阶段的拟压力导数曲线斜率逐渐增大，外区拟径向流阶段的拟压力及其导数曲线位置也越高。

图7 内外区导压系数比对拟压力动态曲线的影响Fig.7 Effect of the pressure conduction coefficient ratio on dynamic pseudo-pressure curve

3 应用实例对比

东海XH气田A2井的测压数据表现出明显的低渗复合气藏斜井的渗流特征，分别应用常规试井模型和非达西试井模型对A2井进行解释(见图8)。由图8可知，常规试井模型仅考虑气藏的径向复合特征，不能体现外区拟径向流阶段的低渗特征；非达西试井模型考虑储层的低渗特性，精细刻画各流动阶段的渗流特征，模型与实测压力数据拟合效果更好，增加非线性因数、启动压力梯度等对拟合结果的约束，有效解决储层参数多解性强的问题。A2井的基础参数及试井解释结果见表1。

图8 A2井压力动态拟合曲线Fig.8 Dynamic pressure fitting curves of well A2

表1 A2井基础参数及试井解释结果

由A2井的压力动态曲线可知，该井的井斜角较小，无早期径向流阶段，类似直井的渗流特征。渗流后期压力传播到断层边界，压力及其导数曲线发生剧烈上翘。考虑低渗气藏非线性渗流和储层应力敏感性，改善该井外区拟径向流的动态压力拟合效果。非达西试井模型解释外区的渗透率为5.9×10-3μm2，表皮因数为8.3；常规试井模型解释的外区渗透率为3.2×10-3μm2，表皮因数为7.9。

基于产能试井实测数据，采用二项式产能方程计算该井的无阻流量为47.85×104m3/d；基于常规试井模型和非达西试井模型解释的参数，计算该井的无阻流量分别为32.49×104m3/d和45.44×104m3/d，非达西试井模型解释的参数更贴近储层实际(见图9)，其中Qg为测试流量。忽略低渗气藏渗流特性的常规试井模型对渗透率的解释与气藏实际的有效渗透率偏差较大；非达西试井模型的建立与应用可有效改进低渗复合气藏参数解释的精度，对低渗气藏的压力与产能评价具有重要意义。

图9 模型解释结果对比Fig.9 Comparison of different model interpretation results

4 结论

(1)基于低渗气藏非达西渗流和应力敏感性的微观机理，完善渗流运动方程，建立低渗复合气藏斜井非达西渗流数学模型；利用等价压力点处理内边界条件，采用有限元方法求解井底压力。

(2)绘制并对比不同模型的井底压力动态曲线，根据曲线特征划分7个流动阶段，即井筒储存、表皮效应过渡、井斜角控制、内区拟径向流、内外区拟径向流过渡、外区拟径向流及边界影响阶段。

(3)非线性因数增大，拟压力动态曲线在外区拟径向流阶段上翘幅度增大；井斜角大于60°后，逐渐出现早期垂向径向流段，且随井斜角增大持续时间延长；内区半径增大，内区拟径向流阶段持续时间延长；内外区导压系数比大于1.0时，过渡流段斜率为正数且随导压系数比的增大而增大。

(4)低渗复合气藏非达西渗流模型相较于传统试井模型压力动态曲线拟合更好，增加非线性因数等对解释结果的约束，降低参数的多解性；两种模型解释的渗透率相差明显，非达西渗流模型可以合理评估低渗气藏的产能。