管道里加一段示踪剂,想看看它多久能被冲到出口,这种事听起来简单得不能再简单——几何是一根长直管,扩散系数小到几乎可以忽略,入口给个浓度阶跃,出口看响应曲线。我最初接手这个算例时,用的是一套最顺手的中心差分离散,残差曲线一路往下掉,看着特别健康。结果把浓度云图调出来,入口附近出现了负浓度,浓度分布沿程还挂着一串衰减的波纹。这就是有限体积法里最经典的翻车现场:对流主导时,中心差分给出的不是"稍微不准",而是物理上根本不可能出现的振荡。
后来把格式换成一阶迎风,振荡立刻消失,代价是阶跃被抹成一条很宽的斜坡,看起来像多加了两个数量级的扩散。再往上做二阶迎风加限制器,斜坡收窄,波纹也没有回来。这套从"崩"到"糊"再到"刚刚好"的折腾过程,基本就是迎风格式(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$ 换成"上风侧单元"和"下风侧单元",形式完全不变。
这个格式为什么稳定?把一维均匀网格的显式更新写出来,$u>0$ 时 $\phi_i^{n+1} = \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)$ 因子,$C=1$ 时数值扩散为零。这不是巧合:$C=1$ 时显式迎风恰好把上一个单元的值原封不动搬到当前单元,网格间距正好等于一个时间步内流过的距离,所以是精确解。很多人拿这个当调试技巧,用 $C=1$ 跑一遍验证输运实现有没有写错,我也常用。
给个量化的直观感受。假设 $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:多取一个上游点
一阶迎风只用了一个上游点,插值精度只能是一阶。想要二阶,最直接的想法是多取一个点。对 $u>0$ 的面 $e$,如果同时用 $\phi_W$、$\phi_P$、$\phi_E$ 三点做二次插值,就得到二阶迎风(SOU):
φ_e = 1.5 φ_P - 0.5 φ_W (二阶迎风,u > 0)QUICK 格式是同一思路的加权变体,权重按到面的距离取,均匀网格下退化成:
φ_e = (6/8) φ_P + (3/8) φ_E - (1/8) φ_W (QUICK,u > 0)QUICK 的权重设计意图是让三点插值在均匀网格上达到三阶截断误差(对面上值的插值而言),实际整体格式仍是二阶。相比二阶迎风,它对下风侧保留了 $3/8$ 的权重,这一点在光滑区能改善相位误差,但在陡峭梯度附近恰恰是最容易出问题的地方。
这两条路线共同的麻烦是:系数不再保证非负。二阶迎风里 $-\phi_W$ 那一项系数是负的,QUICK 里 $-\phi_W$ 也是负的。只要梯度足够陡,插值出来的面值就可能超出 $[\phi_W, \phi_E]$ 的范围,振荡就回来了。所以直接用纯二阶迎风的求解器,在激波、阶跃、自由面这类地方几乎一定需要额外的抑制手段,否则要么发散,要么出现波纹。
高精度带来的第二个代价是计算模板变宽。并行分区时,二阶格式需要多一层 halo 单元;非结构网格上,三点模板在畸变网格上可能根本找不到合适的三点排列。这些都是实现层面躲不开的活。
| 格式 | 均匀网格截断误差 | 面值公式(u>0) | 是否保正 | 典型用途 |
|---|---|---|---|---|
| 一阶迎风 | 一阶 | $\phi_P$ | 是 | 兜底、初场、强源项问题 |
| 二阶迎风 | 二阶 | $1.5\phi_P-0.5\phi_W$ | 否 | 光滑区、结构化网格 |
| QUICK | 二阶(面值三阶) | $0.75\phi_P+0.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$ 反映梯度的光滑程度:$r>0$ 说明相邻梯度同号(光滑过渡),$r\le 0$ 说明出现极值(可能是不连续),此时把 $\psi$ 压到零,格式局部退化成一阶迎风。
常见的几个限制器:
| 限制器 | 表达式 | 特点 |
|---|---|---|
| Minmod | $\max(0,\min(1,r))$ | 最保守,耗散最大,绝不振荡 |
| Van Leer | $(r+ | r |
| Superbee | $\max(0,\min(2r,1),\min(r,2))$ | 最激进,间断最锐,可能削平光滑极值 |
| Van Albada | $(r^2+r)/(r^2+1)$ | 连续光滑,迭代稳定性好 |
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/t=0$ 处的值。对线性标量方程,一阶 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 | 三道波近似 | 鲁棒、保正性好 | 接触面精度略逊 Roe |
| AUSM 类 | 对流项与压力项分别迎风 | 形式简洁、低速友好 | 低马赫数下需额外处理 |
5.2 压力-速度耦合的迎风处理与低速退化
可压缩求解器里有个绕不开的麻烦:马赫数很低时,对流项的迎风耗散正比于 $|u|+c$($c$ 是声速),而物理对流只正比于 $|u|$。比值是 $1+1/M$ 量级,$M=0.01$ 时迎风耗散是被需要的上百倍。表现就是低速区被"糊"得非常厉害,压力场出现异常波动,收敛变慢。这就是为什么很多可压缩代码算不了低速流动。
应对办法有几类:时间预处理(preconditioning),改特征值让它们在低马赫数下趋于一致;通量分裂类的格式里显式对耗散项做马赫数缩放;或者干脆换到不可压/低速求解器框架。这些内容超出了标量迎风格式的范围,但有个认知值得带上:迎风耗散不是固定的,它会随当地特征速度变化,在极端条件下可能远超物理,这一点在任何格式选型里都要先想一想。
6. 实测踩坑记录与验证建议
6.1 收敛了不等于算对了:三个自检动作
有限体积计算里最容易犯的错误,是把"残差降下去了"当成"算对了"。迎风格式的数值扩散恰好会让收敛性变好——耗散越大,迭代矩阵越接近对角占优,收敛越快。所以你用一阶迎风得到的漂亮收敛曲线,很可能掩盖着一个被抹平的解。
我现在固定做三件事。第一件,检查解的极值有没有越界,比如质量分数是否在 $[0,1]$ 内、浓度是否非负。一旦越界,先不用管是格式还是源项的问题,先把范围压回来再说。第二件,做网格无关性检查,加密一档看关键量的变化。第三件,对比至少两个格式的结果,一阶和二阶各跑一遍,差异大的地方就是需要重点确认的区域。这三件事加起来增加不了多少机时,但能挡掉大部分"结果看起来对但其实是数值产物"的情况。
还有个细节:迎风格式在初场的建立阶段特别好用。我一般先跑几百步一阶迎风让场稳定下来,再切到高阶。直接在零初场上高阶开局,限制器会因为初始的强梯度疯狂削峰,很容易把中间产物削成非物理解,然后整场跟着崩。这个"先养后切"的习惯用了很多年,几乎每次都能省掉一堆收敛性上的麻烦。
6.2 边界处的迎风:入口出口不对称处理
边界上的迎风处理经常被写成"入口给 Dirichlet,出口给零梯度",然后就出问题了。问题在于迎风的取法在入口边界上并不对称:入口边界外没有单元,面值只能由给定的边界值确定(这是正确的,因为信息从外部流入);出口边界上信息是从内部流出的,如果用边界处的内部值做迎风,格式才一致。
具体一点:在出口面上,正确的做法是把内部上风单元的值拿去做迎风,而不是拿一个外推值。如果误用零梯度边界条件覆盖掉迎风取值,等价于在出口加了一个反射或阻尼,会出现明显的出口反射波,在云图上表现为出口附近的一圈波纹。我碰到过一次算浓度射流的问题,出口浓度出现了非物理的振荡,翻来翻去发现是出口面用了固定值边界而不是迎风外推。
还有一种常见问题:入口给了强梯度(比如阶跃入口),限制器在入口附近会把斜率削得很厉害,导致入口附近的解精度退化。解决办法是入口附近加密网格,或者用较宽松的限制器。这件事不太可能靠改格式根治,因为入口处的间断是物理存在的,限制器必然要反应。
6.3 用一维和二维标准算例量误差阶
要确认自己的迎风实现对不对,跑标准算例比看云图有效得多。一维纯对流是最基本的:给一个平滑初值(正弦或高斯),常速度输运一个周期,看波形是否回到原处、峰值衰减多少。一阶迎风的峰值衰减会非常明显,二阶加限制器会好很多但仍有削峰。这个算例的另一个用途是量误差阶:在 $L_1$ 或 $L_2$ 范数下,一阶迎风应该能看到一阶收敛率,二阶格式在光滑解上应该是二阶。
收敛率测出来不对,八成是下面几类问题:时间推进和空间离散阶数不匹配(时间用一阶、空间用二阶,整体收敛率会被时间项拉成一阶);限制器在光滑解上也频繁触发(比如 $r$ 的计算写错了符号,或者邻居访问越界取了零);边界处理降阶(入口阶跃附近本来就是一阶)。我建议至少跑一遍三套网格,用相邻两套的误差比取对数算收敛率,比单看一套的误差可靠。
二维可以上纯对流的对角输运或旋转流场(比如把一个高斯团绕着中心转一圈),这类算例能暴露非结构网格上的重构和限制器问题,尤其看限制器有没有在旋转过程中持续削峰。再往上就是经典的对流扩散测试(有解析解的那类),用来看 $Pe_\Delta$ 不同取值下格式的表现。这几个算例都不复杂,但把它们都跑通了,对格式的行为会有一个非常具体的认识,比读十篇格式综述管用。
我个人在实操中的体会是,迎风格式这件事的难点从来不在公式,公式都是三五行的东西,难在知道什么时候该用它、什么时候该从它退出来。一个算例里往往同时存在对流主导区和扩散主导区,一套格式打天下是不可能的。多阶段推进、分区用不同格式、延迟修正、限制器的松紧调节,这些手段的组合才是实际干活的样子。至于数值扩散这位老朋友,我的态度是:先把它当成保证能跑下去的代价接受下来,再通过网格加密和格式升级一点点把它换掉,而不是指望一上来就找到既无耗散又绝对稳定的完美格式——那种东西不存在。