CFD有限体积法迎风格式：从一阶到TVD限制器

做CFD求解器的人大概都有过这种经历扩散项写得漂漂亮亮网格也画得规规矩矩结果一加上对流项残差曲线就开始抽风出口附近冒出一片不该存在的过冲和下冲压力场跟着抖最后整个算例发散。绝大多数情况下问题不在你的线性求解器也不在网格质量而是对流项的离散方式没选对。有限体积法里扩散项是椭圆型的、信息四面八方传播中心差分天然合适对流项是双曲型的、信息只沿特征线单向传播用中心差分就等于让下游的信息逆流而上参与决策物理上就不成立。这一篇专门讲迎风格式——它是有限体积法中对流项离散的地基也是后面所有高阶重构、限制器、通量分裂绕不开的原点。写完这一篇你应该能自己判断什么时候一阶迎风够用什么时候必须上二阶以及为什么加了限制器之后精度反而没那么高但结果更可信。1. 为什么对流项非要迎风不可1.1 中心差分在对流主导问题上的崩盘先把场景摆清楚。一维线性对流方程速度 a 恒定且为正控制方程写成有限体积法习惯的守恒形式∂u/∂t ∂(a·u)/∂x 0均匀网格单元中心记作 i单元宽度 Δx左右两个面分别是 i-1/2 和 i1/2。有限体积法最核心的一步是对每个控制体积分把体积分转成面积分一维情况下就退化成两个面通量之差d(u_i)/dt -(F_{i1/2} - F_{i-1/2}) / Δx到这里为止没有任何近似是精确的守恒写法。所有的学问都集中在面通量 F 怎么用单元平均值去算这一件事上。最朴素的想法是算术平均也就是中心差分F_{i1/2} a · (u_i u_{i1}) / 2代回去展开空间算子变成 -a(u_{i1} - u_{i-1})/(2Δx)二阶精度泰勒展开干净利落看起来很美。问题是它的修正方程里主导误差项是三阶导数色散型误差。色散误差的表现不是把解抹平而是让不同波数的分量以不同速度传播陡峭的阶跃附近会激起波长约 2Δx 的寄生振荡。这就是你看到的那些过冲——不是数值耗散不够是色散在作祟。更狠的判据来自冯·诺依曼稳定性分析。显式时间推进下中心差分的放大因子模长满足 |G|² 1 (aΔt/Δx)² sin²(θ)永远大于等于 1无条件不稳定。也就是说纯中心差分的对流项配上显式时间格式压根不可能收敛无论你把时间步缩到多小它只是振荡得慢一点。要压住它必须额外加人工黏性或者用隐式格式配上强耗散的迭代等于绕了一大圈又把耗散塞回来——那不如一开始就选一个自带合理耗散的格式。注意很多人第一次遇到这个问题时会怀疑是边界条件或者网格拉伸导致的反复折腾进出口设置。一个快速的判别方法是把对流项单独拎出来一维均匀网格、常速度、纯对流、周期边界跑一下如果仍然振荡那百分之百是格式问题跟边界和网格无关。1.2 Godunov 定理与精度—单调性的死结想要不振荡就得要求格式满足单调性初值单调的解演化过程中不能产生新的极值。这个约束比想象中苛刻。关于格式精度和单调性的关系有一个经典结论任何线性格式如果保持单调其精度最高只能到一阶。这就是绕不过去的门槛。想同时拿到二阶精度和单调性只能放弃线性这个前提——让格式系数依赖解本身也就是后来所有限制器类方法、通量分裂方法的出发点。理解了这一点前面那些为什么二阶迎风还会振荡为什么加了限制器在光滑区掉到一阶的困惑就都顺了。一阶迎风单调、稳定、不会产生新极值但它数值扩散大得离谱二阶格式精度高但在间断附近必然过冲。限制器的本质是在二者之间做一个由局部解的光滑度指标驱动的开关光滑区放行二阶重构陡峭区退回一阶或更低的重构从而在整体上既不抹平太多也不激出振荡。我在写求解器的时候习惯把这条逻辑当作选格式的第一原则先问这个算例里有没有激波、接触间断、自由面这类强梯度结构。有就必须走限制器路线别心存侥幸没有比如低速内流、纯扩散主导的传热中心差分或者高阶迎风都可以放心用甚至精度更好。1.3 迎风的物理直觉信息沿特征线传播抛开数学迎风这个词本身就是物理。速度为正流体从左边来那么界面上的状态应该由上游决定下游单元对界面状态没有话语权。这跟人站在河边观察水流是一个道理你看到的水来自上游下游发生什么不会倒着传回来亚临界流动除外那属于另一套处理逻辑。把这个直觉写成公式就是迎风格式的全部内容。速度为正时F_{i1/2} a · u_i (a 0) F_{i-1/2} a · u_{i-1} (a 0)速度变号就换成右侧单元的值。合起来写成统一形式用一个上风指示函数F_{i1/2} a⁺ · u_i a⁻ · u_{i1} a⁺ max(a, 0)a⁻ min(a, 0)这个写法在多维和非结构网格上特别有用后面会展开。它还有个隐藏的好处a⁺ 和 a⁻ 分别非负和非正保证了通量对左右两个单元状态的单调依赖这是构造更复杂格式比如通量分裂、Roe 格式时的基本要求。顺带说一句迎风这个概念在有限体积法里和有限差分里长得不太一样。有限差分里迎风是用上游节点的值代替界面值有限体积法里则是用上游单元的平均值计算界面通量。前者的守恒性要靠额外的技巧保证后者天然守恒——因为每个面通量被相邻两个单元共享一进一出严格抵消。这也是我一直推荐初学阶段直接用有限体积法框架理解迎风的原因省掉很多守恒性验证的麻烦。2. 一阶迎风最稳但最糊2.1 从面通量到代数方程的完整推导把一阶迎风代回半离散方程速度为正的情况下d(u_i)/dt -(a·u_i - a·u_{i-1}) / Δx -a·(u_i - u_{i-1}) / Δx时间上用显式欧拉记 C aΔt/Δx库朗数u_i^{n1} u_i^n - C·(u_i^n - u_{i-1}^n) (1-C)·u_i^n C·u_{i-1}^n这个式子极其直观新时刻的值是旧时刻本单元和上游单元的加权平均权重 (1-C) 和 C只要 0 ≤ C ≤ 1权重都非负加权平均不可能产生超出输入范围的值。这就是单调性的代数证明也是为什么一阶迎风无条件不振荡。C 1 时新值完全等于上游旧值物理上对应一个时间步恰好把信息推进一个网格这是纯对流问题的精确解数值扩散为零。如果是隐式后向欧拉u_i^{n1} C·(u_i^{n1} - u_{i-1}^{n1}) u_i^n整理成三对角形式用追赶法TDMA一个 O(N) 的扫描就解完而且对任意 C 都无条件稳定。工程上的稳态问题基本都是这个结构把时间导数项去掉得到a·(u_i - u_{i-1}) / Δx 0 → u_i u_{i-1}方程退化成简单的上游取值跟纯对流无源项的物理完全一致。这正是很多商业求解器在一阶迎风选项下默认干的事情。2.2 数值扩散的量级估算与它的真实杀伤力一阶迎风的问题不在稳定性在精度。做修正方程分析把 u_{i-1} 在 i 处泰勒展开u_{i-1} u_i - Δx·u_x (Δx²/2)·u_xx - (Δx³/6)·u_xxx ...代回半离散方程d(u_i)/dt -a·(Δx·u_x - (Δx²/2)·u_xx ...) / Δx -a·u_x (aΔx/2)·u_xx O(Δx²)对比原方程 ∂u/∂t a·u_x 0多出来的一项 (aΔx/2)·u_xx 是标准的扩散项形式对应的等效数值扩散系数ν_num a·Δx/2加上时间离散的贡献之后完整的表达式是ν_num (a·Δx/2)·(1 - C)这两个式子非常值得记住。第一个告诉你空间离散的数值扩散与网格尺度成正比、与速度成正比第二个告诉你时间步的影响——C 越接近 1数值扩散越小C 1 时归零。这就解释了为什么有些人明明用的是一阶迎风结果却还行他们把库朗数顶到了接近 1数值扩散被时间离散的负贡献抵消掉了一部分。但这么做风险很大一旦流场速度分布不均某些区域 C 会超过 1格式立刻失稳。数值扩散的杀伤力要用无量纲数衡量。定义网格佩克莱数Pe_Δ a·Δx / ν其中 ν 是真实的物理黏性。Pe_Δ 远大于 2 时数值扩散完全淹没物理扩散你的计算结果跟真实黏性没什么关系相当于在解一个黏性大得多的假流体。举个具体的数空气绕流a 30 m/sΔx 1 mmν 1.5e-5 m²/s那么 ν_num 30 × 0.001 / 2 0.015 m²/s是物理黏性的整整一千倍。这意味着你用一阶迎风算出来的边界层厚度、分离点位置、换热系数全都是假的。这就是为什么一阶迎风在工程上只能用来起场——给迭代一个不发散的初值收敛到一定程度后必须切到高阶格式继续算。提示判断数值扩散是否可接受一个务实的做法是看网格佩克莱数。Pe_Δ 2 时中心差分是安全且精度更好的选择Pe_Δ 在 2 到 10 之间可以用混合格式或者高阶迎风Pe_Δ 超过 10 就必须上二阶以上重构否则结果没有参考价值。2.3 一阶迎风仍然不可替代的三个场景尽管精度堪忧一阶迎风在实际代码里出现的频率依然极高原因很实在。第一它是稳定的起手式。复杂几何、多相流、燃烧这类强非线性问题初始场往往给得很粗糙直接上高阶格式几乎必然在头几十步就崩。先用一阶迎风跑几百步把场抹平残差下降到 1e-3 量级再切高阶这是行业里通行的做法。我在写多相流求解器时就是这个流程切换时机一般选在残差曲线进入稳定下降段之后太早切会振荡太晚切收敛速度被数值扩散拖得很慢。第二它是限制器退化后的保底分支。所有 TVD 限制器在检测到强梯度时限制函数 ψ 会趋于 0重构退化成 u_{i1/2} u_i也就是一阶迎风。所以一阶迎风从来没有真正离开过计算它藏在每一处激波和间断里。第三它适合做代码验证的基准。一阶迎风有解析的数值扩散系数你可以直接用它来校验自己的通量计算、边界处理、时间推进有没有写对。如果一维纯对流算出来的扩散宽度与理论值对不上说明代码里有 bug不用怀疑格式。这个技巧我用了很多次比拿复杂算例对答案高效得多。3. 高阶迎风重构精度是怎么买回来的3.1 二阶迎风的外插思路与它的过冲一阶迎风只用了上游一个单元的信息重构出来的界面值是分段常数精度自然只有一阶。想提高精度最直接的办法是多用几个上游单元做外插。速度为正时用 i 和 i-1 两个单元线性外插到面 i1/2u_{i1/2} u_i (Δx/2)·(u_i - u_{i-1})/Δx (3/2)·u_i - (1/2)·u_{i-1}这就是标准的二阶迎风Second-Order UpwindSOU。系数 3/2 和 -1/2 来自把外插点放在单元中心右侧半个网格处。可以用泰勒展开验证截断误差是 O(Δx²)。如果用三个单元做二次外插就得到三阶精度的格式但系数里会出现负值较大的项振荡风险更高工程上很少单独使用。二阶迎风的问题立刻显现系数 -1/2 是负的。重构出来的界面值可能落在 u_{i-1} 和 u_i 构成区间之外也就是产生了新的极值。在解光滑的地方这没什么误差很小但在阶跃附近u_i - u_{i-1} 很大外插会冲出去很远形成明显的过冲然后过冲在下一步又被外插放大最终演变成剧烈的振荡。这就是为什么单纯把格式换成二阶迎风残差反而可能先升后降甚至直接发散。用前面的单调性判据来检查新时刻值 u_i^{n1} (1-C)u_i C·(3/2 u_i - 1/2 u_{i-1})展开成 u_i 和 u_{i-1} 的线性组合u_{i-1} 的系数是 -C/2是负的。负权重意味着加权平均失效单调性被破坏。这跟前面一阶迎风的代数证明形成了鲜明对照也再次印证了那个结论线性二阶格式不可能单调。3.2 QUICK 格式曲率修正的巧劲QUICKQuadratic Upstream Interpolation for Convective Kinematics是另一个常被提到的三阶迎风类格式。它的思路是在迎风侧用二次多项式插值但曲率项不是用上游三点拟合而是用界面两侧各一点加一个上游点具体形式速度为正u_{i1/2} (3/8)·u_{i1} (3/4)·u_i - (1/8)·u_{i-1}系数 3/8、3/4、-1/8加起来正好等于 1保证常数场能被精确重构一致性条件。这个格式在光滑区是三阶精度比二阶迎风更准尤其在扩散主导或者剪切层这类需要高精度的场合表现好。但它的缺点同样明显系数里带负的 -1/8同样不满足单调性而且它用到了下游单元 u_{i1} 的信息在强对流问题上相当于引入了微弱的逆流影响对库朗数更敏感。我的经验是QUICK 适合用在结构化网格、流动比较平滑的内流问题比如管道内的充分发展流动、换热器通道一旦涉及激波或者自由面必须配限制器否则过冲会很难看。顺带提醒一个实现细节QUICK 在边界附近缺下游点时通常退化成二阶迎风或者一阶迎风这个降级处理如果写得不一致会在入口附近留下一条精度明显偏低的窄带看起来像边界层实际是格式降级造成的假象。第一次遇到这个问题时我查了很久还以为是入口边界条件设错了。3.3 TVD 限制器在光滑与陡峭之间切换限制器的思路是把高阶重构写成一阶迎风 受限修正量的形式。速度为正时u_{i1/2} u_i (1/2)·ψ(r)·(u_i - u_{i-1})其中 ψ(r) 是限制函数r 是相邻梯度之比r (u_i - u_{i-1}) / (u_{i1} - u_i)r 的含义是上游坡度和下游坡度谁更陡。解光滑时 r ≈ 1ψ(1) 取 1重构退化成完整的二阶迎风把 ψ1 代进去正好得到 (3/2)u_i - (1/2)u_{i-1}解在极值点附近时 r ≤ 0上下游坡度反号ψ 取 0格式退回一阶迎风不产生新极值在强梯度但单调的区域ψ 取中间值做适度限制。常用的限制函数有这么几个特性差异不小限制函数表达式光滑区精度耗散特性适用场景Minmodmax(0, min(1, r))二阶较大最稳健强激波、初值粗糙、易发散算例Van Leer(r |r|) / (1 |r|)二阶中等光滑可微通用首选收敛性好Superbeemax(0, min(2r,1), min(r,2))二阶最小压缩性强接触间断、需要锐利界面的场合Van Albada(r² r) / (r² 1)二阶中等偏大稳态问题迭代更稳Korenmax(0, min(2r, (12r)/3, 2))三阶中等追求精度又不想要 QUICK 的振荡选择哪个限制器本质是在锐度和稳健之间调档。Superbee 压缩性最强能把接触间断压得很薄但它在光滑极值点附近会把解削平出现所谓的削波现象长时间推进后波形会畸变。Minmod 最保守激波干净但光滑区耗散偏大。我自己的默认选择是 Van Leer理由是它连续可微隐式迭代时雅可比不会跳变收敛曲线漂亮只有在需要严格捕捉物质界面的多相流算例里才换 Superbee。注意限制器的 r 在分母趋零时会爆掉实现时必须加保护通常是分母加一个极小数 eps或者判断 |u_{i1} - u_i| 小于阈值时直接把 ψ 置 0。这个 bug 很隐蔽表现为偶发的 NaN只在特定的均匀场附近出现调试起来非常费劲。3.4 延迟修正让高阶格式能进隐式框架高阶迎风有个工程上的麻烦它的代数方程模板变宽了。一阶迎风的速度为正时只用到 i 和 i-1系数矩阵是三对角二维五点、三维七点求解器简单高效。二阶迎风用到 i、i-1、i-2 甚至 i1矩阵带宽增加结构化网格上还能忍非结构网格上就麻烦得多。延迟修正Deferred Correction是解决这个矛盾的标准手法。核心思想是隐式部分仍然用一阶迎风保证矩阵结构简单、对角占优、一定收敛把高阶格式与一阶格式的差值当作显式源项用上一轮迭代的值算出来加到右端项a_P·u_P Σ a_nb·u_nb S_u [F_high(u^old) - F_low(u^old)]方括号里那一项就是修正量。收敛之后 u^new ≈ u^old修正量和隐式部分的高阶效应叠加最终解等价于纯高阶格式。这个技巧的价值在于你既能享受高阶精度又能复用一阶迎风的成熟求解器还保持了矩阵的对角占优性迭代稳定性好得多。代价是收敛速度会受影响。修正量是显式处理的当它比较大的时候比如强对流、粗网格收敛会变慢有时候需要欠松弛。我的经验是欠松弛因子取 0.7 到 0.9 之间具体看残差曲线的抖动程度抖得厉害就往 0.7 靠。另外要注意延迟修正的收敛判据不能只看残差还要看修正量本身的量级——修正量没有降到足够小说明高阶效应还没完全体现这时候的收敛解其实还是偏一阶的。4. 从一维到多维迎风在真实网格上的落地4.1 面法向通量与多维迎风的写法律一维的推导很干净但真实的求解器跑的是二维三维、任意多面体网格。这时候迎风要重新表述不再是左边单元还是右边单元而是沿着面法向速度是从哪个单元指向哪个单元。对任意内部面 f设其单位法向为 n从单元 P 指向邻居 N面上的法向速度v_n V_f · nV_f 是面上的速度矢量通常由动量方程插值得到。v_n 0 表示流动从 P 流向 N上风侧是 Pv_n 0 则上风侧是 N。面质量通量ṁ_f ρ_f · v_n · A_fA_f 是面面积。一阶迎风的质量通量输运项写成F_f ṁ_f · φ_upwindφ 是被输运的标量可以是速度分量、温度、组分浓度、湍流量。这个写法跟一维完全同构只是上风的判定换成了法向通量的符号。我强烈建议在代码里把这一步封装成一个函数输入面左右单元的值和面通量输出上风值所有标量输运方程共用避免每个方程各写一套判定逻辑。多维情况下有个细节容易被忽略迎风只对法向速度做判定切向分量怎么处理严格来说多维迎风应该考虑整个速度矢量与面法向的夹角但在有限体积法的框架下通量只跟法向分量有关切向分量不参与面的质量输运所以按法向符号判定是自洽的。这一点跟通量分裂格式不同后者会把整个通量矢量按特征方向分解那是另一套体系本篇不展开。4.2 非结构网格上的上风判定与重构非结构网格四面体、多面体上一阶迎风依然是找上风单元取它的值简单可靠。麻烦的是高阶重构结构化网格上可以用 i-1、i-2 的索引直接取上游点非结构网格没有固定方向只能用梯度重构。思路是这样先给每个单元估算梯度 ∇φ_P然后用一阶泰勒展开把单元值外插到面心φ_f φ_P ∇φ_P · (x_f - x_P)x_f 是面心坐标x_P 是单元中心坐标。这个外插在光滑区是二阶精度但在梯度剧烈变化的地方同样会过冲所以还需要限制。非结构网格上的限制器通常是缩放式的不是前面那种 ψ(r) 形式。最常用的是基于邻居极值的限制算出外插到所有面上的值检查是否超出了本单元与所有邻居单元值的范围如果超了就按比例把梯度缩回去φ_lim min(1, (φ_max - φ_P) / (φ_f_max - φ_P)) ∇φ_P_limited φ_lim · ∇φ_P其中 φ_max 是 P 及所有邻居的最大值。这个做法保证外插值不产生新极值等价于一维 TVD 限制器在多维的推广。它的缺点是会把光滑极值点也削掉跟 Superbee 削波一个道理以及实现上需要一个额外的遍历所有面的循环计算量增加。Barth-Jespersen 限制器是这套方法里最经典的一个实现思路就是上述缩放。我在非结构求解器里用它的时候踩过一个坑限制因子必须对所有分量统一缩放不能每个分量各算各的。如果对速度的 x 分量缩了 0.6、y 分量缩了 0.9重构出来的速度方向就变了等于凭空引入了一个虚假的涡。这个错误在简单剪切流里看不出来一旦算有旋流场就会出现莫名其妙的旋涡结构。4.3 入口出口与内部界面的差异化处理边界上的迎风处理跟内部不一样值得单独说。入口边界如果给定的是Dirichlet值速度入口、质量流量入口那么上风侧在域外面上值直接取给定值没有迎风判定可言。这里容易出的问题是回流某些工况下入口局部出现回流此时上风侧其实在域内但代码仍然强加给定值导致质量不守恒、残差下不去。稳妥的做法是加一个回流判定检测到面法向通量为负时改用内部单元的值外插。出口边界通常用零梯度或者充分发展条件面值取上游单元的值这本身就是迎风处理比较自然。但如果出口附近有回流同样需要特殊处理否则会把出口的扰动反射回计算域形成驻波。内部界面比如多孔介质挡板、内部风扇的逻辑是两个方向都可能成为上风侧需要根据通量符号分别处理。这块的经典错误是只写了一个方向的通量另一个方向忘了镜像导致界面两侧质量不守恒表现是残差收敛到一个平台就下不去了。提示所有边界和内部界面的通量最好在代码里做一次全局质量平衡检查——把所有面的净通量求和理论上应该等于源项的总和。这个检查只要几十行代码但能抓住绝大部分通量符号写反、面积漏乘、边界漏算的问题比逐行 debug 高效得多。5. 实操从零搭一个能跑的一维迎风求解器5.1 时间推进与库朗数的选择一节内容力求把参数选择的逻辑讲透。显式一阶迎风的稳定条件前面已经推过0 ≤ C ≤ 1。但实际操作中没人会贴着 C1 跑原因有两个一是非线性问题里速度场在变局部 C 可能超过 1二是多维情况下稳定性条件更严格。二维显式对流的稳定条件近似是C_x C_y ≤ 1三维则是三个方向库朗数之和不超过 1。这个求和约束在非结构网格上写成对面通量绝对值求和的形式Δt ≤ V_P / Σ_f |v_n,f| · A_fV_P 是单元体积。这个式子对所有网格类型通用实现起来只需要遍历一次面把绝对通量加起来取全局最小值再乘一个安全系数。安全系数我一般取 0.5 到 0.8即实际时间步是理论极限的一半到八成。为什么留这么多余量因为理论极限是针对线性常系数问题的真实问题里有源项、有非线性、有边界反射贴着极限跑很容易在某个局部网格单元上炸掉。而且显式格式一旦局部失稳扰动会迅速扩散到全场抢救成本很高。我的习惯是先用 0.5 跑通确认稳定后再逐步往上加观察残差曲线的抖动情况抖了就退回去。隐式格式理论上无条件稳定但实际受限于非线性迭代的收敛性。稳态问题常用的伪时间步方法局部时间步按 CFL 数给Δt_local CFL · V_P / Σ_f |v_n,f| · A_fCFL 从 5 起步收敛顺利可以逐步加到 50 甚至 100遇到残差上翘就降回来。这个自适应策略比固定 CFL 省时间尤其是从粗糙初场起步的时候。5.2 一个最小可跑的实现骨架下面这段代码把前面所有讨论落到纸面上一维、均匀网格、周期边界、纯对流四个格式可切换一阶迎风、二阶迎风、QUICK、带 Van Leer 限制器的二阶迎风。可以直接跑跑完对比扩散宽度和过冲幅度。import numpy as np def limiter(r, kindvanleer): eps 1e-12 if kind minmod: return np.maximum(0.0, np.minimum(1.0, r)) if kind vanleer: return (r np.abs(r)) / (1.0 np.abs(r) eps) if kind superbee: return np.maximum(0.0, np.maximum(np.minimum(2*r, 1), np.minimum(r, 2))) return np.ones_like(r) def rhs(u, a, dx, schemeupwind1): n len(u) F np.zeros(n 1) # 面 i1/2 位于 u[i] 与 u[i1] 之间索引 face k 对应单元 k-1 与 k for k in range(n 1): iL (k - 1) % n iR k % n iLL (k - 2) % n iRR (k 1) % n if a 0: if scheme upwind1: F[k] a * u[iL] elif scheme upwind2: F[k] a * (1.5*u[iL] - 0.5*u[iLL]) elif scheme quick: F[k] a * (0.75*u[iL] 0.375*u[iR] - 0.125*u[iLL]) elif scheme tvd: d_up u[iL] - u[iLL] d_dn u[iR] - u[iL] if abs(d_dn) 1e-12: r 0.0 else: r d_up / d_dn psi limiter(np.array([r]), vanleer)[0] u_face u[iL] 0.5 * psi * d_up F[k] a * u_face return -(F[1:] - F[:-1]) / dx def run(scheme, n200, cfl0.6, t_end0.5, a1.0): L 1.0 dx L / n dt cfl * dx / a x (np.arange(n) 0.5) * dx # 初始条件一个方波 u np.where((x 0.2) (x 0.4), 1.0, 0.0) t 0.0 while t t_end: u u dt * rhs(u, a, dx, scheme) t dt return x, u跑完之后对比会发现一阶迎风的方波顶多剩下一个低矮的小鼓包扩散宽度约为 sqrt(2·ν_num·t)代进去 v_num a·dx/2·(1-C)误差在 15% 以内就算对二阶迎风的方波边缘锐利但两侧有明显过冲幅值随网格加密应该逐渐减小QUICK 类似但过冲位置略有不同TVD 版本的边缘比一阶锐利同时没有过冲代价是波形顶部的平台区被略微削低。这个算例的验证价值在于如果你的一阶迎风跑出来的扩散宽度跟理论值差很多说明时间推进、面通量索引、边界环绕这几处有 bug。我在写新求解器时第一个跑通的永远是这个算例第二个是纯扩散的解析解两个都对上了才敢上复杂几何。5.3 网格分辨率的经验取值迎风格式的精度是网格相关的实际工程里怎么定网格密度靠的不是公式而是经验加收敛性检查。对一阶迎风网格佩克莱数要控制在 2 以内才不至于把物理扩散完全淹没这是硬指标。反推网格尺度Δx ≤ 2ν / a空气 30 m/s、ν 1.5e-5 的话Δx ≤ 1 微米完全不现实。所以一阶迎风在大佩克莱数问题里根本没有实用价值只能起场。对二阶迎风或者 TVD 格式判断标准换成网格无关性验证。做法是同一算例用三套网格粗、中、细网格量各差约 4 倍算同一个关注的物理量比如阻力系数、最大温度、分离点位置看相邻两套之间的差异是否小于某个阈值常用 2% 到 5%。差异收敛了才能说网格够了。这个过程很花时间但比起拿着一张精美但错误的云图去汇报还是划算的。还有一个常被忽略的点迎风格式的精度损失在边界附近和界面附近最明显因为那里有降级处理。所以网格加密时边界层区域和相界面区域的加密要优先于远场效果差好几倍。6. 常见问题与排查实录6.1 问题速查表现象可能原因排查与解决残差下降后突然上翘局部库朗数超限格式失稳打印全局最大库朗数降时间步或加密网格阶跃附近出现过冲/下冲用了无限制的高阶格式切换到带限制器的格式或改用 Minmod解被抹得过度平滑一阶迎风数值扩散过大提高格式阶数检查网格佩克莱数收敛到一个非零残差平台内部界面或边界通量不守恒做全局质量平衡检查核对通量符号偶发 NaN 且只在某处出现限制器分母趋零、梯度计算出错分母加 eps检查梯度重构的邻居列表波形长时间推进后畸变压缩性限制器的削波效应换 Van Leer 或降到 Minmod减小时间步入口附近有条精度偏低的窄带边界处高阶重构降级检查边界处理是否用了与内部一致的格式稳态解在不同网格上差异很大数值扩散主导网格无关性未达到做三套网格验证优先加密高梯度区6.2 几个踩过的坑和对应的解法限制器写得对但用错了位置是我早期最常犯的错误。有些实现会先把所有面通量算出来存进数组再统一更新单元这没问题但如果限制器函数在算每个面的时候用了最新的、已经部分更新的单元值比如原地更新 u[i] 然后继续算下一个面就会出现方向依赖的结果一维下表现为波形不对称二维下表现为结果跟网格遍历顺序有关。解决办法很简单先完整算一遍所有面通量存到数组再统一更新别偷懒省内存。另一个坑是延迟修正收敛判据。前面提过修正量没降到足够小解其实是偏一阶的。我吃过一次亏残差已经降到 1e-6看着挺收敛但出口温度跟实验数据差了三度回头一查发现修正量还有 1e-3 量级也就是说高阶效应根本没体现出来。后来把判据改成两个都要满足残差小于阈值且修正量的最大绝对值小于某个物理量的特征量乘以 1e-4。加了这条之后同类问题再没出现过。还有一点关于稳态和瞬态的格式选择。稳态问题用高阶迎风配延迟修正收敛后精度是实打实的瞬态问题如果时间步比较大空间离散的阶数再高也没用因为时间离散的误差会主导。这时候要么改用二阶时间格式Crank-Nicolson 或者二阶后向差分要么就把时间步压下来做网格无关性验证。我见过不少人在瞬态算例里死磕空间格式调了半天限制器其实瓶颈在时间精度上。最后分享一个实用的小技巧验证迎风格式实现是否正确除了方波算例还可以跑一个高斯脉冲的纯对流。高斯脉冲有解析解扩散后仍然是高斯标准差按 sqrt(σ0² 2ν_num·t) 增长。这个算例比方波更敏感因为它的波形处处光滑一阶和二阶的差异一眼可见限制器在光滑区的行为也能看出来。我一般用它做第二个验证两个算例都过了才认为对流项的实现可以信了。