第 11 章：对流项的离散化（Discretization of the Convection Term）

11.1 引言（Introduction）

本章开始处理守恒方程中的对流项（即散度算子 \(\nabla\cdot(\rho\mathbf{v}\phi)\)）。在已经建立正交、非正交、结构与非结构网格上一般稳态扩散项离散化方法的基础上，作者转向这个"看似简单"但实际上带来很大麻烦的项。

对流项离散化所遭遇的困难在过去三十多年里一直是被研究的主题，相关文献数量庞大，因此本书用整整两章来讨论它：当前这一章处理基础概念以及高阶（High Order, HO）迎风偏置（upwind biased）格式；下一章再讨论通过有界化（bounding）对流通量而发展出来的高分辨率（High Resolution, HR）非振荡高阶格式。

为叙述清晰起见，新的概念先在一维网格上引入，然后再推广到多维非正交网格。本章首先对一维对流-扩散问题进行离散化，借助一个稳定性判据指出中心差分（Central Difference, CD）格式的缺陷；接着介绍能通过该稳定性测试的迎风格式（Upwind Scheme）；同时解释伴随低阶格式出现的数值扩散（numerical diffusion）误差以及伴随高阶格式出现的数值色散（numerical dispersion）误差。高分辨率（HR）格式族留待下一章处理。

11.2 稳态一维对流与扩散（Steady One Dimensional Convection and Diffusion）

为避免计算的复杂性淹没主要物理思想，作者首先考察一个非常简单的稳态一维对流-扩散问题。所讨论的守恒方程为

\[ \frac{\mathrm{d}(\rho u \phi)}{\mathrm{d} x} - \frac{\mathrm{d}}{\mathrm{d} x}\left(\Gamma_\phi \frac{\mathrm{d}\phi}{\mathrm{d} x}\right) = 0. \]

幸运的是该问题存在解析解，因此可作为参照基准，用来比较各种数值离散化方案的精度。

11.2.1 解析解（Analytical Solution）

对于常截面积的一维稳态问题，连续性方程 \(\mathrm{d}(\rho u)/\mathrm{d} x=0\) 意味着 \(\rho u\) 为常数。记这个常数为 \(\dot m = \rho u\)，则守恒方程可积分为 \(\dot m \phi - \Gamma_\phi\, \mathrm{d}\phi/\mathrm{d} x = c_1\) ，其中 \(c_1\) 为依赖于边界条件的积分常数。改写为 \(\mathrm{d}\phi/\mathrm{d} x = (\dot m/\Gamma_\phi)\phi - c_1/\Gamma_\phi\) ，再作变量替换 \(U = \phi - c_1/\dot m\) 使之化为 \(\mathrm{d}U/\mathrm{d} x = (\dot m/\Gamma_\phi)U\)，分离变量后积分得到

\[ \phi(x) = \frac{c_2\, \Gamma_\phi\, \exp(\dot m x/\Gamma_\phi) + c_1}{\dot m}. \]

把两点 W、E 处的边界值 \(\phi_W\)、\(\phi_E\) 代入，并定义佩克莱数 \(P_{eL} = \dot m L/\Gamma_\phi\)（其中 \(L = x_E - x_W\)）作为对流输运率与扩散输运率之比，最终的解析解为

\[ \frac{\phi - \phi_W}{\phi_E - \phi_W} = \frac{\exp\!\big(P_{eL}\,\tfrac{x-x_W}{L}\big) - 1}{\exp(P_{eL}) - 1}. \]

作者指出，对不同的 \(P_{eL}\) 值，\(\phi\) 在 W、E 之间的形状从纯扩散时的近似直线，逐渐过渡到高 \(P_{eL}\) 时的近似阶跃（step）轮廓。

11.2.2 数值解（Numerical Solution）

数值离散化从把守恒方程在所示一维单元上做体积分开始：

\[ \int_{V_C} \big[\nabla\cdot(\rho\mathbf{v}\phi) - \nabla\cdot(\Gamma_\phi\nabla\phi)\big]\,\mathrm{d} V = 0. \]

把对流与扩散通量分别记为 \(\mathbf{J}_{\phi,\mathrm{C}} = \rho\mathbf{v}\phi\) 与

\[ \mathbf{J}_{\phi,\mathrm{D}} = -\Gamma_\phi\nabla\phi \]

，用散度定理把体积分转化为表面积分，再用面矢量在相对两侧互为反向的事实把面积分化为对面通量的求和。对于常截面积的情况，结果为

\[ (\dot m\Delta y\,\phi)_e - \Gamma_\phi\frac{\mathrm{d}\phi}{\mathrm{d} x}\bigg|_e\,\Delta y \;-\; \big[(\dot m\Delta y\,\phi)_w - \Gamma_\phi\frac{\mathrm{d}\phi}{\mathrm{d} x}\big|_w\,\Delta y\big] = 0. \]

在面 \(e\)、\(w\) 上的速度为已知，\(\mathrm{d}\phi/\mathrm{d} x\) 沿用之前章节的离散化方法，关键问题在于如何由相邻节点的 \(\phi\) 值表达出 \(\phi_e\)、\(\phi_w\)。这一选择面值的做法在文献中被称为"对流格式"（advection scheme）。

11.2.3 预备推导：中心差分（CD）方案（A Preliminary Derivation: The Central Difference Scheme）

"最直观"的选择是采用与扩散项处理相同的线性插值型分布：假定 \(\phi(x) = k_0 + k_1(x-x_C)\)，用穿过面 \(e\) 的两个节点 \(C\)、\(E\) 拟合，得到

\[ \phi_e = \phi_C + \frac{\phi_E - \phi_C}{x_E - x_C}(x_e - x_C). \]

在均匀网格上化简为 \(\phi_e = (\phi_C + \phi_E)/2\)。这就是中心差分（CD）方案：本质上是把 \(\phi\) 在节点 \(x_C\) 附近做 Taylor 展开、舍去二阶及以上项得到的，因此具有二阶精度。

接着把 \(\phi_e = (\phi_C+\phi_E)/2\) 代入一维对流-扩散离散化方程，与相应的 \(\phi_w\) 表达式合并，可以把方程改写为 \(a_C\phi_C + a_E\phi_E + a_W\phi_W = 0\)，其中（对一维常截面积、连续性已保证 \(\dot m_e - \dot m_w = 0\)、并设 \(\Gamma_\phi\) 为常数的情形）

\[ a_E = -\frac{\Gamma_\phi}{x_E - x_C} + \frac{\dot m_e}{2},\quad a_W = -\frac{\Gamma_\phi}{x_C - x_W} - \frac{\dot m_w}{2},\quad a_C = -(a_E+a_W). \]

由此解出 \(\phi_C\)，并用 \(P_{eL}\) 表达为

\[ \frac{\phi_C - \phi_W}{\phi_E - \phi_W} = \frac{1}{2}\Big(1 - \frac{P_{eL}}{2}\Big). \]

将该 CD 数值解与解析解在 \(P_{eL}\) 从 \(-10\) 到 \(+10\) 的范围内比较（Fig. 11.4）：低 \(P_{eL}\) 时两者接近；但当 \(P_{eL}\) 增大到一定程度时，CD 解严重偏离解析解并变得无界——解析解在 \(P_{eL}\to\pm\infty\) 时分别趋近 0 和 1，CD 解却从 \(+\infty\) 单调线性地降到 \(-\infty\)，出现明显的非物理行为。

作者用扩散与对流在"影响域"（zone of influence）上的差异来解释这一现象：扩散是各向同性的，节点 C 的扩散同时受上下游条件的影响；而对流是高度方向性的，只沿流动方向输运性质。线性对称分布对扩散是合理的近似，但对对流而言却忽视了方向偏好——在大 \(P_{eL}\)、对流主导的情形下，这种忽视导致非物理结果。数学上，"有界性"破坏发生于 \(a_E \ge 0\)（流动沿 \(+x\) 时；沿 \(-x\) 时破坏点出现在 \(a_W \ge 0\)），这要求单元佩克莱数 \(P_e = \dot m\,\Delta x/\Gamma_\phi\) 满足

\[ P_e \ge 2. \]

减小网格使 \(P_e < 2\) 即可避免此问题，但很多实际问题中所需的存储和计算代价过大；对纯对流问题（如 Euler 流）则根本不可行。因此需要一个直接的修正方案。

11.2.4 迎风方案（The Upwind Scheme）

作者指出，造成正系数 \(a_E\) 的根源是所用的线性对称分布：它赋予面的两侧节点相同的权重、无方向偏好，适用于扩散这样的椭圆型项，但对流项并非如此。

与对流物理更一致的方案是迎风格式（Fig. 11.6）：单元面值取自流动方向上"上游"节点的 \(\phi\) 值，即

\[ \phi_e = \begin{cases} \phi_C & \dot m_e \ge 0 \\ \phi_E & \dot m_e < 0 \end{cases},\qquad \phi_w = \begin{cases} \phi_C & \dot m_w \ge 0 \\ \phi_W & \dot m_w < 0 \end{cases}. \]

把 \(\dot m_e\phi_e\)、\(\dot m_w\phi_w\) 写成"通量格式" \(a_C\phi_C + a_E\phi_E + a_W\phi_W\) 的形式（使用 \(\max(\cdot,\cdot)\) 记号），再叠加扩散贡献，可得

\[ a_E = -\max(-\dot m_e, 0) - \Gamma_\phi \frac{S_e}{d_{xe}},\quad a_W = -\max(-\dot m_w, 0) - \Gamma_\phi \frac{S_w}{d_{xw}}, \]

而

\[ a_C = \sum_f \mathrm{Flux}_{C_f} = -(a_E + a_W) + (\dot m_e + \dot m_w) \]

；连续性 \(\dot m_e + \dot m_w = 0\) 时，\(a_C = -(a_E+a_W)\)。作者指出，迎风格式保证邻居系数为负、主系数 \(a_C = -(a_W+a_E)\) 满足有界性条件。

在均匀网格、常扩散系数的假设下，\(\phi_C\) 的显式表达式为

\[ \frac{\phi_C - \phi_W}{\phi_E - \phi_W} = \frac{2 + \max(-P_{eL}, 0)}{4 + \max(-P_{eL}, 0) + \max(P_{eL}, 0)} = \frac{2 + \max(-P_{eL}, 0)}{4 + |P_{eL}|}. \]

与解析解及 CD 数值解的对比（Fig. 11.7）显示：低 \(P_{eL}\) 时迎风格式的精度不如 CD 格式（因为它是一阶精度而 CD 是二阶），但在高 \(P_{eL}\) 时，CD 解无界、物理上错误，迎风格式则保持物理上正确。因此迎风格式在精度与稳定性之间作出权衡：解的行为在所有佩克莱数下都有界，代价是低精度；而 CD 方案在某个 \(P_{eL}\) 之后失稳。两种方案都被误差所"感染"，只是表现不同。下一节在介绍下风格式之后，作者将对这些"误差"作出更具体的解释。

11.2.5 下风方案（The Downwind Scheme）

为对照，作者介绍下风格式：面 \(\phi\) 值取自流动下风侧节点的 \(\phi\)（即 Fig. 11.8 中下风节点代表面值），形式化为

\[ \phi_e = \begin{cases} \phi_E & \dot m_e \ge 0 \\ \phi_C & \dot m_e < 0 \end{cases},\qquad \phi_w = \begin{cases} \phi_W & \dot m_w \ge 0 \\ \phi_C & \dot m_w < 0 \end{cases}. \]

按同样方法把对流与扩散通量合并后得到系数

\[ a_E = \max(\dot m_e, 0) - \Gamma_\phi \frac{S_e}{d_{xe}},\quad a_W = \max(\dot m_w, 0) - \Gamma_\phi \frac{S_w}{d_{xw}}, \]

而 \(a_C = -(a_E+a_W) + (\dot m_e+\dot m_w)\)，连续性条件下 \(a_C = -(a_E+a_W)\)。

作者指出，\(\phi_C\) 沿网格的显式表达式变为

\[ \frac{\phi_C - \phi_W}{\phi_E - \phi_W} = \frac{2 - \max(P_{eL}, 0)}{4 - \max(-P_{eL}, 0) - \max(P_{eL}, 0)} = \frac{2 - \max(P_{eL}, 0)}{4 - |P_{eL}|}, \]

当 \(|P_{eL}| \to 4\) 时完全无界，再次验证了上述分析。下风格式单独使用意义不大，但在与其他格式混合预测锐利界面时可能有用。本节引入下风格式的主要目的是为下一节关于稳定性的讨论提供额外的洞见。

11.3 截断误差：数值扩散（Truncation Error: Numerical Diffusion and Anti-Diffusion）

由于离散化过程是近似的，必然引入截断误差；这部分在一维笛卡尔网格上分析最为直接。下面分别给出迎风、下风和中心差分三种格式的扩散/反扩散特性。

11.3.1 迎风方案（The Upwind Scheme）

设流动沿 \(+x\) 方向，迎风格式给出 \(\phi_e = \phi_C\)、\(\phi_w = \phi_W\)。对 \(\phi_C\) 关于面 \(e\) 的 \(\phi_e\) 做 Taylor 展开（均匀网格）：

\[ \phi_C = \phi_e - \frac{\mathrm{d}\phi}{\mathrm{d} x}\bigg|_e \frac{\Delta x}{2} + \frac{1}{2}\frac{\mathrm{d}^2\phi}{\mathrm{d} x^2}\bigg|_e \left(\frac{\Delta x}{2}\right)^2 + \cdots \]

对 \(\phi_W\) 关于 \(\phi_w\) 做类似展开。把这些二阶以上项截去、带入离散化方程的左边，并整理可得

\[ \dot m_e\phi_C - \dot m_w\phi_W - \Gamma_\phi\frac{\mathrm{d}\phi}{\mathrm{d} x}\bigg|_e \Delta y - \Gamma_\phi\frac{\mathrm{d}\phi}{\mathrm{d} x}\bigg|_w \Delta y = \dot m_e\phi_e - \dot m_w\phi_w - \left(\Gamma_\phi + \frac{\dot m\,\Delta x}{2}\right)\frac{\mathrm{d}\phi}{\mathrm{d} x}\bigg|_e \Delta y - \left(\Gamma_\phi + \frac{\dot m\,\Delta x}{2}\right)\frac{\mathrm{d}\phi}{\mathrm{d} x}\bigg|_w \Delta y. \]

显式地，迎风格式离散化时实际求解的方程多出一个额外扩散项，其扩散系数增量为

\[ \Gamma_{\phi,\mathrm{truncation}} = \frac{\dot m\,\Delta x}{2}, \]

作者称之为沿流向（streamwise）数值扩散。这一误差虽降低了精度、改变了有效扩散系数大小，却恰好使解保持有界、物理上正确。

作者指出，要降低这种沿流向的数值扩散必须采用更高阶的对流项近似；但下一节将说明，这种做法必须以保持解的有界性为前提。

11.3.2 下风方案（The Downwind Scheme）

设流动沿 \(+x\)，下风格式给出 \(\phi_e = \phi_E\)、\(\phi_w = \phi_C\)。对 \(\phi_E\) 关于 \(\phi_e\)、\(\phi_C\) 关于 \(\phi_w\) 做 Taylor 展开（均匀网格），截去二阶以上项并代入离散化方程，整理后可得

\[ \dot m_e\phi_E - \dot m_w\phi_C - \Gamma_\phi\frac{\mathrm{d}\phi}{\mathrm{d} x}\bigg|_e \Delta y - \Gamma_\phi\frac{\mathrm{d}\phi}{\mathrm{d} x}\bigg|_w \Delta y = \dot m_e\phi_e - \dot m_w\phi_w - \left(\Gamma_\phi - \frac{\dot m\,\Delta x}{2}\right)\frac{\mathrm{d}\phi}{\mathrm{d} x}\bigg|_e \Delta y - \left(\Gamma_\phi - \frac{\dot m\,\Delta x}{2}\right)\frac{\mathrm{d}\phi}{\mathrm{d} x}\bigg|_w \Delta y. \]

即下风格式引入的截断扩散系数为

\[ \Gamma_{\phi,\mathrm{truncation}} = -\frac{\dot m\,\Delta x}{2}, \]

符号为负，作用是减小实际扩散系数，称为反扩散（anti-diffusion）误差。下风格式的解会"削顶"被对流的廓形，并产生比 CD 方案更剧烈的振荡。

11.3.3 中心差分（CD）方案（The Central Difference Scheme）

CD 方案的截断误差分析较为复杂，因为在有限体积法中梯度计算需要面插值而非直接使用节点值。设速度已知、网格均匀且网格大小为 \(\Delta x\)，则 \(\phi_e - \phi_w\) 的近似式

\[ \tfrac{1}{2}(\phi_E + \phi_C) - \tfrac{1}{2}(\phi_C + \phi_W) = (\phi_e - \phi_w) + \mathrm{TE} \]

为 \(\phi_e - \phi_w\) 的数值近似减去精确值，余项 TE 即为截断误差。

为计算 TE，对 \(\phi_W\)、\(\phi_C\) 关于 \(\phi_w\) 展开，对 \(\phi_C\)、\(\phi_E\) 关于 \(\phi_e\) 展开（网格步长 \(\Delta x\)），取二阶近似后得到面平均值

\[ \tfrac{1}{2}(\phi_E + \phi_C) = \phi_e + \tfrac{\Delta x^2}{8}\phi''_e + \cdots,\qquad \tfrac{1}{2}(\phi_C + \phi_W) = \phi_w + \tfrac{\Delta x^2}{8}\phi''_w + \cdots. \]

两式相减后用 \(\phi''_e - \phi''_w\) 在节点 \(C\) 处的展开，最终在除以 \(\Delta x\) 后得到梯度的截断误差为

\[ \mathrm{TE} = -\frac{\Delta x^2}{8}\phi'''_C + \frac{\Delta x^4}{128}\phi^{\mathrm{v}}_C + \cdots, \]

主项为 \(-\Delta x^2 \phi'''_C/8\)，表明 CD 方案对梯度计算是二阶精度的。

11.4 数值稳定性（Numerical Stability）

作者指出，截断误差 CD（二阶）与迎风（一阶）的精度差异曾让很多研究者推断：既然 CD 用于扩散项很精确，那么用于对流项也应如此。然而前面的分析显示，CD 用于对流项在某些条件下会产生非物理解。根本原因在于：CD 方案在应用于像对流项这样的奇数阶导数时，并不具备内在的对流稳定性（convective stability）。

Leonard 提出对流稳定性概念来解释这一现象。他以一维常速度、含时对流-扩散方程

\[ \frac{\partial(\rho\phi)}{\partial t} + \frac{\partial(\rho u \phi)}{\partial x} - \frac{\partial}{\partial x}\Big(\Gamma_\phi \frac{\partial \phi}{\partial x}\Big) = Q_\phi \]

为例。该方程在以 C 为形心的单元上积分后，左边表示 \(\phi_C\) 在控制体中的时间变化率，右边则表示通过面流入流量的净增量与源项。稳定方案应满足如下要求：若 \(\phi_C\) 出现小的数值偏差，右边的净流入量应对其作负反馈（自校正）。这要求 RHS 关于 \(\phi_C\) 的偏导数为严格负：

\[ \frac{\partial(\mathrm{RHS})}{\partial \phi_C} < 0. \]

作者强调，稳定性 ≠ 有界性 ≠ 精度：稳定方案可能在某些条件下无界（出现超调/欠调或振荡），也可能很扩散、精度低。稳定性的核心是控制数值误差使其不至于无限增长——正如 CD 方案在 \(P_e\) 变化时 \(\phi_C\) 的归一化值从 \(+\infty\) 变到 \(-\infty\) 所表现的那样。

将一般离散化形式代入 CD 方案，作者发现扩散项对 \(\phi_C\) 的灵敏度为

\[ \frac{\partial(\mathrm{RHS}^\mathrm{Diff}_\mathrm{CD})}{\partial \phi_C} = -2\Gamma_\phi \frac{\Delta y}{\Delta x} < 0, \]

是严格负的（\(\Gamma_\phi > 0\)），故 CD 方案的扩散部分是稳定的。对流项的灵敏度为

\[ \frac{\partial(\mathrm{RHS}^\mathrm{Conv}_\mathrm{CD})}{\partial \phi_C} = -\tfrac{1}{2}(\dot m_e + \dot m_w), \]

对稳态流 \(\dot m_e + \dot m_w = 0\)，灵敏度恰好为零。零意味着没有自校正能力；对非稳态流，当流动减速时此项可能为正，成为摆动源，在大 \(P_e\) 下可能引发彻底的数值崩溃。Eq. (11.78) 还表明：对稳态流，用 CD 方案算出的净对流通量与 \(\phi_C\) 的具体取值无关——因而 Fig. 11.9 中无论 \(\phi_C\) 取何值，单元上的净对流通量都相同。

对迎风格式，代入类似分析得到

\[ \frac{\partial(\mathrm{RHS}^\mathrm{Conv}_\mathrm{Upwind})}{\partial \phi_C} = -\max(\dot m_e, 0) - \max(\dot m_w, 0) \le 0, \]

对所有流场非正（一维常截面积下两项不会同时为负）。叠加它本身引入的虚假扩散，作者总结：迎风格式极为稳定，但稳定性的取得以精度为代价。

对下风格式，灵敏度为

\[ \frac{\partial(\mathrm{RHS}^\mathrm{Conv}_\mathrm{Downwind})}{\partial \phi_C} = \max(-\dot m_e, 0) + \max(-\dot m_w, 0) \ge 0, \]

恒非负，叠加其反扩散效应，下风格式是高度不稳定的。

11.5 高阶迎风格式（Higher Order Upwind Schemes）

前面几节表明迎风与中心差分都有严重局限：前者因数值扩散精度差；后者因数值色散（numerical dispersion）误差而失稳。这推动了通过使用高阶迎风偏置插值型分布来同时提高精度与稳定性的研究。本节讨论的高阶格式旨在达到至少二阶精度且无条件稳定。

为强调"通量绑定在单元面上"（而非节点上）的守恒性，作者引入 D（Downwind）、C（Current）、U（Upwind），并在某些情况下再扩展到 DD、UU。Fig. 11.10 给出流动方向为正（→）和负（←）两种情形下这五类节点的相对位置。

11.5.1 二阶迎风方案（Second Order Upwind Scheme）

要构造二阶格式，需采用线性分布（如 CD 方案），但用迎风偏置的节点模板（图 11.11）。线性型 \(\phi(x) = k_0 + k_1(x-x_C)\) 由节点 \(C\)、\(U\) 的 \(\phi\) 值拟合——这意味着面值实际是外推而非插值得到的。

11.5.2 插值分布（The Interpolation Profile）

设 \(\phi(x) = k_0 + k_1(x-x_C)\)，并满足 \(\phi(x_C) = \phi_C\)、\(\phi(x_U) = \phi_U\)，得

\[ \phi(x) = \phi_C + \frac{\phi_C - \phi_U}{x_C - x_U}(x - x_C). \]

把面 \(f\) 的位置 \(x_f\) 代入，并在均匀网格下化简为

\[ \phi_f = \tfrac{3}{2}\phi_C - \tfrac{1}{2}\phi_U. \]

11.5.3 离散化方程（The Discretized Equation）

把该面值代入一维对流-扩散方程 Eq. (11.16)，可得面 \(e\)、\(w\) 上的对流通量

\[ \dot m_e \phi_e = \big(\tfrac{3}{2}\phi_C - \tfrac{1}{2}\phi_W\big)\max(\dot m_e, 0) - \big(\tfrac{3}{2}\phi_E - \tfrac{1}{2}\phi_{EE}\big)\max(-\dot m_e, 0), \]

\[ \dot m_w \phi_w = \big(\tfrac{3}{2}\phi_C - \tfrac{1}{2}\phi_E\big)\max(\dot m_w, 0) - \big(\tfrac{3}{2}\phi_W - \tfrac{1}{2}\phi_{WW}\big)\max(-\dot m_w, 0). \]

代回 Eq. (11.16) 并整理为五对角形式

\[ a_C \phi_C + a_E \phi_E + a_W \phi_W + a_{EE} \phi_{EE} + a_{WW} \phi_{WW} = 0, \]

其中

\[ a_E = \mathrm{Flux}_{F_e} = -\Gamma_\phi \frac{S_e}{d_{xe}} \frac{3}{2} - \max(-\dot m_e, 0)\frac{1}{2} - \max(\dot m_w, 0)\frac{1}{2}, \]

\[ a_{EE} = \mathrm{Flux}_{F_{ee}} = \tfrac{1}{2}\max(-\dot m_e, 0), \]

\[ a_W = \mathrm{Flux}_{F_w} = -\Gamma_\phi \frac{S_w}{d_{xw}}\frac{3}{2} - \max(-\dot m_w, 0)\frac{1}{2} - \max(\dot m_e, 0)\frac{1}{2}, \]

\[ a_{WW} = \mathrm{Flux}_{F_{ww}} = \tfrac{1}{2}\max(-\dot m_w, 0), \]

\[ a_C = -(a_E + a_W + a_{EE} + a_{WW}) + (\dot m_e + \dot m_w). \]

11.5.4 截断误差（Truncation Error）

采用与 CD 方案相同的 Taylor 展开过程，可得二阶迎风（SOU）方案的截断误差为

\[ \mathrm{TE} = -\frac{3}{8}\Delta x^2 \phi'''_C - \frac{1}{4}\Delta x^3 \phi^{\mathrm{iv}}_C + \cdots, \]

主项为 \(-3\Delta x^2 \phi'''_C/8\)，是二阶精度。

11.5.5 稳定性分析（Stability Analysis）

代入 SOU 格式的对流通量到 RHS 中，对 \(\phi_C\) 求偏导：

\[ \frac{\partial(\mathrm{RHS}^\mathrm{Conv})}{\partial \phi_C} = -\tfrac{3}{2}\max(\dot m_e, 0) - \tfrac{3}{2}\max(\dot m_w, 0). \]

对所有流场该值恒为非正，因此 SOU 是稳定格式。但作者强调，该结论在推导条件下（速度恒定）成立；当速度不恒定时一般不成立。

11.5.6 QUICK 方案（The QUICK Scheme）

QUICK（Quadratic Upstream Interpolation for Convective Kinematics）由 Leonard 提出。它用二次多项式迎风偏置插值得到单元面上的 \(\phi\) 值（图 11.12），再用来计算对流项。原始形式需要一般多维二次多项式；但在实践中通常把流动视为局部一维，在每个坐标方向上用一维二次型即可。

11.5.7 QUICK 方案的插值分布（The Interpolation Profile）

一维情形下设 \(\phi = k_0 + k_1 x + k_2 x^2\)，由三个节点 \(U\)、\(C\)、\(D\) 的 \(\phi\) 值定出三个系数，得到

\[ \phi = \phi_U + \frac{(x-x_U)(x-x_C)}{(x_D-x_U)(x_D-x_C)}(\phi_D - \phi_U) + \frac{(x-x_U)(x-x_D)}{(x_C-x_U)(x_C-x_D)}(\phi_C - \phi_U). \]

在均匀网格上，面 \(f\) 上的 \(\phi\) 简化为

\[ \phi_f = \frac{\phi_C + \phi_D}{2} - \frac{\phi_D - 2\phi_C + \phi_U}{8}. \]

把 \(\phi_f\) 代入对流通量公式，并整理为五对角形式 \(a_C \phi_C + a_E \phi_E + a_W \phi_W + a_{EE} \phi_{EE} + a_{WW} \phi_{WW} = 0\) ，其系数为

\[ a_E = -\Gamma_\phi \frac{S_e}{d_{xe}}\tfrac{3}{4} - \tfrac{3}{8}\max(-\dot m_e, 0) + \tfrac{1}{8}\max(\dot m_e, 0) - \tfrac{1}{8}\max(\dot m_w, 0), \]

\[ a_W = -\Gamma_\phi \frac{S_w}{d_{xw}}\tfrac{3}{4} - \tfrac{3}{8}\max(-\dot m_w, 0) + \tfrac{1}{8}\max(\dot m_w, 0) - \tfrac{1}{8}\max(\dot m_e, 0), \]

\[ a_{EE} = \tfrac{1}{8}\max(-\dot m_e, 0),\quad a_{WW} = \tfrac{1}{8}\max(-\dot m_w, 0), \]

\[ a_C = \Gamma_\phi \frac{S_e}{d_{xe}} + \Gamma_\phi \frac{S_w}{d_{xw}} + \max(\dot m_e, 0) - \tfrac{3}{8}\max(-\dot m_e, 0) + \max(\dot m_w, 0) - \tfrac{3}{8}\max(-\dot m_w, 0) = -\sum_F a_F + (\dot m_e+\dot m_w). \]

11.5.8 QUICK 方案的截断误差（Truncation Error）

经同样 Taylor 展开过程，QUICK 方案的截断误差为

\[ \mathrm{TE} = \tfrac{1}{16}\Delta x^3 \phi^{\mathrm{iv}}_C - \tfrac{3}{128}\Delta x^4 \phi^{\mathrm{v}}_C + \cdots, \]

主项为 \(\Delta x^3 \phi^{\mathrm{iv}}_C/16\)，是三阶精度的。

11.5.9 QUICK 方案的稳定性分析（Stability Analysis）

代入 QUICK 对流通量后，对 \(\phi_C\) 求偏导得

\[ \frac{\partial(\mathrm{RHS}^\mathrm{Conv})}{\partial \phi_C} = -\tfrac{3}{8}\max(\dot m_e, 0) - \tfrac{3}{8}\max(\dot m_w, 0) - \tfrac{3}{8}(\dot m_e + \dot m_w). \]

对均匀速度场该值恒为负，QUICK 是稳定格式。但作者强调，这并不保证一般非均匀速度场下的有界性。

11.5.10 FROMM 方案（The FROMM Scheme）

FROMM 方案在跨面的远上游 \(U\) 与下游 \(D\) 节点之间拟合一条线性分布（图 11.13），其模板相对于单元是非对称的迎风偏置——假定 \(U\)、\(C\)、\(D\) 共线。

11.5.11 FROMM 方案的插值分布（The Interpolation Profile）

设 \(\phi(x) = k_0 + k_1(x-x_C)\)，并由 \(\phi(x_D) = \phi_D\)、\(\phi(x_U) = \phi_U\) 拟合：

\[ \phi(x) = \phi_U + \frac{\phi_D - \phi_U}{x_D - x_U}(x - x_U). \]

在 \(x = x_C\) 处得 \(\phi_C = \phi_U + \tfrac{\phi_D-\phi_U}{x_D-x_U}(x_C - x_U)\) ，均匀网格下化简为 \(\phi_C = (\phi_D + \phi_U)/2\)。把 \(x_f\) 代入后再利用 \(x_C\) 处的表达，最终面 \(f\) 值为

\[ \phi_f = \phi_C + \frac{x_f - x_C}{x_D - x_U}(\phi_D - \phi_U), \]

均匀网格下简化为

\[ \phi_f = \phi_C + \frac{\phi_D - \phi_U}{4}. \]

11.5.12 FROMM 方案的离散化方程（The Discretized Equation）

把 FROMM 面值代入对流通量，可得

\[ \dot m_e \phi_e = \big(\phi_C - \tfrac{1}{4}\phi_W + \tfrac{1}{4}\phi_E\big)\max(\dot m_e, 0) - \big(\phi_E - \tfrac{1}{4}\phi_{EE} + \tfrac{1}{4}\phi_C\big)\max(-\dot m_e, 0), \]

\[ \dot m_w \phi_w = \big(\phi_C - \tfrac{1}{4}\phi_E + \tfrac{1}{4}\phi_W\big)\max(\dot m_w, 0) - \big(\phi_W - \tfrac{1}{4}\phi_{WW} + \tfrac{1}{4}\phi_C\big)\max(-\dot m_w, 0). \]

整理为 \(a_C \phi_C + a_E \phi_E + a_W \phi_W + a_{EE}\phi_{EE} + a_{WW}\phi_{WW} = 0\) ，系数为

\[ a_E = -\Gamma_\phi \frac{S_e}{d_{xe}}\tfrac{1}{4} + \tfrac{1}{4}\max(\dot m_e, 0) - \tfrac{1}{4}\max(-\dot m_e, 0) - \tfrac{1}{4}\max(\dot m_w, 0), \]

\[ a_{EE} = \tfrac{1}{4}\max(-\dot m_e, 0), \]

\[ a_W = -\Gamma_\phi \frac{S_w}{d_{xw}}\tfrac{1}{4} + \tfrac{1}{4}\max(\dot m_w, 0) - \tfrac{1}{4}\max(-\dot m_w, 0) - \tfrac{1}{4}\max(\dot m_e, 0), \]

\[ a_{WW} = \tfrac{1}{4}\max(-\dot m_w, 0), \]

\[ a_C = -(a_E + a_W + a_{EE} + a_{WW}) + (\dot m_e + \dot m_w). \]

11.5.13 FROMM 方案的截断误差（Truncation Error）

用相同方法推导，FROMM 的截断误差为 \(O(\Delta x^2)\)，即二阶精度。

11.5.14 FROMM 方案的稳定性分析（Stability Analysis）

把 FROMM 对流通量代入 RHS 后求导：

\[ \frac{\partial(\mathrm{RHS}^\mathrm{Conv})}{\partial \phi_C} = -\tfrac{3}{4}\max(\dot m_e, 0) - \tfrac{3}{4}\max(\dot m_w, 0) - \tfrac{1}{4}(\dot m_e + \dot m_w). \]

对常速度场恒为负，FROMM 是稳定格式；但速度变化时该结论一般不成立。

11.5.15 各方案比较（Comparison of the Various Schemes）

作者把上述各方案关于 \(\phi_C\) 的灵敏度并列比较：SOU 方案最负（系数 \(-3/2\)）；迎风其次（\(-1\)）；FROMM 居中（\(-3/4\)）；QUICK 较小（\(-3/8\)）；CD 方案为零（即对 \(\phi_C\) 变化"中性"）。真正的"自校正"能力是对流与扩散贡献的总和；由于迎风格式引入了额外的虚假扩散，其总体稳定度最高。

Fig. 11.14 在 \(L=1\)、\(\phi(0)=1\)、\(\phi(1)=0\) 的 Dirichlet 条件下比较了各格式的解：低 \(P_e\)（如 \(P_e=1\)，图 a）时所有解都稳定，CD 与 QUICK 精度相当且接近精确解，迎风最差；FROMM 精度高于 SOU，SOU 高于迎风，QUICK 比 FROMM 更精确。高 \(P_e\)（图 b，\(P_e=10\)）时，CD、FROMM、QUICK 出现明显摆动；只有迎风与 SOU 保持平滑且几乎一致——这表明 SOU 仍是高度扩散的。作者指出这些摆动是出流边界强加值与对流主导场上游依赖性冲突所引起的过/欠调。迎风与 SOU 仅基于上游值，因而对出口强加值不敏感；但 SOU 在大梯度区（如激波）下预期会出现振荡。

11.5.16 均匀与非均匀网格上的函数关系（Functional Relationships for Uniform and Non-uniform Grids）

以上各种插值型分布都是在一维笛卡尔网格上推导的。沿曲线坐标轴 \(\xi\)，只要把 \(x\) 替换为 \(\xi\) 即可使用同一函数关系（图 11.15）：均匀网格下函数关系与笛卡尔情形完全一致；非均匀网格下自变量应替换为沿坐标轴的距离 \(\xi\)。若 \(O\) 为原点，\(\xi_U\) 可由下式计算

\[ \xi_U = \xi_{UU} + (\xi_U - \xi_{UU}),\quad \xi_{UU} = \sqrt{(x_{UU}-x_O)^2 + (y_{UU}-y_O)^2 + (z_{UU}-z_O)^2}, \]

\[ \xi_U - \xi_{UU} = \sqrt{(x_U - x_{UU})^2 + (y_U - y_{UU})^2 + (z_U - z_{UU})^2}. \]

更一般地，沿曲线距离由折线累加得到

\[ \xi_1 - \xi_2 = \sqrt{(x_1-x_2)^2 + (y_1-y_2)^2 + (z_1-z_2)^2}. \]

作者把上述各方案在均匀与非均匀网格上的函数关系整理为表（Table 11.1），其中上栏（均匀）一列与前文一致，下栏（非均匀）则把 \(x\) 替换为 \(\xi\)。例如，QUICK 在非均匀网格上的面值为

\[ \phi_f = \phi_U + \frac{(\xi_f - \xi_U)(\xi_f - \xi_C)}{(\xi_D - \xi_U)(\xi_D - \xi_C)}(\phi_D - \phi_U) + \frac{(\xi_f - \xi_U)(\xi_f - \xi_D)}{(\xi_C - \xi_U)(\xi_C - \xi_D)}(\phi_C - \phi_U). \]

11.6 稳态二维对流（Steady Two Dimensional Advection）

稳态二维对流方程为 \(\nabla\cdot(\rho\mathbf{v}\phi) = 0\)。在图 11.16 所示的二维单元 \(V_C\) 上做体积分、用散度定理转化为面积分、再用面求和代替，得

\[ \sum_{f\sim\mathrm{nb}(C)} \big(\rho\mathbf{v}\phi\big)_f\cdot\mathbf{S}_f = 0. \]

对面 \(f\) 上的积分采用单一高斯点近似，式子变为

\[ \sum_f (\rho\mathbf{v}\phi)_f\cdot\mathbf{S}_f = 0 \]

。在笛卡尔网格上把单元面分到 e、w、n、s 四个方向，全离散形式为

\[ (\rho u \Delta y\, \phi)_e - (\rho u \Delta y\, \phi)_w + (\rho v \Delta x\, \phi)_n - (\rho v \Delta x\, \phi)_s = 0. \]

若在每一坐标方向上把流动视为局部一维而采用迎风格式，则代数方程变为 \(a_C\phi_C + a_E\phi_E + a_W\phi_W + a_N\phi_N + a_S\phi_S = 0\) ，其中

\[ a_E = \mathrm{Flux}_{F_e} = -\max(-\dot m_e, 0),\quad a_W = -\max(-\dot m_w, 0), \]

\[ a_N = \mathrm{Flux}_{F_n} = -\max(-\dot m_n, 0),\quad a_S = -\max(-\dot m_s, 0), \]

\[ a_C = \max(\dot m_e, 0) + \max(\dot m_w, 0) + \max(\dot m_n, 0) + \max(\dot m_s, 0) = -\sum_F a_F + \sum_f \dot m_f. \]

若采用 QUICK 方案，离散化方程改为

\[ a_C \phi_C + a_E \phi_E + a_W \phi_W + a_{EE}\phi_{EE} + a_{WW}\phi_{WW} + a_N \phi_N + a_S \phi_S + a_{NN}\phi_{NN} + a_{SS}\phi_{SS} = 0, \]

各系数为

\[ a_E = -\tfrac{3}{4}\max(-\dot m_e, 0) + \tfrac{3}{8}\max(\dot m_e, 0) - \tfrac{1}{8}\max(\dot m_w, 0), \]

\[ a_W = -\tfrac{3}{4}\max(-\dot m_w, 0) + \tfrac{3}{8}\max(\dot m_w, 0) - \tfrac{1}{8}\max(\dot m_e, 0), \]

\[ a_{EE} = \tfrac{1}{8}\max(-\dot m_e, 0),\quad a_{WW} = \tfrac{1}{8}\max(-\dot m_w, 0), \]

\[ a_N = -\tfrac{3}{4}\max(-\dot m_n, 0) + \tfrac{3}{8}\max(\dot m_n, 0) - \tfrac{1}{8}\max(\dot m_s, 0), \]

\[ a_S = -\tfrac{3}{4}\max(-\dot m_s, 0) + \tfrac{3}{8}\max(\dot m_s, 0) - \tfrac{1}{8}\max(\dot m_n, 0), \]

\[ a_{NN} = \tfrac{1}{8}\max(-\dot m_n, 0),\quad a_{SS} = \tfrac{1}{8}\max(-\dot m_s, 0), \]

\[ a_C = -\sum_F a_F + (\dot m_e + \dot m_w + \dot m_n + \dot m_s). \]

作者随后把两种格式应用到一个"斜向速度场中纯对流阶跃廓形"问题：单位方域，速度 \(\mathbf{v} = (1,1)\)，左侧 \(\phi = 1\)、底部 \(\phi = 0\)（图 11.17a）。无扩散的精确解为 \(\phi = 1\)（在图示对角线以上）与 \(\phi = 0\)（以下）。在 \(x=0.5\) 处的廓形对比（图 11.17b）表明：迎风格式输出模糊且精度差但很平滑；QUICK 输出更锐利、精度更高，但在大梯度附近出现超/欠调。

迎风方案精度差的根源是一种新的误差——横向（cross-stream）数值扩散，其本质是把流动当作局部一维处理时引入的"一维插值型"造成。Patankar 和 Stubley 将其识别为多维现象，仅在速度场不与网格对齐时出现。de Vahl Davis 和 Mallinson 给出的二维近似表达式为

\[ \Gamma_{\phi,\mathrm{false}} = \frac{\rho |\mathbf{v}| \Delta x \Delta y \sin(2\theta)}{4\big[\Delta y \sin^3\theta + \Delta x \cos^3\theta\big]}, \]

其中 \(\theta\) 是速度矢量与 \(x\) 轴的夹角。提高插值格式的阶数（如改用 QUICK）可以显著减少这一误差；但 QUICK 仍会因大梯度而出现超/欠调——这正是高阶格式的特征性色散误差。

11.6.1 误差来源（Error Sources）

综合前面的讨论，作者把对流通量离散化中的数值误差分为数值扩散与数值色散两大类。

数值扩散造成锐利梯度被抹平（图 11.18a），又可细分为沿流向（stream wise）与横向（cross stream）两类。沿流向扩散可通过提高插值阶数（图 11.18b）来减小——更锐利的廓形随之而来，但在大梯度区会引入超/欠调。横向扩散（图 11.17b）由"一维插值型"的本质造成，可通过沿流动方向插值（多维型）或采用一维高阶插值型（图 11.18b）来减小。

数值色散误差在大梯度附近表现为振荡，使解变得无界。作者指出，除迎风格式外，所有插值型分布都存在色散误差；其本质是所选插值型的非物理行为（图 11.18b、图 11.19、图 11.14b）。

作者随后给出色散误差的解析度量方法。从忽略扩散与源项、常速度常密度的简化方程

\[ \frac{\partial \phi}{\partial t} + u \frac{\partial \phi}{\partial x} = 0 \]

出发，假定形如 \(\phi(x,t) = \phi(t)\mathrm{e}^{\mathrm{i}kx}\) 的精确解（其中 \(\mathrm{i}^2 = -1\)），其导数为

\[ \partial\phi/\partial x = \mathrm{i}k\phi(x,t) \]

。若数值格式把梯度近似为

\[ \frac{\partial \phi}{\partial x} \approx \frac{1}{\Delta x}\sum_{n=-M}^{N} a_n \phi(x+n\Delta x, t) = \frac{1}{\Delta x}\sum_{n=-M}^{N} a_n \mathrm{e}^{\mathrm{i}kn\Delta x}\phi(x,t), \]

则数值波数 \(k\) 满足

\[ k = \frac{\mathrm{i}}{\Delta x}\sum_{n=-M}^{N} a_n \mathrm{e}^{\mathrm{i}kn\Delta x}. \]

一般 \(k\) 是复数，\(k = \mathrm{Re}(k) + \mathrm{i}\mathrm{Im}(k)\)，于是近似解为

\[ \phi(x,t) = \phi(t)\mathrm{e}^{\mathrm{i}kx} = \phi(t)\underbrace{\mathrm{e}^{\mathrm{i}\mathrm{Re}(k)x}}_{\text{相位/色散}}\underbrace{\mathrm{e}^{-\mathrm{Im}(k)x}}_{\text{幅度/耗散}}. \]

若 \(k\) 是实数，则只有色散误差；若 \(k\) 是复数，则色散与耗散误差同时存在。

对迎风格式，梯度近似为 \((\phi_C - \phi_W)/\Delta x\)，代入精确解后得

\[ k = \frac{\sin(k\Delta x)}{\Delta x} - \mathrm{i}\,\frac{1-\cos(k\Delta x)}{\Delta x}, \]

即 \(k\) 为复数，迎风格式同时具有耗散与色散误差。对CD 方案，梯度为 \((\phi_E - \phi_W)/(2\Delta x)\)，代入得

\[ k = \frac{\sin(k\Delta x)}{\Delta x}, \]

\(k\) 纯虚（按原文："纯实"对应实数轴上的 \(\mathrm{i}\sin\)），故只有色散误差——这正是 CD 方案在锐利梯度处出现振荡与超/欠调的根因。

最后，作者总结：理解了色散误差的本质之后，下一步自然就是发展对流非振荡且高阶的格式；这一目标直到"对流通量有界化"（bounding of the convective flux）方法被理解后才得以实现，相关发展将在下一章中详述。

11.7 非结构网格上的高阶方案（High Order Schemes on Unstructured Grids）

与结构网格类似，HO 方案的函数关系也由 \(U\)、\(C\)、\(D\) 三个节点的 \(\phi\) 值定义。在非结构网格中，\(C\)、\(D\) 节点对任何内部面都容易确定（图 11.20a），但 \(U\) 节点的定义并不直接（图 11.20b、c）。最直接的解决办法是把 HO 方案改写为基于 \(C\)、\(D\) 节点梯度的形式；另一种方法是重建一个伪 \(U\) 节点，这将在下一章与 HR 方案一起详述。

11.7.1 以梯度形式重写 HO 方案（Reformulating HO Schemes in Terms of Gradients）

这一方法建立在结构网格上发展出的"型"之上，最清晰的解释方式是用 QUICK 方案为例，把它的"型"用新术语重写，然后推广到所有用三个点构造的二阶型。

QUICK 方案的函数关系可以写为

\[ \phi_f = \phi_C + \frac{1}{4}\Big(\frac{\phi_D - \phi_U}{2}\Big) + \frac{1}{4}(\phi_D - \phi_C). \]

若用 \(d_{Cf}\) 方向（图 11.19）上 \(C\) 与面 \(f\) 处的梯度近似

\[ \frac{\partial \phi_C}{\partial \xi} = \frac{\phi_D - \phi_U}{2\Delta \xi_f},\qquad \frac{\partial \phi_f}{\partial \xi} = \frac{\phi_D - \phi_C}{\Delta \xi_f}, \]

则 QUICK 方案可改写为

\[ \phi_f = \phi_C + \frac{1}{2}\frac{\partial \phi_C}{\partial \xi}\frac{\Delta \xi_f}{2} + \frac{1}{2}\frac{\partial \phi_f}{\partial \xi}\frac{\Delta \xi_f}{2}, \]

或以向量形式（\(d_{Cf}\) 为 \(C\) 与 \(f\) 之间的向量）写为

\[ \phi_f = \phi_C + \frac{1}{2}\nabla\phi_C\cdot\mathbf{d}_{Cf} + \frac{1}{2}\nabla\phi_f\cdot\mathbf{d}_{Cf}. \]

这非常适合在非结构网格中应用，因为只需要 \(C\) 与 \(f\) 处的梯度信息；只要这些梯度的计算是二阶精度的，具体计算方式无关紧要。这一表达进一步暗示：任何基于三个点的型都可以写为

\[ \phi_f = a\phi_C + b\,\nabla\phi_C\cdot\mathbf{d}_{Cf} + c\,\nabla\phi_f\cdot\mathbf{d}_{Cf}, \]

其中常数 \(a\)、\(b\)、\(c\) 由"型"在结构网格上的形式确定。把 \(\nabla\phi_C\)、\(\nabla\phi_f\) 的近似代入后可整理为

\[ \phi_f = \Big(a - \frac{c}{2}\Big)\phi_C + \Big(\frac{b}{4} + \frac{c}{4}\Big)\phi_D - \frac{b}{4}\phi_U, \]

其中

\[ -\frac{c}{2} = \frac{3}{2} \]

、

\[ \frac{b}{4} + \frac{c}{4} = 0 \]

、

\[ a - \frac{c}{2} = \frac{3}{2} \]

解出 \(a=1\)、\(b=1/2\)、\(c=-1\)。于是 SOU 方案的等价梯度形式为

\[ \phi_f = \phi_C + (2\nabla\phi_C - \nabla\phi_f)\cdot\mathbf{d}_{Cf}. \]

按同样方法，所有前述方案的梯度形式分别为

迎风：\(\phi_f = \phi_C\)
CD：

\[ \phi_f = \phi_C + \nabla\phi_f\cdot\mathbf{d}_{Cf} \]

SOU：

\[ \phi_f = \phi_C + (2\nabla\phi_C - \nabla\phi_f)\cdot\mathbf{d}_{Cf} \]

FROMM：

\[ \phi_f = \phi_C + \nabla\phi_C\cdot\mathbf{d}_{Cf} \]

QUICK：

\[ \phi_f = \phi_C + \tfrac{1}{2}(\nabla\phi_C + \nabla\phi_f)\cdot\mathbf{d}_{Cf} \]

下风：

\[ \phi_f = \phi_C + 2\nabla\phi_f\cdot\mathbf{d}_{Cf} \]

作者提示，关于 \(\nabla\phi_C\)、\(\nabla\phi_f\) 的具体计算方法已在第 9 章详述。

11.8 延迟修正方法（The Deferred Correction Approach）

延迟修正（Deferred Correction, DC）方法由 Khosla 和 Rubin 提出。它是一种压缩（compacting）技术，使原本只写低阶方案的代码能够使用 HO 方案而不违反任何稳定性规则，且适用于任何结构或非结构网格。方法的核心是把面 \(f\) 上用 HO 方案计算的对流通量写为

\[ \dot m_f \phi_f^{\mathrm{HO}} = \underbrace{\dot m_f \phi_f^{\mathrm{Upwind}}}_{\text{隐式}} + \underbrace{\dot m_f \big(\phi_f^{\mathrm{HO}} - \phi_f^{\mathrm{Upwind}}\big)}_{\text{显式}}. \]

第一项按迎风格式隐式表达为节点值的函数，第二项（HO 与迎风之差）用最近一次迭代的 \(\phi\) 值显式计算。以节点值表示时，

\[ \dot m_f \phi_f^{\mathrm{HO}} = \max(\dot m_f, 0)\phi_C - \max(-\dot m_f, 0)\phi_F + \dot m_f\big(\phi_f^{\mathrm{HO}} - \max(\dot m_f, 0)\phi_C + \max(-\dot m_f, 0)\phi_F\big) = \mathrm{Flux}_{C_f}\phi_C + \mathrm{Flux}_{F_f}\phi_F + \mathrm{Flux}_{V_f}, \]

其中 \(\mathrm{Flux}_{C_f} = \max(\dot m_f, 0)\)、\(\mathrm{Flux}_{F_f} = -\max(-\dot m_f, 0)\)、

\[ \mathrm{Flux}_{V_f} = \dot m_f \phi_f^{\mathrm{HO}} - \mathrm{Flux}_{C_f}\phi_C - \mathrm{Flux}_{F_f}\phi_F \]

。代回离散化方程并整理得

\[ a_C \phi_C + \sum_{F\sim\mathrm{NB}(C)} a_F \phi_F = b_C, \]

其中

\[ a_F = \mathrm{Flux}_{F_f} = -\max(-\dot m_f, 0),\quad a_C = \sum_f \mathrm{Flux}_{C_f} = \sum_f \big[\max(\dot m_f, 0) + \max(-\dot m_f, 0)\big] = -\sum_F a_F + \sum_f \dot m_f, \]

而

\[ b_C = Q_\phi V_C - \sum_{f\sim\mathrm{nb}(C)} \mathrm{Flux}_{V_f} = Q_\phi V_C - \sum_{f\sim\mathrm{nb}(C)} \dot m_f\big(\phi_f^{\mathrm{HO}} - \phi_f^{\mathrm{Upwind}}\big). \]

\(b_C^{\mathrm{DC}}\) 体现了 DC 过程引入的额外源项；DC 技术的关键性质是系数矩阵总是对角占优的，因为它是用迎风格式构造的。

Example 2 演示了对 QUICK 方案进行 DC 实施时的系数推导：迎风系数作为 \(a_F\)、\(a_C\)，差值 \(\phi_f^{\mathrm{QUICK}} - \phi_f^{\mathrm{UPWIND}}\) 进入源项

\[ b_C = S_\phi V_C - \sum_{f\sim\mathrm{nb}(C)} \dot m_f\Big(\tfrac{3}{4}\phi_C - \tfrac{9}{8}\phi_U + \tfrac{3}{8}\phi_D\Big). \]

这一压缩过程实现简单，但当 HO 与迎风之间的面值差变大时收敛速度会下降。

11.9.1 uFVM

在 uFVM 中，对流项的装配由 cfdAssembleConvectionTerm 完成：先为内部面装配（cfdAssembleConvectionTermInterior），再对各边界 patch 循环装配（cfdAssembleConvectionTermInletBC、cfdAssembleConvectionTermOutletBC 等）。

内部面装配核心（Listing 11.1）取出内部面的 owner、neighbour 索引、定义迎风索引 pos（pos = 1 当 \(\dot m_f > 0\)），并按迎风格式计算系数：

theMdotName = ['Mdot' theFluidTag];
mdotField  = cfdGetMeshField(theMdotName,'Faces');
mdot_f     = mdotField.phi(iFaces);
iOwners    = [theMesh.faces(iFaces).iOwner];
iNeighbours= [theMesh.faces(iFaces).iNeighbour];
pos        = zeros(size(mdot_f));
pos((mdot_f>0)) = 1;
theFluxes.FLUXC1f(iFaces,1) = mdot_f.*pos;
theFluxes.FLUXC2f(iFaces,1) = mdot_f.*(1-pos);
theFluxes.FLUXVf(iFaces,1)  = 0;

HO 方案在迎风装配之后通过 DC 方法实现。QUICK 装配（cfdAssembleConvectionTermDCQUICK，Listing 11.2）核心为：取迎风梯度 phiGradCf、把梯度插值到面 phiGradf、计算 \(\mathbf{d}_{Cf} = \mathbf{r}_f - \mathbf{r}_C\)，并把修正量

\[ \mathrm{corr} = \dot m_f \cdot \big(\nabla\phi_C + \nabla\phi_f\big)\cdot\mathbf{d}_{Cf} \cdot \tfrac{1}{2} \]

加到 FLUXTf 上。该 DC 值对所有内部面为

\[ \dot m_f \big(\phi_f^{\mathrm{HR}} - \phi_f^{\mathrm{UPWIND}}\big) = \dot m_f\Big[\phi_C + \tfrac{1}{2}\big(\nabla\phi_C + \nabla\phi_f\big)\cdot\mathbf{d}_{Cf} - \phi_C\Big] = \tfrac{1}{2}\dot m_f\big(\nabla\phi_C + \nabla\phi_f\big)\cdot\mathbf{d}_{Cf}. \]

11.9.2 OpenFOAM®

在 OpenFOAM 中，对流项可显式用 fvc::div(mDot, phi)，或隐式用 fvm::div(mDot, phi)。fvc::div(...) 返回一个在每个单元上计算 \(\phi\) 散度的场，加到方程右边；fvm::div(...) 返回一个 fvMatrix，其系数由面值通量的线性化得到，加到方程左边。

两函数的脚本位于 FOAM_SRC/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.C。对流类基于"在面上使用的插值格式类型"声明，Listing 11.3 给出 gaussConvectionScheme 类的声明，其私有成员 tinterpScheme_ 描述面插值类型。该类还定义两个辅助函数 interpolate、flux（Listing 11.4），把 surfaceInterpolation 类包装为从体场到面场的插值。

显式离散化（fvc::div，Listing 11.5）通过 fvc::surfaceIntegrate(flux(faceFlux, vf)) 把 faceFlux 与插值场相乘后做面积分。surfaceIntegrate 函数（Listing 11.6）遍历所有面，对 owner 单元做加、对 neighbour 单元做减，最后除以单元体积 \(V\)。对显式对流项，divergence 算子由以下三步组成：(1) 在单元面上求 \(\phi_f\)；(2) 将面值乘以面质量通量 faceFlux；(3) 把所有面的贡献求和并除以单元体积。

隐式离散化（fvm::div，Eq. 11.163 / Listing 11.7）把面值隐式表达为 owner 与 neighbour 节点值的加权平均

\[ \phi_f = \gamma\phi_O + (1-\gamma)\phi_N, \]

其中下标 O、N 分别指 owner 与 neighbour，\(\gamma\) 为 owner 的权重。隐式算子通过

fvm.lower() = -weights.internalField()*faceFlux.internalField();
fvm.upper() =  fvm.lower() + faceFlux.internalField();
fvm.negSumDiag();

填充系数矩阵。作者指出，OpenFOAM 隐式实现对应"下风加权"——所有格式都以全隐式方式离散，与阶数无关、无需 DC；这与下一章要讲的"下风加权因子法"思路一致。

如前所述，fvc::div 与 fvm::div 都通过 surfaceInterpolationScheme 基类完成面值计算与权重计算。该基类派生出大量插值格式（UML 图 Fig. 11.21 列出 clippedLinear、CoBlended、fixedBlended、downwind、harmonic、limitedSurfaceInterpolationScheme、limiterBlended、limitWith、linear、localBlended、localMax、localMin、midPoint、outletStabilised、reverseLinear、skewCorrected、surfaceInterpolationScheme、weighted 等），其 UML 层次由 Fig. 11.21 给出。

基类中两个主要函数在 surfaceInterpolationScheme.H 中定义：interpolate（Listing 11.11）显式计算面插值；weights（Listing 11.12）是纯虚函数，由每个派生类必须实现。

interpolate 在 surfaceInterpolationScheme 中以 Eq. (11.164) 的改写形式

\[ \phi_f = \gamma\phi_O + (1-\gamma)\phi_N = \phi_N + \gamma(\phi_O - \phi_N) \]

实现——Listing 11.14 给出双参数版本，遍历所有面计算 sfi[fi] = lambda[fi]*(vfi[P[fi]] - vfi[N[fi]]) + vfi[N[fi]]。由于 weights 是纯虚，每个派生类都必须重写它。CD 方案对应 linear<Type>，其 weights 函数（Listing 11.15）直接返回 mesh().surfaceInterpolation::weights()，即基于网格的距离权重。下风格式（Eq. 11.44）的 weights（Listing 11.16）返回 neg(faceFlux_)，其中 neg 在正通量时返回 0、负通量时返回 1——正好与迎风格式所用 pos 函数相反。

11.10 小结（Closure）

本章讨论了对流项的离散化。线性对称型分布的困难被分析清楚，并提出了若干修正方法。迎风格式虽然稳定、给出物理上合理的结果，但高度扩散——它抹平锐利梯度且仅有一阶精度。迎风偏置的 HO 方案提高了精度，但带来一种新的误差——色散误差——表现为大梯度附近的振荡、超调与欠调。本章未试图处理这一误差；它将在下一章（继续讨论对流项的离散化）中作为核心议题被进一步展开。

本章个人批注

本章是 Moukalled《FVM in CFD》中"对流项离散化"的核心章节，篇幅巨大、细节繁多。最值得记的逻辑链条是：CD 失败 → 迎风补救 → 迎风精度差（沿流向虚假扩散）→ HO 方案（SOU/QUICK/FROMM）提高精度 → HO 引入色散误差 → 下一章 HR 格式处理色散。这条主线是整个两章对流离散化叙事的骨架。

稳定性分析（11.4）是我个人读到的最有洞见的部分。Leonard 的"对流稳定性"概念——RHS 关于 \(\phi_C\) 的偏导数为负则误差自校正——是一个把"稳定性"从纯代数（Diagonally dominant 之类）拉回到物理层面的视角。它把迎风为什么稳、CD 为什么在大 \(P_e\) 下崩、为什么下风稳/不稳三件事统一到同一判据下。这个判据在第 12 章谈 TVD/HR 格式时还会再次出现（极限器本质上是约束数值通量关于 \(\phi_C\) 的灵敏度），所以现在把 Eq. (11.74)–(11.82) 这一组等式记牢会有复利。

Taylor 截断误差的物理读法也很有意思：迎风的沿流向扩散 \(\Gamma_{\phi,\mathrm{trunc}} = \dot m \Delta x/2\)（正）解释了它"抹平"的特性；下风的 \(-\dot m \Delta x/2\)（负）解释了它"削顶"且更振荡；CD 的 \(-\Delta x^2 \phi'''_C/8\) 解释了它在锐利梯度（二阶导不连续或高阶导大）处的振荡。这三类误差在 Fourier 分析（11.6.1）里被统一为"色散 vs 耗散"两部分——\(\mathrm{Re}(k)\) 控相位畸变、\(\mathrm{Im}(k)\) 控幅度衰减。迎风 \(k\) 为复（双误差都有），CD \(k\) 纯虚（只有相位畸变 = 振荡）。这个 11.6.1 节是连接"工程格式"与"波动传播理论"的桥梁。

QUICK 截断误差为三阶（\(\Delta x^3 \phi^{\mathrm{iv}}_C/16\)）这一事实意味着它在光滑解上精度最高，但也意味着它对四阶导数（曲率）非常敏感——尖锐界面、激波、接触间断等处 \(\phi^{\mathrm{iv}}\) 大，QUICK 必振荡。这为下一章 NVD/TVD 类"限制器"提供了动机：纯 QUICK 在工程中常用作基准，但实际应用时几乎总需要某种 limiter。

非结构网格上的梯度形式重写（11.7.1）是另一个工程价值很高的内容。Eqs. (11.150)–(11.155) 把所有一维推导出的方案压缩为 \(\phi_f = \phi_C + (\text{梯度组合})\cdot\mathbf{d}_{Cf}\) 的形式，这正是 OpenFOAM 和 uFVM 中高阶方案在非结构网格上实现时所采用的统一形式。理解这个等式族后，回看 OpenFOAM surfaceInterpolationScheme 的 interpolate 与 weights 函数（Listing 11.9–11.16），就明白它们本质上是这一组梯度形式的代码实现。

实际工程上的判断：当对流主导、锐利梯度存在时，单纯高阶格式（QUICK、FROMM）会摆动，但迎风又太平。在没有 limiter 的情况下，工程上常用 SOU 配加密网格的折中——它至少有界、二阶精度、且实现简单。HR 方案（TVD 家族）则是"既不抹平、也不振荡"的最终方案，留给下一章。

与上下章的衔接（一段话）

第 8 章建立了扩散项的离散化方法，第 9 章系统讲解了梯度的各种计算方式，第 10 章求解代数方程组——这三章为对流项的离散化准备了所有"工具箱"零件：扩散项离散化范式、梯度插值方法、稀疏矩阵求解器。第 11 章正是把这些零件组装起来处理对流项这个"看似简单实则棘手"的核心算子。本章的结论是对流项离散化面临精度与稳定性的两难——CD 失稳、迎风扩散、HO 引入色散——这一两难自然把读者带到第 12 章。下一章将引入对流通量有界化（bounding of the convective flux）以及高分辨率（HR）格式族（TVD/NVD），既保留 HO 格式在光滑区的精度、又通过 limiter 把大梯度处的超/欠调约束掉，是"既不抹平、也不振荡"的最终解法。第 13 章则把瞬态项的离散化纳入时间维度，与本章的空间对流项处理互为补充。