基于水平集方法的原状土三维水气两相渗流特性数值研究

蔡沛辰，阙 云

(福州大学 土木工程学院，福州 350108)

1 研究背景

自然界中的土是一种三相体，不仅有固相(土体基质)和液相(水溶液)，而且还存在气相(空气)[1]。水在土体中入渗的过程实际上是水气两相驱替问题，气体的存在会对水流入渗产生一定影响。传统的两相流研究多采用宏观试验的方法[2-3]，但试验成本较高，操作繁琐，不利于开展重复性研究。

近年来，一些代表性体积单元(Representative Elementary Volume，REV)尺度下的两相数值方法得到了发展，其中多相流的数值模拟方法主要包括：流体体积法(Volume of Fluid，VOF)、相场法(Phase-Field Method，PFM)、格子玻尔兹曼法(Lattice Boltzmann Method，LBM)等。如：①基于VOF方法，罗书靖等[4]对宽尾墩三维紊流流场进行了水气两相流场模拟，并完成了对流场紊流特征的正确模拟；魏祥祥等[5]通过分析典型孔隙结构的微观渗流特征和剩余油受力情况，揭示了不同类型微观剩余油的动用机制和规律；②基于PFM，冯其红等[6]构建了砂岩油藏的三维孔隙结构模型，研究了驱替速度、流体性质、润湿性对剩余油分布以及采出程度的影响；③基于LBM，Zakirov等[7]研究了两相排水过程中多孔介质的无序性及其与流速、黏度比、表面张力的耦合对界面发展动力学的影响；张一夫等[8]研究了相间力对微孔通道内水-气两相渗流规律的影响。上述两相流研究已成果硕丰，但多关注于岩石或煤层等对象；PFM需要高精度的解析度来求解两相界面的精度，这将严重影响时间步长的计算[9]；LBM对计算机的运算内存需求巨大，尤其是复杂的真实三维孔隙模型运算时间漫长[10]。水平集(Level Set)法由Osher等[11]提出，之后众多学者已将其用于两相流仿真领域[12-14]，但由于量化技术、地域性差异等原因，多集中在二维孔隙结构模型层面，对REV尺度原状花岗岩残积土的真实三维孔隙两相渗流特性研究鲜有提及。

鉴于此，本文基于CT扫描与AVIZO-COMSOL交互对接技术，并借助Level Set法直观展现三维土体孔隙内水气相界面动态迁移过程、速度场及驱替阻力场分布情况，重点研究孔隙空间结构特征对两相渗流过程的影响。

2 试样与方法

2.1 取样及CT扫描试验

本次试验原状土取样选取福州一处植被茂密的坡地，现场取样如图1所示，最终获取4个15 cm×15 cm×40 cm(长×宽×高)的土柱试样，并试验测试土样的基本物理参数，见表1[15]。

图1 现场取样Fig.1 Sampling of undisturbed soil

表1 试验测试土样参数[15]Table 1 Test parameters[15]

参考课题组之前的研究[16]，经对比分析后选取孔隙连通性较好的试样3作为后续的研究对象，试样3的典型扫描切片如图2所示。

图2 典型CT扫描图像Fig.2 Typical CT scanning image

2.2 3D重构及REV选取

三维重构前，采用拟合法[16-17]得到所选试样的阈值为19 920，图3为不同深度Z处三维图像切面采用中值滤波降噪后的阈值分割(Z分别为9、18、27、36 cm)，其中蓝色为孔隙区域，灰色为基质区域。

图 3 阈值分割Fig.3 Threshold segmentation

后再通过AVIZO可视化软件自带的插值运算将一系列连续的CT扫描图像生成三维立体结构[18]，自此三维重构完成。

代表性体积单元(REV)，指能够有效表征土体宏观物理特性的最小模型单元体[19]，可有效减少计算机内存容量，加快计算速度。在众多学者[18,20-21]REV选取方法的基础上，本文在试样原始尺寸重构模型的基础上，以中心点位置向四周依次扩增尺寸，来获得对应的REV。分析REV尺寸与孔隙度之间的关系，结果如图4所示。从图4可以发现：试样单元尺寸>300体素时，REV孔隙度几乎不发生变化。综合孔隙度和运算效率的因素，最终选取300体素×300体素×300体素的REV作为后续原状土水气两相渗流特性研究的对象。

图4 REV选取示意图Fig.4 Schematic diagram of REV selection

原状土三维重构及连通孔隙模型构建可采用AVIZO软件实现[22]。对试样所选的REV借助AVIZO阈值分割功能，进行三维重构，最后对其进行孔隙连通测试，删除孤立、不连通孔隙，得到原状土的三维连通孔隙模型，其结果如图5(a)所示，实际尺寸为4.5 cm×4.5 cm×4.5 cm。

2.3 AVIZO-COMSOL交互对接

参考文献[23]中AVIZO与多物理场耦合软件COMSOL进行数据交互对接的方法，输出STL模型文件，见图5。至此COMSOL与AVIZO三维模型可视化数据对接过程完成，反之，需重新对网格进行划分。

图5 REV土体孔隙模型(300 体素)Fig.5 REV pore model (300 voxels)

2.4 计算理论

2.4.1 控制方程

采用不可压缩纳维-斯托克斯方程组描述流体流动[24]，控制方程为

式中：u为流体速度(m/s)；ρ为流体密度(kg/m3)；p为压力(Pa)；I为单位向量；μ为流体动力黏度(Pa·s)；∇为梯度算子；F为体积力(N/m3)；g为重力加速度(m/s2)。

界面张力Fst定义为：

Fst=∇·T,

(2)

T=σ[I+(-nnT)]δ,

(3)

(4)

δ=6|Φ(1-Φ)||∇Φ| 。

(5)

式中：σ为表面张力系数；n为界面单位法向量；δ为狄拉克函数；Φ为水气两相界面等值线。

2.4.2 水平集

水平集(Level Set)是一种用于界面追踪和形状建模的数值技术，主要用于描述运动界面的动态变化[11]，并通过跟踪水平集函数的等值线来确定流体界面(Φ=0.5)。控制Φ的传递和重新初始化方程为

(6)

式中γ(m/s)、ε(m)为重新初始化参数。为使Level Set方程计算时稳定性较高，通常可采用ε=hc/ 2，其中hc是界面区域的网格大小，γ值是模型中出现的最大速度。

由于水平集函数是一个平滑阶跃函数，因此可通过式(7)确定全局密度ρ和动力黏度μ。

(7)

式中：ρw、μw分别为水相的恒定密度和黏度；ρg、μg分别为气相的恒定密度和黏度。

2.5 计算参数及边界条件

采用水作为驱替相，空气作为被驱替相，具体材料属性设置如表2所示。

表2 模拟参数设定Table 2 Settings of simulation parameters

为使模拟结果贴近于现实，假设模型初始时刻孔隙中充满空气，沿土体深度方向进行模拟，并设定孔隙壁面为中性湿润。考虑到运算效率及后期驱替可观察度等因素，故压差设置较大，模型初始压力pin=1.1 kPa流入，出口压力pout=0.1 kPa流出。详细边界条件设置如图6所示。

图6 边界条件设置Fig.6 Setting of boundarycondition

3 结果与分析

为研究REV尺度下水气两相流特性，借助AVIZO-COMSOL交互对接技术对选取的REV孔隙模型进行两相渗流模拟研究。其中，选取的REV孔隙模型经测试，其连通孔隙度为0.484，孔隙网络模型(Pore Network Model，PNM)测试渗透率为0.416 μm2。

3.1 数值模型验证

为验证文中数值模型的正确性，采用体积守恒法[26]验证仿真模型中水气两相驱替结果是否正确。由于多孔介质孔隙区域体积一定，故每一时刻气相与水相的总体积是恒定不变的，孔隙模型的总体积为2.25×10-9m3。图7为水相、气相及总体积随时间的变化曲线。从图7可知：水气两相的体积之和随时间恒定不变，且都等于总体积，符合体积守恒定律，也验证了该数值模型的正确性。

图7 体积随时间的变化情况Fig.7 Change of volume over time

3.2 水气两相驱替

3.2.1 孔隙分布特征分析

图8列出了不同时刻渗流过程可视化分布，其中红色为气相，蓝色为液相。图8(a)、图8 (b)分别对REV模型孔道出口处较大孔隙区域进行了圈注，可知水气两相自进口端进入，优先会沿笔直、较大孔隙存在的孔道到达出口端，大孔隙优先流效应较为显著。从图8(c)可以看到，渗流结束后孔隙中仍存在气相，残余气主要呈块状分布于模型孔隙边、死角和孔喉突变部位。分析原因主要是水气相不可压缩特性致使水相无法侵入。

图8 水-气两相渗流过程可视化Fig.8 Visualization of the water-air two-phaseseepage process

3.2.2 速度可视化分析

图9列出了渗流稳定后水气两相渗流速度场的分布。如图9(a)所示，渗流稳定后，REV模型总体流速较小，计算知孔隙渗流平均流速仅为0.876 m/s；同时众多孔隙中存在少数主渗流孔道，普遍有孔道宽且较为笔直的特征，如图9(b) 、图9 (c)所示。此外，由图9中红色圈注部分可以看出，主渗流通道相比周围区域速度更大，计算得平均流速可达1.29 m/s。

图9 水气两相渗流速度场及流线分布示意图Fig.9 Velocity field and streamline distribution ofwater-air two-phase seepage

上述分析表明：水气两相渗流过程中速度受孔道迂回程度控制，且存在明显的“优势通道”，其中主渗流通道平均流速比模型整体高47.26%。

3.2.3 驱替阻力分析

为探究不同时刻水气两相驱替阻力分布情况，绘制如图10所示的阻力场分布。从图10可知：最大驱替阻力出现在孔隙壁面处，且孔道越窄，阻力越大，如图10中红色圈所示；随驱替时间的延长，驱替阻力最大值呈减小趋势(6.41 Pa·s/m2→3.29 Pa·s/m2→2.36 Pa·s/m2)，减小幅度较为明显，驱替阻力最小值却呈增大趋势(0.14 Pa·s/m2→0.25 Pa·s/m2→0.37 Pa·s/m2)，增大幅度较小。分析原因是：水气两相渗流过程中，土体孔隙中空气的存在，势必对流体驱替过程造成阻碍，但随着气体被逐渐驱替出孔隙，这种阻碍也逐渐减小，直至驱替过程结束，即土体基本趋于饱和状态。

图10 不同时刻水气两相驱替阻力场Fig.10 Water-air two-phase displacement resistancefield at different instances

4 结 论

本文基于CT扫描图像与AVIZO-COMSOL交互对接技术，采用Level Set方法研究了REV尺度的原状土3D水气两相渗流特性。主要结论如下：

(1) 渗流结束后孔隙中仍存在气相，残余气主要呈块状分布于模型孔隙边、死角和孔喉突变部位。REV模型孔隙中存在少数主渗流孔道，普遍有孔道宽且笔直的特征，孔隙优先流效应显著。

(2) 水气两相渗流过程中速度受孔道迂回程度控制，且存在明显的“优势通道”，其中主渗流通道平均流速比模型整体高47.26%。

(3) 最大驱替阻力出现在孔隙壁面处，孔道越窄，阻力越大。随驱替时间的延长，驱替阻力最大值呈减小趋势，而驱替阻力最小值呈增大趋势。

本研究中模拟了REV尺度的土体水气两相渗流动态过程，重点研究了中性湿润壁条件下孔隙结构特征对渗流过程的影响情况，仍未考虑到两者相互溶合及孔隙壁面湿润性(亲疏水壁面)的问题，故将此作为今后的研究方向。