考虑压敏的双重介质分形油藏非线性渗流模型

姜瑞忠，张春光，崔永正，张伟，张福蕾，沈泽阳

（中国石油大学（华东）石油工程学院，山东 青岛 266580）

0 引言

在储层开采进程中，地层压力下降使得渗透率与孔隙度发生改变，储层流动原理也随之改变。鉴于双重介质油藏的非均质性和复杂性，国内外许多专家用流固耦合的渗流理论对其进行研究，但流固耦合方法复杂且求解困难。考虑介质随压力的变化建立模型并求得解析解，是另外一种研究方法，但流体在孔隙内的非线性渗流规律被忽略。众多的局限性导致双重介质分形油藏开始被广泛研究。J.Chang等［1］首先提出分形油藏中的试井方法；S.Aprilian等［2］采用干扰试井的解释方法，拟合得到地层的储存系数以及传导率，从而推广了J.Chang所提出的模型；J.Acuna等［3］将分形油藏的试井模型应用到地热矿的2口井，即将分形油藏原理投入了实际应用并得到了矿藏分维；葛家理等［4］用分形维数Df和分形指数来描述自然裂缝的复杂程度与连通状况，建立了分形油藏低速非达西渗流模型，最后求得实空间的精确解；李其深等［5］以储层的均质性为基础建立了分形油藏的双渗模型；李顺初等［6-8］研究了封闭边界且考虑井筒储集系数和表皮效应影响下，存在变流率问题的井底压力拉氏空间解。

本文在前人研究的基础上，建立了考虑压敏的双重介质分形油藏非线性渗流模型；基于有限元方法对模型进行求解，并根据井底压力解绘制渗流特征曲线；对比各模型渗流规律，对参数进行敏感性分析。研究成果对试井分析、油气藏开发动态预测具有实际意义［9］。

1 分形理论适用性分析与应用

经典的渗流理论基于欧式空间并具有一定特征尺度，储层渗透率K是与特征长度有关的值，所以，达西定律在宏观尺度上不能直接应用于分形油藏。如果储层体系在统计意义下具有分形特性，即不再具有特征长度时，经典渗流力学就不能很好地表征实际储层的非均质性。与经典的渗流理论类似，分形网络中的渗透率一定反映分形结构的几何性质，所以，渗透率与分形维数以及异常扩散系数θ必然有关。综上所述，将分形参数应用到具有分形特性的油藏体系中，可以更精确地描述实际油藏的复杂多变性。

首先对具有分形特性的裂缝孔隙度φf和裂缝渗透率Kf进行分析。设每个节点体积相同，分形体内流体储集在体积为Vs的节点处，节点与生产井的径向距离为 r，节点密度为 N（r）。N（r）与 Df的关系［10］为

式中：α为渗流系统渗流节点的数量；rw为井筒半径，m；Kw为距离生产井 rw处的渗透率，10-3μm2；ds为谱维数；φw为距生产井rw处的孔隙度；Ct为综合压缩系数，MPa-1；μ 为流体黏度，mPa·s。

在双孔单渗的情况下，裂缝空间压缩系数定义为

式中：Cf为裂缝空间压缩系数，MPa-1；p为地层压力，MPa。

在欧式条件下，α与φf有关，用B表示相应对称性，即：

式中：B为描述对称性的几何常量。

如果整数维d=2，即考虑平面内嵌入分形体，则根据孔隙度定义可得到：

多孔介质孔隙度与其有效上覆压力呈指数关系。设在r=rw处的孔隙度为φfw，压力敏感性孔隙度φf（r）可表示为

式中：pi，pf分别为原始地层压力、 裂缝压力，MPa；Kfw为岩心分析得到的地层渗透率，10-3μm2；CK为渗透率模数，MPa-1。

分形渗透率与分形孔隙度的定义类似，与扩散半径也呈幂律关系。在生产进程中，储层流体压力的变化会导致微裂缝的张合，储层的绝对渗透率会发生改变。孔隙度会随基岩上覆压力的变化而改变，孔隙度的改变会使基岩通过流体的能力产生变化，这就是变形介质中渗透率与压力有关的原因。参考分形孔隙度定义，加上异常扩散系数就可得到分形渗透率Kf（r）：

2 新模型建立

2.1 模型条件

建立双重介质分形油藏物理模型。建模时对双重介质油藏渗流系统进行假设［11-12］：1）裂缝网络为分形维数的分形储层，该网络被嵌入到二维欧式空间岩块中；2）流体微可压缩，定产量生产；3）不考虑重力、毛细管力以及储层温度变化；4）流体被驱替是指基质渗流到分形裂缝网络，没有孔隙系统的对流，只有裂缝系统中流体直接流向井筒；5）裂缝网络在压力作用下产生微小形变，不影响分形特征，不考虑分形维数变化；6）考虑表皮效应、井筒储集效应和压力敏感性。

目前，研究低渗透储层渗流时，常用拟启动压力梯度模型、分段模型、连续模型以及低速非线性渗流新模型（C模型）进行表征。研究表明，C模型处理低渗非线性流更加灵活，拟合程度更高［13］，所以，本文将C模型中的低渗非线性流运动方程运用到双重介质分形油藏中，可得到：

式中：v为渗流速度，cm/s；c1为流体屈服应力与边界层对渗流共同影响下的参数；c2为边界层对渗流影响下的参数；▽p为启动压力梯度，105Pa/cm。

c1，c2可以通过函数非线性拟合得到。当c1=0时，该模型可简化为达西模型，当c2=0时，该模型简化为拟启动压力梯度模型。最小启动压力梯度▽pmin=c1+c2，大于0。

本文模型是在考虑表皮效应、井筒储集效应以及压敏的基础上建立的。在考虑表皮效应的条件下，引入等效半径rwe=rwe-S（S为表皮系数），将井筒外可能受污染或激励的储层考虑为井筒的一部分，则有：

式中：φwe为等效半径处裂缝孔隙度；Kwe为等效半径处裂缝渗透率，10-3μm2。

2.2 建立模型

对基质和裂缝系统的运动方程、状态方程、连续性方程以及恒定产量和无限大外边界条件进行组合变换，双重介质分形油藏非线性渗流模型［14-15］为

内边界条件：

外边界条件：

初始条件：

式中：λ为窜流系数；δLVD为无因次中间变量；pm为基质压力，MPa；ω为弹性储能比；C为井筒储集系数，m3/MPa；t为时间，s；lc为特征长度，m；h 为油层厚度，m；下标D表示无因次。

2.3 模型求解

将裂缝和基质系统分别用有限元表示。运用Galerkin法得到裂缝有限元方程：

运用格林公式，分步积分可以得到内部单元和封闭条件边界单元有限元方程：

式中：Ni为形函数的分量（i=1，2，3，…，n）。

则有限元方程矩阵形式可简化为

进一步简化单元平衡方程，得：

其中：

同理，对于基质系统可得到有限元单元平衡方程：

式中：A 为单元域面积，m2；Ke为单元刚度矩阵；Fe为单元载荷向量；上标m表示时刻变量。

以上为双重介质分形油藏非线性渗流试井模型的有限元单元平衡方程。求解时需要由基质场的上一时刻值求解裂缝场的压力解，再将该时刻裂缝场压力解引入基质场求基质压力解。

3 渗流规律

图1为不同模型的双重介质油藏压力及其导数曲线的对比图（CD=100，S=0，ω=0.02，λ=10-9）。达西模型中，CKD=0，c1D=c2D=0，θ=0，Df=2.0；非线性流模型中，c1D=c2D=0.5，其他条件与达西模型一样；分形非线性流模型中，CKD=0，c1D=c2D=0.5，θ=0.3，Df=1.9。在考虑压敏的条件下，可以根据实际情况调整CKD的大小。

图1 不同模型试井曲线对比

由图1可以看出，任何模型的试井曲线均存在4个特征流动阶段，即井筒储集段I、裂缝系统径向流段Ⅱ、基质向裂缝的窜流段Ⅲ以及系统总径向流段Ⅳ。基质向裂缝的非线性流影响，表现为压力导数曲线下凹段右移，且下凹程度变浅，该现象为非线性流典型特征。分形非线性流表现为压力及其导数曲线均整体上移，曲线末端发生上翘。在裂缝性低渗透储层中，压力敏感性与分形对试井曲线产生类似的影响，但分形对多个流动阶段均有影响，而压力敏感性通常只对曲线末端产生影响。

4 敏感性分析

图2为不同c1D，c2D时试井曲线的变化情况。

图2 非线性参数对试井曲线的影响

由图2可以看出，从基质向裂缝系统的非线性渗流主要影响窜流段压力导数曲线波动段下凹深度及位置。窜流的非线性越强，窜流强度就越弱，压力导数曲线波动段右移且下凹程度减弱。

CKD对压力及其导数曲线的影响如图3所示。随着CKD的增加，压力敏感性逐渐增强，曲线整体上移幅度变大，但CKD不影响窜流出现的时间与强度。当CKD＞0.05时，中后期压力及压力导数曲线上翘加剧。

图3 渗透率模数对试井曲线的影响

Df对试井曲线的影响如图4所示。Df为描述整个分形油藏内部孔隙结构复杂程度的几何标量。其对压力曲线的影响表现为：Df越小，无因次压力变化越快，且无因次压力曲线发生上翘的幅度越大，上翘出现的时间提前。

图4 分形维数对试井曲线的影响

θ对试井曲线的影响如图5所示。θ越大，裂缝网络连通性与孔隙结构越差，渗流阻力越大。初期θ影响较小，但随着时间增加，θ的影响越来越大，导致压力导数曲线上移幅度增大。分形特征使压力波在储层中逐渐传播并形成动态边界，导致总径向流阶段压力导数曲线上翘幅度增大，体现了分形油藏的典型特征。

图5 异常扩散系数对试井曲线的影响

5 结论

1）达西模型、非线性模型以及本文分形非线性模型的试井曲线均存在4个流动特征阶段，即井筒储集段、裂缝系统径向流段、基岩向裂缝的窜流段以及总系统径向流段。

2）非线性参数越大，压力导数曲线下凹段右移，且下凹程度变小，即窜流强度减弱。θ越大，裂缝网络连通性越差，且渗流阻力越大，初期阶段影响较小，但压力导数曲线上移幅度随着时间的增大而增大；Df越小，压力变化越快且曲线上翘幅度越大，同时上翘出现的时间也提前。