基于XFEM-MBEM 的嵌入式离散裂缝模型流固耦合数值模拟方法1)

杜旭林 程林松 牛烺昱 方思冬 曹仁义 ,3)

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

† (中国石化石油勘探开发研究院,北京 100083)

引言

非常规裂缝性储层中,岩石基质提供油气储存空间,人工裂缝、诱导裂缝和天然裂缝提供其流动通道[1].目前油田实际历史数据已表明,地质力学特性对裂缝性油藏的产量起重要作用,会导致一些关键渗流参数发生改变,如基岩的孔隙度和渗透率;此外,致密油藏以“缝控储量”为主[2],岩石形变会引起裂缝开度变化,这对于压裂水平井产能评价的影响极大.传统的应力敏感模型无法揭示深层次的科学道理[3-4],为了准确模拟此类裂缝性多孔介质中的裂缝动态行为和流体流动机理,需要综合考虑流体和岩石两方面的力学性质.

岩石力学中有效应力的概念最早由Terzaghi[5]和Biot[6]提出,并由Biot[7]发展为较完善的三维固结理论,这是研究应力场与渗流场耦合的早期基础.为了将上述的理论基础应用于传统的渗流模型,Barenblatt和Zheltov[8]首次提出了双重孔隙度/双重渗透率的概念,基岩与裂缝为两种重叠的、不同的连续介质,其中基岩完全被裂缝所围绕,用于表征裂缝性多孔介质中的流体和岩石形变.后续,Bai[9]提出了关于双重介质更为严格的物理解释即双孔有效应力理论,用于分析不同尺度的裂缝和基质的形变.针对多孔介质流固耦合数学模型的建立,Wilson 和Aifantis[10]、Beskos 和Aifantis[11]和Khaled 等[12]建立了基于有限元法的双重介质流固耦合模型,该模型的主要问题在于其忽视了基质与裂缝之间压差导致的耦合流动作用;Lewis 和Schrefler[13]对于油藏中的流固耦合问题进行了研究,分别考虑了单相流、多相流、线弹性变形以及弹塑性变形,并给出了相应的有限元计算格式.在上述早期的研究中,采用连续介质力学理论表征裂缝性油藏,将不连续的裂缝系统粗化为连续介质,能在一定程度上刻画基岩和裂缝的流体流动及力学特征,但对于裂缝的几何形状及传导率具有较强的约束,更适用于描述小尺度裂缝,对于大尺度裂缝处理精度较低,无法满足以体积压裂为开发手段的非常规油气藏的矿场实际需求.

近年来,基于离散裂缝的流固耦合理论受到了广泛关注,此类模型将裂缝离散化显式表征,能准确刻画大尺度裂缝对渗流场和应力场的影响,根据前处理算法的不同,可划分为离散裂缝模型[14](DFM)与嵌入式离散裂缝模型[15-16](EDFM),其中嵌入式裂缝模型最早由Li 和Lee[17]开发,旨在提高天然裂缝建模的精确程度.与使用复杂非结构网格在空间上匹配裂缝几何结构的DFM 方法相比,EDFM 只需使用一组固定的结构化矩形网格,裂缝在前处理中与网格分开并显式精确刻画,这是EDFM 相对于DFM的关键优势[18].流固耦合研究方面,由于裂缝开度与油藏尺度存在极大差异,需对裂缝进行降维处理避免裂缝周围局部加密以减小网格量,而降维后裂缝两侧的位移场存在间断性[19].由于EDFM 采用非匹配性网格,使得其在表征裂缝力学间断效应方面极具挑战性.目前主要处理方法有: 有限元法、位移不连续法和扩展有限元法[20](XFEM).其中,XFEM 与EDFM 具有相似的优点,两者均基于结构化网格,分别相较于传统的基于有限元法的应力场、渗流场计算,避免了网格重构和非结构化网格剖分带来的困难.基于此,Ren 等[21-22]、Yan 等[23]、Zeng 等[24]先后实现了将EDFM 与用于应力场计算的扩展有限元法的耦合研究;Ren 等[25]进一步开发了pEDFM-XFEM流固耦合数值模拟器,能更精确地模拟多相流情况下的岩石和流体的力学性质.

在上述前人的研究中,仍亟待解决的主要问题有两点: 致密基质的渗透率极低,大规模体积压裂之后储层中的大尺度水力压裂缝与大量被激活的、小尺度的天然裂缝共存,对此类“人工裂缝性油藏”,单一数学模型无法进行准确表征,需要更一步发展综合数学模型;目前现有数值模拟方法在计算包含裂缝单元的基质网格内的压力分布时,采用了线性分布假设,这导致了早期及非稳态流动计算精度的不足.本文以上述突出的科学问题为导向,建立了致密油藏基质-天然裂缝-水力压裂缝综合数学模型,基于单孔模型表征基质渗流参数随着压力、应变场的改变,基于双重介质模型表征天然裂缝渗流参数的变化,基于嵌入式离散裂缝模型表征水力压裂缝的动态特征;并提出了XFEM-MBEM 的混合数值离散化方法,用于计算基岩和裂缝间的非稳态窜流,提高了模拟的早期精度,可准确表征致密油藏开采过程中的裂缝变形及流体流动机理,此研究为致密油气开发领域提供了理论基础.

1 流固耦合数学模型

1.1 物理模型及基本假设

本文采用了黑油模型假设,考虑基质和流体的压缩性,油水两相流动均符合达西定律,等温渗流,岩石变形考虑为微小线弹性变形,力学方程中采用常规弹性力学符号约定,即拉伸为正,压缩为负.

1.2 准静态力学平衡方程

图1 为弹性静力平衡状态下的多孔介质,中间存在一条水力压裂缝,pF为作用裂缝内表面的缝内孔隙压力,pp为支撑剂颗粒的支撑力,表示裂缝表面两侧 ΓF上的两个单位法向量,nt表示基岩应力外边界 Γt上的单位外法向量,t为应力外边界上的牵引力,为位移外边界 Γu上的平均固体位移.由此可建立准静态应力平衡方程式(1)～ 式(4)

图1 弹性静力平衡状态下的多孔介质Fig.1 Porous medium in elastostatic equilibrium state

式中,σ 为总应力张量,Pa;u为岩石固体位移,m;Pm为基岩孔隙压力,Pa;ρb为多孔介质基岩密度,kg/m3;g为重力加速度,N/kg.

1.3 油水渗流控制方程

基质体积单元中的油水两相质量平衡方程

离散裂缝单元中的油水两相质量平衡方程

式中,k为渗透率,m2;kro为油相相对渗透率,小数;krw为水相相对渗透率,小数; μ 为流体黏度,mPa·s;s为流体饱和度,小数;B为流体体积系数,小数;p为流体压力,Pa;要补充的是,上述式中下标o 和w 分别表示油相和水相,上标m和F分别表示基岩和水力压裂缝,上标w表示井筒表示基岩与裂缝之间的窜流量表示基岩与井筒之间的流量;φ*表示为Lagrangen 孔隙度,后续进行详细阐述.

1.4 基质-裂缝传质方程

EDFM 中存在4 类非相邻连接关系 (NNRs): 相邻基质网格之间的连接;位于不同基质网格中的相邻裂缝单元之间的连接;位于同一基质网格中的相邻裂缝单元之间的连接;基质网格与所包含裂缝单元之间的连接.其中,传导率可表示为流度与其几何因子的乘积,以油相为例,流度和几何因子可表示为

式中,λij表示单元i和j之间的油相流度;kro(upstream)表示利用上游权格式取值的油相相对渗透率;Gij为单元i和j之间的几何因子;TIij表示单元i和j之间的传导系数.

基质与裂缝之间传质项的处理是基于EDFM 中的第4 类连接关系,原始EDFM 方法[26]在计算包含裂缝单元的基质网格内的压力分布时,采用了Lee 等[27]及Praditia 等[28]提出的压力与到裂缝面的垂向距离成正比的线性分布假设如下

式中ki表示单元i的渗透率,m2;A代表接触面积,m2;〈d〉代表基质网格中心点与裂缝面之间的距离,m;S为所在基质网格的面积,m2.

但在实际物理情景中由于裂缝与基质的渗透率级别相差极大,这种求取“平均法向距离”的近似方法在处理非稳态传质时会导致一定的误差.因此,本文基于Fang 等[29]与Cao 等[30]所提出的混合边界元法,利用稳态渗流控制方程的边界积分方程推导得到双孔基质网格向裂缝网格传质量的新近似格式,以油相为例见式(16),EDFM 中MBEM 流量处理的连接关系示例见图2 所示.

图2 双孔系统EDFM 中MBEM 流量项处理的连接关系示例Fig.2 Example of connection list for MBEM handling flux in EDFM with dual-porosity system

式中,qomF为从裂缝单元流入基质网格的油相流量,m3;下标nF代表裂缝单元个数;h为基质网格厚度,等矩阵形式见附录.

1.5 井方程

经典Peaceman 公式[31]可用于计算井筒与基质或裂缝之间的流量.以油相为例,对于压裂水平井井筒与裂缝单元相连接

式中,wF为水力压裂缝的开度,m;re为等效补给井径,m;rw为油井井径,m;l1和l2分别为裂缝矩形单元的两个边长,m.

2 多孔介质本构关系

2.1 修正的Biot 方程

为准确描述非常规储层中基质与小尺度裂缝间的流动和岩石形变,可采用由Bai[32]提出的双孔有效应力原理,对基质和裂缝分别建立具有严格物理意义的有效应力关系

式中,σm表示基质总应力张量,Pa; σf为裂缝总应力张量为基质有效应力表示裂缝有效应力,Pa; αm表示毕渥系数,小数; ε 为总应变张量,m; εm为基质应变张量,m; εf为裂缝应变张量,m;上述式中下标f 表示小尺度裂缝.

Ren 根据式(19)～ 式(22)及胡克定律,指出双重介质有效应力的表达式为[22]

式中,Dmf为弹性刚度张量,Pa;Cm为基质柔度张量,Pa;Cf为裂缝柔度张量,Pa.

2.2 孔渗参数表征方法

孔隙度是关于体积应变和流体压力的函数[25].考虑线弹性微小形变下的基岩Lagrangen 孔隙度和大尺度裂缝Lagrangen 孔隙度可分别由下式表征

本文中将小尺度天然裂缝和基岩视为双重连续介质,因此小尺度裂缝表征方法与基质相同.对于大尺度水力压裂缝,生产过程中由于缝内流体压力损失会导致裂缝面的闭合,然而当存在支撑剂时,支撑剂发生形变会对裂缝面施加支撑力阻止其闭合.根据胡克定律,由支撑剂引起的支撑力可表示为

式中,Ep为支撑剂弹性模量,Pa;w为裂缝开度的位移,为水力压裂缝的原始开度,m.

此过程中,水力压裂缝的实际导流能力依赖于裂缝渗透率和开度的更新迭代[33]

3 耦合模型的数值离散

3.1 有限体积法离散流动方程

XFEM 与EDFM 具有相似的优点,两者均基于结构化网格,网格划分容易且物理意义清晰.因此本文采用正交矩形网格对基质进行几何离散,针对EDFM 前三类非相邻连接关系,采用满足局部物质守恒且简单易行的块中心有限体积方法来获取渗流控制方程的离散格式,仍以油相为例,此时为公式简洁,方程中去除可能存在的点源/汇项,对式(5)的两侧对时间和控制体积进行积分

式中左侧利用散度定理将体积分转换为面积分,并用符合物理意义的相邻网格之间的法向流量之和近似;右侧对时间的积分精确求解,对体积分采用矩形法估计,由此式(29)可以写为

对于相邻基质网格如图3,式(30)取隐式格式可改写为

图3 矩形网格离散示意图Fig.3 Sketch of discretization for rectangle grids

式中,n表示与该基质网格相邻的网格数量,ΔVi为该基质网格的体积,Δt为相邻时间步.

根据式(29)～ 式(31)的推导过程可以获得油水相的全隐式离散方程,但对于裂缝与基质连接的第4 类非相邻连接关系即含传质项,需要进行单独考虑.这里将已推导的多相流情况下EDFM 第4 类连接中基质网格与裂缝网格之间传质量的近似格式(16),耦合到基质与裂缝之间多相流方程的有限体积离散格式中,这也是混合边界元EDFM 的核心思想.仍以油相为例,对于基质系统

同理,对于裂缝系统

3.2 扩展有限元离散力学方程

XFEM 是在常规的有限元位移模式中根据单位分解的思想加引进一个跳跃函数和裂尖渐进位移场以反映位移不连续性的数值方法.XFEM 中的不同节点和单元的类型见图4,单元位移的近似解可以写成

图4 混合模型离散化方法示意图Fig.4 Sketch of hybrid model discretization strategy

式中,I表示网格中所有节点的集合,L表示包含所有裂缝段所在单元的节点集合,M为裂缝与裂缝相交的节点集合,K表示裂缝尖端所在单元的节点集合;Ni为节点位移形函数,ui为节点位移矢量,ai,j,ci,j,为单元增强节点的额外自由度,H表示Heaviside 函数,J表示裂间连接函数,F表示裂尖渐进函数,三者的表达式如下所示

将式(19)和式(20)代入式(1),考虑边界条件式(2)～ 式(4),得到方程(1)的等效积分弱形式见式(38),弱形式方程包含不同尺度裂缝以及基岩流体压力的相互作用

式中,δε 表示试应变,δu表示试位移,δu=δu+-δu-.

裂缝开度位移可以表示为

采用标准伽辽金有限元法对系统整体进行离散,可获得一个线性平衡方程组如下

式中,K表示整体刚度矩阵,f表示整体力向量.

3.3 整体计算格式及求解策略

本文分别采用有限体积法和扩展有限元法求解渗流控制方程及力学平衡方程,主要变量包括基岩及裂缝中的流体压力、含水饱和度以及分别设置在网格中心和角点的节点位移,对时间项的离散采用一阶向后欧拉全隐式差分格式,对完全耦合的全局方程组采用Newton-Raphson 迭代求解,并利用自动微分算法计算迭代过程中的雅可比矩阵

式中系数矩阵是非对称的,其中反对角项是主变量之间的耦合项.

4 模型验证分析

4.1 连续性介质流固耦合算例

此算例选取了经典的一维弹性土固结的孔弹性问题,将本文模型与Terzaghi[5]所提出的解析解进行对比.假设等温多孔介质由单相流体和固体组成,表现为线性孔弹性特征,模型的宽度为5 m,深度为10 m,顶部为排水边界,其余边界均不排水,底部为零位移边界,左右边界位移在水平方向上约束,于顶部施加20 MPa 的恒定载荷.岩土的孔隙度为0.4,渗透率为100 mD,弹性模量为20 GPa,泊松比为0.2,孔隙比为0.4,一维模型数值解与解析解对比结果如图5 所示,从图中可看出不同时间下本文提出的流固耦合数值解与Terzaghi 解析解的结果较为吻合,误差在合理范围之内,能证明该模型的准确性.

图5 一维模型数值解与解析解的对比Fig.5 Comparison between numerical solution and analytical solution of one dimensional model

4.2 离散裂缝单相流算例

此算例是通过多物理场耦合模拟的商业模拟器COMSOL 的标准有限元 (SFEM) 验证本文建立的流动模块.物理模型为矩形油藏,油藏中心有一条水力压裂缝,设置裂缝半长为基础长度将整个区域无因次化,裂缝无因次长度为4,油藏在X方向的长度为10 且Y方向长度为10,见图6.基于拉氏空间下无因次渗流方程,将本文修正后的EDFM 和有限元法与Blasingame 精确解[34]进行对比,其中裂缝无因次导流能力为200,有限元网格裂缝采用局部加密如图7 所示,裂缝加密的个数为20,如图8 所示.本文模型采用了非匹配性的矩形网格和混合边界元法,不对裂缝进行局部加密,分析离散裂缝压力动态曲线中不同数值算法的适应性.

图6 单裂缝矩形油藏模型Fig.6 Single fracture in rectangular reservoir model

图7 对裂缝网格加密示意图Fig.7 Sketch of fracture mesh refinement

图8 不同数值解与解析解的对比Fig.8 Comparison of different numerical solutions and analytical solution

模拟结果如图7 所示,模拟无因次时间从10-9～10,覆盖了有限导流能力的几个阶段,可以看出修正后的EDFM 方法与精确解无论从早期还是晚期都较为吻合,而有限元模型在早期与解析解的差距很大,中后期曲线重合,这在图9 的误差分析中也可以看出.原因在于有限元法处理网格累积项时,采用线性插值积分,在一个时间步内是稳态过程,而本文模型采用混合边界元法处理时采用了非稳态渗流的基本解,能精确描述初期的流动特征,根据对比结果可以说明本文模型的流动模块具有较高的早期精度.

图9 误差分析Fig.9 Error analysis

4.3 裂缝性介质应力强度因子算例

此算例是通过经典的Tada 解析解[35]验证本文建立的力学模块计算裂缝尖端应力强度因子 (SIFs)的准确性.考虑一个平面应变模型 (10 m×10 m),上边界为施加均匀的拉应力,下边界为固定位移的滚轮边界,模型内部中心存在一条倾斜裂缝,裂缝半长l为1 m,如图10 所示.解析解中应力强度因子是外部应力和裂缝几何尺寸和倾角的函数,因此为方便对比,定义无因次应力强度因子:在裂缝尺寸保持不变的条件下,计算不同裂缝倾角无因次应力强度因子数值解与解析解的结果见图11,对比显示本文模型的数值解与经典解析解的结果几乎保持一致,具有较高的可信度.

图10 裂缝性介质示意图Fig.10 Schematic of fractured medium

图11 不同裂缝倾角无因次应力强度因子数值解与解析解的对比Fig.11 Comparison of numerical and analytical solutions for dimensionless SIFs with different fracture dip angles

5 典型实例分析

如图12 所示的致密油藏,考虑分为压裂改造区(SRV) 与非压裂改造区 (USRV),SRV 内部存在一口水平井和7 条水力压裂缝,水力压裂缝中考虑支撑剂的作用.采用平面应变模型,油藏边界条件如图12所示.SRV 内考虑微裂缝的影响,采用双孔模型;USRV 不考虑微裂缝的影响,采用单孔模型.致密油藏的基本参数见表1,设置该水平井为定压生产1000 d,井底流压为10 MPa,并将模拟的结果与商业模拟器COMSOL 的标准有限元 (SFEM) 作对比,图13 为两种数值模拟器的网格剖分,其中COMSOL-SFEM 中对SRV 内部及裂缝尖端进行了局部加密处理,本文模型XFEM-MBEM 中采用非匹配性的矩形网格.图14 为两种数值模拟器在开发1000 d 时刻的压力场对比,图15 和图16 分别为X-位移场与Y-位移场的对比,从上述场图可看出两种数值模拟器的结果基本一致.以图12 中SRV 内部最中间的那条水力压裂缝为研究目标,分析其不同时刻的开度分布见图17,图中能看出XFEM-MBEM 与COMSOL-SFEM不同时刻裂缝开度的计算结果较为吻合,可以说明XFEM-MBEM 流固耦合模拟的准确性.此外,从图中可以随着油藏的不断开发,裂缝呈不断闭合的趋势,定压生产下靠近井筒附近的裂缝段闭合更加剧烈.图18 为该模型在表1 参数条件下考虑流固耦合作用与不考虑流固耦合作用的生产动态曲线对比,可以看出应力场所引起的渗流参数改变及裂缝开度降低对产能的影响很大,因此对致密油藏压裂水平井进行产能评价时,地质力学效应的影响不可忽略.

图12 致密油藏示意图Fig.12 Sketch of tight oil reservoir

表1 致密油藏基本参数Table 1 Basic parameters of tight oil reservoir

图15 两种数值模拟器开发1000 d 的X-位移场对比(单位: mm)Fig.15 Comparison of X-displacement distribution maps over 1000 days of two simulators (unit: mm)

图16 两种数值模拟器开发1000 d 的Y-位移场对比(单位: mm)Fig.16 Comparison of Y-displacement distribution maps over 1000 days of two simulators (unit: mm)

图17 不同时刻裂缝开度分布Fig.17 Fracture aperture distribution at different time

图18 考虑与不考虑流固耦合作用的生产动态曲线对比Fig.18 Comparison of production dynamic curves with and without flow-geomechanics coupled effect

6 结论

(1)本文提出了基于XFEM-MBEM 的嵌入式离散裂缝模型流固耦合数值模拟方法,采用正交的非匹配型网格进行几何离散,其优势在于裂缝独立于计算网格,与常规离散裂缝模型基于裂缝形状的匹配型网格相比,网格划分难度大大降低,能避免当裂缝间距或夹角非常小时,网格划分质量不佳等问题.本文通过几个实例分别验证了其力学模块、流动模块、耦合模块的准确性.

(2)传统EDFM 方法对于第4 类非相邻连接关系在计算包含裂缝单元的基质网格内的压力分布时采用了线性分布假设,导致了开发早期及非稳态流动计算精度的不足.局部网格加密技术可以减小包含裂缝网格的基质网格尺寸,从而有效解决这个问题,但会对前处理算法带来巨大的复杂性.本文采用MBEM 方法获取了双孔基质网格向裂缝网格传质量的新近似格式,能在不加密情况下提高早期计算精度.

(3) SRV 改造区内的微裂缝对提高储层流体的动用程度有重要影响.然而,用EDFM 和XFEM 进行显式表征是较为困难的,且会增加计算耗时.本文采用双孔有效应力原理耦合双重介质模型隐式表征SRV 区内部的天然裂缝,能有效模拟具有复合分区特征的此类裂缝性油藏问题.

(4)致密油藏开发过程中随着缝内流体被采出,裂缝会出现不同程度的闭合,从而导致了整体导流能力降低,油藏采收率将降低.考虑到非常规储层原油的产量对裂缝的依赖性如此之高.因此,考虑开发过程中流固耦合作用引起的参数变化对产量的评价具有重要意义.

附录