有限体积法迎风格式实战：从数值振荡到数值扩散的取舍

管道里加一段示踪剂想看看它多久能被冲到出口这种事听起来简单得不能再简单——几何是一根长直管扩散系数小到几乎可以忽略入口给个浓度阶跃出口看响应曲线。我最初接手这个算例时用的是一套最顺手的中心差分离散残差曲线一路往下掉看着特别健康。结果把浓度云图调出来入口附近出现了负浓度浓度分布沿程还挂着一串衰减的波纹。这就是有限体积法里最经典的翻车现场对流主导时中心差分给出的不是稍微不准而是物理上根本不可能出现的振荡。后来把格式换成一阶迎风振荡立刻消失代价是阶跃被抹成一条很宽的斜坡看起来像多加了两个数量级的扩散。再往上做二阶迎风加限制器斜坡收窄波纹也没有回来。这套从崩到糊再到刚刚好的折腾过程基本就是迎风格式Upwind Scheme在有限体积法里的全部性格写照。这篇东西想讲清楚的就是迎风到底迎的是什么一阶迎风为什么要背数值扩散这笔债高阶迎风又是怎么补精度同时不把振荡放回来的以及在非结构网格和可压缩流场里这套思路怎么落地。内容偏实战适合已经能跑通一个简单有限体积求解器、正在被振荡或过度耗散折磨的人。1. 中心差分在一维纯对流上翻车的那次经历1.1 残差曲线很漂亮浓度场全是负值先把现场还原清楚。控制方程是标量输运方程$\partial \phi/\partial t u,\partial \phi/\partial x \Gamma,\partial^2\phi/\partial x^2$取 $u$ 为常数、$\Gamma$ 小到接近零。有限体积离散时面通量用相邻两单元的中心平均也就是 $F_e u(\phi_P \phi_E)/2$。这个写法对扩散项天然友好因为扩散本来就是中心性质的过程通量正比于梯度用两侧差商天经地义。麻烦在于把同样的中心平均直接搬去算对流项得到的离散方程里非对角系数会变成负的。非对角系数为负意味着什么意味着当邻居单元的浓度升高时本单元的浓度反而被推低。这在物理上无从解释——示踪剂只能由上游带进来不可能越靠近高浓度区浓度越低。数学上的后果是离散矩阵丢掉了 M 矩阵性质不再满足离散极值原理解可以在没有任何物理原因的情况下越过初边值范围。我那次算例里负浓度的量级大概在 $10^{-3}$看起来不起眼但如果这是个浓度场负浓度在后续的化学反应源项里会被放大成完全错误的反应速率。判断这件事有个特别省事的判据格子 Peclet 数 $Pe_\Delta |u|\Delta x/\Gamma$。一维稳态对流扩散用中心差分时只有当 $Pe_\Delta 2$ 时离散系数才全部为正解才不振荡。注意这个判据跟时间步、跟格式精度都没关系只看网格尺度和当地的对流-扩散比。实际工程里稍微像个对流问题的算例$Pe_\Delta$ 动辄几百上千中心差分必炸只是炸法不同定常问题表现为解发散或出现非物理解非定常问题表现为波纹沿下游传播。离散方式非对角系数符号振荡风险主要代价中心差分可正可负$Pe_\Delta 2$ 时必然振荡无非物理耗散但可能完全不可用一阶迎风恒非正基本无强数值扩散精度一阶二阶迎风可正光滑区可用间断附近振荡需要限制器或足够细的网格混合格式视权重而定一般可控权重的选取带经验性1.2 从特征线看信息只从上游来要理解为什么中心差分错了换个角度想稳态对流方程 $u,\partial\phi/\partial x 0$ 的通解是 $\phi \text{const}$但更本质的结构是双曲型方程的解沿特征线 $dx/dt u$ 传播。$u 0$ 时$x$ 处解的值只由 $x$ 上游的历史决定下游的信息对当前点没有任何影响。这是双曲型方程的方向性也叫依赖域的单侧性。中心差分把左右两侧的邻居等权重地掺和进来等于默认下游也在影响我。在 $Pe_\Delta$ 很大时这种双向依赖与方程真实的单侧依赖严重冲突冲突的表现就是数值解在冲击附近产生高频振荡。反过来说如果格式的依赖域与方程的依赖域一致——只用上游信息——就不会产生这种非物理振荡。这就是迎风格式的出发点也是它名字的由来面向来流方向取值。这个逻辑有个容易被忽略的推论迎风是有方向性的不是一种通用更好的离散。对于扩散主导的问题$Pe_\Delta$ 很小双向依赖才是真实物理这时候上迎风反而引入不必要的偏差。所以真正稳妥的做法不是无脑全用迎风而是按当地 $Pe_\Delta$ 或按特征方向决定用哪种。混合格式、混合迎风-中心格式动机都在这里。还有个细节值得说判断谁是上游在多维、非结构网格上并不总那么显然。一维时看 $u$ 的符号就够二维时可能要同时考虑面法向速度分量可压缩流里还要看各个特征波的传播方向。这部分放到第 4 章和第 5 章细讲。2. 一阶迎风的构造与它背着的数值扩散债2.1 面通量的取值规则两个候选值挑一个一阶迎风的构造简单到有点粗暴面上需要 $\phi$ 的值我不做插值直接取上游单元的中心值。若 u_e 0: F_e u_e * phi_P 若 u_e 0: F_e u_e * phi_E 统一写法: F_e u_e * [ max(u_e,0)*phi_P min(u_e,0)*phi_E ] / u_e写成一个更常用的紧凑形式引入正负分解 $u^ \max(u,0)$、$u^- \min(u,0)$面通量就是 $F_e u^ \phi_P u^- \phi_E$。这个写法可以直接推广到多维和非结构网格把 $u$ 换成面法向速度 $v_f \mathbf{u}_f\cdot\mathbf{n}_f$把 $P$ 和 $E$ 换成上风侧单元和下风侧单元形式完全不变。这个格式为什么稳定把一维均匀网格的显式更新写出来$u0$ 时 $\phi_i^{n1} \phi_i^n - C(\phi_i^n - \phi_{i-1}^n)$其中 $C u\Delta t/\Delta x$ 是 Courant 数。当 $0 \le C \le 1$ 时新值是两个非负旧值的凸组合系数全部非负且和为 1。凸组合的意义在于新值永远落在参与组合的那些值的最大最小值之间不会凭空越界。这条性质在数值上等价于单调性也是它绝不振荡的根本原因。用有限体积的正式语言说一阶迎风离散得到的系数矩阵是 M 矩阵对角元为正非对角元非正且行和对角占优。M 矩阵保证解满足离散极值原理这对浓度、组分质量分数、相分数这类有物理上下界的量至关重要。我后来做多相流和组分输运时只要涉及这类必须保持正性的标量一阶迎风甚至是先保证跑得动再谈精度的兜底选项。2.2 修正方程算账数值扩散到底有多大稳定性有了精度问题立刻浮出来。用泰勒展开把离散格式还原成等效的偏微分方程修正方程一阶迎风给出的不是原方程而是∂φ/∂t u ∂φ/∂x (u Δx / 2)(1 - C) ∂²φ/∂x² O(Δx²)多出来的那项是二阶导数形式上跟物理扩散一模一样所以叫数值扩散numerical diffusion扩散系数 $\Gamma_{num} \dfrac{u\Delta x}{2}(1 - C)$。这个系数有几点很值得琢磨。第一它正比于 $\Delta x$所以网格加密一倍数值扩散减半——收敛还是收敛的一阶收敛。第二它正比于 $u$对流越强、污染越重这跟物理扩散的性质完全不同物理扩散与速度无关。第三它带着 $(1-C)$ 因子$C1$ 时数值扩散为零。这不是巧合$C1$ 时显式迎风恰好把上一个单元的值原封不动搬到当前单元网格间距正好等于一个时间步内流过的距离所以是精确解。很多人拿这个当调试技巧用 $C1$ 跑一遍验证输运实现有没有写错我也常用。给个量化的直观感受。假设 $u 1\ \text{m/s}$$\Delta x 1\ \text{mm}$$\Gamma_{num} \approx 5\times10^{-4}\ \text{m}^2/\text{s}$。如果这是水中的分子扩散量级在 $10^{-9}$数值扩散直接把扩散效应放大了五个数量级。这就是为什么一阶迎风算出来的阶跃看起来糊——它实际上在以一个人为设定的巨大扩散系数做平滑。稳态情形下 $\partial\phi/\partial t 0$数值扩散没有 $(1-C)$ 因子直接是 $u\Delta x/2$。这个差别容易搞混非定常的数值扩散随 Courant 数变化稳态的对角占优判断则与时间步无关。做稳态计算时即便用隐式推进把 $C$ 取得很大数值扩散也不会因此减小。2.3 正性保持与离散极值原理一阶迎风的隐藏优点这里说个我早年没想明白、后来才体会到的点数值扩散虽然讨厌但它带来的正性保持并不是白给的。很多高阶格式之所以要在间断附近退回一阶就是因为高阶插值可能造出超出邻居范围的值而一阶迎风天然不会。具体到操作层面有个判断方式值得记住看单元的更新式里旧时刻各单元值的系数是否全部非负。只要全部非负就不会出现新的极值。任何重构reconstruction只要在某处产生了负系数那一点就有振荡风险。限制器的作用本质就是在重构里做事后修正把可能产生负系数的插值斜率削回去。工程上还有一个实际问题源项。如果源项是消耗性的比如反应速率、衰变量级大、时间步又不小一阶迎风加上显式源项都可能把标量推成负值。这时候光靠格式不够还得对源项做线性化处理把消耗项放到隐式侧用当前值做系数。我见过不止一次把负浓度归罪于格式、结果排查半天发现是源项线性化没做。3. 想提精度就得付出代价高阶迎风的三条路线3.1 二阶迎风与QUICK多取一个上游点一阶迎风只用了一个上游点插值精度只能是一阶。想要二阶最直接的想法是多取一个点。对 $u0$ 的面 $e$如果同时用 $\phi_W$、$\phi_P$、$\phi_E$ 三点做二次插值就得到二阶迎风SOUφ_e 1.5 φ_P - 0.5 φ_W 二阶迎风u 0QUICK 格式是同一思路的加权变体权重按到面的距离取均匀网格下退化成φ_e (6/8) φ_P (3/8) φ_E - (1/8) φ_W QUICKu 0QUICK 的权重设计意图是让三点插值在均匀网格上达到三阶截断误差对面上值的插值而言实际整体格式仍是二阶。相比二阶迎风它对下风侧保留了 $3/8$ 的权重这一点在光滑区能改善相位误差但在陡峭梯度附近恰恰是最容易出问题的地方。这两条路线共同的麻烦是系数不再保证非负。二阶迎风里 $-\phi_W$ 那一项系数是负的QUICK 里 $-\phi_W$ 也是负的。只要梯度足够陡插值出来的面值就可能超出 $[\phi_W, \phi_E]$ 的范围振荡就回来了。所以直接用纯二阶迎风的求解器在激波、阶跃、自由面这类地方几乎一定需要额外的抑制手段否则要么发散要么出现波纹。高精度带来的第二个代价是计算模板变宽。并行分区时二阶格式需要多一层 halo 单元非结构网格上三点模板在畸变网格上可能根本找不到合适的三点排列。这些都是实现层面躲不开的活。格式均匀网格截断误差面值公式u0是否保正典型用途一阶迎风一阶$\phi_P$是兜底、初场、强源项问题二阶迎风二阶$1.5\phi_P-0.5\phi_W$否光滑区、结构化网格QUICK二阶面值三阶$0.75\phi_P0.375\phi_E-0.125\phi_W$否结构化网格、需低相位误差TVD/MUSCL二阶光滑区$\phi_P\frac{1}{2}\psi(r)(\phi_P-\phi_W)$是通用间断问题首选3.2 延迟修正隐式矩阵留在低阶高阶项进源项这里有个实现上的现实问题一阶迎风的系数矩阵是稀疏、对角占优、好解换成二阶迎风模板变宽非结构网格上更麻烦而且负系数会破坏对角占优迭代求解器的收敛性变差。早年很多代码用了一个很聪明的绕法叫延迟修正deferred correction。思路是把高阶通量拆开F_high F_low (F_high - F_low) ↑隐式 ↑显式用上一迭代步的值前一项用一阶迎风的形式参与系数矩阵组装保持矩阵性质良好后一项是高阶与低阶之差当作源项处理用上一轮迭代的已知值算出来。收敛时括号里那一项把低阶的数值扩散正好抵消掉最终的离散精度是高阶的。这个做法我觉得是整个迎风格式体系里最实用的工程技巧之一。它的好处很直接矩阵结构不变可以直接复用一阶迎风的求解器实现上只需要额外算一次高阶面值和一次低阶面值做差内存开销几乎为零。坏处也同样直接收敛过程被削弱了因为显式那部分降低了迭代矩阵的对角占优程度大 Courant 数下可能收敛变慢甚至振荡。用延迟修正的算例我一般会把亚松弛因子调低一些或者先跑几百步一阶把场养起来再切高阶。有个坑要说清楚延迟修正的收敛判据不能只看残差。因为高阶项是挂在源项上的残差下降只说明线性系统解好了不代表高阶项和低阶项已经互相抵消到位。稳妥的做法是同时监控几个物理量的积分值或监测点值看着它们不再漂移才算真收敛。我就吃过这个亏残差降到 $10^{-6}$出口浓度还在缓慢爬升。3.3 TVD与MUSCL让限制器在激波处刹车要既保住二阶精度、又不在间断处振荡主流做法是通量限制器flux limiter。用 MUSCL 的写法把面值写成中心值加上一个受限制的修正φ_e φ_P (1/2) ψ(r) (φ_P - φ_W) r (φ_P - φ_W) / (φ_E - φ_P)$\psi(r)$ 就是限制器函数。$r$ 反映梯度的光滑程度$r0$ 说明相邻梯度同号光滑过渡$r\le 0$ 说明出现极值可能是不连续此时把 $\psi$ 压到零格式局部退化成一阶迎风。常见的几个限制器限制器表达式特点Minmod$\max(0,\min(1,r))$最保守耗散最大绝不振荡Van Leer$(rrSuperbee$\max(0,\min(2r,1),\min(r,2))$最激进间断最锐可能削平光滑极值Van Albada$(r^2r)/(r^21)$连续光滑迭代稳定性好TVD 条件给出的可行域是 $0\le\psi(r)\le 2$ 且 $0\le\psi(r)\le 2r$Sweby 图里的那两块区域。选哪个限制器本质上是在间断锐度和光滑区精度之间选一个位置越靠近可行域上边界间断越锐但越容易把光滑极值削成平台越靠下边界越稳但耗散越接近一阶。我实际选型的经验是通用问题先上 van Leer它对光滑极值相对友好不容易把本该有的物理峰值抹平如果算的是强激波、需要极锐的间断再考虑 superbee如果不确定性很大或者网格质量差minmod 最保险牺牲点精度换不折腾。非结构网格上很多代码会用 Venkatakrishnan 限制器它的特点是连续可微不会因为限制器本身的切换造成迭代抖动。还有个必须提醒的地方限制器在光滑极值附近会把精度降回一阶。这是 TVD 类格式的固有性质不是实现 bug。如果你算的是个平滑的波峰值位置恰好是极值点限制器在那里把斜率削掉长期积分会看到峰值缓慢衰减。这叫削峰clipping。解决办法是上更宽松的限制器或者改用 WENO 这类不依赖 TVD 条件的重构。我在算声波传播这类问题时被这个现象坑过波形越传越矮一度以为是数值耗散其实是限制器在削峰。4. 把迎风搬到非结构网格与多维流场4.1 面法向速度决定谁在上游一维时判断上游只看 $u$ 的符号二维、三维非结构网格上这一步变成了算面法向速度、判断正负、选上风单元。具体流程是先由动量方程得到面上的速度 $\mathbf{u}_f$可能是插值得到的也可能是 Rhie-Chow 类方法构造的再算面法向分量 $v_f \mathbf{u}_f\cdot\mathbf{n}_f$其中 $\mathbf{n}_f$ 是面的外法向统一按 owner 单元指向 neighbor 单元的方向定义然后若 v_f 0: 上风单元 owner下风单元 neighbor 若 v_f 0: 上风单元 neighbor下风单元 owner 一阶: φ_f φ_upwind这里面最容易被低估的是面速度的获取。不可压求解器里速度存储在单元中心面上必须插值。如果只是简单线性插值得到的 $v_f$ 和压力梯度之间会产生棋盘失耦这就是著名的压力-速度耦合问题需要 Rhie-Chow 插值来引入压力梯度的紧致性。这件事跟迎风不是一回事但两者经常一起出现混淆的人不少。简单说Rhie-Chow 解决的是面上的法向速度怎么算迎风解决的是面上的标量怎么取别混为一谈。另一个实操细节面法向方向的定义必须在整个网格里统一。如果一部分面按 owner→neighbor 定义另一部分反过来迎风选边就会随面而异得到的结果看起来像是随机的耗散分布。这种 bug 特别隐蔽因为残差照样下降只是结果不对。我见过一个案例排查了两天才发现是网格处理脚本里某个面的法向符号反了。4.2 梯度重构、限制器与网格质量的纠缠二阶迎风在非结构网格上的标准写法是重构reconstruction用上风单元的梯度把中心值外推到面上。φ_f φ_up ψ_f (∇φ)_up · (x_f - x_up)其中 $(\nabla\phi)_{up}$ 是上风单元的梯度用最小二乘或格林-高斯方法算出来$x_f - x_up$ 是从上风单元中心到面中心的向量$\psi_f$ 是作用于该面的限制器注意这跟 3.3 节里基于 $r$ 的限制器不完全一样非结构网格上常用的是基于邻居极值的限制比如 Venkatakrishnan 或 Barth-Jespersen。这里的坑集中在两处。第一处是梯度重构的质量直接受网格质量影响。高度畸变的四面体网格上最小二乘梯度可能非常不准外推到面上时误差被放大即使加了限制器也会看到局部精度掉到一阶甚至更低。第二处是限制器的作用范围。限制器成功的前提是知道邻居的极值范围因此并行分区时必须保证有足够厚的 halo 层否则限制器在分区边界上失效接缝处会出现可见的精度突变。这个现象在多块并行计算里挺常见云图上表现为规整的方块状纹理。还有个更隐晦的问题非结构网格上上风单元唯一确定但如果面法向速度接近零滞止点、回流区边界上风方向会在两个时间步之间来回跳。跳变时面值从 $\phi_{owner}$ 突变成 $\phi_{neighbor}$相当于在解里注入噪声。处理办法通常是给 $v_f$ 加一个判断阈值接近零时退化成两侧平均或者直接保证该处的 $|v_f|$ 不会因为插值精度而改变符号。4.3 对流扩散混合时的格子Peclet数与混合格式实际问题很少是纯对流。对流扩散并存时迎风的必要性要看当地 $Pe_\Delta$。有个做法叫混合格式hybrid scheme在中心差分和迎风之间按 $Pe_\Delta$ 加权$Pe_\Delta 2$ 用中心差分否则用迎风。这本质上是通过引入数值扩散把 $Pe_\Delta$ 压到 2 以下思路和用扩散换稳定完全一致。这个做法在今天看来有点老派但它有个很实用的性质实现简单、无条件稳定、不需要限制器。如果算的是低雷诺数流动、传热问题、或者网格已经足够细$Pe_\Delta$ 本来就接近 2混合格式的效果和二阶迎风差别不大而后者可能还要担心限制器和收敛性。我处理一些热传导为主、对流只是伴随效应的算例时会直接用混合格式省事。但对流主导时不建议靠混合格式硬扛。它引入的数值扩散是实实在在的等于解了一个物理上不存在的高扩散问题得到的边界层、剪切层、温度梯度都会被显著抹厚。这类算例值得老老实实上二阶迎风加限制器再配合网格加密。判断什么时候该加密网格有个粗略的经验算完之后取几个关键位置的解把网格加密一倍再算一遍看结果变化多少。如果变化明显说明数值扩散还在主导继续加密如果变化很小说明已经进入渐近区。这个方法比看残差靠谱得多也适用于判断限制器有没有把解削过头。5. 通量分裂迎风思想的另一条分支5.1 Godunov型格式与黎曼问题的关系前面讲的都是标量输运。放到可压缩流动里迎风这件事就变成了通量分裂flux splitting。核心思想把通量按特征方向拆开各走各的迎风。Godunov 的思路最干净假设每个单元内的解是分片常数那么在两个单元之间就形成一个局部的黎曼问题——左边状态 $U_L$、右边状态 $U_R$中间会生成波系。面上的数值通量取黎曼解在 $x/t0$ 处的值。对线性标量方程一阶 Godunov 格式退化成正好的一阶迎风这印证了两者的等价性。问题在于欧拉方程的黎曼解需要迭代代价太高。于是有了各种近似黎曼解。Roe 格式是最经典的一个通量写成F 0.5 (F_L F_R) - 0.5 |A| (U_R - U_L)$|A|$ 是 Roe 平均雅可比矩阵的绝对值形式上就是把每个特征波按传播方向拆开、分别取迎风侧。看到这个结构你应该觉得眼熟第一项是中心部分第二项是迎风修正修正的方向由特征值符号决定。标量情形下这就是 $u^\phi_L u^-\phi_R$。Roe 格式的著名弱点是膨胀波处不满足熵条件会产生非物理的膨胀激波实际代码里都要加一个熵修正entropy fix。类似思路还有 HLL/HLLC 系列它们不分解完整波系只假设两三道波代价低、鲁棒性好在工程代码里用得极广。AUSM 系列则把通量拆成对流部分和压力部分分别做迎风处理形式简单、不需要特征分解对低速和跨声速都还比较友好。这几种格式的共同点是都保留了按特征方向取迎风侧的核心思想区别在于怎么近似特征结构、怎么处理退化情形。选型的经验是强激波、高超问题优先考虑 HLLC 或带熵修正的 Roe一般工程流动 AUSM 类足够对鲁棒性要求极高、能接受多一些耗散就用 HLL。通量格式迎风思想体现优点短板Godunov精确黎曼直接取黎曼解最符合物理计算量大Roe按 $A$ 分解特征波HLLC三道波近似鲁棒、保正性好接触面精度略逊 RoeAUSM 类对流项与压力项分别迎风形式简洁、低速友好低马赫数下需额外处理5.2 压力-速度耦合的迎风处理与低速退化可压缩求解器里有个绕不开的麻烦马赫数很低时对流项的迎风耗散正比于 $|u|c$$c$ 是声速而物理对流只正比于 $|u|$。比值是 $11/M$ 量级$M0.01$ 时迎风耗散是被需要的上百倍。表现就是低速区被糊得非常厉害压力场出现异常波动收敛变慢。这就是为什么很多可压缩代码算不了低速流动。应对办法有几类时间预处理preconditioning改特征值让它们在低马赫数下趋于一致通量分裂类的格式里显式对耗散项做马赫数缩放或者干脆换到不可压/低速求解器框架。这些内容超出了标量迎风格式的范围但有个认知值得带上迎风耗散不是固定的它会随当地特征速度变化在极端条件下可能远超物理这一点在任何格式选型里都要先想一想。6. 实测踩坑记录与验证建议6.1 收敛了不等于算对了三个自检动作有限体积计算里最容易犯的错误是把残差降下去了当成算对了。迎风格式的数值扩散恰好会让收敛性变好——耗散越大迭代矩阵越接近对角占优收敛越快。所以你用一阶迎风得到的漂亮收敛曲线很可能掩盖着一个被抹平的解。我现在固定做三件事。第一件检查解的极值有没有越界比如质量分数是否在 $[0,1]$ 内、浓度是否非负。一旦越界先不用管是格式还是源项的问题先把范围压回来再说。第二件做网格无关性检查加密一档看关键量的变化。第三件对比至少两个格式的结果一阶和二阶各跑一遍差异大的地方就是需要重点确认的区域。这三件事加起来增加不了多少机时但能挡掉大部分结果看起来对但其实是数值产物的情况。还有个细节迎风格式在初场的建立阶段特别好用。我一般先跑几百步一阶迎风让场稳定下来再切到高阶。直接在零初场上高阶开局限制器会因为初始的强梯度疯狂削峰很容易把中间产物削成非物理解然后整场跟着崩。这个先养后切的习惯用了很多年几乎每次都能省掉一堆收敛性上的麻烦。6.2 边界处的迎风入口出口不对称处理边界上的迎风处理经常被写成入口给 Dirichlet出口给零梯度然后就出问题了。问题在于迎风的取法在入口边界上并不对称入口边界外没有单元面值只能由给定的边界值确定这是正确的因为信息从外部流入出口边界上信息是从内部流出的如果用边界处的内部值做迎风格式才一致。具体一点在出口面上正确的做法是把内部上风单元的值拿去做迎风而不是拿一个外推值。如果误用零梯度边界条件覆盖掉迎风取值等价于在出口加了一个反射或阻尼会出现明显的出口反射波在云图上表现为出口附近的一圈波纹。我碰到过一次算浓度射流的问题出口浓度出现了非物理的振荡翻来翻去发现是出口面用了固定值边界而不是迎风外推。还有一种常见问题入口给了强梯度比如阶跃入口限制器在入口附近会把斜率削得很厉害导致入口附近的解精度退化。解决办法是入口附近加密网格或者用较宽松的限制器。这件事不太可能靠改格式根治因为入口处的间断是物理存在的限制器必然要反应。6.3 用一维和二维标准算例量误差阶要确认自己的迎风实现对不对跑标准算例比看云图有效得多。一维纯对流是最基本的给一个平滑初值正弦或高斯常速度输运一个周期看波形是否回到原处、峰值衰减多少。一阶迎风的峰值衰减会非常明显二阶加限制器会好很多但仍有削峰。这个算例的另一个用途是量误差阶在 $L_1$ 或 $L_2$ 范数下一阶迎风应该能看到一阶收敛率二阶格式在光滑解上应该是二阶。收敛率测出来不对八成是下面几类问题时间推进和空间离散阶数不匹配时间用一阶、空间用二阶整体收敛率会被时间项拉成一阶限制器在光滑解上也频繁触发比如 $r$ 的计算写错了符号或者邻居访问越界取了零边界处理降阶入口阶跃附近本来就是一阶。我建议至少跑一遍三套网格用相邻两套的误差比取对数算收敛率比单看一套的误差可靠。二维可以上纯对流的对角输运或旋转流场比如把一个高斯团绕着中心转一圈这类算例能暴露非结构网格上的重构和限制器问题尤其看限制器有没有在旋转过程中持续削峰。再往上就是经典的对流扩散测试有解析解的那类用来看 $Pe_\Delta$ 不同取值下格式的表现。这几个算例都不复杂但把它们都跑通了对格式的行为会有一个非常具体的认识比读十篇格式综述管用。我个人在实操中的体会是迎风格式这件事的难点从来不在公式公式都是三五行的东西难在知道什么时候该用它、什么时候该从它退出来。一个算例里往往同时存在对流主导区和扩散主导区一套格式打天下是不可能的。多阶段推进、分区用不同格式、延迟修正、限制器的松紧调节这些手段的组合才是实际干活的样子。至于数值扩散这位老朋友我的态度是先把它当成保证能跑下去的代价接受下来再通过网格加密和格式升级一点点把它换掉而不是指望一上来就找到既无耗散又绝对稳定的完美格式——那种东西不存在。