考虑膜态水渗流的非饱和土渗透系数计算方法研究

翟钱; 沈天伦; 田刚; 戴国亮; 赵学亮; 龚维明; 蔡建国

doi:10.11779/CJGE20221322

考虑膜态水渗流的非饱和土渗透系数计算方法研究 English Version

翟钱^{1, 2,},
沈天伦²,
田刚²,
戴国亮²,
赵学亮²,
龚维明²,
蔡建国²

1.
东南大学混凝土及预应力混凝土结构教育部重点实验室, 江苏南京 210096
2.
东南大学土木工程学院, 江苏南京 210096

基金项目:

国家自然科学基金面上项目 52078128

国家自然科学基金面上项目 52378328

国家自然科学基金面上项目 52178317

国家自然科学基金面上项目 52378329

详细信息

作者简介:
翟钱(1982—)，男，博士，副教授，主要从事非饱和土力学相关方面的科研和教学工作。E-mail: 101012332@seu.edu.cn

中图分类号: TU43
计量
- 文章访问数: 299
- HTML全文浏览量: 57
- PDF下载量: 88
出版历程
- 收稿日期: 2022-02-28
- 网络出版日期: 2024-02-05
- 刊出日期: 2024-01-31

Prediction of the coefficient of permeability for unsaturated soil by considering the film flow

1.
Key Laboratory of Concrete and Prestressed Concrete Structure of Ministry of Education, Southeast University, Nanjing 210096, China
2.
School of Civil Engineering, Southeast University, Nanjing 210096, China

摘要

摘要: 当前采用土-水特征曲线计算非饱和土渗透系数的方法主要分为经验模型和统计模型两大类。经验模型的计算结果取决于模型中的经验系数，统计模型基于毛细模型采用统计学相关理论推导得出，因此统计模型计算结果相对可靠。研究发现，统计模型主要研究毛细水迁移规律，并未考虑土体中的结合水。因物理化学作用，部分土颗粒表面会吸附水膜，该水膜可作为液态水分的传输媒介。水膜厚度主要取决于土颗粒表面与土中水的相互作用，水膜厚度决定了土颗粒表面膜态水的迁移速率。首先计算单个颗粒表层吸附水膜的厚度、膜态水迁移速率；其次考虑土体中不同粒径颗粒的随机接触，利用级配曲线计算不同粒径颗粒的连接概率；最终依据单颗粒膜态水的迁移速率并结合颗粒间的随机接触概率，提出非饱和土膜态水的渗透系数预测模型。为验证模型的可靠性，选取其他文献中的试验数据对比模型预测结果，对比结果表明本文预测结果与文献试验结果吻合度较好，比传统统计模型表现出一定的优越性。
- 非饱和土 /
- 毛细水渗流 /
- 膜态水渗流 /
- 随机连接 /
- 预测模型
Abstract: It has become a common practice that the hydraulic conductivity of unsaturated soils is estimated from the soil-water characteristic curve by using either empirical method or statistical method. The estimated results by the empirical method are dependent on the empirical parameters while those by the statistical method are determined based on the theory of statistics. As a result, the results by the statistical method are more reliable than those by the empirical method. It is observed that the hydroscopic water is misused as the capillary water in the conventional statistical method. Because of the adsorptive force, the water film (hydroscopic water) is attached around the soil particle, and it acts as the transferring medium for the water migration in soils. The flowing rate of the water in the adsorbed water film is a function of the thickness of the water film. Both the thickness of water film around a single soil particle and the flowing rate through the water film are firstly computed. Subsequently, the probability of the connecting between particles with different sizes is calculated by using the grain size distribution (GSD) data. Consequently, a new equation is proposed for the estimation of the hydraulic conductivity of unsaturated soils by considering the film flows in soils. The proposed method is verified using the experimental data from literatures. It is indicated that results by the newly proposed method provide better agreement with the experimental data as compared with those by the conventional statistical method.
- unsaturated soil /
- capillary flow /
- film flow /
- random connection /
- statistical method

HTML全文

0. 引言

DDA(discontinuous deformation analysis)即非连续变形分析，是一种基于非连续介质力学的数值分析方法。该方法以块体作为基本对象进行分析，块体边界可表示实际岩体裂隙，将各块体的位移视作基本未知量，通过块体间的接触关系将独立块体组合成一个块体系统，在处理块体大变形、大位移问题上具有一定的优势。因此DDA被广泛地应用在岩体这一块体离散介质的研究当中。众多学者将DDA应用于边坡稳定、岩体锚固、矿山开挖、地震等众多工程问题当中^[1-6]。

原始DDA程序主要考虑块体应力应变场，难以直接应用到岩体渗流问题模拟当中。考虑复杂渗流场-应力应变场耦合最终建立能够真实模拟岩体赋存条件的DDA模型将大幅提高该算法的模拟精度和实用性，成为国内外学者关注的热点^[7]。例如，张国新等^[8]基于能量原理将水压荷载加入到原始DDA程序中，分析了水位变化条件下岩质边坡的倾倒变形机理。虞松等^[9]提出了一种矩阵搜索方法，识别并删除DDA生成的裂隙岩体中的不连通裂隙，以生成更精确的裂隙水通网络，并对地下油库水封条件开挖的工程案例进行模拟，发现考虑流固耦合所获得的裂隙水压相对较小而涌水量则有所增加。王知深^[10]在DDARF基础上推导达西-非达西流裂隙岩体渗流场和应力场控制方程，并建立DDARF平台裂隙岩体渗流应力耦合模型，将其应用在碳酸盐岩水力压裂过程分析中。江巍等^[11]采用DDA对藕塘滑坡进行了失稳模式分析，采用块体在饱和状态下的强度参数来进行降雨工况的模拟，但在分析过程中没能考虑流固耦合的影响。Hu等^[12]利用离散裂隙网络模型模拟页岩的自然裂隙，结合DDA提出复杂水力裂隙网络耦合流体力学模型，并将其应用在页岩气压裂开采设计中。Morgan等^[13]利用有限体积裂隙网络模型模拟裂缝中的可压缩流体，研究了在内部有开口加压的格里菲斯裂缝在无渗漏无补给条件下压强和裂缝开度的变化，通过与半解析解进行对比验证了算法的正确性，并对高黏性压裂液在不渗透介质中的传播进行了数值分析。

上述研究成果对于DDA流固耦合的处理多是基于饱和稳态渗流理论。在实际工程中降雨入渗、库水位变化等工程案例均属于非饱和瞬态渗流问题，而目前运用DDA处理这一类问题通常是将降雨工况下材料强度参数进行折减，用饱和强度参数进行模拟，该方法并不能反映降雨入渗过程。针对这一问题，本文对DDA进行二次开发，推导DDA瞬态整体渗流矩阵，并将其加入到总体平衡方程中，实现DDA非饱和瞬态渗流流固耦合计算，通过土柱入渗、砂槽入渗试验验证该计算方法的正确性，为实际工程提供新的数值模拟手段和借鉴参考。

1. DDA稳态渗流计算

DDA以岩块和裂隙组成的岩体作为研究对象，由于岩块渗透能力较弱，岩体渗流通常只考虑裂隙渗流。采用立方定律来描述岩体单裂隙渗流：

。

$q = \frac{{g{e^3}}}{{12\mu }}J 。$

(1)

式中： $e$ 为岩体裂隙开度；q为岩体裂隙流量；μ为流体的动力黏滞系数；J为光滑平板两端的水力梯度。

在实际工程当中，岩体裂隙具有一定的起伏和粗糙度。裂隙的起伏和粗糙度会使流体在裂隙中的流线变长，因此在相同情况下单位时间内通过粗糙裂隙的渗流量也会小于光滑裂隙。为了使计算结果更接近于真实值，引入粗糙度修正系数C，可得

。

$q = \frac{{g{e^3}}}{{12C\mu }}J 。$

(2)

根据相关的水力学试验，粗糙度修正系数C可由岩体裂隙内表面凹凸粗糙程度Δ与裂隙开度e计算得到^[14]

$C = 1 + {\left( {\frac{\Delta }{e}} \right)^{1.5}} 。$

(3)

为将其写入DDA程序当中，将其改写成类似于DDA平衡方程的系数矩阵和未知数矩阵乘积的形式，即

$q = T\Delta H \text{，}$

(4)

式中，T为岩体裂隙的导流系数。

对于岩体裂隙网络中的某一节点，其总流量为各条单裂隙流量对该节点流量的影响之和。对一个由n个节点组成的岩体模型来说，其整体渗流方程可写为

$\left[ {\begin{array}{*{20}{c}} {{T_{11}}}&{{T_{12}}}& \ldots &{{T_{1n}}} \\ {{T_{21}}}&{{T_{22}}}& \ldots &{{T_{2n}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{T_{n1}}}&{{T_{n2}}}& \ldots &{{T_{nn}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{H_1}} \\ {{H_2}} \\ \vdots \\ {{H_n}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{Q_1}} \\ {{Q_2}} \\ \vdots \\ {{Q_n}} \end{array}} \right\} 。$

(5)

式中：裂隙导流系数矩阵T的对角线元素T_ii为与节点i相连接的所有裂隙的导流系数之和；T_ij为由节点i和节点j组成裂隙的导流系数的相反数；水头矩阵H中的H_i为节点i对应的水头；流量矩阵Q中的Q_i为节点i对应的总流量，由连接该点的所有裂隙流量叠加得到，且以流体流入该点为负流出该点为正。

在渗流的模拟过程当中，为了求解式（5）还需要在迭代前给出边界条件，即流量边界和水头边界。当某一节点的水头为已知量时，在迭代时将该点导流系数矩阵T的对角线元素T_ii设置为1，其余元素T_ij设置为0，以保证水头不变。当某一节点的流量为已知量时，将流量矩阵Q中对应的元素Q_i设为已知量即可^[15]。

2. DDA瞬态渗流计算

2.1 瞬态整体渗流矩阵推导

为了运用DDA进行瞬态渗流场模拟，将用体积含水量表示的二维非饱和Richards方程引入到DDA程序当中：

$\frac{\partial }{{\partial x}}\left( {{D_x}(\theta )\frac{{\partial \theta }}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {{D_y}(\theta )\frac{{\partial \theta }}{{\partial y}}} \right) = \frac{{\partial \theta }}{{\partial t}} 。$

(6)

式中：θ为体积含水量，D_x(θ)，D_y(θ)分别为x方向和y方向上的水力扩散系数。

推导符合DDA形式的瞬态整体渗流矩阵，需要将二维Richards方程进行离散化处理。DDA中的块体边界代表了岩体中的裂隙，具有实际的物理意义，因此直接用DDA中块体划分形成的网格作为基础，将 ${D_x}(\theta )\frac{{\partial \theta }}{{\partial x}}$ 视为整体采用中心差分将式（6）中的第一项偏导数离散化，可得

$\frac{\partial }{\partial x}\left({D}_{x}(\theta )\frac{\partial \theta }{\partial x}\right)=\frac{D\left({\theta }_{i+\frac{1}{2}hx}^{j+1}\right)\frac{\partial {\theta }_{i+\frac{1}{2}hx}^{j+1}}{\partial x}-D\left({\theta }_{i-\frac{1}{2}hx}^{j+1}\right)\frac{\partial {\theta }_{i-\frac{1}{2}hx}^{j+1}}{\partial x}}{\Delta x}。$

(7)

式中：θ的上标j+1为时间差分，表示在第j+1时间步对应的体积含水量，下标 $i + \frac{1}{2}hx$ 为空间差分，表示对第i个节点在x方向上取1/2个空间步长对应的体积含水量；Δx为x方向上两个节点的间距。

对式（7）中的偏微分式继续离散化处理，可得

$\left.\begin{array}{l}D\left({\theta }_{i+\frac{1}{2}hx}^{j+1}\right)\frac{\partial {\theta }_{i+\frac{1}{2}hx}^{j+1}}{\partial x}=D\left({\theta }_{i+\frac{1}{2}hx}^{j+1}\right)\frac{{\theta }_{i+hx}^{j+1}-{\theta }_{i}^{j+1}}{\Delta x}\text{ }\text{，}\\ D\left({\theta }_{i-\frac{1}{2}hx}^{j+1}\right)\frac{\partial {\theta }_{i-\frac{1}{2}hx}^{j+1}}{\partial x}=D\left({\theta }_{i-\frac{1}{2}hx}^{j+1}\right)\frac{{\theta }_{i}^{j+1}-{\theta }_{i-hx}^{j+1}}{\Delta x}\text{ }。\end{array}\right\}$

(8)

对于y方向，也有类似于式（8）的计算结果，只需要将空间差分的下标和节点间距分别记为y方向上的对应值即可。

对式（6）的等式右侧同样进行离散化处理，有

$\frac{{\partial \theta }}{{\partial t}} = \frac{{\theta _i^{j + 1} - \theta _i^j}}{{\Delta t}} 。$

(9)

由于差分在渗流模拟过程中具有一定的滞后性，因此对于水力扩散系数D(θ)可以近似的用第j时步的计算结果去表达第j+1时步^[16]。为了简化算式表达进行如下处理：

$\left.\begin{array}{l}\alpha {x}_{i}^{j}=\frac{D\left({\theta }_{i+\frac{1}{2}hx}^{j}\right)}{\Delta {x}^{2}}\text{ }\text{，}\\ \beta {x}_{i}^{j}=\frac{D\left({\theta }_{i-\frac{1}{2}hx}^{j}\right)}{\Delta {x}^{2}}\text{ }\text{，}\\ \alpha {y}_{i}^{j}=\frac{D\left({\theta }_{i+\frac{1}{2}hy}^{j}\right)}{\Delta {y}^{2}}\text{ }\text{，}\\ \beta {y}_{i}^{j}=\frac{D\left({\theta }_{i-\frac{1}{2}hy}^{j}\right)}{\Delta {y}^{2}}\text{ }。\end{array}\right\}$

(10)

对式（7）~（10）进行整理可以得对某一节点体积含水量离散形式的表达式：

$\begin{aligned} \frac{{\theta _i^{j + 1} - \theta _i^j}}{{\Delta t}} =& \alpha x_i^j\theta _{i + hx}^{j + 1} + \beta x_i^j\theta _{i - hx}^{j + 1} + \alpha y_i^j\theta _{i + hy}^{j + 1} +\\ &\beta y_i^j\theta _{i - hy}^{j + 1} - (\alpha x_i^j + \beta x_i^j + \alpha y_i^j + \beta y_i^j)\theta _i^{j + 1} 。 \end{aligned}$

(11)

式中：Δt为瞬态渗流的时间步长，上标为j表示本时步的体积含水量，均为已知量；j+1表示下一时步的体积含水量，通过本时步的状态求解；下标i+hx表示在x方向上与第i个节点间隔hx的下一个节点对应的体积含水量。

由此可知式（11）中每节点在x方向和y方向上均有两个相邻节点连接，对应了采用较为规则的四边形网格划分的情况。在DDA中块体可以是任意形状，因此网格的划分也会随之变得不规则，对网格中某一点i来说，其体积含水量仍然受到与之相连的其他节点的影响，将式（11）整理成矩阵形式，可得

$\left[ {\begin{array}{*{20}{c}} {{A_{11}}}&{{A_{12}}}& \ldots &{{A_{1n}}} \\ {{A_{21}}}&{{A_{22}}}& \ldots &{{A_{2n}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{A_{n1}}}&{{A_{n2}}}& \ldots &{{A_{nn}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {\theta _1^{j + 1}} \\ {\theta _2^{j + 1}} \\ \vdots \\ {\theta _n^{j + 1}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{B_1}} \\ {{B_2}} \\ \vdots \\ {{B_n}} \end{array}} \right\} 。$

(12)

与式（5）相似，式（12）的系数矩阵A是一个只在相邻节点连接的位置上有对应数值的稀疏矩阵。其中，非对角元素A_ij表示由i，j两点构成裂隙对应的系数，即式（10）中的4个代数式之一；对角元素A_ii表示在第i行中所有非对角元素A_ij的和再加上时间步长 $\Delta t$ 的倒数；而等式右侧的常数矩阵B为每一个节点本时步的体积含水量θ与时间步长 $\Delta t$ 的比值。由此可知，对于系数矩阵A由于有时间项的存在，其每一行的对角线元素都大于该行非对角线元素的总和，是具有严格对角占优特点的矩阵。为了保证其收敛性，采用Gauss-Seidel迭代法求解。

对式（12）进行求解需要已知初始条件下每个节点的体积含水量，以及计算模型的边界条件。对DDA模型先进行稳态渗流计算，以稳态渗流的结果作为瞬态渗流的初始条件。对体积含水量已知的边界，将对角线元素A_ii设为1，非对角线元素A_ij设为0即可。

得到节点体积含水量后，通过下式转化成有效饱和度：

${S_{\text{e}}} = \frac{{\theta - {\theta _{\text{r}}}}}{{{\theta _{\text{s}}} - {\theta _{\text{r}}}}} 。$

(13)

式中：S_e为有效饱和度；θ_s为饱和体积含水量；θ_r为残余体积含水量。

根据Van Genuchten^[17]模型和Mualem模型^[18]来表示有效饱和度和相对渗透系数、裂隙压力间的关系：

${P_{\text{c}}} = \frac{{{{({S_{\text{e}}}^{ - \frac{1}{m}} - 1)}^{\frac{1}{n}}}}}{\alpha } \text{，}$

(14)

${T_{\text{r}}} = {S_{\text{e}}}^{\frac{1}{2}}{\left[ {1 - {{\left( {1 - {S_{\text{e}}}^{\frac{1}{m}}} \right)}^m}} \right]^2} 。$

(15)

式中：m，n，α均为土水特征曲线经验拟合参数，m可由 $m = 1 - \frac{1}{n}$ 计算得到；P_c为压力水头；T_r为相对渗透系数，加入到导流系数矩阵T中。

2.2 流固耦合计算

将渗流场计算得到的水头矩阵块体外荷载加入到DDA总体平衡方程中。假设有一块体i，该块体由点1和点2组成的边界受到水压力的作用。将水压力视作在两点之间呈梯形分布的荷载，如图 1所示。

图 1 块体水压力

Figure 1. Water pressures on block

下载: 全尺寸图片幻灯片

点1上的水压力h分别向x，y方向分解，并规定与坐标轴同向为正，与坐标轴反向为负。

$\left.\begin{array}{l}{f}_{x}{}_{1}=-{h}_{1}\mathrm{sin}\alpha \text{ }\text{，}\\ {f}_{y}{}_{1}={h}_{1}\mathrm{cos}\alpha \text{ }。\end{array}\right\}$

(16)

式中， $\alpha$ 角为在点1处块体边界线与水平线的夹角。

对于点2上的水压力沿x，y方向分解可得到相似的结果。类似于DDA的体积力荷载子矩阵，点1和点2的水压力对块体i产生的势能为

${\varPi _i} = - \int\limits_0^{{l_{12}}} {({f_{{\text{xs}}}}{u_{\text{s}}} + {f_{{\text{ys}}}}{v_{\text{s}}})} {\text{d}}s 。$

(17)

式中：l₁₂为点1到点2的距离；s为在该裂隙长度方向上选取的某一微元段；u_s和v_s为该微元段在x，y方向上的平移分量。

对式（17）进行势能最小化处理：

$- \frac{{\partial {\Pi _i}(0)}}{{\partial {d_{ri}}}} = {\boldsymbol T_i}^{\text{T}}({x_1},{y_1})\left[ E \right] + {\boldsymbol T_i}^{\text{T}}({x_2},{y_2})\left[ G \right]$

(18)

式中： ${\boldsymbol T_i}^{\text{T}}$ (x₁, y₁)和T_i^T(x₂, y₂)分别为点1和点2的位移转换矩阵，[E]和[G]由下式表示：

$\begin{array}{l}\left[E\right]=\left(\begin{array}{cc}-\frac{{l}_{12}}{3}\mathrm{sin}\alpha & -\frac{{l}_{12}}{6}\mathrm{sin}\alpha \\ \frac{{l}_{12}}{3}\mathrm{cos}\alpha & \frac{{l}_{12}}{6}\mathrm{cos}\alpha \end{array}\right)\left\{\begin{array}{c}{p}_{1}\\ {p}_{2}\end{array}\right\}\text{ }\text{，}\\ \left[G\right]=\left(\begin{array}{cc}-\frac{{l}_{12}}{6}\mathrm{sin}\alpha & -\frac{{l}_{12}}{3}\mathrm{sin}\alpha \\ \frac{{l}_{12}}{6}\mathrm{cos}\alpha & \frac{{l}_{12}}{3}\mathrm{cos}\alpha \end{array}\right)\left\{\begin{array}{c}{p}_{1}\\ {p}_{2}\end{array}\right\}\text{ }。\end{array}$

(19)

式（18）的计算结果为一个6×1的子矩阵，将其加入到DDA总体平衡方程的荷载项F_i中，即式(20)的右侧：

$\left[\begin{array}{cccc}{K}_{11}& {K}_{12}& \dots & {K}_{1n}\\ {K}_{21}& {K}_{22}& \dots & {K}_{2n}\\ ⋮& ⋮& ⋮& ⋮\\ {K}_{n1}& {K}_{n2}& \dots & {K}_{nn}\end{array}\right]\left\{\begin{array}{c}{D}_{1}\\ {D}_{2}\\ ⋮\\ {D}_{n}\end{array}\right\}=\left\{\begin{array}{c}{F}_{1}\\ {F}_{2}\\ ⋮\\ {F}_{n}\end{array}\right\}\text{ }。$

(20)

式中，K_ij为6×6的子矩阵，用于描述块体之间的相互作用，D_i和F_i为6×1的子矩阵用以描述块体上的位移和荷载。

通过式（18）可以将渗流场计算的结果作为外荷载施加到应力应变场计算中，实现渗流场对应力应变场的影响。外荷载的施加则会改变下一时步岩体裂隙的开度，反向施加影响于渗流场，进而实现渗流场与应变场的双向耦合。

3. 算例验证

3.1 土柱入渗

该算例选取长1 m土柱计算水体从其顶部入渗的压力水头分布^[19]。在模型顶部初始施加-7 m的水头，底部施加-8 m的水头。当降雨开始后模型底部水头不发生变化，顶部施加1 m的水头作为边界条件，如图 2所示。

图 2 计算模型及边界条件

Figure 2. Computational model and boundary conditions

下载: 全尺寸图片幻灯片

该模型的具体计算参数如表 1所示，将46800 s的压力水头数值解与解析解进行对比。

表 1 模型参数

Table 1. Model parameters

初始裂隙开度/mm	θ_r	θ_s	α/m^-1	n	Δt/s
0.111	0.363	0.168	1	1.53	180

下载: 导出CSV

| 显示表格

如图 3所示，数值解与解析解的压力水头均呈现随高程下降而逐渐降低的趋势，且模型上下两端压力水头随高程变化较慢而模型中部压力水头随高程变化较快，这是由于模型上下两端设置了固定水头作为边界条件，而模型中部受边界条件影响较小。与解析解相比，数值解模型顶部压力水头上升较快，原因可能为当块体体积含水量接近饱和体积含水量时，压力水头会近似处理成饱和状态压力水头，从而导致压力水头上升较快。数值解与解析解的压力水头具有相似的变化规律，验证了该数值方法的正确性。

图 3 46800 s压力水头分布

Figure 3. Distribution of pressure head at 46800 s

下载: 全尺寸图片幻灯片

3.2 砂槽入渗试验

该算例采用砂槽入渗试验对程序的正确性进行验证，该试验常用于检验饱和-非饱和非稳定渗流算法的性能^[20]。试验区域为一个高2 m宽6 m的矩形区域。边界条件为两侧边界自由排水、底部边界不透水、上部边界降雨。地下水位在0.65 m处，在试验开始后，对该试验区域中心1 m的范围内均匀施加0.148 m/h的降雨强度，降雨时长为8 h。初始裂隙宽度为0.0003 m，饱和渗透系数为9.72×10^-5 m/s，其数值模型如图 4所示。

图 4 算例模型

Figure 4. Example model

下载: 全尺寸图片幻灯片

由于该模型是完全对称的，因此截取模型的右半段进行分析。完成建模后按照表 2参数对模型进行赋值并计算。

表 2 模型参数^[19]

Table 2. Model parameters

初始裂隙开度/mm	θ_r	θ_s	α/m^-1	n	Δt/s
1.29	0.01	0.3	3.3	4.1	3600

下载: 导出CSV

| 显示表格

降雨开始2 h后的压力水头等值线如图 5所示。

图 5 2 h压力水头等值线

Figure 5. Contours of pressure head at 2 h

下载: 全尺寸图片幻灯片

由于在模型的左上角施加降雨边界，在该边界条件下的地层压力水头较大，随着右侧边界条件变为自由边界，压力水头迅速下降。符合在降雨区域地层饱和度上升负孔隙水压消散，在未降雨区域地层状态与初始状态接近这一规律。整体模型压力水头等值线呈左高右低的趋势，且随着深度降低等值线越发趋于平缓。符合浅表地层受降雨影响较大，体积含水量上升，压力水头上升，而深处地层受到影响较小的规律。与NMM的计算结果主要在降雨边界附近存在差异，其主要原因是NMM可以通过数学片对块体进一步细分，而DDA只在块体的边界有计算点。因此DDA模拟结果的等值线分布相NMM要更为稀疏，但总体分布规律与NMM相似。该模型主要研究降雨引起地下水位的变化，对于地下水位以下部分未进行插值。

降雨开始8 h后的压力水头等值线如图 6所示。随着降雨的进行，模型的压力水头整体呈上升趋势，尤其是在降雨边界处，由2 h的-0.49 m上升到了-0.36 m。地下水位线也随降雨进行而有所升高，且呈现左高右低的趋势，符合离降雨边界较近处饱和度上升较快的规律。等值线出现偏折是由于受到边界条件影响，在降雨边界和自由边界的分界处产生计算数值上的突变。DDA计算结果压力水头的分布规律与NMM计算结果相似。选取试验区域中央截面，分析地下水位随时间上升的规律如图 7所示。

图 6 8 h压力水头等值线

Figure 6. Contours of pressure head at 8 h

下载: 全尺寸图片幻灯片

图 7 地下水位线变化

Figure 7. Variation of groundwater level line

下载: 全尺寸图片幻灯片

图 7将Vuaclin的室内试验结果和DDA数值模拟的结果进行对比。DDA数值模拟水位线上升的高度稍小于室内试验得到的结果，与试验结果相比最大误差出现在降雨4 h后，为3.74%。因此可以认为数值模拟结果与室内试验结果较为接近，进而认为该程序可以应用在非饱和瞬态渗流相关的模拟当中。

4. 工程案例

将改进后的DDA程序应用到重庆市巫山县三溪乡大梯子岩危岩体中。根据工程地质剖面图建立如图 8所示的计算模型。模型的力学参数如表 3所示。

图 8 计算模型

Figure 8. Model

下载: 全尺寸图片幻灯片

表 3 模型力学参数

Table 3. Mechanical parameters for model

材料	密度/(kg·m^-3)	弹性模量/GPa	泊松比	抗拉强度/MPa	黏聚力/MPa	摩擦角/(°)
危岩体	2400	40	0.3	10	10	35
基岩	2400	50	0.25	20	50	50

下载: 导出CSV

| 显示表格

在自然工况下对危岩体进行数值计算，结果如图 9所示。如图 9所示，在自然工况下，位移最大值点出现在危岩体的底部，其值为4.35 mm。可能由于差异风化导致危岩体底部出现凹腔，在上部自重等荷载的作用下导致危岩体在底部产生较大位移。其此时危岩体的整体稳定性系数为1.2127，危岩体整体处于相对稳定状态。

图 9 自然工况位移云图

Figure 9. Nephogram of displacement under natural condition

下载: 全尺寸图片幻灯片

对于危岩体的非饱和渗流采用Van Genuchten模型描述，具体取值如表 4所示。

表 4 渗流参数

Table 4. Seepage parameters

初始裂隙开度/mm	Δt/s	n	α/m^-1	${\theta _{\text{s}}}$	${\theta _{\text{r}}}$
0.3	1800	2	0.1	0.39	0.05

下载: 导出CSV

| 显示表格

根据当地降雨统计数据模拟20 h强度为7.07 mm/h的较长时间持续降雨，得到结果如图 10所示。

图 10 降雨工况体积含水量云图

Figure 10. Nephogram of volume water content under rainfall condition

下载: 全尺寸图片幻灯片

经过20 h降雨作用后，危岩体顶部体积含水量接近饱和体积含水量，产生饱和度大于0.95的暂态饱和区^[21]，而危岩体下部体积含水量则几乎不受到降雨的影响。由于降雨边界只在该模型顶部设置，因此底部节点的体积含水量受到降雨影响小，降雨从危岩体顶部入渗且无法渗透整个危岩体。危岩体右侧降雨入渗深度相对于左侧较深：一方面由于危岩体顶部有一定坡度，右侧高度低于左侧更有利于雨水入渗；另一方面，降雨入渗计算网格与切割岩体裂隙网格相同，岩体右侧裂隙长度较左侧更短，因此在相同的降雨条件下，危岩体右侧体积含水量更容易受到影响。而岩体顶部最右侧临空块体由于没有联通的裂隙对其进行切割，属于孤立块体，因此不在降雨入渗计算网格当中。危岩体的总位移云图如图 11所示。

图 11 降雨工况总位移云图

Figure 11. Nephogram of total displacement under rainfall condition

下载: 全尺寸图片幻灯片

危岩体最大位移为2.6 cm，出现在模型顶部的临空块体，发生较大位移的区域位于危岩体底部，其位移为1.8 cm。与自然工况相比，危岩体的位移有明显增长。危岩体顶部受到降雨影响，块体的强度有一定弱化，且顶部受到降雨形成的水头压力使其产生较大的位移。底部块体由于顶部块体的重力和变形增大，也产生了较大的位移。此时危岩体的整体稳定性系数为0.6304，危岩体整体处于不稳定状态。危岩体顶部和底部块体出现裂缝扩展、错动等大变形现象，如图 12所示。

图 12 块体大变形

Figure 12. Large deformations of block

下载: 全尺寸图片幻灯片

该模拟结果符合工程实际规律，且能用于描述块体大变形过程，验证了本文提出的方法在工程应用中的可行性。

5. 结论

针对DDA难以进行瞬态渗流场模拟这一问题对其进行二次开发，将渗流场计算引入到原始DDA当中，形成一种流固耦合的DDA方法并通过相关经典算例和工程案例进行验证。提出了一种可以反映降雨入渗过程的DDA方法，为实际工程提供新的数值模拟手段，得到以下2点结论。

（1）对土柱入渗和砂槽入渗的模拟结果表明，本文建立的非饱和瞬态流固耦合计算模型与试验结果和NMM计算模型相比具有较好的一致性。能够对降雨造成的地下水位抬升、负压水头消散等变化过程进行模拟，拓展了DDA的适用范围。

（2）对大梯子岩危岩体的模拟结果表明，降雨导致危岩体顶部出现暂态饱和区，造成岩体强度参数弱化以及裂隙水压力上升，从而产生较大的位移。在顶部块体的重力和位移作用下，危岩体底部块体产生相对较大的位移，危岩体整体处于欠稳定状态。在危岩体治理中应做好排水措施，防止降雨造成危岩体失稳崩塌。

图 1 平面和曲面土粒吸附水膜示意图

Figure 1. Illustration of adsorbed water films around flat surface and curved surface

下载: 全尺寸图片幻灯片

图 2 土单元和土粒碎片示意图

Figure 2. Illustration of soil matrix and divided fraction

下载: 全尺寸图片幻灯片

图 3 截面A上不同土粒碎片间的连接示意图

Figure 3. Illustration of connection between different fractions on section A

下载: 全尺寸图片幻灯片

图 4 5种验证土的性质

Figure 4. Properties of five verification soils

下载: 全尺寸图片幻灯片

图 5 最佳拟合土-水特征曲线（单峰和双峰）

Figure 5. Best fitted SWCCs for soils (including both unimodal and bimodal SWCCs)

下载: 全尺寸图片幻灯片

图 6 预测渗透系数与实测渗透系数之间的比较

Figure 6. Comparisons between estimated results and measured hydraulic conductivities for soils

下载: 全尺寸图片幻灯片

表 1 不同种类土质的有效粒径

Table 1 Effective particle sizes of soils collected from different literatures

土的类型	有效粒径/mm
土的类型	d₉₅	d₈₅	d₇₅	d₆₅	d₅₅	d₄₅	d₃₅	d₂₅	d₁₅	d₅
亚砂土	0.230	0.188	0.151	0.122	0.099	0.081	0.067	0.053	0.017	0.00025
亚黏土	0.500	0.397	0.345	0.197	0.075	0.021	0.009	0.005	0.004	0.00020
UNSODA 4010	0.358	0.191	0.161	0.133	0.111	0.090	0.067	0.049	0.031	0.00250
UNSODA 4031	0.210	0.091	0.081	0.068	0.058	0.018	0.008	0.017	0.010	0.00800
UNSODA 4650	1.140	0.561	0.481	0.402	0.345	0.296	0.251	0.211	0.114	0.00350
SM 6-62	0.075	0.060	0.055	0.048	0.040	0.030	0.020	0.010	0.006	0.00200
SM 35-119	0.060	0.050	0.040	0.030	0.020	0.010	0.008	0.003	0.002	0.00100

下载: 导出CSV

表 2 单峰土-水特征曲线的拟合参数

Table 2 Fitting parameters for estimated unimodal SWCCs

土的类型	${a_{\text{f}}}$ /kPa	${n_{\text{f}}}$	${m_{\text{f}}}$	C_r/kPa
亚砂土	12.12	1.13	1.36	1500
亚黏土	4.29	4.33	0.51	1500
UNSODA 4031	5.06	0.84	0.85	1500
UNSODA 4650	1.97	3.03	0.91	1500
SM 6-62	5.04	1.31	0.65	1500
SM 35-119	3.82	1.22	0.49	1500
注： ${a_{\text{f}}}$ ， ${n_{\text{f}}}$ ， ${m_{\text{f}}}$ 为Fredlund和Xing公式中的拟合参数。

下载: 导出CSV

表 3 双峰土-水特征曲线的拟合参数

Table 3 Fitting parameters for bimodal SWCC

土的类型	Fredlund和Xing式中的拟合参数
土的类型	θ_s1	${a_{{\text{f1}}}}$ /kPa	${n_{{\text{f1}}}}$	${m_{{\text{f1}}}}$	θ_s2	${a_{{\text{f2}}}}$ /kPa	${n_{{\text{f2}}}}$	${m_{{\text{f2}}}}$
UNSODA 4010	1	2.91	2.78	0.26	0.725	6.45	9.22	0.36
注：θ_s1， ${a_{{\text{f1}}}}$ ， ${n_{{\text{f1}}}}$ ， ${m_{{\text{f1}}}}$ ，θ_s2， ${a_{{\text{f2}}}}$ ， ${n_{{\text{f2}}}}$ ， ${m_{{\text{f2}}}}$ 为Fredlund和Xing式双峰SWCC中的拟合参数。

下载: 导出CSV

表 4 土的残余吸力ψ_r和ψ_film

Table 4 Residual suctions ψ_r and ψ_film for soils

参数	亚砂土	亚黏土	UNSODA 4010	UNSODA 4031	UNSODA 4650
ψ_r/kPa	124.63	13.38	12.58	20	7.9
ψ_film/kPa	500	70	50	50	70

下载: 导出CSV

参考文献(26)

[1]	WEBER T K D, DURNER W, STRECK T, et al. A modular framework for modeling unsaturated soil hydraulic properties over the full moisture range[J]. Water Resources Research, 2019, 55(6): 4994-5011. doi: 10.1029/2018WR024584
[2]	翟钱, 朱益瑶, 叶为民, 等. 全吸力范围非饱和土水力渗透系数的计算[J]. 岩土工程学报, 2022, 44(4): 660-668. doi: 10.11779/CJGE202204008 ZHAI Qian, ZHU Yiyao, YE Weimin, et al. Estimation of hydraulic conductivity of unsaturated soils under entire suction range[J]. Chinese Journal of Geotechnical Engineering, 2022, 44(4): 660-668. (in Chinese) doi: 10.11779/CJGE202204008
[3]	GARDNER W R. Mathematics of isothermal water conduction in unsaturated soil[J]. Highway Research Board Special Report 1958, 40: 78-87
[4]	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
[5]	FREDLUND D G, XING A Q, HUANG S Y. Predicting the permeability function for unsaturated soils using the soil-water characteristic curve[J]. Canadian Geotechnical Journal, 1994, 31(4): 533-546. doi: 10.1139/t94-062
[6]	ZHAI Q, RAHARDJO H. Estimation of permeability function from the soil-water characteristic curve[J]. Engineering Geology, 2015, 199: 148-156. doi: 10.1016/j.enggeo.2015.11.001
[7]	叶为民, 钱丽鑫, 白云, 等. 由土-水特征曲线预测上海非饱和软土渗透系数[J]. 岩土工程学报, 2005, 27(11): 1262-1265. doi: 10.3321/j.issn:1000-4548.2005.11.005 YE Weimin, QIAN Lixin, BAI Yun, et al. Predicting coefficient of permeability from soil-water characteristic curve for Shanghai soft soil[J]. Chinese Journal of Geotechnical Engineering, 2005, 27(11): 1262-1265. (in Chinese) doi: 10.3321/j.issn:1000-4548.2005.11.005
[8]	蔡国庆, 盛岱超, 周安楠. 考虑初始孔隙比影响的非饱和土相对渗透系数方程[J]. 岩土工程学报, 2014, 36(5): 827-835. doi: 10.11779/CJGE201405004 CAI Guoqing, SHENG Daichao, ZHOU Annan. Approach for predicting the relative coefficient of permeability of unsaturated soils with different initial void ratios[J]. Chinese Journal of Geotechnical Engineering, 2014, 36(5): 827-835. (in Chinese) doi: 10.11779/CJGE201405004
[9]	TOKUNAGA T K, WAN J M. Water film flow along fracture surfaces of porous rock[J]. Water Resources Research, 1997, 33(6): 1287-1295. doi: 10.1029/97WR00473
[10]	TOKUNAGA T K, WAN J M, SUTTON S R. Transient film flow on rough fracture surfaces[J]. Water Resources Research, 2000, 36(7): 1737-1746. doi: 10.1029/2000WR900079
[11]	LIU H H. A constitutive-relationship model for film flow on rough fracture surfaces[J]. Hydrogeology Journal, 2004, 12(2): 237-240.
[12]	TOKUNAGA T K. Hydraulic properties of adsorbed water films in unsaturated porous media[J]. Water Resources Research, 2009, 45(6): W06415.
[13]	LEBEAU M, KONRAD J M. A new capillary and thin film flow model for predicting the hydraulic conductivity of unsaturated porous media[J]. Water Resources Research, 2010, 46(1): 1-15.
[14]	TULLER M, OR D. Hydraulic conductivity of variably saturated porous media: film and corner flow in angular pore space[J]. Water Resources Research, 2001, 37(5): 1257-1276. doi: 10.1029/2000WR900328
[15]	IWAMATSU M, HORII K. Capillary condensation and adhesion of two wetter surfaces[J]. Journal of Colloid and Interface Science, 1996, 182(2): 400-406. doi: 10.1006/jcis.1996.0480
[16]	BIRD R B, STEWART W E, LIGHTFOOT E N. Transport phenomena[M]. New York: John Wiley & Sons, Inc, 1960.
[17]	PACHEPSKY Y A, SHCHERBAKOV R, VARALLYAY G, et al. On obtaining soil hydraulic conductivity curves from water retention curves [J]. Pochvovedenie, 1984, 10: 60-72.
[18]	MUALEM Y. Hydraulic conductivity of unsaturated soils: prediction and formulas[M]//SSSA Book Series. Madison, WI, USA: Soil Science Society of America, American Society of Agronomy, 2018: 799-823.
[19]	NEMES A, SCHAAP M G, LEIJ F J, et al. Description of the unsaturated soil hydraulic database UNSODA version 2.0[J]. Journal of Hydrology, 2001, 251(3/4): 151-162.
[20]	SCHINDLER U G, MÜLLER L. Soil hydraulic functions of international soils measured with the Extended Evaporation Method (EEM) and the HYPROP device[J]. Open Data Journal for Agricultural Research, 2017, 3: 1-7.
[21]	FREDLUND D G, XING A Q. Equations for the soil-water characteristic curve[J]. Canadian Geotechnical Journal, 1994, 31(4): 521-532. doi: 10.1139/t94-061
[22]	ZHAI Q, RAHARDJO H, SATYANAGA A. Effect of bimodal soil-water characteristic curve on the estimation of permeability function[J]. Engineering Geology, 2017, 230: 142-151. doi: 10.1016/j.enggeo.2017.09.025
[23]	ZHAI Q, RAHARDJO H. Determination of soil-water characteristic curve variables[J]. Computers and Geotechnics, 2012, 42: 37-43. doi: 10.1016/j.compgeo.2011.11.010
[24]	ZHANG F X, FREDLUND D G. Examination of the estimation of relative permeability for unsaturated soils[J]. Canadian Geotechnical Journal, 2015, 52(12): 2077-2087. doi: 10.1139/cgj-2015-0043
[25]	ZHAI Q, RAHARDJO H, SATYANAGA A, et al. Framework to estimate the soil-water characteristic curve for soils with different void ratios[J]. Bulletin of Engineering Geology and the Environment, 2020, 79(8): 4399-4409. doi: 10.1007/s10064-020-01825-8
[26]	ZHAI Q, RAHARDJO H, SATYANAGA A. Effects of residual suction and residual water content on the estimation of permeability function[J]. Geoderma, 2017, 303: 165-177. doi: 10.1016/j.geoderma.2017.05.019

施引文献(5)

期刊类型引用(2)

1.	邱军领，贾玎，赖金星，唐琨杰，强磊. 含裂隙土质隧道降雨入渗双通道渗流模型. 岩土工程学报. 2025(03): 548-558 . 本站查看
2.	张小东，路喆津，刘敏，杨峰，赵炼恒. 强降雨作用下引入分形理论的浅层边坡入渗及稳定性分析. 哈尔滨工业大学学报. 2024(11): 140-150 . 百度学术