考虑接触角及粒径级配影响的土水特征曲线滞回模型

伏映鹏; 廖红建; 吕龙龙; 柴小庆

doi:10.11779/CJGE202203012

考虑接触角及粒径级配影响的土水特征曲线滞回模型 English Version

1.
西安交通大学人居环境与建筑工程学院, 陕西西安 710049
2.
长安大学地质工程与测绘学院, 陕西西安 710054

基金项目:

国家自然科学基金项目 51879212

国家自然科学基金项目 41630639

中央高校基本科研业务费专项资金项目 xzy022019009

陕西省国际科技合作计划重点项目 2019KWZ-09

详细信息

作者简介:
伏映鹏(1992—)，男，博士研究生，主要从事非饱和土力学方面的研究。E-mail: cumtfyp@outlook.com

通讯作者:
廖红建, E-mail: hjliao@xjtu.edu.cn

中图分类号: TU431
计量
- 文章访问数: 328
- HTML全文浏览量: 39
- PDF下载量: 131
出版历程
- 收稿日期: 2021-02-23
- 网络出版日期: 2022-09-22
- 刊出日期: 2022-02-28

Hysteretic model for fitting soil-water characteristic curves considering contact angle and grain-size distribution

1.
School of Human Settlements and Civil Engineering, Xi'an Jiaotong University, Xi'an 710049, China
2.
College of Geology Engineering and Geomatics, Chang'an University, Xi'an 710054, China

摘要

摘要: 土水特征曲线描述了非饱和土中含水率（或饱和度）与基质吸力（或土水势）之间的关系。土水特征曲线存在滞回效应，即在相同的吸力下，土体在蒸发或重力排水的脱湿过程中赋存的水量比土体在入渗、毛细上升的增湿过程中所赋存的水量多。针对影响土水特征曲线滞回特性的两个关键因素——接触角及由颗粒粒径级配（孔隙分布）引起的墨水瓶效应，首先，基于势能最小原理采用数值模拟方法研究了斜面液滴行为，提出了增湿、脱湿过程微细观界面关键参数接触角与含水率的演化关系；其次，通过修正颗粒间液桥体积表达式，使Young-Laplace方程可用于计算w-lns平面上边界影响区、过渡区及部分残余区的基质吸力；最后，通过引入Rosin-Rammler颗粒粒径分布函数，建立了考虑接触角及粒径级配影响的土水特征曲线滞回模型，该模型参数物理意义明确，能够较好地描述接触角及颗粒粒径对土水特征曲线滞回特性的影响。以西安曲江某基坑压实黄土土水特征曲线为例进行模型验证，结果表明该模型能够很好地反映基质吸力与含水率的非线性关系，该模型可进一步用于非饱和土抗剪强度、本构模型等的研究。
- 土水特征曲线 /
- 滞回模型 /
- 接触角 /
- 粒径级配 /
- Young-Laplace方程 /
- Rosin-Rammler分布
Abstract: The soil-water characteristic curve (SWCC) describes the relationship between water content (or saturation) and matric suction (or soil-water potential) in unsaturated soils. It is worth noting that the amount of water stored in soils during the drying process caused by evaporation or gravity drainage is more than that during the wetting process caused by infiltration and capillary rise under the same matric suction level. This phenomenon is described by the hysteretic behavior of the SWCCs. Based on existing knowledge, the contact angle and the ink-bottle effect caused by the grain-size distribution (or pore-size distribution) are the two main factors that affect the hysteretic behaviors of the SWCCs. First, the behaviors of a liquid drop sitting on an inclined surface are investigated by the numerical method, and a relationship between the contact angle and the water content during the wetting and drying processes is proposed based on the principle of the minimum potential energy. Second, the expression for the volume of the liquid bridge is modified to make the Young-Laplace equation keep valid in calculating the matric suction in the boundary affected zone, transition zone and the part of the residual zone on the w-lns plane. Finally, an analytical model for analyzing the hysteretic behaviors of SWCCs is proposed after introducing the Roshin-Rammler grain-size distribution function. The model has a clear physical meaning, and the parameters can be obtained easily. Moreover, the model can describe the hysteretic behaviors of the SWCCs. Taking the SWCC of recompacted loess in a foundation pit in Qujiang, Xi'an as an example, the model is verified. The results show that the proposed model can well reflect the nonlinear relationship between the matric suction and the water content. The model can also be used for studying the shear strength and the constitutive model for unsaturated soils.
- soil-water characteristic curve /
- hysteretic model /
- contact angle /
- grain-size distribution /
- Young-Laplace equation /
- Rosin-Rammler distribution

HTML全文

0. 引言

淤堵按其产生的原因可分为机械淤堵、化学淤堵和生物淤堵等3种类型。岩土工程中,机械淤堵多被称为“渗透淤堵”,它指的是细小土颗粒在渗流作用下进入岩土多孔介质中滞留下来进而降低其渗透性的现象^[1]。渗透淤堵在某些工况下具有积极的作用,如坝基淤堵能减少水库的渗漏量、基坑止水帷幕淤堵能提高其截水效果等^[2-3]。但另一方面,渗透淤堵也常有不利影响甚至造成严重灾害,如减压井等工程排渗系统因淤堵而失效、斜坡体内淤堵引发滑坡、堰塞坝淤堵诱发漫顶溃决等^[4-6]。因此,研究岩土多孔介质的渗透淤堵特性对于工程实践中“趋其利、避其害”具有重要意义。

截至目前,国内外学者大多通过物理模拟试验研究岩土体或土工织物的渗透淤堵作用,如径向渗透模型试验、梯度比试验、砂槽模型试验等^[3-4,7-8]。物理模拟试验可获得试样淤堵过程中的流量、水力梯度和淤堵后的渗透系数等宏观水力参数,但无法洞悉淤堵后滞留颗粒在多孔介质中的空间分布特征及其与孔-喉结构的关系,从而也难以明晰渗透淤堵的内在机制。

岩土多孔介质中的渗透淤堵过程本质上是渗流-颗粒的相互作用即流-固耦合,数值模拟技术是研究此类问题的重要手段。Dabir等^[9]基于孔隙网络模拟了多孔介质中颗粒运移和滞留造成的淤堵。但孔隙网络模拟技术将多孔介质复杂的孔-喉结构简化为规则的球体-圆柱网络,与实际情况差别较大。李识博等^[10]基于商业离散元软件PFC中的流体动力学-离散单元法耦合（CFD-DEM）模块模拟了坝基粗粒土中的渗透淤堵作用。CFD-DEM中的流体计算基于宏观连续介质假设并采用经验或半经验的渗流模型（如达西定律）进行求解,且采用的是粗糙网格和局部平均的方案,难以进行更为精细的孔隙尺度的模拟^[11-12]。

将基于介观的格子Boltzmann方法（LBM）与离散单元法耦合形成的LBM-DEM能有效克服上述CFD-DEM的不足,更适合从孔隙尺度模拟复杂的流体-颗粒系统^[13-15]。近年来,Wang等^[16-17]采用LBM-DEM耦合算法成功地模拟了接触侵蚀和水力劈裂等复杂的岩土流-固耦合问题。为此,本文拟采用LBM模拟岩土多孔介质中的渗流,采用DEM模拟土颗粒在多孔介质中的相互作用与运动,并采用浸入运动边界法（IMB）处理移动的固体颗粒与流体之间的相互作用,基于LBM-IMB-DEM对岩土多孔介质的渗透淤堵过程开展孔隙尺度的模拟,探究多孔介质骨架颗粒分布、骨架颗粒粒径、骨架颗粒定向性、孔隙率及水压力等因素对其渗透淤堵特性的影响。

1. 数值算法

1.1 LBM基本原理

本文基于介观的LBM^[18-19]模拟岩土多孔介质中的渗流。在LBM中,流体的计算域由边长为h（格子步长）的立方体格子所划分,在每个节点上定义一个粒子分布函数 $f_{α} (r, t)$ 表示在 $r$ 位置处 $t$ 时刻以 $e_{α}$ 速度运动的流体粒子的数目。在离散的每个时间步 $Δ t$ 内,流体粒子以离散的速度 $e_{α}$ 运动到相邻的节点上,紧接着在该节点处与其它的粒子发生碰撞,碰撞后粒子的分布函数将更新,该过程称为迁移与碰撞,可用公式分别表示为

$f_{α} (r + e_{α} Δ t, t + Δ t) = f_{α} (r, t)$ ,

(1)

$f_{α} (r, t) = Ω (f_{α} (r, t))$ ,

(2)

式中, $Ω$ 为碰撞算子,它规定了流体粒子的碰撞行为。通过迁移与碰撞的不断循环,流体粒子逐步趋向平衡态分布 $f_{α}^{eq}$ 。

本文采用最为常用的碰撞算子即单松弛时间BGK碰撞模型^[20]：

$Ω (f_{α} (r, t)) = f_{α} (r, t) - \frac{Δ t}{τ} (f_{α} (r, t) - f_{α}^{eq} (r, t))$ ,

(3)

式中, $f_{α}^{eq} (r, t)$ 为局部平衡态分布函数,τ为无量纲弛豫时间。

本文LBM采用的离散速度模型为常用的D2Q9模型^[21],其离散速度、平衡态分布函数及相关计算参数为

$e_{α} = {\begin{array}{l} (0,0) (α = 0) \\ c {\cos [(α - 1) \frac{π}{2}], \sin [(α - 1) \frac{π}{2}]} (α = 1,2,3,4) \\ \sqrt{2} c {\cos [(2 α - 1) \frac{π}{4}], \sin [(2 α - 1) \frac{π}{4}]} (α = 5,6,7,8) \end{array},$

(4)

$f_{α}^{eq} = ρ ω_{α} [1 + \frac{e_{α} u}{c_{s}^{2}} + \frac{{(e_{α} u)}^{2}}{2 c_{s}^{4}} - \frac{u^{2}}{2 c_{s}^{2}}]$ ,

(5)

$c = h / Δ t$ ,

(6)

$c_{s}^{2} = c^{2} / 3$ ,

(7)

$ω_{α} = {\begin{cases} 4 / 9 (α =0) \\ 1 / 9 (α =1,2,3,4) \\ 1 / 36 (α =5,6,7,8) \end{cases}$ ,

(8)

式中,c为格子速度, $c_{s}^{}$ 为格子声速,ω_α为权系数,ρ为流体密度,u为宏观速度。

根据质量和动量守恒定律,模型的宏观密度ρ,速度u和压力p分别定义为

$ρ = \sum_{α = 0}^{8} f_{α}$ ,

(9)

$u = \frac{1}{ρ} \sum_{α = 0}^{8} f_{α} e_{α}$ ,

(10)

$p = ρ c_{s}^{2}$ 。

(11)

此外,流体的运动黏滞系数υ与无量纲弛豫时间τ及格子参数间关系为

$υ = c_{s}^{2} (τ - \frac{1}{2}) Δ t$ 。

(12)

在LBM中,还需要处理好流体的边界,如压力边界和速度边界等。另外,为计算方便一般将物理量做无量纲处理,即采用格子单位。关于边界处理及单位转换可参考文献[13,18,19]。

采用LBM计算低黏度的流体如空气和水时常出现数值不稳定,将工程中常用的大涡模拟（LES）引入LBM中可有效解决该问题,具体可参见文献[13]。

1.2 IMB方法

本文采用IMB方法^[22]处理流体流动与移动的固体颗粒之间的相互作用。在IMB中,需要先将移动的固体颗粒用流体节点进行离散并分别标记出固体内部节点、固体边界节点和流体边界节点,具体标记方法可参考文献[13]。针对这些节点,在LBM的碰撞算子 $Ω$ 中要引入1个附加碰撞项 $Ω_{α}^{S}$ 来表征固体颗粒对流体的作用,其公式为：

$Ω = - \frac{Δ t}{τ} (1 - B) [f_{α} (r, t) - f_{α}^{eq} (r, t)] + B Ω_{α}^{S}$ ,

(13)

$Ω_{α}^{S} = [f_{- α} (r, t) - f_{- α}^{eq} (ρ, u)] - [f_{α} (r, t) - f_{α}^{eq} (ρ, u_{P})]$ ,

(14)

式中,B为附加碰撞项的权函数, $- α$ 为方向 $α$ 的反方向, $u$ , $u_{P}$ 分别为节点 $r$ 处流体和固体的速度。B值由节点的固体占比 $ε$ （节点处控制单元被固体覆盖的体积分数）计算得到,

$B = \frac{ε (τ / Δ t - 0.5)}{(1 - ε) + (τ / Δ t - 0.5)}$ 。

(15)

另一方面,流体对移动的固体颗粒的作用由流体流动施加在颗粒上的力和力矩表征,

$F_{f} = c h \sum_{k = 1}^{n} (B_{k} \sum_{α = 0}^{8} Ω_{α}^{S} e_{α})$ ,

(16)

$T_{f} = c h \sum_{k = 1}^{n} [(r_{k} - r_{P}) \times (B_{k} \sum_{α = 0}^{8} Ω_{α}^{S} e_{α})]$ ,

(17)

式中,n为该固体颗粒所标记节点的数量, $(r_{k} - r_{P})$ 为从颗粒中心指向第k个节点的位置向量。

1.3 DEM基本理论

本文采用基于软球模型的DEM^[23]模拟固体颗粒之间的相互作用,即允许相接触的两个颗粒间有少量重叠。采用线性接触刚度模型和摩擦滑移模型模拟相接触的两颗粒间的本构行为（力-位移法则）,

$F_{c} = F_{c}^{n} + F_{c}^{t} = F_{n} n + F_{t} t$ ,

(18)

$F_{n} = k_{n} δ_{n}$ ,

(19)

$F_{t} = k_{t} δ_{t} \leq μ F_{n}$ ,

(20)

$T_{c} = R \times F_{c}^{t}$ ,

(21)

式中, $F_{c}$ 为颗粒间接触力, $n$ , $t$ 为接触的法向和切向单位向量, $F_{c}^{n}$ , $F_{n}$ 为法向接触力矢量及其大小, $F_{c}^{t}$ , $F_{t}$ 为切向接触力矢量及其大小, $k_{n}$ , $k_{t}$ 分别为接触法向刚度和切向刚度, $δ_{n}$ , $δ_{t}$ 分别为接触法向重叠量和切向相对位移, $μ$ 为颗粒间摩擦系数, $T_{c}$ 为接触力产生的力矩, $R$ 为从颗粒中心指向接触点的位置向量。

固体颗粒可能受到的力包括颗粒间的接触力、颗粒与固体边界墙的接触力、流体的作用力以及重力等,在这些力的作用下颗粒的运动遵循牛顿定律（运动法则）,

$m a = F_{f} + \sum_{i} F_{c i} + F_{b}$ ,

(22)

$I \ddot{θ} = \sum_{i} T_{c i} + T_{f}$ ,

(23)

式中, $m$ , $I$ 分别为颗粒的质量和惯性矩, $a$ , $\ddot{θ}$ 分别为颗粒的线加速度和角加速度, $F_{c i}$ , $T_{c i}$ 分别为第i个接触产生的接触力和力矩, $F_{b}$ 表示其它作用力如重力等。颗粒的运动方程采用显式中心有限差分法求解,以获得每个时刻颗粒的速度和位置。

1.4 LBM-IMB-DEM数值计算平台及其验证

基于1.1—1.3节所述算法的基本理论,运用Fortran语言开发了LBM-IMB-DEM耦合数值计算平台,其中LBM模拟流体流动、DEM模拟固体颗粒间的相互作用与运动、IMB模拟移动的固体颗粒与流体之间的相互作用,耦合程序的计算流程如图1所示。

图 1 LBM-IMB-DEM耦合程序计算流程

Figure 1. Flow chart of coupled LBM-IMB-DEM

下载: 全尺寸图片幻灯片

需要说明的是,LBM与DEM均采用显式时步循环算法且有各自的时间步长,一般来说DEM的时间步长 $Δ t_{D}$ 比LBM的时间步长 $Δ t$ 要小,因此可将DEM的时间步长适当缩短使得 $Δ t / Δ t_{D}$ 为整数,进而可将DEM循环作为LBM的子循环（如图1所示）以达到耦合计算时间的统一。

本文通过一系列经典算例对LBM-IMB-DEM数值计算平台进行了多方面的验证,主要包括：模拟二维Poiseuille流和圆柱绕流验证了LBM的准确性；模拟砂堆和双轴试验并与PFC模拟的结果对比验证了DEM的准确性；模拟单颗粒沉降和两个颗粒沉降DKT的过程验证了IMB的准确性。限于篇幅,这里只给出单颗粒下沉模拟的参数和结果。

单颗粒下沉算例常被用来验证颗粒-流体耦合算法的准确性。本文建立的单颗粒下沉模型中流体计算域为宽2 cm、高6 cm的长方形区域,将1个直径为0.25 cm的圆形固体颗粒放置于坐标为（1 cm,4 cm）的位置（原点位于计算域的左下角点）。固体颗粒的密度和法向刚度分别为1.25 g/cm³和1.0×10⁷ N/m,计算域四周为固体边界墙且其刚度为1.0×10⁸ N/m,不考虑颗粒与边界墙接触时的切向作用。流体的密度与运动黏度分别为1.0 g/cm³和1.0×10^-5 m²/s,流体计算域的格子步长取为0.01 cm,弛豫时间为0.65,区域四周为不透水边界并由标准反弹准则^[18]进行处理。LBM-IMB-DEM计算过程中,固体颗粒在重力作用下开始向下运动并受到向上的流体力作用,随后颗粒一直下落直到底边界并逐渐稳定下来,将该过程中颗粒竖向位置随时间的变化与Wan等^[24]基于CFD模拟的结果进行对比,如图2所示。可以看到,两者吻合较好,这说明LBM-IMB-DEM耦合算法能较准确地模拟颗粒-流体耦合问题。

图 2 单颗粒下沉模拟结果对比

Figure 2. Comparison of simulated results of single particle sedimentation

下载: 全尺寸图片幻灯片

2. 渗透淤堵的LBM-IMB-DEM模拟

实际的岩土体中颗粒粒径、颗粒形态、颗粒定向性与孔隙等内部结构异常复杂,为了便于系统地研究其淤堵特性,本文建立了简化的岩土多孔介质渗透淤堵二维模型,如图3所示。模型总长50 mm、宽20 mm,其中多孔介质区域长为40 mm（图3中2条虚线之间）,它由一定尺寸的骨架颗粒按一定孔隙率随机分布所构成,这些骨架颗粒在后续的数值计算过程中保持固定。图3中左侧绿色线框内的灰色圆盘为添加的拟淤堵颗粒,其直径为0.8 mm。模型的上、下边界不透水且为固体颗粒的墙边界,水流从左侧进入、右侧流出,左、右边界为压力边界并由非平衡态反弹法^[18]处理,上、下不透水边界及固定的骨架颗粒由标准反弹法^[18]处理。基于LBM-IMB-DEM数值计算平台对模型进行计算以模拟水流作用下拟淤堵颗粒在多孔介质中的迁移和沉积过程,主要计算参数列入表1。

图 3 岩土多孔介质渗透淤堵二维模型

Figure 3. Two-dimensional model for seepage clogging in rock and soil porous media

下载: 全尺寸图片幻灯片

表 1 渗透淤堵LBM-IMB-DEM模拟参数

Table 1. LBM-IMB-DEM simulation parameters for seepage clogging

固体参数	参数值	流体参数	参数值
颗粒密度/(kg·m^-3)	2650	流体密度/(kg·m^-3)	1000
法向刚度/(N·m^-1)	1×10⁷	运动黏度/(m²·s^-1)	1×10^-6
切向刚度/(N·m^-1)	5×10⁶	弛豫时间	0.51
摩擦系数	0.3	格子步长/m	8×10^-5
DEM子循环步数	4	LBM时间步长/s	2.13×10^-5

下载: 导出CSV

| 显示表格

本文定义“淤堵渗透系数比”来表征多孔介质渗透淤堵的严重程度,其定义为多孔介质淤堵后渗透系数与淤堵前渗透系数的比值,该值越小表示淤堵越严重。另外,还定义了“滞留体积百分比”来表征淤堵后多孔介质中滞留颗粒数量的多少,其定义为多孔介质中滞留颗粒的总体积与孔隙体积之比（百分数）。

待颗粒运动与流场都稳定后停止淤堵过程的模拟。随后,针对淤堵前、后的多孔介质模型（图3中2条红色线之间区域）,在保持其中颗粒（包括淤堵颗粒）位置不变的条件下,将其入口与出口密度分别设为1000.1和1000.0（雷诺数较小,基本为层流）并进行迭代计算,当流场收敛到稳态后停止迭代并由达西定律计算其渗透系数,据此可进一步得到该模型的淤堵渗透系数比。需要指出的是下文中所涉及的流场图也均是在此条件下获得的。

3. 模拟结果分析

本文主要模拟了不同骨架颗粒分布、骨架颗粒粒径、多孔介质孔隙率、骨架颗粒定向性及水流入口压力（水流入口密度）等条件下多孔介质的渗透淤堵,下面就这些条件下的模拟结果进行分析。需要说明的是,格子单位与物理单位的比值为一个固定的参考量,也就是说格子单位的变化规律与物理单位的变化规律是完全一致的,一般为了计算方便在计算过程中统一取格子单位。为方便起见,在下文的分析中也直接采用了格子单位（l.u.）下得到的模拟结果。

3.1 骨架颗粒分布的影响

在骨架颗粒是粒径（直径）为37.5（3 mm）的圆形,孔隙率为0.7,水流入口与出口密度分别为1001.0和1000.0的条件下,通过改变随机数种子生成了5种不同骨架颗粒随机分布的多孔介质模型并进行渗透淤堵过程的LBM-IMB-DEM数值模拟。模拟结果如下：5个模型的淤堵渗透系数比分别为0.41,0.36,0.43,0.37和0.43,平均值为0.40；5个模型的滞留体积百分比分别为3.82%,2.53%,2.51%,2.82%和2.55%,平均值为2.85%。其中第1,2种骨架颗粒分布的模型淤堵前后的流场与滞留颗粒分布情况如图4所示。

图 4 不同骨架颗粒分布时多孔介质淤堵前、后的流场与滞留颗粒分布

Figure 4. Distribution of retained particles and rate of fluid flow in porous media versus distribution of skeleton particles

下载: 全尺寸图片幻灯片

由模拟结果及图4可知：渗透淤堵发生后,尽管多孔介质中滞留颗粒体积只占总孔隙体积很小一部分（小于5.0%）,但这些颗粒多驻留在孔喉处,致使流动通道数显著减少,因而淤堵后的渗透系数仅为淤堵前的40%；骨架颗粒分布1的模型中滞留颗粒相对于分布2时的较多,但部分滞留颗粒位于离入口较近的骨架颗粒表面和非流动区域,而分布2的滞留颗粒基本都位于孔喉处,所以分布2的淤堵渗透系数比反而较小,即淤堵更为严重。

由上述分析可得,多孔介质中骨架颗粒的空间分布不同,其孔隙结构也不同,渗透淤堵特性也有所区别。为此,在下文分析不同骨架颗粒粒径等条件下（除不同水流入口压力的情形）多孔介质渗透淤堵时,淤堵前渗透系数、淤堵渗透系数比和滞留体积百分比均采用5个随机分布模型的平均值。

3.2 骨架颗粒粒径的影响

在孔隙率为0.7,骨架颗粒是圆形,水流入口与出口密度分别为1001.0和1000.0的条件下,分别生成骨架颗粒粒径为25,37.5和50的多孔介质模型并进行渗透淤堵过程的LBM-IMB-DEM数值模拟。各模型淤堵后的淤堵渗透系数比、滞留体积百分比及淤堵前的渗透系数如图5所示,淤堵前、后的流场和滞留颗粒分布情况如图6所示。

图 5 不同骨架颗粒粒径时多孔介质的淤堵渗透系数比、滞留体积百分比和淤堵前渗透系数

Figure 5. Variation of permeability ratio of clogging, percentage of retention volume and permeability coefficient before clogging with size of skeleton particles

下载: 全尺寸图片幻灯片

图 6 不同骨架颗粒粒径时多孔介质淤堵前、后的流场与滞留颗粒分布

Figure 6. Distribution of retained particles and rate of fluid flow in porous media versus size of skeleton particles

下载: 全尺寸图片幻灯片

由图5,6可知,随着骨架颗粒粒径的增大,淤堵前多孔介质中的流动通道数虽然减少,但其宽度显著增加,透水性也随之增大,细颗粒在水流携带下也更容易通过,故淤堵后滞留体积百分比逐渐减小,淤堵渗透系数比不断增大。因此,骨架颗粒粒径越小,多孔介质越容易发生淤堵,淤堵的程度也更严重,如图6（a）中细颗粒很快就堵塞了入口附近的孔喉,后续颗粒继续堆积在入口处,最终所有流动通道都被完全堵塞。

3.3 多孔介质孔隙率的影响

在骨架颗粒是37.5的圆形、水流入口与出口密度分别为1001.0和1000.0的条件下,分别生成孔隙率为0.6,0.7和0.8的多孔介质模型并进行渗透淤堵过程的LBM-IMB-DEM数值模拟。各孔隙率的模型淤堵前渗透系数及淤堵后的淤堵渗透系数比和滞留体积百分比如图7所示,淤堵前、后的流场和滞留颗粒分布情况如图8所示。

图 7 不同孔隙率时多孔介质的淤堵渗透系数比、滞留体积百分比和淤堵前渗透系数

Figure 7. Variation of permeability ratio of clogging, percentage of retention volume and permeability coefficient before clogging with porosity of porous media

下载: 全尺寸图片幻灯片

图 8 不同孔隙率时多孔介质淤堵前、后的流场与滞留颗粒分布

Figure 8. Distribution of retained particles and rate of fluid flow in porous media versus porosity

下载: 全尺寸图片幻灯片

由图7,8可见,随着孔隙率的增加,多孔介质中孔隙变大,连通性也更好,故其透水性增强,细颗粒在水流携带下也更容易通过,滞留体积百分比随之减小,淤堵渗透系数比逐渐增大。因此,孔隙率越小,多孔介质越容易发生淤堵。

3.4 骨架颗粒定向性的影响

为了探究多孔介质中骨架颗粒定向性对其渗透淤堵特性的影响,本文建立了如图9所示的条形颗粒,其由3个同面积的圆组合而成且颗粒面积等于直径为37.5的圆形颗粒的面积。在骨架颗粒为条形、多孔介质孔隙率为0.7、水流入口与出口密度分别为1001.0和1000.0的条件下,分别生成颗粒长轴方向平行于流动方向、垂直于流动方向和随机的3种多孔介质模型。

图 9 颗粒形状示意图

Figure 9. Schematic diagram of particle shapes

下载: 全尺寸图片幻灯片

通过渗透淤堵过程的LBM-IMB-DEM数值模拟得到这3种模型淤堵后的淤堵渗透系数比和滞留体积百分比以及淤堵前渗透系数,如图10所示,淤堵前、后的流场和滞留颗粒分布情况如图11所示。

图 10 不同骨架颗粒定向性时多孔介质的淤堵渗透系数比、滞留体积百分比和淤堵前渗透系数

Figure 10. Variation of permeability ratio of clogging, percentage of retention volume and permeability coefficient before clogging with orientation of skeleton particles

下载: 全尺寸图片幻灯片

图 11 不同骨架颗粒定向性时多孔介质淤堵前、后的流场与滞留颗粒分布

Figure 11. Distribution of retained particles and rate of fluid flow in porous media versus orientation of skeleton particles

下载: 全尺寸图片幻灯片

由图10,11可见：骨架颗粒长轴平行于流动方向时,多孔介质中流动通道较宽、连通性较好且迂曲度较小,故其透水性最强,细颗粒在水流携带下也最容易通过,淤堵渗透系数比最大；骨架颗粒长轴垂直于流动方向时,多孔介质中流动方向的孔喉直径较小且渗流通道的迂曲度较大,故其淤堵前透水性最弱,最后被细颗粒完全堵塞；骨架颗粒长轴方向随机时多孔介质的淤堵渗透系数比居于前两者之间。

从图5,7,10可以看出,表征多孔介质渗透淤堵程度的淤堵渗透系数比的变化规律与淤堵前渗透系数的变化是一致的,这也表明透水性强的多孔介质更不容易发生淤堵。然而,淤堵渗透系数比和滞留体积百分比的变化规律并不完全一致,如前所述的图4中两种情形,图10中也出现了局部的不一致情况,这表明多孔介质中滞留颗粒越多并不一定代表其淤堵越严重。

3.5 水压力的影响

为了探究水压力对多孔介质渗透淤堵特性的影响,将图4（b）中所示多孔介质模型分别在水流入口密度为1000.5,1000.75和1001.0（出口密度均为1000.0）的条件下进行渗透淤堵过程的LBM-IMB-DEM模拟。这3种条件下多孔介质淤堵后的淤堵渗透系数比和滞留体积百分比如图12所示。

图 12 不同水流入口密度时多孔介质的淤堵渗透系数比、滞留体积百分比

Figure 12. Variation of permeability ratio of clogging and percentage of retention volume with inlet density

下载: 全尺寸图片幻灯片

由图12可见,当水流出口密度保持不变,随着水流入口密度的增加,多孔介质淤堵后的滞留体积百分比逐渐增加,淤堵渗透系数比逐渐减小。究其原因,水流入口的密度越大即水压力越大（式（11））,流速也就越大,在较大的流速下入口处的细颗粒运动速度较快并趋于在多孔介质中分散开来,这样其驻留在孔喉处的概率更大,最终滞留的颗粒更多,淤堵也更为严重。如图13所示,入口密度为1001.0时细颗粒堵塞了多孔介质上下边界附近的流动通道,但入口密度为1000.5时这些部位并未发生堵塞。

图 13 不同水流入口密度时多孔介质淤堵前、后的流场与滞留颗粒分布

Figure 13. Distribution of retained particles and rate of fluid flow in porous media versus inlet density

下载: 全尺寸图片幻灯片

4. 结论

本文基于LBM-IMB-DEM耦合算法从孔隙尺度模拟了岩土多孔介质渗透淤堵的过程,分析了多孔介质骨架颗粒分布、骨架颗粒粒径、孔隙率、骨架颗粒定向性以及水压力等因素对多孔介质淤堵后滞留体积百分比和淤堵渗透系数比的影响,并通过滞留颗粒在多孔介质中的分布及相应流场的变化情况探究了各因素影响的内在机制。

（1）LBM-IMB-DEM耦合算法能较准确地模拟颗粒-流体动力耦合问题,在岩土流-固耦合精细模拟等研究领域具有广阔的应用前景。

（2）岩土多孔介质中骨架颗粒的空间分布不同,其孔隙结构特征也不同,渗透淤堵特性也有所区别。

（3）岩土多孔介质的骨架颗粒粒径或孔隙率越小,渗透淤堵后的滞留颗粒越多,淤堵也更为严重。

（4）岩土多孔介质中骨架颗粒的定向性对其渗透淤堵特性有一定的影响,当骨架颗粒为条形且其长轴与流动方向垂直时淤堵最为严重。

（5）随水压力增大,渗透淤堵后滞留颗粒在多孔介质中更趋于分散,淤堵更为严重。

（6）岩土多孔介质渗透淤堵的程度与其淤堵前的渗透性密切相关,表现为透水性强的多孔介质相对不易发生淤堵。

后续研究还应进一步考虑其他因素如拟淤堵颗粒的粒径与形状、骨架颗粒的粒径分布,骨架颗粒的形状、流体的黏度等对多孔介质渗透淤堵特性的影响,并应进一步开发三维的LBM-DEM耦合计算程序以便开展渗透淤堵过程的三维模拟。

图 1 SE数值模型点、线、面示意图

Figure 1. Schematic diagram of points, lines and surfaces in SE numerical model

下载: 全尺寸图片幻灯片

图 2 不同迭代步下的液滴能量及其形状

Figure 2. Energy and shapes of liquid drop under different iteration steps

下载: 全尺寸图片幻灯片

图 3 液滴位于不同倾角时的形状及接触角变化规律

Figure 3. Shapes of liquid drop and contact angles under different inclined angles

下载: 全尺寸图片幻灯片

图 4 接触角演化模型

Figure 4. Model for evolution of contact angle

下载: 全尺寸图片幻灯片

图 5 液桥形态的变化示意图

Figure 5. Change of water morphology with water content increasing

下载: 全尺寸图片幻灯片

图 6 含水率与结合水膜厚度及颗粒平均半径之间的关系（N_c = 9，G=2.7）

Figure 6. Relationship among water content, water film thickness and mean particle radius(N_c = 9, G=2.7)

下载: 全尺寸图片幻灯片

图 7 增湿过程前进接触角与土水特征曲线关系

Figure 7. Relationship between advancing contact angles and SWCCs during wetting

下载: 全尺寸图片幻灯片

图 8 脱湿过程后退接触角与土水特征曲线关系

Figure 8. Relationship between receding contact angles and SWCCs during drying

下载: 全尺寸图片幻灯片

图 9 4种不同的粒径级配曲线

Figure 9. Four different grain-size distribution curves

下载: 全尺寸图片幻灯片

图 10 粒径级配参数对土水特征曲线的影响（λ=0.00002）

Figure 10. Effects of parameters of grain-size distribution curves on soil-water characteristic curves (λ=0.00002)

下载: 全尺寸图片幻灯片

图 11 颗粒粒径级配（a=0.9236, b=0.0142）

Figure 11. Grain-size distribution curve(a=0.9236, b=0.0142)

下载: 全尺寸图片幻灯片

图 12 模型验证结果

Figure 12. Results of model validation

下载: 全尺寸图片幻灯片

表 1 黄土的物理性质指标

Table 1 Physical properties of test loess

相对质量密度^*	含水率	孔隙比	黏粒^**	粉粒^**	砂粒^**
2.71	16.61%	0.79	18.40%	81.30%	0.03%
注：^实验室所测相对质量密度范围为2.70~2.72，本文计算时使用值为2.71，^*粒组划分标准为：黏粒（d＜0.002 mm），粉粒（0.002~0.075 mm），砂粒（d＞0.075 mm），采用文献[49，50]计算时，应按照文献中的标准重新确定。

下载: 导出CSV

参考文献(55)

[1]	高游, 孙德安, 张俊然, 等. 考虑孔隙比和水力路径影响的非饱和土土水特性研究[J]. 岩土工程学报, 2019, 41(12): 2191–2196. http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract18065.shtml GAO You, SUN De-an, ZHANG Jun-ran, et al. Soil-water characteristics of unsaturated soils considering initial void ratio and hydraulic path[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(12): 2191–2196. (in Chinese) http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract18065.shtml
[2]	ANANDARAJAH A, AMARASINGHE P M. Microstructural investigation of soil suction and hysteresis of fine-grained soils[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2012, 138(1): 38–46. doi: 10.1061/(ASCE)GT.1943-5606.0000555
[3]	刘艳, 于建涛. 动态土水特征曲线滞后模型研究[J]. 岩土工程学报, 2021, 43(1): 62–68. http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract18420.shtml LIU Yan, YU Jian-tao. Hysteresis model for soil-water characteristic curve under dynamic conditions[J]. Chinese Journal of Geotechnical Engineering, 2021, 43(1): 62–68. (in Chinese) http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract18420.shtml
[4]	ARYA L M, PARIS J F. A physicoempirical model to predict the soil moisture characteristic from particle‐size distribution and bulk density data[J]. Soil Science Society of America Journal, 1981, 45(6): 1023–1030. doi: 10.2136/sssaj1981.03615995004500060004x
[5]	FREDLUND M D, WILSON G W, FREDLUND D G. Use of the grain-size distribution for estimation of the soil-water characteristic curve[J]. Canadian Geotechnical Journal, 2002, 39(5): 1103–1117. doi: 10.1139/t02-049
[6]	ASSOULINE S, TESSIER D, BRUAND A. A conceptual model of the soil water retention curve[J]. Water Resources Research, 1998, 34(2): 223–231. doi: 10.1029/97WR03039
[7]	郑方, 刘奉银, 王磊. 粒度对非饱和土土水特征曲线滞回特性的影响[J]. 水利与建筑工程学报, 2019, 17(5): 19–24. https://www.cnki.com.cn/Article/CJFDTOTAL-FSJS201905003.htm ZHENG Fang, LIU Feng-yin, WANG Lei. Influence of unsaturated soil granularity on hysteretic behavior of soil-water characteristic curve[J]. Journal of Water Resources and Architectural Engineering, 2019, 17(5): 19–24. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-FSJS201905003.htm
[8]	陈宇龙, 内村太郎. 粒径对土持水性能的影响[J]. 岩石力学与工程学报, 2016, 35(7): 1474–1482. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201607018.htm CHEN Yu-long Uchimura Taro. Influence of particle size on water retention of soils[J]. Chinese Journal of Rock Mechanics and Engineering, 2016, 35(7): 1474–1482. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201607018.htm
[9]	朱青青, 苗强强, 陈正汉, 等. 考虑基质势影响的非饱和土水分运移规律测试系统研制[J]. 岩土工程学报, 2016, 38(增刊2): 240–244. http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract16687.shtml ZHU Qing-qing, MIAO Qiang-qiang, CHEN Zheng-han, et al. Development of test system for unsaturated soil water movement law considering influence of matrix potential[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(S2): 240–244. (in Chinese) http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract16687.shtml
[10]	姚仰平, 王琳. 影响锅盖效应因素的研究[J]. 岩土工程学报, 2018, 40(8): 1373–1382. http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract17471.shtml YAO Yang-ping, WANG Lin. Influence factors for "pot-cover effect"[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(8): 1373–1382. (in Chinese) http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract17471.shtml
[11]	VAN GENUCHTEN M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980, 44(5): 892–898. doi: 10.2136/sssaj1980.03615995004400050002x
[12]	FREDLUND D G, XING A. Equations for the soil-water characteristic curve[J]. Canadian Geotechnical Journal, 1994, 31(4): 521–532. doi: 10.1139/t94-061
[13]	BROOKS R H, COREY A T. Hydraulic properties of porous media and their relation to drainage design[J]. Transactions of the ASAE, 1964, 7(1): 26–0028. doi: 10.13031/2013.40684
[14]	BRUTSAERT W. Some methods of calculating unsaturated permeability[J]. Transactions of the ASAE, 1967, 10(3): 400–404. doi: 10.13031/2013.39683
[15]	ZAPATA C E, HOUSTON W N, HOUSTON S L, et al. Soil–water characteristic curve variability[J]. Advances in Unsaturated Geotechnics, 2000: 84–124.
[16]	ZHAI Q, RAHARDJO H. Quantification of uncertainties in soil–water characteristic curve associated with fitting parameters[J]. Engineering Geology, 2013, 163: 144–152. doi: 10.1016/j.enggeo.2013.05.014
[17]	LIKOS W J, LU N, GODT J W. Hysteresis and uncertainty in soil water-retention curve parameters[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2013, 140(4): 04013050.
[18]	LI P, LI T L, VANAPALLI S K. Prediction of soil-water characteristic curve for Malan loess in Loess Plateau of China[J]. Journal of Central South University, 2018, 25(2): 432–447. doi: 10.1007/s11771-018-3748-1
[19]	BAYAT H, MAZAHERI B, MOHANTY B P. Estimating soil water characteristic curve using landscape features and soil thermal properties[J]. Soil and Tillage Research, 2019, 189: 1–14. doi: 10.1016/j.still.2018.12.018
[20]	TOMASELLA J, HODNETT M G. Estimating soil water retention characteristics from limited data in Brazilian Amazonia[J]. Soil Science, 1998, 163(3): 190–202. doi: 10.1097/00010694-199803000-00003
[21]	PHAM K, KIM D, YOON Y, et al. Analysis of neural network based pedotransfer function for predicting soil water characteristic curve[J]. Geoderma, 2019, 351: 92–102. doi: 10.1016/j.geoderma.2019.05.013
[22]	WEI C, DEWOOLKAR M M. Formulation of capillary hysteresis with internal state variables[J]. Water Resources Research, 2006, 42(7): 260–273.
[23]	LI X S. Modelling of hysteresis response for arbitrary wetting/drying paths[J]. Computers and Geotechnics, 2005, 32(2): 133–137. doi: 10.1016/j.compgeo.2004.12.002
[24]	NIMMO J R. Semiempirical model of soil water hysteresis[J]. Soil Science Society of America Journal, 1992, 56(6): 1723–1730. doi: 10.2136/sssaj1992.03615995005600060011x
[25]	GAN Y, MAGGI F, BUSCARNERA G, et al. A particle–water based model for water retention hysteresis[J]. Géotechnique Letters, 2013, 3(4): 152–161. doi: 10.1680/geolett.13.00046
[26]	MUALEM Y. Extension of the similarity hypothesis used for modeling the soil water characteristics[J]. Water Resources Research, 1977, 13(4): 773–780. doi: 10.1029/WR013i004p00773
[27]	MUALEM Y. Prediction of the soil boundary wetting curve[J]. Soil Science, 1984, 137(6): 379–390. doi: 10.1097/00010694-198406000-00001
[28]	刘艳, 赵成刚, 李舰, 等. 相间交界面对非饱和土应力状态的影响[J]. 力学学报, 2017, 49(2): 335–343. https://www.cnki.com.cn/Article/CJFDTOTAL-LXXB201702010.htm LIU Yan, ZHAO Cheng-gang, LI Jian, et al. The influence of interfaces on the stress state in unsaturated soils[J]. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(2): 335–343. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-LXXB201702010.htm
[29]	张鹏程, 汤连生, 姜力群, 等. 基质吸力与含水量及干密度定量关系研究[J]. 岩石力学与工程学报, 2013, 31(增刊1): 2792–2797. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2013S1027.htm ZHANG Peng-cheng, TANG Lian-sheng, JIANG Li-qun, et al. Research of quantitative relations of matric suction with water content and dry density[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 31(S1): 2792–2797. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2013S1027.htm
[30]	ZHOU A N. A contact angle-dependent hysteresis model for soil–water retention behaviour[J]. Computers and Geotechnics, 2013, 49: 36–42. doi: 10.1016/j.compgeo.2012.10.004
[31]	LIKOS W J, LU N. Hysteresis of capillary stress in unsaturated granular soil[J]. Journal of Engineering Mechanics, 2004, 130(6): 646–655. doi: 10.1061/(ASCE)0733-9399(2004)130:6(646)
[32]	WANG X, LI J. A novel liquid bridge model for estimating SWCC and permeability of granular material[J]. Powder Technology, 2015, 275: 121–130. doi: 10.1016/j.powtec.2015.01.044
[33]	SOLIGNO G, DIJKSTRA M, VAN ROIJ R. The equilibrium shape of fluid-fluid interfaces: derivation and a new numerical method for Young's and Young-Laplace equations[J]. The Journal of Chemical Physics, 2014, 141(24): 244702. doi: 10.1063/1.4904391
[34]	BRAKKE K A. The surface evolver[J]. Experimental Mathematics, 1992, 1(2): 141–165. doi: 10.1080/10586458.1992.10504253
[35]	CHOU T H, HONG S J, SHENG Y J, et al. Drops sitting on a tilted plate: receding and advancing pinning[J]. Langmuir, 2012, 28(11): 5158–5166. doi: 10.1021/la300257t
[36]	BEREJNOV V, THORNE R E. Enhancing drop stability in protein crystallization by chemical patterning[J]. Acta Crystallographica Section D, Biological Crystallography, 2005, 61(12): 1563–1567. doi: 10.1107/S0907444905028866
[37]	THAMPI S P, GOVINDARAJAN R. Minimum energy shapes of one-side-pinned static drops on inclined surfaces[J]. Physical Review E, Statistical, Nonlinear, and Soft Matter Physics, 2011, 84(4): 046304. doi: 10.1103/PhysRevE.84.046304
[38]	SANTOS M J, VELASCO S, WHITE J A. Simulation analysis of contact angles and retention forces of liquid drops on inclined surfaces[J]. Langmuir, 2012, 28(32): 11819–11826. doi: 10.1021/la3019293
[39]	FISHER R A. On the capillary forces in an ideal soil; correction of formulae given by WB Haines[J]. The Journal of Agricultural Science, 1926, 16(3): 492–505. doi: 10.1017/S0021859600007838
[40]	NISHIYAMA N, YOKOYAMA T. Water film thickness in unsaturated porous media: effect of pore size, pore solution chemistry, and mineral type[J]. Water Resources Research, 2021, 57, e2020WR029257.
[41]	SMITH W O, FOOTE P D, BUSANG P F. Packing of homogeneous spheres[J]. Physical Review, 1929, 34(9): 1271–1274. doi: 10.1103/PhysRev.34.1271
[42]	BUTT H J, KAPPL M. Normal capillary forces[J]. Advances in Colloid and Interface Science, 2009, 146(1/2): 48–60.
[43]	黎澄生, 孔令伟, 柏巍, 等. 土–水特征曲线滞后阻塞模型[J]. 岩土力学, 2018, 39(2): 59–604. LI Cheng-sheng, KONG Ling-Wei, Bai Wei, et al. Hysteresis model of soil-water characteristic curve[J]. Rock and Soil Mechanics, 39(2): 598–604. (in Chinese)
[44]	KASANGAKI G J, MEDERO G M, OOI J Y. Factors Influencing Water Retention Characteristics of Granular Materials[M]// Multiphysical Testing of Soils and Shales, Berlin: Springer, 2013.
[45]	ESMAEELNEJAD L, SIAVASHI F, SEYEDMOHAMMADI J, et al. The best mathematical models describing particle size distribution of soils[J]. Modeling Earth Systems and Environment, 2016, 2(4): 1–11.
[46]	ROSIN P, RAMMLER E. The laws governing the fineness of powdered coal[J]. Journal of the Institute of Fuel, 1933, 7: 29–36.
[47]	井彦林, 王昊, 陶春亮, 等. 非饱和黄土的接触角与孔隙特征试验[J]. 煤田地质与勘探, 2019, 47(5): 157–162. doi: 10.3969/j.issn.1001-1986.2019.05.022 JING Yan-lin, WANG Hao, TAO Chun-liang, et al. Experimental study on contact angle and pore characteristics of unsaturated loess[J]. Coal Geology & Exploration, 2019, 47(5): 157–162. (in Chinese) doi: 10.3969/j.issn.1001-1986.2019.05.022
[48]	陶高梁, 朱学良, 胡其志, 等. 黏性土压缩过程临界孔径现象及固有分形特征[J]. 岩土力学, 2019, 40(1): 81–90. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201901005.htm TAO Gao-liang, ZHU Xue-liang, HU Qi-zhi, et al. Critical pore-size phenomenon and intrinsic fractal characteristic of clay in process of compression[J]. Rock and Soil Mechanics, 2019, 40(1): 81–90. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201901005.htm
[49]	SHIRAZI M A, BOERSMA L, HART J W. A unifying quantitative analysis of soil texture: improvement of precision and extension of scale[J]. Soil Science Society of America Journal, 1988, 52(1): 181–190. doi: 10.2136/sssaj1988.03615995005200010032x
[50]	SHIRAZI M A, BOERSMA L. A unifying quantitative analysis of soil texture[J]. Soil Science Society of America Journal, 1984, 48(1): 142–147. doi: 10.2136/sssaj1984.03615995004800010026x
[51]	刘建维. 面向润湿作用的黄土结构性及吸附特性试验研究[D]. 西安: 长安大学, 2017. LIU Jian-wei. Experimental Study on the Structure and Adsorption Characteristics of Loess oriented to Wetting Function[D]. Xi'an: Changan University, 2017. (in Chinese)
[52]	栾茂田, 李顺群, 杨庆. 非饱和土的理论土-水特征曲线[J]. 岩土工程学报, 2005, 27(6): 611–615. http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract11688.shtml LUAN Mao-tian, LI Shun-qun, YANG Qing. Theoretical soil-water characteristic curve for unsaturated soils[J]. Chinese Journal of Geotechnical Engineering, 2005, 27(6): 611–615. (in Chinese) http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract11688.shtml
[53]	刘星志, 吴悦, 潘诗婷, 等. 颗粒级配对非饱和红土土-水特征曲线的影响[J]. 水利水运工程学报, 2018(5): 103–110. https://www.cnki.com.cn/Article/CJFDTOTAL-SLSY201805015.htm LIU Xing-zhi, WU Yue, PAN Shi-ting, et al. Influences of different grain size contents on soil-water characteristic curve of unsaturated laterite based on fractal theory[J]. Hydro-Science and Engineering, 2018(5): 103–110. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLSY201805015.htm
[54]	穆青翼, 党影杰, 董琪, 等. 原状和压实黄土持水特性及湿陷性对比试验研究[J]. 岩土工程学报, 2019, 41(8): 1496–1504. http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract17975.shtml MU Qing-yi, DANG Ying-jie, DONG Qi, et al. Water-retention characteristics and collapsibity behaviors: comparison between intact and compacted loesses[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(8): 1496–1504. (in Chinese) http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract17975.shtml
[55]	伏映鹏, 廖红建, 刘雪刚, 等. 非线性连续卸荷路径下黄土的强度与变形特性研究[J]. 重庆大学学报, 2021, 44(5): 26–37. https://www.cnki.com.cn/Article/CJFDTOTAL-FIVE202105004.htm FU Ying-peng, LIAO Hong-jian, LIU Xue-gang, et al. Study on strength and deformation characteristics of loess under nonlinear continuous unloading path[J]. Journal of Chongqing University, 2021, 44(5): 26–37. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-FIVE202105004.htm