第 12 章：高分辨率格式（High Resolution Schemes）
12.1 归一化变量格式（The Normalized Variable Formulation, NVF）
归一化变量格式（NVF）由 Leonard 引入并因 Gaskell 与 Lau 简化的对流有界性准则（CBC）而广为流传；归一化变量图（NVD）则是分析高阶（HO）与高分辨率（HR）格式的有力工具。NVF 是一种面（face）格式构造方法：在面 \(f\) 上要构造 \(\phi_f\) 的值时，先对因变量在局部做归一化，并用上游节点 C、下游节点 D 以及远上游节点 U 的值来表达归一化变量 \(\tilde\phi_f\)。对均匀网格作图时，U、C、D 沿同一条线排列；对于非结构网格，则把 C、D 与外推出来的 U 放在一起参考（图 12.1）。归一化定义本身是
把 \(\phi_f = f(\phi_U, \phi_C, \phi_D)\) 改写成 \(\tilde\phi_f = \tilde f(\tilde\phi_C)\)，是因为归一化后 \(\tilde\phi_U = 0\)、\(\tilde\phi_D = 1\)；此时 \(\tilde\phi_C\) 成为场 \(\phi\) 平滑性的指示器：\(\tilde\phi_C \in (0, 1)\) 对应单调剖面，\(\tilde\phi_C < 0\) 或 \(\tilde\phi_C > 1\) 意味着 C 点处出现极值，而 \(\tilde\phi_C = 0\) 或 \(\tilde\phi_C = 1\) 则指示梯度跳跃（图 12.2）。
归一化的另一项作用是把各种高阶格式转化为 \(\tilde\phi_f\) 与 \(\tilde\phi_C\) 之间的线性关系。例如第 11 章讨论的各 HO 格式在归一化后变成
对所有基于三个节点值的高阶格式，\(\tilde\phi_f\) 都可以写成 \(\tilde\phi_f = \ell \tilde\phi_C + k\)，其中 \(\ell\) 与 \(k\) 是与格式有关的常数；于是把这些直线一起画在 \((\tilde\phi_C, \tilde\phi_f)\) 平面上便得到 NVD（图 12.3）。例 1 与例 2 是这一节的两个算例：例 1 演示如何从
（QUICK）出发，对两端做归一化并利用
这一关键等式，得到
；例 2 证明对在均匀笛卡尔网格上发展的二阶格式，Taylor 展开给出
，代入一般形式后只在 \(\tilde\phi_C = 1/2\)、\(\tilde\phi_f = 3/4\) 这一个公共点上恒等地成立，因此所有二阶格式必过 NVD 上的点 \(Q(1/2, 3/4)\)。
NVD 的几何意义还有更多：Upwind 格式非常扩散（diffusive），Downwind 格式非常压缩（compressive，即反扩散），任何接近 Upwind 直线的格式就偏扩散，接近 Downwind 直线的就偏压缩。本书第 11 章展示过 HO 格式虽然把 Upwind 的截断误差大幅降下来、保持稳定，但本身是无界的——它们在突变或陡梯度附近会产生欠/过冲甚至振荡；在某些场景里小幅过冲可以容忍，但当被输运变量本身是湍流黏性系数等关键量时，振荡会带来灾难性后果。这种振荡行为对所有 HO 线性对流格式都存在：它们不是极值保持的（extrema preserving），最大不会单调不增、最小不会单调不减。Godunov 与 Ryabenki 证明线性单调格式至多一阶，所以任何二阶及以上的线性格式都做不到极值保持；要想得到极值保持的格式必须使用非线性限制子函数（limiter function）。研究中形成了两大类方法：第一类向一阶迎风加一个受限的反扩散通量（flux limiter），第二类在无界 HO 格式上叠一层平滑扩散（flux blending）。前者多步且要在两个通量之间精细平衡，计算代价较高，所以本书只介绍属于第一类的两种 HR 构造：把 HO 格式与有界低阶格式复合，以一个判据切换（composite scheme），或对一阶迎风加一个乘以限制子的反扩散通量（TVD 格式）。复合方案先在 NVF/NVD 框架下给出，CBC 准则正是在此框架下为高阶插值的有界性而引入。
12.2 对流有界性准则（The Convection Boundedness Criterion, CBC）
数值格式应当保留所描述物理现象的固有属性。考虑对流的物理本质——把流体属性从上游输运到下游——可以理解一个对流近似必须具有的"输运性"：数值对流格式应该迎风偏置，否则会失去对流稳定性。因此分析对流格式时除了包围界面的 C、D 两个节点外，还要看远上游节点 U 的值；NVF 把 U 的影响也纳入分析，对识别格式是否单调极有帮助。Spekreijse、Barth 与 Jespersen 的"单调/有界"定义涉及界面的所有邻居，而 Leonard、Gaskell 与 Lau 的定义只考虑沿局部坐标方向的相邻点，其条件为
归一化后写成 \(\min(\tilde\phi_C, 1) \le \tilde\phi_f \le \max(\tilde\phi_C, 1)\) 。Gaskell 与 Lau 为隐式稳态流场给出的 CBC 准则要求：格式函数 \(\tilde f(\tilde\phi_C)\) 连续；在 \(0 \le \tilde\phi_C \le 1\) 的单调区里以 \(\tilde\phi_C\) 为下界、以 1 为上界，并经过 \((0, 0)\) 与 \((1, 1)\) 两点；而在 \(\tilde\phi_C < 0\) 或 \(\tilde\phi_C > 1\) 的非单调区里函数值恒等于 \(\tilde\phi_C\)（即 Upwind）。CBC 的物理直觉（图 12.4、12.5）在于：当 \(\tilde\phi_C\) 处在单调剖面内时，界面插值不应在面两侧节点之间引入新的极值，于是被夹在 \(\phi_C\)、\(\phi_D\) 之间；当 \(\tilde\phi_C\) 接近 \(\phi_D\)（仍在单调区）时，\(\phi_f\) 也跟着趋近 \(\phi_D\)，等 \(\tilde\phi_C = 1\) 时 \(\phi_f = \phi_D\)，这就是 \(\tilde\phi_f\) 过 \((1, 1)\) 点的来源。当 \(\tilde\phi_C > 1\) 时，\(\phi_f\) 被强制取 Upwind 值 \(\phi_C\)——这是"在节点之间能达到的、且仍被面两侧节点夹住的最大出流值"，使过大的出流把可能的过冲压低；只要没有源项之类的外部机制产生极值，极值就会被衰减。\(\tilde\phi_C < 0\) 时情况镜像对称：\(\phi_f\) 取 \(\phi_C\)，直到 \(\phi_C = \phi_U\) 时落到 \((0, 0)\)。在 \(\tilde\phi_C < 0\) 或 \(\tilde\phi_C > 1\) 的区域，对流占主导，Upwind 是一个极好的近似。
12.3 高分辨率（HR）方案（High Resolution (HR) Schemes）
用 NVD 来构造一个有界的高阶（HR）方案相对直接：取任意一个高阶基础方案，再用一组事先定好的曲线把它"夹"在 NVD 的有界区域内即可。具体地，单调区 \(0 \le \tilde\phi_C \le 1\) 内的剖面要经过 \((0,0)\) 与 \((1,1)\)、同时落在 NVD 上面的上三角阴影区（图 12.4）之内；非单调区 \(\tilde\phi_C < 0\) 或 \(\tilde\phi_C > 1\) 内则直接采用 Upwind 剖面。图 12.6 展示了一些著名的 HR 格式：MINMOD、BOUNDED CD、OSHER、MUSCL、SMART、STOIC、SUPERBEE。复合 HR 格式在剖面连接点、水平段与垂直段处要避免硬转折以改善收敛行为——例如 SMART（基于 QUICK）把 \(0 \le \tilde\phi_C \le 1/6\) 段改成 \(\tilde\phi_f = 3\tilde\phi_C\)（即 SMART 原式），STOIC 则把 \(0 \le \tilde\phi_C \le 1/5\) 段改写；为了更陡的收敛，SMART 与 STOIC 的水平段也可以微调，例如在 \(9/10 \le \tilde\phi_f \le 1\)（对应 \(7/10 \le \tilde\phi_C \le 1\)）上把剖面改写为 \(\tilde\phi_f = \tilde\phi_C / 3 + 2/3\)，或者在 \(0.95 \le \tilde\phi_f \le 1\) 段对 Bounded CD 也做类似处理。修改后的 SMART、STOIC、SUPERBEE 的 NVD 见图 12.7。
把这些 HR 方案的函数关系写成数学形式，可以列出若干经典格式（为简洁只写主支与"其余"分支）：
文献中按这一思路发展的方案还有 CLAM、UTOPIA、SHARP、ULTRA-SHARP 等。例 3 给出了如何由 QUICK 出发构造 SMART：以 QUICK 直线 \(\tilde\phi_f = 3\tilde\phi_C/4 + 3/8\) 为黑线，再用"经过 \((0,0)\)、\((1,1)\)、落在 NVD 上三角区"的蓝线把过陡、过缓的部分切掉；最终 SMART 的三段函数就是 Modified SMART 形式。
12.4 TVD 框架（The TVD Framework）
总变差减小（Total Variation Diminishing, TVD）是发展 HR 对流格式的另一种主流框架。离散一维空间上节点解 \(\phi_i\) 的总变差（TV）定义为
；一个数值方法称为 TVD 是指解的 TV 不随时间增加，即 \(\mathrm{TV}(\phi^{t+\Delta t}) \le \mathrm{TV}(\phi^t)\) 。Harten 证明单调格式必 TVD、TVD 必极值保持；本节不重述 TVD 的完整数学推导，只按 Sweby 的方法给出 TVD 格式的构造思路。回到第 11 章所用的非稳态一维对流方程（在无扩散、无源时简化为
），其 RHS 用五点模板离散可以写成
Sweby 与 Harten 证明上式所表示的格式 TVD（单调）的充分条件是三个不等式
其中 \(a\)、\(b\) 取决于所采用的具体对流格式。第 11 章已经指出：一阶 Upwind 过度扩散、二阶 CD 高度色散，所需要的是一个介于两者之间、"既稳定又精确"的格式。一个直接的构造是把 CD 写为
即"Upwind + 反扩散通量"——CD 本身色散，所以右端第二项具有反扩散性质，是 CD 二阶精度的来源；副作用是数值扩散被削弱，因而在梯度大处出现非物理振荡。一个更好的办法是只取该反扩散通量的一部分加到 Upwind 上，使二阶精度得以保留而不致振荡。做法是把这一通量乘以一个限制子（limiter / flux limiter）\(\psi(r)\)，其中 \(r\) 通常取相邻梯度之比，得到
为保持反扩散通量的符号不变，\(\psi(r_f)\) 取非负值。构造 TVD 格式就转化为：找一个满足 TVD/单调条件的限制子 \(\psi(r)\)。把 \(\phi_f\) 写成图 12.8 上一维域中单元面上的对流通量形式，分别代入 RHS 后再利用 \(\dot m_e + \dot m_w = 0\)（连续性），把 \(\mathrm{RHS}\) 整理为
与上式比较得
、\(b = 0\)；TVD 条件 \(0 \le a + b \le 1\) 化为
。为对任意相邻单元同时成立，限制子须满足
把三条合起来就得到 Sweby 的 TVD-CBC：\(\psi(r) = \min(2r, 2)\) 当 \(r \ge 0\)，否则 \(\psi(r) = 0\)。这个条件可以在 \(r\)–\(\psi\) 平面上（图 12.9）画成 TVD 单调区；任何落在这个蓝色区域内的限制子都给出 TVD 格式，这就是 Sweby 图（也叫 \(r\)–\(\psi\) 图）。Sweby 图与 NVD 非常类似：把已有 HO 格式写成
的形式就能给出对应限制子。CD 由 \((\phi_D + \phi_C)/2\) 立即得 \(\psi_{CD}(r_f) = 1\)；SOU 由 \(\phi_f = \phi_C - (\phi_C - \phi_U)/2\) 解出 \(\psi_{SOU}(r_f) = r_f\)。两者在 Sweby 图上都过 \((1, 1)\)，二阶精度的限制子也必须过此点（图 12.10）。Van Leer 已证明任何二阶格式都可以写成 CD 与 SOU 的加权平均，所以二阶限制子既过 \((1, 1)\) 又要落在由 CD 与 SOU 两条线围成的区域内（图 12.10 的蓝色区）；这个区域在 NVD 上对应一段从 \((0.5, 0.75)\) 出发、夹在 CD 与 SOU 直线之间的子区域（图 12.11）。Sweby 指出当 \(r_f < 0\) 时 \(\psi(r_f) = 0\)，这意味着解的极值点处会丢失二阶精度。
把以上方法推广，可以得到各 HO 方案的限制子形式：Upwind \(\psi = 0\)、Downwind \(\psi = 2\)、FROMM \(\psi = (1 + r_f)/2\)、SOU \(\psi = r_f\)、CD \(\psi = 1\)、QUICK \(\psi = (3 + r_f)/4\)。FROMM 是 CD 与 SOU 的算术平均；它的 TVD 形式为 \(\psi(r_f) = (1 + r_f)/2\)。除了 Upwind 外，其余限制子都不完全落在单调区内——这正是它们无界的原因。把它们各自截到 TVD 单调区里（图 12.9、12.12），就得到 TVD-HR 格式（图 12.13a–d），函数式为
12.5 NVF 与 TVD 的关系（The NVF-TVD Relation）
NVF 与 TVD 走的是不同的有界化路线，但可以证明二者几乎等价。首先建立 \(r_f\) 与 \(\tilde\phi_C\) 的关系：
即 \(\tilde\phi_C = r_f/(1 + r_f)\)。以 Upwind 限制子 \(\psi(r_f) = 0\) 为例：\(\psi = 0\) 意味着 \(\phi_f = \phi_C\)，所以 \(\tilde\phi_f = \tilde\phi_C\)，与 NVF 中 Upwind 的 \(\tilde\phi_f = \tilde\phi_C\) 完全一致；TVD-CBC 在 \(r_f < 0\) 时强制 Upwind，对应到 NVF 条件是 \(\tilde\phi_C < 0\) 或 \(\tilde\phi_C > 1\)（分母为负或分子大于分母）；NVF-CBC 要求在 \(0 \le \tilde\phi_C \le 1\) 内函数单调递增，Sweby 图上对应 \(0 \le r_f < +\infty\)，而当 \(\tilde\phi_C \to 1\) 时 \(r_f \to +\infty\)，两区间确实重合。TVD-CBC 的另一条 \(\psi(r_f) \le 2\) 等价于 \(\phi_f = \phi_C + (\phi_D - \phi_C) = \phi_D\)，即 \(\tilde\phi_f = 1\)，与 NVF-CBC 一致。TVD-CBC 中 \(\psi(r_f) \le 2r_f\) 等价于
归一化得 \(\tilde\phi_f = 2\tilde\phi_C\)，比 NVF-CBC 在该处的条件更严格——这是两种格式之间唯一真正的差异（图 12.14a：TVD-CBC 的单调区被压缩到 Upwind 线与蓝色区）。把 NVF-CBC 中 \(\tilde\phi_f = 0\) 对应到 TVD 上相当于 \(r_f = 0\)（图 12.14b）。对二阶精度，TVD 要求 \(w(1) = 1\)；代入 \(r_f = 1\) 得 \(\tilde\phi_C = 1/2\)；再代回 \(\phi_f\) 公式得 \(\tilde\phi_f = 3/4\)——正好是 NVD 上的 \(Q(0.5, 0.75)\)。Van Leer 已经证明任何二阶格式都是 CD 与 SOU 的加权平均，因此二阶 TVD 格式的 NVD 必然落在 CD 与 SOU 直线所夹的蓝色区内（图 12.14c）。
把以上对应推广，可以从任一 TVD 格式得到等价的 NVF 格式，反之亦然。从 NVF 框架出发，把 \(\phi_f\) 写成
而 TVD 框架给出
令两式相等，整理得
化简后写成
用这个变换可以验证：Upwind 限制子 \(\psi = 0\) 对应 NVF 的 \(\tilde\phi_f = \tilde\phi_C\)；Downwind 限制子 \(\psi = 2\) 对应 \(\tilde\phi_f = (2 + 2r_f)/[2(1 + r_f)] = 1\)；SOU 限制子 \(\psi = r_f\) 在 NVF 中是 \(\tilde\phi_f = 3\tilde\phi_C/2\)。例 4 是 Van Leer 限制子 \(\psi(r_f) = (r_f + |r_f|)/(1 + |r_f|)\) 转换到 NVF 下的推导：分 \(0 \le \tilde\phi_C < 1\)、\(\tilde\phi_C \ge 1\)、\(\tilde\phi_C < 0\) 三种情况，得到 Van Leer 在 NVF 下的最终形式是
12.6 非结构网格中的 HR 方案（HR Schemes in Unstructured Grid Systems）
第 11 章提到过非结构网格上构造 HR 方案的一个常见困难：没有清晰的"上游节点 U"，而 \(\tilde\phi_C\) 或 \(r_f\) 的计算又需要 U。一个简单的做法是人为造一个虚拟 U：让 U 落在 C、D 两点的连线上，并且令 C 恰好是 U、D 段的中点（图 12.15）。于是 \(\phi_D - \phi_U = 2 r \phi_C \cdot d_{UD} = 2 r \phi_C \cdot d_{CD}\) ，从而
其中 \(d_{CD}\) 是 C 到 D 的向量，\(d_{UD}\) 是 U 到 D 的向量。算出 \(\phi_U\) 之后，无论是 NVF 还是 TVD 流程都可以按之前讨论的方式进行。
12.7 HR 方案的延迟修正（Deferred Correction for HR Schemes）
HR 方案的数值实现最好用一个例子来展示。考虑带源项的多维对流方程 \(\nabla \cdot (\rho \mathbf{v} \phi) = Q_\phi\)，已知速度场。在图 12.16 的体积 \(V_C\) 上积分、应用散度定理、再把面积分化为对单元各面的求和，得到
把面 \(f\) 上的质量流率记为 \(\dot m_f = (\rho \mathbf{v})_f \cdot \mathbf{S}_f\)，方程化为
\(\phi_f\) 由前面任一对流格式给出，但为使代数方程对主节点上的未知数可解，离散方程应该形如
把 \(\phi_f\) 直接用相邻节点值表示会带来不稳定，下一小节会说明。
12.7.1 直接用节点值带来的困难（The Difficulty with the Direct Use of Nodal Values）
把 TVD 格式直接套到对流通量上（图 12.16 的一维化形式）会得到
其中 \(\dot m_f^+ = \max(\dot m_f, 0)\)、\(\dot m_f^- = -\min(\dot m_f, 0)\)。代回离散方程
可得
退化到图 12.8 的一维情形并取正流向，方程化为
其中
由于 \(0 \le \psi(r) \le 2\)，\(a_E\) 与 \(a_W\) 的符号相反（Upwind 情形 \(\psi = 0\) 除外），违反了稳定性的基本规则之一，使迭代收敛困难。NVF 走同样的路也会遇到同一种缺陷。第 11 章介绍过一种补救——延迟修正（Deferred Correction, DC）程序：代数方程中的系数基于 Upwind，而 HR 与 Upwind 之间的差异作为源项加入右端。DC 实现简单、结构与非结构网格上都能用，但 Upwind 与 HR 面值的差异越大收敛越慢；这一现象可以在 NVD 上直观地看出来——Upwind 线与所选 HR 格式线之间的归一化差就是收敛速度的"距离"，差越大收敛越慢。为此研究者发展了更隐式的实现方法，下一节介绍其中两个。
12.8 DWF 与 NWF 方法（The DWF and NWF Methods）
DC 程序的代价是 HR 方案实现中的收敛速度下降。本节介绍两种为绕开这一代价而设计的实现：Leonard 与 Mokhtari 的下游权重因子法（Downwind Weighing Factor, DWF，OpenFOAM® 中所采用）与 Darwish 与 Moukalled 的归一化权重因子法（Normalized Weighing Factor, NWF）。两种方法都在图 12.16 的三维非结构网格上、求解守恒方程
这样的背景下展开；面值 \(\phi_f\) 由 HR 格式给出，下面的讨论关心的是如何把这些面值"最有效"地塞进离散方程。
12.8.1 下游权重因子（DWF）法（The Downwind Weighing Factor (DWF) Method）
DWF 的定义是
利用它可以把面值改写为
即把 HR 格式的估计 \(\phi_f\)（或归一化值 \(\tilde\phi_f\)）按权重分摊到 Upwind 和 Downwind 两端，效果是缩小离散系数的模板。HR 格式给出的 \(\phi_f\) 一定落在 \(\phi_C\) 与 \(\phi_D\) 之间，所以 \(\text{DWF}_f\) 始终在 0 与 1 之间。表 12.1 列出 Upwind、SOU、CD、FROMM、QUICK、Downwind、MINMOD、Bounded CD、OSHER、SMART、STOIC、MUSCL 等若干 HO/HR 格式在均匀网格上 \(\text{DWF}_f\) 的函数关系——例如 Upwind 给 \(\text{DWF}_f = 0\)、CD 给 \(\text{DWF}_f = 1/2\)、QUICK 给 \(\text{DWF}_f = 1/4 + 1/8(1 - \tilde\phi_C)^{-1}\) 等。比较 TVD 形式
与上式可以立即看出
但上节已经论证，这种实现会破坏对角占优，因此格式不稳定；这正是 OpenFOAM® 在求解对流主导问题时经常遇到数值困难的根源。
为完整起见，仍按 DWF 推导一般方程。把面值写成
其中 \(\text{DWF}_f^+ = (\phi_f - \phi_C)/(\phi_F - \phi_C)\) 、 \(\text{DWF}_f^- = (\phi_f - \phi_F)/(\phi_C - \phi_F)\) 。于是
离散方程
中
退化到一维（\(\dot m_w = -\dot m_e\)）并取正流，得到
\(a_E\) 与 \(a_W\) 符号相反，又一次违反基本系数规则。当 \(\text{DWF} > 0.5\)（即 \(\phi_f > 0.5(\phi_C + \phi_D)\)）时，\(a_C\) 还会变负，迭代根本无法进行；这种情况在所有 HR 格式下只要 \(\tilde\phi_C > 0.5\) 都会出现——DWF 把大量 HR 通量"压"到 Downwind 一端，行为上很像中心差分。
12.8.2 归一化权重因子（NWF）法（The Normalized Weighing Factor (NWF) Method）
NWF 正是为弥补 DWF 的缺陷而提出的。它的思路是把归一化插值剖面线性化为
其中 \(\ell\)、\(k\) 是各区间上的常数（斜率与截距），区间数取决于所用的 HR 方案。这对几乎所有 HR 格式都是精确的。以 MINMOD 为例，把 \(\tilde\phi_f = \ell \tilde\phi_C + k\) 与 MINMOD 的 NVF 形式对比，可得
再把上式改写为
整理得
其中 \(\phi_U\)、\(\phi_D\)、\(\phi_C\) 的位置取决于流动方向。表 12.2 给出 Upwind、SOU、CD、FROMM、QUICK、MINMOD、OSHER、MUSCL、SMART 等 HO/HR 方案在均匀网格（非结构/结构）上的 \([\ell, k]\) 值。在非结构网格上 U 是虚拟节点，\(\phi_U\) 这一项按延迟修正处理；其修正量 \((1 - \ell - k)\phi_U\)（其中 \(\ell < 0\) 或 \(k < 0\)）的绝对值小于标准 DC 中直接用 \(\phi_U\) 时的值，所以 NWF 比标准 DC 需要的欠松弛更少、收敛更快。
把上式代入对流通量，并在正反两个方向上同时展开，得到
代回代数方程
得到
其中
。NWF 原本是为结构网格发展的，在结构网格上 U 是真实节点，可以被直接解出，于是方程的模板扩大为包含 EE 和 WW 等远节点——在图 12.8 的一维结构网格上，代数方程展开为
其中
按 NWF 的重写，由表 12.2 可见在大多数区间上 \(\ell > k\)，因此 \(a_C\) 始终为正、不会失稳；只有在 Downwind 线附近 \((\ell, k) = (0, 1)\) 时 \(a_C\) 才为零，此时把 \((\ell, k)\) 改设为 \((L, 1 - L)\)，其中 \(L\) 通常取前一个区间的 \(\ell\) 值。这使 NWF 比 DWF 健壮得多——它能保证 \(a_C\) 始终为正。
12.8.2.1 TVD 框架下的 NWF 方法（The NWF Method in the Context of the TVD）
除 MUSCL（Van Leer 限制子）以外，前面所有 HR TVD 格式的限制子在 Sweby 图上都呈直线，因此可以统一写成
其中 \(m\)、\(n\) 是只与几何量有关的常数，在 \(\psi(r_f)\) 的每个区间上取一组值；区间数取决于具体的 HR TVD 方案。例如与 MINMOD（其 TVD 形式为 \(\psi(r_f) = \max(0, \min(1, r_f))\)）对比，可得
把 \(\psi(r_f) = m r_f + n\) 代入 TVD 插值公式
整理为
其中 \(\phi_U\)、\(\phi_D\)、\(\phi_C\) 的位置仍取决于流向。表 12.3 列出 Upwind、SOU、CD、FROMM、QUICK、OSHER、MUSCL、SUPERBEE 等格式在均匀网格上的 \([m, n]\) 值。把上式与 NVF-NWF 的形式 \(\phi_f = \ell \phi_C + k \phi_D + (1 - \ell - k)\phi_U\) 对照，可得
因此可以照搬 NVF-NWF 的实现思路来构造 TVD-NWF。
12.9 边界条件（Boundary Conditions）
对流项的边界条件一般比扩散项要简单。本节按 Inlet、Outlet、Wall、Symmetry 四种类型分别给出实现细节。典型边界单元见图 12.17：边界单元 C 有一个面位于边界上，该面形心记为 b，外法向面积矢量为 \(\mathbf{S}_b\)。在多维纯对流问题的边界单元上离散得到
内部面通量按前述格式离散；无论哪种边界条件类型，边界通量都可以用边界面值写成
因此边界单元的离散方程是
其中下标 \(f\) 指内部面，\(b\) 指边界面。具体边界条件就是指定边界上的 \(\phi_b\) 或者边界通量 \(\mathbf{J}_{\phi, C}|_b\)。
12.9.1 入口边界条件（Inlet Boundary Condition）
入口处（图 12.18）\(\phi\) 通常被指定；速度场已知，对流通量也就已知。把该项移到方程右端作为源项，(12.107) 化为
若内部面对流通量用 HR 格式（以 DC 实现），则边界单元的代数方程为
其中
这里 \(F\) 指 C 点的内部邻居，\(f\) 指边界单元的内部面。
12.9.2 出口边界条件（Outlet Boundary Condition）
出口处（图 12.19）下游没有信息，但作为方向性现象，\(\phi\) 在边界上的值强烈依赖于上游——Upwind 和 SOU 格式就不需要任何下游信息，因为面值已经是上游节点的函数。一种最有效的处理是假设 \(\phi\) 在出口充分发展（fully developed），即 \(\phi\) 在边界法向上的梯度为零
；常见做法是在出口采用 Upwind \(\phi_b = \phi_C\)，自动得到零法向梯度。HR 格式在内部面上以 DC 实现时，出口单元的代数方程为
其中
这里 \(f\) 指边界单元的内部面，\(C\)、\(F\) 分别指 owner 与 neighbor。
12.9.3 壁面边界条件（Wall Boundary Condition）
如图 12.20，壁面上法向速度为零，因此对流通量为零、不出现在代数方程中。HR 格式以 DC 实现时，壁面单元的代数方程为
其中
这里 \(f\) 仍指内部面，\(C\)、\(F\) 仍指 owner 与 neighbor。
12.9.4 对称边界条件（Symmetry Boundary Condition）
对称边界上没有穿越边界的流动，处理方式与壁面边界类似——把法向对流通量设为零即可。
12.10 计算实现提示（Computational Pointers）
12.10.1 uFVM
与第 11 章的 HO 格式类似，uFVM 中 HR 格式通过延迟修正（DC）方法实现。实现上的主要差别是用 NVF 或 TVD 关系式代替直接用梯度来计算面值。Listing 12.1（cfdAssembleConvectionTermDCSTOIC）展示了用 NVF 形式实现 STOIC HR 格式的过程：先取所需场、再为所有内部面设定 Upwind/Downwind 索引；接着对每个面识别 \(\phi_C\)、\(\phi_D\)，构造 \(\phi_U\)，计算 \(\tilde\phi_C\) 并由 HR 格式的 NVF 关系式得到 \(\tilde\phi_f\)；最后由 \(\tilde\phi_f\) 重建面值 \(\phi_f = \tilde\phi_f(\phi_D - \phi_U) + \phi_U\) 供 DC 程序使用。Listing 12.1 的关键 MATLAB 片段包括：
- 取出面质量流率 \(\dot m_f\)、owner/neighbor 索引、Upwind/Downwind 标志；
- 由虚拟上游
构造 \(\phi_U\)； - 算 \(\tilde\phi_C = (\phi_C - \phi_U)/(\phi_D - \phi_U)\)； - 分段赋值：\(\tilde\phi_C \le 0\) 走 Upwind、\(0 < \tilde\phi_C \le 0.2\) 走 \(3\tilde\phi_C\)、\(0.2 \le \tilde\phi_C < 0.5\) 走 \(0.5\tilde\phi_C + 0.5\)、\(0.5 \le \tilde\phi_C < 5/6\) 走 \(0.75\tilde\phi_C + 3/8\)、\(5/6 \le \tilde\phi_C < 1\) 走 1、\(\tilde\phi_C \ge 1\) 走 Upwind； - 重构 \(\phi_f = \tilde\phi_f(\phi_D - \phi_U) + \phi_U\)； - 修正通量 \(\text{corr} = \dot m_f (\phi_f - \phi_C)\)，并以 DC 方式叠加到原通量上。
12.10.2 OpenFOAM®
如第 11 章所述，OpenFOAM® 通过基类 surfaceInterpolationScheme 完成对流项离散；HR 格式都是特殊的面插值算法。OpenFOAM® 把所有 TVD 方案统一在 limitedSurfaceInterpolationScheme 类中（图 12.21），它继承自 surfaceInterpolationScheme。如 Listing 12.2 所示，该类主要包含三个成员函数：表示基类虚函数的 weights、带三个参数的局部 weights、以及新的虚基函数 limiter。基类虚 weights（Listing 12.3）内部会调用局部 weights，把通量 phi、CD 权重、限子返回值作为参数。局部 weights（Listing 12.4）的实现是
pWeights[face] = pWeights[face] * CDweights[face] + (1.0 - pWeights[face]) * pos(faceFlux_[face]);
看似不直观，但当通量为正时化简为
\(\text{weight} = \text{CDweight} \cdot \psi + (1 - \psi) = 1 - \psi/2\)
（均匀网格 \(\text{CDweight} = 1/2\)），代回第 11 章的 \(\phi_f = \phi_N + \text{weight}(\phi_O - \phi_N)\) 正好得到 TVD 公式 (12.30)。OpenFOAM® 把限子组织为 LimitedScheme 类（Listing 12.5）——它通过 Limiter 模板参数继承自 limitedSurfaceInterpolationScheme，并把 limiter（Listing 12.6）特化为虚函数。其 limiter 函数体（Listing 12.7）实际上调用了辅助函数 calcLimiter，后者负责为 TVD 限子计算 \(\psi\) 值并存入 tlimiterField。calcLimiter（Listing 12.8）依次：把被插值场保存为 lPhi、用 fvc::grad(lPhi) 计算其梯度、取出 CD 权重、收集单元中心坐标、最后调用嵌套模板类 Limiter::limiter 计算每个面的限子值
\(\text{pLim}[\text{face}] = \text{Limiter::limiter}(\cdots)\)
。
以 SUPERBEE 限子函数（Listing 12.9）为例，它接收 CD 权重、流通量、owner/neighbor 上的 \(\phi\) 与梯度、以及两点间向量 \(\mathbf{d}\)，内部先用 LimiterFunc::r 算出 \(r\)，再返回 \(\max(\max(\min(2r, 1), \min(r, 2)), 0)\)——与公式 (12.44) 完全一致。\(r\) 的定义由 Listing 12.10 给出：当 \(\dot m_f > 0\) 时取
、否则取
，最终 \(r = 2(\text{grad}_f / \text{grad}_f^{\phi}) - 1\)（\(\text{grad}_f^{\phi} = \phi_F - \phi_C\)）——这与第 11 章 (12.66) 一致。
12.11 小结（Closure）
本章讨论 HO 对流格式的有界化问题：通过施加对流有界性准则（CBC）把 HO 格式约束到有界范围内，所得的 HO 有界格式记为 HR 格式。归一化变量格式（NVF）与总变差减小（TVD）是发展 HR 格式的两大框架——尽管走的是不同路径，但可以证明二者几乎等价。本章还介绍了两种在结构与非结构网格上实现 HO/HR 格式的技术：下游权重因子法（DWF）与归一化权重因子法（NWF）。下一章将转向非稳态项的离散化。
12.12 习题（Exercises）
本章习题均为源文给读者的练习题（NVF/TVD 格式互相转换、非均匀网格上的 \(Q\) 点、OSHER/SMART 的 DWF/NWF 推导与一维数值实验、Smith-Hutton 与斜流场步态廓形问题的 OpenFOAM®/uFVM 求解、Doxygen 派生类列举等），按精读规约第 13 条不纳入内容概述。
本章个人批注
第 12 章在结构上把"什么是有界 HR 格式"切成两条路——NVF 与 TVD——然后再论证二者实质上等价，这条论述是整章的脊柱。我自己读下来的几个关键心得。
第一，NVD 上所有二阶格式都必过 \(Q(1/2, 3/4)\)——这并不是一个偶然的几何现象，而是 Taylor 展开后 \(b = (c - a)/3\) 这一个约束在均匀网格上的反映。例 2 把这条逻辑展开得相当清楚，给出了一个"看图说话"格式背后真正的代数理由。QUICK 因为斜率正好是 0.75 所以是三阶精度的，这是一个非常干净的"二阶必须过 Q、过 Q 还不够、要斜率 0.75 才三阶"的几何论据。
第二，TVD 单调区的推导看上去很技术（把 Sweby 的 \(r\)–\(\psi\) 图一层层剥开），但物理直觉其实只有一句话：在梯度大处少用反扩散通量、在光滑区多保留。各种限制子（MINMOD、OSHER、MUSCL、SUPERBEE、Van Leer）的差异就在于"少用"和"多用"之间的折中。SUPERBEE 最激进（最接近 Downwind），MINMOD 最保守（最接近 Upwind），Van Leer 和 MUSCL 居中。例 4 把 Van Leer 限制子翻译成 NVF 形式 \(\tilde\phi_f = 2\tilde\phi_C^2\) 的过程特别有助于把这两个框架在脑子里对齐。
第三，DWF 与 NWF 这条线索本质上解决的是"在迭代求解时，HR 格式怎么塞进代数方程而不破坏对角占优"。DWF 的失败很说明问题：\(\text{DWF}_f\) 把 HR 通量分到 Upwind/Downwind 两端，但分给 Downwind 那一份导致 \(a_E\) 与 \(a_W\) 符号相反——这是结构性的，不可能靠调参解决。NWF 用 \([\ell, k]\) 的分解把 \(\phi_f\) 写成 \(\ell\phi_C + k\phi_D + (1 - \ell - k)\phi_U\)，再让 \(\phi_U\) 这一项以 DC 方式进右端，关键观察是 \(\ell\) 在大多数区间上都大于 \(k\)，所以 \(a_C\) 始终为正（只有 Downwind 线附近 \((\ell, k) = (0, 1)\) 时 \(a_C\) 才为零，因此用前一个区间的 \(\ell\) 值做修补）。这套思路与 Sweby 图上的 TVD 单调区是相通的——都是在不破坏对角占优的前提下让格式尽量激进。
第四，虚拟 U 节点的处理是 12.6 节里最值得记的一点：非结构网格没有天然的 U，于是让 U 落在 C–D 延长线上并取 C 为中点，得到
。这样所有 NVF/TVD 的工具都可以照搬，但代价是 \(U\) 实际上是 C 的"延长"，其精度依赖于 \(\nabla\phi_C\) 的计算精度——这又把第 9 章的梯度计算与第 12 章的 HR 格式串起来了。
最后，边界条件那一节实际是把内部面公式套到边界单元上时对 \(\dot m_f\)、\(\phi_f\)、\(\text{Flux}\) 的重新安排。Inlet 直接把 \(\dot m_b \phi_b\) 移到右端，Outlet 约定 \(\phi_b = \phi_C\)（充分发展），Wall 和 Symmetry 直接置对流通量为零——读到这里我对这套离散化在边界上的处理才真正踏实下来。
与上下章的衔接（一段话）
第 12 章在全书的位置是"对流项"的第二阶段：第 11 章建立了对流项的离散化基础（CD 失败、Upwind 的数值扩散、SOU/QUICK/FROMM 等 HO 格式、Leonard 稳定性准则、DC 程序），第 12 章则把所有 HO 格式提升为有界的 HR 格式——通过 NVF/NVD 或 TVD/Sweby 图两条等价的路径，把"在光滑区用 HO、在梯度大处退化到 Upwind"这个核心思想表达成可操作的数学形式，并在 12.6–12.8 节解决了"如何在结构/非结构网格上把这些格式实际装进代数方程且不破坏对角占优"的问题（DWF → NWF）。第 12 章末尾给出的 OpenFOAM®/uFVM 实现提示（12.10 节）也对应着第 11 章同节的对流项基类，为下一章的时间项离散化做了一个清晰的收尾——第 13 章将把迄今为止的稳态讨论拓展到非稳态对流问题，时间项的离散化会让 \(\phi\) 随时间推进，于是 HR 格式的"不破坏对角占优"性质会在非稳态迭代中变得更加关键。