第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：*******************
文章编号：1674-2974（2023）09-0119-12DOI：10.16339/jki.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 细观结构表征qnmlgb
作为一种多孔介质材料，土体通常可划分为土体颗粒（固相）和孔隙区域，其空间分布函数如下［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 多孔介质渗透能力的影响，采用两个依次递进的生长过程来构建3
D 多孔介质模型，当各方向的生长速度不同时，即可获得各向异性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
发布评论