边界条件选型:结果诊断与可信度验证

用入流特征数、质量与能量闭合、域长外推三类证据判断边界条件是否欠定或被污染,给出可核对的判据阈值、OpenFOAM 诊断配置与失败模式对照表。

CFD结果诊断与可信度验证边界条件与初始化边界条件选型入流特征数质量闭合专题

边界条件选型:结果诊断与可信度验证

边界条件选错极少直接报错。它更常见的表现是残差在某个平台不再下降、目标量随计算域尺寸单调漂移,或者出口截面上出现本该流出的流体反向灌入。本文给出一套只用求解器自带输出即可完成的诊断顺序:先用入流特征数判断边界条件是否欠定,再用质量与能量闭合、域长外推定位污染来源,最后与独立基准对照,给出可引用的误差条。

残差平台抬高不等于离散格式有问题

求解器残差停在一个高位平台,第一嫌疑往往被归给对流格式或松弛因子,但边界欠定的典型指纹是连续性型残差(不可压缩流的 pp 方程、可压缩流的 ρ\rho 方程)单独抬高,而动量残差已经很低。原因是边界上没有足够的约束把进入域内的质量流"接住":若入口给定质量流量,出口又给定总压,则流量由两侧共同决定,方程组的定解性被破坏,压力修正会在边界附近反复补偿。

判据是质量收支。对任意稳态不可压缩算例,全域质量不平衡率应满足

εm˙=∣m˙in−m˙out∣m˙in×100%\varepsilon_{\dot m} = \frac{\left| \dot m_{in} - \dot m_{out} \right|}{\dot m_{in}} \times 100\%

工程上 εm˙<0.1%\varepsilon_{\dot m} < 0.1\% 视为闭合良好,0.1%∼1%0.1\% \sim 1\% 需要解释来源(通常是松弛未收敛或时间平均),超过 1%1\% 说明边界组合存在结构性冲突,继续迭代不会改善。举例:内径 D=0.05 mD = 0.05\ \mathrm{m} 的圆管,水在 U=2 m/sU = 2\ \mathrm{m/s} 下流动,ρ=998 kg/m3\rho = 998\ \mathrm{kg/m^3},则截面面积 A=πD2/4=1.9635×10−3 m2A = \pi D^2/4 = 1.9635 \times 10^{-3}\ \mathrm{m^2},入口质量流量

m˙=ρAU=998×1.9635×10−3×2=3.92 kg/s\dot m = \rho A U = 998 \times 1.9635\times 10^{-3} \times 2 = 3.92\ \mathrm{kg/s}

若出口报告 3.94 kg/s3.94\ \mathrm{kg/s},则 εm˙=0.02/3.92=0.51%\varepsilon_{\dot m} = 0.02/3.92 = 0.51\%,落在"需要解释"的区间,此时应先把松弛因子收敛到机器精度再看,而不是先改格式。该工况的雷诺数为 Re=ρUD/μ=998×2×0.05/10−3=9.98×104Re = \rho U D/\mu = 998 \times 2 \times 0.05 / 10^{-3} = 9.98\times 10^4,属于充分湍流,任何边界污染都会在数十倍管径后才衰减。

入流特征数决定边界上必须给定几个量

双曲型系统的边界条件个数由指向域内的特征线数量决定,这是选型的第一性依据。设边界外法向为 nn,局部法向马赫数为 Mn=un/aM_n = u_n / a,则亚声速入流需要给定 4 个量(三维速度三分量与一个热力学量),亚声速出流只允许给定 1 个量(压力),超声速入流给定全部 5 个量,超声速出流一个都不给:

nin={4,Mn<1 入流1,Mn<1 出流5,Mn>1 入流0,Mn>1 出流n_{in} = \begin{cases} 4, & M_n < 1 \ \text{入流} \\ 1, & M_n < 1 \ \text{出流} \\ 5, & M_n > 1 \ \text{入流} \\ 0, & M_n > 1 \ \text{出流} \end{cases}

把这条规则用在可压缩喷管上很容易暴露错误:出口局部 Mn=0.85M_n = 0.85 时若同时指定静压与速度,就多给了一个条件,求解器只能丢掉其中一个;反之出口 Mn=1.6M_n = 1.6 时若指定压力,等于在超声速段强行反射信息,压力场会出现非物理的下游扰动。不可压缩流是 Mn→0M_n \to 0 的极限,出流只保留"压力水平"这一个自由度,这正是压力出口的物理含义。

域长外推:把人为边界推离关注区

任何有限计算域都在某处施加了人为条件。判断它是否污染目标量,不能靠肉眼看流场,而要靠目标量对域长的敏感性:

SL=q(L2)−q(L1)q(L1)×100%S_L = \frac{q(L_2) - q(L_1)}{q(L_1)} \times 100\%

其中 LL 以特征尺寸 DD 为单位,qq 是关注量(压降、升力系数、出口平均温度等)。工程判据:域长加倍时 ∣SL∣<1%|S_L| < 1\% 才认为出口位置合格。某弯管算例在 L=20DL = 20D 得到压降 1250 Pa1250\ \mathrm{Pa},L=40DL = 40D 得到 1235 Pa1235\ \mathrm{Pa},L=80DL = 80D 得到 1229 Pa1229\ \mathrm{Pa}:第一次加倍 $S_L = (1235-1250)/1250 = -1.2%,第二次加倍,第二次加倍S_L = (1229-1235)/1235 = -0.49%。敏感性在缩小但尚未低于。敏感性在缩小但尚未低于1%,说明,说明40D仍不足,报告里应写明该量在仍不足,报告里应写明该量在80D下的不确定度约为下的不确定度约为0.5%$,而不是声称"已收敛"。

对钝体绕流,尾迹恢复距离的经验下限为 Lwake≈10D∼20DL_{wake} \approx 10D \sim 20D,层流管流入口段为 Ldev≈0.06 Re DL_{dev} \approx 0.06\,Re\,D;本例 Re=9.98×104Re = 9.98\times 10^4,代入得 Ldev≈0.06×9.98×104×0.05=299 mL_{dev} \approx 0.06 \times 9.98\times 10^4 \times 0.05 = 299\ \mathrm{m},即约 6000D6000D,这就是为什么工程上从不用层流入口段公式去估湍流管流,而直接取 L≈10DL \approx 10D 量级并在直段中给定充分发展剖面。

用函数对象把闭合检查自动化

质量收支与壁面通量应作为每次运行的固定输出,而不是事后手工统计。OpenFOAM 中在 system/controlDict 末尾追加:

functions
{
    flowIn
    {
        type            surfaceFieldValue;
        libs            (fieldFunctionObjects);
        fields          (phi);
        regionType      patch;
        name            inlet;
        operation       sum;
        writeFields     false;
    }
    flowOut
    {
        type            surfaceFieldValue;
        libs            (fieldFunctionObjects);
        fields          (phi);
        regionType      patch;
        name            outlet;
        operation       sum;
        writeFields     false;
    }
    wallFlux
    {
        type            wallHeatFlux;
        libs            (fieldFunctionObjects);
        patches         (wall);
    }
}

运行时用 postProcess -func flowIn -latestTime 与 postProcess -func flowOut -latestTime 分别取值,两者之差即为 εm˙\varepsilon_{\dot m} 的分子。对共轭传热算例再加 wallHeatFlux,把界面两侧的通量同时打出,若两侧符号相同或量级差 10 倍以上,说明界面拓扑或热物性接错了。

失败模式对照

现象 根因 判定试验
压力方程残差卡在 10−310^{-3} 不再下降 入口给流量、出口给总压,流量被双重约束 把出口改为固定静压,观察连续性残差是否立刻下降一个量级
压降随域长单调减小且未收敛 出口距弯头或扩张段不足 出口后移 10D10D 重算,∣SL∣>1%\lvert S_L \rvert > 1\% 即判未收敛
出口截面出现逆向质量流 回流未给定温度/湍流量,求解器用默认值填充 输出出口 ϕ<0\phi < 0 的面积占比,稳态 > 2% 或瞬态峰值 > 10% 需处理
对称面出现法向速度 对称面被误设为壁面或滑移壁面 输出对称面 unu_n 峰值,应小于参考速度的 0.1%0.1\%
全域压力随时间缓慢漂移 参考压力单元落在回流区或大梯度区 更换 referenceCell 到静止区,检查域平均压力是否变平
壁面热流符号与预期相反 界面两侧法向定义不一致 分别输出两侧 wallHeatFlux,检查是否同号

与独立基准对照

闭合与敏感性只证明"自洽",不证明"正确"。最后一步必须引入外部证据:解析解、基准实验或文献关联式。层流圆管可用 Hagen-Poiseuille 关系 f=64/Ref = 64/Re 校核摩擦系数,平板边界层可用 Blasius 解校核 CfC_f,绕流可用实验升阻力系数校核。对照时应报告相对误差与不确定度来源,例如"L=80DL = 80D 下压降 1229 Pa1229\ \mathrm{Pa},域长不确定度 0.5%0.5\%,与文献关联式偏差 2.3%2.3\%",而不是给一个没有误差条的单一数字。若偏差随网格加密单调减小,属于离散误差;若随域长变化而变,属于边界污染;若两者都不变,则要怀疑湍流模型或物性输入。

参考文献

  1. Hirsch C., Numerical Computation of Internal and External Flows, 2nd ed., Butterworth-Heinemann, 2007.
  2. Ferziger J.H., Perić M., Street R.L., Computational Methods for Fluid Dynamics, 4th ed., Springer, 2020.
  3. Roache P.J., Verification and Validation in Computational Science and Engineering, Hermosa Publishers, 1998.
  4. AIAA, Guide for the Verification and Validation of Computational Fluid Dynamics Simulations, AIAA G-077-1998, 1998.
  5. The OpenFOAM Foundation, OpenFOAM User Guide, Release v2312, 2023.