黑油模型广泛用于模拟油藏中的三相流动, 而EDFM可以准确模拟裂缝的几何结构, 两者结合可以很好地模拟裂缝性油藏中的水驱过程^[18]。使用EDFM模拟复杂的裂缝系统已经被证明是可行的^[19]。通过将裂缝离散为小的分段裂缝并嵌入到待计算区域中, EDFM融合了双孔双渗方法和离散裂缝方法的优势。因此, EDFM可以避免非结构网格带来的高计算成本, 同时又能保留油藏模拟器本身的优势。

流线模拟技术已经被广泛地应用于计算油藏流体流动轨迹和流动时间。Muskat^[21]首先利用流线和流管研究流体运移的问题, 在此基础上, 许多研究人员使用流线方法来模拟简单的水驱问题^{[26, 27, 28]}, 并对二维空间不可压缩流体适用的势函数和流线函数进行定义。流线的“ 飞行时间” 作为空间变量引入之后^[29], 饱和度的运移方程便可以沿着流线进行数值求解, 从此流线模型逐渐代替了流管模型。流线算法随着非常规油气储集层中裂缝模拟技术的进步而逐步改进。Pollock^[23]建立了针对矩形网格计算流线轨迹和流动时间的经典方法, 随后该方法被扩展到复杂的三维油藏结构中^{[30, 31]}, 进而通过引入基于角点网格的拟飞行时间对该方法做进一步简化^[32]。对于流线方法在非结构网格、断层网格和局部网格中的应用, 有许多学者通过引入局部边界层的方法来处理交界面上的流速不连续问题^{[33, 34, 35]}。

流线模拟在天然裂缝性油藏中的应用存在2个棘手的问题：①由于裂缝的存在, 流线模拟需要计算两侧具有不同精度交界面上的流线轨迹; ②油藏模拟器计算得到的基质和裂缝网格上的流量对于流线模拟来说并不够。新的局部边界层的方法可以解决第1个问题^{[32, 33, 34, 35]}。为了解决第2个问题, 本文利用已知基质和裂缝流量, 基于质量守恒定律来计算三维网格6个面的流量。

二维情形下只有1条裂缝的基质网格中体积流量的关系见图1。基质被裂缝分成2个小的基质网格, 假设基质网格的流量为f₁— f₆、裂缝网格的流量为f₇— f₁₀, 其中f₇和f₈是基质网格和裂缝网格共享的。假设f₃+f₄+f₉=F₁、f₅+f₆+f₁₀=F₂、f₇+f₈=F₃, 引入EDFM方法的油藏模拟器可以计算基质网格每个面上的流量（f₁、f₂、F₁、F₂）、裂缝网格上的流量（f₉、f₁₀）以及裂缝和基质网格之间的总流量F₃。对于流线模拟来说, 每个新网格所有6个面的流量也需要额外进行计算。

图1

	Figure Option ViewDownloadNew Window
	图1 二维网格中裂缝流量与基质流量关系图

在本文研究的所有案例中, 最下方一面裂缝的流量占那一侧流量的大部分, 因此可以假设f₃=f₄、f₅=f₆, 联立所有流量关系得到如下方程组：

$\left\{ \begin{matrix} {{f}_{3}}=\frac{1}{2}({{F}_{1}}-{{f}_{9}}) \\ {{f}_{4}}=\frac{1}{2}({{F}_{1}}-{{f}_{9}}) \\ -{{f}_{5}}-{{f}_{7}}=-{{f}_{1}}-{{f}_{3}} \\ \begin{align} & {{f}_{5}}-{{f}_{6}}=0 \\ & -{{f}_{6}}+{{f}_{8}}={{f}_{2}}-{{f}_{4}} \\ & {{f}_{7}}-{{f}_{8}}={{f}_{10}}-{{f}_{9}} \\ \end{align} \\\end{matrix} \right.$ （1）

在（1）式中, F₁、f₁、f₂、f₉、f₁₀均为已知, 可以直接计算f₃=f₄=1/2（F₁– f₉）, （1）式可以改写为：

$\left\{ \begin{matrix} {{f}_{5}}+{{f}_{7}}={{f}_{1}}+{{f}_{3}} \\ {{f}_{5}}-{{f}_{6}}=0 \\ {{f}_{6}}-{{f}_{8}}=-{{f}_{2}}+{{f}_{4}} \\ {{f}_{7}}-{{f}_{8}}={{f}_{10}}-{{f}_{9}} \\\end{matrix} \right.$ （2）

将它写成线性系统如下,

$A\overset{\to }{\mathop{f}}\, =\overset{\to }{\mathop{F}}\, $ （3）

其中

$A=\left( \begin{matrix} 1 & 0 & 1 & 0 \\ 1 & -1 & 0 & 0 \\ 0 & 1 & 0 & -1 \\ 0 & 0 & 1 & -1 \\\end{matrix} \right)$$\overset{\to }{\mathop{f}}\, =\left( \begin{matrix} {{f}_{5}} \\ {{f}_{6}} \\ {{f}_{7}} \\ {{f}_{8}} \\\end{matrix} \right)$$\overset{\to }{\mathop{F}}\, =\left( \begin{matrix} {{f}_{1}}+{{f}_{3}} \\ 0 \\ -{{f}_{2}}+{{f}_{4}} \\ {{f}_{10}}-{{f}_{9}} \\\end{matrix} \right)$

求解（3）式便可得到流线模拟所需要的全部基质网格和裂缝网格流量。对于与裂缝相交的井筒网格, 在油藏模拟器计算出井筒网格的流量之后, 流线方法会把井筒网格与其他网格一样处理。

假设有一天然裂缝性油藏, 注水井位于油藏西南角, 生产井位于油藏东北角。为了集中分析流体在天然裂缝中的流动机理, 避免裂缝和基质相交造成的复杂结构的影响, 假设天然裂缝互不相交, 走向角为0° 或者90° 。假设裂缝渗透率与基质渗透率相差很大, 裂缝与油藏边界不相交, 流体为不可压缩流体, 流动遵循达西定律。

油藏模型尺寸为152.4 m× 152.4 m× 6.1 m, 划分为20× 20× 1共400个网格, 所有网格具有相同尺寸7.62 m× 7.62 m× 6.10 m。以0.00018 m³/s速度恒定注水, 生产井的压力为6.89 MPa。模型基本参数见表1。

表1

表1

表1 模型基本参数 参数 符号 单位 数值 注水速度 Q m3/s 0.00018 油藏温度 T K 288.706 油藏渗透率 Ki m2 1.97× 10-14 油藏厚度 h m 6.1 油藏孔隙度 ϕ % 10 原油黏度 μ Pa∙ s 0.001125 原油密度 ρ kg/m3 784.5 裂缝高度 hf m 6.1

表1 模型基本参数

为了显示基于EDFM的流线模拟的准确性, 将流线轨迹和飞行时间值与含水饱和度结果进行了比对（见图2）。生产600 d后的含水饱和度结果（见图2a）与流线轨迹和飞行时间值（见图2b）具有很好的一致性。

图2

	Figure Option ViewDownloadNew Window
	图2 含水饱和度结果（a）与流线轨迹和飞行时间值（b）对比图

在重建几何关系和相邻关系的基础上, 利用流线算法计算每条流线的轨迹以及流线上流体流动时间。每条流线都是独立计算, 为将来的并行计算节省时间。计算时间的长短主要取决于网格数量, 随着额外裂缝网格的引入, 会增加计算时间。对于无裂缝的均质油藏, 油藏模拟器的计算时间为75 s, 流线模拟的计算时间为240 s; 对于裂缝性油藏, 油藏模拟器的计算时间为120 s, 流线模拟的计算时间为300 s。可以看出, 额外裂缝网格的增加会增加计算成本, 但是在同等精确度的情况下, 本文方法的计算效率较高。

2.4.1 无裂缝均质油藏

无裂缝均质油藏（0号模型）的计算结果显示（见图3）, 流线模拟可以将无裂缝均质油藏水驱流体的分布情况清楚地刻画出来, 因为油藏是均质的, 所以水驱前沿均匀推进。

图3

	Figure Option ViewDownloadNew Window
	图3 无裂缝均质油藏流线和飞行时间值计算结果

2.4.2 裂缝导流能力对流动轨迹和流动时间的影响

为了研究裂缝导流能力对流动轨迹和流动时间的影响, 设定了高、低2种裂缝导流能力（见图4）。在油藏左边存在1条长度为83.8 m的垂直裂缝, 裂缝导流能力分别设定为7.52× 10^-9 m²· m（1号模型）和7.52× 10^-11 m²· m（2号模型）, 生产1 000 d时的流线轨迹和飞行时间值的分布如图4所示。

图4

	Figure Option ViewDownloadNew Window
	图4 不同裂缝导流能力下的流线轨迹和飞行时间结果

当裂缝导流能力为7.52× 10^-9 m²· m时, 水驱前沿到达裂缝之前, 驱替流体均匀推进, 当水驱前沿到达裂缝时, 流速突然增加, 在生产1 000 d的时候, 经过裂缝流动的水驱前沿已经接近到达生产井, 而其他部分刚刚推进到了离注水井118 m的地方（见图4a）。图4a显示一方面裂缝可以作为流体的快速通道, 另一方面流体只能在裂缝附近流动而不能到达远离裂缝的地方, 裂缝又起到了“ 屏障” 的作用。与图3相比, 图4a中出现被油藏工程师称为“ 死区” 或者“ 阻塞区” 的空白区域, 该区域的存在会减少驱扫体积, 降低水驱采收率。为了避免过早见水, 可以采取注入凝胶和纳米颗粒等方法来延缓通过裂缝快速流动的水流。

在生产1 000 d的时候, 高裂缝导流能力的水驱前沿到达了距离裂缝上端38.1 m的位置（见图4a）, 而低裂缝导流能力的水驱前沿只到裂缝上端15.2 m的位置（见图4b）。随着裂缝导流能力的减小, 裂缝的“ 快速通道” 作用和“ 屏障” 作用也相应减弱, 有更多的水流可以穿过裂缝到达另一侧, 从而减小了“ 死区” 体积。

2.4.3 裂缝数目对流动轨迹和流动时间的影响

为了研究裂缝数目对流动轨迹和流动时间的影响, 设定了3条裂缝（3号模型）和12条裂缝（4号模型）2种裂缝数目模型（见图5）。其中3条裂缝模型的裂缝长度分别为83.8 m、76.2 m和68.6 m, 相应的流线轨迹和飞行时间分布显示, 一部分水流被第1条裂缝改变路径的同时, 仍有一部分水流穿过第1条裂缝, 随后被第2条裂缝改变轨迹。第2条裂缝的“ 屏障” 功能比第1条裂缝弱一些, 不过“ 快速通道” 的功能还是很明显。第3条裂缝的“ 屏障” 功能相比第2条裂缝弱一些（见图5a）。12条裂缝模型的流线轨迹和飞行时间分布显示, 随着裂缝离注水井的距离越来越远, “ 快速通道” 作用越来越弱（见图5b）。

图5

	Figure Option ViewDownloadNew Window
	图5 不同裂缝数目的流线轨迹和飞行时间结果

2.4.4 裂缝位置对流动轨迹和流动时间的影响

为了研究裂缝位置对流动轨迹和流动时间的影响, 设定了3种裂缝位置的模型（见图6）。第1种裂缝位置模型（5号模型）包含垂直和水平2条裂缝, 长度均为83.8 m, 流线和飞行时间结果显示, 水平裂缝与垂直裂缝具有同样的“ 快速通道” 作用和“ 屏障” 作用（见图6a）。第2种裂缝位置模型（6号模型）包含8条裂缝, 流线和飞行时间结果显示, 由于第1条垂直裂缝被分成了2条短垂直裂缝, 流线轨迹发生了明显的变化, 更多水流得以穿过裂缝。另外, 水平裂缝的存在提高了见水速度（见图6b）。第3种裂缝位置模型（7号模型）在第2种裂缝位置模型的基础上将上边第1条水平裂缝的位置下移9.1 m, 流线和飞行时间结果显示, 由于这条裂缝的位置下移, 导致裂缝的“ 屏障” 作用更加明显, 进一步减少了驱扫体积（见图6c）。关于死区的定量分析对压裂方案设计和井位优化有所帮助。

图6

	Figure Option ViewDownloadNew Window
	图6 不同裂缝位置的流线轨迹和飞行时间结果

设计水驱方案时主要考虑见水时间和驱扫体积这2个因素。见水时间越晚, 驱扫体积越大, 水驱效果越好。目前见水时间是通过分析产水量的突变时间来确定, 驱扫体积是通过研究压力或者含水饱和度分布来粗略确定, 这种分析方法准确性较差。利用本文提出的流线技术, 见水时间可以通过计算所有流线上水驱前沿到达生产井的最小时间来确定, 这个步骤通过比较所有流线上的飞行时间值就可以完成, 大大提高了准确性; 驱扫体积可以通过累加至少有1条流线通过的网格体积而近似得到。本文研究采用均匀网格, 所有网格的体积相同, 可以用驱扫网格数来代替驱扫体积, 驱扫网格数除以总网格数即可得到驱扫效率。由此可以定量分析裂缝对水驱效果的影响（见表2）。

表2

表2

表2 见水时间和驱扫体积统计表 模型编号 见水时间/d 驱扫网格数/个 驱扫效率/% 0 3 912 400 100.0 1 2 169 348 87.0 2 3 532 400 100.0 3 2 608 378 94.5 4 2 387 376 94.0 5 2 259 324 81.0 6 2 056 375 93.8 7 2 047 370 92.5

表2 见水时间和驱扫体积统计表

表2结果显示, 随着裂缝导流能力的增高以及裂缝位置复杂程度的增加, 见水时间加快、驱扫体积减少。裂缝数目对于流体流动轨迹和流动时间的影响不能一概而论, 1号模型确实比0号模型见水时间快、驱扫体积少, 但是5号模型比1号模型有更多的裂缝, 见水时间却慢了4%、驱扫体积减少了7%。3号模型比1号模型多了2条裂缝, 见水时间却慢了20%、驱扫体积增加了9%。在驱扫体积接近的情况下, 4号模型的见水时间比3号模型快了8%, 表明并不是垂直裂缝越多越好。在驱扫体积接近的情况下, 4号模型的见水时间比6号和7号模型慢了近17%, 表明两组裂缝会比单组裂缝具有更快的见水时间。7号模型和6号模型的唯一差别是1条水平裂缝的位置, 但是驱扫体积却差了1%, 表明裂缝的位置对流体流动轨迹和流动时间的影响较大。