基于连接元法的多尺度裂缝性油藏数值模拟1)

湛文涛 赵 辉 饶 翔 刘 伟 徐云峰

(长江大学石油工程学院,武汉 430100)

1 引言

随着人工压裂技术相关材料、方法和工艺取得新的进展,水平井压裂技术成为超低渗油藏、致密油藏等非常规油气资源商业化开发中的主要技术手段[1].在压裂过程中,地层中会形成许多导流裂缝,相比于天然裂缝,导流缝通常尺度较大,渗透率很高,但对于油藏区域而言,裂缝的尺度相对小,水力压裂缝网储层具有多尺度的特征.裂缝几何(形状,长度以及方向)对流体流动具有显著的影响[2-3],分析这些影响因素,裂缝油藏数值模拟模型是必要的.

目前裂缝性油藏数值模拟方法主要分为两类,连续介质模型和离散介质模型[4],其中,连续介质模型主要包括等效连续介质模型[5]和双重介质模型[6-7].在连续介质中,裂缝和基质被处理成等维度的离散单元.在离散介质模型中,n维油藏系统中的裂缝被作为n-1 维处理,即在二维的油藏系统中,裂缝被认为是一维的线段,而三维问题中裂缝则是二维的平面.这是这两种方法在处理裂缝形态上本质的区别.在连续介质模型和离散介质模型下,一些学者采用有限差分[8-11]、有限元[12-14]和有限体积[15-17]等数值模拟方法进行流动模拟计算,这些数值模拟计算方法都是基于网格剖分离散表征计算域.在连续介质模型中,裂缝和复杂边界几何的精细描述需要更小的网格,这会造成数百万计甚至数亿计的模型自由度,极大增加计算耗时.在离散介质模型中,虽然将裂缝降维处理,能够降低裂缝表征的难度,但高质量的匹配性网格难以生成.嵌入式离散裂缝模型基于正交网格的基质离散方案无法适用于实际油藏复杂边界,针对该问题程林松等[18]提出基于两套节点的格林元法耦合嵌入式裂缝模型能适用于任意形态的基质网格,可实现复杂油藏边界问题的求解,但该计算过程相对复杂.

近年来,无网格法在计算力学以及流体动力学领域得到了广泛应用[19-22],常用的无网格法有加权最小二乘[23]、广义有限差分[24]、无单元Galerkin 法[25]以及光滑粒子流体运动学法[26].无网格法通过点云的方式离散计算域,能够更加灵活地刻画裂缝和复杂边界,对于油藏复杂几何的描述相比于网格类方法更加简单.一些学者将无网格运用于裂缝性油藏,比如将无网格法应用到裂缝扩展建模中,结合应力分析裂缝几何形状演变过程[27].基于无网格广义有限差分法和嵌入式离散裂缝模型,建立一种裂缝性油藏数值模拟无网格方法[28],验证其方法对于复杂几何形态刻画的优势性,但该方法还存在着计算效率问题和不能直观反映井节点间的连通关系.基于井间连通性思想[29-33],赵辉等[34]引入广义有限差分理论,建立非欧物理连通网络等效模型,推导满足物质守恒且具有物理意义的渗流特征参数,提出一种新的油藏数值模拟计算方法—无网格连接元法(connection element method,CEM),为油藏数值模拟提供了新思路.基于连接元核心思想,定义无网格节点控制体积,Rao 等[35]发展了扩展有限体积法,本质上它是连接元法的全隐式格式.

目前对于复杂裂缝网络的刻画,传统方法存在着匹配性网格自适应生成困难和结构化网格难以适用于裂缝网络和油藏边界的复杂几何等问题.在应用于实际油藏时,数值模型的网格数将十分巨大,导致计算效率低、历史拟合难,同时难以获取注采井间的动态连通性以及流线追踪.本文将连接元法应用于裂缝性油藏的开发动态模拟,不同于网格体系(匹配网格或者结构化网格),针对多尺度裂缝储层的表征,连接元法沿着裂缝的走势配置相应的不同尺度的连接单元,精细刻画裂缝.该方法具有无网格特征,离散表征更加自由和灵活,能够更加容易刻画裂缝和复杂油藏边界.此外,物理连通网络的连接关系相比于网格体系和无网格点云能够更加直观揭示井间的连通关系.

2 裂缝油藏连接元法基本原理

区别于网格体系,该方法离散油藏计算域时不需要网格剖分,能够避免网格体系表征油藏复杂几何(裂缝、断层以及不规则边界)的困难.从流动的角度,将油藏计算域离散成为一系列连接单元构成的连接网状拓扑结构,采用全隐式的离散方案基于牛顿迭代法可实现渗流控制方程的高效求解.以油水两相流为例,介绍裂缝性油藏连接元方法的基本原理.

2.1 裂缝油藏连接元单元体系构建

在油藏数值模拟方法中,基于网格体系的离散方式,其主要目的是获取离散后的网格单元体间的连接拓扑结构,按照其拓扑连接关系分析势场和流场的相互作用关系.换而言之,网格剖分的本质是建立离散单元体之间的连接关系,但基于网格体系的离散会限制其连接关系,它们只能建立相邻网格间的连接关系,比如正交网格、三角剖分网格、PEBI网格离散.连接单元体系是一种物理网络结构,它是由节点间的连接单元所构成.如图1,利用连接单元将某个油藏离散为连接单元体系.显然,基于连接单元体系离散避免了网格剖分对于不规则油藏边界和裂缝的匹配性网格离散的困难,尤其对于尺度差异大的裂缝表征.这种方案几何离散更加简单,不再受限于网格相邻关系,而是更加灵活直观地表征节点间的相互关系.

图1 油藏的连接单元离散Fig.1 Reservoir discretization based on connection element

在实际油藏开发过程中,往往需要在地层中布置合理的生产井和注入井.在油藏数值模拟中,部署的井往往是作为源汇项处理.基于网格体系的传统数值模拟对于源汇项问题,往往是在背景网格的控制域上采用积分的方式处理.基于连接单元体系的节点没有实质控制域,为方便处理带有源汇项的物理问题,在节点处需要给出一个控制体积域的概念.下面基于物质守恒原理给出节点控制域的定义,所有节点间的控制域不相交,且所有控制域之和为整个油藏区域

其中,Ω 是整个油藏控制域;Ωi是节点i的控制域;VΩ是油藏总体积;Vi是节点控制体积.

2.2 渗流特征参数计算

连接单元体系的本质将实际三维流场转换为由一维连接单元构成的物理连通网络,为在连接单元体系上对渗流方程进行计算,需要建立表征连接单元的渗流特征参数.定义节点间连接传导率来表征连接单元的流体渗流能力,连接传导率与势场(如渗流问题的压力场)中对势的计算(压力)相关.下面以两相流控制方程为例,推导裂缝性油藏连接单元渗流特征参数的计算方法.

2.2.1 基质间渗流参数表征

首先给出双重介质基质系统两相流渗流控制方程

式中,km和krσ分别是基质系统绝对渗透率和相对渗透率,mD;μo和 μw是油和水黏度,mPa·s;∇ 是哈密尔顿梯度算子;pσ,m是基质层压力,MPa;sσ,m是基质层饱和度,1;t是时间,d;qosi是地面标况下的体积流量(源汇项),d-1;δ是狄拉克函数,1;τmf是基质与裂缝之间物质交换的量.

对于式(2),压力扩散项在节点控制域 Ωi内积分,得到

式中,i为影响域中心节点,该影响域包含着其他ni个节点,代表以i节点为中心节点计算的连接单元(i,η) 的传导率.考虑以节点j为中心的影响域(包含节点i),该影响域包含着其他nj个节点,可以得到

根据连接单元物质平衡原理及所有节点的控制体积之和等于油藏体积的原则,可以求得每个节点的控制体积[28,34-35],进而可以定义基质节点间的连接单元 (i,j)的传导率为

2.2.2 显式裂缝间渗流参数表征

首先给出双重介质裂缝系统两相流渗流控制方程

式中,kf和krσ分别是裂缝绝对渗透率和相对渗透率,mD;μo和 μw是油和水黏度,mPa·s;∇ 是哈密尔顿梯度算子;pf,σ是基质层压力,MPa;sf,σ是基质层饱和度,1;t是时间,d;qσ,well是地面标况下的体积流量(源汇项),d-1;δ是狄拉克函数,1;τmf是基质与裂缝之间物质交换的量.

对于裂缝系统而言,类似于网格体系的传导系数定义裂缝层两裂缝节点间的传导率为

2.2.3 两相流基质与显式裂缝间渗流参数表征

对于基质与裂缝间物质交换的刻画,Moinfar等[36]将嵌入式离散裂缝模型(EDFM)运用到三维问题,定义基质网格向裂缝网格的传导系数为

式中,A是裂缝与基质的界面面积,m2;k是渗透率张量,n为 法向向量,〈d〉为f与m之间的平均法向距离,m.

然而本文方法对于节点控制体积没有一个具体的形状描述,针对裂缝与基质的界面面积A的刻画,采取节点控制域内裂缝总长度Lf,i与裂缝高度hf,i的乘积.对于某中心节点i,Lf,i被取作所有以节点i为端点的裂缝连接单元的一半的总和

式中,Λ表示以节点i为端点的裂缝连接单元的其他端节点构成的集合;Lf,iζ是以节点i和节点 ζ为端点的裂缝连接单元 (i,ζ)的长度,m.

在处理基质节点与裂缝节点间的物质交换量的表达式为

本文考虑裂缝贯穿整个油藏储层厚度,即裂缝高度与油藏储层厚度h相同.是节点i的控制域内基质与裂缝窜流的调和平均渗透率,mD.d是裂缝与基质窜流的等效法向距离,m

因此定义基质层与裂缝层窜流的传导率为

2.3 渗流控制方程离散求解

基于局部双重介质的裂缝油藏两相流渗流控制方程如下.

基质系统

裂缝系统

结合离散节点控制域 Ωi可将方程离散,其渗流控制方程的离散形式如下

对于油相和水相的相对渗透率,采用油藏数值模拟中常用的迎风格式,即

边界条件主要是第一类边界条件(Dirichlet 边界)和第二类边界条件(Neumann 边界).处理第二类边界条件时,需要在边界处加入虚拟点,辅助计算边界处的导数[35].虚拟点的加入方式在参考文献里有详细介绍,这里不再赘述.

2.4 连接单元体系路径追踪

基于网格体系的方法难以直观获取各井之间的相互作用关系.目前,针对裂缝油藏尚未形成简单实用的连通性定量表征方法,而基于连接体系的连接元法可通过INSIM-FPT[31]中的路径搜索算法获取各节点之间的所有流动路径及各流动路径的劈分系数.劈分系数的数学描述如下,假设在第n个时间步,上游节点i与下游节点j之间直接相连(即存在连接单元),则节点i到节点j的劈分系数为[31]

式中,是与节点i相连的下游节点数.

劈分系数反映了上游节点处的流体流到下游节点的比例,体现了节点之间的流动相互作用,从而可以直观揭示注水受效、水窜识别等矿场实际问题.此外,在某一时间步计算得到各节点控制体积的平均压力后,在无网格连接体系的基础上会形成压力高低判别的有向图,即对于某连接单元i-j,如果在第n时间步,节点i的压力高于节点j,则连接单元i-j的方向定义为由i指向j.在这样一个有向图的基础上,可以采用图论中的路径搜索算法获取各节点之间的所有路径.

3 应用

本节主要的目的是探索裂缝性油藏连接元的计算性能.下面给出几个数值算例来验证本文方法的有效性和优越性,包括平行多缝单相流,复杂边界平行多缝两相流.引入L2范数误差函数,以精细网格剖分的EDFM 作为参考解,对比压力、含水饱和度场图以及含水率曲线.此外引入路径追踪方法,可以在连接单元体系下获取流动路径和分析井节点间的连通性,以此体现本文方法独特的优越性.下面给出两个算例相同的物性参数见表1

表1 油藏物性参数Table 1 Physical parameters of reservior

其中,u(i) 是计算值,uref(i)是参考值,T是节点(网格)总数.

3.1 倾斜平行缝网单相流

油藏尺寸为960 m×520 m×10 m,油藏中心一口水平井P1,在地层中射开6 条平行倾斜裂缝,裂缝纵向上贯穿油藏,裂缝厚度与油藏厚度相同.基质渗透率0.1 mD,以每天10 方的产量进行衰竭开发,模拟计算生产500 天.油藏网格剖分为120×65×1,网格大小Dx=Dy=8 m,Dz=10 m,网格总数7800,采用嵌入式离散裂缝模型精细解作为参考解.如图2 (a)所示,给出了24×13 粗化网格的嵌入式离散裂缝油藏模型,网格大小Dx=Dy=40 m,Dz=10 m,网格总数312.如图2 (b)所示,以等间距Dx=Dy=40 m,Dz=10 m 的24×13×1 的构建连接元模型,影响域半径为,总节点数312,连接单元总数2173.

图2 裂缝性油藏模型离散Fig.2 The discretization of fractured reservoir model

根据渗流控制方程(19),对于单相流问题,通过EDFM 方法和CEM 方法求解压力,计算整个油藏的压力场分布.图3～图5 分别是EDFM 精细网格剖分、EDFM 粗化网格剖分和CEM 法在第100 天和第500 天的压力场图.结果表明,本文方法与参考解是一致的,说明了该方法的有效性和正确性.在24 ×13 配点模型下,统计计算机CPU 耗时(计算机型号:12 th Gen Intel(R) Core(TM) i5-12400 F),基于精细网格的EDFM 计算耗时138.14 s,粗化网格的嵌入式离散裂缝模型计算耗时29.76 s,压力场图的计算精度96.2%,而连接单元法的计算耗时31.24 s,压力场图计算精度99.1%.在计算耗时相当的情况下,计算精度提高了2.9%.

图3 EDFM 计算压力,8 mFig.3 The pressure calculation of EDFM,8 m

图4 EDFM 计算压力,40 mFig.4 The pressure calculation of EDFM,40 m

图5 CEM 计算压力,40 mFig.5 The pressure calculation of CEM,40 m

下面从影响域半径和配点间距两个方面分析本文方法的稳定性.首先采用均匀离散配点的方式离散油藏计算域,分别取如图6 所示的3 种不同节点影响半径构建连接单元体系,需要说明是,从物质流动的角度出发,共线的3 个点,只取相邻两点构建连接单元,计算不同模型的误差如图7 所示,结果表明,过小或者过大的影响域半径会降低计算精度,这个结果也与无网格法的影响域半径对计算精度的敏感性一致.因此,为获取相对高的计算精度,本文在算例验证中影响域半径取.

图6 不同影响域半径示意图Fig.6 Diagram of different influence radius

图7 不同影响域半径压力计算误差Fig.7 Calculation error of pressure with different radius of influence

图2 中已经以等间距Dx=Dy=40 m,影响域半径re=构建连接元模型,下面分别给出以等间距Dx=Dy=20,10 m,两种布点方案,影响域半径+0.01构建连接元模型,计算100 d 和500 d 的压力分布如图8 所示.图9 是3 种不同配点间距连接元模型的压力计算误差,结果表明,随着配点间距越小,计算精度越高,说明本文方法具有良好的收敛性.

图8 不同间距布点CEM 计算压力Fig.8 The pressure is calculated by CEM at different collocation intervals

图9 不同间距布点CEM 压力计算误差Fig.9 Calculation error of CEM pressure at different collocation intervals

3.2 不规则边界平行缝网油水两相流

设计一个不规则边界油藏的概念算例,在油藏中布置有3 口水平生产井以及7 口注水井,裂缝纵向上贯穿油藏,裂缝厚度与油藏厚度相同,油藏区域如图10 (a)所示.油藏总体积为4.212 ×105m3,基质渗透率10 mD.3 口水平井均以350 m3/d 产液量,7 口注水井以150 m3/d 注入量的生产制度进行开发.如图10 (b),以网格大小为Dx=Dy=4 m,Dz=10 m 离散油藏计算域,网格总数33325,以此精细网格剖分的EDFM 数值模型作为参考解.图10(c)是连接元油藏离散模型,节点总数1053,连接单元总数7868.连接单元法沿着油藏边界配置相应的节点,相比于网格离散,能更灵活地匹配实际油藏边界.为了对比边界处的渗流场的计算精度,在油藏左下边界上取一个观测点A(如图10 (b)),对比注采过程中压力和含水饱和度的计算结果.

图10 不规则油藏模型Fig.10 Irregular reservoir model

以精细网格剖分EDFM 计算压力、含水饱和度、含水率以及产油速度等参数作为参考解.结合两相流的控制方程离散格式,采用连接单元法计算油藏的压力和饱和度分布,求解井点的含水率以及产油速度等参数.如图11,分别给出了第500 d EDFM 与CEM 计算的压力结果.图12 是第100 d和500 d EDFM 与CEM 计算的饱和度,本文方法计算的结果与参考解一致.图13 是两种方法计算观察点处的压力和含水饱和度动态曲线.结果表明,本文方法在边界处具有较高的计算精度.通过压力和饱和度的计算,说明了裂缝性油藏连接元法具有较高的计算精度.此外,对比了3 口水平生产井的含水率和产油速度等动态参数,如图14 所示,连接元的井点参数计算结果与参考解结果一致.基于精细网格的嵌入式离散裂缝模型计算耗时621.41 s,连接单元法的计算耗时126.53 s,3 口生产井P1,P2 和P3 的含水率的计算精度分别为90.41%,88.26% 和91.89%.连接单元法在较少的节点体系下取得了较高的计算精度,相比于精细的嵌入式离散裂缝模型而言,计算速度提高了近5 倍.因此,本文方法能够取得计算精度和计算效率的更优平衡.此外,引入路径追踪方法,计算7 口注水井的劈分系数,如图15 所示.

图11 压力场图Fig.11 Pressure profile

图12 含水饱和度场图Fig.12 Water saturation profile

图13 边界观测点结果对比Fig.13 Comparison of boundary observation points

图14 含水率与产油速度曲线Fig.14 Water cut and oil rate curves

图15 注采连通示意图Fig.15 Injection-production connectivity profile

4 结论

(1)本文构建一种基于无网格连接单元体系的裂缝性油藏数值模拟方法,该方法具有无网格特征,相比于网格离散,对于裂缝几何形态、分布以及不规则油藏边界的精细刻画非常自由灵活,能够更加适用于裂缝性油藏复杂几何特征的刻画.

(2)相较于传统方法,本文方法在离散渗流控制方程时,构造了未知函数导数的多点差分近似格式,具有更高的估计精度,即使在相对粗化的模型下,也能具有更加丰富的节点(网格)间连接拓扑结构.因此,该方法能够通过粗化模型,降低油藏数值模型计算自由度,提高计算效率的同时,保证了较高的计算精度,能够取得计算精度和计算效率的更优平衡.

(3)相较于传统方法,本文方法能够在粗化模型下,利用路径追踪算法高效直观地获取节点间的相互作用,计算无网格连接体系下各流动路径的注水劈分系数,实现裂缝性油藏井间连通关系的定量识别.