三维压裂缝网不稳定压力半解析求解方法1)

王志凯 程林松,2) 曹仁义,3) 王 进 贾 品 王选茹

* (中国石油大学(北京)石油工程学院,北京 102249)

† (中国石油天然气股份有限公司长庆油田分公司,西安 710000)

引言

随着常规油气资源储量不断降低,全球非常规油气资源勘探开发进入活跃期,致密油、页岩油等非常规油气资源探明储量不断上升[1-3].特别是近5年,中国非常规油气进入了规模发现与开发上产阶段,致密、页岩油气展现出了巨大的资源潜力[4-7].与常规储层不同,致密储层多为源储一体,储集层致密,未经压裂储层难以形成工业油流.开发中需采用体积压裂技术将储层“打碎”,形成“人工油气藏”[8,9].通常认为压裂缝沿最大主应力方向延伸[10-12],但受原始地应力及压裂工艺影响,压裂阶段储层应力状态非常复杂,现场微地震数据及耦合应力场压裂模拟结果表明,体积压裂后储层缝网展布复杂,存在不同的倾角、形态及穿透比[13-16].

对于上述储层,其渗流规律较为复杂,单一线性流难以表征其流动状态[17-19],樊冬艳等[20]、Al Rbeawi 等[21-22]基于点源函数叠加原理研究了相同倾斜方向下多段压裂无限导流倾斜缝典型试井曲线特征.贾品等[23-25]、Teng 和Li[26]、万义钊等[27]通过将裂缝离散为微元的方法实现不同裂缝形态的表征,同时Teng 和Li[28-29]、王苏冉[30]通过引入空间倾斜面源实现单条倾斜缝压力、产量瞬态分析.但上述研究所考虑的裂缝只能沿单一方向倾斜,在此条件下,可通过调整坐标系的方法降低模型复杂度实现模型有效求解.但唐慧莹等[14]研究表明,即使在常规多段压裂过程中,压裂缝也会因缝间应力干扰发生裂缝面弯曲,出现裂缝倾斜方向不同的现象.因此为更准确表征水力压裂缝空间展布,建立可表征任意方向倾斜、任意形态展布的多段压裂井压力瞬态分析模型[31-34]是非常必要的.

本文通过离散裂缝网络表征不同裂缝形态,结合点源函数面积分与裂缝微元局部坐标系,构建了不受空间坐标系约束的封闭盒式储层任意平面缝源函数,并验证了其准确性.针对三重积分面源积分函数存在的计算效率较低的问题,提出了应用点源或特殊线源代替面源的方法实现了倾斜裂缝快速求解,探讨了此类求解方法的准确性及可行性,基于此模型划分了单段、多段倾斜压裂缝流态,分析了裂缝导流能力、裂缝倾角、裂缝高度以及裂缝段间距等因素的影响.

1 多段压裂井三维缝不稳定渗流模型

1.1 物理模型及模型假设

斜缝,裂缝半缝长为lf/2,m;沿裂缝面延伸高度为hf,m;裂缝宽度为w,m;各裂缝倾斜方向不同,每条裂缝内采用局部坐标 (eεI,eξI)进行定位.

如图1 所示,均质封闭储层 (xe,ye,h) 中存在一口体积压裂大斜度井,沿井筒延伸方向存在M条倾储层及裂缝孔隙度为 φm,φf;渗透率为km,kf,mD;综合压缩系数为Ctm,Ctf,MPa−1.流体黏度为μ,mPa·s;体积系数为B;油井定产qw生产,m3/d;井储系数为C,m3/MPa;不考虑表皮影响.

图1 多段压裂井三维缝物理模型Fig.1 Physical model of multi-fractured well with 3D fractures

模型基本假设如下:

(1) 将渗流阶段划分为基质向裂缝渗流及裂缝向井筒流动两部分,不考虑基质向井筒的直接渗流;

(2) 不考虑井筒沿程压力降,忽略渗流过程中的重力影响;

(3) 裂缝有限导流,单条裂缝内孔渗数据保持一致;

(4) 储层内流体为油单相,无毛管力影响.

1.2 裂缝内部流动模型

水力裂缝形态受储层物性及压裂工艺影响[35],常见裂缝形态有矩形缝、椭圆缝、菱形缝等.本文将裂缝离散为矩形微元,实现不同形态裂缝有效表征.如图2 所示,井筒穿过椭圆缝中心,结合椭圆裂缝形态将其离散为17 个矩形微元,其中中心点微元(编号9)为井筒所在微元.

图2 椭圆缝离散表征Fig.2 Discrete characterization of elliptical fracture

模型推导与计算过程中,裂缝微元相关参数定义如下:第I条裂缝微元总量为NfI,其中井筒所在微元编号NfwI;裂缝系统微元总量为Nf;第i个微元的长、宽为 Δ εi, Δ ξi,m;定义每个裂缝微元的渗透率为Ki,相应的可确定裂缝微元有限导流能力为Kiwi,通过给定不同的Ki值即可实现裂缝不同有限导流能力求解.

基于裂缝内二维不稳定渗流,建立裂缝内流动综合控制方程.

考虑井筒与椭圆缝相交于微元9,以该微元为研究对象,采用有限差分方法将公式(1)离散处理

式中, β为单位转换系数,取值0.0864;为第n时间步基质向第i个微元的供给,m3/d;为第n时间步第i条裂缝向井筒的 供给,m3/d;为第n时间步第i个微元的压力,MPa;Δd(i,j)为微元i,j中心点连线的距离,m.

为简化模型推导过程,引入无量纲参数

将无量纲参数代入式(2)

化简可得

式中,TD(i,j)为i,j两裂缝微元间的无因次传导率,αDi为相邻时间步间的影响.

将式(4)应用到离散裂缝系统所有微元上,可得包含Nf个裂缝微元的流动方程有限差分方程组

式中,矩阵T大小为Nf×Nf,T(i,j) =TD(i,j),表征微元间不稳定渗流影响;INf为Nf×Nf的单位矩阵,表征基质源汇项影响;矩阵0 大小为Nf×M的零矩阵;矩阵a大小为Nf×M,表征井筒源汇项影响;矩阵R1大小为Nf× 1,由第n−1 时间步参数决定

1.3 油藏流动模型

考虑上述离散方法所得第I个微元中心点坐标为 (xDI,yDI,zDI),其局部坐标系由矩形微元相邻边所在方向的单位向量 (eεI,eξI)构成,基于中心点坐标及局部坐标系可确定微元内任意点坐标

结合Newman 乘积与封闭边界无限大平面源经典解[36]可得封闭储层(x,y,z) 处单一点源以定产q生产时在点(x0,y0,z0)处产生压力降

将无量纲参数代入式(14)

由叠加原理可知,空间面源可由空间点源沿微元面积分获得,因此空间面源函数表达式

基于空间面源压力分布,结合压力叠加原理及杜哈美原理得第i个微元第n时间步在第j个微元处产生的压力降

为简化方程形式,定义式(16)中的积分项为参数

从而

基于式(19)建立基质向裂缝流动矩阵

式中,矩阵B大小为Nf× 2(Nf+M),表征基质向各微元流动的影响;矩阵R2大小为Nf× 1,由n− 1 时间步中压力及流量参数构成参数求解过程中计算效率较低,因此考虑以下

由式(19)可知,随着Nf及n的增大,式(18)的调用次数迅速升高,而式(18)中的三重积分会导致两种求解方法.

(1) 线源代替面源

当多条裂缝法向量共面时,可建立坐标轴(x,y,z)使其中一条坐标轴与裂缝微元局部坐标系的一条向量平行,实现模型简化,如图3 所示.

图3 裂缝微元局部坐标系示意图Fig.3 Schematic diagram of local coordinate system of fracture elements

以局部坐标轴eε与坐标轴x平行为例,此时xξ,yε,zε均为0,式(12)中yD(ε,ξ),zD(ε,ξ)将转化成ξ的函数

相应的参数可简化为双重积分形式

考虑参数由线源沿微元斜边积分所得,通过调整ξ 方向上微元数量使单个微元内该方向上线源函数变化较小,引入积分中值定理,以微元中心点(xDI,yDI,zDI)处线源与Δξ 之积近似代替源函数积分实现式(24)的进一步简化,获得参数的一重积分形式

(2) 点源代替面源

参照方案一的思路,直接对三重积分形式下参数使用积分中值定理,以微元中心点(xDI,yDI,zDI)处点源与微元面积Δε × Δξ 之积近似代替源函数积分实现式(18)的简化,该条件下需保证裂缝微元的尺寸在一定范围内,后面将对其精度进行讨论

将上述方案中的代入式(21)可实现简化求解 条件下油藏流动模型求解矩阵的构建.

1.4 模型耦合及求解

由Peaceman 公式[37]可知井筒所在微元内边界条件可由稳态径向流公式处理

基于Hurst 等人[38]的研究,在模型中考虑井储系数CwD的影响.

同时,井筒无限导流条件下,各裂缝处井筒压力相等,即

本研究选择武汉市洪山区光谷地区部分城乡结合部2006年、2009年、2011年3年的遥感影像(3期遥感影像购买自中国遥感数据网，都为6月份的数据)作为初始数据，空间分辨率分别为2.4 m、2.4 m和2.0 m.为了满足本研究需要，把研究区用地简单分为城市用地和非城市用地，具体详见表1，土地利用类型分类参考参考文献[1]中的分类标准.

结合式(27)−(29)建立与井筒内流动相关渗流方程组

式中,矩阵C大小为2M× 2(Nf+M),表征了式(27)−(29)中第n时间步压力及产量的系数项,矩阵R3大小为Nf× 1,由n−1 时间步相关参数构成.

该方程组中存在2(Nf+M)个方程,其中未知量包含Nf个裂缝微元的压力pfD,Nf个裂缝微元的产量qfD,M条裂缝处的井底压力pfwD以及M条裂缝向井筒的供给qfwD共计2(Nf+M)个,未知量数与方程数相等 ,方程组可封闭求解.

2 模型验证分析

2.1 三重积分模型验证

本模型中,考虑多段压裂大斜度井压裂缝为有限导流,且裂缝面沿任意方向展布,现有文献相关研究较少.为验证模型准确性,分别验证有限导流单一倾斜缝与无限导流多段压裂缝求解的准确性.针对图4 所示部分贯穿倾斜缝(partially penetrating inclined fracture,PPIF),Teng 和Li[28]通过建立坐标轴1 (Coordinate 1)实现模型点源函数求解简化.为验证本文模型三重积分求解的准确性,将坐标轴1 沿z 轴旋转45°建立新坐标轴2,增加模型复杂度.

图4 坐标轴转换示意图Fig.4 Schematic diagram of coordinate axis transformation

模型数据与Teng 图4 数据相同,具体无因次参数如下:xeD=yeD= 20,hD= 0.5,θ= 45°,hfD=CfD= 5,Cs= 10,CwD= 0,wD= 1.0 × 10−5,γ= 1.3 × 10−8.

坐标轴2 条件下微元内局部坐标系如下

图5 本文模型与Teng 模型[28]结果对比Fig.5 Comparison of the result of our model with Teng model[28]

为研究多缝条件下模型求解准确性,选取Al Rbeawi 和Tiab[22]2013年提出的无限大地层多段压裂倾斜缝图版(A-1)进行对比验证,选取裂缝倾角为15°,hD= 1,无因次段间距DD= 1.如图6 所示.

图6 本文模型与Al Rbeawi 模型结果对比Fig.6 Comparison of the result of our model with Al Rbeawi model

可以看出,不同时间下两模型计算结果基本吻合,由于本模型采用封闭边界,当tD>10 压力响应到达封闭边界后的压力及压力导数变化与Al Rbeawi模型不同由于该模型采用三重积分进行源函数计算,求解时间较长(数小时),适用性较差,考虑对点源函数求解进行进一步调整.

2.2 源函数调整准确性验证

线源、点源简化模型求解主要误差来自于积分中值的选取,当裂缝微元沿积分方向的长度足够小时,点源函数沿积分方向上单调变化,此时选取积分中点作为积分中值可取得较好的积分效果.图7 为不同网格划分下线源求解所得曲线变化对比.可以看出,当网格划分较为粗糙时,线源法所得初期压力及压力导数曲线误差较大.特别的压力导数曲线存在较为明显的回流特征,考虑是由于中心线所在线源与积分中值相差较大,各线源之间存在流向各微元中心线源的径向流动所致(图8).

图7 网格划分对线源求解效果影响Fig.7 Effect of meshing on line source solution

图8 回流段流动示意图Fig.8 Diagram of flow in flowback period

随着划分微元数目增多,简化模型所得压力及压力导数曲线与面源积分所得压力及压力导数曲线趋于一致.大量模拟表明当所划分微元无因次边长在0.15 左右时,采用点源或线源代替面源的方法可获得较高的精度.选取的单缝,划分微元数目15 × 3,模拟获得相应典型试井曲线如图9所示.可以看出在tD大于10−7时间内采用点源或线源函数代替面源积分可取得较高计算精度,相同模型下,面源积分计算时长为435 443 s,而线源积分计算时长仅为17 s,点源积分计算时长仅为22 s,极大地提升了模型计算效率.

图9 简化模型求解效果示意图Fig.9 Diagram of simplified model

3 典型曲线分析

3.1 流动阶段划分

通过分析典型试井曲线中的压力及压力导数曲线特征,可实现不同流动阶段识别.选取无因次参数xeD=yeD= 20,hD= 0.5,θ= 45°,hfD=50,Cs= 10,CwD= 5,wD= 1.0 × 10−5,γ= 1.3 × 10−8.绘制考虑井储的压力及压力导数曲线.从图10 可以看出,单一倾斜缝典型试井曲线可划分为8 个阶段.

图10 单一倾斜缝流动阶段划分Fig.10 Identification of the flow regimes of incline fracture

(1)井筒储集及过渡阶段:井储阶段压力及压力导数曲线重合,均为斜率等于1 的直线,该阶段主要流体供给来自井筒;

(2)裂缝线性流阶段:压力导数曲线斜率为1/2,该阶段基质暂未动用,与井相连的裂缝内部流体近线性流入井筒;

(3)双线性流阶段:压力导数曲线斜率为1/4,该阶段除裂缝内部流动外,近缝基质向裂缝的线性流动开始出现;

(4)基质线性流阶段:压力导数曲线斜率为1/2,该阶段裂缝内压力趋于稳定,近缝基质向裂缝的流动仍以线性流为主;

(5)早期径向流阶段:压力导数曲线斜率为0,此时波及范围扩大,基质向裂缝的供给由线性流转变为径向流;

(6)椭圆流阶段:压力导数曲线斜率近似为0.36,通常认为该阶段为缝周基质径向流向裂缝系统径向流过渡的阶段;

(7)晚期径向流阶段:压力导数为斜率等于0 的直线,流体围绕裂缝整体成径向流;

(8)边界控制流阶段:流体流动到封闭边界,压力及压力导数响应主要受边界影响,压力及压力导数均为斜率等于1 的直线.

为研究裂缝条数对流动阶段影响,取缝间距DD=1,对比裂缝条数分别为1,2,3 条件下压力及压力导数曲线(图11)可以看出,随着裂缝条数增加,主要流动阶段仍为8 个,但与单条缝相比存在一定的差别.随着裂缝段数增多,裂缝及近缝区域流动影响增强,井储阶段提前进入过渡阶段,相应的双线性流阶段出现时间略有提前;由于缝间干扰的存在,裂缝条数越多,早期径向流阶段越不太明显,过渡阶段斜率更大,晚期径向流及边界控制流阶段,流体相对于裂缝系统进行流动,不同裂缝条数下压力及压力导数曲线基本重合.

图11 多段压裂倾斜缝流动阶段Fig.11 Flow regimes of multi-stage incline fractures

3.2 影响因素分析

本文所建模型在处理裂缝过程中主要进行了以下几点考虑:多段压裂缝、裂缝有限导流、裂缝以任意形态展布以及裂缝向任意方向倾斜,因此针对上述几点考虑分别进行裂缝导流能力、裂缝倾角、井筒倾角、裂缝纵向高度以及裂缝段间距等参数敏感性分析.

(1)裂缝导流能力敏感性分析

图12 给出了无因次裂缝导流系数分别为0.5,5,50 及500 时压力及压力导数曲线响应,可以看出随着裂缝导流能力的升高,井储阶段及过渡阶段更早结束,早期及晚期径向流阶段时间变短,特别是早期径向流阶段基本消失.

图12 裂缝导流能力敏感性分析Fig.12 Sensitivity analysis of well storage coefficient

(2)裂缝倾角敏感性分析

考虑裂缝纵向高度保持一致,对比裂缝面倾角分别为15,30,45,60 以及75 时压力及压力导数曲线变化,由图13 可以看出,裂缝面倾角主要影响双线性流阶段至椭圆流阶段的压力导数曲线,当裂缝倾角小于45°时,倾角影响较小,随着倾角进一步增大,双线性流作用时间增长,早期径向流更为明显,相应的椭圆流阶段倾角增大.

图13 裂缝倾角敏感性分析Fig.13 Sensitivity analysis of well storage coefficient

(3)裂缝高度敏感性分析

取裂缝倾角为45°,考虑裂缝无因次纵向高度在0.2–1.0 范围内变化.由图14 可以看出,当裂缝高度较小时,裂缝内线性流即双线性流阶段较不明显,且基质线性流后期存在一定的回流现象,随着裂缝高度的升高,随着裂缝高度的增大,各阶段流态逐渐出现.

图14 裂缝高度敏感性分析Fig.14 Sensitivity analysis of fracture height

(4)裂缝段间距敏感性分析

取裂缝倾角为15°,考虑无因次裂缝段间距在1–4 范围内变化,由图15 可以看出,随着裂缝段间距增大,早期径向流至晚期径向流阶段的特征更为不明显,表明随着裂缝段间距的增大,缝间干扰越晚出现且系统径向流出现时间越晚.

图15 段间距敏感性分析Fig.15 Sensitivity analysis of fracturing interval

4 结论

(1)对于任意方向展布裂缝可通过裂缝离散表征,耦合裂缝内数值解与基质内面源函数解析解的方法进行求解,但存在计算效率低的问题.借助点源、特殊线源替代面源的方法可大大提升计算效率,且在裂缝微元划分较为精细(微元无因次边长小于0.15)时精度较高.

(2)倾斜缝典型试井曲线可划分为井筒储集及过渡流、裂缝线性流、双线性流、基质线性流、早期径向流、椭圆流、晚期径向流及边界控制流八个流动阶段.随着裂缝段数增多,缝间影响增强,井储、过渡阶段及径向流阶段特征减弱.

(3)裂缝倾角及裂缝段间距主要影响裂缝线性流到早期径向流阶段压力及压力导数曲线,特别的当裂缝高度较小时压力导数存在回流现象.