1D Richards 方程隐式差分中如何实现 Neumann–Dirichlet 上边界切换？

2026-11-04数据阿明

可以先以降雨通量作为 Neumann 边界求解。如果收敛后的表面压力水头超过 20 cm，就将边界改为 h = 20 cm 的 Dirichlet 条件，重新求解同一个时间步。切换后还需检查实际入渗通量是否超过可用降雨量。

这套方法适用于“不显式计算地表储水、超额降雨直接排走”的限压边界。如果 20 cm 表示真实的最大积水深度，还需要联立地表储水方程。仅修改边界的 if 条件，无法描述雨停后积水继续入渗的过程。

1. 统一向下为正的通量约定

以下将压力水头记为 h，坐标 z 和水通量 q 都以向下为正，降雨强度 P ≥ 0。忽略其他源汇项时：

[ q=K(h)\left(1-\frac{\partial h}{\partial z}\right), \qquad \frac{\partial\theta}{\partial t}=-\frac{\partial q}{\partial z}. ]

可以这样检查符号：当压力水头沿深度不变时，重力应产生向下的 q = K。

采用原代码的边界节点与算术平均导水率：

[ q_{1/2}=K_{1/2}\left(1-\frac{h_2-h_1}{\Delta z}\right). ]

因此，指定入渗通量 qT 的边界行应为：

B(1,1) = -1;
B(1,2) = 1;
d(1) = dz * (1 - qT/K12);

原代码写成了 dz * (qT/K12 + 1)，其中通量项的符号不符合上述约定。原代码输出通量时还使用了负的 dz，得到的是向上为正的通量，也应一并修改。

内部节点右端项中的重力项：

-(dt/dz) * (Kjp - Kjm)

在这里采用的约定下是正确的，不应随上边界公式一起反号。

2. 用边界约束判断是否切换

不显式计算地表储水时，限压边界可写成：

[ h_s\le h_{\max},\qquad q_s\le P,\qquad (h_{\max}-h_s)(P-q_s)=0. ]

其中 h_max = 20 cm，q_s 为实际进入土壤的通量。这意味着：

通量边界：q_s = P，求解后必须满足 h_s ≤ h_max。
水头边界：h_s = h_max，求解后必须满足 q_s ≤ P，剩余的 P − q_s 作为排走的超额水量。

不能仅根据上一次 Picard 迭代的 h_it(1) 判断边界。应先让固定边界条件下的非线性问题收敛，再检查是否满足另一项约束。这样也便于区分“非线性迭代失败”和“边界需要切换”。

这里没有引入蒸发或禁止渗出的额外限制；若土体向表面排水，q_s 可以为负。

3. 混合形式需要修正储水项

原代码使用：

d(j) = CC(j)*h_old(j) - (dt/dz)*(Kjp - Kjm);

这相当于用迭代中的 C 近似整个时间步的储水变化，收敛后一般不等于真正的 θ(h_new) − θ(h_old)。

若要离散混合形式中的含水量时间导数，可以使用 modified Picard 线性化：

[ \theta(h^{k+1})\approx\theta(h^k)+C(h^k)(h^{k+1}-h^k). ]

内部节点右端项相应改为：

d(j) = C(j)*h_it(j) ...
     - (theta_it(j) - theta_old(j)) ...
     - (dt/dz)*(Kjp - Kjm);

下面的实现采用了这项修正，也会检查离散水量平衡残差。

4. MATLAB 实现：固定边界求解，再检查约束

以下示例保留原有参数、网格、Kfun、Cfun、Qfun 和结果数组，替换原来的主时间循环，并在脚本末尾添加辅助函数。它实现的是上述无显式地表储水的限压模型。

降雨使用每个时间区间内的平均强度，避免用闭区间判断脉冲端点时多计一个时间步。对于其他 P(t)，应替换为相应的时间区间平均值。

h_max = 20;
tolH = 1e-6;       % cm
tolQ = 1e-6;       % cm/h
tolTheta = 1e-8;   % volumetric water content

qSurface = zeros(1,N);
qExcess = zeros(1,N);

for n = 1:N-1
    h_old = h(:,n);
    h_guess = h_old;

    % Mean rainfall over [t(n), t(n+1)].
    P = 0;
    for k = 1:size(events,1)
        overlap = max(0, min(t(n+1),events(k,2)) ...
                         - max(t(n),events(k,1)));
        P = P + q_pulse*overlap/dt;
    end

    mode = 'flux';
    acceptedStep = false;
    totalIter = 0;

    for activeIter = 1:4
        [h_try,qFace,ok,it] = fixedBoundaryStep( ...
            h_old,h_guess,mode,P,h_max,h_bot,dt,dz, ...
            Ks,Kfun,Cfun,Qfun,maxIter,tolH,tolQ,tolTheta);
        totalIter = totalIter + it;

        if ~ok
            error('Picard did not converge at step %d; reduce dt and retry.',n);
        end

        qIn = qFace(1);

        if strcmp(mode,'flux')
            if h_try(1) <= h_max + tolH
                acceptedStep = true;
                break;
            end
            mode = 'head';
        else
            if qIn <= P + tolQ
                acceptedStep = true;
                break;
            end
            mode = 'flux';
        end

        % Retry the SAME time step with the other boundary condition.
        h_guess = h_try;
    end

    if ~acceptedStep
        error('Boundary switching failed at step %d; reduce dt and retry.',n);
    end

    h(:,n+1) = h_try;
    theta(:,n+1) = Qfun(min(h_try,0));
    qSurface(n+1) = qIn;
    qExcess(n+1) = max(0,P-qIn);
    iterate(n+1) = totalIter;

    % Face fluxes are used in the balance equations.
    % Node values below are interpolated for output only.
    q(1,n+1) = qFace(1);
    q(M,n+1) = qFace(end);
    q(2:M-1,n+1) = 0.5*(qFace(1:end-1)+qFace(2:end));
end

辅助函数采用稀疏矩阵和 MATLAB 的线性求解器，不调用原来的 ThomasAlgorithm。

function [u,qFace,ok,it] = fixedBoundaryStep( ...
    hOld,u,mode,P,hMax,hBot,dt,dz, ...
    Ks,Kfun,Cfun,Qfun,maxIter,tolH,tolQ,tolTheta)

M = numel(hOld);
thetaOld = Qfun(min(hOld,0));
ok = false;
qFace = nan(M-1,1);

for it = 1:maxIter
    hEval = min(u,0);
    K = Kfun(hEval);
    C = Cfun(hEval);
    thetaIt = Qfun(hEval);
    K(u >= 0) = Ks;
    C(u >= 0) = 0;
    Kf = 0.5*(K(1:end-1)+K(2:end));

    if any(~isfinite(Kf)) || any(Kf <= 0)
        error('Invalid face conductivity; inspect units and constitutive functions.');
    end

    B = spalloc(M,M,3*M);
    d = zeros(M,1);

    if strcmp(mode,'flux')
        B(1,1) = -1;
        B(1,2) = 1;
        d(1) = dz*(1-P/Kf(1));
    else
        B(1,1) = 1;
        d(1) = hMax;
    end

    for j = 2:M-1
        Km = Kf(j-1);
        Kp = Kf(j);
        B(j,j-1) = -dt*Km/dz^2;
        B(j,j) = C(j)+dt*(Km+Kp)/dz^2;
        B(j,j+1) = -dt*Kp/dz^2;
        d(j) = C(j)*u(j) ...
             - (thetaIt(j)-thetaOld(j)) ...
             - dt*(Kp-Km)/dz;
    end

    B(M,M) = 1;
    d(M) = hBot;
    uNew = B\d;

    if any(~isfinite(uNew))
        return;
    end

    % Recompute nonlinear fluxes at the candidate solution.
    Knew = Kfun(min(uNew,0));
    Knew(uNew >= 0) = Ks;
    KfNew = 0.5*(Knew(1:end-1)+Knew(2:end));
    qFace = KfNew.*(1-diff(uNew)/dz);
    thetaNew = Qfun(min(uNew,0));

    % Interior mixed-form residual, in water-content units.
    residual = thetaNew(2:M-1)-thetaOld(2:M-1) ...
             + (dt/dz)*(qFace(2:end)-qFace(1:end-1));

    if strcmp(mode,'flux')
        boundaryOK = abs(qFace(1)-P) <= tolQ;
    else
        boundaryOK = abs(uNew(1)-hMax) <= tolH;
    end

    converged = max(abs(uNew-u)) <= tolH ...
             && max(abs(residual)) <= tolTheta ...
             && boundaryOK;
    u = uNew;

    if converged
        ok = true;
        return;
    end
end
end

示例中的容差可以调整，但不能保证在当前参数下收敛。代码不会自动缩小时间步；如果不收敛或边界反复切换，程序会停止。此时应从同一个旧时间层减小 dt 后重试，不能保存未收敛的解并继续计算。

5. 如果 20 cm 是真实积水深度，需要增加地表储水

令 S 为地表积水深度，R 为排走的径流强度，单位分别为 cm 和 cm/h。地表水量平衡为：

[ \frac{S^{n+1}-S^n}{\Delta t} =P^{n+1}-q_s^{n+1}-R^{n+1}. ]

存在积水时，表面压力水头应满足：

[ h_s^{n+1}=S^{n+1},\qquad 0<S^{n+1}\le20. ]

若超过 20 cm 的水立即形成径流，还需满足：

[ R^{n+1}\ge0,\qquad (20-S^{n+1})R^{n+1}=0. ]

这时有三种状态：

地表干燥：S = 0，且求解得到的 h_s ≤ 0；无蒸发条件下，降雨可以作为入渗通量。
积水未满：0 < S < 20，R = 0，需要联立储水方程和土壤方程求解 S 与 q_s。
积水达到上限：S = 20，使用 h_s = 20，并由储水方程计算 R，检查 R ≥ 0。

从有积水转为干燥的时间步，还需要计入旧时间层的存水 S^n。上面的限压代码没有包含这些状态，不能用于预测积水形成和消退的时间。

雨停后，不要仅凭 h_it(1) > 0 就设置 h_s = 0：这会把正的表面压力水头立即归零，却没有说明对应的地表水去了哪里。真实积水应随入渗逐渐减少，耗尽后再转为零降雨通量边界。

6. 使用前核对单位和离散水量平衡

alpha_v*h 必须无量纲。若 alpha_v = 0.076 的单位确实是 m⁻¹，而 h 使用 cm，则应先转换为：

alpha_v = 0.076/100; % cm^-1

若参数本来就是 0.076 cm⁻¹，则只需修正注释。这个选择会影响初始压力水头、导水率和入渗过程，不能根据代码注释确定。

原有 Qfun、Kfun 和 Cfun 的非饱和表达式也不应直接用于正压力水头。示例通过 min(h,0) 求值，并在 h ≥ 0 时使用 K = Ks、C = 0，使 θ = θ_s，沿用无额外饱和储水项的假设。

示例保留了原来的边界节点差分结构，以第一条网格边的通量近似表面通量。若要严格统计整个 80 cm 土柱的储水量，应为边界半网格建立相应的守恒离散，不能将所有节点的 θ 等权相加后当作精确水量平衡。验证时应检查时间步和网格加密后的变化，以及累计降雨、实际入渗、底部出流、储水变化和排走水量是否相互一致。