缝洞型碳酸盐岩油藏水平井分形非线性渗流

姜瑞忠，张春光，郜益华，耿艳宏，余 辉，李昊远

（1.中国石油大学（华东）石油工程学院，山东青岛 266580；2.中海油研究总院有限责任公司，北京 100028；3.西安石油大学地球科学与工程学院，西安 710000）

0 引言

缝洞型油藏的孔洞、基质及裂缝系统相互影响，为了描述其储层复杂的孔隙结构，出现了三重介质模型的概念［1-2］。描述缝洞型碳酸盐岩油藏的渗流特征，须充分考虑裂缝系统的分形特征、致密基质块的非线性渗流特征，以及储层的应力敏感性［3-5］。

国内外学者针对三重介质模型进行了深入的研究。姚军等［6］创新变井筒储存三重介质油藏试井理论，对介质间窜流、弹性储容比及变井储参数对压力响应的影响进行了研究。Camacho-Velazquez等［7］对所建立的三重孔隙模型考虑了溶洞系统与井筒连通和不连通2 种情况下的解，并对样板曲线进行了详细的研究。李成勇等［8］建立了三重介质油藏的水平井各向异性渗流数学模型，并对水平井的压力动态进行了系统的研究。张冬丽等［9］建立了三重介质数值试井模型，很好地反映出了缝洞型储层的非均质性特征，该数值试井模型可进一步考虑缝洞型油藏中具体的裂缝和溶洞，进而推广至两相流试井。Gomez 等［10］利用全局最优化解法对缝洞型储层试井曲线解释方法进行研究，提高了计算效率和试井解释精度。Li 等［11］对不同边界条件下缝洞型储层的动态特征进行了研究，形成了一套系统的刻画缝洞型储层动态特征的方法。不仅如此，分形理论在储层渗流领域的应用也取得了较大的进展。Chang 等［12］提出了分形裂缝网络嵌入基质的表征公式，并通过对扩散方程系数的适当修正来描述分形系统中的单相流动。Acuna 等［13］将分形理论应用于裂缝性储层试井中，并基于试井解释明确了分形参数的物理意义。Razminia 等［14］利用分形几何和分数阶微积分的概念分析了两区复合储层的数学模型，并推导了径向复合系统中井底瞬态压力解析解。张本健等［15］针对裂缝性气藏试井渗流问题，使用分形网络模拟气藏裂缝系统并将分形网络嵌入到气藏中，计算出试井模型的流动动态特征典型曲线，结果表明分形网络能够很好地模拟裂缝性气藏的渗流特征。

虽然目前对三重介质模型的研究较多，但仍存在无法精确表征储层非均质性特征的问题，而分形理论的应用使得非均质储层，尤其是裂缝的表征更灵活且更贴近实际，但该理论在缝洞型储层中的应用研究较少且未成体系。此外，对缝洞型储层致密基质块非线性渗流特征的研究也不可忽视。研究低渗储层最常用的渗流模型主要有拟启动压力梯度模型、分段模型及连续模型，这3 种模型均由不同数学函数拟合实验数据所得，缺乏物理背景，表征不够灵活。

针对缝洞型碳酸盐岩储层渗流模型目前存在的问题，综合考虑裂缝分形特征和基质系统非线性渗流机理，完善非线性渗流理论，建立三重介质油藏水平井分形非线性渗流模型；基于有限元原理求解水平井动态压力，并对比各模型的压力动态曲线，总结渗流规律。在分析非线性参数、分形指数及水平井长度等相关参数的敏感性之后，将所建立模型进行实际井动态压力拟合并解释相关渗流参数，以期验证模型的准确性及实用价值。

1 缝洞型碳酸盐岩油藏渗流机理

1.1 基质非线性渗流

将低渗储层微观渗流机理、边界层流动理论和毛细管渗流模型相结合，简化并改进非线性渗流模型。考虑边界层厚度与流体屈服应力值，在Hagen-Poisseuille 定律的基础上进行修正［16］可得

式中:v为通过岩心的流速，cm/s；K为岩心渗透率，mD；μ为流体黏度，mPa·s；r0为毛管半径，mm；δ为边界层厚度，mm；τ0为流体屈服应力，MPa为压力梯度，MPa/m。

对于相同毛管，δ/r0与压力梯度成反比［17］，故假设；对于同一种流体，屈服应力值保持不变，假设8τ0/3r0=a2。将式（1）变形并简化可得

式中:a1，a2均为实验拟合得到的常数值。

式中:c1和c2是通过对实验数据进行非线性拟合得到的参数，c1反映了流体的屈服应力以及边界层对渗流的影响，c2反映了边界层对渗流的影响。

当c1=0 时，式（3）简化为达西模型；当c2=0 时，式（3）简化为拟启动压力梯度模型。相当于启动压力梯度项。当v＞0 时最小启动压力梯度为启动压力梯度的存在与否以及取值大小均取决于c1，c2的值，即取决于矿场实际［18］。

1.2 裂缝分形特征

利用分形维数及异常扩散系数来表征分形裂缝中的渗透率和孔隙度，能够精确地反映地下储层复杂的裂缝形态。设分形体内流体聚集在某节点处，且节点体积处处相等，则渗流节点的数量可结合节点密度求出，即

式中:Df为分形维数；α为节点的数量；ds为谱维数；rw为井筒半径，m；Kfw为生产井处的渗透率，mD；φfw为生产井处的孔隙度；Cft为裂缝综合压缩系数，MPa-1。

考虑维数为d的欧氏空间内嵌入分形体，可得到用分形维数表征的孔隙度表达式为

式中:B为描述对称性的几何常量；V为储层总体积，m3；Vf为储层孔隙体积，m3；Vs为储集流体的节点体积，m3；r为距生产井的距离，m。

裂缝孔隙度主要受孔隙空间聚集方式（用分形维数表征）的影响，且与有效上覆压力呈指数关系，所以考虑应力敏感性的裂缝分形孔隙度表达式［19-20］为

式中:Cf为裂缝孔隙压缩系数，MPa-1；pi为原始地层压力，MPa；pf为裂缝系统压力，MPa。

裂缝渗透率不仅取决于孔隙空间的聚集方式，还与裂缝之间的连通性（用异常扩散系数表征）密不可分。储层压力变化会引起裂缝张合状态的变化，从而影响储层渗透率，在对应力敏感地层不稳定试井的研究中发现，渗透率与基岩上覆压力之间也成指数关系。在分形孔隙度定义的基础上引入异常扩散系数便可得到考虑应力敏感性的分形渗透率［19-20］，即

式中:αK为渗透率模数，MPa-1；θ为异常扩散系数。

2 缝洞型碳酸盐岩分形油藏水平井非线性渗流模型

2.1 物理模型

图1 物理模型示意图Fig.1 Schematic diagram of the physical model

建立圆形三重介质碳酸盐岩分形油藏水平井物理模型［21］（图1），并设该模型符合以下假定条件:①模型由基质、裂缝和溶洞系统组成，基质块向裂缝和溶洞系统发生拟稳态窜流，溶洞系统向裂缝发生拟稳态窜流；②外边界为封闭或定压条件，储层厚度为h，原始压力为pi，水平井半长为L，水平井纵向深度为zw，油藏半径为re；③考虑储层的表皮系数、井筒储集系数以及裂缝渗透率的各向异性，裂缝系统水平渗透率为Kfh，垂直渗透率为Kfv，基质系统渗透率为Km，溶洞系统渗透率为Kv；④水平井位于油藏平面中心，产量为Q）；⑤忽略重力与毛管压力的影响；⑥考虑裂缝系统的应力敏感性及裂缝的分形特征，基质向裂缝系统的窜流为非线性渗流。

2.2 数学模型

在笛卡尔坐标系中，将流体流动的质量守恒方程和状态方程与非线性运动方程［式（3）］相结合，可得到缝洞型碳酸盐岩油藏水平井分形非线性渗流模型［22］，即

式中:rD为无因次渗流半径；LD为无因次水平井半长；hD为无因次储层厚度；tD为无因次时间；pfD为无因次裂缝压力；pmD为无因次基质压力；pvD为无因次溶洞压力；δLVDf为裂缝的无因次中间变量；δLVDv为溶洞的无因次中间变量；λmf为基质向裂缝窜流系数；λvf为溶洞向裂缝窜流系数；λmv为基质向溶洞窜流系数；ωf为裂缝弹性储能比；ωv为溶洞弹性储能比；下标D 代表无因次；下标m，f，v 分别代表基质、裂缝和溶洞系统。

2.3 模型求解

对于缝洞型碳酸盐岩油藏水平井分形非线性渗流模型，首先采用隐式求解裂缝系统，应用Galerkin法，得到裂缝系统的有限元方程为

式中:Ni为形函数的分量，i=1，2，…，n。

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

由此可得，单元有限方程的矩阵形式为

假设E为单位矩阵，其他各矩阵的表达形式为

对矩阵方程进一步简化可得

其中Ke，Fe的表达形式为

式中:Ke为非线性系数矩阵；Fe为单元载荷向量；上标n表示时刻变量。

式（18）为三重介质分形非线性模型裂缝系统的有限元单元平衡方程，同理可采用显式求解方法得到基质和溶洞系统的有限元单元平衡方程，且计算结果表明隐式求解裂缝系统、显式求解基质和溶洞系统的方法符合计算精度。

考虑到无限导流能力模型比均匀流量模型求解困难，通常选取0.7L处为水平井在这2 种模型下的等价压力点，从而可以借助均匀流量模型评价水平井井底压力。通过把水平井考虑成单元源汇项并积分，应用Delta 函数将水平井划分成n个节点，并进行无因次化得到单元源汇项有限元方程。

所求压力解经过Laplace 变换后，通过Duhamel原理引入井筒储集与表皮系数:

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

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

3 渗流规律

图2 为碳酸盐岩油藏水平井的达西、非线性以及分形非线性模型的压力动态曲线。各基本参数取值如下:CD=100，S=0，ωf=0.01，ωv=0.1，λmf=0.001，λmv=0.000 5，λvf=0.1，LD=5，hD=200，αKD=0.02。假设非线性参数cD=c1D=c2D，分形指数β=2+θ-Df，则各模型其他参数的取值如下:达西模型cD=0，θ=0，Df=2；非线性流模型cD=0.3，θ=0，Df=2；分形非线性模型cD=0.5，θ=0.3，Df=1.9。在考虑储层应力敏感性的基础上，可以根据实际情况适当调整各模型的αKD值。

不同模型在不同流动阶段曲线的形态各异，但压力动态曲线仍然呈9 个特征流动阶段（图2）。以达西模型为例，具体流动阶段划分及不同流动阶段压力动态曲线特征［23-24］如表1 所列。通过分析和对比非线性、分形非线性模型的压力动态曲线与达西模型的压力动态曲线的异同点（图2）可发现，非线性参数主要影响基质-裂缝窜流强度，而分形指数使得溶洞-裂缝窜流后的流动阶段的压力动态曲线发生上翘。

图2 不同模型压力动态曲线对比1.井筒储存阶段；2.表皮效应过渡阶段；3.早期径向流阶段；4.早期线性流阶段；5.中期径向流阶段；6.溶洞-裂缝窜流阶段；7.溶洞-裂缝拟径向流阶段；8.基质-裂缝窜流阶段；9.系统拟径向流阶段Fig.2 Comparison of dynamic pressure curves of different models

表1 不同流动阶段压力动态曲线特征Table 1 Characteristics of dynamic pressure curves in different flow stages

4 敏感性分析

首先对非线性参数和分形指数进行敏感性分析，各基本参数取值如下:CD=100，S=0，ωf=0.1，ωv=0.5，λm=0.001，λv=0.1，LD=5，hD=200，αKD=0.02。图3 与图4 分别为非线性参数cD和分形指数β对压力动态曲线的影响。

非线性参数主要对基质向裂缝的窜流段有显著影响，基本不影响溶洞向裂缝系统的窜流段。随着非线性参数的增大，基质向裂缝的窜流强度减弱，持续时间变短。在压力导数曲线上表现为第2个下凹段延迟出现，下凹程度减弱，但窜流段结束的时间不受影响。

图3 非线性参数对压力动态曲线的影响Fig.3 Effect of nonlinear parameters on dynamic pressure curves

图4 分形指数对压力动态曲线的影响Fig.4 Effect of fractal index on dynamic pressure curves

分形指数对压力动态曲线的影响主要体现在溶洞向裂缝窜流后的流动阶段，而对早期径向流和早期线性流阶段只产生微小的影响。从溶洞-裂缝拟径向流阶段开始，压力及其导数曲线随着时间的增加而逐渐上翘，且上翘幅度随着分形指数的增大而加剧。分形指数并不影响各个窜流段出现的时间以及窜流强度，所以压力导数曲线上下凹段的深度和持续时间均保持不变。

图5 与图6 分别为水平井长度LD和纵向位置zwD对压力动态曲线的影响，各基本参数取值如下:CD=100，cD=0，β=0.1，S=0，ωf=0.01，ωv=0.1，λm=0.001，λv=0.1，hD=200，αKD=0.02。研究表明，水平井越长，压力及其导数曲线越低，早期径向流越明显，从而后续流动阶段出现的时间延迟。当水平井长度小到一定程度时，早期径向流易被井筒储集效应所掩盖。考虑储层分形和应力敏感性特征，不同长度水平井的压力动态曲线并不在后期汇聚在一起，而是呈平行分布。在其他参数不变的情况下，水平井越靠近油层中部，即zwD的值越接近0.5，早期径向流阶段出现得越早，同时早期线性流持续时间越长，但水平井纵向上位置的变化并不影响其他流动阶段。

图5 水平井长度对压力动态曲线的影响Fig.5 Effect of horizontal well length on dynamic pressure curves

图6 水平井纵向位置对压力动态曲线的影响Fig.6 Effect of horizontal well position on dynamic pressure curves

5 应用实例分析

为了验证本文模型的正确性，进一步认识具有分形特征的缝洞型碳酸盐岩油藏的渗流特征，对四川CN 油田S63 井的渗流参数进行了解释。为了便于试井解释，假定裂缝垂向渗透率与水平渗透率相等，即β=1。S63井的基础参数取值如表2 所列。

各个模型的异同已在渗流规律及敏感性分析部分进行了详细对比，但为方便对比、分析本文分形非线性模型的合理性，同时应用达西模型和本文模型进行试井拟合（图7）。2 个模型的拟合结果均显示压力导数与压力曲线分开。实际压力导数曲线和理论曲线在前期拟合程度相对较低，这是由于实际储层的井筒储集系数发生变化所导致的，在实际试井过程中普遍存在类似的问题。在渗流的中后期，虽然实测数据不完整，压力导数曲线仍明显出现反映孔洞缝三重介质储层渗流特征的2 个下凹段，且相对于达西模型，本文模型的拟合结果更加精确，理论计算值和实测值拟合得更好。选取本文的分形非线性模型进行S63 井的试井解释，表3为该模型所解释的渗流参数。

表2 S63 井基础参数取值Table 2 Values of basic parameters of well S63

实际生产曲线的基质向裂缝的窜流段明显比达西模型所呈现的下凹程度小，且试井解释下的非线性参数为0.15，基质渗透率为1.52 mD，这都体现了基质确实存在非达西渗流。拟合结果表明，采用本文模型进行渗流参数解释更加精确，这说明储层考虑分形特征和非线性渗流后更符合油井的实际情况，可应用本文的分形非线性模型进行实际试井解释及渗流规律分析。

图7 实际算例压力动态拟合曲线Fig.7 Dynamic pressure fitting curves of the actual example

表3 S63 井渗流参数解释Table 3 Explanation of seepage parameters of well S63

6 结论

（1）引入分形理论描述裂缝系统的分形特征，同时应用非线性渗流新模型描述致密基质块的渗流特征，在此基础上建立缝洞型油藏水平井分形非线性渗流模型，并基于有限元方法求解水平井动态压力。

（2）通过对比达西模型与分形非线性模型的异同点，总结缝洞型碳酸盐岩油藏水平井的渗流规律，分析各个流动阶段的特点，将渗流过程划分为9 个流动阶段，即井筒储存、表皮效应过渡、早期径向流、早期线性流、中期径向流、溶洞向裂缝窜流、溶洞-裂缝拟径向流、基质向裂缝窜流以及系统拟径向流。

（3）非线性参数主要影响基质向裂缝系统的窜流段，随着非线性参数的增大，基质向裂缝系统窜流强度减小，持续时间变短。分形指数不影响窜流出现的时间及强度，它对压力动态曲线的影响主要体现在溶洞向裂缝窜流后的阶段，压力及导数曲线随时间增加而逐渐上翘，且上翘幅度随分形指数的增大而加剧。

（4）利用本文的分形非线性模型进行实际井动态压力拟合并解释相关渗流参数，考虑了储层分形特征和非线性渗流后，模型更贴近矿场实际，可用其进行实际试井解释及渗流规律分析。