非饱和花岗岩残积土水-气两相驱替过程数值模拟

2021-11-20蔡沛辰

水文地质工程地质 2021年6期

蔡沛辰，阙云，李显

（福州大学土木工程学院，福建福州 350108）

非饱和土广泛存在于自然界中[1]。与饱和土不同，非饱和土是一种三相土，固相（土体基质）、液相（水溶液）和气相（空气）同时存在。事实上，水流在非饱和土体中的入渗过程实质上是水在下渗过程中驱替空气的两相流问题。孔隙细观尺度下的两相流研究具有十分广泛的工程背景，石油开采、废弃物处理、地下水有机污染及泄水建筑中的掺气水流等都与细观多相渗流机理有着密切的联系[2−4]。传统的两相流研究多采用试验方法，如玻璃刻蚀、薄片平面驱替等二维方法[5]，岩心驱替物理试验、宏观驱替CT 扫描等三维方法[6−7]，但试验方法研究成本较高，操作较繁琐困难，且精度一般无法满足工程需要，更不利于开展重复性研究。因此，为克服试验方法的不足，细观尺度下的数值模拟方法得到了迅速发展和应用。

数值模拟方法利用计算机刻画细观孔隙结构，模拟水驱气的整个动态过程，具有重复、方便和可视化等优点。目前，细观尺度下多相流的数值模拟方法主要包括：格子Boltzmann、孔隙网络模型、相场等方法。如：毛欢[8]基于格子Boltzmann 方法，利用Shan-Chen 模型在孔隙尺度上对两相不混溶驱替过程进行了模拟；Qren 等[9]基于孔隙网络模型方法，通过构建大量砂岩孔隙网络模型，进行了多相渗流模拟，并通过模拟结果预测了多相流的相对渗透率；陈民锋等[10]建立油水两相三维孔隙网络模型，模拟了储层岩心的初次油驱和二次水驱过程；朱光谱等[11]基于相场方法，对孔隙尺度下表面活性剂两相流进行了模拟研究。上述相关研究已取得丰硕成果，但多集中于岩石多相流的研究，且研究对象与实际孔隙结构存在一定差异；应用较广泛的孔隙网络模型方法无法直观展现任一时刻各质点流速变化情况及无法呈现相界面移动状态。Level Set 方法最早由Osher 等[12]提出，现已用于两相流的数值模拟领域[13−15]，但受限于孔隙结构理想化或油气/水两相驱替，且未涉及土体细观孔隙内的水驱气两相模拟。鉴于此，本文选用花岗岩残积土作为研究对象，对其进行工业CT 扫描，建立二维细观孔隙模型，借助COMSOL Multiphyscis 有限元软件，采用Level Set 方法进行数值模拟，刻画土壤内部各孔隙区域水-气两相驱替的动态可视化过程，并分析驱替过程中的驱替速度、阻力及效率的规律。

1 CT 扫描及细观模型构建

1.1 CT 扫描及图像处理

试验原状土选自福州市某地山坡，现场取样如图1所示。为保证取样过程尽量不破坏原状土孔隙结构，首先将取样位置周围土壤挖出，保留核心土，之后将中间预留的土体上部削成与模型盒相当的大小，再将模型盒慢慢向下压入，依次重复削压，直至土体到达模型盒底后再向下超削约5 cm 深度，以便将土柱从根部切断，最终获取尺寸为15 cm×15 cm×40 cm 的土柱试样。

图1 现场取样Fig.1 Field sampling

为保证原状土体孔隙CT 扫描图像的真实性，对其进行如下处理：首先，将获取的土柱试样用保鲜防水材料封装，防止土体内部水分过分蒸发；其次，将封装后的试样放置于定制的木质模板箱中，防止搬运途中的扰动对其孔隙结构造成较大影响；最后，运至实验室后，拆除模板和封装材料，进行工业CT 扫描试验。

CT 扫描试验设备为C450KV 高能量工业CT，工作电压为450 kV，电流为63 mA，扫描最低分辨率为0.15 mm。对获取的图像，通过Image J 中的Threshold功能进行二值化处理，形成只包含黑白两色的图像，再利用中值滤波对其进行降噪，除去图像中孤立的噪点，最终得到一系列15 mm×15 mm 的二维扫描切片(图2）。考虑到计算时间的因素，从原始图像中选取孔隙连通性较好的一部分区域作为模拟对象，并将其等分为4 部分，局部放大如图3所示，模型a—d 尺寸均为3.0 mm×1.6 mm。

图2 二维扫描切片Fig.2 2D scanning slice

图3 模型a—d 二值化图像Fig.3 Model a—d binarized image

1.2 计算几何模型

基于上述步骤，类似处理其他图像以获得不同孔隙几何模型。最后，将这些图像以数组的形式存储在MATLAB 中，借助COMSOL-MATLAB 接口[16]将它们转换为计算域，并作为计算几何模型导入COMSOL中进行仿真模拟。以图3 中模型b 为例，构建的计算几何模型如图4（a）所示。此外，为清楚查看网格剖分质量分布情况，将图4（a）中红色圈出部位放大，得到图4（b），可以看出：模型b 的网格质量大多处于0.8 左右，质量分布较为均匀，也为后续数值模拟结果的精确度奠定了基础。

图4 模型b 的计算几何模型及网格质量分布图Fig.4 Computational geometric model and mesh quality distribution map of model b

2 水-气两相驱替仿真模拟

2.1 控制方程

2.1.1 纳维尔·斯托克斯方程

采用不可压缩的纳维尔·斯托克斯（N-S）方程描述流体流动[17]，控制方程如下：

式中：ρ—流体密度/（kg·m−3）；

u—流速/（m·s−1）；

p—压力/Pa；

I—单位矩阵；

μ—流体动力黏度/（Pa·s）；

F—体积力/（N·m−3）；

g—重力加速度/（m·s−2）。

2.1.2 界面张力

界面张力Fst由下式定义：

式中：σ—表面张力系数；

n—界面单位法向量；

δ—狄拉克函数；

Φ—水-气两相界面等值线。

2.1.3 Level Set 方法

“Level Set”接口通过跟踪水平集函数的等值线来确定流体界面，等值线Φ= 0.5 决定界面的位置。控制Φ的传递和重新初始化的方程为：

式中：γ—重新初始化参数/（m·s−1）；

ε—界面厚度控制参数/m。

由于水平集函数是一个平滑阶跃函数，因此可通过下式确定全局密度和动力黏度：

式中：ρw—水相密度/（kg·m−3）；

μw—水相动力黏度/（Pa·s）；

ρg—气相密度/（kg·m−3）；

μg—气相动力黏度/（Pa·s）。

2.2 材料属性

采用水作为驱替相，空气作为被驱替相，具体材料属性见表1，其中界面张力的设定见文献[18]。

表1 材料属性Table 1 Material properties

2.3 边界及初始条件

水-气两相驱替模拟中，依据构建的几何模型，进行边界条件设定：上下边界为进出口，左右边界为不透水层，壁为无滑移条件，初始界面设置为t= 0 时刻水相和气相的接触面。以模型b 为例(图5)，其他模型的设定与模型b 类似。其中，Pin为入口水压，Pout为出口水压。

图5 边界条件设定Fig.5 Boundary condition setting

3 结果与讨论

3.1 数值模型验证

为验证文中数值模型的正确性，采用体积守恒方法[19]验证几何模型中水-气两相驱替结果是否正确。由于多孔介质孔隙区域体积一定，故每一时刻气相与水相的总体积是恒定不变的。对于二维多孔介质孔隙模型而言，应采用孔隙区域面积来验证，即每一时刻气相与水相所占的孔隙面积是恒定不变的，模型a—d 孔隙总面积依次为8.56×10−7，1.83×10−6，1.75×10−6，8.87×10−7m2。

t= 0 时刻，气相存在于初始界面以下区域，水相存在于初始界面以上区域。水相从上边界进入孔隙区域，以初始界面为界限，通过水压力驱替气相，到下边界流出，并在孔隙区域内保持恒定的水压差。其中，考虑到运算效率及后期驱替可观察度等因素，以入口水压1.1 kPa、出口水压0.1 kPa 为例进行驱替运算仿真。图6 为模型a—d 水相、气相占孔隙区域面积随时间变化的曲线。从图中可以观察到：4 个模型中，水-气两相在孔隙区域中的面积之和恒定不变，且都等于各自模型的孔隙总面积，符合守恒定律，也验证了该数值模型的正确性。

图6 水、气相占孔隙区域的面积随时间变化曲线Fig.6 The area curve of water and gas in the pore area with time

3.2 水-气两相驱替

3.2.1 动态驱替过程分析

采用Level Set 方法模拟不同时刻水驱气的动态可视化过程，图7 为模型a—d 不同时刻下水-气两相驱替过程示意图，其中绿色为气相，蓝色为液相。

从图7（a）可以发现：流体会优先选择较大孔隙进行驱替，之后才会选择较窄孔隙，且在成圆度较高的孔隙处易出现“绕流”现象。计算知此处成圆度高达0.83，成圆度表征孔隙形状趋近于圆的程度，定义为：4π 倍的孔隙面积与孔隙周长平方的比值[20]。由图7（b）可看出：驱替流体存在“回流”现象，且部分死角孔隙流体不能进入，其原因为两种流体均不可压缩且不相溶，水相驱替使其形成了只含有空气的封闭空间。由图7（c）可知：笔直且较大的孔道，驱替速度相对较高，如通道1 明显快于2 和3。再结合图7（d）可见：在孔道相互交错区域，流体并不直接相融合，而是经过一段时间后才汇合。

综上，对于细观尺度水-气两相驱替模拟，Level Set 方法能很好地捕捉两种不混溶流体之间的界面位置。水-气两相驱替过程存在大孔隙优先流特征，且“绕流”现象一般易于出现在孔隙成圆度较高处。受水-气相参数限制，部分死角孔隙形成了只含气相的封闭空间。笔直且较宽的孔道，驱替速度相对较高，孔道相互交错区域，多条流体需经过一段时间后才能融合。

3.2.2 驱替速度可视化分析

运用后处理模块中的速度可视化功能，以模型b 为例，绘制不同时刻水-气两相驱替过程速度可视化图(图8)。可以看出：多孔介质模型众多孔隙中仅存在少数组成的主要渗流通道，且主要通道普遍有着孔隙半径大、孔道笔直的特征。由各不同时刻速度分布云图可知，流速较大处多出现在孔道相对较窄或相交汇处，且在驱替压差不变的条件下为5×10−4，1×10−3，1.5×10−3，2×10−3s 时，最大流速分别为0.878，0.984，1.03，1.05 m/s，即随驱替时间延长，流速有逐渐增大趋势。上述现象表明：水-气两相驱替速度受孔道迂回程度控制，且流体优先选择主要通道进行渗流，存在明显的“优势通道”。驱替速度与驱替时间以先快后缓的特征呈正相关变化，最大增长率为10.77%，最小仅为1.90%。

图9 为不同孔道横截面速度分布图，其中孔道a、b 位置如图7（b）所标注。从图中可以看出：孔道b 中速度最大位置出现在中心位置x=3.745 mm 处，而孔道a 中速度最大位置在4.815 mm 处，并不在中心位置4.775 mm，存在右偏现象，与前人研究结果有所不同[13]。分析原因是：孔道a 所处位置的上方流体流入的方向为右上到左下，相比孔道b，流体会偏向孔道右侧通过。此外，由图9（a）可知驱替持续时间越长，驱替速度越大，与图8 相对应，但孔道b 流速最大出现在5×10−4s 处，到1×10−3s时流速反而有所降低，不符合前文所得结论。原因分析如下：孔道b 在1×10−3s 时，出现了“回流”现象，如图7（b）所示，使得通过孔道b 的流量骤降，进而使其速度减小，降低幅度为21.62%。

图7 不同时刻水-气两相驱替过程动态示意图Fig.7 Dynamic schematic diagram of the water-gas two-phase displacement process at different times

图8 不同时刻水-气两相驱替过程速度可视化图Fig.8 Visualized diagram of water-gas two-phase displacement process speed at different times

图9 不同孔道横截面速度分布图Fig.9 Velocity distribution diagram of different tunnel cross sections

3.2.3 驱替阻力场分析

图10 描述了不同时刻水-气两相驱替阻力分布情况。可以发现：驱替阻力随驱替时间的延长而增加，且越靠近孔壁区域，红色程度越大，尤其是狭窄区域，表明驱替阻力最大出现在孔壁处，且孔道越窄，阻力越大。

图10 不同时刻水-气两相驱替阻力分布情况Fig.10 Resistance distribution of water-gas two-phase displacement at different times

3.2.4 驱替效率分析

水-气两相驱替效率Et定义为：t时刻下，孔隙区域内水相饱和度与水-气两相饱和度的比值。

式中：Swt—t时刻水相饱和度；

Sgt—t时刻气相饱和度。

为探究不同驱替压差对驱替过程的效率影响情况，以0.1 kPa 为初始驱替压差，分三阶段加压，每阶段加压0.3 kPa，进行驱替模拟，并得结果如图11（a）所示。从图中可以发现：0.1～0.4 kPa、0.4～0.7 kPa、0.7～1.0 kPa 驱替效率分别增加了25.49%、6.25%和1.47%，说明驱替效率随驱替压差的增大而增大，且第一阶段加压增速效果更加显著。

图11（b）为相同驱替压差1 kPa 下不同孔隙结构

图11 不同情况下驱替率曲线Fig.11 Displacement rate curve under different conditions

模型驱替效率的对比曲线。计算可知，模型a—d 孔隙率分别为26.59%、38.14%、36.43%和18.48%。由图可知：初期，孔隙率最小的模型d，驱替效率最为显著；之后，模型c 急剧上升，而模型a、b、d 增速明显减缓。结合图7 分析原因：模型c 虽孔隙率较大，但孔隙结构并不复杂，其孔隙均质性较好，无明显狭窄、不连通的孔隙。除此之外，根据文献[21]可知，驱替效率还受到孔隙结构形态、孔径、连通率等因素影响。

上述分析表明：驱替效率与驱替压差成正比关系，且初期加压增速效果显著，可达25.49%，而后期仅为1.47%。驱替初期，相比孔隙率大的模型而言，孔隙较小模型的驱替效率较大，而随驱替时间持续，驱替效率还将受孔隙结构特征的影响。