缝洞型油藏储集空间以大型溶洞、溶孔及裂缝为主，裂缝为主要的导流通道，基块基本不具有渗透性^[1-2]。储集空间尺度差异大，从几微米到几十米不等，具有很强的非均质性^[2-3]。同时，油藏中存在多种缝洞组合模式，不同组合模式下的流体流动形式与剩余油分布不尽相同。由于这种复杂性的存在，许多学者提出了诸多不同的数学模型，从理论上揭示在不同尺度上的流动规律，但是难以用于大规模油藏的数值模拟^[4-8]。另有学者根据此类油藏微观缝洞结构特点进行物理仿真实验来认识不同缝洞组合模式下流动规律及剩余油分布，但仿真实验模型设计难度大，成本较高，费时费力，同时物理实验中的观察和测量也存在很多问题^[9-12]。

拟颗粒方法力求刻画微观尺度上单粒子的运动，将宏观上连续的流体采用拟颗粒性质的流体微团来表征，流体间各部分的力用拟颗粒之间的作用力来代替，发挥其在界面追踪上的天然优势，能够较好地模拟远离平衡态的系统^[13-16]。本文基于拟颗粒方法进行缝洞型油藏油水两相流动的模拟，力求刻画油藏中多尺度油水流动特征，从微观角度了解水驱油的真实细节，进而研究流体流动的宏观特征，判断不同缝洞组合模式下剩余油的分布状况。文中首先建立多相流动的拟颗粒数学模型；然后建立了其表面张力模型，通过在两相界面处引入一种斥力的方式实现表面张力的模拟；通过模拟两种典型缝洞组合模式的不混溶驱替过程，并与实验结果进行对比，从理论及现象上验证了该方法的正确性；最后基于该方法对缝洞型油藏的简化模型进行了水驱油模拟，并以此为基础研究了不同缝洞组合模式下剩余油的分布状况，为进一步的挖潜奠定了理论基础。

1 多相流动拟颗粒数学模型 1.1 拟颗粒方法基本原理

拟颗粒方法的核心思想是插值理论，对于任意一个宏观变量f(r)(密度、压力、温度等)的值都可以通过它周围一系列邻近点的权重插值得到，具体表达形式如式(1)所示

$ f(r) = \int {f(r')W(r-r', h){\rm{d}}r'} $

(1)

式中：W—权重函数；r，r'—某粒子到其附近粒子的距离，m；h—光滑长度，是权重函数作用范围的一种度量，m。

图 1为权重系数示意图，其中：m_k—粒子k的质量，kg；v_k—粒子k的速度，m/s；r_ij—粒子i和粒子j间的距离，m；κ—与r处权重函数类型相关的常数，它确定权重函数的有效作用范围。可以看出权重函数通常是偶函数，值随两个粒子之间距离的变大而变小。其在整个作用范围内满足归一化条件，必须满足Dirac的δ函数，并符合紧致性条件式(2)中的

$ W(r-r', h)\; = \;0, \;\;\;\;\;|r-r'| > \kappa h $

(2)

权重函数近似法具有关于光滑长度的二阶精度，权重函数有各种形式，包括四次方函数、高斯型函数%被认为、三次样条函数以及五次样条函数等。拟颗粒方法的稳定性强烈依赖权重函数的二阶导数，对于低雷诺数流动，当使用拟不可压缩状态方程时，三次样条函数会使压力场和速度场产生显著的噪音，五次样条函数就相对比较稳定，但计算量相对较大。采用式(3)所示的五次样条函数计算^[13]。

$ W(s) = \frac{7}{{478\pi }}\left\{ {\begin{array}{*{20}{l}} {{{(3-s)}^5}-6{{(2-s)}^5} + 15{{(1 - s)}^5}\;\;\;\;0 \le s < 1}\\ {{{(3 - s)}^5} - 6{{(2 - s)}^5}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;1 \le s < 2}\\ {{{(3 - s)}^5}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;2 \le s < 3}\\ {0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;s \ge 3} \end{array}} \right. $

(3)

式中：s=|r_i −r_j|/h；r_i，r_j—到粒子i，j的距离，m。

假设知道f(r)在一系列邻近离散点距离r₁，r₂，…，r_N处的值，则f(r)可以近似表示为

$ f(r) \approx \sum\limits_{j = 1}^N {f({r_j})} W(r-{r_j}, h){V_j} $

(4)

式中：V_j—粒子j所占的体积，m³。

由于V_j=m_j/ρ_j，式(4)进一步变为

$ f(r) = \sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}f({r_j})} W(r-{r_j}, h) $

(5)

式中：m_j—粒子j的质量，kg；ρ_j—粒子j的密度，kg/m³。

式(5)的误差主要取决于这N个无序位置的分布状态。一般认为，为了使这种近似具有较高的精度，二维情况下必须保证权重函数作用范围内有20个左右的粒子。类似的，可以得到变量一阶导数和二阶导数的表示形式，就可以基于周围粒子的信息得到相关变量，避免了计算时对网格的使用。

1.2 互不相容多相流拟颗粒模型

采用连续性密度法，将二阶导数形式应用到流体连续性方程，经变换可得到式(6)所示的连续性方程的另一种常见形式

$ \frac{{{\rm{d}}{\rho _i}}}{{{\rm{d}}t}} = {\rho _i}\sum\limits_j {\frac{{{m_j}}}{{{\rho _j}}}{v_{ij}} \cdot {{W'}_{ij}}} $

(6)

式中：ρ_i—粒子i的密度，kg/m³；t—时间，s；v_ij—粒子i相对粒子j的速度，m/s；W_ij^'—权重系数的导数，无因次。

式(6)不存在边界缺陷，完全适用于自由表面问题，同时适用于多相流体交界面处的密度计算，但不能严格保证质量守恒。因此，为了保证系统的质量守恒，在使用连续性密度法的同时，可以周期性地使用密度求和法。

N-S方程中压力项有许多不同的离散形式，为保证动量守恒，采用如式(7)所示的对称格式。

$ \frac{1}{{{\rho _i}}}{p'}_i = \sum\limits_j {{m_j}} (\frac{{{p_i}}}{{\rho _i^2}} + \frac{{{p_j}}}{{\rho _j^2}})W_{ij}^{'} $

(7)

式中：p_i^′ —粒子i对时间的压力导数，MPa/s；p_i，p_j—粒子i、j处的压力，MPa。

N-S方程中黏性项处理相对比较复杂，为了模拟激波问题，人工黏性的思想也被引入到拟颗粒中，后来在弹性力学的拟颗粒模拟中也得到了广泛的应用。拟颗粒中最常见的一种人工黏性为Monaghan型人工黏性，具体如式(8)所示

$ {\mathit{\Pi} _{ij}} = \left\{ {\begin{array}{*{20}{l}} {\frac{{-\alpha c{{\tilde \mu }_{ij}} + \beta \tilde \mu _{ij}^2}}{{{{\bar \rho }_{ij}}}}\;\;\;\;{v_{ij}} \cdot {r_{ij}} < {\rm{0}}}\\ {0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他\;} \end{array}} \right. $

(8)

其中

$ {\tilde \mu _{ij}} = \frac{{h{v_{ij}} \cdot {r_{ij}}}}{{r_{ij}^2 + 0.01{h^2}}} $

(9)

式中：v_ij—粒子i到j的速度矢量，m/s；α，β—无因次相关系数；c—流体的声速，m/s；ρ_ij—粒子i和j的平均密度，kg/m³。

Monaghan型人工黏性不仅提供了冲击波面必不可少的耗散，而且防止了粒子相互接近时的非物理穿透。

理论上，低雷诺数流动的模拟要求流体的真实黏性力。Morris把权重函数的一阶导数和有限差分结合起来，提出了如式(10)的表达式

$ \frac{1}{{{\rho _i}}}{\mu _i}^\prime {v_i} = \sum\limits_j {\frac{{{m_j}({\mu _i} + {\mu _j}){v_{ij}}}}{{{\rho _i}{\rho _j}}}\left( {\frac{1}{{{r_{ij}}}}\frac{{\partial {W_{ij}}}}{{\partial {r_i}}}} \right)} $

(10)

式中：µ_i^′—i粒子的黏度导数；v_i—i粒子的速度，m/s；µ_i，µ_j—i、j粒子的黏度，Pa·s。

这种表达式严格保证了动量的守恒，得到了广泛的应用。基于式(7)和式(10)，N-S方程的拟颗粒表达式为

$ \frac{{{\rm{d}}{v_i}}}{{{\rm{d}}t}} = - 2\sum\limits_j {{m_j}} \frac{k}{\rho }W_{ij}^{'} + 2\frac{\mu }{{{\rho ^2}}}\sum\limits_j {{m_j}{v_{ij}}} \frac{1}{{{r_{ij}}}}\frac{{\partial {W_{ij}}}}{{\partial {r_i}}} + {\rm{g}} $

(11)

式中：ρ—密度，kg/m³；µ—黏度，Pa·s；g—重力加速度，g=9.8 m/s²。

为了避免应用于不互溶液-液多相流时，界面处密度统计的失真导致可能出现人工表面张力问题，Tartakovsky和Meakin提出粒子数密度的概念。考虑到液-液界面的表面张力以及壁面处的润湿力等一些其他作用力，不互溶液-液多相流体系动量演化方程的最终通用表达式为

$ \frac{{{\rm{d}}{v_i}}}{{{\rm{d}}t}} =-2\sum\limits_j {{m_j}} \frac{k}{{{n_i}}}{W_{i{j^\prime }}} + \frac{{2\mu }}{{{n_i}^2}}\sum\limits_j {\frac{{{m_j}{v_{ij}}}}{{{r_{ij}}}}\frac{{\partial {W_{ij}}}}{{\partial {r_i}}} + {\rm{g}} + \frac{{{F_i}}}{{{m_i}}}} $

(12)

式中：n_i—i粒子的数密度，kg/m³；m_i—i粒子的质量，kg；F_i—i粒子所受其他作用力，N。

1.3 状态方程与积分算法

模拟不可压缩流体流动需要有效计算动量方程中的压力项，基于拟颗粒技术可以利用拟不可压缩状态方程来模拟不可压缩流动，常用的状态方程为

$ p = {c^2}\rho $

(13)

c必须足够大以保证流体的基本不可压缩，但又不能太大，以避免时间步长太小。Morris根据对N-S方程各项大小的估算，认为c²应该与三者中的最大值有可比性，即

$ {c^2} \sim \max \left( {\frac{{v_0^2}}{\delta },\frac{{\nu {v_0}}}{{{L_{\rm{0}}}\delta }},\frac{{{\rm{g}}{L_{\rm{0}}}}}{\delta }} \right) $

(14)

式中：v₀—初始速度，m/s；L₀—控制长度，m；ν—运动黏度，m²/s；$ \delta = \frac{{\Delta \rho }}{{{\rho _{\rm{0}}}}} $，密度变化率，无因次；∆ρ—密度变化量，kg/m³；${\rho _{\rm{0}}}$—初始密度，kg/m³。

直接把式(15)应用于多相流时，由于界面处密度的不连续，会导致相应压力的跳跃。因此，同样采用基于数密度的人工状态方程

$ p = {c^2}n $

(15)

式中：n—数密度，kg/m³。

拟颗粒的积分方式有多种，比如，欧拉算法、预估校正算法、蛙跳算法以及四阶龙格库塔算法等。欧拉算法相对简单，当前的位置、速度和作用力完全由上一时刻的各个物理量所决定，具体计算公式如下

$ {v_i}(t + \Delta t) = {v_i}(t) + \Delta t{f_i}(t)/{m_i} $

(16)

式中：v_i(t + ∆t)—t + ∆t时刻粒子i的速度，m/s；v_i(t)—t时刻粒子i的速度，m/s；∆t—时间变化量，s；f_i(t)—t时刻粒子i上的作用力，N。

$ {r_i}(t + \Delta t) = {r_i}(t) + \Delta t{v_i}(t + \Delta t) $

(17)

式中：r_i(t+∆t)—t+∆t时刻粒子i的距离，m；r_i(t)—t时刻粒子i的距离，m。

$ {f_i}(t + \Delta t) = {f_i}\left( {r(t + \Delta t), {v_i}(t + \Delta t)} \right) $

(18)

式中：f_i(t+∆t)—t+∆t时刻粒子i上的作用力，N。

为了保证计算的稳定性，显式积分必须满足以下几个时间判据。第一个判据为CFL条件，如式(19)所示

$ \Delta t \le 0.25\frac{h}{c} $

(19)

第二个判据为粒子的加速度f_a带来的额外限制，为

$ \Delta t \le 0.25\mathop {\min }\limits_a {\left( {\frac{h}{{{f_a}}}} \right)^{1/2}} $

(20)

最后一个是黏性耗散带来的时间限制

$ \Delta t \le 0.125\frac{{{h^2}}}{\nu } $

(21)

当精度足够高(h较小)或黏度很大时，式(21)通常是最主要的一项。实际上，权重函数的类型和粒子的排布方式都会影响式(19)~式(21)中的系数。

2 表面张力模型

流体的宏观性质是流体分子微观运动的统计平均结果，对分子之间相互作用的具体细节并不敏感。拟颗粒是基于对流体宏观行为的连续描绘，需要对分子层次上的成对粒子间的相互作用采用粗粒化的方法进行光滑。对于流体相的内部分子，分子之间的相互作用相互抵消，从而微观作用的净力为零。由于这种作用在相界面处不能相互抵消，从而导致了表面张力的产生。对于不互溶流体来说，两相在界面处的相互作用应该可以被理解为某种形式的斥力。因此，拟颗粒中不采用在所有粒子对之间引入不同形式的引力的做法，而通过在两相界面处的不同种类粒子间引入一种斥力来实现表面张力的模拟，大大简化了表面张力模型。

图 2为表面张力模型的示意图。可看出每一相流体的内部粒子并不受到额外的作用力，只有界面处的粒子才受到周围另一相粒子的作用力。而界面粒子的确定，宏观粒子方法无需其他数值方法所采用的复杂界面追踪技术，具有独特的优势。在实际处理不同粒子间的相互作用时，只进行邻近网格的搜索。当搜索到的一对粒子分别属于不同的流体时，这两个粒子都被确定为界面粒子。

计算过程中采用的表面张力表达式为

$ {\boldsymbol{F}_{ij}} = S\frac{{{\boldsymbol{r}_{ij}}}}{{r_{ij}^2}} $

(22)

式中：F_ij—粒子j作用在粒子i上的力，N；S—表征表面张力强度的参数，N·m，r_ij—i粒子到j粒子的距离矢量，m。

F_ij是一个斥力，它随着粒子之间距离的变大而变小，它仅仅在界面处的不同种类粒子间存在。表面张力系数是由两相性质决定的一个常数，它表现了两相的界面性质。从以上部分可以看出，模型中参数本质上对应于不同的表面张力系数。测定表面张力系数的方法很多，但基本上都是根据受力平衡原理，常用的有毛细上升法、最大泡压法、滴重法及躺滴法等。

3 算例分析 3.1 正确性验证

为了验证上述理论方法的正确性，设计双向连通型缝洞组合与上缝下洞组合两种模式，将拟颗粒计算结果与实际物理实验结果对比。

双向连通型缝洞组合模型结构示意图及几何尺寸如图 3所示。

在初始时刻，假设仅洞下面的裂缝充满水，其余部分全部充满油。初始化粒子排布时，首先在一个刚好覆盖整个计算区域的矩形里按三角形排布方式填满粒子。然后，根据缝洞的构型，把边界外但距离边界在光滑函数的作用范围内(3h)的粒子指定为壁面粒子。同时，把位于水和油分布区域的粒子分别指定为水粒子和油粒子。最后，删除所有其他剩余粒子。图 3b给出了图 3a中方框部分对应的初始粒子排布，其中的水粒子、油粒子及壁面粒子分别用不同的颜色表示。

在中性壁面条件下研究竖直裂缝与水平裂缝对驱替行为的影响。水和油的密度分别取ρ_w=1.0 × 10³ kg/m³和ρ_o=0.8 × 10³ kg/m³，黏度分别是µ_w=5× 10⁻⁴ Pa· s和µ_o=1.5× 10⁻² Pa· s。入口处控制区的长度取0.44 mm，来流速度为5 × 10⁻⁴ m/s，表面张力参数取3.12 N·m。整个驱替过程从时刻t=0开始。

对于竖直裂缝，水从底端注入系统，随着时间的延长，越来越多的水流入洞中，同时油也被持续不断地从上端的裂缝驱出。整个驱替过程中，油和水之间的界面非常清晰，且水基本以活塞流的方式不断地向上推进。最后，几乎洞里所有的油都被水驱出(图 4a)。图 4a左边为数值模拟结果，右边为物理实验结果，二者现象基本相似。物理实验中油水界面更加平坦，主要是由于不同的缝洞几何比引起的。模拟中的几何比仅为10，但实验中却达到1.67×10³。随着比率的增加，洞中流体的实际速度以级数的方式快速下降，使得界面变得更加平坦。

对于水平裂缝，水从左端注入系统。由于水的密度大于油的密度，重力对于这种状况下的驱替过程应该发挥支配作用。当水从左端流入洞后，它沿着壁面以小细流的方式快速沉降。在水流到达洞底部的过程中，它的速度逐渐降为零。随着水流的持续注入，越来越多的水在底部聚集，使得油水界面不断上升。然而，在界面到达裂缝高度后，流出左边裂缝的水几乎以水平的方式直接流入右边裂缝，从而高于裂缝的油不再被驱出(图 4b)。图 4b中，左边为拟颗粒计算结果，右边为物理实验结果，二者在现象上比较接近。与模拟结果相比，实验中沿壁面沉降的水流是如此的细小，以至于几乎看不到。实际上，这种差别可以从Froude数进行解释。Froude数计算公式为

$ {F_r} = \frac{v}{{\sqrt {{\rm{g}}l} }} $

(23)

式中：v—流速，m/s；l—特征长度，m。

随着孔洞直径的增大，F_r快速减小，重力的影响变得更加显著，从而实验中的水流更易于沿壁面沉降(图 5)。

对于如图(5)所示的上缝下洞组合模式，在中性壁面下进行拟颗粒计算，当竖直裂缝宽度达到1.0 mm时，水也不能沉降。只有当裂缝宽度进一步增大到1.5 mm时，水才开始沿着竖直裂缝的左壁面沉降，并且在裂缝的底端以液滴的形式不断掉落。液滴再穿过洞里的油，沉降到洞的底部。部分油被水压在了洞的底部。随着时间的延长，液滴持续不断地形成并继而沉降，最后在洞的底部合并。在这个合并过程中，最初的油水界面被强烈的扰动，一些油甚至喷溅出来，最终以小油滴的形式分散于水中。随着界面的不断上升，水变得越来越不容易沿着壁面沉降。水甚至从竖直裂缝的中部沉降。最后，当竖直裂缝中的水达到一定高度后，水直接穿过水平裂缝流出体系，导致竖直裂缝中的一大部分油最终被残留。

图 5b为裂缝宽度为50 µm时的实验结果。此时，水直接从水平裂缝流出体系，从而不能驱出洞中的油。而当裂缝增大到1 800 µm，注入水就能进入洞中把油全部驱出。该实验结果进一步说明了模拟结果的正确性。

3.2 缝洞油藏多尺度流动模拟

缝洞型油藏中孔、洞、缝内的流动状态差异很大，为典型的多尺度多相复杂流动，不同尺度上的流动特征也不相同，如细缝内的流动属于渗流，而大的空洞尺寸上的流动则属于管流或空腔流。针对这种复杂流动特性，采用拟颗粒模拟技术建立研究该体系的多尺度离散模拟方法，能够表征缝洞型油藏的多尺度结构，并考虑了体系中存在油水两相界面，以及两相与基质的相互作用。

通过对实际缝洞型油藏进行简化，设计一个包含缝、洞和缝的交叉点的算例，该算例对应物理尺度约3.88×1.20 m，初始时刻缝洞中充满油，水从左上角注入，油从右上角采出，驱替持续时间约6.93 h。其中油黏度为15 mPa·s，进口速度为5×10^-5 m/s。

图 6是模拟过程中几个时刻的图像，其中蓝色部分代表水，红色部分代表原油，可以看出水沿裂缝向前推进，在不同的缝洞组合模式下，注入水并不能能够进入所有的溶洞，缝洞的连接点不同造成了溶洞驱替程度的不同，在溶洞顶部形成大量的剩余油，有些溶洞受大裂缝的影响，不能被水波及到，从而形成了高导流通道圈闭的剩余油。

拟颗粒方法计算量巨大，基于该方法对真实的缝洞型油藏地质模型进行多尺度流动模拟，需要采用更高效的算法，性能更强的硬件系统，对油藏的生产动态进行准确的预测，为油藏的开发过程提供参考与指导。

4 结论

(1) 基于拟颗粒方法建立了互不相容多相流动数学模型，发挥其在界面追踪上的优势，通过在两相界面处的不同种类粒子间引入一种斥力的方式，实现了表面张力的模拟，并据此建立了表面张力模型，能够有效实现缝洞型油藏多尺度流动模拟。

(2) 通过将不同缝洞组合模式下的水驱油数值模拟结果与室内物理实验结果做对比，从理论分析与实验现象对比方面验证了该方法的正确性。同时指出，该方法能够取代室内物理实验，实现不同缝洞组合模式下的油水流动模拟，大大节省了物力财力，提高了效率。

(3) 该方法用于简化的缝洞型油藏，能够从微观上判断水驱油过程的真实细节，较好地模拟复杂油藏油水流动情况，判断不同缝洞组合类型下剩余油的分布状况，为实施剩余油挖潜措施奠定理论基础。

(4) 该方法计算量巨大，有待于进一步探索高效的算法，发展高性能计算机，提高模拟计算能力，从而进行大规模的油藏数值模拟。