基于微观渗流实验和仿真的致密砂岩气–水两相微观渗流机理研究

Integrated Experimental and Numerical Study on Gas-Water Two-Phase Flow Mechanisms in Tight Sandstone at Microscale

摘要: 致密砂岩微观渗流过程中，多孔介质内的气水两相驱替是重要的基础流动现象。在微观孔隙尺度下，由于孔隙尺度微小，毛细管力作用显著，导致非均匀运移，这使得传统方法难以准确追踪束缚水的动态行为。为此，需要寻找一种新型的微观孔道内气水两相运动轨迹记录技术，以实现对气水运移的效率和分布规律进行研究，从而为油气开采、CO₂地质封存等领域提供可靠技术支撑。本文针对气–水在孔隙尺度下的流体运移，利用高压微观渗流模型实验系统和微观气水两项渗流仿真开展孔隙尺度下超临界气水驱替的可视化机理研究。结合实验观测和仿真模拟，揭示气–水界面动态行为，为非常规油气开发和碳封存具有指导意义。

Abstract: Gas-water two-phase displacement in porous structures is a fundamental and critical flow behavior in the microscopic seepage of tight sandstone. At the microscale, the dominant capillary forces due to tiny pore sizes lead to heterogeneous flow, making it challenging for conventional methods to accurately track the dynamic behavior of bound water. Therefore, it is essential to develop a novel technique for recording the motion trajectories of gas-water two-phase flow in microchannels to investigate displacement efficiency and distribution patterns. This will provide reliable technical support for oil and gas extraction and CO₂ geological sequestration. In this study, the pore-scale fluid migration of gas-water under supercritical conditions is visually investigated using a high-pressure microscopic seepage experimental system and numerical simulations of gas-water two-phase seepage. By integrating experimental observations and simulation results, the dynamic behavior of gas-water interfaces is revealed, offering valuable insights for unconventional hydrocarbon development and carbon sequestration.

1. 引言

我国低渗砂岩气藏分布广泛、资源储量丰富且开发潜力巨大。自21世纪以来，低渗砂岩气藏的探明储量和产量呈现快速上升趋势，截至2020年，低渗砂岩气在全国天然气产量中的占比已达30% [1]。砂岩气藏储层普遍存在非均质性强、微观孔隙结构复杂、渗透率低等特征[2]，致使渗流规律复杂化、单井控制储量偏低，开发难度显著增加。其中，微观孔隙结构特征作为储层评价的关键参数，主要包括孔隙和喉道的几何形态、空间分布、连通性以及孔隙–骨架相互作用等微观特性[3]-[5]，这些特征直接决定了储层的储集性能、渗透能力以及流体赋存状态和运移效率。

针对致密砂岩储层非均质性强、孔隙结构复杂的特点，众学者从实验和模拟两种方法开展了对气–水两相渗流规律的研究和优选。微观可视化实验是一种基于仿真的物理模拟技术，其核心工艺是通过激光微加工或光刻技术，将真实岩心的孔隙网络结构高精度地复制在透明玻璃芯片上，从而构建可光学观测的微观孔隙模型[6]。孔令荣和曲志浩[7] 1991年开发的砂岩孔隙模型被用于水驱油实验，成功观测了不同流体(气体、煤油、正辛醇)的渗吸过程。2018年王璐等[8]采用气驱水微观可视化方法，探究了缝洞型碳酸盐岩储层中水相的生成机制与分布特征。2021年肖文联等[9]选取鄂尔多斯盆地的多块岩心进行了水驱油的微观可视化实验，分析了驱替类型与岩心渗透率和孔隙结构的关系，对剩余油的形态进行了整理。与传统的室内实验相比，计算机辅助数值模拟典型地克服了实验观测尺度的缺点，可以更准确地描述孔隙尺度的流体流动动力学。目前广泛应用的流体流动孔隙尺度建模方法主要有孔隙网络模型(PNM)、晶格玻尔兹曼方法(LBM)两大类[10]-[13]。孔隙网络模型在计算相对渗透率和残余油饱和度方面具有较好的适用性，但其采用的孔喉结构过于简化，难以准确反映两相渗流过程中相界面的动态演变及相体积分数的变化规律。虽然玻尔兹曼方法能够对岩石孔隙内的两相渗流过程进行定性模拟，但该方法同样存在局限性，无法精确表征孔隙空间内相界面的形态变化特征。水平集方法由Osher和Sethian提出，已广泛应用于气液两相流自由界面追踪研究。国内外许多学者采用该方法探究了两相渗流过程中的相态演变规律。例如，高亚军等(2016)结合Levelset方法，研究了孔隙尺度两相窜流的临界喉道直径比特征[14]。吴丰等(2019)基于玻璃刻蚀实验与数值模拟的对比分析，证实了微观气水渗流的动态特性[15]。近年来，随着计算技术的发展，有限元数值模拟方法被广泛应用于研究多孔介质内单相及多相流体的流动特性，并逐渐成为该领域的重要研究手段。

在前人的研究基础上，研究通过高温高压微观驱替实验与数值模拟相结合的方法，研究致密砂岩储层气水两相渗流特征。实验方面搭建了压力8 MPa，温度140℃的微观玻璃刻蚀模型实验系统，通过高精度CCD相机实时观测气–水两相在多孔结构中的动态分布特征；模拟方面基采用水平集法对气液两相压差下的相界面进行动态追踪，揭示气–水界面动态行为，为致密砂岩气藏开发提供指导意义。

2. 高温高压条件下微观渗流机理实验

多孔结构内超临界气–水两相流体的驱替是致密砂岩微观渗流中关键与基本的流动运移行为。孔隙尺度超临界气水驱替的流动特征、流体分布及其对运移能力的影响机制，是理解宏观两相流动的基础。本章拟采用自主研发的高压微观可视化渗流实验装置，重点研究超临界状态下气水两相流体的动态驱替过程，通过实时显微观测与图像分析技术，深入揭示孔隙尺度下气水两相运移的微观机理及其渗流规律。

2.1. 高压微观渗流实验系统与方法

取用真实岩心薄片YQ6，埋深2137.14 m，孔隙度为14.4%，渗透率4.96 mD。本研究采用刻蚀工艺加工技术制备的二维孔隙结构芯片(图1)，结合精密流体控制与显微观测系统，开展超临界气水两相流动实验。实验系统具有以下技术特征：(1) 可维持气相的临界压力条件；(2) 实现宽范围压力调控以研究降压过程中的相态变化；(3) 通过实时显微成像解析孔隙尺度下的流体分布与运移机制。

Figure 1. Real rock thin section and pore network-etched micromodel

图1. 真实岩心薄片和孔道刻蚀的微观模型

2.2. 实验装置与实验流程

本研究基于自主搭建的超临界气水驱排采综合实验平台(图2)，该平台由四个功能模块构成：微观模型观测子系统、流量控制子系统、温度控制子系统和气水反应器子系统四大部分组成。具体过程如下：

(1) 样品制备。实验体系在15 MPa和60℃条件下恒温搅拌24小时，确保气液两相达到充分饱和平衡状态。

(2) 模型处理。500 μL/min流速去离子水冲洗10分钟；120℃恒温干燥24小时。

(3) 驱替实验。采用去离子水饱和微观模型后，逐步加压至8 MPa。

(4) 观察分析。实验完成后，通过显微成像分析毛细滞留气液相的形貌特征及空间分布。

Figure 2. Micromodel displacement system

图2. 微观模型驱替系统

2.3. 实验结果

微观渗流实验在8 MPa，40℃下进行，实验采用CO₂作为驱替气体，此时气体为超临界状态，同时采用罗丹明B对地层水染色，呈玫瑰红(浓度为50 PPm)，在图3中，模型饱和地层水，饱和水过程以200 uL/min流量注入微观模型，同时建立实验压力、温度，超临界气驱水(排驱)过程都以100 uL/min流量注入。实验观测气相在入口区及多孔介质内部的局部驱替动态特征，一共进行3组实验。

Figure 3. Water saturation in the model

图3. 模型饱和水

(1) 不同压差下气驱水实验

开展不同压差条件下的气驱水实验，以探究气体在不同孔隙结构中的渗流规律。设置0.5 MPa、1.0 MPa、1.5 MPa三个压差梯度，通过微观渗流驱替实时监测气水两相分布及运移特征。重点分析驱替效率、残余水饱和度的变化规律，见图4。

Figure 4. Gas-water mutual displacement results under different injection pressure differentials

图4. 不同驱替压差下气水互驱结果

以0.5 MPa驱替压差，0.6 s气前缘突破，0.8 s开始，气向四周发散，3 s左右，达到稳定，剩余少量束缚水。超临界气突破时间晚，驱替速度慢，驱替效率低。以1.5 MPa的压差驱替，0.5 s气前缘突破，0.7s开始，气向四周发散，2.8 s左右，达到稳定，剩余少量束缚水。超临界气首先突破优势孔道，再由连接孔道向四周发散驱替结束后束缚水较少，驱替效率较高。

气体驱替过程表现出典型的非均匀推进特征，其动态行为主要呈现以下规律：1) 初始阶段优先突破高渗透通道；2) 次级阶段通过连通孔隙网络发生径向扩展；3) 驱替动力学参数间存在显著耦合关系，具体表现为驱替速率与突破时间呈反比，且低速率情况下整体驱替效率明显降低。部分孔道能清晰观察到残余水膜，部分残余水占据着孔道；驱替过后，大量地层水被捕获，在孔道内形成大量聚集，聚集形状主要分为：间断型、渐进段塞型、水膜型、段塞型，储层物性越好，残余水越少；通过实验结果并进行测量，统计了40处不同地方水锁情况；出现水膜现象处的喉道宽度在1.33~5.12 μm，平均喉道宽度在2.87 μm (图5(a))；出现非连通孔隙死体积剩余水现象处的喉道宽度在0.51~2.0 μm，平均喉道宽度在0.93 μm (图5(b))。

Figure 5. Throat width under different phenomena

图5. 不同现象下吼道宽度

3. 微观气水两相渗流仿真

微尺度流体流动数值模拟的核心在于正确选取物理场模型，需要依据流动特性筛选适用的控制方程作为求解约束条件。从介质类型角度划分，流体运动主要可分为多孔介质内流动与非受限流动两大类，其本质差异体现在流动约束条件的不同：前者受复杂孔隙结构限制，后者仅受宏观固体边界约束。在对流动体系进行类型判别时，需综合评估空间尺度、时间特性、计算域特征、边界条件、介质属性、流体密度、雷诺数、马赫数、温度场分布以及粘性耗散等多维参数指标，具体分类体系参见图6。

Figure 6. Classification of fluid flow [16]

图6. 流体流动分类[16]

本研究采用多物理场耦合方法模拟储层孔隙内的气液两相流动，流动控制基于Navier-Stokes方程，采用层流物理场接口，相界面追踪应用水平集方法实现动态捕捉，通过数值求解获得孔隙尺度下的速度场和压力场分布特征。具体方程表达式如下：

质量守恒：(连续性方程)

$\frac{\partial \rho }{\partial t}+\rho \nabla \left(v\right)=0$ (1)

动量守恒：

$\frac{\partial \left(\rho v\right)}{\partial t}+\rho \left(v\cdot \nabla \right)\left(v\right)=\nabla \left[-pl+\mu \left(\nabla v+{\left(\nabla v\right)}^T\right)\right]+F$ (2)

其中：ρ为密度，kg/m³；V为速度向量，m/s；μ为粘度，mPa·s；p为压力，MPa；F为体积力，N/m³。

本研究采用水平集方法追踪气液两相界面动态演化过程。该方法通过构建高维空间中的水平集函数，以其零等值面表征相界面位置。在数值求解过程中，通过对水平集方程进行迭代计算直至收敛，可准确描述气液两相界面的动态变化特征。相界面演化方程表述如下：

$\frac{\partial \phi }{\partial t}+v\nabla \phi =\gamma \nabla \cdot \left({\epsilon }_{ls}\nabla \phi -\phi \left(1-\phi \right)\frac{\nabla \phi }{|\nabla \phi |}\right)$ (3)

其中，方程左边表达式$\frac{\partial \phi }{\partial t}+v\nabla \phi$ 用于描述相界面的移动。式中，$\phi$ 为水平集变量，一般取值0~1，$\phi <0.5$ 为流体1，$\phi >0.5$ 为流体2。${\epsilon }_{ls}$ 为界面厚度控制参数，一般取$h_{\max }/2$ ，$h_{\max }$ 为最大网格尺寸，单位m。$\gamma$ 水平集函数重新初始化参数，单位m/s，一般取值流体流动的最大速度。

流体属性控制方程：

$\rho ={\rho }_1+\left({\rho }_2-{\rho }_1\right)\phi$ (4)

$\mu ={\mu }_1+\left({\mu }_2-{\mu }_1\right)\phi$ (5)

式中，${\rho }_1$ 为流体1的密度，kg/m³；${\rho }_2$ 为流体2的密度，kg/m³；${\mu }_1$ 为流体1的粘度，mPa·s；${\mu }_2$ 为流体2的粘度，mPa·s。

3.1. 几何模型的建立

本研究基于环境扫描电子显微镜(SEM)表征结果，系统构建了致密砂岩有机质孔隙网络的数值模型。具体建模方法如下：首先通过图像处理算法对典型有机质孔隙的SEM图像进行二值化处理和形态学分析，精确提取孔隙–喉道系统的结构参数；继而采用非结构化网格生成技术，在保持孔隙几何特征的前提下，对复杂孔隙空间进行自适应网格剖分；最后基于计算流体力学原理，建立考虑表面润湿性和毛细力效应的多相流控制方程。整个建模过程实现了从微观结构表征到宏观流动模拟的多尺度耦合，建模流程如图7所示。

Figure 7. Flowchart of physical-geometric model construction

图7. 物理几何模型构建流程图

3.2. 网格剖分

网格剖分作为数值建模的关键环节，直接影响计算结果的收敛性。本研究采用混合网格策略：1) 基础网格采用自由三角形单元进行全域自动剖分；2) 基于曲率自适应算法优化局部网格密度；3) 在边界区域设置多层渐变式网格。针对二维模型特性，通过控制最大/最小单元尺寸和边界层参数，实现了计算精度与效率的平衡。

图8展示了最终的网格剖分结果，表1统计了相应的网格质量参数。该几何体中网格顶点数为48,554个，三角形网络数为90,461个，边单元数为6679个，顶点单元数为956个。

Figure 8. Mesh generation results

图8. 网格剖分结果图

Table 1. Statistics of flow model mesh parameters

表1. 流动模型网格参数统计表

流动模型	网格顶点数	三角形网格数	边单元数	顶点单元数
两相流动	48,554	90,461	6679	956

3.3. 模型边界条件

几何模型建立后，需合理设置物理场参数和边界条件以约束渗流模拟过程。

3.3.1. 边界条件的设定

边界条件的设定主要分为流体力学计算(CFD)和水平集两个部分，而CFD部分包括单相流的Darcy定律以及两相流中的层流，见图9。

对于临界气驱替地层流动模型，做如下假设：

1) 孔隙空间饱和盐水，左侧注入临界态气体驱替；

2) 左侧入口边界：施加定压Dirichlet条件，约束法向流动(无回流)及切向速度；

3) 右侧出口边界：设定定压Dirichlet条件，允许气液混相流动及存在切向速度分量；

4) 模型上下两条边为对称边界。

3.3.2. 求解算法的设定

本研究采用瞬态求解算法模拟砂岩储层微观孔隙内的流体流动动力学过程。基于有限元方法，通过时间离散化处理，实现流动过程的动态仿真。求解过程中采用自适应时间步长策略，其中当前时间步的计算结果将作为下一时间步的初始条件。时间步长参数设置详见图10所示。

Figure 9. Boundary condition configuration

图9. 边界条件设置

Figure 10. Transient solver timestep settings

图10. 瞬态求解器时间步长设定

数值模拟中时间参数的优化配置对计算性能具有重要影响。本研究综合考虑模型几何特征和流动特性，设置了三个关键时间参数：初始时刻(t₀)、时间步长(Δt)和终止时刻(t_end)。通过合理配置这些参数，在保证计算精度的同时显著提升了求解效率。

3.2. 流体驱替仿真结果分析

根据获得的完全驱替后平衡状态的仿真结果，观察气体驱替后剩余的残余饱和盐水的体积分数变化。压差分别为10 KPa、30 KPa、50 KPa、70 KPa。根据所得数据绘制了不同压差条件下的驱替效率图和束缚水体积分数图。从压力场的分布可以看出，在不同驱替压差下的束缚水体积分布存在差异，具体结果如图11所示。

Figure 11. The distribution of irreducible water volume under varying pressure differentials

图11. 不同压差下束缚水体积分布

通过添加探针，记录了不同驱替压差下束缚水体积变化情况如下图12所示，随着驱替压差从10 Mpa增加到30 Mpa，束缚水体积减小明显，由最初的36.8%降低为35.1%，随后驱替压差增加减少较为缓慢，这是由于部分束缚水受毛细管力和表面张力的共同作用被锁定在微小孔隙和颗粒表面，形成不可动水相，即使继续增大驱替能量也难以将其有效去除。

Figure 12. Irreducible water saturation under different displacement pressure gradients

图12. 不同驱替压差束缚水体积分数

4. 结论

本文采用高温高压微观驱替实验和仿真模拟的方法对多孔结构内气–水两相的驱替流动运移行为进行研究。得出以下结论：

(1) 在8 MPa、140℃条件下，气体驱替过程首先突破优势孔道，再由连接孔道向四周发散，驱替速度越小，突破时间越晚，驱替速度、效率越低。

(2) 孔道内残余水分布呈现两种典型形态：喉道宽度1.33~5.12 μm (平均2.87 μm)区域倾向于形成水膜；喉道宽度0.51~2.0 μm (平均0.93 μm)区域更易在死端孔隙中捕集水分。

(3) 基于CFD方法构建的考虑毛细管力–粘滞力–惯性力耦合作用的水平集模型能有效模拟气–水界面演化过程。

(4) 驱替压差由10 MPa增至70 MPa时，束缚水饱和度呈现两阶段下降特征：初期压差由10 MPa增加到30 MPa，束缚水体积分数由36.8%下降至35.09%，降幅达到4.65%；随着后期压差继续增大到50 MPa和70 MPa，束缚水体积分数下降至34.65%和34.13%，后期降幅减缓至1.25和1.5%。这种非线性变化源于毛细管力与表面张力协同作用导致的微孔隙水相锁定效应。

此外，研究所采取二维模型无法完全反映真实岩心的三维孔隙连通性和流体捕集机理，尤其在表征孔隙喉道空间分布、迂曲度及多相流体竞争运移等复杂过程时可能存在偏差，未来研究需结合三维成像技术以更精确地模拟实际储层条件。同时，本研究的结论基于致密砂岩的岩样数据，其孔隙结构、矿物组成和润湿性可能与其他岩石类型存在显著差异。因此，推广结论至其他储层(如页岩或碳酸盐岩)时需谨慎，建议通过多类型岩样的对比实验进一步验证。

综上所述，过合理调控驱替参数可有效改善气水两相渗流通道，但对纳米级孔隙中的不可动水相需采取特殊解堵措施，本研究为致密砂岩气藏高效开发提供了关键理论支撑。

NOTES

^*通讯作者。