岩石气水两相渗流的玻璃刻蚀驱替实验与有限元数值模拟对比

吴 丰 ，姚 聪 ，丛林林 ，袁 龙 ，闻 竹 ，张凤生 ，习研平

（1.西南石油大学地球科学与技术学院，成都610500；2.油气藏地质及开发工程国家重点实验室·西南石油大学，成都610500；3.延长油田股份有限公司，陕西延安716000；4.中国石油集团测井有限公司测井应用研究院，西安710077；5.中国石油集团渤海钻探工程有限公司第一录井公司，天津300280）

0 引言

多孔介质中流体的微观渗流（尤其是两相渗流）对地下油气资源开发、河流堤坝防渗、环境保护、化学工程等领域均具有重要意义［1-5］。在油气资源开发领域，微观两相渗流对流体渗流机理、提高采收率等研究具有重要影响［6-8］。国内外学者针对岩石微观两相渗流已开展了大量的研究，所采用的研究手段分为物理实验观测方法［9-15］和数值模拟方法［16-21］两大类。

为了便于观察，目前的岩石微观两相渗流物理实验观测方法都是基于二维孔隙结构的，主要包括宏观驱替CT扫描实验法［9-12］、岩心切片微观驱替实验法［2］、玻璃刻蚀微观驱替实验法［3，22-32］等。宏观驱替CT扫描实验法是利用CT扫描技术对柱塞岩心的流体驱替过程进行观察，可以获得流体在柱塞岩心中的分布及驱替优势通道，但驱替过程中所采用的CT扫描技术分辨率往往较低，无法清晰地观察到孔隙尺度的微观驱替过程及渗流特征，难以满足油气开发精细研究的需要。随着岩石两相渗流机理研究的不断深入，有学者针对宏观驱替实验的不足，提出了岩心切片微观驱替实验法和玻璃刻蚀微观驱替实验法等微观驱替实验方法［2-3，22-32］。这2种微观实验方法采用的驱替和观测手段相似，主要区别在于驱替和观测的对象（物理模型）不同。岩心切片微观驱替实验法直接将岩心切成薄片再进行驱替和观测，能尽量保留岩石的孔隙结构和润湿性特征，且制作技术难度较小，成本低，但模型透明程度差，无法清晰地观测流体的微观渗流特征［2］。玻璃刻蚀微观驱替实验法则是基于岩心铸体薄片制作透明的微观孔隙网络模型再进行驱替和观测，具有较强的可视化能力，便于直接观测微观孔隙中流体的动态渗流特征，因此在水驱油的机理研究中运用非常广泛［31-32］。近年来，众多学者在微观驱替实验的理论和技术上取得了不少成果［2，27-29，31］，针对常规岩心驱替实验和玻璃刻蚀微观驱替实验，不仅实现了显微镜下直接观察两相流体渗流动态特征，还进行了一些两相流体渗流参数定量分析的有益尝试。

近些年，随着微观渗流理论的进步及计算机性能的提升，岩石微观渗流数值模拟方法也迅速发展。目前使用较多的岩石微观渗流数值模拟方法主要包括格子玻尔兹曼方法［16-17］和基于Navier-Stokes方程的数值模拟方法［34-36］。这2种方法均以岩石微观孔隙网络为基础，对流体的运动进行模拟。格子玻尔兹曼方法起源于20世纪80年代中期，是一种以格子气自动机为基础发展起来的流场模拟方法，该类方法采用离散介质模型来模拟流体的运动。基于Navier-Stokes方程的数值模拟方法则将流体视为不间断的整体，用微积分方法来计算流体运动的参数，具体的计算方法则包括有限差分法（FDM）［18-19］和有限元法（FEM）［34-37］。一般来说，有限差分法是以微分方程为基础，通过离散化处理求解区域后，用差分和差商近似代替微分和微分商［18］，因此求解微分方程和边界条件，可以归结为求解一个线性方程组以得到数值解。有限差分法具有直观、网格划分简单等优点［18］，但在某些情况下对于复杂的边界是不适用的，由于有限差分法的缺点，它几乎被网格划分更加灵活、边界适应性更强的有限元法所取代，因此有限元数值模拟在流体渗流模拟中也逐渐得到广泛的应用［18］。

现有的岩石微观两相渗流研究主要针对油水两相流体，相对而言，针对气水两相流体的岩石微观渗流研究还比较少，并且到底气水两相流体的岩石微观渗流研究更适合物理实验观测方法还是数值模拟方法尚无定论。笔者利用相同的孔隙结构模型，分别采用玻璃刻蚀微观驱替实验法和基于Navier-Stokes方程的有限元数值模拟方法对气水两相流体的微观渗流特征进行研究，并对比分析物理实验观测方法和数值模拟方法的优缺点，以期为气水两相流体的岩石微观渗流机理研究和气藏的精细高效开发提供参考。

1 玻璃刻蚀驱替实验方法

1.1 玻璃刻蚀模型制备

在油气田开发实验中所采用的微观模型种类很多，微观玻璃刻蚀模型是微观模型中常用的模型之一。本文基于岩心铸体薄片照片构建玻璃刻蚀模型，如图1所示：第一，将岩石铸体薄片［图1（a）］在ImageJ软件中转换为灰度图［图1（b）］；第二，将灰度图在Avizo软件中作分割处理，划分出孔隙和岩石颗粒［图1（c）］；第三，基于分割图在Coreldraw软件中绘制刻蚀模型图像［图1（d）］；第四，将刻蚀模型图像刻蚀到玻璃上，制作玻璃刻蚀模型［图1（e）］。微观玻璃刻蚀模型的制备工艺（尤其是光刻技术）直接决定了玻璃刻蚀模型的精度，进而影响渗流实验的结果。由于目前光刻技术主要应用于光电技术领域，在油气资源开发中应用较少，一些学者［38-40］在微观玻璃刻蚀模型制作理论与技术上做了深入研究，为微观渗流研究提供了理论指导。本文所采用的微观玻璃刻蚀模型制备工艺的具体操作步骤分3步：①清洗匀胶。打开真空泵，吸住清洗后的玻璃片；先在玻璃片上滴打底胶作为增粘剂；烘烤10 min后滴上光刻胶，匀胶30 s，得到一定厚度的光刻胶膜。此处采用底胶和刻蚀胶配合使用的办法，能有效防止刻痕浅和孔隙“V”形槽变形。②曝光显影。在烤胶板上进行前烘10 min，然后进行曝光；利用光刻机进行曝光，曝光时间为80～120 s；将曝光后的玻璃片，放入显影液中，显影时间为5 min，后烘时间为10 min。③刻蚀去胶。在烤胶板上后烘并封蜡，在酸蚀液中腐蚀5～15 min；溶解经光照后的光刻胶，获得刻蚀模型。

图1 微观玻璃刻蚀模型制作流程图Fig.1 Flow chart for fabrication of micro glass etching model

1.2 玻璃刻蚀驱替实验

目前，基于微观玻璃刻蚀模型的驱替实验主要针对油水两相流开展研究，油和水的颜色差异较大， 可以比较清晰地对实验过程进行观测［24，27，41-42］。对于气水两相流研究，气体和水的颜色相近，因此气水两相流体的界面不易观察。为了更好地观测驱替实验过程，采用亚甲基蓝染料对水进行染色，以区分水和气体（氮气）。准备好驱替模型和驱替流体后，即可开始驱替实验。驱替实验装置如图2所示，具体实验步骤如下：①将微观玻璃刻蚀模型安装固定好，给玻璃刻蚀模型施加合适的围压，并对模型抽真空，再从入口端注入水（亚甲基蓝染色水）直至模型完全饱和水。②在24℃的恒温条件下，在入口端启动微量恒速泵，将初始驱替速度设置为 0.05 mL/min，将气体（氮气）持续（20 min）注入玻璃刻蚀模型孔隙中。③当气水分布达到稳定状态后，逐渐加大入口端流量，分别以0.1 mL/min流量驱替 20 min，0.2 mL/min流量驱替10 min，0.5 mL/min流量驱替5 min等，直至出口基本无水流出，驱替结束。整个实验过程中利用微摄像机拍摄实验过程的图像。

1.3 驱替图像处理与分析

不同驱替时刻的驱替图像可以由图像处理系统得到，为了使两相流界面更易区分及分析渗流特征，需对图片作进一步处理（图3）：首先将实验得到的驱替图片导入ImageJ软件，转换成灰度图；再将灰度图导入Avizo软件，进行分割处理，在气体和水所占孔隙空间填充不同颜色加以区分。

图2 气水两相流玻璃刻蚀驱替实验装置示意图Fig. 2 Schematic diagram of glass etching displacement experiment device for gas-water two-phase flow

图3 驱替图像处理流程示意图Fig.3 Flow diagram of displacement image processing

2 有限元数值模拟方法

2.1 有限元模型构建

二维有限元模型构建（图4）与玻璃刻蚀模型一样，也是基于岩心铸体薄片照片［图4（a）］，先经过灰度图转化［图 4（b）］和分割处理［图 4（c）］这 2 个步骤。有限元模型制作的第三步，是利用AutoCAD软件，基于分割后的孔隙网络图按照真实比例构建孔隙结构模型［图4（d）］；第四步是将孔隙结构模型导入COMSOL Mutilphysics软件来定义求解域［图4（e）］，模型左端为入口，右端为出口，闭合曲线所围区域为岩石颗粒，其余连通区域为孔隙和喉道，连通区域设置为求解域；第五步是采用自由三角形对求解域进行网格划分［图4（f）］，小喉道网格划分较为密集，大孔隙中间网格划分相对稀疏，越靠近孔隙壁网格划分越密集。

2.2 单相流模拟方法

二维单相流体渗流数值模拟实验采用缓慢的蠕动流接口求解Stokes方程。由于模型尺寸较小，所以在数值模拟过程中认为孔隙内流体是恒温等密度流。假设流体不可压缩，且不考虑流体的重力影响。模型入口在左端，出口在右端。在恒定压力条件下，流体从模型入口流入，单相流体渗流Stokes方程和连续性方程可表示为

图4 有限元模型构建流程示意图Fig.4 Flow chart of finite element model construction

式中：P为压力，Pa；I为单位矩阵；μ为流体动态黏度，Pa·s；u 为流体速度向量，m/s；F 为体积力向量，N/m3，ρ为流体密度，kg/m3；T 为绝对温度（293.15 K）。

边界条件：模型左端各孔隙入口压力相等且恒定不变，右端各孔隙出口压力也相等且恒定不变；除出入口外的其他边界均无流体的流入和流出；孔隙壁为中性润湿，流体流速在孔隙壁方向的法向分量为零。

2.3 两相流模拟方法

气水两相流数值模拟实验采用Navier-Stokes方程并结合Level set方法构建微观两相流数学模型，并利用有限元方法求解方程。模拟实验采用层流两相流接口，在模型中，指定材料1（润湿相）为水，饱和孔隙和喉道，材料2（非润湿相）为气体，气体从模型左端入口注入，驱替孔隙内的水，水由出口流出。用水平集函数φ定义气水两相流各自的体积分数，为水时，φ=0，为气体时，φ=100%，并且利用重新初始化的水平集函数φ值追踪两相流界面，研究气水两相流界面动态特征。气水两相流界面的动态方程可表示为

式中：φ为水平集函数；t为两相流驱替时间，s；γ是初始化参数（默认设置为1 m/s）；ε为控制界面厚度的参数（设置为hmax/2，hmax为两相流界面流经区域的最大网格单元的大小）。

除了能定义流体界面，水平集函数φ也被用于描述流体的特性（密度和黏度跃变），即

式中：ρ为密度，kg/m3；ρwater和 ρair分别为水和气体的密度，kg/m3；μwater和 μair分别是水和气体的动态黏度，Pa·s。

流体的流速与声音的速度相比特别小，可以忽略不计，所以可以假设水和气体均为不可压缩的。Navier-Stokes方程可以用来描述不可压缩两相流体的质量以及动量的转输，并且考虑了界面张力，即气水两相流动的Navier-Stokes方程和连续性方程为

式中：g为重力，N/m；Fst为气水界面张力，N/m。

两相流界面上的表面张力可表示为

式中：σ为界面张力系数，N/m；n为垂直于界面的单位向量；δ为界面上的狄拉克函数，1/m。

垂直于界面的单位向量为

狄拉克函数δ可以根据如下光滑函数近似计算得到

3 实验与数值模拟结果

3.1 研究区储层特征

川西地区蓬莱镇组致密砂岩储层按照物性及孔隙结构类型可分为3类（图5）：Ⅰ类储层，孔隙度＞11%，渗透率＞0.8 mD，核磁共振呈明显双峰分布，可动峰幅度较高，压汞排驱压力低于0.3 MPa，孔喉类型为中孔-细喉道［图5（a）］。Ⅱ类储层，孔隙度为9%～11%，渗透率为0.1～0.8 mD，核磁共振呈明显双峰分布，可动峰幅度中等，压汞排驱压力为0.3～1.0 MPa，孔喉类型为中孔-大喉道［图 5（b）］。Ⅲ类储层，孔隙度为6%～9%，渗透率为0.03～0.10 mD，核磁共振双峰分布不明显，可动峰幅度较低，压汞排驱压力为1～2 MPa，孔喉类型为小孔-大喉道［图5（c）］。针对Ⅰ类储层、Ⅱ类储层和Ⅲ类储层各选取了一张典型铸体薄片照片，同时开展玻璃刻蚀驱替实验、单相流渗流模拟和气水两相流渗流模拟。

图5 川西地区蓬莱镇组致密砂岩储层孔隙结构特征Fig.5 Pore structure characteristics of tight sandstone reservoir of Penglaizhen Formation in western Sichuan

3.2 玻璃刻蚀驱替实验结果

微观可视化玻璃刻蚀驱替实验选用3张具有不同孔隙结构类型的铸体薄片制作玻璃刻蚀模型进行气驱水实验，研究气水两相流渗流特征及驱替后的气、水分布规律。玻璃刻蚀气驱水实验初始阶段，水饱和孔隙空间，3个模型均是左端为入口，右端为出口；将气体由左端入口注入模型进行驱替实验，随着实验的进行，被驱替的水由右端出口流出，气体与水的基本参数如表1所列。

表1气体和水的参数Table 1 Parameters of gas and water

在气驱水实验中，由于孔隙和喉道的半径存在微观非均质性，不同孔喉半径具有不同的毛管压力，气体驱替时所受毛管阻力也不同，驱替前缘速度不同，存在明显的微观指进现象。即将气体由左端入口注入模型，当驱替压力大于喉道的毛管阻力时，气体率先通过喉道；当驱替压力小于毛管阻力时，气体先停止向前移动一段时间，然后发生运移跳跃；驱替前缘沿着毛管阻力较小的连通路径向前推进并率先到达出口，形成优势通道。当气体到达模型出口后，注入的气体多半沿着优势通道流动。从图6可看出，Ⅰ类储层驱替前缘向前推进较为均匀，指进现象较弱，驱替效果较好［图6（a）］；Ⅱ类储层模型出口端见气时，含水饱和度为56.87%，最终含水饱和度为52.29%［图6（b）］；Ⅲ类储层模型出口端见气时，含水饱和度为57.80%，最终含水饱和度为46.25%，含水饱和度降低幅度较小，驱替效果差，指进现象较为明显［图6（c）］。实验结束后，剩余水主要分布在小喉道、小喉道所控制的大孔隙及孔隙盲端。

图6 玻璃刻蚀驱替实验过程中不同驱替状态的两相分布图Fig.6 Two-phase flow distribution of different displacement states in glass etching displacement experiment

3.3 有限元数值模拟结果

3.3.1 单相流模拟结果

构建的Ⅰ类、Ⅱ类、Ⅲ类储层模型的尺寸分别为11.70 mm×10.98 mm，11.76 mm×10.97 mm，12.21 mm×10.93 mm。单相流数值模拟时，流体在恒压条件下从模型左端入口被注入到孔隙空间，气体首先占据大喉道和大喉道控制的孔隙，接着气体逐渐填充小喉道和小喉道控制的孔隙，随着时间的推移，模拟结果逐渐接近于稳态。如图7（a）—（c）所示，模型中流体的压力从左（入口）向右（出口）逐渐降低，由于3类储层模型孔隙结构不同，毛管阻力不同，压力等值线分布也就不同。Ⅰ类储层模型中压力大于0.67 kPa的区域较大，Ⅱ类储层模型中压力为0.57～0.67 kPa的区域较大，Ⅲ类储层模型中压力为0.26～0.57 kPa的区域较大。Ⅰ类储层模型比Ⅱ类储层模型和Ⅲ类储层模型具有更大的孔隙半径，因此，它的压力沿着驱替的方向降低的最慢。如图 7（d）—（f）所示，在流速图中，蓝色表示低流速，红色表示高流速。通道中心处的流速比入口处的更大，越靠近孔隙壁流速越小。Ⅰ类储层模型的孔喉半径最大，所以毛管阻力更小，流动速度更大。

3.3.2 两相流模拟结果

为了便于实验结果对比分析，气水两相流数值模拟实验采用了和微观玻璃刻蚀气驱水实验相同的孔隙结构模型，数值模拟的模型出入口条件、气体和水的相关参数也都与玻璃刻蚀驱替实验保持一致，气体和水的相关参数如表1所列。

图7 3类储层模型单相流数值模拟结果Fig.7 Numerical simulation results of one-phase flow of three types of reservoir models

图8 Ⅰ类储层模型气水两相流数值模拟实验驱替过程Fig.8 Displacement process of gas-water two-phase flow numerical simulation experiment of typeⅠreservoir model

气驱水数值模拟实验过程中存在明显的微观指进现象，即气体沿着一条或多条阻力较小、由大喉道和连通孔隙组成的连续路径向前运移并率先到达出口（图8—图10）。模拟实验初始阶段，气体由左端注入，当驱替前缘到达模型左端时，模型中Ⅰ类、Ⅱ类、Ⅲ类储层的含水饱和度分别为34.69%［图8（c）］、56.87%［图 9（e）］和57.80%［图 10（e）］。Ⅰ类储层模型比Ⅱ类和Ⅲ类储层模型具有更大的孔喉半径，驱替时间更快，当两相流界面前缘到达模型出口后，驱替范围继续变大（图8）。由于Ⅱ类储层和Ⅲ类储层模型孔喉半径较小，毛管阻力较大，模型右端出口见气后，流动通道中的驱替阻力减小，被注入的气体基本沿着优势通道向前运移，波及范围不再扩大，驱水效率基本不再增加，指进现象更为明显（图 9，图 10）。

图9 Ⅱ类储层模型气水两相流数值模拟实验驱替过程Fig.9 Displacement process of gas-water two-phase flow numerical simulation experiment of typeⅡreservoir model

图10 Ⅲ类储层模型气水两相流数值模拟实验驱替过程Fig.10 Displacement process of gas-water two-phase flow numerical simulation experiment of typeⅢreservoir model

4 实验与数值模拟结果对比

为了对比玻璃刻蚀驱替实验结果和有限元数值模拟结果的异同点，对玻璃刻蚀实验和数值模拟实验的驱替时间进行归一化处理。玻璃刻蚀驱替实验和有限元数值模拟结果对比表明：

（1）归一化后的数值模拟结果与玻璃刻蚀驱替实验结果基本一致（图11）。相同含水饱和度条件下，数值模拟与玻璃刻蚀驱替实验的两相流体界面前缘基本相同［图 6（h），图 10（e）］，数值模拟和玻璃刻蚀驱替实验都存在明显的指进现象。当驱替前缘到达模型入口端时，Ⅰ类、Ⅱ类和Ⅲ类储层模型的含水饱和度分别为34.69%，56.87%和57.80%，这与图11中曲线的拐点基本一致。数值模拟和玻璃刻蚀驱替实验中，Ⅰ类储层模型比Ⅱ类和Ⅲ类储层模型的孔喉半径更大，驱替时间更快。当两相流界面前缘到达模型出口后，Ⅰ类储层模型驱水效率继续增加，驱替范围继续变大，Ⅱ类和Ⅲ类储层模型驱水效率基本不再增加，波及范围不再扩大。

（2）数值模拟和玻璃刻蚀驱替实验均为微观可视化实验。玻璃刻蚀驱替实验可直接观测整个连续的驱替过程，比如气体在通过狭窄喉道时的运移跳跃现象。数值模拟则能较好地观测不同时刻的动态驱替特征。在驱替过程的初始阶段，气体率先驱替分布在大喉道及与大喉道相连的孔隙内的水，在驱替结束后，残余水主要分布在小喉道、小喉道控制的大孔隙和孔隙盲端。

（3）微观玻璃刻蚀驱替实验方法能较好地观测气体在通过狭窄喉道时的运移跳跃现象，但玻璃刻蚀驱替实验受仪器和刻蚀工艺精度不够等不利因素的影响，容易导致驱替压力不稳定、刻蚀成型的孔隙尺寸与设计尺寸有一定误差等现象。整体而言，玻璃刻蚀模型制备的制作成本高，模型精度较低，孔喉容易变形、成功率太低，其工艺技术还需进一步提高。

（4）数值模拟不受仪器和刻蚀工艺精度不够等不利因素的影响，可完美地重现岩石的微观孔隙结构特征，但数值模拟只考虑了有限的边界条件，因此无法准确模拟一些特殊的微观渗流现象。整体而言，有限元数值模拟方法则具有操作简单、可重复实验、运算能力强、成本低的特点，为微观可视化渗流研究提供了一种新方法。

（5）目前，玻璃刻蚀驱替实验和有限元数值模拟各有优缺点，结合二者开展气水两相渗流研究的效果更佳。

图11 玻璃刻蚀驱替实验与有限元数值模拟的含水饱和度变化规律对比Fig.11 Comparison of water saturation variation between glass etching displacement experiment and finite element numerical simulation

5 结论

（1）微观玻璃刻蚀驱替实验方法与有限元数值模拟方法均可用于岩石孔隙中气水两相流体的渗流研究，在采用相同孔隙结构模型的前提下，这2种方法所得到的流体驱替前缘形态和渗流规律大致相似。

（2）气水两相流玻璃刻蚀驱替实验过程中存在明显的指进现象；驱替结束后，分布在大孔喉和与大孔喉相连的孔隙中的水被驱替，未被驱替的水主要分布在小孔喉、小孔喉控制的大孔隙和孔隙盲端。

（3）在单相渗流模拟实验达到稳态时，流动通道中心的流速高于入口端，孔隙半径越大，流动通道中心的速度越高，压力下降越慢。两相流数值模拟实验过程中也存在明显的指进现象，气体沿着由大喉道及大喉道控制的孔隙组成的连续通道向前驱替，达到驱替稳态后，残余水主要分布在小孔喉、与小孔喉连接的大孔隙及孔隙盲端。

（4）微观玻璃刻蚀驱替实验方法能较好地观测微观指进现象和气体在通过狭窄喉道时的运移跳跃现象，但玻璃刻蚀模型制备时，孔喉容易变形、成功率太低。有限元数值模拟方法则具有操作简单、可重复实验、运算能力强、成本低的特点，为微观可视化渗流研究提供了一种新的选择。