基于QSGS 法3D 重构土体渗流场的LBM 数值模拟

阙云 ，邱婷，蔡沛辰，马宏岩，谢秀栋，薛斌

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

在基坑、边坡、堤坝、隧道等工程项目中，渗流作用往往会引起土体或结构的变形失稳，严重威胁工程安全［1-6］，明晰土体内部渗流机理是预防和治理结构渗流破坏的关键一环.目前，国内外研究者多通过宏观渗流试验、数值模拟方法研究土体的渗流特性，但传统渗流试验方法不可避免地会对原有结构产生扰动，试验数据将存在一定测量误差，且反映的是宏观渗流表现，无法从根本上揭示多孔介质孔隙内部的渗流过程，以及孔隙参数对渗流特性的影响机理.因此，从细观尺度出发探究土体渗流特性尤为必要.

当前细观渗流研究主要集中在两方面：重构多孔介质模型和细观渗流数值模拟.1）多孔介质模型重构方法.其主要包括球体沉降法［7］、硬球Monte-Carlo法［8］、分数布朗运动法［9］、随机生长方法［10-11］、模拟退火法（Simulate Anneal Arithmetic，SAA）［12］和CT扫描重构技术.采用前三种方法构建多孔介质模型相对复杂［13］，故常采用后三种方法.SAA法重构模型与真实土体孔隙结构较为接近［14］，随机法具有高效便捷且节约计算机内存资源的优势，但上述方法存在生成孔隙结构的尺寸和位置难以控制的缺点.为解决此问题，Wang 等［11］提出了四参数随机生长法（Quartet Structure Generation Set，QSGS）.研究土体孔隙参数对渗透特性影响情况时采用QSGS 法比CT 扫描法更为方便、快捷，CT 扫描法更适合研究真实土体孔隙结构内部的渗流特性.2）细观渗流数值模拟方法.其主要包括基于连续介质模型的CFD方法［15］、基于分子动力学模型的格子Boltzmann 方法（Lattice Boltzmann Method，LBM）［16］.LBM 方法具有可准确求解流体流动Navier-Stokes 方程以及处理复杂几何边界的优势，广泛应用于多相流和多组分流动中，如多相流密度对比模拟、非混相和部分混相驱替过程模拟、两相渗流模拟等［17-18］.

由于QSGS 能与LBM 方法灵活对接实现土体细观渗流场仿真模拟，诸多研究者已采用QSGS-LBM法研究了重构土体微细观渗流问题，如周潇等［19］基于QSGS 法重构土体微观结构，通过LBM 方法建立饱和土体渗流模型，探讨了微观土体结构中水的渗流变化规律.申林方等［20］对传统四参数随机生长法进行改进，并采用LBM 方法研究了恒定流速下重构土体的细观渗流场.Jun 等［21］基于QSGS-LBM 方法分析了蒸压加气混凝土砌块的平面渗流场特性.蔡沛辰等［22］采用QSGS 法构建土体模型，采用LBM 方法研究重构土在不同条件下的细观渗流机理.李滔等［23］通过QSGS-LBM 方法分析了多孔介质渗透率与孔隙尺度各向异性、孔隙分布非均质性的关系.

虽然QSGS-LBM 方法在细观渗流仿真中取得了较为丰硕的成果，但仍存在以下不足：一方面，传统QSGS 法重构土体模型存在土颗粒团大小、形状较为均一的缺陷，并不符合真实土体细观结构中土颗粒团分布情况.另一方面，LBM渗流场仿真中模型尺寸大小与计算机内存容量及运行速度之间存在矛盾，而最优模型尺寸的选取是解决问题的关键.在细观渗流研究领域，关于孔隙参数、孔隙结构各向异性等因素对3D渗流场特性影响机理的研究相对欠缺.

鉴于此，本文采用传统与改进QSGS法等效重构各向同性/各向异性3D土体孔隙模型，综合孔隙连通性和运算时间因素选取最优仿真模型，并结合LBM进行渗流场数值模拟，直观展现各孔隙区域内渗流速度及流线分布情况，分析模型尺寸、孔隙结构、孔隙率、土体颗粒大小及模型各向异性等因素下重构土体渗流场的细观渗流特性.

1 土体细观结构表征及构建

1.1 细观结构表征

作为一种多孔介质材料，土体通常可划分为土体颗粒（固相）和孔隙区域，其空间分布函数如下［20］：

式中：G（x）为随机变量，反映孔隙分布情况，且G（x）的平均值为.

1.2 四参数随机生长法

四参数随机生长法［11，13］可通过控制分布概率pc、生长概率pd、概率密度pirs（在i方向上第r相在第s相上的生长概率）和孔隙率n四个参数重构土体细观孔隙结构.令固相（土体颗粒）作为生长相，孔隙作为非生长相，则3D多孔介质模型的重构流程如下：

1）在3D 空间中，令生长相以一定的分布概率pc随机布设，pc必须小于设定的孔隙率n.

2）在3D 空间中，令分布的生长相单元以一定的生长概率pd朝19个运动方向生长，如图1所示.

3）重复上述步骤，当生长相达到设定的孔隙率n时终止生长，3D多孔介质模型重构完毕.

控制分布概率pc=0.01（大颗粒），孔隙率n=0.36、0.50 和0.64，各向生长概率为一定值0.01（假设土体为各向同性），经过生长，得到不同孔隙率下60×60×60格点大小的3D模型，如图2所示.

图2 不同孔隙率的 QSGS 模型（以pc=0.01为例）Fig.2 QSGS model with different porosity（taking pc=0.01 as an example）

1.3 改进四参数随机生长法

为更好地研究孔隙分布各向异性对3D 多孔介质渗透能力的影响，采用两个依次递进的生长过程来构建3D 多孔介质模型，当各方向的生长速度不同时，即可获得各向异性3D模型，详细步骤如下：

步骤1：大团颗粒生成较大的固相，固相生长核概率为pt，生长概率为pd，其中d代表大团颗粒的19个生长方向，孔隙率为n1（pt

步骤2：小颗粒团生成较小的固相，固相生长核概率为ps，生长概率为pk，其中k代表小颗粒团的19个生长方向，且ps远小于pt.

两个步骤相互结合，直至达到设定的孔隙率，最终生成具有各向异性的3D重构土体细观模型.

各向异性的参数可通过原状土体CT 切片中孔隙的典型结构和几何形态特征确定，若同时调整各个方向的pd可生成与CT 扫描模型良好吻合的各向异性多孔介质模型.本文模型通过设定孔隙率n=0.597（真实土体孔隙率平均值），生长概率p1～7=0.01、p8～11=0.2、p12～18=0.01、p19=0.005，所生成的3D重构土体细观模型如图3所示，图中灰色区域为孔隙、黑色区域为基质.

图3 考虑各向异性的 3D 重构模型Fig.3 Three-dimension anisotropic soil reconstruction model

2 格子 Boltzmann 模型

2.1 格子 Boltzmann 方程

格子Boltzmann 方程是离散形式的Boltzmann-BGK 方程［24］，其可从细观层面仿真模拟多孔介质中的流体渗流过程.不含外力项时，粒子分布函数的演化方程可表示为：

式中：Fα（ω，t）为格点ω位置处在t时刻沿α方向的粒子分布函数；eα为离散速度；δt为离散时间；τ为弛豫时间；Feqα(ω，t)为局部平衡态粒子分布函数.

2.2 格子 Boltzmann 基本模型

格子Boltzmann 模型通常由三部分组成，即离散速度模型（格子）、平衡态分布函数和分布函数的演化方程［25］.本文选择的模型是3D 层面最为常用的D3Q19模型，如图4所示.

图4 D3Q19 模型Fig.4 D3Q19 model

D3Q19 模型一共有如图4 所示的19 个运动方向，离散速度eα满足式（3）.

式中：c=δx/δt，其中δx为网格步长，δt为时间步长，均取1.

D3Q19模型的权系数wα如下：

DdQm模型的局部平衡态粒子分布函数可表示为：

式中：cs为格子声速，取值为；ρ为密度；wα为D3Q19模型的权系数；u为宏观速度.

同时基于LBE 基本模型，采用Chapman-Enskog展开方法［26］推导所对应的Navier-Stokes 方程.密度ρ、宏观速度u、压力差p及运动黏度系数μ与弛豫时间τ（无量纲）之间关系如下：

2.3 真实单位与格子单位转换

针对实际问题的数值仿真，常用的编程思路有两种：程序代码中物理量直接采用实际物理单位；程序代码中物理参数都做无量纲化处理［24］.本文采用的是后者，即无量纲化的格子单位，参考Succi 著作［27］中的方法对LBM 单位换算部分进行阐述.模型中基本参数（格子单位）：长度L、密度ρ、时间t、压力差p和运动黏度系数μ.模型对应的参数（物理单位）：长度为L′、密度为ρ′、时间为t′、压力差为p′、运动黏度系数为μ′.为实现上述两者的转换，需引入部分参考量，即参考长度Lr、参考密度ρr和参考速度ur［24］.定义为：

式中：c′和c分别为物理单位和格子单位下的声速.

对于一个特定的问题，模拟中的长度L、密度ρ、声速c和运动黏度系数μ是已知的；实际物理量也可以通过相关公式获得.故而参考密度ρr、参考速度ur可以确定，但L′与Lr仍无法确定.鉴于此，补充如下关系式：

此外，时间t、压力p与t′、p′之间的转换，可基于下式进行求解：

至此，格子与实际物理单位间的换算完成［24］.通常，可采用下式将格子单位下的δx=δy=δz=1，δt=1 及c2=1/3，转换成物理单位.

2.4 边界条件处理

本文采用标准反弹格式［28］模拟渗流土体颗粒与流体间无滑移的渗流行为，其基本思想为：流体节点从（α-1，t）入射到（α，t）的粒子分布函数Fα1，迁移碰撞后沿原路径返回，据此获得节点（α，t）处的粒子分布函数Fα0，其余节点类似，粒子分布函数Fα0可表示为：

式中：α1 为指向壁面的方向；α0 为α1 的相反方向；xb为固体壁面处的格点；xf为流体格点，xf=xb-eαδt.

2.5 收敛性判别

参考文献［28］中收敛判别方法判断结果是否收敛.给定一个小量ε=10-6，若Error＜ε，则计算结果收敛，否则返回计算.3D判别式如下：

2.6 LBM算例验证

通过Poiseuille 流验证自编LBM 程序的正确性.选取60×60×100 的网格区域作为验证计算的3D 模型，具体计算参数如表1所示.

表1 LBM验证算例的计算参数表（格子单位）Tab.1 Calculation parameter of LBM verification example（grid unit ）

本文所有计算参数单位均为格子单位［20，30］，边界条件处理同2.4节.

对比自编LBM 程序计算结果与Poiseuille 流理论值，如图5 所示.从图5 中可知，LBM 计算值与Poiseuille流理论值基本吻合，流速最大误差仅为1.47%.

图5 Poiseuille 理论值与 LBM 计算值对比（格子单位）Fig.5 Comparison of Poiseuille theoretical value and LBM calculated value（grid unit ）

3 结果与讨论

基于QSGS重构技术和LBM对不同类型的3D土体模型进行渗流仿真模拟，设定模型初始压差：pin=1.000 6 流入，pout=1.000 0 流出，边界条件详细设定如图6所示，其他参数设置参见2.6节算例.

图6 计算模型边界条件Fig.6 Boundary conditions of the calculation of model

3.1 最优仿真模型尺寸确定方法

本文选用pc=0.01，n=0.64，尺寸为30、40、50、60、80、100、120、150 和200 的立方体模型，通过分析模型尺寸大小对渗流稳定时间和孔隙连通情况的影响来确定数值仿真最优模型尺寸的大小.渗流稳定状态判断标准［19］：当格点离散速度不再随时步增加而发生变化时，可认为渗流已经达到稳定状态，图7为渗流稳定后不同典型尺寸模型的速度场图像.

图7 不同典型尺寸模型渗流稳定后速度场切片Fig.7 Velocity field slices after seepage stability of different typical size models

采用Matlab 软件编程，统计3D 模型的孔隙参数，统计结果如表2所示.从表2中可看出：随模型尺寸扩大，孔隙连通程度明显增大，模型尺寸扩大至100×100×100 格点大小后，继续扩大模型尺寸时发现连通孔隙率nc增加并不明显，再结合图7（b）分析可知，随模型尺寸增大渗流稳定所需时间呈直线形增大.鉴于此，综合孔隙连通性和运算时间因素，本文最终选取100×100×100格点大小的多孔介质模型作为最优仿真单元模型，用于后续研究.该研究方法可为其他类似3D岩土体最优仿真单元选取提供一定的借鉴.

表2 不同尺寸模型孔隙参数表（格子单位）Tab.2 Pore parameters of different size models（grid units ）

图8 为典型模型尺寸对渗流速度和稳定时间的影响情况.从图8 中可发现：模型尺寸越大，格点速度u先随之增大后减小，而稳定所需时间t始终呈增长趋势，统计知尺寸30、40、50、60、80、100、120、150和200的模型所需时间分别为212 min、657 min、853 min、1 091 min、1 957 min、2 661 min、4 582 min、6 283 min和103 043 min（运算所采用计算机配置为 Intel Core i5-10400@3.60GHz 8G RAM，64 bit operating system，Windows10）.分析原因是：由表2 结论可知，模型尺寸越大，孔隙间的连通率随之增大，连通孔隙增多，即速度也将增大，但当模型尺寸增大到一定界限时，连通孔隙率nc将会保持在一个稳定值左右；同时随模型尺寸扩大，模型的孔隙结构也更加复杂，流体在孔隙结构中来回碰撞、迁移的次数随之增加，这将使得渗流达到稳定所需的时间增加，故速度在一定程度上会降低.

图8 3D模型尺寸对渗流时间和速度的影响Fig.8 Influence of 3D model size on seepage time and velocity

3.2 孔隙结构对渗流场的影响

以pc=0.01，n=0.64，100×100×100格点大小的模型为例进行渗流模拟，探讨孔隙结构特征对渗流过程的影响.

图9（a）～（d）分别列出了渗流稳定后不同切面的渗流速度及流线分布示意图.从图9（a）～图9（c）中可明显地发现：流体粒子自入口到出口，整个过程中格点速度呈逐渐减小趋势.同时结合三维模型孔隙结构特征分析可知，土体的孔隙结构特征对流速影响较大.如：格子坐标（1，12，46）处，格点速度为2.12×10-3，远高于其他孔隙区域速度.再观察图9（d）流线图可知，QSGS 重构的三维孔隙模型待渗流稳定后，速度分布较为集中，且存在x向主渗流通道.

图9 3D渗流模拟计算结果（pc=0.01，n=0.64）Fig.9 Results of three-dimensional seepage simulation（pc=0.01，n=0.64）

3.3 孔隙率对渗流场的影响

以pc=0.01，n=0.36、0.50、0.64、0.78，100×100×100 格点大小的模型为例，分析孔隙率变化情况下，渗流稳定后渗流速度及渗透率变化情况.

图10（a）～图10（d）分别列出了不同孔隙率的三维模型渗流稳定后典型切面的渗流速度分布示意图.从图10 中可看出：随着模型孔隙率逐渐增大，渗流流体在孔隙中的分布范围更加广泛且均匀，但不论模型孔隙率大小，速度最大值都出现在孔隙较大处.

图10 不同孔隙率的3D渗流模拟结果（pc=0.01）Fig.10 Results of 3D seepage simulation with different porosity（pc=0.01）

图11 为不同孔隙率的模型沿x向截面的平均渗流速度变化曲线.从图11 中不难发现，总体上孔隙率越大，平均渗流速度就越大，但临近出口位置处（x=80～100）变化较为杂乱，无规律可循.此外，为探究孔隙率大小与渗透率k之间的关系，绘制了孔隙率n=0.36～0.78 与渗透率的关系，并对其进行拟合，得到两者的关系表达式，如图12 所示.从图12 中可发现：计算渗透率随模型孔隙率的增大逐渐增大，且近似呈线性关系，拟合表达式为y=0.148x-0.012 2，相关系数为0.938 4.

图12 孔隙率与渗透率的拟合关系Fig.12 Fitting relationship between porosity and permeability

3.4 土颗粒大小对渗流场的影响

以pc=0.01（大颗粒）、0.05（中颗粒）、0.1（小颗粒），n=0.64，100×100×100格点大小的3D模型为例，分析土颗粒大小变化情况下，渗流场速度分布情况.

图13 为不同土颗粒大小的渗流模拟结果.从图13（a）～图13（c）中可以看出：当分布概率为pc=0.01 时，渗流场稳定后的流速主要分布范围为0～2×10-3；当分布概率为pc=0.05 时，渗流场稳定后的流速主要分布范围为0～1×10-3；当分布概率为pc=0.1时，渗流场稳定后的流速主要分布范围为0～1×10-3.并且随着土颗粒粒径减小，对应的速度场分布更为均匀，但是速度呈减小的趋势.从图13（d）～图13（f）中可以看出：当土颗粒粒径越小时，渗流场流线越密集，相应的渗流孔道也越细窄，渗流流体粒子在土体孔隙中的分布越集中.

图13 不同土颗粒大小的3D渗流模拟结果Fig.13 Results of three-dimensional seepage simulation with different soil particle sizes

为分析土体颗粒大小与速度间定量化关系，绘制不同土颗粒大小截面出口位置渗流速度分布曲线，如图14 所示.由图14 可知：当分布概率pc=0.01时，出口位置渗流速度波动在0.25×10-3～0.45×10-3区间，最大波动幅度可达1.18×10-3；当分布概率pc=0.05时，出口位置渗流速度在0.25×10-3～0.40×10-3间波动，最大波动幅度可达0.89×10-3；而当分布概率pc=0.1 时，出口位置渗流速度波动在0.25×10-3～ 0.35×10-3区间，最大波动幅度为0.88×10-3.分析上述原因是，3D重构土体颗粒越小，土体内部孔隙分布越均匀，流体渗流过程中的流速分布也更加集中.由上可见，土体中的土颗粒大小对渗流特性具有重要的作用.

图14 不同土颗粒大小出口渗流速度分布曲线Fig.14 Distribution curve of outlet seepage velocity with different soil particle size

图15 为不同土颗粒大小与渗透率的关系曲线.由图15可知，孔隙率n相同情况下，土体渗透率随土颗粒尺寸减小而减小，且分布概率pc=0.05 时，渗透率降低特别显著，降低幅度达到42.86%，当pc＞0.05时，渗透率基本稳定在0.078上下.

图15 不同土颗粒大小与渗透率的关系曲线Fig.15 Relationship curves between soil particle size and permeability

3.5 各向异性对渗流场的影响

对1.3 节各向异性3D 重构模型施加如图6 所示的边界条件，并借助LBM进行渗流场仿真计算，模拟结果如图16 和图17 所示.从图16 中可发现：各向异性模型孔隙渗流速度均小于各向同性的渗流速度，且整个过程中格点速度呈逐渐减小趋势.分析原因是：各向异性孔隙模型重构过程中，考虑因素更多，重构结果更加贴近于真实的土体孔隙结构，孔隙结根据能量守恒定律可知，相比各向同性孔隙模型渗构复杂无序，这将对流体渗流过程造成一定的阻碍，渗流速度势必会更小.

图17 3D各向异性模型出口边界渗流速度分布情况Fig.17 Seepage velocity distribution at exit boundary of 3D anisotropic soil reconstruction model

再观察图16（d）可知，各向异性重构孔隙模型待渗流稳定后，流线分布相互交错，分布均匀.此外，对比各向同性模型流线分布可以明显发现，各向异性模型的流线较为稀疏，不集中，且x向主渗流通道不明显.

图17 为各向异性3D 孔隙模型出口边界渗流速度分布情况.从图17 中可发现：考虑各向异性的模型出口边界渗流速度波动幅度较大、变化趋势规律性不明显；出口边界平均渗流速度为0.24×10-3，且存在多处格点速度为0 的区域（图17 中格点20～ 30 范围）.分析其主要原因可能是：相较于各向同性模型，各向异性模型重构所得孔隙分布复杂无序，且孔隙空间结构不均一，具有非均匀性.

综上，3D各向同/异性模型渗流场的格点速度均存自渗流入口到出口呈逐渐减小的趋势.但各向异性重构孔隙模型渗流场流线稀疏、分布相互交错，在x方向的主渗流通道不明显，且出口各个位置的流速变化无明显规律，平均流速为0.24×10-3，存在多处流速为0 的区域.而各向同性重构孔隙模型速度分布较为集中，且存在明显的x向主渗流通道，最大流速出现在孔隙较大处，可达2.12×10-3.

此外，对比文献［31］中提出的2D 各向异性模型，其在孔隙连通性好的孔道中形成了明显的主渗流通道，而本文中3D 各向异性模型在x向主渗流通道不明显，分析其主要原因可能是：模型尺寸、孔隙率、初始压差、各向异性等参数设置不同，并且3D 模型相较于2D 模型在空间结构上更加复杂、无序，因此，对流体渗流过程造成更大阻碍.

4 对比分析

周潇等［19］、申林方等［20］、Jun等［21］、蔡沛辰等［22］、李滔等［23］基于LBM模拟细微观尺度下流体在多孔介质中渗流过程，表3列出了本文与其他学者关于多孔介质重构模型微细观渗流研究结果对比情况.

表3 多孔介质渗流研究结果对比Tab.3 Comparison of research results of porous media seepage

从表3中可以发现：

1）研究者对多孔介质渗流研究一般都直接给出了研究对象的尺寸，而对于最优仿真模型尺寸选取鲜有提及.本文综合孔隙连通性和渗流稳定时间等因素选取3D 最优仿真模型尺寸为100×100×100 格点大小，可见3D 模型尺寸比2D 模型［22］尺寸小1/3左右.

2）目前渗流研究分析对象多为速度场、孔隙率和计算渗透率，本文则更加深入地研究了3D 孔隙结构与速度间的关系、土颗粒大小对渗流场影响等，发现多孔介质渗流过程中存在明显的“指进”和优先流效应；土颗粒越小，速度场分布越均匀，但速度更小.在2D 模型中最大格点速度为3.24×10-2，而3D 模型孔隙渗流最大格点速度仅为2.12×10-3，可见3D 渗流速度明显小于2D 渗流速度，分析原因可能是：3D 模型孔隙相比于2D 模型孔隙结构更为复杂，对渗流过程造成一定程度阻碍.

3）在土体细微观渗流研究领域，仅有少部分学者在重构土体模型时考虑各向异性，并且结论多为定性表述.本文基于3D 各向异性模型分析了不同位置处渗流速度大小分布情况，结论为定量化表征，并发现渗流速度u在出口边界处波动幅度较大，平均渗流速度u=0.24×10-3，且存在多处u=0的区域.

4）重构细观多孔介质模型，已被广泛用于研究多孔介质的孔隙结构及预测其水力学特征等方面，但模型能否真实反映实际对象的孔隙结构，目前尚无定论［32］.本文在前人基础上，采用改进QSGS 法生成了具有各向异性的3D 重构土体孔隙结构模型，更接近真实原状土体CT扫描模型，具有较好的适用性.

5 结论

1）3D 模型随尺寸扩大，孔隙连通程度明显增大，最大增长率为19.23%，出现在100×100×100格点大小处，继续扩大模型尺寸发现连通孔隙率nc增加并不明显.

2）流体粒子易在孔隙连通性好、孔径大的区域形成主渗流通道，且存在“指进效应”，孔道中间流速最大，可达格点速度0.002 1，越靠近孔壁渗流速度越小，接近0.

3）相同土颗粒大小下，随模型孔隙率逐渐增大，渗流流体在孔隙中的分布范围愈加广泛，平均渗流速度也越大.3D 模型计算渗透率随模型孔隙率的增大逐渐增大，且近似呈线性关系，拟合表达式为y=0.148x-0.012 2，相关系数为0.938 4.

4）孔隙率n相同的情况下，3D 土体渗透率随着土颗粒尺寸减小而减小，且分布概率pc=0～0.05 时渗透率显著降低，降低幅度达42.86%，当pc＞0.05时，渗透率基本稳定在0.078上下.

5）考虑各向异性模型渗流场的格点速度自渗流入口到出口呈逐渐减小的趋势，且在出口边界处的渗流速度波动幅度较大、变化趋势规律性不明显，平均渗流速度为0.24×10-3，同时存在多处格点速度为0的区域.此外，相较于各向同性模型流线分布，各向异性模型的流线相互交错，较为稀疏、不集中，且主渗流通道不明显.