基于孔隙网络模型的非稳态两相渗流模拟研究

任广磊 吴 倩 李 闽 孙华超

（1.中国石化华北油气分公司勘探开发研究院，河南 郑州 450006；2.西南石油大学，四川 成都 610500）

0 引言

研究岩石多孔介质孔隙空间内的两相流动，在地质、地球物理、石油工程和化学工程等领域都具有重要的意义［1-3］。国内外许多学者采用连续介质理论对多孔介质多相渗流进行研究［4-5］，但是由于孔隙尺度毛细管压力、黏滞力的非连续性，目前对孔隙尺度油水两相渗透率机理缺乏足够的认识。岩石的单／两相流的渗流特性在很大程度上取决于岩石孔隙尺度的非均质性。利用Bernabé 等人［6-7］提出的逾渗网络模拟和归一化方法能够定量地描述孔隙尺度非均质性，即孔隙连通性和孔径分布对绝对液体渗透率、气体表观渗透率和电传输特性的影响［8-9］。然而，孔隙尺度非均质性对两相渗流的影响尚不明确。

非混相驱替模式通常可以分为毛管指进、黏性指进和稳定驱替3 种类型［10］。对于典型的多孔介质水／油驱替，油水界面的演化依赖于黏性力与毛细管压力之间的竞争［3］795。在这一过程中，如果某一油水界面处两端的压差大于毛细管压力，则油水界面可以向前推进；反之，某一油水界面处两端的压差小于毛细管压力，但是两端仍存在压差，则该处的油水界面不发生移动，进而不同部位的油水界面移动产生差异形成指进现象［3］796。国内外许多学者对动态驱替过程和动态网络模拟算法进行了大量的研究，通过实验和动态网络模拟划分了不同流体侵入过程中注入流体波及的区域，并分析了注入压力／速度、流体黏度等对驱替过程的影响。

传统的孔隙网络模拟中，根据基尔霍夫定律，流入和流出节点的流量之和为零。流体在孔隙网络中流动满足拉普拉斯方程，认为流体压力可以瞬间从入口传递到出口，即渗流力学中的稳态渗流模型。这一假设后来进一步扩展到动态网络模型中，形成了动态网络模拟技术［11］。但是实际的流体和岩石通常具有（微）可压缩性，同时流体压力无法瞬间传播到无穷远处，尤其是多孔介质尺度较大时，传统动态网络模拟技术无法描述压力传播过程，进而不能准确地描述多相流体渗流过程。因此，需要在传统孔隙网络模型中引入非稳态渗流理论。通过将稳态渗流的拉普拉斯方程转换为非稳态渗流的压力扩散方程，在分析流体黏度、毛细管数和多孔介质非均质性分布等因素对驱替过程的影响基础上，利用动态网络深入研究油润湿性非均质岩石中黏性指进的动力增长机制，以理清油湿多孔介质中的黏性指进现象，并提出抑制黏性指进的措施。

1 孔隙网络模拟

1.1 网络模型构建

孔隙网络模型作为研究流体在孔隙中流动规律的模拟手段，可以用于研究流体从上一级多孔介质渗流到更大一级多孔介质渗流中的流动状态［12］。笔者采用二维正方形网络模型构建孔隙网络。在二维正方形网络中，每个节点连接4个相邻节点，因此二维正方形网络的最大配位数zmax＝4。网络模型的连通概率把1～P倍数量的节点相互间连接的通道半径设为0，通过这种方法将一部分管束去除可以生成部分连通的不规则孔隙网络模型。例如：若网络配位数z为3，则是通过随机将二维正方形网络中25%的连接通道去除后生成的模型。二维正方形网络模型的渗流临界连通概率为Pc＝50%，对应的平均配位数为 2［6］5，即当z＝2 时，模型可能发生断开，流体无法在其中流动。为了探讨不同因素的影响，构建了节点数为22 500，油润湿角均设置为160°。通过随机分配网络中各孔道半径ri，模型的孔喉半径服从均匀分布［rmin，rmax］，其中，rmin＝180 μm，rmax＝280 μm。网络模型的平均孔喉长度

1.2 单相非稳态渗流模拟方法

在真实岩心中，由于流体的黏性作用，流体质点黏附在物体表面上，形成流体不滑移现象，即相对速度为零，因而产生摩擦阻力和能量耗散。因此假设孔隙网络中的流体流动遵循能量耗散最低的原则，且黏性流动仅在流体窜流分支方向发生，流动过程中孔隙网络模型遵循的质量守恒定律通过基尔霍夫定律描述，即流入的流体体积等于流出的流体体积，由此将真实岩心基质流动简化为孔隙网络模型流动［6，13］。

根据基尔霍夫定律，引入无流动边界条件，解得各节点的流动压力，由此求得各截面的平均流速。模拟过程中整体流动方向设定为水平方向。该模型中对于单独的节点i和j，设两点压力分别为pi和pj，两点间连接孔道的半径和长度分别为rij和lij，流体黏度为μ，水力传导系数为则两节点间的体积流量qij用泊肃叶原理可以表达为：

传统的孔隙网络模拟中，驱替相和被驱替相的压力默认相等，根据基尔霍夫定律，两节点间的总流量为零。则流体在孔隙网络中流动满足拉普拉斯方程：

通过泰勒展开、有限差分法以及质量守恒定律可得：

根据式（3）对孔隙网络中的所有节点进行方程构建，则可构建出孔隙网络的单相流线性方程组或矩阵方程：

实际渗流过程中，通常认为油水两相和岩石骨架具有微可压缩性。若考虑流体和岩石微可压缩性，可得到非稳态渗流方程：

通过泰勒展开以及Crank-Nicolson 隐式有限差分法，可将式（5）转换为矩阵方程：

可使用超松弛迭代法进行求解，迭代公式如下：

式中，α为松弛因子（经试算取1.75）。

在得到压力场后，即可求得整个孔隙网络的流量通量q。

1.3 两相流体界面移动模拟方法

通过调研，孔隙网络模型的多相流数值模拟方法有（准）静态网络模拟和动态网络模拟。其中，（准）静态网络模拟只考虑毛细管压力对渗流的影响，而忽略流体的黏性压降和时间等影响因素。动态网络模拟技术可研究孔隙尺度多相流的瞬态流动变化现象，可同时考虑毛细管压力、黏滞力、重力和时间等对多相流的影响。

孔隙网络中的多相流数值模拟存在以下假设条件：①所有的流体都被认为包含在孔隙管道和节点中，但是所有的压降都发生在节点之间的管道中；②孔隙网络中两相流体之间只存在一个界面；③网络中有两种不可混溶、不可压缩的流体流动；④两种流体在管口界面上的毛细管压差与管道半径成反比；⑤流体流量可以用泊肃叶方程计算；⑥孔隙空间中只有活塞式驱替（图1）［14］。

图1 活塞式驱替示意图

在黏度为μnw的非湿相侵入前，孔隙网络被黏度为μw的湿相流体占满。模拟的驱替过程开始后，侵入流体以一定的速率自孔隙网络左端注入，毛细管压力pcij使用杨－拉普拉斯方程求解：

单个管道的流量采用Washburn 管流方程［12，15］求解，当圆柱形的管道中存在两相流的弯液面时，qij就可以用Hagen-Poiseulle延伸方程表示：

在管道中只存在一个两相凹液面的情况下，μeff可以用下式计算［3］797。

当管道中只有单相流体时，pc＝0，此时式（8）可以简化为式（3）（μeff＝μnw或μw）。在两相流中，每个节点的总体积流量依然满足守恒定律，即∑jqij＝0 ，据此可以由式（6）构建线性方程组，通过逐次超松弛迭代法求解，同时需要根据时间步长Δt进行校正。驱替过程一直进行直到整个孔隙网络空间被侵入流体占满。然后重新计算并更新压力场和饱和度，进行下一步计算。整个过程保持注入速度Q不变。该研究中流体注入速度的取值范围介于5 × 10-9～1 ×10-6m3／s。流体黏度比M＝μw／μnw＝200（水的黏度为1 mPa·s，油的黏度为200 mPa·s）。

由于流动过程中两相界面随时间推进，因而必须选定一个时间步长，使该步长内每个两相界面均发生了适量的位移Δx，这里的“适量”是指必须在保证精度的前提下尽量减少运算次数，为此，引入最小时间步长和修正时间步长。修正时间步长是指将渗流瞬态划分为长度相等的时间间隔Δtk（k＝1，2，…），当Δtk足够小时，相邻时间步长的压力场变化可以认为是稳定、线性的，此时可以使用与单相流部分相同的方法，即超松弛迭代法对压力场进行求解。最小时间步长Δti（i＝1，2，…）是指在所有管道中的凹液面到达下一节点的时间中，选取Δti作为这一步运算总体的时间步长，此时除到达下一节点凹液面外的其他凹液面的位移为Δxij＝vij·Δtmin，由此可以求得该时间点的水力传导系数gij和两相在孔隙网络中的分布。显然最小时间步长法中Δti与压力降Δp和凹液面位置有关，因此每次迭代的步长取决于该次计算的具体情况，而不是都相等。Δti的引入使得此模型可以使用尽量少的迭代次数得出模拟结果，在保证计算精度的同时大大提高了计算效率。

该模拟技术与传统黑油模型模拟方法中的IMPES模拟方法有类似之处，通过隐式有限差分方法计算孔隙流体的压力场分布，然后显示计算两相流体界面移动。孔隙网络模型可以通过GPU 加速计算技术扩展至超过1×108个网络节点，实现基于同一渗流物理模型下的从孔隙网络、室内岩心和物理模型实验到井筒再到单井油藏模型尺度的多尺度跨越。笔者重点从孔隙尺度下非稳态两相渗流机理出发开展研究，通过模拟研究提出压制黏性指进的技术方法。

2 模拟结果及讨论

结合孔隙网络模型与模拟方法，首先对非稳态网络模拟算法进行了验证，在保证模拟结果准确性的基础上对黏性指进的影响因素进行了模拟，进而得到压制黏性指进的技术方法。

2.1 非稳态网络模拟验证

笔者建立了一种以均质网络模型验证动态网络模拟算法可靠性的方法。构建的孔隙网络圆盘形模型的直径为7.5 cm，该模型中心为流体的注入端，周围圆形径向方向为出口端（图2）。如图2所示，当注入低黏度流体驱替高黏度流体（不利黏度比）时，在均质多孔介质中，注入流体呈现较强的十字型黏性指进，这与前人的实验和模拟结果基本一致［16］。

图2 模拟结果验证图

传统孔隙网络模拟是基于稳态渗流理论建立的数值模型，在孔隙尺度稳态渗流理论的基础上，借鉴黑油模型考虑流体的可压缩性和压力传播，建立了孔隙尺度非稳态两相渗流理论和模拟方法。实际渗流过程中，流体的流动和孔隙压力传播同时进行，非稳态渗流理论能够较好地描述这一过程。模拟过程中压力传播过程如图3 所示，从图3 可以看出，随着流体从入口端注入，流体压力随着时间的增加从注入端逐渐向远端波及。这一模拟结果表明非稳态渗流模拟方法和计算程序能够正确模拟流体压力的传播过程。

图3 不同时刻压力传播场图

2.2 渗透率梯度压制黏性指进

从孔隙尺度到油藏尺度易发生黏性指进导致注入流体过早的突破［17］，基于此认识，笔者构建了具有一定渗透率梯度的圆盘形孔隙网络模型研究黏性指进与渗透率梯度的直接关系。定义渗透率梯度为：

构建的孔隙网络模型中，中间注入端口的孔喉半径最大，孔喉半径沿径向方向随着距离的增大而线性减小。笔者得到模拟结果（图4），根据模拟结果可以发现，渗透率梯度对黏性指进具有较好的抑制作用，具体应用条件为当注入速度较快时，需要较大的渗透率梯度才能起到抑制黏性指进的作用；当注入速度较慢时，较小的渗透率梯度也能够起到较好的黏性指进抑制作用。产生上述条件的原因可归结为压力传播、黏滞力和毛细管压力三者共同作用的结果：当没有渗透率梯度时，压力传播速度较快，因此注入速度越快，黏性指进越强；当有渗透率梯度时，由于渗透率梯度减小的方向上，压力传播速度逐渐减慢，从而使得注入流体能够在压力波及的范围内进行较为均匀地推进，从而能够提高注入流体的波及效率和采收率。模拟得到的渗透率梯度系数开始发生作用的取值界限为1×10-3。根据注入速度和渗透率梯度系数之间的关系，把图3整理为相图，得到图5。根据图5 中红色直线，拟合得到临界注入流速Qc与渗透率梯度系数λ之间有如下关系：

式中，a和b为拟合系数（a＝7× 10-10，b＝2.059）。

图4 不同注入速度和渗透率梯度系数下的水驱油过程剖面图

图5 黏性指进区与均衡驱替区划分剖面图

从图5可知，当注入速度小于临界注入流速Qc时（红色线左端）为均衡驱替区，当注入速度大于临界注入流速Qc时（红色线右端）为黏性指进区。从图4可知，均质地层本身不利于油气的采出，合理利用储层平面非均质性可以进一步提高油气田的采出程度。这一研究可为后续结合实际储层平面非均质性开展亿级大规模孔隙网络模拟研究，估算实际储层的渗透率梯度和临界注入流速Qc用于油气田开发实际应用奠定理论基础。

3 结论

1）考虑非稳态过程的动态网络模型能够精确描述孔隙级两相流动过程。

2）孔隙尺度渗流过程中，孔隙流体压力传播影响压力波及区域内流体的黏滞力和毛细管压力之间的平衡，进而影响两相流体界面移动。

3）在渗透率梯度减小的方向上，压力传播速度逐渐降低，注入流体能够在压力波及范围内进行较为均匀地推进，从而有助于提高流体波及效率。

4）黏性指进是决定驱替波及效率的主要因素，降低注入速度以及沿孔喉半径减小的方向驱替可以有效地压制黏性指进，从而有助于提高波及效率并最终提高采收率。

物理量注释：z、zmax分别为配位数、最大配位数，无量纲；P、Pc分别为连通概率、临界连通概率，%；ri、rmax、rmin分别为孔道半径、最大孔道半径、最小孔道半径，μm；为平均孔喉长度，μm；rij、lij分别为i节点与j节点间的孔喉半径和长度，μm；gij为i节点与j节点间的水力传导系数，μm3· mPa-1· s-1；qij为体积流量，m3／s；μ、μw、μo分别为流体黏度、水黏度、油黏度，mPa · s；pi、pj分别为i节点、j节点的压力，Pa；Ct为地层综合压缩系数，Pa-1；Q为注入速度，m3／s；Δt为时间步长，s；gij的下标ij表示与第i节点相连的第j个节点；q为孔隙网络的流量通量，m3／s；pcij为毛细管压力，MPa；γ为界面张力，N／m；θ为润湿接触角，（°）；μeff为两相流体的有效黏度，mPa·s；xij为与凹液面位置有关的无量纲数，无量纲；M为流体黏度比，无量纲；Δxij为凹液面从第i节点到第j节点的位移，μm；Δti为最小时间步长，s；Qc为临界注入流速，m3／s；λ为渗透率梯度系数，无量纲；li为模型空间任意一点到模型中心的径向距离，μm。