A NUMERICAL SIMULATION APPROACH FOR EMBEDDED DISCRETE FRACTURE MODEL COUPLED GREEN ELEMENT METHOD BASED ON TWO SETS OF NODES
-
摘要: 对于原始嵌入式离散裂缝模型(EDFM), 在计算包含裂缝单元的基质网格内的压力分布时采用了线性分布假设, 这导致了油藏开发早期对非稳态窜流量的计算精度不足. 因此, 本文提出了一种两套节点格林元法的EDFM数值模拟方法. 两套节点格林元法的核心思想是将压力节点与流量节点区分开, 一套压力节点设置在单元顶点, 另一套流量节点设置在网格边的中点, 满足局部物质守恒、具有二阶精度的同时, 可适用于任意网格类型. 本文将两套节点格林元法与EDFM耦合, 采用了非稳态渗流控制方程的边界积分形式推导了基质网格与裂缝网格之间传质量的新格式, 代替了线性分布假设以提高模拟精度; 此外, 修正后的EDFM能适应任意形态的基质网格剖分, 拓展了原始EDFM仅适用于矩形基质网格、难以考虑复杂油藏边界的局限性. 研究表明: 通过对比商业模拟软件tNavigator® LGR模块与原始EDFM, 验证了本文模型具有较高的早期计算精度; 以复杂油藏边界−缝网−SRV分区模型为例, 通过对比SFEM-COMSOL商业模拟软件, 验证了本文模型处理复杂问题的适应性. 本文研究可用于裂缝性油藏开发动态的精确模拟.Abstract: For the original embedded discrete fracture model (EDFM), the linear-distribution assumption is adopted in calculating the pressure distribution in matrix grids containing fracture elements, which leads to the lack of accuracy in solving unsteady interflux in the early stage of oil reservoir development. Therefore, this paper proposes a numerical simulation approach for EDFM coupled Green element method based on two sets of nodes. The main idea of the Green element method with two sets of nodes is to distinguish pressure nodes from flux nodes, in which one set of pressure nodes is set at the vertex of grids and another set of flux nodes is set at the edge-midpoint of grids. It not only meets the local material conservation and has second-order accuracy, but also can be applied to any grid type. In this paper, the Green element method based on two sets of nodes is coupled with EDFM, and a new scheme of mass transfer between matrix cell and fracture elements is derived by adopting the boundary integral form of the unsteady flow control equation, which replaces the linear distribution assumption to improve the simulation accuracy. In addition, the modified EDFM adapts to any form of matrix mesh generation, which extends the limitations of the original EDFM which is only suitable for rectangular matrix mesh and difficult to consider complex reservoir boundaries. The research shows that the proposed model has high accuracy in the early stage and it is verified by the LGR module of commercial software tNavigator® and the original EDFM. Taking the SRV-zoning model considering fracture networks and complex reservoir boundaries as an example, the flexibility of the proposed model for solving complicated problems is demonstrated by comparing the business simulation software named SFEM-COMSOL. This study can be used for the accurate simulation of dynamic production performance in fractured oil reservoirs.
-
引 言
对非常规裂缝性储层的压力和流体分布进行精确模拟, 是许多能源工程技术人员面临的一个热点和挑战性工作[1]. 目前为止, 已有许多数值方法被用于解决这一物理问题, 如等效连续介质模型[2]、离散裂缝模型[3]和嵌入式离散裂缝模型(embedded discrete fracture method, EDFM)[4-5]. 由于具有较高的关键地质特征可视化水平, 其中EDFM的应用最为广泛. 但由文献[6-8]所提出的原始EDFM在计算基质网格与裂缝单元间传导系数的几何因子时, 采用了压力与到裂缝面的垂向距离成正比的线性假设, 由于非常规储层裂缝与基质的渗透率极差较大, 这种理想假设无法准确描述早期的非稳态窜流, 会导致一定的误差; 此外, 原始EDFM多采用仅适用于矩形网格剖分的有限差分或有限体积法, 难以考虑复杂油藏边界, 技术应用层面存在局限性.
针对EDFM早期精度不足、网格适用性不佳的问题, 目前国内外有两个主要研究思路: 其一是优化EDFM网格剖分的前处理算法[9-11]; 其二是寻求一种更加高效、适用性更强的数值解法. 针对基质网格和裂缝单元之间传质量的准确计算问题, 由于边界元法相较于有限元法具有半解析精度, 可能会是一种较为有效的思路[12]. 然而, 由于实际储层的非均质性极强, 原始边界元法无法有效处理这一问题. 格林元法(Green element method, GEM)本质上是边界元法的一种变体, 在剖分的每个网格中建立边界积分方程, 结合了有限元法的变分原理, 更适用于解决非均匀介质的非线性问题. 原始GEM由Taigbenus等[13-17]提出, 核心思想是将计算域由多边形单元划分, 单元顶点被视为未知节点, 通过用网格内部压力近似表达式的法向导数估计边界上的法向流量, 但仅显式考虑了节点的压力值, 因此总体精度不高, 误差会随着多边形单元尺寸的减小而增大. 后续的学者对原始GEM进行了改进, Archer等[18-19]指出, 采用较大网格范围内构造的插值基函数能减少原始GEM的误差; Pecher等[20]和Lorinczi等[21-24]提出了基于流量向量的格林元法; Taigbenus[25-26]后续发展了修正流量项的GEM, 但始终未解决无法显式求解的根本问题; 方思冬等[27]结合了混合有限元法和边界元法二者的优势, 将压力节点和流量节点设置在网格边的中心点, 建立了物理意义明确的混合边界元法; 为提高GEM的鲁棒性, Rao等[28]提出了模拟格林元法, 该方法耦合了模拟有限差分法的核心思想, 可满足局部质量守恒, 能稳定地求解二阶偏微分方程; Rao等[29]还提出了两套节点格林元的基本思想, 但没有实现与EDFM的耦合及相关应用. 综合上述研究来看, 目前大部分格林元法在精度、收敛性及适应的网格类型上各有优劣性, 但若要将格林元的优势引入EDFM的发展中, 就必须建立一种同时保持高精度、强鲁棒性且能适用于非结构网格的格林元法.
本文从原始EDFM的适用局限性出发, 利用格林元方法在处理非定常和非均质问题方面具有高精度的独特优势, 建立一种能与EDFM合理匹配的改进格林元法, 能精确求解基质与裂缝间的非稳态窜流问题, 拓展原始EDFM技术应用层面的局限性, 使其既能求解各类复杂油藏渗流问题, 又能使求解结果具有明确的物理意义.
1. 格林元法雏形及其局限性
本文研究的是二维孔隙空间中的单相等温渗流问题, 流动遵守Darcy定律, 并考虑封闭外边界条件. 原始GEM[13]的基本思想是将边界积分方程应用至各个子单元中, 离散格式与直接法推导的有限元法较为相似, 被广泛应用在如式(1)的二阶偏微分方程中
$$ \nabla \cdot (K\nabla p) = c\frac{{\partial p}}{{\partial t}} + f $$ (1) 式中, p是计算域内基质的孔隙流体压力, MPa; K和c为渗流介质的属性参数函数, 均为标量; f是计算域内的源汇项强度, 若渗流介质某基质网格内含有裂缝单元, 则有
$f = \dfrac{{{q_{{\rm{omf}}}}}}{V}$ , 其中${q_{{\rm{omf}}}}$ 表示该基质网格与裂缝单元间的传质速度, m3/d;$ V $ 是基质网格体积, m3.上式可以改写为
$$ {\nabla ^2}p = - \nabla \psi \cdot \nabla p + \sigma \frac{{\partial p}}{{\partial t}} + \nu f $$ (2) 式中,
$\psi = {\rm{ln}}K$ ,$ \nu = 1/K $ ,$ \sigma = cv $ .拉氏方程在无限大计算域的基本解
$G\left( {{\boldsymbol{M}},{{\boldsymbol{M}}_i}} \right)$ 为$$ {\nabla ^2}G\left( {{\boldsymbol{M}},{{\boldsymbol{M}}_i}} \right) = \delta ({\boldsymbol{M}},{{\boldsymbol{M}}_i}) $$ (3) 式中,
${{\boldsymbol{M}}_i}$ 为空间中的源点坐标,${\boldsymbol{M}}$ 为空间中任意一点的坐标,$\delta ({\boldsymbol{M}},{{\boldsymbol{M}}_i})$ 为Delta函数.根据格林第二公式和式(3), 可以得到边界积分方程如下
$$ \begin{split} & - \lambda {p_i} + \int_{\partial \varOmega } (p\nabla G \cdot {\boldsymbol{n}} - G\nabla p \cdot {\boldsymbol{n}}){\text{d}}s + \\ &\qquad \iint_\varOmega {G\left( - \nabla \psi \cdot \nabla p + \sigma \frac{{\partial p}}{{\partial t}} + \nu f\right)} {\text{d}}A = 0 \end{split} $$ (4) 式中,
$ \lambda $ 为边界积分方程中所选取源点所处位置的特征角度,$ \varOmega $ 和$ \partial \varOmega $ 分别表示计算域及计算域的边界, n是计算域边界上的外法向量,$ i $ 是所选取源点的序号.Taigbenu等[13-17]利用了法向流量的概念, 即令
$ q = - K\nabla p \cdot {\boldsymbol{n}} $ , 则上式改写为$$ \begin{split} & - \lambda {p_i} + \int_{\partial \varOmega } \left(p\nabla G \cdot {\boldsymbol{n}} + G\frac{q}{K}\right){\text{d}}s +\\ &\qquad \iint_\varOmega {G\left( - \nabla \psi \cdot \nabla p + \sigma \frac{{\partial p}}{{\partial t}} + \nu f\right)} {\text{d}}A = 0 \end{split} $$ (5) 从式(5)可以看出, 式中存在与流体压力和介质属性参数有关的域积分, 由此体现了GEM与边界元法的显著差异性, 其可以有效处理计算域非均质及参数非线性等更为复杂的问题. 以前学者提出过的原始GEM[13-17]、基于流量向量的GEM[21-24]及流量修正的GEM[25], 都是在上述GEM雏形基础上的进一步发展, 此类经典方法存在的突出共性问题在于子单元顶点处的外法向流量
$({q_n}=\nabla p \cdot {\boldsymbol{n}}) $ 不具备整体连续性, 使得局部物质不守恒.2. 压力流量两套节点格林元法
为此, 本文提出了一种压力流量两套节点的改进格林元法, 用一组节点包含表示压力值的多边形单元的所有顶点作为压力节点; 用另一组包含表示法向流量值的网格边的中心点作为流量节点. 其中心思想是通过将流量节点设置在网格边的中点, 并作为常单元, 以此满足局部物质守恒. 该方法可适用于任意形状的网格剖分, 下述推导以三角形网格为例, 如图1所示, 改进的格林元法具有压力和流量两套不同节点, 其中压力节点分布在三角形单元的三个顶点, 流量节点分布在三角形单元三条边的中心点. 压力节点分别标记为1, 2, 3, 流量节点分别标记为a, b, c.
![]()
图 1 压力和流量两套节点三角形单元示意图Figure 1. Sketch of pressure and flux two sets of nodes in triangle cell对三角形单元获取如式(4)的边界积分方程, 其包含有与压力有关的积分项, 需对压力节点值进行插值来估计压力函数在整个单元上的取值, 通常压力节点仍然选择在多边形网格的顶点, 以保证插值函数之和在整个单元内的取值均为1. 每个压力节点相应的基函数表示为
$ {\phi _i} $ , 三角形单元内的压力值可由节点值和基函数的加权平均得到, 即$$ p = \sum\limits_{i = 1}^3 {{p_i}} {\phi _i} $$ (6) 但对于每个流量节点, 可根据法向流量的分段连续性, 将网格边上的法向通量视为一个常量, 这是与原始格林元法截然不同的地方. 流量节点所在的边记作Γz (z = a, b, c), 按上述思想, 式(4)的离散格式为
$$ \begin{split} & \sum\limits_{j = 1}^3 {{R_{ij}}{p_j}} + \sum\limits_{z = a}^c {{L_{iz}}{{\left(\frac{q}{K}\right)}_z}} - \sum\limits_{j = 1}^3 {\sum\limits_{m = 1}^3 {{U_{imj}}{\psi _m}{p_j}} } + \\ &\qquad \sum\limits_{j = 1}^3 {\sum\limits_{m = 1}^3 {{W_{imj}}\left[{\sigma _m}\frac{{{\text{d}}{p_j}}}{{{\text{d}}t}} + {\upsilon _m}{f_j}\right]} } = 0 \end{split} $$ (7) 式中
$$\begin{split} &{R_{ij}} = \displaystyle\int_{{\varGamma _e}} {{\phi _j}\nabla G({\boldsymbol{M}},{{\boldsymbol{M}}_i}) \cdot {\boldsymbol{n}}{\text{d}}s - {\delta _{ij}}\lambda } \\ & {L_{iz}} = \displaystyle\int_{{\varGamma _z}} G({\boldsymbol{M}}, {{\boldsymbol{M}}_i}) {\text{d}}s \\ & {U_{imj}} = \displaystyle\iint_{{\varLambda _e}} {G({\boldsymbol{M}},{{\boldsymbol{M}}_i})\dfrac{{\partial {\phi _m}}}{{\partial x}}}\dfrac{{\partial {\phi _j}}}{{\partial x}}{\text{d}}A + G({\boldsymbol{M}},{{\boldsymbol{M}}_i}) \dfrac{{\partial {\phi _m}}}{{\partial y}} \dfrac{{\partial {\phi _j}}}{{\partial y}}{\text{d}}A\\ & {W_{imj}} = \displaystyle\iint_{{\varLambda _e}} {G({\boldsymbol{M}},{{\boldsymbol{M}}_i}){\phi _m}{\phi _j}}{\text{d}}A \end{split}$$ 则单元内的局部方程组的矩阵形式为
$$ {\left( {{B_{iz}}} \right)_{3 \times 3}}{{\boldsymbol{q}}_n} = - {\left( {{E_{ij}}} \right)_{3 \times 3}}{\boldsymbol{p}} - {\left( {{C_{ij}}} \right)_{3 \times 3}}\frac{{{\text{d}}{p_j}}}{{{\text{d}}t}} - {\left( {{F_i}} \right)_3} $$ (8) 式中,
${E_{ij}} = \displaystyle\sum\limits_{m = 1}^3 {\left( {{R_{ij}} - {U_{imj}}{\psi _m}} \right)}$ ,${B_{iz}} = {L_{iz}}{\left(\dfrac{1}{K}\right)_z}$ ,${C_{ij}} = \displaystyle\sum\limits_{m = 1}^3 {{W_{imj}}{\sigma _m}}$ ,${F_i} = \displaystyle\sum\limits_{j = 1}^3 {\displaystyle\sum\limits_{m = 1}^3 {{W_{imj}}{\upsilon _m}} } {f_j}$ , 在不含有裂缝单元或井筒的基质网格中(无源汇项), 此时$ {F_i} = 0 $ ,${\boldsymbol{p}} = {\left( {{p_1},{p_2},{p_3}} \right)^{\rm{T}}}$ ,$ {\boldsymbol{q}} = {\left( {{q_a},{q_b},{q_c}} \right)^{\rm{T}}} $ .计算矩阵B的逆矩阵, 将上式变形为
$$ {{\boldsymbol{q}}_n} = - {{\boldsymbol{B}}^{ - 1}}{\boldsymbol{Ep}} - {{\boldsymbol{B}}^{ - 1}}{\boldsymbol{C}}\frac{{{\text{d}}{\boldsymbol{p}}}}{{{\text{d}}t}} - {{\boldsymbol{B}}^{ - 1}}{\boldsymbol{F}} $$ (9) 利用权重系数θ控制显隐式程度, 得到
$$ \begin{split} {{\boldsymbol{q}}_n} = & - \left(\theta {{\boldsymbol{B}}^{ - 1}}{\boldsymbol{E }}+ \frac{{{{\boldsymbol{B}}^{ - 1}}{\boldsymbol{C}}}}{{\Delta t}}\right){{\boldsymbol{p}}^{\left( {n + 1} \right)}} + \\ & \left[ { - (1 - \theta ){{\boldsymbol{B}}^{ - 1}}{\boldsymbol{E }}+ \frac{{{{\boldsymbol{B}}^{ - 1}}{\boldsymbol{C}}}}{{\Delta t}}} \right]{{\boldsymbol{p}}^{\left( n \right)}} - {{\boldsymbol{B}}^{ - 1}}{\boldsymbol{F}} \end{split} $$ (10) 简写为
$$ {{\boldsymbol{q}}_n} = {\boldsymbol{M}}{{\boldsymbol{p}}^{\left( {n + 1} \right)}} + {\boldsymbol{b}} $$ (11) 式中
$$ \left\{\begin{split} &{\boldsymbol{M}} = - \left(\theta {{\boldsymbol{B}}^{ - 1}}{\boldsymbol{E}} + \dfrac{{{{\boldsymbol{B}}^{ - 1}}{\boldsymbol{C}}}}{{\Delta t}}\right)\\ & {\boldsymbol{b}}=\left[ { - (1 - \theta ){{\boldsymbol{B}}^{ - 1}}{\boldsymbol{E}} + \dfrac{{{{\boldsymbol{B}}^{ - 1}}{\boldsymbol{C}}}}{{\Delta t}}} \right] {{\boldsymbol{p}}^{\left( n \right)}} - {{\boldsymbol{B}}^{ - 1}}{\boldsymbol{F}} \end{split}\right.$$ 3. 改进格林元法与EDFM的耦合求解
基质与裂缝之间的窜流量被认为是基质单元中的源项或汇项, 因此有必要明确哪个单元中存在源项或汇项, 而这和基质与裂缝之间的几何性质有关. 本文所提出的改进格林元法可适用于任意单元类型的网格剖分, 如矩形单元和三角形单元, 其中三角形单元的前处理, 本文借助了开源程序Distmesh三角形非结构化网格剖分的Matlab代码[30-31]. 相比于矩形结构化网格, 三角形非结构化网格没有规则的拓扑结构, 网格节点的分布更具灵活性, 在诸多工程力学领域均有较好的应用[32-33]. 如图2所示的二维EDFM, 其中两条离散裂缝分别以红色和绿色的线段绘制, 通过基质单元的网格边界将与其相交的裂缝离散为多个裂缝单元, 并对每个裂缝单元进行编号.
针对渗流介质某基质网格内含有裂缝单元的情况, 可将源点取在裂缝网格中心, 由于该点是基质网格内点, 其特征角度为2π, 则得到相应的边界积分方程离散格式如下
$$ \begin{split} & \sum\limits_{j = 1}^3 {{R_{ij}}{p_f}} + \sum\limits_{z = a}^c {{L_{iz}}{{\left(\frac{q}{K}\right)}_z}} - \sum\limits_{j = 1}^3 {\sum\limits_{m = 1}^3 {{U_{imj}}{\psi _m}{p_j}} } + \\ &\qquad \sum\limits_{j = 1}^3 {\sum\limits_{m = 1}^3 {{W_{imj}}\left({\sigma _m}\frac{{{\text{d}}{p_f}}}{{{\text{d}}t}} + {\upsilon _m}{f_j}\right)} } = 0 \end{split} $$ (12) 式中
$$\begin{split} &{R_{ij}} = \displaystyle\int_{{\varGamma _e}} {{\phi _j}\nabla G({\boldsymbol{M}},{{\boldsymbol{M}}_{fi}}) \cdot {\boldsymbol{n}}{\text{d}}s - 1}\\ & {L_{iz}} = \displaystyle\int_{{\varGamma _z}} G({\boldsymbol{M}}, {{\boldsymbol{M}}_{fi}}) {\text{d}}s\\ & {U_{imj}} = \displaystyle\iint_{{\varLambda _e}} {G({\boldsymbol{M}},{{\boldsymbol{M}}_{fi}})\dfrac{{\partial {\phi _m}}}{{\partial x}}}\dfrac{{\partial {\phi _j}}}{{\partial x}}{\text{d}}A + G({\boldsymbol{M}},{{\boldsymbol{M}}_{fi}})\dfrac{{\partial {\phi _m}}}{{\partial y}} \dfrac{{\partial {\phi _j}}}{{\partial y}} {\text{d}}A\\ & {W_{imj}} = \displaystyle\iint_{{\varLambda _e}} {G({\boldsymbol{M}},{{\boldsymbol{M}}_{fi}}){\phi _m}{\phi _j}}{\text{d}}A \end{split}$$ 此时
${{\boldsymbol{M}}_{fi}}$ 是裂缝网格中心坐标.将式(12)变形后, 得到$$ {\left( {{B_{iz}}} \right)_{3 \times 3}}{{\boldsymbol{q}}_{{\rm{omf}}}} = {\left( {{R_{ij}}} \right)_{3 \times 3}}{{\boldsymbol{p}}_f} - {\left( {{{E'}_{ij}}} \right)_{3 \times 3}}{\boldsymbol{p}} - {\left( {{C_{ij}}} \right)_{3 \times 3}}\frac{{{\text{d}}{{\boldsymbol{p}}_f}}}{{{\text{d}}t}} - {\left( {{F_i}} \right)_3} $$ (13) 式中,
${E'_{ij}} = \displaystyle\sum\limits_{m = 1}^3 { {{U_{imj}}{\psi _m}} }$ , 即有$$ {{\boldsymbol{q}}_{{\rm{omf}}}} = {{\boldsymbol{B}}^{ - 1}}{\boldsymbol{R}}{{\boldsymbol{p}}_f} - {{\boldsymbol{B}}^{ - 1}}{\boldsymbol{Ep}} - {{\boldsymbol{B}}^{ - 1}}{\boldsymbol{C}}\frac{{{\text{d}}{{\boldsymbol{p}}_f}}}{{{\text{d}}t}} - {{\boldsymbol{B}}^{ - 1}}{\boldsymbol{F}} $$ (14) 式(14)即为考虑压力瞬态效应的基质网格与裂缝单元之间的传质表达式, 式中
$ {{\boldsymbol{B}}^{ - 1}} $ ,$ {\boldsymbol{R}} $ ,$ {\boldsymbol{E}} $ ,$ {\boldsymbol{C}} $ 仅仅与基质网格和裂缝单元的几何参数有关, 表征了裂缝单元分布在基质网格中的几何特征. 与传统线性分布假设不同, 该表达式中除了包含基质孔隙压力和裂缝孔隙压力, 还包含了裂缝孔隙压力对时间的导数项, 即体现了裂缝孔隙压力的瞬时变化率. 因此, 相较于传统两点线性流量估计格式, 可提高针对低渗透/致密储层中由于基质与裂缝两者间渗透率差异极大所产生的早期非稳态窜流量的数值计算精度.对于相邻的基质网格, 存在连续性条件: 共有的压力节点具有相同的压力, 共有的流量节点处的流量相同, 即两个相邻基质网格在共有流量节点处的外法线流量的代数和等于0, 这也是两套节点格林元法的基本原理. 如图3所示, 以两个相邻三角形基质单元为例, 介绍耦合过程的细节.
![]()
图 3 两套节点格林元方法中相邻单元方程组耦合示意图Figure 3. Sketch of coupling adjacent element equations in two sets of nodes-based GEM假设基质单元e1内的方程组为
$$ \left[ {\begin{array}{*{20}{c}} {q_a^{(2)}} \\ {q_b^{(2)}} \\ {q_c^{(2)}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {H_{11}^{e1}}&{H_{12}^{e1}}&{H_{13}^{e1}} \\ {H_{21}^{e1}}&{H_{22}^{e1}}&{H_{23}^{e1}} \\ {H_{31}^{e1}}&{H_{32}^{e1}}&{H_{33}^{e1}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {p_1^{(2)}} \\ {p_2^{(2)}} \\ {p_3^{(2)}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {b_1^{e1}} \\ {b_2^{e1}} \\ {b_3^{e1}} \end{array}} \right]{\text{ }} $$ (15) 单元e2内的方程组为
$$ \left[ {\begin{array}{*{20}{c}} { - q_a^{(2)}} \\ {q_d^{(2)}} \\ {q_e^{(2)}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {H_{11}^{e2}}&{H_{12}^{e2}}&{H_{13}^{e2}} \\ {H_{21}^{e2}}&{H_{22}^{e2}}&{H_{23}^{e2}} \\ {H_{31}^{e2}}&{H_{32}^{e2}}&{H_{33}^{e2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {p_2^{(2)}} \\ {p_3^{(2)}} \\ {p_4^{(2)}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {b_1^{e2}} \\ {b_2^{e2}} \\ {b_3^{e2}} \end{array}} \right] $$ (16) 则两个相邻基质单元方程组耦合为
$$ \begin{split} & \qquad\qquad \left[ \begin{gathered} 0 \\ q_b^{(2)} \\ q_c^{(2)} \\ q_d^{(2)} \\ q_e^{(2)} \\ \end{gathered} \right] - \left[ {\begin{array}{*{20}{c}} {b_1^{e1} + b_1^{e2}} \\ {b_2^{e1}} \\ {b_3^{e1}} \\ {b_2^{e2}} \\ {b_3^{e2}} \end{array}} \right] =\\ & \left[ {\begin{array}{*{20}{c}} {H_{11}^{e1}}&{H_{12}^{e1} + H_{12}^{e2}}&{H_{13}^{e1} + H_{13}^{e2}}&{H_{13}^{e2}} \\ {H_{21}^{e1}}&{H_{22}^{e1}}&{H_{23}^{e1}}&0 \\ {H_{31}^{e1}}&{H_{32}^{e1}}&{H_{33}^{e1}}&0 \\ 0&{H_{21}^{e2}}&{H_{22}^{e2}}&{H_{23}^{e2}} \\ 0&{H_{31}^{e2}}&{H_{32}^{e2}}&{H_{33}^{e2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {p_1^{(2)}} \\ {p_2^{(2)}} \\ {p_3^{(2)}} \\ {p_4^{(2)}} \end{array}} \right] \end{split} $$ (17) 此外, 对于离散裂缝网络中的流体流动问题, 可采用有限差分法对流动方程进行离散. 如图4所示, 单条裂缝可离散成多个裂缝单元, 离散裂缝交汇处的流量交换可采用星三角变换方法[34-35]进行处理.
![]()
图 4 裂缝交汇的星三角变换原理图Figure 4. Schematics of star-delta transformation for fracture segments intersecting得到相应的有限差分隐式格式, 共
$ {n_f} $ 个方程$$ \begin{split} {B'_i}p_{f,i - 1}^{\left( {n + 1} \right)} + &\left( {{{C'}_i} - \frac{{{{F'}_i}}}{{\Delta t}}} \right)p_{f,i}^{\left( {n + 1} \right)} + {D'_i}p_{f,i + 1}^{\left( {n + 1} \right)} + {E'_i}{q_{{\rm{omf}},i}} + \\ & {E'_i}{q_{{\text{well}}}} + \frac{{{{F'}_i}}}{{\Delta t}}p_{f,i}^{\left( n \right)} = 0 \end{split} $$ (18) 式中,
$ {B'_i} $ ,$ {C'_i} $ ,$ {D'_i} $ ,$ {E'_i} $ ,$ {F'_i} $ 为相关系数,${q_{{\rm{well}}}}$ 表示从裂缝单元流入井筒的流量, 对于不与井筒相连的裂缝单元,${q_{{\rm{well}}}} = 0$ , 未知量为$ {p_f} $ ; 对于与井筒相连的裂缝单元, 假设油井在恒定的井底压力$ {p_{wf}} $ 下进行生产, 此时$ {p_f} = {p_{wf}} $ , 未知量为${q_{{\rm{well}}}}$ . 由此, 上式结合井方程[36]即可以封闭求解.假设基质单元节点数为
$ {n_p} $ 、基质单元边的个数为${n_{{\rm{ed}}}}$ 、裂缝单元的个数为$ {n_f} $ . 由此可知, 全局方程组由两部分构成: 一部分是基于式(11), 采用式(17)耦合方法获取的表征基质网格间渗流的方程组; 第二部分是由式(18)构成的表征裂缝网格间流动的方程组, 且这两部分方程组中的$ {{\boldsymbol q}_{{\rm{omf}}}} $ 用式(14)表示.整体上, 全局方程组共有
$ {n_e} + {n_f} $ 个方程(其中第一部分方程有$ {n_e} $ 个, 第二部分方程有$ {n_f} $ 个)和$ {n_p} + {n_f} $ 个未知量(包含$ {n_p} $ 个基质网格节点压力和$ {n_f} $ 个来源于裂缝单元的未知量). 因为$ {n_e} $ 大于$ {n_p} $ , 可得知由改进格林元法获得的该全局方程组为超定方程组, 在封闭外边界条件下, 此时方程组可表示为$$ {\left( {{b_i}} \right)_{{n_e} + {n_f}}} = {\left( {{H_{ij}}} \right)_{\left( {{n_e} + {n_f}} \right) \times \left( {{n_p} + {n_f}} \right)}}{\left( {p_i^{\left( {n + 1} \right)}} \right)_{{n_p} + {n_f}}} $$ (19) 基于泛函分析中的正交投影定理, 该超定方程组的解可等价于计算式(20)的解. 由此, 所有的未知量均可获得.
$$ \left( {{{\boldsymbol{H}}^{\rm{T}}}{\boldsymbol{H}}} \right){{\boldsymbol{p}}^{\left( {n + 1} \right)}} = {{\boldsymbol{H}}^{\rm{T}}}{\boldsymbol{b}} $$ (20) 4. 模型验证及讨论
4.1 模型验证
将本文提出的两套节点格林元法EDFM与原始EDFM[6]及商业模拟软件tNavigator® LGR模块进行对比, 其中tNavigator®是由俄罗斯Rock Flow Dynamics (RFD)公司开发的高性能商业模拟软件, 其结果可以参考为精确解. 该实例是采用单水平井和单压裂缝的机理模型验证本文模型对长期产能预测的准确性, 低渗透储层物性参数值见表1. 图5展示了矩形储层域的示意图, 地层中心有一口生产井和一条压裂210 m的主裂缝, 水平井筒穿过裂缝中点. 对储层域进行离散化, tNavigator® LGR模块局部加密网格处理的裂缝必须沿网格线方向展布, 而在EDFM中裂缝走向不受限制. 水平井采用10 MPa恒定井底流压制度进行生产, 图6展示了三种不同模型计算1000天的压力分布场图, 图7展示了三种不同模型所获得的日产油量曲线. 结果表明, 三种不同模型的结果吻合较好, 验证了该方法对产能预测的整体精度.
表 1 数值算例的输入参数Table 1. Input parameters of numerical caseProperties Value Properties Value reservoir thickness/m 10 initial reservoir pressure/MPa 20 matrix porosity 0.15 matrix permeability/mD 1 rock compressibility coefficient/MPa−1 1.08 × 10−4 fracture aperture/m 0.005 fracture porosity 0.45 fracture permeability/mD 70000 oil viscosity/(mPa·s) 2 oil compressibility coefficient/MPa−1 3 × 10−3 ![]()
图 6 三种不同模型压力分布对比图(1000 d)Figure 6. Comparison of pressure distribution maps for three various simulators (1000 d)![]()
图 6 三种不同模型压力分布对比图(1000 d)(续)Figure 6. Comparison of pressure distribution maps for three various simulators (1000 d) (continued)![]()
图 7 三种不同模型产油速度的对比图(1000 d)Figure 7. Comparison of oil rate for three various simulators (1000 d)4.2 基质−裂缝间瞬态流动对比
该算例旨在验证两套节点格林元法EDFM能够较好地反映基质−裂缝间的瞬态流动, 算例的基本属性与表1中的参数值相同, 讨论中考虑了两个重要因素: 网格大小和基质渗透率, 将网格尺寸设置为2 m × 2 m和10 m × 10 m, 基质渗透率分别为0.1 mD, 1 mD和10 mD. 计算时间设置为10 d. 为了实现对早期瞬态流动的准确描述, 与上一个算例相比, 在模拟过程中使用了更小的时间步长. 图8展示了三种模型(网格步长10 m × 10 m)模拟10 d压力分布的比较, 从压力场的结果来看, 不同模型之间没有显著差异. 图9比较了在不同网格尺寸(2 m × 2 m和10 m × 10 m)和基质渗透率(0.1 mD, 1 mD, 10 mD)条件下, 由三种模型(原始EDFM、本文EDFM和tNavigator®)计算所得的日产油量曲线. 与tNavigator®结果相比, 从生产早期可以明显看出, 本文EDFM比原始EDFM具有更高的精度, 这是因为本文提出的两套节点格林元法能够有效地反映局部基质网格和裂缝单元之间的瞬态流动, 取代了原始EDFM的线性流动假设, 瞬态效应的程度与网格尺寸和基质渗透率密切相关. 平均相对误差可由式(21)获得, 计算结果如图10所示. 结果表明, 本文EDFM的平均相对误差比原始EDFM小得多. 同时, 随着基质渗透率的降低和网格尺寸的增大, 瞬态效应的影响增大, 平均相对误差增大. 本文所提出的两套节点格林元法实现了对基质−裂缝间瞬态流动更准确的表征, 能提高模拟的早期精度.
![]()
图 8 三种不同模型压力分布对比图(10 d)Figure 8. Comparison of pressure distribution maps for three various simulators (10 d)![]()
图 10 原始EDFM和修正EDFM之间的误差比较Figure 10. Error comparison between the original EDFM and the modified EDFM$$ {A_e} = \frac{1}{N}\sum\limits_{i = 1}^N {\left| {\frac{{{q_i} - q_i^{({\rm{ref}})}}}{{q_i^{({\rm{ref}})}}}} \right|} \times 100\text{%} $$ (21) 式中,
$ {q_i} $ 表示由EDFM计算所得的井筒流量, m3/d;$q_i^{\left( {{\rm{ref}}} \right)}$ 为由tNavigator®计算所得的井筒流量, m3/d; N为时间步数.4.3 复杂油藏边界−缝网−SRV分区模型
原始EDFM仅支持矩形结构化网格剖分, 在网格适用性方面存在局限性. 由于格林元法的引入, 本文改进后的EDFM可支持非匹配性的非结构化网格剖分, 能较完善地考虑复杂油藏边界−缝网及SRV分区的非常规油气储层体积压裂井产量分析的复杂问题. 本算例考虑基质渗透率为0.1 mD的致密储层, 人工裂缝渗透率为5000 mD, 其他参数值同表1, 图11(a)展示了储层域的基本模型, 图11(b)和图11(c)分别考虑了圆形和矩形的SRV改造区. COMSOL商业模拟软件以标准有限元(SFEM)为基础, 对DFM的模块化处理极为成熟. 本文以SFEM-COMSOL数值解为参考, 以图11(a)中的基本模型为例, 对本文模型进行验证, 不同模型的网格剖分对比见图12. 水平井采用10 MPa恒定井底流压制度进行生产, 图13对比了本文模型与SFEM-COMSOL两种模型模拟1000 d的日产油量曲线, 结果对比本文模型的解与参考解较为吻合, 验证了本文模型在考虑复杂油藏边界条件下产量预测的准确性.
![]()
图 12 不同模型的网格剖分对比图Figure 12. Comparison diagram of mesh division for different simulators![]()
图 12 不同模型的网格剖分对比图(续)Figure 12. Comparison diagram of mesh division for different simulators (continued)![]()
图 13 修正EDFM与COMSOL产油速度的对比图(1000 d)Figure 13. Comparison of oil rate between modified EDFM and COMSOL (1000 d)图14展示了本文模型对图11(a)的基本模型不同开发阶段得到的压力分布情况. 图15展示了本文模型对图11(b)和图11(c)两种不同SRV分区模拟100天的压力分布对比, 从图中能看出开发早期以缝网内部流体流动为主, 生产井产量完全由裂缝供给; 当考虑SRV分区时, 开发中期离散裂缝周围的内区流体以垂直于裂缝面的方向流入裂缝, 由于内外区基质渗透率级别存在差异, 开发后期压力波才慢慢向外区传播. 本文所提出的方法采用非稳态渗流控制方程的边界积分形式推导了基质网格与裂缝网格之间传质量的新格式, 因此也可以将EDFM的适用性拓展至油藏压裂井早期生产动态特征分析. 图16对比了不考虑分区、矩形SRV和圆形SRV三种SRV形态的生产动态特征, 其中矩形SRV对应的外区流动形态为复合线性流, 圆形SRV对应的为拟径向流, 不考虑分区的产量动态特征位于二者之间.
![]()
图 15 不同SRV分区模型的压力分布对比(100 d)Figure 15. Comparison of pressure distribution of different SRV zoning models (100 d)![]()
图 16 不同SRV对应的产量动态特征对比Figure 16. Comparison of dynamic production characteristics corresponding to different SRV5. 结 论
(1)本文提出了满足局部物质守恒且具有高精度、适用于任意网格剖分的压力流量两套节点的格林元方法, 能用于二阶偏微分方程的稳定求解; 将两套节点格林元法与EDFM进行耦合求解, 能更精确地获得瞬态压力和流量分布, 通过实例验证了本文模型的准确性.
(2)针对原始EDFM仅适用于矩形网格剖分、难于考虑复杂油藏边界的局限性, 本文提出的修正EDFM能适应任意多边形单元的网格剖分, 可适用于复杂油藏边界、缝网及SRV分区等复杂渗流模型的计算, 是对原始EDFM在网格适用性方面的拓展和升级.
-
![]()
图 1 压力和流量两套节点三角形单元示意图
Figure 1. Sketch of pressure and flux two sets of nodes in triangle cell
![]()
图 3 两套节点格林元方法中相邻单元方程组耦合示意图
Figure 3. Sketch of coupling adjacent element equations in two sets of nodes-based GEM
![]()
图 4 裂缝交汇的星三角变换原理图
Figure 4. Schematics of star-delta transformation for fracture segments intersecting
![]()
图 6 三种不同模型压力分布对比图(1000 d)
Figure 6. Comparison of pressure distribution maps for three various simulators (1000 d)
![]()
图 6 三种不同模型压力分布对比图(1000 d)(续)
Figure 6. Comparison of pressure distribution maps for three various simulators (1000 d) (continued)
![]()
图 7 三种不同模型产油速度的对比图(1000 d)
Figure 7. Comparison of oil rate for three various simulators (1000 d)
![]()
图 8 三种不同模型压力分布对比图(10 d)
Figure 8. Comparison of pressure distribution maps for three various simulators (10 d)
![]()
图 10 原始EDFM和修正EDFM之间的误差比较
Figure 10. Error comparison between the original EDFM and the modified EDFM
![]()
图 12 不同模型的网格剖分对比图
Figure 12. Comparison diagram of mesh division for different simulators
![]()
图 12 不同模型的网格剖分对比图(续)
Figure 12. Comparison diagram of mesh division for different simulators (continued)
![]()
图 13 修正EDFM与COMSOL产油速度的对比图(1000 d)
Figure 13. Comparison of oil rate between modified EDFM and COMSOL (1000 d)
![]()
图 15 不同SRV分区模型的压力分布对比(100 d)
Figure 15. Comparison of pressure distribution of different SRV zoning models (100 d)
![]()
图 16 不同SRV对应的产量动态特征对比
Figure 16. Comparison of dynamic production characteristics corresponding to different SRV
表 1 数值算例的输入参数
Table 1 Input parameters of numerical case
Properties Value Properties Value reservoir thickness/m 10 initial reservoir pressure/MPa 20 matrix porosity 0.15 matrix permeability/mD 1 rock compressibility coefficient/MPa−1 1.08 × 10−4 fracture aperture/m 0.005 fracture porosity 0.45 fracture permeability/mD 70000 oil viscosity/(mPa·s) 2 oil compressibility coefficient/MPa−1 3 × 10−3 
下载: 导出CSV
-
[1] 曹仁义, 程林松, 杜旭林等. 致密油藏渗流规律及数学模型研究进展. 西南石油大学学报(自然科学版), 2021, 43(5): 113-136 (Cao Renyi, Cheng Linsong, Du Xulin, et al. Research progress on fluids flow mechanism and mathematical model in tight oil reservoirs. Journal of Southwest Petroleum University (Science &Technology Edition) , 2021, 43(5): 113-136 (in Chinese) [2] Obembe AD, Hossain ME. A new pseudosteady triple-porosity model for naturally fractured shale gas reservoir//SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, 2015
[3] 苏皓, 雷征东, 李俊超等. 储集层多尺度裂缝高效数值模拟模型. 石油学报, 2019, 40(5): 587-593, 634 (Su Hao, Lei Zhendong, Li Junchao, et al. An effective numerical simulation model of multi-scale fractures in reservoir. Acta Petrolei Sinica, 2019, 40(5): 587-593, 634 (in Chinese) [4] Rao X, Cheng L, Cao R, et al. A modified projection-based embedded discrete fracture model (pEDFM) for practical and accurate numerical simulation of fractured reservoir. Journal of Petroleum Science and Engineering, 2020, 187: 106852
[5] 杜旭林, 程林松, 牛烺昱等. 基于XFEM-MBEM 的嵌入式离散裂缝模型流固耦合数值模拟方法. 力学学报, 2021, 53(12): 3416-3427 (Du Xulin, Cheng Linsong, Niu Langyu, et al. Numerical simulation for coupling flow and geomechanics in embedded discrete fracture model based on XFEM-MBEM. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(12): 3416-3427 (in Chinese) [6] Li L, Lee SH. Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media. SPE Reservoir Evaluation & Engineering, 2008, 11(4): 750-758
[7] Lee SH, Lough MF, Jensen CL. Hierarchical modeling of flow in naturally fractured formations with multiple length scales. Water Resources Research, 2001, 37(3): 443-455
[8] Hajibeygi H, Karvounis D, Jenny P. A hierarchical fracture model for the iterative multiscale finite volume method. Journal of Computational Physics, 2011, 230(24): 8729-8743
[9] 朱大伟, 胡永乐, 崔明月等. 局部网格加密与嵌入式离散裂缝模型耦合预测压裂改造井产能. 石油勘探与开发, 2020, 47(2): 341-348 (Zhu Dawei, Hu Yongle, Cui Mingyue, et al. Productivity simulation of hydraulically fractured wells based on hybrid local grid refinement and embedded discrete fracture model. Petroleum Exploration and Development, 2020, 47(2): 341-348 (in Chinese) doi: 10.11698/PED.2020.02.12 [10] Xue X, Yang C, Onishi T, et al. Modeling hydraulically fractured shale wells using the fast-marching method with local grid refinements and an embedded discrete fracture model. SPE Journal, 2019, 24(6): 2590-2608
[11] Liu L, Huang Z, Yao J, et al. An efficient hybrid model for 3D complex fractured vuggy reservoir simulation. SPE Journal, 2020, 25(2): 907-924
[12] Cao Y, Killough JE. An improved boundary element method for modeling fluid flow through fractured porous medium//SPE Reservoir Simulation Conference, Society of Petroleum Engineers, 2017
[13] Taigbenu AE. A more efficient implementation of the boundary element theory//Proc. 5th International Conference on Boundary Element Technology (BETECH 90), Newark, Delaware, 1990
[14] Taigbenu AE. The Green element method. International Journal for Numerical Methods in Engineering, 1995, 38: 2241-2263 doi: 10.1002/nme.1620381307
[15] Taigbenu AE, Onyejekwe OO. Green element simulations of the transient nonlinear unsaturated flow equation. Applied Mathematical Modelling, 1995, 19(11): 675-684
[16] Taigbenu AE, Onyejekwe OO. A mixed Green element formulation for the transient Burgers equation. International Journal for Numerical Methods in Fluids, 1997, 24(6): 563-578
[17] Taigbenu AE, Akpofure E. The Green Element Method. Springer, 1999
[18] Archer R, Horne RN, Onyejekwe OO. Petroleum reservoir engineering applications of the dual reciprocity boundary element method and the Green element method//Proceedings of Twenty First International Conference on Boundary Element Method Southampton, 1999
[19] Archer R. C1 continuous solutions from the Green element method using Overhauser elements. Applied Numerical Mathematics, 2006, 56(2): 222-229 doi: 10.1016/j.apnum.2005.04.001
[20] Pecher R, Harris SD, Knipe RJ, et al. New formulation of the Green element method to maintain its second-order accuracy in 2D/3D. Engineering Analysis with Boundary Elements, 2001, 25(3): 211-219
[21] Lorinczi P, Harris SD, Elliot L. Modified flux-vector–based Green element method for problems in steady-state anisotropic media-Generalisation to triangular elements. Engineering Analysis with Boundary Elements, 2011, 35: 495-498
[22] Lorinczi P, Harris SD, Elliott L. Unsteady flux-vector-based Green element method. Transport in Porous Media, 2011, 87(1): 207-228
[23] Lorinczi P, Harris SD, Elliott L. Modified flux-vector-based Green element method for problems in steady-state anisotropic media. Engineering Analysis with Boundary Elements, 2009, 33(3): 368-387
[24] Lorinczi P, Harris SD, Elliott L. Modelling of highly-heterogeneous media using a flux-vector-based Green element method. Engineering Analysis with Boundary Elements, 2006, 30(10): 818-833
[25] Taigbenu AE. The flux-correct Green element formulation for linear, nonlinear heat transport in heterogeneous media. Engineering Analysis with Boundary Elements, 2008, 32(1): 52-63
[26] Taigbenu AE. Enhancement of the accuracy of the Green element method: Application to potential problems. Engineering Analysis with Boundary Elements, 2012, 36(2): 125-136
[27] 方思冬, 程林松, 饶翔等. 一种改进格林元方法及在渗流问题中的应用. 计算力学学报, 2021, 38(6): 787-795 (Fang Sidong, Cheng Linsong, Rao Xiang, et al. An improved Green element method and its application in seepage problems. Chinese Journal of Computational Mechanics, 2021, 38(6): 787-795 (in Chinese) doi: 10.7511/jslx20200907001 [28] Rao X, Cheng L, Cao R, et al. A mimetic Green element method. Engineering Analysis with Boundary Elements, 2019, 99: 206-221
[29] Rao X, Cheng L, Cao R, et al. A novel Green element method based on two sets of nodes. Engineering Analysis with Boundary Elements, 2018, 91: 124-131
[30] Persson PO, Strang G. A simple mesh generator in MATLAB. SIAM Review, 2004, 46(2): 329-345
[31] Persson PO. Mesh generation for implicit geometries. [PhD Thesis]. Cambridge: MIT, 2004
[32] Yan CZ, Guo H, Tang ZC. Three-dimensional continuous-discrete pore-fracture mixed seepage model and hydromechanical coupling model to simulate rock fracture driven by fluid. Journal of Petroleum Science and Engineering, 2022, 215: 110510
[33] Yan CZ, Gao Y, Guo H. A FDEM based 3D discrete mixed seepage model for simulating fluid driven fracturing. Engineering Analysis with Boundary Elements, 2022, 140: 447-463
[34] Jia P, Cheng L, Huang S, et al. Transient behavior of complex fracture networks. Journal of Petroleum Science and Engineering, 2015, 132: 1-17
[35] 贾品, 程林松, 黄世军等. 水平井体积压裂缝网表征及流动耦合模型. 计算物理, 2015, 32(6): 685-692 (Jia Pin, Cheng Linsong, Huang Shijun, et al. Characterization of fracture network by volume fracturing in horizontal wells and coupled flow model. Chinese Journal of Computational Physics, 2015, 32(6): 685-692 doi: 10.3969/j.issn.1001-246X.2015.06.008 [36] Peaceman DW. Interpretation of well-block pressure in numerical reservoir simulation with nonsquare gridblocks and anisotropic permeability. SPE Journal, 1983, 18(3): 183-194
计量
- 文章访问数: 1100
- HTML全文浏览量: 398
- PDF下载量: 95
