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

第50 卷第 9 期2023年9 月
Vol.50，No.9
Sept. 2023湖南大学学报（自然科学版）
Journal of Hunan University（Natural Sciences）
基于QSGS法3D重构土体渗流场的LBM数值模拟
阙云†，邱婷，蔡沛辰，马宏岩，谢秀栋，薛斌
（福州大学土木工程学院，福建福州，350108）
摘要：基于四参数随机生长法（QSGS）以及改进QSGS方法等效重构三维（3D）各向同性/各向异性土体内部孔隙结构，结合格子玻尔兹曼方法（LBM），采用MATLAB自编程序进行细
观渗流场数值模拟，分析不同孔隙参数及模型各向异性等因素对渗流特性的影响规律. 结果
表明：3D模型尺寸100×100×100格点时，连通孔隙率n c达到最大增幅19.23%. 流体易在连通性
好、孔喉直径大的部位形成主渗流通道，且存在“指进效应”，流体越靠近孔道中轴线流速越
快. 3D重构土体的计算渗透率k随模型孔隙率n增大而增大，随土颗粒尺寸减小而减小，当分
布概率p c = 0~0.05（大、中颗粒）时渗透率显著降低，降低幅度达42.86%. 考虑各向异性的重构
模型出口边界渗流速度波动幅度较大、变化趋势规律性不明显，且流线分布稀疏、相互交错、
主渗流通道不明显. 该研究能揭示流体在3D土体孔隙结构中的流动规律，并为3D土体重构
和细观渗流模拟研究方法提供借鉴.
关键词：格子玻尔兹曼；四参数随机生长法；细观渗流；各向异性；计算渗透率
中图分类号：TU42；O357.3 文献标志码：A
LBM Numerical Simulation of Seepage Field of 3D Reconstructed Soil
Based on QSGS Method
QUE Yun†，QIU Ting，CAI Peichen，MA Hongyan，XIE Xiudong，XUE Bin
（School of Civil Engineering， Fuzhou University， Fuzhou 350108， China）
Abstract：The quartet structure generation set （QSGS）， as well as improved QSGS， were used to equivalently reconstruct the three-dimension（3D）isotropy and anisotropy pore structure of the real soil. Combined with the lattice Boltzmann method （LBM）， the numerical simulation of the seepage field employing the MATLAB automatic program was carried out to explore the influence of different pore parameters and model anisotropy on seepage characteristics. The results showed that when the mesh size of the 3D model was 100×100×100，the maximum increase of pore connectivity rate n c reached 19.23 %. The fluid preferred to form the main seepage channel in the position with good connectivity and large pore throat diameter，and there was a “finger-in”effect. The closer the fluid was to the axis of the channel，the greater the flow velocity was. The calculated permeability k of 3D reconstructed soil increased with the increase of model porosity n and decreased with the decrease of soil particle size. When the distribution probability P c was 0~0.05， the permeability k decreased significantly， and the reduction was 42.86%. The seepage velocity at the outlet boundary of the reconstructed model considering anisotropy fluctuated
∗收稿日期：2022-07-24
基金项目：国家自然科学基金资助项目（41772297）， National Natural Science Foundation of China（41772297）
作者简介：阙云（1980―），男，江西黎川人，福州大学教授，博士
† 通信联系人，E-mail：*******************.cn
文章编号：1674-2974（2023）09-0119-12DOI：10.16339/ki.hdxbzkb.2023108
湖南大学学报（自然科学版）2023 年greatly and irregularly. The streamline distribution streamline sparse and interlaced，with no prominent main seepage channel. This study can better reveal the flow law of fluid in the 3D pore structure and provide some references for the research methods of 3D soil reconstruction and meso-seepage simulation.
Key words：Lattice Boltzmann；quartet structure generation set；meso-seepage；anisotropy；calculate permeability
在基坑、边坡、堤坝、隧道等工程项目中，渗流作用往往会引起土体或结构的变形失稳，严重威胁工程安全［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)=
ì
í
î
1，x位于土体固相中；
0，x位于土体孔隙中.（1）式中：G（x）为随机变量，反映孔隙分布情况，且G（x）的平均值为<G（x）>.
120
第 9 期
阙云等：基于QSGS 法3D 重构土体渗流场的LBM 数值模拟
1.2 四参数随机生长法
四参数随机生长法［11，13］
可通过控制分布概率p c 、
生长概率p d 、概率密度p i rs （在i 方向上第r 相在第s 相
上的生长概率）和孔隙率n 四个参数重构土体细观孔隙结构. 令固相（土体颗粒）作为生长相，孔隙作为非生长相，则3D 多孔介质模型的重构流程如下：
1）在3D 空间中，令生长相以一定的分布概率p c
随机布设，p c 必须小于设定的孔隙率n .
2）在3D 空间中，令分布的生长相单元以一定的生长概率p d 朝19个运动方向生长，如图1所示.
3）重复上述步骤，当生长相达到设定的孔隙率n 时终止生长，3D 多孔介质模型重构完毕.
控制分布概率p c = 0.01（大颗粒），孔隙率n =
0.36、0.50和0.64，各向生长概率为一定值0.01（假设土体为各向同性），经过生长，得到不同孔隙率下60×
60×60格点大小的3D 模型，如图2所示.1.3 改进四参数随机生长法
为更好地研究孔隙分布各向异性对3D 多孔介质渗透能力的影响，采用两个依次递进的生长过程来构建3D 多孔介质模型，当各方向的生长速度不同时，即可获得各向异性3D 模型，详细步骤如下：
步骤1：大团颗粒生成较大的固相，固相生长核概率为p t ，生长概率为p d ，其中d 代表大团颗粒的19个生长方向，孔隙率为n 1（p t <n 1）.步骤2：小颗粒团生成较小的固相，固相生长核
概率为p s ，生长概率为p k ，其中k 代表小颗粒团的19
个生长方向，且p s 远小于p t .
两个步骤相互结合，直至达到设定的孔隙率，最
终生成具有各向异性的3D 重构土体细观模型.
各向异性的参数可通过原状土体CT 切片中孔隙的典型结构和几何形态特征确定，若同时调整各个方向的p d 可生成与CT 扫描模型良好吻合的各向异性多孔介质模型. 本文模型通过设定孔隙率n = 0.597（真实土体孔隙率平均值），生长概率p 1~7=0.01、p 8~11=0.2、
p 12~18=0.01、p 19=0.005，所生成的3D 重构土体细观模型如图3所示，图中灰色区域为孔隙、黑色区域为基质.
2 格子 Boltzmann 模型
2.1 格子 Boltzmann 方程
格子Boltzmann 方程是离散形式的Boltzmann-BGK 方程［24］，其可从细观层面仿真模拟多孔介质中的流体渗流过程. 不含外力项时，
粒子分布函数的演
图1 QSGS 生长相生长方向
Fig.1 Growth direction of QSGS growth phase
（a ） 3D 模型 （b ） 典型切片
图3 考虑各向异性的 3D 重构模型
Fig.3 Three-dimension anisotropic soil reconstruction model
（a ） n =0.36 （b ） n =0.50
（c ） n =0.64
图2 不同孔隙率的 QSGS 模型（以p c = 0.01为例）
Fig.2 QSGS model with different porosity
（taking p c = 0.01 as an example ）
121
湖南大学学报（自然科学版）2023 年
化方程可表示为：
F α(ω+e αδt ，t +δt )=F α(ω，t )-F α(ω，t )-F eq
α(ω，t )
τ
.
（2）
式中：F α（ω，t ）为格点ω位置处在t 时刻沿α方向的粒子分布函数；e α为离散速度；δt 为离散时间；τ为弛豫
时间；F eq
α(ω，t )为局部平衡态粒子分布函数. 2.2 格子 Boltzmann 基本模型
格子Boltzmann 模型通常由三部分组成，即离散
速度模型（格子）、平衡态分布函数和分布函数的演化方程［25］. 本文选择的模型是3D 层面最为常用的D3Q19模型，如图4所示.
D3Q19模型一共有如图4所示的19个运动方
向，离散速度e α满足式（3）.
e α=c éëêêê000100-100 0 0 0 1 0-1001 0 1-1 1-11-1 1-10 0 0 00-110-1-1 10 0 0 0 1-1 1-1 0 0 01-1-1 1 1-1-11ùû
úúú
.
（3）
式中：c =δx /δt ，其中δx 为网格步长，δt 为时间步长，均取1.
D3Q19模型的权系数w α如下：
w α=ìíîïïïïïï13，
α=0；118，α=1，2，⋯，6；
1
36，α=7，8，⋯，18.（4）
D d Q m 模型的局部平衡态粒子分布函数F eq
α可表
示为：
F eq
α
=ρw α[1+e α⋅u c 2s +(e α⋅u )22c 4-u 22c 2s ].
（5）
式中：c s 为格子声速，取值为
ρ为密度；w α为
D3Q19模型的权系数；u 同时基于LBE 基本模型，采用Chapman-Enskog
展开方法［26］推导所对应的Navier-Stokes 方程.密度ρ、宏观速度u 、压力差p 及运动黏度系数μ与弛豫时
间τ（无量纲）之间关系如下：
ρ=∑α=0
18F α.
（6）u =1ρ∑α=018
F αe α.
（7）p =ρc 2s .
（8）τ=μc 2s δt +12
.
（9）
2.3 真实单位与格子单位转换
针对实际问题的数值仿真，常用的编程思路有两种：程序代码中物理量直接采用实际物理单位；程序代码中物理参数都做无量纲化处理［24］. 本文采用的是后者，即无量纲化的格子单位，参考Succi 著作［27］中的方法对LBM 单位换算部分进行阐述. 模型中基本参数（格子单位）：长度L 、密度ρ、时间t 、压力差p 和运动黏度系数μ. 模型对应的参数（物理单位）：长度为L ′、密度为ρ′、时间为t ′、压力差为p ′、运
动黏度系数为μ′. 为实现上述两者的转换，需引入部分参考量，即参考长度L r 、参考密度ρr 和参考速度
u r
［24］
. 定义为：L r =L′L .（10）
ρr =ρ′ρ.（11）
u r =c′c
.（12）
式中：c ′和c 分别为物理单位和格子单位下的声速.
对于一个特定的问题，模拟中的长度L 、密度ρ、声速c 和运动黏度系数μ是已知的；实际物理量也可以通过相关公式获得. 故而参考密度ρr 、参考速度u r 可以确定，但L ′与L r 仍无法确定. 鉴于此，补充如下关系式：
L r u r =
μ′
μ
.（13）
此外，时间t 、压力p 与t ′、p ′之间的转换，可基于
图4 D3Q19 模型Fig.4 D3Q19 model
122
第 9 期阙云等：基于QSGS法3D重构土体渗流场的LBM数值模拟下式进行求解：
L r u r=t′
t=t r.（14）
p=p′t r 2
L r2ρr.（15）至此，格子与实际物理单位间的换算完成［24］. 通常，可采用下式将格子单位下的δx=δy=δz=1，δt=1及c2 =1/3，转换成物理单位.
δ′
x=δ′y=δ′z=L r.（16）δ′t=L r u r.（17）
c′=u r
3.（18）2.4 边界条件处理
本文采用标准反弹格式［28］模拟渗流土体颗粒与流体间无滑移的渗流行为，其基本思想为：流体节点从（α-1，t）入射到（α，t）的粒子分布函数Fα1，迁移碰撞后沿原路径返回，据此获得节点（α，t）处的粒子分布函数Fα0，其余节点类似，粒子分布函数Fα0可表示为：
F
α0(x b，t)=Fα1(x f，t).（19）式中：α1为指向壁面的方向；α0为α1的相反方向；x b 为固体壁面处的格点；x f为流体格点，x f= x b-eαδt.
采用非平衡态外推格式［29］模拟土体渗流压力边界处渗流行为，其基本思想为：将粒子分布函数分解为两部分，即平衡态部分F eqα(O，t)和非平衡态部分
F neq
α(O，t)，O为土体渗流压力边界处格点，平衡态F eq
α(O，t)可以由边界条件得到，而非平衡态F neqα(O，t)需通过非平衡外推确定，粒子分布函数可表示为：ì
í
î
F
α(O，t)=F eqα(O，t)+F neqα(O，t)，
F neq
α(O，t)=Fα(B，t)-F eqα(B，t).
（20）2.5 收敛性判别
参考文献［28］中收敛判别方法判断结果是否收敛. 给定一个小量ε =10-6，若Error＜ε，则计算结果收敛，否则返回计算. 3D判别式如下：
Error =)-u(ω，t)]+[u(ω，t+δ)-u(ω，t)]+[u(ω，t
∑[u x(ω，t+δt)2+u y(ω，t+δt)2+u z(ω，t+δt)2]（21）
2.6 LBM算例验证
通过Poiseuille流验证自编LBM程序的正确性. 选取60×60×100的网格区域作为验证计算的3D模型，具体计算参数如表1所示.
本文所有计算参数单位均为格子单位［20，30］，边界条件处理同2.4节.
对比自编LBM程序计算结果与Poiseuille流理论值，如图5所示. 从图5中可知，LBM计算值与Poi‑
seuille流理论值基本吻合，流速最大误差仅为1.47%.
3 结果与讨论
基于QSGS重构技术和LBM对不同类型的3D土体模型进行渗流仿真模拟，设定模型初始压差：p in= 1.000 6流入，p out=1.000 0流出，边界条件详细设定如图6所示，其他参数设置参见2.6节算例.
3.1 最优仿真模型尺寸确定方法
本文选用p c = 0.01，n = 0.64，尺寸为30、40、50、60、80、100
、120、150和200的立方体模型，通过分析模型尺寸大小对渗流稳定时间和孔隙连通情况的影响来确定数值仿真最优模型尺寸的大小. 渗流稳定状态判断标准［19］：当格点离散速度不再随时步增加而发生变化时，可认为渗流已经达到稳定状态，图7
表1 LBM验证算例的计算参数表（格子单位）
Tab.1 Calculation parameter of LBM verification example （ grid unit ）
长度L 100模型尺寸D
60
网格步长δx
1.0
时间步长δt
1.0
雷诺数Re
100
入口水压p in
1.000 6
出口水压p out
1.0
图5 Poiseuille 理论值与 LBM 计算值对比（格子单位）
Fig.5 Comparison of Poiseuille theoretical value and LBM
calculated value （ grid unit ）
123
湖南大学学报（自然科学版）2023 年
为渗流稳定后不同典型尺寸模型的速度场图像.采用Matlab 软件编程，统计3D 模型的孔隙参数，统计结果如表2所示. 从表2中可看出：随模型尺寸扩大，孔隙连通程度明显增大，模型尺寸扩大至100×100×100格点大小后，继续扩大模型尺寸时发
现连通孔隙率n c 增加并不明显，再结合图7（b ）分析可知，随模型尺寸增大渗流稳定所需时间呈直线形增大. 鉴于此，综合孔隙连通性和运算时间因素，本文最终选取100×100×100格点大小的多孔介质模型作为最优仿真单元模型，用于后续研究. 该研究方法可为其他类似3D 岩土体最优仿真单元选取提供一定的借鉴.
图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 ****************8G RAM ， 64 bit operating sys‑tem ， Windows10）. 分析原因是：由表2结论可知，模型尺寸越大，孔隙间的连通率随之增大，连通孔隙增
多，即速度也将增大，但当模型尺寸增大到一定界限时，连通孔隙率n c 将会保持在一个稳定值左右；同时随模型尺寸扩大，模型的孔隙结构也更加复杂，流体在孔隙结构中来回碰撞、迁移的次数随之增加，这将使得渗流达到稳定所需的时间增加，故速度在一定程度上会降低.
3.2 孔隙结构对渗流场的影响
以p c =0.01，n = 0.64，100×100×100格点大小的模
型为例进行渗流模拟，探讨孔隙结构特征对渗流过程的影响.
图9（a ）~（d ）分别列出了渗流稳定后不同切面的渗流速度及流线分布示意图. 从图9（a ）~图9（c ）中可明显地发现：流体粒子自入口到出口，整个过程中格点速度呈逐渐减小趋势. 同时结合三维模型孔隙结构特征分析可知，土体的孔隙结构特征对流速影响较大. 如：格子坐标（1，12，46）处，格点速度
为
图6 计算模型边界条件
Fig.6 Boundary conditions of the calculation of model
（a ） 30×30×30 速度场切片 （b ） 40×40×40 速度场切片
（c ） 50×50×50 速度场切片 （d ） 60×60×60 速度场切片
（e ） 80×80×80 速度场切片 （f ） 100×100×100 速度场切片
图7 不同典型尺寸模型渗流稳定后速度场切片
Fig.7 Velocity field slices after seepage stability of different
typical size models
表2 不同尺寸模型孔隙参数表（格子单位）Tab.2 Pore parameters of different size models
（ grid units ）
模型尺寸D
3040506080
100120150200孔隙率n
0.64
孔隙数量/个103241478
826
1 9903 9116 89910 80119 802平均孔径
3.423.443.423.433.403.393.373.633.95
连通孔隙率 n c
0.260.290.320.340.360.410.400.410.42
124
第 9 期
阙云等：基于QSGS 法3D 重构土体渗流场的LBM 数值模拟
2.12×10-3，远高于其他孔隙区域速度. 再观察图9（d ）流线图可知，QSGS 重构的三维孔隙模型待渗流稳定后，速度分布较为集中，且存在x 向主渗流通道.
3.3 孔隙率对渗流场的影响
以p c = 0.01，n = 0.36、0.50、0.64、0.78，100×100×
100格点大小的模型为例，分析孔隙率变化情况下，渗流稳定后渗流速度及渗透率变化情况.
图10（a ）~图10（d ）分别列出了不同孔隙率的三维模型渗流稳定后典型切面的渗流速度分布示意图. 从图10中可看出：随着模型孔隙率逐渐增大，渗流流体在孔隙中的分布范围更加广泛且均匀，但不论模型孔隙率大小，速度最大值都出现在孔隙较大处.
图11为不同孔隙率的模型沿x 向截面的平均渗流速度变化曲线. 从图11中不难发现，总体上孔隙率越大，平均渗流速度就越大，但临近出口位置处（x = 80~100）变化较为杂乱，无规律可循. 此外，为探
究孔隙率大小与渗透率k 之间的关系，绘制了孔隙率n = 0.36~0.78与渗透率的关系，并对其进行拟合，得到两者的关系表达式，如图12所示. 从图12中可发现：计算渗透率随模型孔隙率的增大逐渐增大，且近似呈线性关系，拟合表达式为y = 0.148x -0.012 2，相关系数为0.938 4.
3.4 土颗粒大小对渗流场的影响
以p c = 0.01（大颗粒）、0.05（中颗粒）、0.1（小颗
粒），n =0.64，100×100×100格点大小的3D
模型为例，
（a ） 十字交叉位置速度分布切片 （b ） 进口位置速度分布切片
（c ） 出口位置速度分布切片 （d ） 3D 流线图
图9 3D 渗流模拟计算结果（p c = 0.01，n = 0.64）
Fig.9 Results of three- dimensional seepage simulation
（p c = 0.01，n = 0.64
）
（a ）步长L 与格子速度u
关系
（b ）模型尺寸D 与渗流稳定时间t 的关系
图8 3D 模型尺寸对渗流时间和速度的影响
Fig.8 Influence of 3D model size on seepage time and velocity
（a ） 典型位置速度切片，n = 0.36 （b ） 典型位置速度切片，n = 0.50
（c ） 典型位置速度切片，n = 0.64 （d ） 典型位置速度切片，n = 0.78
图10 不同孔隙率的3D 渗流模拟结果（p c = 0.01）
Fig.10 Results of 3D seepage simulation with different porosity
（p c = 0.01）
125
湖南大学学报（自然科学版）2023 年
分析土颗粒大小变化情况下，渗流场速度分布情况.
图13为不同土颗粒大小的渗流模拟结果.从 图13（a ）~图13（c ）中可以看出：当分布概率为p c = 0.01时，渗流场稳定后的流速主要分布范围为0~2×10-3
；当分布概率为p c = 0.05时，渗流场稳定后的流
速主要分布范围为0~1×10-3
；当分布概率为p c =0.1时，渗流场稳定后的流速主要分布范围为0~1×10-3
. 并且随着土颗粒粒径减小，对应的速度场分布更为均匀，但是速度呈减小的趋势. 从图13（d ）~图13（f ）
中可以看出：当土颗粒粒径越小时，渗流场流线越密集，相应的渗流孔道也越细窄，渗流流体粒子在土体孔隙中的分布越集中.
为分析土体颗粒大小与速度间定量化关系，绘制不同土颗粒大小截面出口位置渗流速度分布曲线，如图14所示. 由图14可知：当分布概率p c =0.01时，出口位置渗流速度波动在0.25×10-3
~0.45×10-3
区
间，最大波动幅度可达1.18×10-3
；当分布概率p c =0.05时，出口位置渗流速度在0.25×10-3
~0.40×10-3
间波
动，最大波动幅度可达0.89×10-3；而当分布概率p c =0.1时，出口位置渗流速度波动在0.25×10-3 ~ 0.35×10-3区间，最大波动幅度为0.88×10-3. 分析上述原因是，3D 重构土体颗粒越小，土体内部孔隙分布越均匀，流体渗流过程中的流速分布也更加集中. 由上可见，
土体中的土颗粒大小对渗流特性具有重要的作用.
图15为不同土颗粒大小与渗透率的关系曲线. 由图15可知，孔隙率n 相同情况下，土体渗透率随土颗粒尺寸减小而减小，且分布概率p c = 0.05时，渗透率降低特别显著，降低幅度达到42.86%，当p c ＞0.05时，渗透率基本稳定在0.078上下.3.5 各向异性对渗流场的影响
对1.3节各向异性3D 重构模型施加如图6所示的边界条件，并借助LBM 进行渗流场仿真计算，模拟结果如图16和图17所示. 从图16中可发现：各向异性模型孔隙渗流速度均小于各向同性的渗流速度，且整个过程中格点速度呈逐渐减小趋势. 分析原因是：各向异性孔隙模型重构过程中，考虑因素更多，重构结果更加贴近于真实的土体孔隙结构，孔隙结根据能量守恒定律可知，
相比各向同性孔隙模型渗图12 孔隙率与渗透率的拟合关系
Fig.12 Fitting relationship between porosity and permeability
（a ）典型位置速度切片，p c = 0.01 （b ）典型位置速度切片，p c = 0.05
（c ）典型位置速度切片，p c = 0.1 （d ）3D 流线，p c = 0.01
（e ）3D 流线，p c = 0.05 （f ）3D 流线，p c = 0.1
图13 不同土颗粒大小的3D 渗流模拟结果
Fig.13 Results of three-dimensional seepage simulation with
different soil particle sizes
图11 平均渗流速度沿x 向截面变化曲线Fig.11 The curve of average seepage velocity
along x -direction section
126
第 9 期
阙云等：基于QSGS 法3D 重构土体渗流场的LBM 数值模拟
构复杂无序，这将对流体渗流过程造成一定的阻碍，渗流速度势必会更小.
再观察图16（d ）可知，各向异性重构孔隙模型待渗流稳定后，流线分布相互交错，分布均匀. 此外，对比各向同性模型流线分布可以明显发现，各向异性模型的流线较为稀疏，不集中，且x 向主渗流通道不明显.
图17为各向异性3D 孔隙模型出口边界渗流速度分布情况. 从图17中可发现：考虑各向异性的模型出口边界渗流速度波动幅度较大、变化趋势规律性不明显；出口边界平均渗流速度为0.24×10-3，且存在多处格点速度为0的区域（图17中格点20 ~ 30范围）. 分析其主要原因可能是：相较于各向同性模型，
各向异性模型重构所得孔隙分布复杂无序，且孔隙空间结构不均一，具有非均匀性.
综上，3D 各向同/异性模型渗流场的格点速度均存自渗流入口到出口呈逐渐减小的趋势. 但各向异性重构孔隙模型渗流场流线稀疏、分布相互交错，
在
（a ）十字交叉位置速度分布切片 （b ）进口位置流速分布切片
（c ）出口位置流速分布切片 （d ）3D 流线图
图16 3D 各向异性模型渗流模拟结果
Fig.16 Seepage simulation results of three-dimension anisotro‑
pic soil reconstruction model
（a ）出口速度分布，p c = 0.01
（b ）出口速度分布，p c = 0.05
（c ）出口速度分布，p c = 0.1
图14 不同土颗粒大小出口渗流速度分布曲线Fig.14 Distribution curve of outlet seepage velocity
with different soil particle size
图15 不同土颗粒大小与渗透率的关系曲线Fig.15 Relationship curves between soil particle size
and permeability
图17 3D 各向异性模型出口边界渗流速度分布情况Fig.17 Seepage velocity distribution at exit boundary of 3D
anisotropic soil reconstruction model
127
湖南大学学报（自然科学版）2023 年
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中可以发现：
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法生
表3 多孔介质渗流研究结果对比
Tab.3 Comparison of research results of porous media seepage
研究
尺度LBM
模型
最优模型尺寸/格子
孔隙
结构
孔隙率
颗粒
尺寸
模型同异性
细观尺度
格子单位
D2Q9
―
―
大孔隙率流速大，
小孔隙率流速小
大颗粒
流速大
各向同性
细观尺度
格子单位
D2Q9
―
流动方向受孔
隙连通性影响
―
―
各向同性
微观尺度
―
D2Q9
―
在孔隙连通性较好
的区域形成“主渗
流通道”，流速较大
多孔介质模型n越
大，渗透系数k越大
颗粒越大
渗透率越大
各向同性
细观尺度
格子单位
D2Q9
2D：300
指进现象，孔道中间
u可达0.032 4；符合
Poiseuille流理论
n越小的土体u波动
更稳定
土颗粒越大
流速越大
各向同性
细观尺度
―
D3Q19
―
―
非均质多孔介质的渗
透率k与n呈正相关
颗粒越大
渗透率越大
模型各向异性越强，
渗流通道迂曲度越小，
流动能量损失越小，
渗透率越大
细观尺度
格子单位
D3Q19
3D：100
渗流稳定后，u分布较集
中，存在x向主渗流通
道，且u可达2.12×10-3
n越大， u平均值越大，
临近出口变化杂乱、
无规律
土颗粒越小
速度越小
各向异性模型出口边
界u波动幅度较大，平
均u= 0.24×10-3，存在
多处u=0的区域
128。