📚 课程目录第4章质量守恒数值校验实战

4. 质量守恒数值校验实战

连续性方程是质量守恒在流体力学中的数学表达。有限体积法（FVM）把计算域切成若干控制体，对每个控制体写出“净流出质量通量 = 内部质量变化率”的积分形式。数值校验的核心，就是检查离散后的通量账本是否平账。本节围绕不可压缩与可压缩两类流动，给出可落地的校验方案、脚本实现与边界通量误差分析方法。

一、先明确符号与参照约定

设控制体体积为 V，第 f 个面面积为 A_f，外法向单位向量为 n_f，面心速度为 u_f，密度为 ρ_f。定义流出为正：面质量通量 ṁ_f = ρ_f (u_f · n_f) A_f。于是离散连续性方程为：

Σ_f ṁ_f  +  d(ρV)/dt  = 0        （流出为正）

不可压缩流动中 ρ 为常数，方程退化为 Σ_f (u_f · n_f) A_f = 0，即体积通量守恒；可压缩流动则必须保留密度与时间项。全文与图示均采用“流出为正”，若某软件采用流入为正，需整体取负，切勿混用。

二、不可压缩流动：体积通量闭合校验

以二维顶盖驱动方腔为例，取 64×64 均匀网格，顶盖速度 U = 1 m/s，其余三壁无滑移。校验分三层：

逐面通量核算：对每个内部面，用两侧单元速度插值得到 u_f，计算 Q_f = (u_f · n_f) A_f。内部面在相邻两单元中符号相反，累加时应严格抵消。
逐单元残差：R_i = Σ_f Q_f，理论上应为 0。实际因插值不满足离散无散度条件，残差非零，需统计其量级。
全域闭合：所有单元残差之和应等于边界净通量，即 Σ_i R_i = Σ_boundary Q_f。方腔四壁均为固壁，边界净通量应为 0，故全域残差和应趋近机器精度。

一个常见的坑是：压力修正方程求解后速度场满足离散连续性，但若在输出时对速度做了额外插值或滤波，再重新计算通量就会出现“假残差”。因此校验脚本必须直接读取求解器内部用于通量计算的面速度，而不是后处理插值结果。

三、可压缩流动：含密度与时间项的通量校验

以一维Sod激波管为例，长度 L = 1 m，200 个单元，初始间断位于 x = 0.5：左态 (ρ, u, p) = (1.0, 0, 1.0)，右态 (0.125, 0, 0.1)，绝热指数 γ = 1.4，CFL 取 0.4。可压缩校验的关键是时间项不能忽略：

R_i = Σ_f ṁ_f  +  (ρ_i^{n+1} - ρ_i^{n}) V_i / Δt

两端均为固壁，边界质量通量为 0，因此全域 Σ_i R_i 应等于全域质量变化率，而全域质量本身应守恒（总质量不变）。校验时同时输出两个量：逐单元残差的最大绝对值（衡量局部平账程度），以及全域总质量随时间的漂移（衡量全局守恒性）。Sod 问题在激波与接触间断处残差较大属正常现象，因为那里梯度陡、重构误差大；但全域总质量漂移应保持在 10⁻¹² 量级（双精度）。

四、边界通量误差分析

边界是误差最容易积累的地方。对出口边界，若采用零梯度外推，面密度与面速度都取自内部单元，通量误差主要来自外推假设；对入口边界，给定总压总温时需反算密度与速度，反算过程引入的误差会直接进入通量账本。建议单独统计边界通量误差：

e_b = | Σ_boundary ṁ_f  -  (M^{n+1} - M^{n}) / Δt | / max(|Σ_boundary ṁ_f|, ε)

其中 M 为全域质量。该相对误差在收敛良好的算例中应低于 10⁻⁶；若明显偏大，优先检查边界面的密度取值是否与内部重构一致、时间推进是否与通量计算使用同一时刻的状态。

五、校验脚本骨架

下面给出一个与求解器解耦的 Python 校验脚本骨架，输入为面通量数组与单元密度、体积、时间步长。它同时适用于两类流动：不可压缩时把密度视为常数、时间项置零即可。

import numpy as np

def check_mass_balance(face_flux, cell_vol, rho_old, rho_new, dt, boundary_mask):
    # face_flux: 每个面的 ṁ_f，流出为正
    # 内部面在两侧单元中符号相反，需按单元-面关联表累加
    residual = np.zeros_like(cell_vol)
    for f, (cL, cR) in enumerate(face_cells):
        residual[cL] += face_flux[f]
        residual[cR] -= face_flux[f]
    residual += (rho_new - rho_old) * cell_vol / dt

    bflux = face_flux[boundary_mask].sum()
    dM_dt = ((rho_new - rho_old) * cell_vol).sum() / dt
    e_b = abs(bflux - dM_dt) / max(abs(bflux), 1e-30)

    return {
        "max_cell_residual": np.abs(residual).max(),
        "global_residual": residual.sum(),
        "boundary_flux": bflux,
        "mass_drift": dM_dt,
        "relative_boundary_error": e_b,
    }

脚本输出五个量：逐单元残差最大值、全域残差和、边界净通量、全域质量变化率、边界相对误差。对不可压缩算例，rho_old = rho_new = 1，时间项自动为零，global_residual 应等于 boundary_flux，二者之差即离散闭合误差。

六、结果判读与常见问题

判读时不要只看单一数字。若逐单元残差最大值很大但全域残差和很小，说明误差在内部相互抵消，属于局部重构问题，通常出现在大梯度区；若全域残差和随时间单调增长，说明边界通量与内部通量使用了不一致的状态量，属于实现层面的错误。可压缩算例中，若总质量漂移超过 10⁻¹⁰，应检查通量计算与时间推进是否采用同一时间层，以及边界 ghost 单元的密度是否每步更新。

建议把校验脚本挂到每个时间步之后运行，只输出五个标量，开销可忽略；一旦 relative_boundary_error 超过阈值就中断并保存现场，比事后翻日志高效得多。

最后强调一点：质量守恒校验不是“跑通即可”，而是要建立可比较的误差基线。同一求解器、同一网格下，先记录一个已知良好算例的残差量级，后续改动格式、边界条件或时间推进时，用同一脚本复测，误差量级是否漂移就是最直接的回归信号。