基于数字岩心与格子Boltzmann方法的致密砂岩自发渗吸模拟研究

汪勇，孙业恒，梁栋，蔡建超

1 中国石化胜利油田分公司勘探开发研究院，东营 257015

2 中国石化胜利油田分公司科技管理部，东营 257015

3 中国石油大学油气资源与探测国家重点实验室，北京 102249

*通信作者, wangyong733.slyt@sinopec.com

0 引言

致密油气藏广泛发育于我国各大含油气盆地，其资源量约为常规油气的3倍[1]。中国致密油藏的勘探与开发起步较晚，但发展速度较快，目前已在鄂尔多斯、松辽、准噶尔、渤海湾等多个含油气盆地发现致密油勘探开发有利区[2-3]。随着常规石油资源量的减少，低渗—致密油藏的重要性日益突出，提高致密油藏开发技术水平[4]，对我国石油工业持续稳定发展具有重要战略意义。

在水力压裂和注水开发过程中，自发渗吸是低渗—致密油藏基质—裂缝间产生油水置换现象的重要机制[5-6]，基质中的油气得以进入高渗裂缝网络从而提高油气产能。致密油藏中自发渗吸是一个复杂过程，受多种因素影响：岩石微观结构、流体特性及边界条件等[7-8]。不同条件下，影响渗吸规律的主导因素不尽相同。精确刻画渗吸规律是提高油藏采收率的关键基础问题之一。

近年来，数值模拟方法已被广泛应用于油藏开发渗流特性的研究[9-10]，而基于数字岩心模型的孔隙尺度流动模拟为研究热点[11-12]。格子Boltzmann方法(Lattice Boltzmann Method，LBM)[13]是一种通过流体粒子的碰撞和迁移来描述流体流动的介观方法。相比于传统计算流体力学，该方法易于处理复杂几何边界和实现并行计算，常用于微观渗流模拟，或直接计算岩石宏观物性参数等[14-15]。

本文基于胜利油田东营凹陷金家—樊家鼻状构造带致密砂岩微观孔隙结构的数字岩心模型，应用基于LBM的多相流模型开展孔隙尺度自发渗吸模拟研究。研究聚焦于岩心微观孔隙结构与润湿性对自发渗吸过程中流体界面变化及采出程度的影响，有助于明确致密油藏水力压裂开发中自发渗吸效应的特征及探明其影响因素。

1 致密砂岩数字岩心模型构建

岩心样品来自胜利油田樊154区块，岩性为致密砂岩，埋藏深度2506.91～3205.01 m不等。由于原始样品具有一定的含油性，所以首先进行洗油处理。之后在洗完油的标准岩心上钻取了直径为3 mm的样品，用于CT扫描构建数字岩心模型。CT扫描的图像数据来自于VerXRM500型微米CT显微镜，光源电压30～160 kV，最高分辨率 0.5 μm。

共有三块岩心进行了微米CT扫描，扫描电压为110 kV，扫描分辨率为1 μm。扫描结束后，扫描出的二维图像由重构软件重构成三维数据体。图1a、b、c为此次研究区编号为1、2和3岩心的二维扫描原始灰度图像。从扫描图像可以看出，1号样品孔隙发育较好，孔隙连通性较高，孤立孔隙较少，孔隙数量较高；2号样品孔隙呈裂缝状分布，同样具有较好的孔隙连通性，孔隙数量较高。3号样品孔隙发育较差，孔隙连通性不高，孤立孔隙较多，孔隙数量较少。为了提升计算效率与排除孤立孔隙对渗吸模拟的影响，三块岩心的大小为1503像素的子体积模型进行一系列包括提取、滤波、二值化等处理手段从原始CT扫描数据中提取出用于自发渗吸模拟(图1d、e、f，孔隙度依次为0.186，0.139，0.103)。从图中可以看出，提取出的子体积模型很好的表征了三块岩心的孔隙结构，1号岩心孔隙较多，连通性较好；2号岩心孔隙呈片状发育，连通性较好；3号岩心连通孔隙较少，分布非均质性较强。

2 孔隙尺度自发渗吸数值模拟

此次研究中主要应用基于LBM的颜色梯度模型开展致密储层数字岩心模型自发渗吸模拟计算。下面简要介绍控制流体输运的LBM基础理论。

LBM利用离散空间网格格点上具有不同速度方向的流体粒子的碰撞和迁移来表征流体流动。一个格点某一离散速度方向的流体粒子占有量即为分布函数fi(x)，其在固定格子上的运动方程如下所示

式中，Ωi为碰撞算子，决定着演化方程表征的流动物理及刻画流动问题的能力[13]。流体流动过程中宏观物理量，如密度和速度可分别由下式计算得到

图1 编号为1、2、3号岩心的二维CT扫描灰度图像与三维子体积图像Fig. 1 2D CT scanning grayscale images and 3D sub-volume images of cores 1, 2 and 3

LBM的介观特征为从底层构建考虑多相流体微观相互作用的宏观流动模型提供了基础理论框架[16]。目前已发展了多种多相流模型，其中颜色梯度模型(Colour Gradient Model)由于其具有界面张力可独立调节的特性和较小的界面厚度等特点，被广泛用于多孔介质多相流体流动研究[17]。该模型以红蓝两色标记非混溶流体，通过在碰撞项中添加微扰动以引入表面张力，即在该模型中表面张力被视为压力局部各向异性。碰撞后非混溶流体粒子通过重新标色驱使流体流向相同颜色的流体区域，达到相分离的目的[16,18]。

由于致密油储层孔隙结构的复杂性及毛细流动过程中较小的毛管数(较低的流速)，为了提高基于三维数字岩心模型渗吸模拟的精度和稳定性，模拟中利用了优化扰动项的颜色梯度模型及多松弛(MRT)碰撞项[19]。优化的颜色梯度模型中包括三个时空演化方程。其中一个为全局分布函数，用于控制压力和速度场的演化，另外两个演化方程用于模拟两相界面的演化。优化的颜色梯度如下[20]：

两个独立的演化方程用于计算红蓝流体密度场(ρr和ρb)对流，其中红色流体的演化方程和平衡分布函数分别为

自发渗吸模拟开始前润湿流体(蓝色)位于入口端20层格子缓冲层内，非润湿流体(红色)饱和整个孔隙空间和出口端20层格子缓冲层。在整个模拟过程中两相流体不受外力作用，完全依靠毛管力驱动润湿流体渗入致密油储层孔隙空间，因此计算域六面均采用周期边界条件。与渗吸方向垂直的两个侧面，在迁移前需要与缓冲层内流体序参数进行匹配。固体壁面采用反弹边界以适应复杂的孔隙空间结构。由于本次研究在孔隙尺度考察孔隙结构特征和润湿性对自发渗吸的影响，并且忽略重力的影响，因此在模拟中润湿和非润湿流体的黏度均设置为0.01，密度设置为1.0(格子单位)。模拟中不同的润湿角是通过设定固体壁面序参数而实现。在不同孔隙结构自发渗吸模拟中设置的润湿角与不同岩心样品的实验测取值进行匹配(三块岩心实验测量的润湿角均在30°～40°范围内)。

3 模拟结果及分析

3.1 微观孔隙结构对自发渗吸的影响

为了直观对比不同微观孔隙结构自发渗吸驱替过程中的渗吸前缘演化特征，图2以一定时间间隔展示了自发渗吸过程中润湿和非润湿流体分布。结果显示不同孔隙结构内润湿和非润湿流体相间界面演化规律整体相似。这是由于在数值模拟中三种孔隙结构的润湿角设置值较为接近，且均为强润湿性。根据孔隙尺寸的相对大小，润湿流体的流动方式可明显分为两类：在大孔隙内润湿流体优先沿孔隙壁面或角隅流动，该现象造成大孔隙内两相界面明显滞后。当大孔隙内的润湿流体流经喉道处时，由于“卡断”效应，造成非润湿流体被捕获在孔隙中央。尺寸相对小的孔隙内润湿流体快速填充整个孔隙，当大孔隙周围有很多小孔隙存在时，大小孔隙内润湿流体填充速率的差异引起“绕流”效应，造成大量的非润湿流体残留在大孔隙中。1号岩心模型的孔隙结构尺寸大小和连通性相对均匀，非润湿流体主要因“卡断”效应而滞留；2号岩心模型孔隙结构片状发育明显及连通性相对较好，非润湿流体主要因“绕流”效应而滞留；3号岩心模型孔隙结构也呈现出片状发育，但由于其连通性差，该模型内非润湿流体主要因“卡断”效应而滞留。

图3为不同岩心采出程度随时间的变化。1号岩心孔隙结构的渗吸速率一直变化，计算终止时刻仍未达到稳定状态，渗吸过程仍在缓慢发生。2号和3号岩心孔隙结构渗吸开始时刻渗吸速率十分迅速，达到拐点后2号岩心后期渗吸速率几乎没有明显变化，而3号岩心渗吸速率减缓，渗吸过程仍在持续，直到达到下一个拐点，渗吸速率基本达到稳定。三种类型岩心孔隙结构自发渗吸的最终采出程度分别为0.60，0.70，0.51。2号岩心孔隙结构内的最终采出程度最大，这也可以从计算停止后两相流体的赋存状态定性判断出(图2)。而较大的最终采出程度与该岩心微观孔隙结构的片状发育及较好的连通性直接相关。

图4更为直观的展示了不同孔隙结构模型渗吸过程中渗吸速率随时间的变化规律。渗吸初始时刻，3号岩心孔隙结构的渗吸速率最大，1号岩心孔隙结构的渗吸速率最小。渗吸初始阶段，三种孔隙结构的自发渗吸速率均迅速减小，随后以一定的速率持续缓慢变化。在0.25×106～1.0×106计算时间内，2号岩心的渗吸速率仍然保持较大的速率，而在此时间段内其它两块岩心孔隙结构的自发渗吸速率相对较小，并以较小的波动缓慢减小。

3.2 岩石润湿性对自发渗吸的影响

为深入认识不同润湿性对致密砂岩自发渗吸过程中孔隙内流体动态填充和整体采收率的影响，以2号岩心的微观孔隙结构为基础，在不同润湿条件下开展自发渗吸模拟计算。模拟中不同润湿角设置值分别为20.0°，37.3°，60.0°，80.0°。

图2 编号为1、2、3号岩心子体积模型内自发渗吸过程中两相界面演化(蓝色为润湿流体，红色为非润湿流体，岩石骨架不显示)Fig. 2 Evolution of two-phase interface during spontaneous imbibition in the sub-volume model of cores 1, 2 and 3 (blue indicates the wetting fluid, red indicates the non-wetting fluid, and the rock skeleton is not shown)

图3 三块岩心自发渗吸采出程度随时间变化关系Fig. 3 The recovery degree of spontaneous imbibition of three cores varies with time

图4 三块岩心自发渗吸速率随时间变化关系Fig. 4 The spontaneous imbibition rates of three cores varies with time

图5以一定时间间隔展示了不同润湿条件下自发渗吸过程中的两相分布。结果显示，润湿性对渗吸前缘的形态及两相流体的空间分布有较大影响。强润湿条件下(润湿角≤37.3°)，润湿流体主要以“角流”形式优先侵入孔隙角隅。渗吸前缘整体杂乱、分散，主终端界面滞后明显。随着润湿角增大，以“角流”形式的流动明显减小，两相界面的形态也相应变得规则、紧凑。当润湿角进一步增大到80°时，润湿流体流动能力明显减弱，仅能侵入入口端附近的部分孔隙。

图5 2号样品不同润湿条件下孔隙内渗吸前缘演化(蓝色-润湿流体，红色-非润湿流体，岩石骨架不显示)Fig. 5 Evolution of imbibition interface in the pores under different wetting conditions of core 2 (blue indicates the wetting fluid, red indicates the non-wetting fluid, and the rock skeleton is not shown)

图6 渗吸初期不同润湿条件下非润湿流体分布(红色-非润湿流体，灰色-岩石骨架，润湿流体不显示)Fig. 6 Distribution of non-wetting fluid under different wetting conditions in the early stage of imbibition (red indicates the non-wetting fluid, gray indicates the rock skeleton, and the wetting fluid is not shown)

渗吸初期，强、弱润湿条件下(润湿角≤60°)模型入口端均有一定体积的非润湿相液滴析出(图6)，定性说明了在该润湿条件下可以发生逆向渗吸现象。但不同润湿条件下逆向渗吸发生程度和发生位置明显不同。润湿性越强(润湿角越小)，入口端面上尺寸较小和较大孔隙均可作为逆向渗吸发生位置。随着润湿角增大，逆向渗吸现象仅能通过尺寸较大的孔隙发生，当润湿性接近中性润湿时，基本不发生逆向渗吸现象。

不同润湿条件下渗吸采出程度随时间的变化如图7所示。可以看出，同一微观孔隙结构模型的自发渗吸驱替效果显著受润湿强度的影响。润湿性接近中性润湿时，自发渗吸现象不明显。当润湿角减小时，渗吸驱替效率明显提高，但润湿角由37.3°进一步减小到20°时，渗吸采收率曲线提升效果并不显著。渗吸采收率曲线的这种差异主要源于不同润湿条件下孔隙尺度渗吸前缘不同的演化特征。中性润湿条件下润湿流体侵入程度低并很快达到稳定，渗吸驱替效果最差。弱润湿(润湿角为60°)和强润湿(润湿角≤37.3°)条件下，渗吸前缘主要沿着片状孔隙和小孔隙形成的优势路径向多孔介质内部流动。由于同一模型内该优势通道较为相似，相应地这三种润湿条件下的最终采收程度较为接近。此外，由于润湿角越小，毛管力越大，渗吸速率越快，因此强润湿条件下的渗吸采收率曲线高于弱润湿条件下的采收率曲线。当微观孔隙内两相流体达到稳定时，强、弱润湿条件下模型从入口到出口范围内均有大量非润湿流体残留。其中，单一孔隙内的非润湿相液滴主要由于“卡断”效应被捕获，而位于模型中部较大体积的“团簇”状非润湿相主要由于“绕流”效应被捕获。计算终止时刻不同润湿条件下对应的最终采出程度如表1所示。润湿角越小，自发渗吸最终采出程度越高。

图7 2号岩心不同润湿条件下采出程度随时间的变化Fig. 7 The recovery degree of core 2 under different wetting conditions varies with time

图8 2号岩心不同润湿条件下自发渗吸速率随时间的变化Fig. 8 The spontaneous imbibition rates of core 2 under different wetting conditions varies with time

表1 2号岩心不同润湿条件下渗吸驱油效果Table 1 The oil displacement effect by imbibition of core 2 under different wetting conditions

图8直观展示了不同润湿条件下渗吸速率随时间的变化。渗吸初始阶段，润湿角越小，渗吸速率越大。强(润湿角≤37.3°)、弱(润湿角为 60°)和中性润湿(润湿角为80°)条件下的渗吸速率差异显著，而强润湿条件下(润湿角≤37.3°)的渗吸速率的大小和变化趋势较为接近。该变化趋势表明，在强润湿条件下进一步减小润湿角并不能显著提高渗吸速率。不同润湿条件下自发渗吸持续时间显著不同。强润湿条件下渗吸速率快，润湿流体快速突破并达到稳定状态；中性润湿条件下渗吸速率最慢，润湿流体侵入程度低，渗吸持续时间最短；而弱润湿条件下渗吸速率居中，渗吸持续时间最长。

4 结论

本研究通过致密砂岩的数字岩心模型进行了自发渗吸模拟计算，重点分析了孔隙结构特征和润湿条件对自发渗吸的影响规律。结论如下：

(1)致密砂岩储层的微观孔隙结构特征显著影响自发渗吸驱替机理和采收规律。形态特征以片状发育且连通性较好的孔隙结构中渗吸速率快且非润湿流体主要以“绕流”形式捕集，最终采出程度高；孔隙尺寸细小且连通性较好的孔隙结构内渗吸速率稳定，无较大波动，渗吸现象持续时间长，非润湿流体主要以“卡断”方式捕集，最终采出程度一般；形态特征以片状发育但连通性较差的孔隙结构渗吸速率波动显著，非润湿流体主要以“卡断”方式捕集，最终采出程度低。

(2)润湿性显著影响自发渗吸过程中渗吸前缘演化。润湿角越小，润湿流体以“角流形式”优先侵入孔隙角隅，两相界面杂乱、分散，主终端液面滞后明显，渗吸前缘后非润湿相滞留明显。

(3)润湿性对渗吸初期的逆向渗吸及之后的顺向渗吸过程影响显著。逆向渗吸的发生程度和发生位置与润湿性直接相关，强润湿条件下，逆向渗吸发生程度高和发生位置多。顺向渗吸过程中，润湿性越强，渗吸作用越明显，渗吸速率越快，渗吸采出程度越高。