有限体积法的基础
理论：把守恒律写成“盒子的收支”
概述
听说CFD软件几乎都用有限体积法。它和结构分析用的有限元法到底有什么不同？
两者的出发点不同。有限元法从微分方程乘以权函数、在整个区域积分的“弱形式”出发。有限体积法把空间划分为小盒子（控制体、单元），直接写出收支：“盒内物理量的增加 = 从面流入量 − 流出量 + 源项”。相邻盒子通过共享面交换完全相同的量，所以质量不会在某处凭空消失或增加。这是它在流体计算中受青睐的最大原因。
积分形式的输运方程
在密度 $\rho$、速度 $\mathbf{u}$、扩散系数 $\Gamma$ 的场中，任意物理量 $\phi$（速度分量、温度、浓度等）在控制体 $V$ 及其表面 $S$ 上的守恒写作：
从左到右依次是时间项、对流项、扩散项和源项。对单元 $P$ 把面积分换成各面 $f$ 的求和，得到离散方程（$N$ 为面另一侧的相邻单元，$F_f$ 为通过该面的质量流量）：
公式看起来很简单，难点在哪里？
未知量位于单元中心，但方程中出现的是面上的值 $\phi_f$和面上的梯度。如何从单元中心插值到面，几乎决定了有限体积法的精度与稳定性，后面几节都围绕这一点展开。
单元中心型与节点中心型
| 方式 | 未知量位置 | 主要应用 |
|---|---|---|
| 单元中心型 | 单元中心；单元即控制体 | OpenFOAM, Ansys Fluent, Simcenter STAR-CCM+ |
| 节点中心型 | 网格节点，在节点周围构造对偶单元 | 如 Ansys CFX |
推广传热与流动计算的一本书
常被认为让有限体积法走进工程实践的，是 S. V. Patankar 的教科书《Numerical Heat Transfer and Fluid Flow》（1980年）。书中以可以手算跟踪的方式，讲解了由控制体收支组装离散方程的步骤，以及交替修正压力和速度的 SIMPLE 算法（Patankar 与 Spalding 于1972年发表）。如今商业软件内部运行的许多思路，都延续自这本书。
离散：如何确定面上的值
扩散项：中心差分与非正交修正
扩散项的面梯度，基本用两侧单元中心之差 $(\phi_N-\phi_P)/|\mathbf{d}|$ 近似。只有连接单元中心的向量 $\mathbf{d}$ 垂直于面时才准确；在斜交（非正交）网格上，需要加上利用单元中心梯度的修正项。修正项多为显式处理，非正交性大时迭代容易不稳定。
对流项：格式的选择
| 格式 | 精度 | 特性 |
|---|---|---|
| 一阶迎风 | 一阶 | 始终有界且稳定；数值扩散大，梯度被抹平 |
| 中心差分 | 二阶 | 网格佩克莱数超过2时振荡 |
| 二阶（线性）迎风 | 二阶 | 用上游梯度外推；实用标准；可能轻微过冲 |
| QUICK | 三阶（均匀网格） | 精度高但不保证有界 |
| TVD（带限制器） | 二阶（突变处降为一阶） | 仅在急剧变化处切换为迎风以防振荡 |
网格佩克莱数
对流与扩散孰强孰弱，可以用以单元宽度 $\Delta x$ 为特征长度的佩克莱数判断：
在中心差分中，连接下游相邻单元的系数为 $D - F/2$（$D=\Gamma/\Delta x$，$F=\rho u$）。当 $Pe_\Delta>2$ 时它变为负值，方程中就混入了“上游升高、下游反而降低”的非物理关系。这就是振荡的根源。
算例：一维对流扩散
长度 $L=1$ m、$\rho=1$ kg/m³、$\Gamma=0.1$ kg/(m·s)、左端 $\phi=1$、右端 $\phi=0$ 的稳态对流扩散问题有精确解，非常适合检验格式的特性：
流速 $u=2.5$ m/s（$Pe=25$）、5个单元（网格佩克莱数为5）的结果如下：
| 位置 x [m] | 精确解 | 中心差分 | 一阶迎风 |
|---|---|---|---|
| 0.1 | 1.0000 | 1.0356 | 0.9998 |
| 0.3 | 1.0000 | 0.8694 | 0.9987 |
| 0.5 | 1.0000 | 1.2573 | 0.9921 |
| 0.7 | 0.9994 | 0.3521 | 0.9524 |
| 0.9 | 0.9179 | 2.4644 | 0.7143 |
同一条件下增加单元数时的最大误差如下：
| 单元数 | 网格佩克莱数 | 中心差分 | 一阶迎风 |
|---|---|---|---|
| 5 | 5 | 1.547 | 0.204 |
| 10 | 2.5 | 0.537 | 0.158 |
| 20 | 1.25 | 0.160 | 0.120 |
| 40 | 0.625 | 0.044 | 0.079 |
| 80 | 0.31 | 0.012 | 0.048 |
| 160 | 0.16 | 0.003 | 0.026 |
5个单元的中心差分，边界值明明是0和1，却算出2.46……完全坏掉了。
对，而且线性方程组本身解得很好，看残差根本发现不了。迎风格式虽不振荡，但单元数加倍误差只减小约1.7倍（一阶精度）。中心差分在网格佩克莱数低于2后，每加密一倍误差约减小4倍（二阶精度），从40个单元起反而更准确。实务中以二阶迎风类格式和TVD为标准，正是为了同时避开这两种缺点。
压力—速度耦合
不可压缩流动中没有直接决定压力的方程。求解器先由动量方程得到预测速度，再求解压力修正方程使连续性方程成立，并修正速度和压力，如此反复。稳态计算常用SIMPLE算法（及改进版SIMPLEC），非稳态计算常用PISO算法。
- 速度与压力同置于单元中心时，压力会出现棋盘状振荡。通过 Rhie–Chow 插值构造面通量可以避免，主要商业软件都在内部完成了这一处理。
- SIMPLE 类算法需要欠松弛。压力0.3、动量0.7是常用的初始值；发散时调低，收敛过慢时调高。
实务：网格与求解器设置
网格质量的判读
| 指标 | 影响 | 参考 |
|---|---|---|
| 非正交性 | 扩散项精度与迭代稳定性 | OpenFOAM 的 checkMesh 在超过70°时警告；需调整修正次数和限制 |
| 偏斜度 | 面值插值精度 | checkMesh 在超过4时警告 |
| 相邻单元体积比 | 梯度评估精度 | 避免尺寸突变，逐步过渡 |
| 长宽比 | 与流动方向不一致时精度下降 | 边界层中可以较大（沿流向拉长的单元） |
各软件的阈值不同，这些只作参考，最终依靠下面的收敛检验来判断。
收敛判定
- 不要因为残差下降三个数量级就结束。要确认压降、升力、出口温度等关心的量不再随迭代变化。
- 确认流入与流出质量流量之差相对于流量足够小（例如约0.1%以下）。
- 最终结果用二阶格式给出。一阶迎风仅用于起步稳定，中途再切换。
主要软件中的设置
| 软件 | 对流格式设置 | 压力—速度耦合 |
|---|---|---|
| OpenFOAM | 在 system/fvSchemes 的 divSchemes 中写入 Gauss upwind、Gauss linearUpwind grad(U)、Gauss linear 等 | 由求解器决定（simpleFoam＝SIMPLE，pimpleFoam＝PIMPLE） |
| Ansys Fluent | 在空间离散中选择 First Order Upwind、Second Order Upwind、QUICK 等 | SIMPLE / SIMPLEC / PISO / Coupled |
| Simcenter STAR-CCM+ | 在物理模型中选择对流离散阶数 | 分离式（Segregated）或耦合式（Coupled）求解器 |
斜向流动中出现的“伪扩散”
众所周知，一阶迎风的数值扩散在流动与网格成45°斜交时最大。即使物理上完全没有扩散，斜向运动的尖锐界面在结果中也会变得模糊、变宽。若这种“伪扩散”大于物理扩散系数，结果反映的就是网格而不是物理。在评估混合或热扩散的计算中尤其要注意。
常见问题与验证
症状与原因
| 症状 | 可能原因 | 对策 |
|---|---|---|
| 超出边界值的数值、条纹 | 网格佩克莱数大却用中心差分 | 加密网格，或改用线性迎风、TVD |
| 几十次迭代即发散 | 高非正交单元、松弛因子过大 | 修正相关单元、降低松弛因子、先用一阶起步 |
| 残差停在某一数值 | 本质非稳态的流动（分离、涡脱落）被当作稳态求解 | 若监测量在振荡，改为非稳态计算（CFD的时间积分法） |
| 更换格式结果变化很大 | 网格过粗，数值扩散占主导 | 先进行网格收敛检验 |
验证方法
- 在像上面一维问题那样有精确解的问题上，确认格式的精度阶数，即网格加密一倍误差减小到几分之一。
- 在实际几何上，至少用三套网格比较关心的量，并用 Richardson 外推或 GCI（网格收敛指数）估算离散误差。
- 在掌握离散误差之后再与实验比较。把依赖网格的结果凑到与实验一致，换个工况就会失准。
我想自己改变佩克莱数，重现这种振荡。
用对流扩散方程模拟器改变流速和扩散系数，就能找到中心差分开始振荡的分界。控制体的组装方法可以用一维有限体积法模拟器查看，由此得到的三对角方程组的解法可以用托马斯算法（TDMA）模拟器跟踪。各格式的详细比较整理在迎风格式一文中。