非饱和土有限元渗流计算：控制方程、离散与非线性迭代

饱和渗流中，渗透系数通常可以视为常数；非饱和土却不同，含水率与渗透系数都会随压力水头变化。降雨入渗时，某个高斯点的压力水头只要略有改变，就可能使该点的渗透系数跨越数个数量级，进而改变单元渗透矩阵。非饱和渗流计算的核心困难不是线性方程求解，而是本构参数与未知压力场之间的强非线性耦合。

本文以二维瞬态渗流为对象，采用压力水头 $\psi$ 作为基本未知量，依次说明质量守恒、Richards 方程、土水特征曲线、有限元弱式、时间离散和高斯点迭代。最后给出适合程序实现的收敛判据与排错清单。

开始前需要统一约定： $x$ 水平向右， $z$ 竖直向上；压力水头 $\psi=u_w/\gamma_w$ ，非饱和区内通常有 $\psi<0$ ；总水头为 $H=\psi+z$ ；源项 $Q>0$ 表示向土体内补水。

一、从微元水量平衡到 Richards 方程

1. 质量守恒

取尺寸为 $dx\times dz$ 、单位厚度的二维微元。令达西通量 $\boldsymbol{q}=(q_x,q_z)^T$ ，通量方向按坐标正方向取正。微元内储水量的增加等于净流入量与源项之和，因此有

$\frac{\partial \theta}{\partial t} +\boldsymbol{\nabla}\cdot\boldsymbol{q}=Q \tag{1}\label{eq:mass}$

其中， $\theta$ 为体积含水率， $Q$ 的量纲为 $\mathrm{s^{-1}}$ 。若程序使用节点流量而不是体积源强度，需要先根据单元面积和厚度完成量纲换算。

2. 达西定律与符号检查

各向异性介质中的达西定律写为

$\boldsymbol{q} =-\boldsymbol{K}(\psi)\boldsymbol{\nabla}H =-\boldsymbol{K}(\psi) \left(\boldsymbol{\nabla}\psi+\boldsymbol{e}_z\right) \tag{2}\label{eq:darcy}$

若材料主轴与坐标轴重合，则

$\boldsymbol{K}= \begin{bmatrix} K_x(\psi) & 0\\ 0 & K_z(\psi) \end{bmatrix}$

负号表示水从高总水头流向低总水头。实现时最容易出错的是重力项：当 $z$ 轴向上时， $\boldsymbol{\nabla}z=\boldsymbol{e}_z$ ；若程序把竖直坐标向下定义，重力项必须同步改变符号。

将式 $\eqref{eq:darcy}$ 代入式 $\eqref{eq:mass}$ ，得到混合形式的 Richards 方程：

$\boxed{ \frac{\partial\theta(\psi)}{\partial t} =\boldsymbol{\nabla}\cdot \left[ \boldsymbol{K}(\psi) \left(\boldsymbol{\nabla}\psi+\boldsymbol{e}_z\right) \right]+Q } \tag{3}\label{eq:richards}$

它同时适用于非饱和区和饱和区。非线性来自两个位置：储水关系 $\theta(\psi)$ 和渗透关系 $\boldsymbol{K}(\psi)$ 。COMSOL 的 Richards 方程接口同样把压力作为未知量，并允许水力性质随饱和状态改变；这也是此类问题需要非线性迭代的根本原因。

二、土水特征曲线与非饱和渗透系数

1. 有效饱和度

定义有效饱和度

$S_e=\frac{\theta-\theta_r}{\theta_s-\theta_r} \tag{4}\label{eq:se}$

其中， $\theta_s$ 和 $\theta_r$ 分别为饱和与残余体积含水率。采用 van Genuchten 模型时，常写为

$S_e(\psi)= \begin{cases} \left[1+(\alpha|\psi|)^n\right]^{-m}, & \psi<0\\ 1, & \psi\ge 0 \end{cases}, \qquad m=1-\frac{1}{n} \tag{5}\label{eq:vg}$

于是

$\theta(\psi)=\theta_r+(\theta_s-\theta_r)S_e(\psi) \tag{6}\label{eq:theta}$

$\alpha$ 控制进气吸力的尺度， $n$ 控制曲线陡峭程度。van Genuchten 模型的连续形式不仅便于拟合试验数据，也能避免不连续斜率给数值迭代带来的额外困难。

2. 容水度

压力水头形式的容水度定义为

$C(\psi)=\frac{d\theta}{d\psi} \tag{7}\label{eq:capacity}$

当 $\psi<0$ 时，式 $\eqref{eq:vg}$ 对应

$C(\psi)= (\theta_s-\theta_r)mn\alpha^n|\psi|^{n-1} \left[1+(\alpha|\psi|)^n\right]^{-m-1} \tag{8}\label{eq:cvg}$

若程序以孔隙水压力 $u_w$ 为未知量，则

$\frac{d\theta}{du_w} =\frac{1}{\gamma_w}\frac{d\theta}{d\psi} =\frac{C(\psi)}{\gamma_w} \tag{9}\label{eq:pressure-capacity}$

因此不能把针对压力水头标定的容水度直接用于孔压方程，否则会遗漏 $\gamma_w$ ，导致储水项出现数量级错误。

3. Mualem–van Genuchten 渗透关系

各向同性时，非饱和渗透系数可写成

$K(\psi)=K_sK_r(S_e) \tag{10}\label{eq:k}$

其中 $K_s$ 为饱和渗透系数。常用的 Mualem–van Genuchten 相对渗透系数为

$K_r(S_e)= S_e^{\ell} \left[ 1-\left(1-S_e^{1/m}\right)^m \right]^2 \tag{11}\label{eq:kr}$

$\ell$ 为孔隙连通参数，常用初值为 $0.5$ ，但最终应由试验或可靠资料确定。原始研究表明，低含水率段的土水特征曲线会显著影响非饱和渗透系数预测，因此不能只关注接近饱和区的拟合效果。

三、有限元弱式

在区域 $\Omega$ 内取试函数 $N_i$ 。把式 $\eqref{eq:richards}$ 乘以 $N_i$ 并分部积分，若边界入流通量 $\bar q_{\mathrm{in}}>0$ 定义为流入区域，可得

$\int_{\Omega}N_i\frac{\partial\theta}{\partial t}\,d\Omega +\int_{\Omega}\boldsymbol{\nabla}N_i^T \boldsymbol{K}(\psi) \left(\boldsymbol{\nabla}\psi+\boldsymbol{e}_z\right)d\Omega =\int_{\Omega}N_iQ\,d\Omega +\int_{\Gamma_q}N_i\bar q_{\mathrm{in}}\,d\Gamma \tag{12}\label{eq:weak}$

采用单元形函数插值

$\psi\approx\boldsymbol{N}\boldsymbol{\psi}_e, \qquad \boldsymbol{\nabla}\psi =\boldsymbol{B}\boldsymbol{\psi}_e \tag{13}\label{eq:interp}$

在一次线性化迭代中，单元容量矩阵和渗透矩阵分别为

$\boldsymbol{C}_e =\int_{\Omega_e}\boldsymbol{N}^T C(\psi)\boldsymbol{N}\,d\Omega \tag{14}\label{eq:capacity-matrix}$

$\boldsymbol{K}_e =\int_{\Omega_e}\boldsymbol{B}^T \boldsymbol{K}(\psi)\boldsymbol{B}\,d\Omega \tag{15}\label{eq:conductivity-matrix}$

重力项为

$\boldsymbol{g}_e =\int_{\Omega_e}\boldsymbol{B}^T \boldsymbol{K}(\psi)\boldsymbol{e}_z\,d\Omega \tag{16}\label{eq:gravity-vector}$

这些积分通常在高斯点计算。对四节点四边形单元， $2\times2$ 积分对应四个高斯点；每个点都要根据当前压力水头重新计算 $S_e$ 、 $\theta$ 、 $C$ 与 $K$ ，再组装到单元矩阵中。

四、时间离散与非线性残量

使用后向欧拉格式，从时刻 $t_n$ 前进到 $t_{n+1}=t_n+\Delta t$ 。相比直接把 $C\,\partial\psi/\partial t$ 离散，混合形式保留含水率增量：

$\frac{\partial\theta}{\partial t} \approx \frac{\theta(\psi^{n+1})-\theta(\psi^n)}{\Delta t} \tag{17}\label{eq:time}$

离散残量可以写为

$\begin{aligned} \boldsymbol{R}(\boldsymbol{\psi}^{n+1}) =&\sum_e\int_{\Omega_e}\boldsymbol{N}^T \frac{\theta(\psi^{n+1})-\theta(\psi^n)}{\Delta t}\,d\Omega\\ &+\sum_e\int_{\Omega_e}\boldsymbol{B}^T \boldsymbol{K}(\psi^{n+1}) \left(\boldsymbol{B}\boldsymbol{\psi}^{n+1} +\boldsymbol{e}_z\right)d\Omega -\boldsymbol{f}^{n+1} \end{aligned} \tag{18}\label{eq:residual}$

目标是求得 $\boldsymbol{R}=\boldsymbol{0}$ 。Celia 等人的研究指出，混合形式比单纯的压力水头形式更有利于保持质量守恒；时间项采用合适的集中处理还可减小非物理振荡。

五、高斯点本构更新与收敛

1. Picard 迭代框架

在时间步 $n+1$ 的第 $m$ 次迭代中，以 $\boldsymbol{\psi}^{(m)}$ 计算所有高斯点本构量，组装近似切线矩阵，然后求解增量：

$\boldsymbol{A}^{(m)}\Delta\boldsymbol{\psi}^{(m)} =-\boldsymbol{R}^{(m)} \tag{19}\label{eq:increment}$

$\boldsymbol{\psi}^{(m+1)} =\boldsymbol{\psi}^{(m)} +\omega\Delta\boldsymbol{\psi}^{(m)}, \qquad 0<\omega\le1 \tag{20}\label{eq:relaxation}$

$\omega$ 为松弛系数。普通时间步可取 $\omega=1$ ；干湿交界面移动剧烈或迭代振荡时，可暂时减小到 $0.3\sim0.8$ 。完整 Newton 法还会把 $dK/d\psi$ 对残量的贡献放入雅可比矩阵，局部收敛更快，但实现和调试成本也更高。

2. 为什么不能只检查节点压力

当土体很干时， $K(\psi)$ 曲线可能非常陡。节点压力增量已经很小，不代表高斯点渗透系数已经稳定；反过来，绝对渗透系数极小时，单纯计算相对误差又容易被浮点噪声放大。因此建议同时监控三类指标。

压力水头增量：

$\eta_{\psi} =\frac{\|\Delta\boldsymbol{\psi}\|_{\infty}} {\max(\|\boldsymbol{\psi}^{(m+1)}\|_{\infty},\psi_{\mathrm{ref}})} \tag{21}\label{eq:head-error}$

高斯点渗透系数变化：

$\eta_K =\max_g \frac{|K_g^{(m+1)}-K_g^{(m)}|} {\max(K_g^{(m+1)},K_{\mathrm{floor}})} \tag{22}\label{eq:k-error}$

平衡残量：

$\eta_R =\frac{\|\boldsymbol{R}^{(m+1)}\|_2} {\max(\|\boldsymbol{f}^{n+1}\|_2,R_{\mathrm{ref}})} \tag{23}\label{eq:residual-error}$

只有三项同时满足

$\eta_{\psi}<\varepsilon_{\psi}, \qquad \eta_K<\varepsilon_K, \qquad \eta_R<\varepsilon_R \tag{24}\label{eq:convergence}$

才接受当前时间步。容差必须结合单位、网格与材料曲线确定；调试初期可从 $10^{-5}\sim10^{-4}$ 的相对容差开始，再通过时间步和网格收敛试验收紧。 $K_{\mathrm{floor}}$ 只用于误差归一化，不应偷偷改变物理本构值。

3. 质量平衡复核

非线性迭代收敛后还应检查整个区域的水量平衡：

$E_m= \Delta W -\Delta t\left(Q_{\mathrm{in}}-Q_{\mathrm{out}}+Q_{\mathrm{source}}\right) \tag{25}\label{eq:mass-error}$

其中

$\Delta W =\int_{\Omega} \left[\theta^{n+1}-\theta^n\right]d\Omega \tag{26}\label{eq:storage-change}$

若残量已经很小而 $E_m$ 仍持续累积，应优先检查时间项是否使用了含水率增量、边界通量符号、积分权重、单元厚度和节点流量分配，而不是盲目增加迭代次数。

六、边界条件的程序表达

1. 定水头边界

给定 $\psi=\bar\psi$ 或 $H=\bar H$ 。必须确认输入量究竟是孔压、压力水头还是总水头，三者不能混用。施加本质边界后，应同步修改矩阵与右端项，避免只覆盖解向量而破坏平衡方程。

2. 流量边界

给定法向入流 $\bar q_{\mathrm{in}}$ 。降雨强度并不总能全部转化为入渗通量：当表面达到饱和或所需入渗能力不足时，边界可能需要从给定通量切换为给定压力，并把剩余水量计为地表径流。

3. 无流边界

有 $\boldsymbol{q}\cdot\boldsymbol{n}=0$ 。在弱式中它属于自然边界，不需要额外装配边界向量，但仍要确认模型底部是否真的可以视为不透水层。

4. 渗出面

渗出面通常同时包含压力约束与单向出流条件，边界范围会随解变化。若直接把整个潜在渗出边界固定为零压力，往往会高估排水范围；更稳妥的做法是在每次非线性迭代后更新活动边界集合。

七、建议的程序流程

读取网格、材料、初始条件和边界条件
计算初始高斯点 θ、C、K

for n = 0 ... nstep - 1
    取 ψ^(0) = ψ^n

    for m = 0 ... iter_max - 1
        for 每个单元
            for 每个高斯点
                插值得到 ψ_g
                更新 Se_g、θ_g、C_g、K_g
                计算容量、渗透、重力和源项贡献
            end for
            组装全局残量 R 和切线矩阵 A
        end for

        施加本质边界
        求解 A Δψ = -R
        更新 ψ^(m+1) = ψ^(m) + ω Δψ
        重新计算 η_ψ、η_K、η_R

        if 三项均收敛
            接受 ψ^(m+1)
            exit
        end if
    end for

    if 未收敛
        减小 Δt，回滚到 ψ^n 后重算
    else
        检查全局水量平衡
        输出节点压力、高斯点饱和度与通量
        根据迭代次数调整下一步 Δt
    end if
end for

时间步自适应通常比单纯提高最大迭代次数更有效：若连续几个时间步只需少量迭代，可以适度放大 $\Delta t$ ；若迭代超过上限、误差振荡或边界状态频繁切换，则应回滚并减小 $\Delta t$ 。

八、常见不收敛原因

现象	常见原因	优先处理方法
第一个时间步立即发散	初始场与边界条件不协调	先求稳态初始场，或缓慢加载边界
压力振荡、渗透系数来回跳变	时间步过大， $K(\psi)$ 过陡	减小时间步并使用松弛更新
节点压力收敛但流量不稳定	只检查了 $\Delta\psi$	增加高斯点 $K$ 与残量判据
水量误差持续累积	时间项、边界通量或单位错误	用式 $\eqref{eq:mass-error}$ 分项核对
饱和线附近难以收敛	本构关系在 $\psi=0$ 附近不光滑	检查分段函数及一阶导数连续性
干燥区矩阵病态	$K$ 过小、网格质量差	改善网格，使用缩放和稳健预条件
降雨通量导致正压力异常增长	未实现通量—压力边界切换	加入积水与径流判定