摘要:缝洞型介质通常具有非均质性强、结构多尺度的特征. 传统数值方法在解决此类多尺度流动问题时, 难以兼顾计算精度与计算效率, 无法实际应用. 对此, 本文提出了多孔介质流体流动的多尺度分解法, 并应用于缝洞介质流动模拟, 能够大幅减小计算的复杂度, 同时可以通过控制均化程度控制计算精度. 该方法将求解空间分为若干个子空间的正交直和, 从而获得一个近线性的计算复杂度; 以分层计算的方式实现了快速计算, 另外这种方法是一种无网格方法, 具有较好的地层适应性. 同时, 采用离散缝洞模型简化缝洞结构, 进一步提高了计算效率. 详细阐述了基于多尺度分解法的多孔介质流体流动数值计算格式的建立, 重点介绍了如何在不同的层次上计算基函数. 数值结果表明, 本文提出的计算方法不仅能够准确捕捉多孔介质中的精细流动特征, 而且具有很高的计算效率, 是一种有效的流动模拟方法.
关键词: 缝洞型介质/
多尺度分解法/
流动模拟/
离散缝洞模型

English Abstract

全文HTML

--> --> -->

缝洞型介质广泛存在于自然界中, 过去几年间, 针对其开展的数值模拟技术引起了广泛关注, 但是缝洞型介质流动模拟给常规数值方法带来很大的挑战. 缝洞型介质除了本身具有的非均质性, 其中所包含的裂缝和溶洞横跨多个尺度. 如果要描述缝洞介质所有尺度上的流动特征, 则需要很精细的网格划分, 使用常规数值方法求解需要巨大的计算量, 即便使用超级计算机也难以满足需求. 因此, 常规数值方法难以描述流体所有尺度上的流动特征. 尺度升级方法被广泛视为能减小计算量的方法^[1,2], 但是由于其忽视了小尺度信息导致不能反映精细的介质非均质性. 近年来, 多尺度方法被应用于流动模拟^[3-7], 这种方法能够在粗网格上通过计算多尺度基函数捕捉小尺度特征, 从而在保证计算精度的同时减少计算量. 然而, 当粗网格和细网格尺度相差较大时, 多尺度方法在差值算子及约束算子的构造方面存在着困难. 多重网格方法是公认的求解偏微分方程的最快的方法之一^[8-11], 因此, 本文首先将多重网格法应用于非均质介质流动问题, 并进一步构建了一种新的缝洞型介质流动模拟的多重网格法(多尺度分解法).
虽然多重网格方法被拓展至求解多种偏微分方程, 但是当方程系数粗糙时, 会导致该方法收敛性差^[12,13]. 因此, 对于强非均质性油藏, 传统多重网格法也存在收敛性差的问题. 近年来, 虽然改进的多重网格方法在一定程度上能够缓和粗糙系数对收敛性的影响^[14-17], 但能够精确适应粗糙系数的多重网格法一直都是数学界的难题^[18].
目前, 主要有三种方法用于缝洞型碳酸盐岩数值模拟: 1)等效介质模型^[19]; 2)多重介质模型^[20]; 3)离散缝洞网络(DFVN)模型^[21,22]. DFVN模型能精细描述缝洞, 基岩和裂缝系统视为渗流区域, 溶洞系统视为自由流区域, 是离散裂缝模型的有效延伸. DFVN模型属于宏观尺度的精细模型, 能克服三重介质模型和等效介质模型的瓶颈, 但受目前计算机硬件和数值模拟技术的制约, 只能处理小规模的缝洞型油藏. 针对此问题, 对于小尺度裂缝和溶洞, 对其进行等效处理, 用等效后的渗透率场描述其造成的非均质性; 使用离散裂缝模型描述大尺度裂缝, 准确刻画裂缝流动状态.
本文基于多尺度分解思想^[23,24], 构建了缝洞型介质多重网格计算格式, 能够准确描述介质的非均质性, 并能准确模拟大尺度裂缝存在时的压力差分布. 文中首先介绍了离散缝洞的数学模型和多尺度分解法的基本原理, 然后构建了缝洞介质的多尺度分解计算格式, 并给出了几种不同的算例. 数值结果表明, 多尺度分解法可以在充分反映裂缝溶洞影响的同时大幅减少计算复杂度, 是一种潜力很大的流动模拟数值方法.

2

-->

2.1.缝洞介质数学模型

本文通过计算等效渗透率来描述介质中的小尺度裂缝和溶洞导致的非均质性; 采用离散裂缝模型表征大尺度裂缝, 两者结合形成多尺度缝洞混合模型(图1).
图 1 缝洞型介质示意图
Figure1. Schematic of fractured-vuggy porous medium.

将基岩和裂缝系统视为渗流区域, 其运动方程符合Darcy定律, 即

$\mu {{{K}}_l}^{ - 1} \cdot {{{v}}_l} + \nabla p - \rho {{f}} = 0, $

式中$\mu$为流体黏度(mPa·s); K_i为渗透率(${\text{μ}}{\rm m}^2$); v_i为渗流速度(m/s); 对于岩块l = m, 对于裂缝l = f; p为多孔介质中的平均压力(Pa); $\rho$是流体密度(kg/m³); f为单位质量力(m/s²).
溶洞的空间尺度较大, 其流动模式视为自由流, 为简单起见本文仅考虑牛顿流体的低雷诺数运动. 为了避免引入复杂的耦合边界条件, 降低自由流-渗流耦合问题的复杂性, 采用Brinkman方程描述溶洞内的自由流:

$\mu {{K}}_v^{ - 1} \cdot {{{v}}_v} - \frac{\mu }{\phi }{\nabla ^2}{{{v}}_v} + \nabla {p_v} - \rho {{f}} = 0, $

式中v_v是流体真实速度(m/s); K_v是溶洞内渗透率(${\text{μ}}{\rm m}^2$); $\phi $为溶洞内孔隙度(无充填时为1); p_v为流体压力(Pa).
假设流体流动过程温度不变, 综合以上方程可以推导出DFVN模型的控制方程组:

$\left\{ \begin{aligned}&\mu {{{K}}_l}^{ - 1} \cdot {{v}} - B\frac{\mu }{\phi }{\nabla ^2}{{v}} + \nabla p - \rho {{f}} = 0, \\&\nabla \cdot {{v}} = 0, \end{aligned} \right.$

式中
${{{K}}_l} = \left\{ \begin{aligned}&{{{K}}_{\rm{m}}}, \;{\text{在岩块系统中}}, \\&{{{K}}_{\rm{f}}}, \;{\text{在裂缝系统中}}, \\&{{{K}}_{\rm{v}}}, \;{\text{在溶洞系统中}};\end{aligned} \right.$
$ B =\left\{ \begin{aligned}&0, \; {\text{在岩块和裂缝系统中}}, \\&1, \; {\text{在溶洞系统中}};\end{aligned} \right.$
${{v}}$为速度场. 结合合适的边界条件可得到完整的离散缝洞网络单相稳态渗流数学模型.
2 -->

2.2.等效渗透率计算

基于流量等效计算等效渗透率的方法, 往往仅考虑目标单元内的非均质性, 这样得到的等效渗透会存在误差. 对此, 本文采用超样本方法, 在包含目标单元的更大区域内建立离散缝洞网络模型并进行有限元数值求解, 从而获取超样本单元内的压力和速度场; 最后, 基于体积平均法计算目标网格内的等效渗透率.
在超样本单元内结合封闭定压边界条件建立离散缝洞网络模型, 采用有限单元法进行求解^[21], 根据求得的压力及速度场, 采用体积平均法计算目标网格上的速度和压力梯度的体积平均值^[25]:

${\left\langle u \right\rangle ^j} = \frac{1}{{{V_{\rm{b}}}}}\int_V {{u^j}{\rm{d}}V} = \frac{1}{{{V_{\rm{b}}}}}\sum\limits_{l = 1}^N {{u_l}{V_l}} , $

${\left\langle {\nabla p} \right\rangle ^j} = \frac{1}{{{V_{\rm{b}}}}}\int_V {{{\left( {\nabla p} \right)}^j}{\rm{d}}V} = \frac{1}{{{V_{\rm{b}}}}}\sum\limits_{l = 1}^N {{{\left( {\nabla p} \right)}_l}{V_l}} , $

式中$j = x, y$(表示定压边界施加的坐标轴方向); ${V_{\rm{b}}}$表示目标原始网格的体积; ${u_l}$和${\left( {\nabla p} \right)_l}$分别表示目标原始网格中第$l$个单元的速度和压力梯度; ${V_l}$表示第$l$个单元的体积; $N$表示目标原始网格中单元的总数. 在计算得到${\left\langle u \right\rangle ^j}$和${\left\langle {\nabla p} \right\rangle ^j}$后, 根据达西定律, 可以求得等效渗透率, 详细计算过程可以参考文献[26].

多尺度分解法的目的是尽可能快的、以近线性的计算复杂度求解偏微分方程. 令q为网格层数, 定义${{\cal{I}}^{\left( q \right)}}$为每一层的细网格单元索引$i =$$ \left( {{i_1}, \cdots, {i_q}} \right)$. 令${{{i}}^{\left( k \right)}} = \left( {{i_1}, \cdots, {i_k}} \right)$, ${{\cal{I}}^{\left( k \right)}}{{ = }}\left\{ {{{{i}}^{\left( k \right)}}:{{i}} \in {{\cal{I}}^{\left( q \right)}}} \right\}$, 其中1 ≤ k ≤ q, ${{i}} = \left( {{i_1}, \cdots, {i_q}} \right) \in {{\cal{I}}^{\left( q \right)}}$; ${{{i}}^{\left( {k, {k'}} \right)}}$为满足$j \in {{\cal{I}}^{\left( {{k'}} \right)}}$和${j^{\left( k \right)}} = i$的单元, 其中${k'} \in \left\{ {k, \cdots, q} \right\}$. ${{ \pi} ^{\left( {k, k + 1} \right)}}$为${{\cal{I}}^{\left( k \right)}} \times {{\cal{I}}^{\left( {k + 1} \right)}}$矩阵, 满足当$j \notin {i^{\left( {k, k + 1} \right)}}$时${ \pi} _{i, j}^{\left( {k, k + 1} \right)} = 0$, 同时${{ \pi} ^{\left( {k, k + 1} \right)}}{\left( {{{ \pi} ^{\left( {k, k + 1} \right)}}} \right)^{\rm{T}}}{{ = }}{{{I}}^{\left( k \right)}}$, 其中${{{I}}^{\left( k \right)}}$为${{\cal{I}}^{\left( k \right)}} \times {{\cal{I}}^{\left( k \right)}}$的单位矩阵.
不同网格层上的测试函数定义为

$\phi _i^{\left( k \right)} = \sum\limits_{j \in {I^{\left( {k + 1} \right)}}} {{\pi} _j^{\left( {k, k + 1} \right)}\phi _j^{\left( {k + 1} \right)}} .$

如图2所示, $\varOmega $被分为若干${2^{ - k}} \times {2^{ - k}}$的子单元${\left( {\tau _i^{\left( k \right)}} \right)_{i \in {I^{\left( k \right)}}}}$. 令$\phi _i^{\left( k \right)} = 1_{\tau _i^{\left( k \right)}}^{}/\sqrt {\left| {\tau _i^{\left( k \right)}} \right|} $, 其中$\left| {\tau _i^{\left( k \right)}} \right|$为$\tau _i^{\left( k \right)}$的体积, $1_{\tau _i^{\left( k \right)}}^{}$为$\tau _i^{\left( k \right)}$的指示函数. 对于边界条件问题, 使用非零边界的测试函数. 通过假设A和B在进行一个博弈游戏来定义基函数, 具体理论证明见文献[24]. 根据博弈理论^[27,28]并结合流动方程, B的最佳选择为
图 2 区域$\varOmega $的网格剖分示意图
Figure2. Schematic of grid partition of solution space.

$p ^{\left( k \right)}\left( x \right){{ = }}\sum\limits_{i \in {{\cal{I}}^{\left( k \right)}}} {\psi _i^{\left( k \right)}\left( x \right)} \int\limits_\varOmega {p\left( y \right)\phi _i^{\left( k \right)}\left( y \right)} \;{\rm{d}}y, $

其中$\psi _i^{\left( k \right)}$为基函数, 定义为

$\psi _i^{\left( k \right)} = \mathbb{E}\left[\!{v|\int_\varOmega \!{v\left( y \right)\phi _j^{\left( k \right)}\left( y \right){\rm{d}}y = {\delta _{i, j}}, \;\;i \!\in\! {{\cal{I}}^{\left( k \right)}}} } \!\right]\!.$

在实际计算中, 可以通过求解以下方程获得

$\left\{ \begin{aligned}&{\text{求最小值}}{\left\| \psi \right\|_\lambda }, \\&{\text{受限于}}\psi \in H_0^1\left( \varOmega \right){\text{和}}\int_\varOmega {{\phi _j}\psi = {\delta _{ij}}, \;\;j = 1, \cdots, m, } \end{aligned} \right.$

其中${\left\| \psi \right\|_\lambda }{{ = }}\int_\varOmega {\nabla {\psi ^{\rm{T}}}} \lambda \nabla \psi \;{\rm{d}}x$为能量范数. 利用以下定理证明(9)式可以用于基函数的计算.
可以证明以下定理成立:
1) 最优问题(9)式可以确定唯一一个由方程(8)定义的${\psi _i}$;
2) 对于$j \in \left\{ {1, \cdots, m} \right\}$, $\sum\nolimits_{i = 1}^m {{w_i}{\psi _i}} $是${\left\| \psi \right\|_\lambda }$的最小值, 满足$\int_\varOmega {\psi \left( x \right){\phi _j}\left( x \right) = {w_j}} $;
3) ${\left\langle {{\psi _i}, {\psi _j}} \right\rangle _\lambda } \!\!=\! \varTheta _{i, j}^{ - 1}$, 其中$ {\left\langle {v, w} \right\rangle _\lambda } \!=\!\!\int_\varOmega {{{\left( {\nabla v} \right)}^{\rm{T}}}\lambda \nabla w}, $ $\varTheta _{i, j} = \int_{{\varOmega ^2}} {{\phi _i}\left( x \right)G\left( {x, y} \right){\phi _j}\left( y \right)} {\rm{d}}x{\rm{d}}y$, $G\left( {x, y} \right)$为协方差函数.
以均匀介质为例, 首次给出了流动问题的多重网格法的基函数$\psi _i^{\left( k \right)}$, 如图3所示, 每一层的基函数的共同特点是这些基函数都是局部化的. 本文没有对基函数进行截取, 当基函数包含所有信息时可视为精确的基函数, 单项流动中也存在其他高精度基函数^[29]. 局部化的思想是目前的热点问题之一, 该思想有利于信息的集中处理, 目前大数据等研究中均有涉及.
图 3 分层基函数示意图
Figure3. Illustration of gamblets $\psi _i^{\left( k \right)}$.

定义空间

${II^{\left( k \right)}} = {\rm{span}}\left\{ {\psi _i^{\left( k \right)}|i \in {{\cal{I}}^{\left( k \right)}}} \right\}, $

其中$k \in \left\{ {1, \cdots, q} \right\}$. 从(6)式可以推出此空间是嵌套的. 因此, 对于$k \in \left\{ {1, \cdots, q - 1} \right\}$,

$\psi _i^{\left( k \right)} = \sum\limits_{j \in {{\cal{I}}^{\left( {k{{ + }}1} \right)}}} {R_{i, j}^{\left( {k, k + 1} \right)}} \psi _j^{\left( {k + 1} \right)}, $

其中, ${{R}}^{\left( {k, k + 1} \right)}$为${{\cal{I}}^{\left( k \right)}} \times {{\cal{I}}^{\left( {k{{ + }}1} \right)}}$的约束矩阵,

$\begin{split}&R_{i, j}^{\left( {k, k + 1} \right)}=\int_\varOmega {\psi _i^{\left( k \right)}\phi _j^{\left( {k + 1} \right)}}\\ = &\;\mathbb{E}\!\left[\!\int_\varOmega \!\!v\left( y \right)\phi _j^{\left( {k + 1} \right)}\!\left( y \right) |\!\int_\varOmega \!v\!\left( y \right)\phi _l^{\left( k \right)}\left( y \right){\rm{d}}y \!=\! {\delta _{i, l}}, l \!\in\! {\cal{I}}^{(k)}\!\right]\!.\end{split}$

令${{A}} ^{\left( k \right)}$是刚度矩阵$A_{i, j}^{\left( k \right)}{{ = }}{\int_\varOmega {\left( {\nabla \psi _i^{\left( k \right)}} \right)} ^{\rm{T}}}\lambda \;\nabla \psi _i^{\left( k \right)}$, 那么

${{A}} ^{\left( k \right)} = {{{R}}^{\left( {k, k + 1} \right)}}{{{A}}^{\left( {k + 1} \right)}}{{{R}}^{\left( {k + 1, k} \right)}}.$

至此, 以上嵌套公式足以进行多重网格计算, 但是本文进一步对其进行多尺度分解. 定义${{\cal{J}}^{\left( k \right)}}$为k-元组, 其中${{j}} = \left( {{j_1}, \cdots, {j_k}} \right)$, $\left\{ {{j^{\left( {k - 1} \right)}}|{{j}} \in {{\cal{J}}^{\left( k \right)}}} \right\} =$$ {{\cal{I}}^{\left( {k{{ - }}1} \right)}}$. 令${{{W}}^{\left( k \right)}}$为${{\cal{J}}^{\left( k \right)}} \times {{\cal{I}}^{\left( k \right)}}$矩阵, 满足当$ {j^{\left( {k - 1} \right)}} \ne$$ {i^{\left( {k - 1} \right)}}$时$W_{j, i}^{\left( k \right)} = 0$, 对于$s \in {{\cal{I}}^{\left( {k - 1} \right)}}$, $ t \in \left\{1, \cdots, {m_s} -1\right\},$ ${t'} \in \left\{ {1, \cdots, {m_s}} \right\}$, $W_{\left( {s, t} \right), \left( {s, {t'}} \right)}^{\left( k \right)} = {\delta _{t, {t'}}} - {\delta _{t + 1, {t'}}}$. 令

$\chi _i^{\left( k \right)} = \sum\limits_{j = {{\cal{I}}^{\left( k \right)}}} {W_{i, j}^{\left( k \right)}} \psi _j^{\left( k \right)}, $

其中$i \in {{\cal{J}}^{\left( k \right)}}$, 定义空间

${III^{\left( k \right)}} = {\rm span}\left\{ {\chi _i^{\left( k \right)}|i \in {{\cal{I}}^{\left( k \right)}}} \right\}, $

${III^{\left( k \right)}}$为${II^{\left( k-1 \right)}}$关于${II^{\left( k \right)}}$的正交补集. 令${ \oplus _\lambda }$表示正交直和, 那么可以得到空间的多尺度分解为

${II^{\left( k \right)}} = {II^{\left( 1 \right)}}{ \oplus _\lambda }{III^{\left( 2 \right)}}{ \oplus _\lambda } \cdot \cdot \cdot { \oplus _\lambda }{III^{\left( k \right)}}.$

进一步可以推导出

$p = {p^{\left( 1 \right)}} + \left( {{p^{\left( 2 \right)}} - {p^{\left( 1 \right)}}} \right) + \cdot \cdot \cdot + \left( {{p^{\left( q \right)}} - {p^{\left( {q - 1} \right)}}} \right), $

此为压力p的正交分解, ${p^{\left( {k + 1} \right)}} - {p^{\left( k \right)}}$为流动方程在空间${III^{\left( {k{\rm{ + }}1} \right)}}$里的解. 流动方程的$\chi _i^{\left( k \right)}$及多尺度分解如图4和图5所示, $\chi _i^{\left( k \right)}$可以理解为在每一层上的基函数.
图 4 基函数$\chi _i^{\left( k \right)}$示意图
Figure4. Illustration of basis functions $\chi _i^{\left( k \right)}$.

图 5 区域$\varOmega $的多尺度解示意图
Figure5. Illustration of solution of the multiple spaces.

当介质内存在大尺度裂缝时, 则采用离散裂缝模型表征大裂缝, 即对裂缝进行降维处理^[8]. 假设裂缝介质流动方程为FEQ, 离散裂缝模型的计算格式为

$\int_\varOmega {{\rm{FEQ}}\;} {\rm{d}}\varOmega = \int_{{\varOmega _{\rm{m}}}} {{\rm{FEQ}}\;} {\rm{d}}{\varOmega _{\rm{m}}}\; + \sum\limits_i {{d_i}} \times \int_{{\varOmega _{\rm{f}}}} {{\rm{FEQ}}\;} {\rm{d}}{\varOmega _{{\rm{f}}, i}}, $

其中d_i为裂缝开度, ${\varOmega _{\rm{m}}}$为基岩区域, ${\varOmega _{\rm{f}}}$为裂缝区域. 在每一个子空间内求解离散裂缝模型, 捕捉裂缝对渗流的影响, 然后将各个子空间的解加起来获取最终缝洞介质数值解, 即

$\begin{split}{p_f} \;= \;& p_f^{\left( 1 \right)} + \left( {p_f^{\left( 2 \right)} - p_f^{\left( 1 \right)}} \right) + \cdot \cdot \cdot +\\& \left( {p_f^{\left( q \right)} - p_f^{\left( {q - 1} \right)}} \right) + \left( {{p_f} - p_f^{\left( q \right)}} \right).\end{split}$

本节首先通过一个数值算例与参考算例的对比验证了程序的正确性, 把使用传统模拟有限差分法求得的数值解作为参考解; 然后通过复杂离散裂缝模型数值算例进一步验证了方法的正确性和程序的鲁棒性. 对于单相流动而言, 着重对比压力差分布, 为了实现压力差的对比, 设介质中一点的压力为已知, 令压力分布为正值.
2

-->

4.1.小尺度缝洞模型算例

考虑图6(a)所示的缝洞型介质物理模型, 模型大小为100 m × 100 m. 模型包含小裂缝和溶洞, 不存在大裂缝, 模型中的参数如表1所列. 溶洞简化为规则球形, 溶洞内不考虑充填(渗透率无充填时为无穷大), 其中的流动视为自由流. 图6(b)是通过均化理论求得的等效渗透率场图. 令流体自上而下流动, 边界为不流动边界.

参数名称	参数值
裂缝渗透率/${\text{μ}}{\rm m}^2$	1 × 104
裂缝开度/m	1 × 10–3
流体黏度/mPa·s	1
基岩渗透率/${\text{μ}}{\rm m}^2$	0.001

表1裂缝性介质模型基本参数
Table1.Parameters of fractured porous medium.

图 6 (a)小尺度缝洞型油藏几何模型; (b)等效渗透率场图; (c) 多层网格系统
Figure6. (a) Geometrical model of a small-scale-fractured vuggy porous oil reservoir; (b) equivalent permeability of fractured-vuggy medium; (c) nested grid system.

令q = 4, 即网格系统包含4层网格, 如图6(c)所示. 采用多尺度分解法求解得到每一层上的解, 并根据(17)式获得最终的压力分布. 小尺度裂缝与溶洞拥有较高的渗透率, 会造成多孔介质的非均质性, 进而影响整体压力分布. 在数值算例中采用L²范数计算压力误差

${p_v} = \frac{{\left\| {{p_{\rm{f}}} - {p_{{\rm{mg}}}}} \right\|_2^2}}{{\left\| {{p_{\rm{f}}}} \right\|_2^2}}, $

式中${p_{\rm{f}}}$为参考压力解, ${p_{{\rm{mg}}}}$ 为多尺度分解法所计算的压力解.
图7给出了参考解和多重网格解的对比. 可以看出, 当k = 3时多重网格法所得压力分布与参考解几乎一致, 因此多重网格解能够精确地反映非均质性, 从而捕捉到小尺度裂缝和溶洞的影响. 表2列出了不同k时的计算误差, 其中当k = 4时为多层分解法的精确解, 当k = 3和k = 2时为均化计算, 计算误差均很小, 体现了该方法的精确性. 小尺度缝洞模型的计算结果表明, 此方法可以有效地进行非均质地层以及小尺度缝洞介质的流动模拟.
图 7 对于小尺度缝洞模型, 参考解和多重网格解对比　(a) 参考解; (b) k = 3时的多重网格解
Figure7. Comparison of reference solution and gamblets solution for a small-scale-fractured vuggy porous medium: (a) Reference solution; (b) gamblets solution with k = 3.

	k
	2	3	4
计算误差	0.0823	0.0132	3.2549 × 10–15

表2对于小尺度缝洞模型, 不同k时的计算误差
Table2.Relative error in different k for a small-scale-fractured vuggy porous medium.

2 -->

4.2.长裂缝模型算例

算例2考虑存在长裂缝的情况. 考虑如图8所示的裂缝型介质, 在小裂缝和溶洞的背景下包含一条长裂缝. 裂缝开度d_f = 1 × 10^–3 m, 裂缝渗透率为$d_{\rm f}^2/12 $ ${\text{μ}}{\rm m}^2$, 基岩孔隙度$\phi =0.2$. 边界均为不流动边界, 流体从上往下流动. 采用离散裂缝模型表征裂缝.
图 8 长裂缝介质模型
Figure8. Geometrical model of a fractured-vuggy porous medium with a long fracture.

本算例考虑长裂缝与小裂缝溶洞共存的情况. 压力分布图9展示了不同k时多尺度分解法与参考解的对比. 当k = 3时计算误差为0.0131, 当k = 2时计算误差为0.0741, 在进行均化的前提下依然能够反映裂缝的存在, 保证较高的计算精度. 参考解与多重网格解的对比结果说明, 本文构建的多尺度分解法能够有效地处理长裂缝的情况, 并能通过均化大幅减少计算量, 同时保持较高的计算精度.
图 9 对于长裂缝模型, 参考解和多重网格解对比　(a)参考解; (b) k = 3时的多重网络解; (c) k = 2时的多重网格解
Figure9. Comparison of reference solution and gamblets solution for a fractured-vuggy porous medium with a long fracture: (a) Reference solution; (b) gamblets solution with k = 3; (c) gamblets solution with k = 2.

2 -->

4.3.大尺度缝洞模型算例

自然介质中, 除了小裂缝, 往往伴随着压裂等增产措施产生大裂缝, 该算例检验本文提出的方法对大裂缝的模拟能力. 图10所示为10 m × 10 m缝洞型介质, 其中包含的小尺度裂缝和溶洞的参数与算例1相同. 在包含小尺度溶洞的同时, 还包含着尺度较大的裂缝网络系统. 裂缝网络由6条长裂缝组成, 长裂缝开度d_f = 1 × 10^–3 m, 裂缝渗透率为$ d_{\rm f}^2/12$ ${\text{μ}}{\rm m}^2 $, 基岩孔隙度$\phi = 0.2$. 流体从上至下流动.
图 10 大尺度缝洞介质几何模型
Figure10. Geometrical model of a large-scale-fractured vuggy porous medium.

多孔介质中的长裂缝作为导流通道, 对压力分布有显著的影响, 并且裂缝之间也会相互影响. 该方法在每一个子空间内求解离散裂缝模型, 图11为离散裂缝模型在各层上的解, 可以很明显地看出, 各个层次上的解都能反映出裂缝的存在. 然后将k个层次上的解加起来, 即得到最终的解.
图 11 大尺度缝洞介质模型在各层上的解
Figure11. Solutions in different levels for a large-scale-fractured vuggy porous medium.

图12给出了参考解和多重网格解的对比, 可以看出, 本文构建的数值算法在捕捉小尺度裂缝和溶洞的同时, 能够精确模拟长裂缝对压力场的影响, 并反映裂缝间的相互作用. 结合表3可知, 当均化程度较大(k = 1)时存在较大误差. 图13给出了沿x = 5 m的参考解与不同k时多重网格解的对比. 裂缝在流动过程中可视为一个等势体, 对应于图11中较平缓的部分, 与实际符合, 随着均化程度加大误差变大. 该算例展示了本文构建的数值方法对长裂缝的精确模拟能力.
图 12 对于大尺度缝洞介质模型, 参考解和多重网格解对比　(a) 参考解; (b) k = 3时的多重网格解; (c) k = 2时的多重网格解; (d) k = 1时的多重网格解
Figure12. Comparison of reference solution and gamblets solution for a large-scale-fractured vuggy porous medium: (a) Reference solution; (b) gamblets solution with k = 3; (c) gamblets solution with k = 2; (d) gamblets solution with k = 1.

实际上, 对于任意形状的一个裂缝性介质, 一旦求出刚度矩阵及对应的W和${\text{π}}$, 就可以有效地进行任意形状的裂缝性介质数值模拟. 因此, 该方法能够处理形状复杂的裂缝系统. 但是, 对于裂缝分布极其复杂的裂缝系统, 尤其裂缝之间距离非常小的情况, 多层网格剖分会存在一定困难, 造成分层的基函数不适应高渗透的裂缝, 进而影响该方法的计算精度.

	k
	1	2	3
计算误差	0.1925	0.0831	0.0141

表3对于大尺度缝洞介质模型, 不同k时的计算误差
Table3.Relative error with different k for a large-scale-fractured vuggy porous medium.

图 13 沿x = 5 m的参考解和不同均化程度下多重网格解对比曲线图
Figure13. Comparison of reference solution and gamblets solution with different k when x = 5 m.

1)随着现代地质建模技术的发展, 缝洞型地质模型可能包含数百万甚至数亿个网格, 采用传统的数值方法对其进行求解, 计算量巨大, 甚至超出当今计算机的计算能力. 本文首次提出缝洞介质的多尺度分解计算模式, 构建局部化的基函数, 该方法能够快速求解流动方程, 非常适合对强非均质型油藏进行精细流动模拟;
2)结合离散裂缝模型和等效介质模型, 既可以精确反映大尺度裂缝的流体流动, 又可以捕捉小尺度溶洞和裂缝的影响, 从而可以对缝洞型介质进行高效流动模拟;
3)每一层网格的基函数是独立计算的, 可以采用并行计算得到, 进一步减少了计算量. 因此, 本文提出的方法对于缝洞型油藏数值模拟有很高的潜在价值. 渗透率场对计算精度有一定影响, 因此提高算法的鲁棒性是研究重点. 下一步可采用压力梯度对比的形式验证方法对多相流动计算的精确性^[30].