多翼裂缝垂直井的瞬态压力分析

2021-12-26张守江乔红军张永飞陈一鸣曾凡华

西南石油大学学报(自然科学版) 2021年5期

张守江，乔红军，张永飞，陈一鸣，曾凡华

1.延长油田股份有限公司勘探开发技术研究中心，陕西延安 716000；2.陕西延长石油（集团）有限责任公司研究院，陕西西安 710065；3.辽宁石油化工大学石油天然气工程学院，辽宁抚顺 113001；4.里贾纳大学研究生院，萨斯喀彻温省里贾纳 SKS4S0A2

引言

水力压裂作为油气井增产、注水井增注的一项重要技术措施，不仅广泛应用于低渗透油气藏，在中高渗透油气藏的增产改造中也取得了很好的效果。在过去20 年中，世界各地非常规油气市场利用该技术均得到了实质性的发展[1-4]。

水力压裂增产增注的原理是通过制造水力裂缝来降低地层内流体的渗流阻力，进而改变流体的渗流状态。如果水力裂缝能连通油气层深处的生产层和天然裂缝，则增产效果会更加明显[5-6]。但水力压裂所产生的水力裂缝与天然裂缝并不是严格正交化的，地层内的裂缝网络通常呈不规则的相间分布和复杂的互联场景。Sierra 等[7]与Molenaar 等[8]利用微震监测、光纤温度传感等技术对压裂增产情况进行了监测，结果表明多处裂缝的生长状态不尽相同，复杂的地应力场和地层非均质性对裂缝的位置和发育情况具有较严重的影响。

为了研究复杂裂缝网络的特性，国内外众多学者对此进行了研究，王晓冬等[9]利用不稳定渗流理论对封闭地层有限导流垂直裂缝井的压力动态特征进行了研究并得到了较为准确的中期径向流公式和晚期拟稳态渗流公式，而后采用数值反演方法给出了无因次导流能力、裂缝壁面表皮以及渗流区域半径等因素的对比图版；Cinco 等[10]建立了有限导流垂直裂缝板状油藏的瞬态数学模型，通过理论推导及实例分析，验证了半对数压力分析方法的准确性，并运用该方法对地层内早期阶段的瞬态压力数据进行了分析并得到了地层特征和裂缝特征；Wang等[11]基于Fick 定律和达西流动规则提出了一种考虑边界效应的半解析模型，并运用该模型对非对称裂缝煤层井的瞬变压力特性进行了分析，通过改变裂缝导流能力、煤地层孔隙度和渗透率等参数，明确了各变量对井底瞬态压力行为造成的影响；任宗孝等[12]在Ozkan[13]平面源函数理论的基础上，通过等效变换和坐标变换推导出了一种新型平面源函数计算方法，并运用该方法建立了复杂裂缝水平井非稳态压力的半解析模型，模型识别了地层内的5 种流动类型以及裂缝数量、裂缝分布、储集比等关键参数对压力变化的影响；Luo 等[14]通过引入裂缝翼与坐标轴夹角概念，利用笛卡尔旋转坐标系方法对裂缝垂直井的瞬态压力特性进行了分析，在此基础上，通过改变裂缝翼交角、裂缝长度、裂缝非对称性等条件，对瞬态压力响应的敏感性因素进行了分析及讨论。

综合分析国内外研究现状可知，鲜有学者针对同一垂直井眼产生的复杂破裂模式对瞬变压力响应的影响进行研究，而对于复杂裂缝网络压力响应特性分析主要通过引入裂缝翼与坐标轴夹角、依靠坐标变化的方法进行研究，该方法对于特性条件下的规则裂缝很适用，但是对于不规则、非对称的裂缝网格计算效果并不是很好。

由于无论是水平井、还是垂直井，都必须对孔隙度、渗透率和压差等因素进行适当的研究，以确保流体有效地向上流动。压差对油井的流量和产量有很大的影响，特别是在垂直井多个非对称裂缝的生产过程中。

基于此，提出一种基于笛卡尔坐标系的坐标计算方法，通过给定地层内裂缝所在的位置信息计算任意段的压力响应特征，简化常规坐标变化的计算过程提高了求解和分析的效率。

1 物理模型及假设

地层水力裂缝通常不单一存在，Germanovich等[15]通过对裂缝岩芯的直接观察发现它们具有不规则的分支，为了便于分析问题，可对实际的复杂裂缝进行一定的简化处理，结果如图1 所示。

图1 水力复杂裂缝示意图Fig.1 Schematic diagram of complex hydraulic fractures

无限大平板油藏中有限导流垂直裂缝井的物理模型如图2 所示。

图2 无限大平板油藏中有限导流垂直裂缝示意图Fig.2 Schematic diagram of finite diversion vertical fractures in an infinitely large flat reservoir

模型需要满足下列假设条件：

（1）假设地层上下为不渗透边界，水平为无限大边界的均质油藏被压裂开一条垂直裂缝，裂缝关于井筒对称，裂缝半长为xf，裂缝宽度为wf。

（2）假设裂缝具有一定的渗透率Kf，沿着裂缝存在不可忽略的压降损失。

（3）忽略重力和毛细管压力的影响。

（4）地层流体流动过程符合达西渗流准则。

2 数学模型

2.1 数学推导

根据叠加理论，将油藏地层系统分为地层-裂缝子系统和裂缝子系统。

对于地层-裂缝子系统，将每个裂缝都离散成Ni段，并且对每段性质进行单独定义。裂缝离散示意图如图3 所示。

图3 裂缝离散示意图Fig.3 Schematic diagram of fracture dispersion

对于裂缝子系统，其内部的流体流动假定为线性流动，离散后的裂缝段利用Cinco 等[16]所提出的线源解决方案进行求解。

首先，基于流量和压力的连续性条件，将流体在两个子系统中的流动耦合到每个单元的界面上；然后，在Laplace 域对裂缝进行离散和求解，并应用Stehfest 反演算法[17]确定压力分布。为了便于该模型的求解，无因次压力和时间定义为

式中：

pD—无因次压力；

K—地层渗透率，mD；

h—地层厚度，m；

pi—原始地层压力，Pa；

p—压力，Pa；

Q—流量，m3·s−1；

B—流体体积系数，无因次；

µ—流体黏度，Pa·s；

tD—无因次时间；

t—时间，s；

xf—裂缝半长，m；

ct—综合压缩系数，Pa−1；

φ—地层孔隙度，%。

为了考虑地层内气体流动，引入拟压力

式中：pa—拟压力，Pa；

pr—参考压力，Pa；

Z—偏差系数，无因次。

求解时，无因次坐标定义为

式中：xD、yD—裂缝无因次坐标；

x、y—裂缝实际坐标，m。

地层内部裂缝的无因次导流能力和裂缝扩散率定义为

式中：CfD—无因次裂缝导流能力；

Cη—裂缝扩散系数，无因次；

Kf—裂缝渗透率，mD；

ctf—裂缝压缩率，%；

φf—裂缝孔隙率，%；

wf—裂缝宽度，m。

2.1.1 地层-裂缝子系统

极坐标下用于描述地层流体流动的方程为

式中：qfD—无因次裂缝流量。

当考虑井筒储集效应时，近井地带内边界条件应受到限制

式中：CD—井筒储存系数，无因次；

rD—源点的无因次极径。

地层子系统解可通过Laplace 域中的线性源解获得

u—Laplace 变量；

K0—修正贝塞尔函数；

2.1.2 裂缝子系统

研究表明[9，16，18-20]，若裂缝无因次导流能力等于或大于300，裂缝具有无限导流能力的假设是有效的。

因此，假定每个水力裂缝均是有限导流裂缝，裂缝每段的扩散系数方程可以描述为

式中：pfD—无因次裂缝内流体压力；

下标i—第i个离散裂缝单元。

通过在Laplace 域求解式（10），可得

基于Van Kruysdijk[21]的结果，xD处的无因次压力在拉普拉斯域可以简化为

式（12）中，无因次压力由管段边界流量和地层流量决定，系数bi，ci由裂缝扩散率决定。

2.2 离散化与耦合

为了计算简便，在大多数解析方法和半解析方法的研究中会假设裂缝的断裂模式为正交拓扑状态[22]，但实际的水力裂缝网络通常是相互交错、复杂多变的。为了解决裂缝非正交的问题，许多研究人员，如Luo 等[14，23-24]引入了表示非正交拓扑中裂缝翼和坐标轴的交角θ；虽然该方法可以解决一般非正交裂缝的计算问题，但与此同时增加了方程求解的难度，并且该方法对于不规则复杂裂缝的计算和分析会产生一定的误差。为了解决上述问题，提出了一种新的离散耦合方法，通过对地层内的裂缝坐标计算方法进行优化设计，直接利用裂缝所在位置的信息进行计算。

为了说明本文所提出的离散耦合方法，以3条规则交叉裂缝网络和两条交叉非规则裂缝网络为例，裂缝分布及离散结果如图4 所示。图4a 和图4b为裂缝在二维平面上的规则曲线，每个裂缝均被分成两部分。例如，裂缝a可分为和。然后将每个裂缝部分进一步离散成Ni段。图4c 与图4d为裂缝在二维平面上的不规则曲线，这些不规则曲线可以用近似直线代替。例如，曲线e可以被分割并接近直线裂缝的和。然后将每裂缝片段进一步离散成Ni段。更为复杂的裂缝形态或裂缝几何结构可以在不考虑两条裂缝夹角的情况下进行处理。

图4 裂缝离散网络示意图Fig.4 Schematic diagram of irregular fracture network

由于裂缝每段的各部分都可以看作地层中的源（汇），因此，任何两源（汇）之间的干扰只取决于其在二维平面上的距离。相邻两段与段中心的界面耦合是组装代数矩阵的基本步骤，需遵循3 条原则。

（1）采用叠加原理计算节段中心压力。

（3）通量的连续性可以表示为

最后，采用Stehfest 算法[17]将Laplace 域的流量和压力分布转化为真实时域。

3 模型验证

利用试井软件对模型进行了验证。计算模型为单一垂直裂缝井，不考虑井筒储存效应的影响，选择两种不同裂缝导流能力的实例说明所提出方法的有效性和准确性。图5 显示了两个方案的比较结果。

由图5 可知，当裂缝无因次导流能力相对较高（如Cfd=750）时，模型与商业软件KAPPA 的计算结果相一致，无因次压力和压力导数的双对数诊断曲线均出现了斜率为0.5 的地层线性流动和数值为0.5 的拟线性流动[10]；当裂缝无因次导流能力相对较低（如Cfd=3.75）时，模型与商业软件KAPPA 的计算结果也吻合较好，无因次压力和压力导数的双对数诊断曲线均出现了斜率为0.25 的双线性流动和数值为0.50 的拟线性流动。

图5 本文提出的模型与KAPPA 提出的模型对比Fig.5 Comparison between the model proposed in this paper and the model proposed by KAPPA

对比分析可知，无论高导流能力还是低导流能力的裂缝，模型均可准确地解释地层内的流动阶段，因此，利用该模型讨论复杂裂缝井的压力瞬态行为是准确可靠的。

4 结果与讨论

为了分析无限大地层有限导流垂直井多翼裂缝的压力瞬变特性，对于油藏地层本身，考虑地层渗透率、地层孔隙度、地层厚度、流体黏度、产量等参数；对于裂缝系统，关注裂缝数量、分布形式、非对称因子、裂缝长度和裂缝导流能力等参数。

4.1 基础数据

为了分析多翼非对称裂缝对垂直井内瞬态压力响应特征及变化规律，设置检查裂缝性垂直井中瞬态压力行为的基本数据如表1 所示。

表1 基础数据表Tab.1 Basic data table

4.2 敏感性分析

4.2.1 裂缝数量

为了研究裂缝数量对地层内瞬态压力特性的影响。选取4 种不同的裂缝数量进行对比分析，其中，裂缝长度均为100 m、导流能力均为200 mD·m；通过代入每个裂缝在平面内的坐标，进而求解得到不同裂缝数量的瞬态压力响应特征。裂缝断裂数设计方案如图6 所示。

图6 裂缝断裂数设计方案Fig.6 Design scheme of fracture and fracture number

不同裂缝数量的瞬态压力响应特征如图7 所示。由图7 可知，在裂缝导流能力均为200 mD·m时，不同裂缝数量垂直井的瞬态压力响应特征及流动阶段近似相同，均包括早期阶段的井筒储集和裂缝双线性流阶段、中期的地层线性流阶和后期的拟径向流阶。

在井筒储集与双线性流两个阶段，裂缝数量对垂直井内瞬态压力响应特性的影响较小，这是由于流动初期并未有地层流体流入裂缝，因此，裂缝的承载阻力相对较小；井筒储集阶段的瞬态压力响应特征具体表现为无因次压力与压力导数双对数曲线上均呈斜率为1.00 的直线，双线性流阶段的瞬态压力响应特征具体表现为无因次压力与压力导数双对数曲线上均呈斜率为0.25 的平行直线。

在地层线性流阶段，裂缝数量对垂直井内瞬态压力响应的影响较大；随着裂缝数量的增加，无因次压力与压力导数在双对数曲线上均出现向下平行移动的现象，即裂缝数量越多，则流体进入到地层线性流阶段的起始时间就越滞后且流动周期越短。在拟径向流阶段，裂缝数量对垂直井内压力效应的影响较大；随着裂缝数量的增加，无因次压力双对数曲线发生了下移现象，而无因次压力导数双对数曲线发生了右移现象，即裂缝数量的增加使得井内拟径向流的发生时间出现延后现象，该阶段瞬态压力响应曲线的无因次压力导数曲线均呈值为0.5 的水平直线。

4.2.2 裂缝分布

为了研究裂缝分布方式对地层内压力变化的影响，选取4 种情况进行对比分析。保证地层内的裂缝总长度始终为400 m，随着裂缝数量的增加，原始裂缝的长度在逐渐减小，例如，当地层内只有一条裂缝时，其长度为400 m，当地层内裂缝数量增加为两条时，其每条裂缝的长度为200 m，以此类推。裂缝分布的设计方案如图8 所示。

图8 裂缝分布方式设计方案Fig.8 Design scheme of fracture distribution

不同分布方式裂缝的瞬态压力响应特性如图9所示。可以看出，地层内不同裂缝分布形式的瞬态压力特征曲线与裂缝数量的敏感性分析结果相似，包括井筒储集阶段、裂缝双线性流阶段、地层线性流阶段和拟径向流阶段。

图9 不同裂缝分布方式的瞬态压力诊断曲线Fig.9 Transient pressure diagnosis curve of different fracture distribution modes

虽然不同裂缝分布方式的压力特征曲线与裂缝数量近似相同，但是在每个流动阶段的内部却又存在明显的差异，尤其是双线性流和地层线性流阶段。通过分析双线性流阶段的压力诊断曲线可知，模型1 的双线性流动特性最不明显，而其他3 种情况的双线性流动特征较为明显且远离了模型1；通过对比可知，该阶段模型3 与模型4 的油藏采集效率远高于模型1 与模型2，这是由于模型3 与模型4的裂缝长度虽然较小，但是裂缝的控制范围明显增加，因此，地层内的线性流动与裂缝内的线性流动更加强烈。通过分析地层线性流阶段的压力诊断曲线可知，该阶段的压力响应特征正好与双线性流阶段相反，该阶段模型3 与模型4 油藏采集效率远低于模型1 与模型2，这是由于在地层线性流动阶段，大部分地层流体会沿着裂缝逐渐的流入裂缝系统内，随着裂缝数量的增加，相邻裂缝之间的间距将会减小，裂缝之间的相互影响会加强。

4.2.3 裂缝导流能力

为了研究裂缝导流能力对地层内瞬态压力变化特性的影响，选取具有3 条多翼裂缝进行分析，其中，裂缝的半长均为100 m。裂缝分布形式如图10所示，裂缝导流能力设置情况如表2 所示。

图10 裂缝导流设计方案Fig.10 Design scheme of fracture diversion

表2 裂缝导流能力设置表Tab.2 Setting of fracture conductivity

不同导流能力裂缝的瞬态压力响应特性如图11所示。

由图11 可知，不同导流能力裂缝的瞬态压力响应曲线也经历了井筒储集阶段、裂缝双线性流阶段、地层线性流及拟径向流阶段等4 个流动阶段。通过对比分析双对数压力响应曲线可知，4 种导流能力裂缝瞬态压力响应曲线的井筒储集、地层线性流及拟径向流阶段基本重合，其中，井筒储集阶段发生的时间较短，即井筒储集阶段并不是很明显；裂缝导流能力的改变对双线性流阶段，也就是裂缝早期的双线性流动阶段影响较为明显，随着裂缝导流能力的增加，裂缝双线性流动阶段的无因次压力及压力导数诊断曲线均发生了下移现象，根据相关学者的研究可知，无因次压力及压力导数诊断曲线下移说明压裂支撑效果逐渐增强。

图11 不同裂缝导流能力的瞬态压力诊断曲线Fig.11 Transient pressure diagnostic curve of different fracture conductivity

当裂缝导流能力为100 mD·m 时，无因次压力导数诊断曲线出现了明显的波动现象，而当裂缝导流能力为200 mD·m 或者更大时，该现象逐渐消失，这是由于裂缝导流能力越高，地层内的流体越容易、越快地通过裂缝，即当裂缝导流能力越大时，裂缝之间干扰作用越低。

4.2.4 裂缝长度

为了探究地层内裂缝长度对瞬态压力变化特性的影响，选取4 种裂缝长度模型进行对比分析，除了横向裂缝的长度按200 m 递增，其余裂缝长度均为100 m，裂缝的导流能力均为200 mD·m，裂缝长度的设计方案见图12。

图12 裂缝长度设计方案示意图Fig.12 Design scheme of fracture length

不同长度裂缝的瞬态压力响应特性如图13 所示。可以看出，4 个不同横向裂缝长度模型的瞬态压力响应曲线都经历了井筒储集阶段、裂缝双线性流阶段、地层线性流阶段及拟径向流阶段等4 个流动阶段。

图13 不同裂缝长度的瞬态压力诊断曲线Fig.13 Transient pressure diagnostic curves of fracture lengths

观察4 种裂缝长度下的无因次压力诊断曲线可知，裂缝长度的变化对地层线性流阶段影响较大，对其他3 个流动阶段影响较小，在该阶段，随着裂缝长度的增加，裂缝的无因次压力诊断曲线发生了明显的下移现象，即相同的条件下，模型4（横向裂缝800 m）的压降最大，更有利于地层内资源的开采，其次是模型3（横向裂缝600 m）和模型2（横向裂缝400 m），而模型1（横向裂缝200 m）的压降最小、效率最低。

通过观察裂缝内各个流动阶段的发生时间和周期可知，随着裂缝长度的增加，模型2、模型3、模型4的无因次压力诊断曲线逐渐偏离了模型1 的无因次压力诊断曲线，这是由于在地层裂缝内的物质具有一定的膨胀作用，裂缝长度越大，则膨胀作用就越强烈，因此，压力响应曲线偏离原始曲线的程度就越大，但是当裂缝长度增加到一定的程度时，裂缝与裂缝对裂缝的影响作用逐渐减弱，因此，无因次压力诊断曲线偏移的现象逐渐变得不明显。

4.2.5 裂缝断裂角

为了研究裂缝的断裂角度对地层内瞬态压力特性的影响，以两条裂缝垂直井为例，设计4 种断裂角度进行对比分析，设计方案如图14 所示。

图14 断裂角度设计方案Fig.14 Fracture angle design plan

模型裂缝长度均为100 m、导流能力均为200 mD·m 时，不同断裂角度的瞬态压力曲线如图15 所示。可以看出，地层内不同断裂角度的瞬态压力响应曲线近似重合，并且均可划分为井筒储集阶段、裂缝双线性流阶段、地层线性流阶段和拟径向流阶段等4 个阶段，其中，井筒储集阶段的无因次压力曲线与无因次压力导数曲线均呈现斜率为1.00 的直线段；裂缝双线性流阶段的裂缝内无因次压力与压力导数曲线均呈现斜率为0.25 的平行直线，其纵坐标相差对数周期（0.602）；地层线性流阶段的地层内无因次压力曲线与无因次压力导数曲线均呈现斜率为0.50 的平行直线，其纵坐标相差约对数周期（0.301）；拟径向流动阶段与无限导流垂直裂缝相类似，其无因次压力导数双对数曲线呈值为0.50 的水平直线。