虚拟仿真APP开发在土体孔隙细观两相渗流研究中的应用

蔡沛辰， 阙 云， 李 显

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

1 研究背景

土体孔隙的渗流问题在坝体工程、边坡防护、基坑开挖等领域具有非常重要的研究意义[1-3]，由于土体渗流涉及水-气两相耦合问题，使得求解过程变得较为棘手[4]。而传统的渗流研究多集中于宏观尺度，用试验手段测得渗透系数来表征土体的渗流特性。但一方面这种研究无法直观准确地描述孔隙结构特征、流体渗流路径以及孔隙结构与渗流特性之间的内在联系[5]，另一方面，试验方法耗时、耗力、耗资，且不利于重复性开展。因此，为克服上述试验方法的不足，一些细观尺度下的数值仿真方法得到了广泛应用并迅速发展。

目前，细观尺度下数值仿真手段主要包括以下两大类：

(1)基于计算机编程的仿真手段。例如苗杰[6]借助Poreflow半代码编程模拟器，对煤岩进行了水-气两相流模拟，并将模拟结果与实测结果进行了对比验证；徐宗恒等[7]采用格子玻尔兹曼方法纯代码编程，对土体孔隙进行了整个渗流过程的定量化研究。

(2)基于专业有限元软件的模拟手段，例如李辉等[8]从土石坝渗流角度出发，借助COMSOL软件建立模型，对不同水位下的坝体浸润线、渗透水压力以及坝体内比降等相关渗流参数进行了研究；戚蓝等[9]、湛文涛等[10]和邢皓枫等[11]从降雨强度角度出发，分别采用ABAQUS和SEEP软件分析了降雨类型和降雨强度对边坡滑移特征时空分布规律的影响，并详细探讨了降雨强度和历时对路基土边坡渗流场的影响情况。上述研究虽然已取得了丰硕的成果，但计算机编程仿真方法要求研究者不仅需具备较强的专业知识，且需熟悉相应的编程语言，适用对象面相对较窄，而应用最为广泛的专业软件模拟方法虽对程序语言能力要求较低，但还是较侧重于针对专业的科研工作者使用，且重复性操作较为繁琐、耗时。

鉴于此，本文以原状花岗岩残积土为研究对象，对其进行工业CT扫描，建立反映真实土体孔隙结构的几何模型，采用虚拟仿真APP开发手段封装土体细观尺度下水-气两相渗流数值仿真平台，对孔隙内水-气两相流动过程进行动态可视化研究，并验证了结果的真实性，以期探寻一种高效、快捷的研究手段，使科研工作者面对友好的自定义界面，从而提高科研效率。

2 CT扫描及土体细观模型构建

试验原状土选自福建省福州市某地山坡，现场取样如图1所示。取样点选取植被茂密的位置，取样前先清理表面腐殖层，清理面积为100 cm×100 cm，厚度为20 cm，最终获取尺寸为15 cm×15 cm×40 cm的土柱试样。测试土样的基本物理参数，结果如表1所示，其中，渗透率定义为：在一定压差下，土体本身允许流体通过的能力，由于花岗岩残积土渗透性较小，渗透率测定一般采用变水头试验[12]。

表1 所选试样的基本物理参数

图1 试验原状土现场取样[5]

对获取的原状土采用英华公司C450KV高能量工业CT扫描仪进行XZ面扫描，CT扫描工作电压为450 kV，工作电流为63 mA，扫描最低分辨率为0.15 mm，之后得到一系列能真实体现土壤孔隙分布的CT扫描切片，如图2所示。

通过MATLAB软件中的imread和im2bw等函数功能，对CT扫描图像进行二值化处理，形成只包含黑白两色的图像，再基于小波变换对二值化图像进行降噪处理，除去图像中孤立的噪点，最后选取中间位置的切片，截取其中孔隙连通性较好的区域作为模拟对象，如图3所示，截取的区域大小为3 mm×1.65 mm。将图像以数组的形式存储在MATLAB中，借助COMSOL-MATLAB接口[13]将它们转换为计算域，并作为计算几何模型导入COMSOL中进行仿真模拟，如图4所示。

图2 试验土样二维扫描切片[5] 图3 本文所采用的土样二值化扫描图像 图4 进行仿真模拟的土样几何模型

3 虚拟仿真平台及APP开发

3.1 虚拟仿真平台

COMSOL Multiphysics被誉为“第一款真正的任意多物理场直接耦合分析软件”，是一款具备完善的理论基础、集功能性、灵活性和实用性于一体，并可通过内部预设的专业求解模块方便用户进行应用和进一步开发拓展的软件[14]。其中，COMSOL APP是通过高度专业化的用户界面与COMSOL模型进行交互的一种直观而有效的研究方式。

3.2 仿真APP构建过程

下面以渗流场研究中经典的水-气两相流问题为例，详细叙述如何通过基于COMSOL Multiphysics构建的虚拟仿真平台封装搭建独立运行的仿真APP，便于科研工作者省时、省力地研究土体孔隙中的渗流规律[15]。

3.2.1 数值模拟平台建立

(1)几何模型导入。将建立的计算模型通过COMSOL-MATLAB接口转换为计算域(图4)，并作为计算几何模型导入COMSOL中进行仿真模拟。

(2)分别采用N-S(Navier-Stokes)方程和Level Set方法描述流体流动状态和水-气两相界面的瞬时动态变化情况[16]。其中Level Set方法接口通过跟踪水平集函数的等值线Φ来确定流体界面，Φ=0表示流体为气相，Φ=1表示流体为水相，一般取Φ=0.5的等值线作为相界面。控制Φ的传递和重新初始化的方程为：

(1)

式中：Φ为水-气两相界面等值线函数；t为渗流时间，s；u为流速，m/s；γ为重新初始化参数，m/s；ε为界面厚度控制参数，m。

为使Level Set方程在计算过程中具有较高的稳定性，公式(1)中的ε通常可取值为ε=hc/2，其中hc为界面区域的网格大小。此外，可通过下式确定全局密度ρ(kg/m3)和动力黏度μ(Pa·s)：

ρ=ρw+(ρg-ρw)Φ

(2)

μ=μw+(μg-μw)Φ

(3)

式中：ρw、ρg分别为水相和气相的密度，kg/m3；μw、μg分别为水相和气相的动力黏度，Pa·s。

计算步骤如下：

第1步，将构建的多孔介质几何模型导入到COMSOL中，设置相应的尺寸参数，并用多边形矢量功能设置两相渗流的初始界面位置，将多孔介质模型分为两部分。

第2步，对多孔介质几何模型两部分分别添加研究者所需要的仿真材料属性，驱替相添加为“水”，被驱替相添加为“空气”。

第3步，采用水平集模块和层流模块进行多物理场耦合计算，驱替相与被驱替相本构关系均为牛顿流，且在动量方程中包含表面张力，边界条件为不透水层，边壁为无滑移条件。水相从上边界进入孔隙区域，以初始界面为界线，通过水压力驱替气相，到下边界流出，整个过程中孔隙区域内保持恒定的水压差。

第4步，使用软件自带的“自由三角形网格”功能对模型进行网格划分，类型和尺寸设置为流体动力学和极细化(若划分不成功，即由于多孔介质模型孔隙结构较为复杂或均质性差，需要对网格调整为粗化)。

第5步，使用相初始化求解器和瞬态求解器共同进行求解，求解器的时间步长为range(0～0.2 s,每10-5s计算1次)，点击计算，经过计算后可得到两相渗流的动态可视化过程以及速度场、压力场等示意图。

(3)交互搭建APP。使用COMSOL中APP开发器功能对模型建立过程进行交互APP封装搭建。自主设计封装完成的APP启动界面如图5所示。该系统主要包括“参数设置区”“边界设定区”“功能操作区”和“结果展示区”等，详细如图6所示。其中，“参数设置区”和“边界设定区”主要包含“模型导入”“材料属性设定”“边界条件选定”“压差参数输入”等；“功能操作区”主要包含“构建模型”“网格剖分”“计算”“创建报告”和“帮助”等按钮；“结果展示区”主要为模型可视化展示以及计算后得到两相渗流过程、速度场等结果示意图的展示。

图5 自主设计封装完成的APP启动界面

3.2.2 仿真APP封装

(1)主窗口设定。在主窗口下菜单栏中添加3个子项，并分别设置为“保存”“另存为”和“退出”，快捷方式链接为“CTRL+S”“CTRL+SHIFT+S”和“CTRL+E”。

(2)建立相应功能的表单，并在主表单中完成APP运行界面的整体布局规划。在表单窗口建立4个表单，其中第1个表单作为主表单用于APP窗口的整体搭建布局设计，为防止表单过多，将“参数设置区”“功能操作区”和“结果展示区”布置在其上；其他3个表单分别为集合表单、信息表单和APP主要功能简介表单，将信息表单和APP主要功能简介表单及主表单复合到几个表单中，以主表单为默认窗口展示。

(3)程序编写，添加相应功能。对于所需输入的参数及仿真过程在之前已完成，现可以直接链接到APP中，并为参数设置区、功能操作区和结果展示区添加相应功能。

“参数设置区”的功能实现选择全局定义命令下的所需参数进行链接，可通过添加Method并输入程序语句实现参数初始化功能；“边界设定区”需要插入“选择输入”，并以“显示”的方式与之前建模过程进行链接，显示于主表单的图像窗口；“功能操作区”的功能实现选择模型ROOT命令下的相关操作；“结果展示区”的功能实现选择结果命令下的相关可视化显示，且可在工具栏命令下添加图像放大、缩小、打印等基本功能。

(4)APP搭建封装完成。至此，所有按钮及相关功能已封装完成，通过compiler功能对其编译处理，并可根据自身需求设置相应的APP图标和启动界面，最后输出可供Windows、Linux和macOS三大系统使用的EXE格式的APP文件。

图6 虚拟仿真APP运行界面

4 结果与分析

基于CT扫描技术构建的土体孔隙细观模型，以水-气两相渗流为例(水相为驱替相，气相为初始相)，借助虚拟仿真APP开发技术，将复杂的理论和物理场封装起来，模拟不同时刻水-气两相渗流的动态可视化过程。为使得模拟结果更好地贴近于真实情况，土体孔隙渗流方向设为沿深度方向进行，并设置模型初始压差以Pin=110 Pa入渗、Pout=10 Pa流出，详细边界条件设置如图7所示，其他相关参数如表2所示。

图7 土体孔隙水-气两相渗流边界条件设定

表2 土体孔隙水-气两相渗流计算参数

4.1 仿真数值模型验证

为验证所构建数值模型的正确性，采用文献[18]的方法，验证几何模型中渗流是否符合Darcy定律。水流通过压力进行驱动，从上边界进入孔隙区域，到下边界流出，并在孔隙区域内保持恒定的水压差100 Pa。整个孔隙区域XZ面的水流流速云图如图8所示。由图8可以发现，在渗流路径上存在一些高速区，经计算得出整个孔隙区域的平均流速为0.041 4 m/s，平均体积流量为8.19×10-8m3/s。

图8 土体孔隙区域XZ面模型流速云图

通过改变模型边界的水压差计算相应的体积流量，对其进行参数化研究，图9为渗流体积流量随水压差的变化趋势。由图9可看出，渗流体积流量随施加的水压差基本呈线性变化，符合Darcy定律，也验证了该计算模型的正确性。

图9 土体孔隙水-气两相渗流体积流量随水压差的变化趋势

4.2 水-气两相渗流过程

图10为土体孔隙水-气两相渗流过程动态示意图。由图10(a)、10(b)可看出，水体会优先选择较大孔隙进行渗流，之后才会选择较窄孔隙。由图10(c)可发现，渗流相在成圆度较高的孔隙处易出现“绕流”现象，通过计算得出此处成圆度达0.76，成圆度表征孔隙形状趋近于圆的程度，定义为4π倍的孔隙面积与孔隙周长平方的比值[19]。再对比观察图10(c)、10(d)可以发现，渗流水体还存在“回流”现象，且部分死角孔隙水流不能进入，其原因是孔隙中的水、气均不可压缩且不相溶，使孔隙中形成了只含有空气的封闭空间。

图10 土体孔隙水-气两相渗流过程动态示意图

综上所述，通过虚拟仿真APP平台对于细观尺度水-气两相渗流进行模拟，能很好地显示两种不混溶流体之间的界面位置；水-气两相渗流过程存在大孔隙优先流特征，且“绕流”现象一般易于出现在孔隙成圆度较高处；受水、气相参数限制，部分死角孔隙形成了只含气相的封闭空间。

4.3 水-气两相渗流速度场

在水-气两相渗流细观速度场的基础上，运用COMSOL后处理模块中的粒子追踪方法绘制不同时刻流速场下的粒子追踪图像，如图11所示(粒子起点控制参数表示为X：range(0,0.006,3)，Z：1.65，单位：mm)，图中粒子轨迹线条的颜色与图例相对应，表示流速的大小；线条密度表示通过孔隙的水流流量大小。

对比图11(a)和 11(b)可知，渗流路径1-1′和2-2′相比其他路径孔道的弯曲程度更低，渗流速度也相应较大。再对比观察图11(c)和11(d)可以发现，t=4×10-3s时刻孔隙中的流量明显大于t=8×10-4s时刻，同时水-气两相渗流存在“绕流”现象，出现“绕流”现象部位的速度相比其周围速度更小，且这种现象随着渗流时间的增加而更加显著，如图11(d)中圆圈标出的“绕流”部位的平均“绕流”流速仅为0.005 8 m/s，而其周围的渗流平均流速为0.087 m/s，两者相差14倍。

图11 土体孔隙水-气两相渗流速度场不同时刻的粒子追踪图像

上述现象表明，水-气两相渗流流速大小受孔道弯曲程度控制，且流体优先选择连通性较好的孔道进行渗流。水-气两相渗流中发生“绕流”部位的流速远小于其周围的渗流流速。

4.4 水-气两相渗流压力场

图12为土体孔隙水-气两相渗流不同时刻压力场分布示意图。由图12可知，在初始阶段，孔隙所受的压力沿两相渗流方向逐渐减小，且孔径越小则压力越大。但随着渗流历时的持续增加，孔隙所受到的压力分布范围沿渗流方向逐渐扩大，且孔径大、连通性好的孔道的压力值越大。此外，由于土体内部的孔隙结构复杂，造成了压力场分布的不均匀性。

图12 土体孔隙水-气两相渗流不同时刻压力场分布示意图

4.5 讨 论

表3列出了本文与其他学者关于土(岩)体孔隙水-气(油)两相流的研究结果对比。由表3可以得出以下结论：

(1)目前国内外研究学者对两相流的研究多集中于岩石和理论模型领域，少有涉及真实原状土体的两相渗流研究。

(2)岩石类的水-气两相流研究一般均可得出优势流和指近现象，而除此之外，本文还发现土体孔隙中存在“回流”和“绕流”现象。

(3)本文与前人研究角度不同，不局限于直接研究土体两相渗流场的整体流速大小，而是研究了“绕流”区域内流速分布与周围区域流速的大小关系。此外，还研究了两相渗流动态过程中压力场的分布变化情况，发现土体孔隙两相渗流压力场分布存在不均匀性的特点。

表3 本文与其他学者关于土(岩)体孔隙水-气(油)两相流的研究结果对比

5 结 论

本文以原状花岗岩残积土水-气两相渗流研究为例，借助数值模型构建虚拟仿真APP，将复杂的理论和物理场进行封装，并验证了虚拟仿真APP开发在土体孔隙细观两相渗流研究中的可行性。得出以下结论：

(1)水-气两相渗流过程存在大孔隙优先流特征，且“绕流”现象一般易出现于孔隙成圆度较高的部位。受水、气相参数的限制，部分死角孔隙形成了只含气相的封闭空间。

(2)两相渗流流速主要受到孔道迂回程度的影响，孔径大、弯曲程度低的孔隙中渗流速度相对较大，存在明显的“优势通道”。发生“绕流”现象部位的流速远小于其周围的渗流流速。

(3)土体孔隙两相渗流压力场分布存在不均匀性。在渗流初期，孔隙所受压力沿渗流方向逐渐减小，且孔径越小则压力越大。但随着渗流历时的增加，孔隙压力沿渗流方向增大，且孔径大、连通性好的孔隙所受压力越大。

(4)土体孔隙两相流是一个复杂的动态过程，本研究中模拟了土体孔隙中水-气两相渗流过程，关注点为土体孔隙结构特征对两相流的影响，但未考虑水、气相互溶合性及基质域渗流的问题，故将此作为今后的研究方向。本文研究成果能够为原状土体细观渗流机理的定量化研究提供一定帮助，所采用的研究思路对其他类似岩土体的渗流研究也具有一定的借鉴意义。