企业官网建设流程全解析

做CFD求解器的人大概都有过这种经历：扩散项写得漂漂亮亮，网格也画得规规矩矩，结果一加上对流项，残差曲线就开始抽风，出口附近冒出一片不该存在的过冲和下冲，压力场跟着抖，最后整个算例发散。绝大多数情况下，问题不在你的线性求解器，也不在网格质量，而是对流项的离散方式没选对。有限体积法里，扩散项是椭圆型的、信息四面八方传播，中心差分天然合适；对流项是双曲型的、信息只沿特征线单向传播，用中心差分就等于让下游的信息逆流而上参与决策，物理上就不成立。这一篇专门讲迎风格式——它是有限体积法中对流项离散的地基，也是后面所有高阶重构、限制器、通量分裂绕不开的原点。写完这一篇，你应该能自己判断：什么时候一阶迎风够用，什么时候必须上二阶，以及为什么加了限制器之后精度反而没那么高但结果更可信。

1. 为什么对流项非要"迎风"不可

1.1 中心差分在对流主导问题上的崩盘

先把场景摆清楚。一维线性对流方程，速度 a 恒定且为正，控制方程写成有限体积法习惯的守恒形式：

∂u/∂t + ∂(a·u)/∂x = 0

均匀网格，单元中心记作 i，单元宽度 Δx，左右两个面分别是 i-1/2 和 i+1/2。有限体积法最核心的一步是对每个控制体积分，把体积分转成面积分，一维情况下就退化成两个面通量之差：

d(u_i)/dt = -(F_{i+1/2} - F_{i-1/2}) / Δx

到这里为止没有任何近似，是精确的守恒写法。所有的学问都集中在"面通量 F 怎么用单元平均值去算"这一件事上。最朴素的想法是算术平均，也就是中心差分：

F_{i+1/2} = a · (u_i + u_{i+1}) / 2

代回去展开，空间算子变成 -a(u_{i+1} - u_{i-1})/(2Δx)，二阶精度，泰勒展开干净利落，看起来很美。问题是它的修正方程里，主导误差项是三阶导数，色散型误差。色散误差的表现不是把解抹平，而是让不同波数的分量以不同速度传播，陡峭的阶跃附近会激起波长约 2Δx 的寄生振荡。这就是你看到的那些过冲——不是数值耗散不够，是色散在作祟。

更狠的判据来自冯·诺依曼稳定性分析。显式时间推进下，中心差分的放大因子模长满足 |G|² = 1 + (aΔt/Δx)² sin²(θ)，永远大于等于 1，无条件不稳定。也就是说纯中心差分的对流项配上显式时间格式，压根不可能收敛，无论你把时间步缩到多小，它只是振荡得慢一点。要压住它，必须额外加人工黏性或者用隐式格式配上强耗散的迭代，等于绕了一大圈又把耗散塞回来——那不如一开始就选一个自带合理耗散的格式。

注意：很多人第一次遇到这个问题时会怀疑是边界条件或者网格拉伸导致的，反复折腾进出口设置。一个快速的判别方法是把对流项单独拎出来，一维均匀网格、常速度、纯对流、周期边界跑一下，如果仍然振荡，那百分之百是格式问题，跟边界和网格无关。

1.2 Godunov 定理与"精度—单调性"的死结

想要不振荡，就得要求格式满足单调性：初值单调的解，演化过程中不能产生新的极值。这个约束比想象中苛刻。关于格式精度和单调性的关系有一个经典结论：任何线性格式，如果保持单调，其精度最高只能到一阶。这就是绕不过去的门槛。想同时拿到二阶精度和单调性，只能放弃"线性"这个前提——让格式系数依赖解本身，也就是后来所有限制器类方法、通量分裂方法的出发点。

理解了这一点，前面那些"为什么二阶迎风还会振荡""为什么加了限制器在光滑区掉到一阶"的困惑就都顺了。一阶迎风单调、稳定、不会产生新极值，但它数值扩散大得离谱；二阶格式精度高，但在间断附近必然过冲。限制器的本质是在二者之间做一个由局部解的光滑度指标驱动的开关：光滑区放行二阶重构，陡峭区退回一阶或更低的重构，从而在整体上既不抹平太多，也不激出振荡。

我在写求解器的时候习惯把这条逻辑当作选格式的第一原则：先问这个算例里有没有激波、接触间断、自由面这类强梯度结构。有，就必须走限制器路线，别心存侥幸；没有，比如低速内流、纯扩散主导的传热，中心差分或者高阶迎风都可以放心用，甚至精度更好。

1.3 迎风的物理直觉：信息沿特征线传播

抛开数学，"迎风"这个词本身就是物理。速度为正，流体从左边来，那么界面上的状态应该由上游决定，下游单元对界面状态没有话语权。这跟人站在河边观察水流是一个道理：你看到的水来自上游，下游发生什么不会倒着传回来（亚临界流动除外，那属于另一套处理逻辑）。

把这个直觉写成公式，就是迎风格式的全部内容。速度为正时：

F_{i+1/2} = a · u_i (a > 0) F_{i-1/2} = a · u_{i-1} (a > 0)

速度变号就换成右侧单元的值。合起来写成统一形式，用一个上风指示函数：

F_{i+1/2} = a⁺ · u_i + a⁻ · u_{i+1} 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^{n+1} = 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^{n+1} + C·(u_i^{n+1} - u_{i-1}^{n+1}) = 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_{i+1/2} = u_i，也就是一阶迎风。所以一阶迎风从来没有真正离开过计算，它藏在每一处激波和间断里。

第三，它适合做代码验证的基准。一阶迎风有解析的数值扩散系数，你可以直接用它来校验自己的通量计算、边界处理、时间推进有没有写对。如果一维纯对流算出来的扩散宽度与理论值对不上，说明代码里有 bug，不用怀疑格式。这个技巧我用了很多次，比拿复杂算例对答案高效得多。

3. 高阶迎风重构：精度是怎么"买"回来的

3.1 二阶迎风的外插思路与它的过冲

一阶迎风只用了上游一个单元的信息，重构出来的界面值是分段常数，精度自然只有一阶。想提高精度，最直接的办法是多用几个上游单元做外插。速度为正时，用 i 和 i-1 两个单元线性外插到面 i+1/2：

u_{i+1/2} = u_i + (Δx/2)·(u_i - u_{i-1})/Δx = (3/2)·u_i - (1/2)·u_{i-1}

这就是标准的二阶迎风（Second-Order Upwind，SOU）。系数 3/2 和 -1/2 来自把外插点放在单元中心右侧半个网格处。可以用泰勒展开验证，截断误差是 O(Δx²)。如果用三个单元做二次外插，就得到三阶精度的格式，但系数里会出现负值较大的项，振荡风险更高，工程上很少单独使用。

二阶迎风的问题立刻显现：系数 -1/2 是负的。重构出来的界面值可能落在 u_{i-1} 和 u_i 构成区间之外，也就是产生了新的极值。在解光滑的地方这没什么，误差很小；但在阶跃附近，u_i - u_{i-1} 很大，外插会冲出去很远，形成明显的过冲，然后过冲在下一步又被外插放大，最终演变成剧烈的振荡。这就是为什么单纯把格式换成二阶迎风，残差反而可能先升后降，甚至直接发散。

用前面的单调性判据来检查：新时刻值 u_i^{n+1} = (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 格式：曲率修正的巧劲

QUICK（Quadratic Upstream Interpolation for Convective Kinematics）是另一个常被提到的三阶迎风类格式。它的思路是在迎风侧用二次多项式插值，但曲率项不是用上游三点拟合，而是用界面两侧各一点加一个上游点，具体形式（速度为正）：

u_{i+1/2} = (3/8)·u_{i+1} + (3/4)·u_i - (1/8)·u_{i-1}

系数 3/8、3/4、-1/8，加起来正好等于 1，保证常数场能被精确重构（一致性条件）。这个格式在光滑区是三阶精度，比二阶迎风更准，尤其在扩散主导或者剪切层这类需要高精度的场合表现好。

但它的缺点同样明显：系数里带负的 -1/8，同样不满足单调性；而且它用到了下游单元 u_{i+1} 的信息，在强对流问题上相当于引入了微弱的"逆流"影响，对库朗数更敏感。我的经验是，QUICK 适合用在结构化网格、流动比较平滑的内流问题，比如管道内的充分发展流动、换热器通道；一旦涉及激波或者自由面，必须配限制器，否则过冲会很难看。

顺带提醒一个实现细节：QUICK 在边界附近缺下游点时，通常退化成二阶迎风或者一阶迎风，这个降级处理如果写得不一致，会在入口附近留下一条精度明显偏低的窄带，看起来像"边界层"，实际是格式降级造成的假象。第一次遇到这个问题时我查了很久，还以为是入口边界条件设错了。

3.3 TVD 限制器：在光滑与陡峭之间切换

限制器的思路是把高阶重构写成"一阶迎风 + 受限修正量"的形式。速度为正时：

u_{i+1/2} = u_i + (1/2)·ψ(r)·(u_i - u_{i-1})

其中 ψ(r) 是限制函数，r 是相邻梯度之比：

r = (u_i - u_{i-1}) / (u_{i+1} - u_i)

r 的含义是"上游坡度和下游坡度谁更陡"。解光滑时 r ≈ 1，ψ(1) 取 1，重构退化成完整的二阶迎风（把 ψ=1 代进去正好得到 (3/2)u_i - (1/2)u_{i-1}）；解在极值点附近时 r ≤ 0（上下游坡度反号），ψ 取 0，格式退回一阶迎风，不产生新极值；在强梯度但单调的区域，ψ 取中间值，做适度限制。

常用的限制函数有这么几个，特性差异不小：

选择哪个限制器，本质是在"锐度"和"稳健"之间调档。Superbee 压缩性最强，能把接触间断压得很薄，但它在光滑极值点附近会把解"削平"，出现所谓的削波现象，长时间推进后波形会畸变。Minmod 最保守，激波干净但光滑区耗散偏大。我自己的默认选择是 Van Leer，理由是它连续可微，隐式迭代时雅可比不会跳变，收敛曲线漂亮；只有在需要严格捕捉物质界面的多相流算例里才换 Superbee。

注意：限制器的 r 在分母趋零时会爆掉，实现时必须加保护，通常是分母加一个极小数 eps，或者判断 |u_{i+1} - u_i| 小于阈值时直接把 ψ 置 0。这个 bug 很隐蔽，表现为偶发的 NaN，只在特定的均匀场附近出现，调试起来非常费劲。

3.4 延迟修正：让高阶格式能进隐式框架

高阶迎风有个工程上的麻烦：它的代数方程模板变宽了。一阶迎风的速度为正时只用到 i 和 i-1，系数矩阵是三对角（二维五点、三维七点），求解器简单高效。二阶迎风用到 i、i-1、i-2 甚至 i+1，矩阵带宽增加，结构化网格上还能忍，非结构网格上就麻烦得多。

延迟修正（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 · n

V_f 是面上的速度矢量，通常由动量方程插值得到。v_n > 0 表示流动从 P 流向 N，上风侧是 P；v_n < 0 则上风侧是 N。面质量通量：

ṁ_f = ρ_f · v_n · A_f

A_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。但实际操作中没人会贴着 C=1 跑，原因有两个：一是非线性问题里速度场在变，局部 C 可能超过 1；二是多维情况下稳定性条件更严格。

二维显式对流的稳定条件近似是：

C_x + C_y ≤ 1

三维则是三个方向库朗数之和不超过 1。这个"求和"约束在非结构网格上写成对面通量绝对值求和的形式：

Δt ≤ V_P / Σ_f |v_n,f| · A_f

V_P 是单元体积。这个式子对所有网格类型通用，实现起来只需要遍历一次面，把绝对通量加起来，取全局最小值再乘一个安全系数。安全系数我一般取 0.5 到 0.8，即实际时间步是理论极限的一半到八成。

为什么留这么多余量？因为理论极限是针对线性常系数问题的，真实问题里有源项、有非线性、有边界反射，贴着极限跑很容易在某个局部网格单元上炸掉。而且显式格式一旦局部失稳，扰动会迅速扩散到全场，抢救成本很高。我的习惯是先用 0.5 跑通，确认稳定后再逐步往上加，观察残差曲线的抖动情况，抖了就退回去。

隐式格式理论上无条件稳定，但实际受限于非线性迭代的收敛性。稳态问题常用的伪时间步方法，局部时间步按 CFL 数给：

Δt_local = CFL · V_P / Σ_f |v_n,f| · A_f

CFL 从 5 起步，收敛顺利可以逐步加到 50 甚至 100，遇到残差上翘就降回来。这个自适应策略比固定 CFL 省时间，尤其是从粗糙初场起步的时候。

5.2 一个最小可跑的实现骨架

下面这段代码把前面所有讨论落到纸面上，一维、均匀网格、周期边界、纯对流，四个格式可切换：一阶迎风、二阶迎风、QUICK、带 Van Leer 限制器的二阶迎风。可以直接跑，跑完对比扩散宽度和过冲幅度。

import numpy as np def limiter(r, kind='vanleer'): 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, scheme='upwind1'): n = len(u) F = np.zeros(n + 1) # 面 i+1/2 位于 u[i] 与 u[i+1] 之间，索引 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, n=200, cfl=0.6, t_end=0.5, a=1.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 问题速查表

6.2 几个踩过的坑和对应的解法

限制器写得对但用错了位置，是我早期最常犯的错误。有些实现会先把所有面通量算出来存进数组，再统一更新单元，这没问题；但如果限制器函数在算每个面的时候用了最新的、已经部分更新的单元值（比如原地更新 u[i] 然后继续算下一个面），就会出现方向依赖的结果，一维下表现为波形不对称，二维下表现为结果跟网格遍历顺序有关。解决办法很简单：先完整算一遍所有面通量存到数组，再统一更新，别偷懒省内存。

另一个坑是延迟修正收敛判据。前面提过，修正量没降到足够小，解其实是偏一阶的。我吃过一次亏，残差已经降到 1e-6，看着挺收敛，但出口温度跟实验数据差了三度，回头一查发现修正量还有 1e-3 量级，也就是说高阶效应根本没体现出来。后来把判据改成两个都要满足：残差小于阈值，且修正量的最大绝对值小于某个物理量的特征量乘以 1e-4。加了这条之后，同类问题再没出现过。

还有一点关于稳态和瞬态的格式选择。稳态问题用高阶迎风配延迟修正，收敛后精度是实打实的；瞬态问题如果时间步比较大，空间离散的阶数再高也没用，因为时间离散的误差会主导。这时候要么改用二阶时间格式（Crank-Nicolson 或者二阶后向差分），要么就把时间步压下来做网格无关性验证。我见过不少人在瞬态算例里死磕空间格式，调了半天限制器，其实瓶颈在时间精度上。

最后分享一个实用的小技巧：验证迎风格式实现是否正确，除了方波算例，还可以跑一个高斯脉冲的纯对流。高斯脉冲有解析解，扩散后仍然是高斯，标准差按 sqrt(σ0² + 2ν_num·t) 增长。这个算例比方波更敏感，因为它的波形处处光滑，一阶和二阶的差异一眼可见，限制器在光滑区的行为也能看出来。我一般用它做第二个验证，两个算例都过了，才认为对流项的实现可以信了。