考虑孔隙尺寸分布演化的持水曲线和相对渗透系数模型

Modeling Water-Retention and Relative Hydraulic Conductivity Curves Considering Pore-Size Distribution Evolution

摘要

孔隙尺寸分布（PSD）显著影响非饱和土的持水和渗流特性。假设土体PSD服从分形定律，采用分形维数定量表征土体体积变形对PSD的影响，通过结合孔隙比和分形维数之间的内在约束，提出了可描述体积变形过程中PSD演变的分形方程，推导了考虑PSD演变的非饱和土持水曲线模型（SWRC）和相对渗透系数模型（RHCC）。模型参数可通过两个不同初始孔隙比下土体的持水曲线数据进行校准，然后用于预测任意初始孔隙比下的SWRC和RHCC。最后，采用7组试验数据对提出的模型进行了验证，结果表明提出的基于分形理论的水力模型可较好地描述体积变形对非饱和土持水和渗流特性的影响。

Abstract

Pore size distribution （PSD） significantly affects the water-retention and seepage characteristics of unsaturated soils. By assuming that the PSD of soil follows a fractal law and that soil volume deformation affects the PSD via changes in the fractal dimension， a fractal-based equation describing the evolution of PSD during soil volume deformation is proposed. This equation combines intrinsic constraints between the pore size ratio and the fractal dimension. Consequently， models for the unsaturated soil water retention curve （SWRC） and the relative hydraulic conductivity curve （RHCC） that consider the evolution of PSD are derived. The model parameters can be calibrated using water retention curve data from soil samples with two different initial porosities， enabling prediction of the SWRC and RHCC under various deformation states. Finally， the proposed model is validated against seven sets of experimental data， demonstrating that the hydraulic model based on fractal theory effectively describes the impact of volume deformation on the water retention and seepage characteristics of unsaturated soils.

Graphical abstract

关键词

非饱和土 / 持水特性 / 渗透系数 / 孔隙尺寸分布

Key words

unsaturated soil / water retention / hydraulic conductivity / pore size distribution

^* 联系作者： 王 琛（1989-），男，博士，副教授，博士生导师，主要从事岩土工程与基础工程领域的教学与科研工作。 E-mail：cwang33@ tongji.edu.cn

0 引 言

非饱和土的水力特性通常采用持水曲线（SWRC）和相对渗透系数曲线（RHCC）来表征，其在模拟液-气流动与变形之间的耦合行为、估算非饱和土的抗剪强度以及解释土壤中污染物的扩散特性等方面发挥着重要作用^［1-2］。研究人员提出了一系列SWRC和RHCC模型来描述非饱和土的水力特性，如VG模型^［3］、BC模型^［4］和Fredlund-Xing模型^［5］等，这些模型通常能够较好地描述特定初始孔隙比下土体的持水和渗流特性。然而，当土体发生显著的体积变形时，土体的孔隙尺寸分布（PSD）也会发生变化，进而显著影响非饱和土的水力特性，但这些早期模型无法描述非饱和土SWRC和RHCC随PSD的演变规律。

鉴于这一原因，研究人员在过去几十年中提出了各种改进的水力本构模型，可概括为经验模型和理论模型。前者通常通过经验函数将进气值与初始干密度（或孔隙比、体积应变）关联起来^［5-10］。如Huang等^［8］在BC模型中引入一个单参数的孔隙比函数，提出了改进的SWRC模型，并推导了可反映初始干密度影响的RHCC模型。Gallipoli等^［9］在VG模型引入了幂函数以描述体积变形对饱和度的影响。基于饱和度和体积变形的内在约束，Zhou等^［10］提出了一个与初始孔隙比相关的单参数函数，其可以结合任何既有的持水曲线模型来量化初始干密度对SWRC的影响。后一种改进模型则基于特定物理机制，如Simms等^［11］基于网络方法研究了体积变化对非饱和土水力特性的影响。Hu等^［12］假设变形状态下的PSD可以通过对参考状态下的PSD进行缩放和平移来获得，并推导出改进的水力模型。然而，由于采用了经验PSD模型或SWRC模型，这些模型的参数没有明确的物理意义，导致适应性和扩展性有限。

分形理论因其具有明确的物理背景已被广泛用于模拟多孔介质的PSD、孔隙粗糙度、迂曲度和颗粒尺寸分布，进而研究多孔介质的水力特性。如，Tyler等^［13］假设颗粒尺寸分布服从分形分布并提出了相应的SWRC模型。Yu等^［14］采用分形理论提出了一个多相流的渗透率模型。Russell等^［15］采用遵循分形缩放定律的非均匀毛细管簇来研究SWRC的滞后行为。尽管当前研究人员提出了不同的分形水力模型，但很少涉及体积变形对水力特性的影响，因此难以描述非饱和土SWRC和RHCC随PSD的演变规律。

本文考虑了土体体积变形引起的PSD演变，基于面积分形维数和PSD之间的内在关系，通过引入体积变形的影响扩展了分形定律。在此基础上，推导了基于分形理论的非饱和土的水力模型，包括SWRC和RHCC模型。在该模型中，所有参数均可通过SWRC数据校准，然后用于预测土壤的RHCC。最后，采用7组不同土的试验数据对提出的水力模型进行了验证。

1 模型建立

1.1　 PSD演变的分形描述

PSD显著影响土体的持水和渗流特性。为描述土体的PSD特征，研究人员提出了多种分布函数，如对数正态分布^［16］、高斯分布^［17］以及分形分布^［13］。由于其明确的物理背景，分形理论在多孔介质的传热传质中得到了广泛应用^［14］。试验结果表明 ^［18-19］，体积变形会改变土体的PSD，进而改变面积分形维数。因此，下面将基于试验数据，分析土体在压缩过程中PSD演变规律，然后通过引入体积变形的影响来扩展当前的分形理论。

假设高度为L₀的代表性体积单元（REV），其中孔隙尺寸分布从r_max（孔隙的最大半径）到r_min（孔隙的最小半径），并遵循分形定律。为简化起见，假设REV内的孔隙系统为毛细管簇。根据分形理论^［14］，孔隙半径大于r的累积孔隙数量N可以表示为

N(\epsilon \ge r)=(r/r_{\mathrm{m}\mathrm{a}\mathrm{x}}{)}^{-D_{\mathrm{f}}}（1）

式中：r为孔隙半径；ε为测量标度；r_max为最大孔隙半径；D_f为面积分形维数（1≤D_f≤2）。由式（1）可知，REV中的孔隙总数量N_t可表示为

N_{\mathrm{t}}={\left(r_{\mathrm{m}\mathrm{i}\mathrm{n}}/r_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)}^{-D_{\mathrm{f}}}（2）

此外，由式（1）可得半径范围为r至r+dr的孔隙数量，即

-\mathrm{d}N=D_{\mathrm{f}}{r}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{D_{\mathrm{f}}}{r}^{-\left(D_{\mathrm{f}}+\mathrm{1}\right)}\mathrm{d}r（3）

式中：-\mathrm{d}N>\mathrm{0}，负号表示孔隙数量随孔径的增大而减少。

对于PSD服从分形分布的多孔介质，孔隙率ϕ与面积分形维数D_f的关系^［20］可表示为

D_{\mathrm{f}}=\mathrm{2}-\frac{\mathrm{l}\mathrm{n}\{e/[c(\mathrm{1}+e)]\}}{\mathrm{l}\mathrm{n}{r}_{\mathrm{m}\mathrm{i}\mathrm{n}}-\mathrm{l}\mathrm{n}{r}_{\mathrm{m}\mathrm{a}\mathrm{x}}}（4）

式中：ϕ=e/（1+e）为孔隙率，e为孔隙比；c是材料常数。Yu等^［20］表示，对于谢尔宾斯基地毯和谢尔宾斯基篮子等严格自相似的多孔介质，c=1。然而，土体的孔隙系统为统计自相似多孔介质，因此c>1。

试验研究表明^{［11，21］}，当土体压缩时，较大的孔隙更容易被压缩，但较小的孔隙几乎不会被压缩。因此，最小孔径r_min在土体压缩过程中基本不发生变化，故可视为常数。结合式（4）可得参考状态与变形状态下土体最大孔径之间的关系，即

\mathrm{l}\mathrm{n}\frac{r_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\text{'}}}{r_{\mathrm{m}\mathrm{a}\mathrm{x}}}=\frac{\mathrm{l}\mathrm{n}\{e/[c(\mathrm{1}+e)]\}}{\mathrm{2}-D_{\mathrm{f}}}-\frac{\mathrm{l}\mathrm{n}\{e^{\text{'}}/[c(\mathrm{1}+e^{\text{'}})]\}}{\mathrm{2}-D_{\mathrm{f}}^{\text{'}}}（5）

式中：D_{\mathrm{f}}^{\text{'}}、e^{\text{'}}和r_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\text{'}}分别为变形状态下的分形维数、孔隙比和最大孔径。

Khoshghalb等^［22］指出分形多孔介质的最大孔隙半径与其孔隙比e的指数函数之间存在线性关系，即

r_{\mathrm{m}\mathrm{a}\mathrm{x}}\propto e^{\eta }（6）

式中：η为材料参数，其与孔隙尺寸分布的特征相关。若忽略土体在压缩过程中的颗粒破碎，η可视为常数。由式（6）可得

\frac{r_{\mathrm{m}\mathrm{a}\mathrm{x}}}{r_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\text{'}}}={\left(\frac{e}{e^{\text{'}}}\right)}^{\eta }（7）

结合式（5）和式（7）可得面积分形维数与孔隙比之间的关系，即

D_{\mathrm{f}}^{\text{'}}=\mathrm{2}-\frac{(\mathrm{2}-D_{\mathrm{f}})\mathrm{l}\mathrm{n}\{e^{\text{'}}/[c(\mathrm{1}+e^{\text{'}})]\}}{\mathrm{l}\mathrm{n}\{e/[c(\mathrm{1}+e)]\}+\eta (\mathrm{2}-D_{\mathrm{f}})\mathrm{l}\mathrm{n}(e/e^{\text{'}})}（8）

由式（2）积分可得半径大于r孔隙的累积体积V，即

\begin{array}{l}V=V_{\mathrm{T}}\frac{{\int }_{r_{\mathrm{m}\mathrm{i}\mathrm{n}}}^r\mathrm{\pi }{r}^{\mathrm{2}}{L}_{\mathrm{0}}\mathrm{d}N(r)}{{\int }_{r_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{r_{\mathrm{m}\mathrm{a}\mathrm{x}}}\mathrm{\pi }{r}^{\mathrm{2}}{L}_{\mathrm{0}}\mathrm{d}N(r)}\\ =V_{\mathrm{T}}\frac{r^{\mathrm{2}-D_{\mathrm{f}}}-{r_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{\mathrm{2}-D_{\mathrm{f}}}}{{r_{\mathrm{m}\mathrm{a}\mathrm{x}}}^{\mathrm{2}-D_{\mathrm{f}}}-{r_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{\mathrm{2}-D_{\mathrm{f}}}}\end{array}（9）

式中：V_T是孔隙总体积。

因此，式（7）~式（9）可用于定量表征体积变形对土体PSD演变规律的影响。

1.2　 SWRC模型

孔隙半径与基质吸力的关系可采用Young-Laplace方程描述^［5］，即

s=\frac{\tau {\sigma }_{\mathrm{s}}\mathrm{c}\mathrm{o}\mathrm{s}\theta }{r}（10）

式中：τ为几何常数；σ_s为孔隙水表面张力；θ为接触角，假设其为常数。

根据式（10）可知，孔隙越小，基质吸力越大。这表明在土壤润湿过程中，小孔隙优先被水填充，而在干燥过程中，大孔隙优先排水。因此，存在一个对应于特定基质吸力的临界孔隙半径r_c，使得半径小于r_c的孔隙被水填充，而其他孔隙则被空气填充。对于给定的临界孔隙半径r_c，由式（2）可得土体的饱和度S_r为

\begin{array}{l}{S}_{\mathrm{r}}=\frac{{\int }_{r_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{r_{\mathrm{c}}}\mathrm{\pi }{r}^{\mathrm{2}}{L}_{\mathrm{0}}\mathrm{d}N(r)}{{\int }_{r_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{r_{\mathrm{m}\mathrm{a}\mathrm{x}}}\mathrm{\pi }{r}^{\mathrm{2}}{L}_{\mathrm{0}}\mathrm{d}N(r)}\\ =\frac{(r_{\mathrm{c}}/r_{\mathrm{m}\mathrm{a}\mathrm{x}}{)}^{\mathrm{2}-D_{\mathrm{f}}}-(r_{\mathrm{m}\mathrm{i}\mathrm{n}}/r_{\mathrm{m}\mathrm{a}\mathrm{x}}{)}^{\mathrm{2}-D_{\mathrm{f}}}}{\mathrm{1}-(r_{\mathrm{m}\mathrm{i}\mathrm{n}}/r_{\mathrm{m}\mathrm{a}\mathrm{x}}{)}^{\mathrm{2}-D_{\mathrm{f}}}}\end{array}（11）

将式（3）和式（10）代入式（11）可得

S_{\mathrm{r}}=\frac{c(s_{\mathrm{a}}{/s)}^{\mathrm{2}-D_{\mathrm{f}}}-e/(\mathrm{1}+e)}{c-e/(\mathrm{1}+e)}（12）

式中：s（=κT_scosθ/r_c）为基质吸力；s_a（=κT_scosθ/r_max）为进气值。式（12）即为给定初始孔隙比下的持水曲线方程。

然而，当土体变形时，结合式（7）、式（8）、式（10）和式（12）可得土体在任意初始孔隙比的持水曲线方程，即

S_{\mathrm{r}}=\frac{c[\alpha /(s{e^{\text{'}}}^{\eta }{)]}^{\mathrm{2}-D_{\mathrm{f}}^{\text{'}}}-e^{\text{'}}/(\mathrm{1}+e^{\text{'}})}{c-e^{\text{'}}/(\mathrm{1}+e^{\text{'}})}（13）

式中：\alpha =s_{\mathrm{a}}{e}^{\eta }，其与土体在参考孔隙比下的进气值相关；D_{\mathrm{f}}^{\text{'}}为变形状态下土体的分形维数，可由式（8）求得；e'和e分别为变形状态和参考状态下土体的孔隙比。

1.3　 RHCC模型

Bousfield等^［23］提出单个孔隙中的体积流率可由Hagen-Poiseuille方程描述，即

q_{\mathrm{p}}(r)=\frac{\mathrm{\pi }{r}^{\mathrm{4}}}{\mathrm{8}\mu }\frac{\mathrm{\Delta }P}{L_{\mathrm{t}}}（14）

式中：q_p为单个孔隙的体积流率；r为孔隙半径；μ为孔隙水的黏度；\mathrm{\Delta }P为水力梯度；L_t为流径的迂曲长度。

Brooks等^［4］提出流径长度与饱和度的关系可表示为

L_{\mathrm{t},\mathrm{s}}={S_{\mathrm{r}}}^{\mathrm{2}}{L}_{\mathrm{t},\mathrm{u}}（15）

式中：L_t,_s和L_t,_u分别为饱和与非饱和状态下的流径长度。

假设土体在饱和状态下流径长度与孔隙长度相等，则饱和状态下的REV的总体积流率Q_{\mathrm{t}}可由式（2）和式（15）得到，即

\begin{array}{l}{Q}_{\mathrm{t}}={\int }_{r_{min}}^{r_{max}}{q}_p(r)dN(r)\\ =\frac{\mathrm{\pi }\mathrm{\Delta }P}{\mathrm{8}\mu L_{\mathrm{0}}}\frac{D_{\mathrm{f}}{r}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{D_{\mathrm{f}}}}{\mathrm{4}-D_{\mathrm{f}}}(r_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{4}-D_{\mathrm{f}}}-r_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{4}-D_{\mathrm{f}}})\end{array}（16）

另一方面，由达西定律可得饱和状态下的REV的总体积流率Q_{\mathrm{t}}为

Q_{\mathrm{t}}=K_{\mathrm{s}}\frac{\mathrm{\Delta }P}{L_{\mathrm{0}}}A=K_{\mathrm{s}}\frac{\mathrm{\Delta }P}{L_{\mathrm{0}}}\frac{A_{\mathrm{p}}}{{\varphi }_{\mathrm{R}}}（17）

式中：K_s是饱和渗透系数；A为REV的横截面积；A_p和ϕ_R分别为REV横截面上的总孔隙体积和孔隙率。

结合式（16）和式（17）可得饱和渗透系数，即

K_{\mathrm{s}}=\frac{\pi \varphi }{\mathrm{8}\mu A_{\mathrm{p}}}\frac{D_{\mathrm{f}}{r}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{4}}}{\mathrm{4}-D_{\mathrm{f}}}[\mathrm{1}-{(\frac{r_{\mathrm{m}\mathrm{i}\mathrm{n}}}{r_{\mathrm{m}\mathrm{a}\mathrm{x}}})}^{\mathrm{4}-D_{\mathrm{f}}}]（18）

式中：2≤4-D_f≤3试验结果表明^{［11，21］}，对于土体而言，r_{\mathrm{m}\mathrm{i}\mathrm{n}}/r_{\mathrm{m}\mathrm{a}\mathrm{x}}\ll {\mathrm{10}}^{-\mathrm{2}}。因此，(r_{\mathrm{m}\mathrm{i}\mathrm{n}}/r_{\mathrm{m}\mathrm{a}\mathrm{x}}{)}^{\mathrm{4}-D_{\mathrm{f}}}可近似视为0，故式（18）可简化为

K_{\mathrm{s}}=\frac{\mathrm{\pi }\varphi }{\mathrm{8}\mu A_{\mathrm{p}}}\frac{D_{\mathrm{f}}}{\mathrm{4}-D_{\mathrm{f}}}{r}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{4}}（19a）

当土体发生体积变形时，饱和渗透系数K_{\mathrm{s}}^{\text{'}}可表示为

K_{\mathrm{s}}^{\text{'}}=\frac{\pi {\varphi }^{\text{'}}}{\mathrm{8}\mu A_{\mathrm{p}}^{\text{'}}}\frac{D_{\mathrm{f}}^{\text{'}}}{\mathrm{4}-D_{\mathrm{f}}^{\text{'}}}{r}^{\text{'}}{}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{4}}（19b）

式中：A_{\mathrm{p}}^{\text{'}}为发生体积变形后REV横截面的孔隙总面积。

结合式（19a）和式（19b）可得

\frac{K_{\mathrm{s}}^{\text{'}}}{K_{\mathrm{s}}}=\frac{A_{\mathrm{p}}}{A_{\mathrm{p}}^{\text{'}}}\frac{{\varphi }^{\text{'}}}{\varphi }\frac{D_{\mathrm{f}}^{\text{'}}(\mathrm{4}-D_{\mathrm{f}})}{D_{\mathrm{f}}(\mathrm{4}-D_{\mathrm{f}}^{\text{'}})}{(\frac{r_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\text{'}}}{r_{\mathrm{m}\mathrm{a}\mathrm{x}}})}^{\mathrm{4}}（20）

REV横截面的总孔隙体积可由式（2）求得，即

\begin{array}{l}{A}_{\mathrm{p}}={\int }_{r_{\mathrm{m}\mathrm{i}\mathrm{n}}}^{r_{\mathrm{m}\mathrm{a}\mathrm{x}}}\mathrm{\pi }{r}^{\mathrm{2}}\mathrm{d}N(r)\\ =\frac{\pi D_{\mathrm{f}}{r}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{D_{\mathrm{f}}}}{\mathrm{2}-D_{\mathrm{f}}}(r_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{2}-D_{\mathrm{f}}}-r_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{2}-D_{\mathrm{f}}})\end{array}（21）

因此，A_{\mathrm{p}}/A_{\mathrm{p}}^{\text{'}}可表示为

\frac{A_{\mathrm{p}}}{A_{\mathrm{p}}^{\text{'}}}=\frac{D_{\mathrm{f}}}{D_{\mathrm{f}}^{\text{'}}}\frac{\mathrm{2}-D_{\mathrm{f}}^{\text{'}}}{\mathrm{2}-D_{\mathrm{f}}}\frac{\mathrm{1}-(r_{\mathrm{m}\mathrm{i}\mathrm{n}}/r_{\mathrm{m}\mathrm{a}\mathrm{x}}{)}^{\mathrm{2}-D_{\mathrm{f}}}}{\mathrm{1}-(r_{\mathrm{m}\mathrm{i}\mathrm{n}}/r_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\text{'}}{)}^{\mathrm{2}-D_{\mathrm{f}}^{\text{'}}}}\frac{r_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{2}}}{r^{\text{'}}{}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{2}}}（22）

结合式（11）、式（20）和式（22）可得土体在变形后的饱和渗透系数为

K_{\mathrm{s}}^{\text{'}}=\frac{\mathrm{2}-D_{\mathrm{f}}^{\text{'}}}{\mathrm{2}-D_{\mathrm{f}}}\frac{(\mathrm{4}-D_{\mathrm{f}})}{(\mathrm{4}-D_{\mathrm{f}}^{\text{'}})}\frac{{\varphi }^{\text{'}}}{\varphi }\frac{c-\varphi }{c-{\varphi }^{\text{'}}}{(\frac{e^{\text{'}}}{e})}^{\mathrm{2}\eta }{K}_{\mathrm{s}}（23）

当土体处于非饱和状态时，REV的体积流率可由式（2）和式（15）确定，即

\begin{array}{l}{Q}_{\mathrm{u},\mathrm{t}}={\int }_{r_{min}}^{r_c}{q}_p(r)dN(r)\\ =\frac{\pi \mathrm{\Delta }P}{\mathrm{8}\mu L_{\mathrm{t},\mathrm{u}}}\frac{D_{\mathrm{f}}{r}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{D_{\mathrm{f}}}}{\mathrm{4}-D_{\mathrm{f}}}(r_{\mathrm{c}}^{\mathrm{4}-D_{\mathrm{f}}}-r_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{4}-D_{\mathrm{f}}})\end{array}（24）

因此，结合式（15）、式（16）和式（24）可得土体在非饱和状态下的相对渗透系数，即

\begin{array}{l}{k}_{\mathrm{r}}=\frac{Q_{\mathrm{u},\mathrm{t}}}{Q_{\mathrm{t}}}\\ =\frac{L_{\mathrm{0}}}{L_{\mathrm{t},\mathrm{u}}}\frac{(r_{\mathrm{c}}/r_{\mathrm{m}\mathrm{a}\mathrm{x}}{)}^{\mathrm{4}-D_{\mathrm{f}}}-(r_{\mathrm{m}\mathrm{i}\mathrm{n}}/r_{\mathrm{m}\mathrm{a}\mathrm{x}}{)}^{\mathrm{4}-D_{\mathrm{f}}}}{\mathrm{1}-(r_{\mathrm{m}\mathrm{i}\mathrm{n}}/r_{\mathrm{m}\mathrm{a}\mathrm{x}}{)}^{\mathrm{4}-D_{\mathrm{f}}}}\end{array}（25）

结合式（11）与式（15）可得

k_{\mathrm{r}}={S_{\mathrm{r}}}^{\mathrm{2}}{\left\{S_{\mathrm{r}}\left[\mathrm{1}-\frac{e}{c(\mathrm{1}+e)}\right]+\frac{e}{c(\mathrm{1}+e)}\right\}}^{\frac{\mathrm{4}-D_{\mathrm{f}}}{\mathrm{2}-D_{\mathrm{f}}}}（26）

因此，结合式（26）和式（8）可得土体在变形后的RHCC，即

k_{\mathrm{r}}^{\text{'}}={S_{\mathrm{r}}}^{\mathrm{2}}{\left\{S_{\mathrm{r}}\left[\mathrm{1}-\frac{e^{\text{'}}}{c(\mathrm{1}+e^{\text{'}})}\right]+\frac{e^{\text{'}}}{c(\mathrm{1}+e^{\text{'}})}\right\}}^{\frac{\mathrm{4}-D_{\mathrm{f}}^{\text{'}}}{\mathrm{2}-D_{\mathrm{f}}^{\text{'}}}}（27）

式中：D_{\mathrm{f}}^{\text{'}}可由式（8）确定。

由式（27）可得土体在变形后的相对渗透系数，即

K_{\mathrm{u}}^{\text{'}}=k_{\mathrm{r}}^{\text{'}}{K}_{\mathrm{s}}^{\text{'}}（28）

式（28）即为提出的考虑PSD演化的RHCC模型，其包含四个参数（α、η、D_f和c）。α和D_f分别为参数状态下土体的进气值和面积分形维数；η为与孔隙尺寸分布相关的材料常数；c为与孔隙迂曲度和异质性相关的材料常数。这些模型参数均可由参考状态下土体的SWRC试验数据进行校准，然后用于预测土体在变形后的SWRC和RHCC。

2 模型验证与讨论

2.1　初始孔隙比对PSD的影响

图1展示了原状黏土^［24］和黏土（soil 1）^［25］的累积孔隙体积在压缩过程中的演变规律，其中实线为式（9）的预测曲线。由该图可知，随着土体体积发生变化，PSD也随着改变，其可采用提出的扩展分形理论来描述。土体的分形维数随着孔隙比减小而减小。如图1所示，提出的扩展分形理论可较好地描述PSD随体积变形的演变规律。

2.2　常吸力下孔隙比对饱和度的影响

在恒定吸力下，饱和度会随体积变形而变化，这可通过所提出的分形SWRC模型来描述。如图2、表1所示，采用了四种不同土的试验数据对提出的模型进行验证。对于Maryland黏土^［26］，采用s= 900 kPa条件下的S_r-e数据用来校准模型参数。参考状态下土壤的分形维数D_f（e=0.7）为1.78，参数c、η和α分别确定为5、4.97和20 kPa。然后，使用这些参数预测其他恒定吸力下的S_r-e曲线。如图2（a）所示，预测曲线（s=200 kPa和100 kPa）与试验数据吻合良好。其他三种土也表现出类似的结果，如图2（b）—图2（d）所示。

2.3　初始孔隙比对RHCC的影响

Laliberte等^［30］对Touchet粉质壤土和Columbia砂质壤土的两种土进行了水力特性测试。Touchet粉质壤土是由15%黏土、53%粉土和35%砂土组成的混合物，而Columbia砂质壤土是由11%黏土、53%粉土和54%砂土组成的混合物。

图3和图4展示了Touchet粉质壤土和Columbia砂质壤土在不同孔隙比下的SWRC和RHCC特征。因此，可以采用从SWRC数据校准的参数来评估提出的水力模型预测RHCC的能力。对于Touchet粉质壤土，使用孔隙比为0.97和0.86的SWRC数据来校准模型参数，如图3（a）所示。参考状态下（e=0.97）的土体的面积分形维数为1.30，其他模型参数列在了表1中。然后，使用这些参数来预测其他孔隙比下土体的SWRC和RHCC。如图3（a）和图3（b）所示，预测曲线与试验数据吻合良好。此外，对于给定基质吸力下的Touchet粉质壤土，相对渗透系数k_r随孔隙比的减小而增加，而饱和渗透系数在土体压缩过程中减小。

由于砂土的孔隙尺寸通常大于粉土，在相同孔隙比下，Columbia砂质壤土（图4）的吸力范围和进气值小于Touchet粉质壤土（图3）。对于参考状态下的Columbia砂质壤土（e=1.19），面积分形维数为1.10。图4（b）展示了使用从SWRC数据获得的参数（表1）预测的曲线，结果表明预测曲线与试验数据吻合良好。

图5所示为Leij等^［31］对Berlin粗砂进行的水力特性测试结果，包含SWRC和RHCC的完整水力特性数据。参考状态下（e=0.54）的分形维数为1.02。图5比较了不同压缩状态下的预测曲线和实验数据。结果再次表明，所提出的模型能够较好地描述孔隙比对水力特性的影响。

上述示例表明，采用扩展的分形定律反映土体压缩过程中PSD的演变，并基于此建立的非饱和土水力模型，可有效描述PSD演变对非饱和SWRC和RHCC的影响。

2.4　 PSD与孔隙结构异质性对RHCC的影响

图6展示了PSD对非饱和土相对渗透系数的影响。不同结构的土在给定参考孔隙比下的D_f值不同。如图6（a）所示，土A和土B是两种不同结构的土，参考状态（e=1.0）下的面积分形维数D_f分别为1.5和1.7。较大的D_f意味着小孔隙的体积比例较大，导致流体阻力较大。因此，如图6（b）所示，在相同饱和度或基质吸力下，土A的K_u大于土B。

另一方面，体积变形会改变土体的PSD，D_f也随之变化，如图6（a）所示。对于土A，随着孔隙比的减小，面积分形维数也减小，且K_u在土体压缩过程中也会随之减小。此外，从图6（a）可观察到，对于给定的土（如土A），分形维数与孔隙比曲线（D_f-e曲线）的斜率随孔隙比的增加而减小。这是因为当孔隙比较大时，土壤压缩会整体压缩所有孔隙，因此最大和最小孔隙尺寸的比例随着孔隙比的减小而缓慢变化。然而，随着孔隙比的减小，较小的孔隙变得越来越难被压缩，压缩的孔隙空间主要由较大的孔隙承担。这导致土体中最大和最小孔隙尺寸的比例迅速变化。因此，D_f-e曲线的斜率随着孔隙比的增加而逐渐减小。

图7展示了孔隙迂曲度和异质性对非饱和土相对渗透系数的影响。c值越大，孔隙的迂曲度和异质性越大，从而导致流体阻力越大。因此，如图7（a）所示，c值较大的土壤其K_u较小。在给定的参考孔隙比e=1.0下，土A（c=20）和土B（c=7）的分形维数D_f相等（D_f=1.7），这意味着两种土在e=1.0时具有相同的孔隙尺寸分布。采用参考孔隙比是为了方便比较两种土在不同c值下D_f的变化模式。若土A和土B的c值也相等，则两种土的D_f随孔隙比的变化应相同。如图7（b）所示，尽管c值有所不同，但两种土的D_f-e曲线差异很小。这是因为参数c主要影响孔隙迂曲度，其对主要受孔隙尺寸分布影响的D_f影响较小。另一方面，从图7（b）可以看出，土壤孔隙的迂曲度越小，孔隙尺寸分布随孔隙比的变化越大。值得注意的是，可以通过侵入-挤出循环汞压入法定量表征孔隙空间的迂曲度和非均匀性。然而，若能够直接以三维形式呈现孔隙的空间结构，将更有利于研究实际土壤的水力特性。

3 结 论

本文基于孔隙尺寸分布（PSD）随土体孔隙比减小而演变的试验证据，通过引入体积变形的影响，提出了扩展的分形缩放定律。然后，采用面积分形维数反映PSD的演变，导出了考虑体积变形影响的非饱和土SWRC模型。在此基础上，结合Hagen-Poiseuille方程，在未引入任何新参数的情况下推导得到了可反映PSD演变影响的RHCC模型。提出的非饱和土水力模型中包含四个具有明确物理意义的参数，其均可通过SWRC数据校准，然后用于预测土体的RHCC。采用了7组从砂土到黏土的试验数据来评估提出的水力模型。结果表明，该模型能够较好地描述孔隙比变化对SWRC和RHCC的影响。

此外，模型参数分析结果表明，随着孔隙比的减小，PSD随孔隙比的演变会使得分形维数减小，进而导致非饱和土的相对渗透系数随之减小。而孔隙的迂曲度和异质性越大则会导致流体阻力增大，进而降低土体的持水和渗流特性，导致土体的相对渗透系数减小。

参考文献

[1] 谢定义.非饱和土土力学［M］.北京：高等教育出版社，2015.

[2] XIE Dingyi.Soil mechanics for unsaturated soils［M］.Beijing：Higher Education Press，2015.（in Chinese）

[3] 师晨翔，杨永平，马际首，考虑土体基质吸力刚性挡墙基坑抗倾覆稳定性分析［J］.结构工程师，2018，34（3）：152-157.

[4] SHI Chenxiang，YANG Yongping，MA Jishou，et al.Anti-overturning stability analysis of rigid wall considering soil matric suction［J］.Structural Engineers，2018，34（3）：152-157.（in Chinese）

[5] 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.

[6] BROOKS R H，COREY A T.Hydraulic properties of porous media［C］//Hydrology Paper No.3.Fort Collins，CO，USA：Colorado State University，1964.

[7] FREDLUND D G，XING A.Equations for the soil-water characteristic curve［J］.Canadian Geotechnical Journal，1994，31（4），521-532.

[8] 周葆春，孔令伟，梁维云，压缩过程中非饱和膨胀土体变特征与持水特性的水力耦合效应［J］.岩土工程学报，2015，37（4）：629-640.

[9] ZHOU Baochun，KONG Lingwei，LAING Weiyun，et al.Hydro-mechanical coupling effects on volume change and water retention behaviour of unsaturated expansive soils during compression［J］.Chinese Journal of Geotechnical Engineering，2015，37（4）：629-640.（in Chinese）

[10] 林志强，钱建固，时振昊.毛细-吸附作用下考虑孔隙比影响的单/双峰土体持水曲线模型［J］.岩土力学，2021，42（9）：2499-2506.

[11] LIN Zhiqiang，Qian Jiangu，SHI Zhenhao.Modeling unimodal/bimodal soil-water retention curves considering the influence of void ratio under capillarity and adsorption［J］.Rock and Soil Mechanics，2021，42（9）：2499-2506.（in Chinese）

[12] HUANG S，BARBOUR S L，FREDLUND D G.Development and verification of a coefficient of permeability function for a deformable unsaturated soil［J］.Canadian Geotechnical Journal，1998，35（3）：411-425.

[13] GALLIPOLI D，BRUNO A W，D'ONZA F，et al.A bounding surface hysteretic water retention model for deformable soils［J］.Géotechnique，2015，65（10）：793-804.

[14] ZHOU A N，SHENG D，CARTER J P.Modelling the effect of initial density on soil-water characteristic curves［J］.Géotechnique，2012，62（8）：669-680.

[15] SIMMS P H，YANFUL E K.A pore-network model for hydromechanical coupling in unsaturated compacted clayey soils［J］.Canadian Geotechnical Journal，2005，42（2）：499-514.

[16] HU R，CHEN Y F，LIU H H，et al.A water retention curve and unsaturated hydraulic conductivity model for deformable soils：consideration of the change in pore-size distribution［J］.Géotechnique，2013，63（16）：1389-1405.

[17] TYLER S W，WHEATCRAFT S W.Application of fractal mathematics to soil water retention estimation［J］.Soil Science Society of America Journal，1989，53（4）：987-996.

[18] YU B，LI J，LI Z，et al.Permeabilities of unsaturated fractal porous media［J］.International Journal of Multiphase Flow，2003，29（10）：1625-1642.

[19] RUSSELL A R.How water retention in fractal soils depends on particle and pore sizes，shapes，volumes and surface areas［J］.Géotechnique，2014，64（5）：379-390.

[20] KOSUGI K I.Three-parameter lognormal distribution model for soil water retention［J］.Water Resources Research，1994，30（4）：891-901.

[21] XU C，TORRES-VERDÍN C.Pore system characterization and petrophysical rock classification using a bimodal Gaussian density function［J］.Mathematical Geosciences，2013，45（6）：753-771.

[22] LAPIERRE C，LEROUEIL S，LOCAT J.Mercury intrusion and permeability of Louiseville clay［J］.Canadian Geotechnical Journal，1990，27（6）：761-773.

[23] NINJGARAV E，CHUNG S G，JANG W Y，et al.Pore size distribution of Pusan clay measured by mercury intrusion porosimetry［J］.KSCE Journal of Civil Engineering，2007，11（3）：133-139.

[24] YU B，LI J.Some fractal characters of porous media［J］.Fractals，2001，9（3），365-372.

[25] TANAKA H，SHIWAKOTI D R，OMUKAI N，et al.Pore size distribution of clayey soils measured by mercury intrusion porosimetry and its relation to hydraulic conductivity［J］.Soils and Foundations，2003，43（6）：63-73.

[26] KHOSHGHALB A，PASHA A Y，KHALILI N.A fractal model for volume change dependency of the water retention curve［J］.Géotechnique，2015，65（2）：141-146.

[27] BOUSFIELD D W，KARLES G.Penetration into three-dimensional complex porous structures［J］.Journal of Colloid and Interface Science，2004，270（2）：396-405.

[28] LAPIERRE C，LEROUEIL S，LOCAT J.Mercury intrusion and permeability of Louiseville clay［J］.Canadian Geotechnical Journal，1990，27（6）：761-773.

[29] GRIFFITHS F J，JOSHI R C.Change in pore size distribution due to consolidation of clays［J］.Géotechnique，1989，39（1）：159-167.

[30] BURTON G J，SHENG D，AIREY D.Experimental study on volumetric behaviour of Maryland clay and the role of degree of saturation［J］.Canadian Geotechnical Journal，2014，51（12）：1449-1455.

[31] SIVAKUMAR R.Effects of anisotropy on the behaviour of unsaturated compacted clay：PhD thesis［D］.Belfast，UK：Queen's University Belfast，UK，2005.

[32] D'ONZA F，GALLIPOLI D，WHEELER S，et al.Benchmark of constitutive models for unsaturated soils［J］.Géotechnique，2011，61（4）：283-302.

[33] RAVEENDIRARAJ A.Coupling of mechanical behaviour and water retention behaviour in unsaturated soils：PhD thesis［D］.Glasgow，UK：University of Glasgow，2009.

[34] LALIBERTE G E，COREY A T，BROOKS R.Properties of unsaturated porous media［M］.Fort-Collins，CO，USA：Colorado State University，1966.

[35] LEIJ F J，ALVES W J，VAN GENUCHTEN M，et al.The UNSODA unsaturated soil hydraulic database，version 1.0，Report EPA/600/R-96/095［R］.Cincinnati，OH，USA：National Risk Management Research Laboratory，US Environmental Protection Agency，1996.

基金资助

国家重点研发计划课题(2022YFC3803000)

国家重点研发计划课题(2016YFC080020)