第 8 章 · 从欧拉方程到完整求解器

前几章的零件都齐了：方程（第 2 章）、差分（第 4 章）、稳定性（第 5 章）、激波处理（第 7 章）。这一节把它们组装成一台完整的多维欧拉求解器，走的是 CFD 开发者的标准路线。

时间推进法的总体结构

可压缩流的定常解常被当作"非定常问题跑到时间趋于无穷"来求——时间推进法。每一时间步： 1. 由当前流场算通量（无粘通量用迎风或中心 + 人工粘性）； 2. 按格式（显式 MacCormack / 隐式 Beam-Warming）更新守恒量 $(\rho, \rho u, \rho v, e)$； 3. 边界条件：远场用特征线理论定入流/出流性质，物面法向速度为零。

多维问题用维分裂（dimensional splitting）降维：一个二维步拆成两个一维步，$x$ 方向推半步、$y$ 方向推半步。格式分析只需要在一维做，工程实现大大简化。

边界条件：一半的坑在这里

• 物面：无粘流法向速度为零 $\mathbf{V}\cdot\mathbf{n}=0$，切向滑移；粘性流全部速度为零（无滑移）。
• 远场：用一维特征理论判断——进入计算域的特征数决定要指定几个物理量，剩下的由相容关系外推。亚声速出流指定压力、超声速出流全部外推，这些规则的根源都在特征值符号。
• 边界处理不当的典型症状：解整体收敛了，但边界附近持续小幅振荡，或升力阻力算错。

收敛判据与加速

定常解的收敛判据用残差（守恒方程右端不平衡量）的量级下降 3-4 个数量级。加速手段：局部时间步（每个网格用自己的最大允许 $\Delta t$，只适用于求定常解）、多重网格、隐式残差光顺。

老师的提醒

• 从一维 Sod 到二维圆柱绕流再到三维，难度不是线性增长：二维以上的边界条件、网格生成、诊断手段各是独立的坑。
• 写求解器第一天就该写的工具：残差打印、质量/动量/能量的全局守恒检查（通量积分差应为零）。

自检

• 时间推进法求定常解的逻辑是什么？
• 远场边界"进入几条特征就指定几个量"的原理？
• 维分裂牺牲了什么？

更深一层：把"零件"组装成求解器的工程学

从一章一个零件到一台完整求解器，中间隔着的是架构决策。第一层：数据结构——结构网格用 $(i,j,k)$ 数组加邻接偏移（缓存友好、快），非结构网格用面-单元拓扑表（几何自由、访存散）；第二层：时间推进框架——显式 Runge-Kutta 简单可并行但受 CFL 拘束，隐式近似因式分解（Beam-Warming）无条件稳定但每步要解块三对角系统；第三层：收敛加速——局部时间步、残差光顺、多重网格、隐式残差平滑，四件套通常叠加使用。每一层的选择互相牵制：选了非结构网格，多重网格就得用聚合式而非几何式；选了显式推进，隐式边界条件就麻烦。

求解器开发者的真实日常不是推公式，而是调试。一套最小调试纪律值得全行业背诵：先在一维均匀网格上复现已知精确解；再加几何复杂性（曲线网格、贴体坐标）；再加物理复杂性（粘性、湍流模型）；每加一层都与上一层的结果做回归对比。求解器是逐层长出来的，不是一次写成的——跳层开发的人最终会在第 N 层发现第 2 层的 bug，然后从头再来。

推导细读：边界条件的特征计数为什么是"铁律"

以一维欧拉方程为例，雅可比矩阵特征值 $u - a, u, u + a$。出流边界处逐个数指向域外的特征：超声速出流三条全出（外推一切）、亚声速出流一条入（指定压力、外推其余）、亚声速入流两条入（指定速度与熵/温度、外推压力）、超声速入流三条全入（指定全部来流条件）。这个计数的深层依据是适定性：指定的量数 = 入流特征数，多指定则超定（无解或边界层病态），少指定则欠定（解漂移）。工程上边界条件错误的典型症状：收敛后升阻力仍缓慢漂移、远场出现反射波纹、残差卡在某个水平锯齿震荡——三条症状对照特征计数表，十分钟内可定位病灶。

在真实世界里：一架机翼算例的完整决策树

商业软件里跑机翼极曲线，决策链长这样：马赫数 < 0.3 → 压力基不可压；0.3–0.8 → 密度基/可压缩修正如需；跨声速 → 密度基 + 激波捕捉格式 + 更细的尾迹网格。湍流模型：附着流为主选 Spalart-Allmaras（一方程、鲁棒），有大分离选 k-ω SST（对逆压梯度更诚实）。网格：物面第一层高度按 $y^+ \approx 1$ 反推（粘性底层解析）、远场 50 倍弦长、尾迹区加密三倍。收敛判据：残差降 4 个量级 + 升力线平台 + 力矩不漂。这条决策树的每个分叉背后都是本书某一节的内容——学完回头再看，决策树就是目录。

常见误区与辨析

• "残差降到底就是对的"：残差只说明离散方程被满足，不能说明方程组（含模型）正确——升力系数还要与网格无关解和实验数据双重对照。
• "边界离得远就不用讲究"：远场边界处理不当的反射误差随距离衰减很慢，压力基求解器常用 20–100 倍弦长；偷近了升力能差百分之几。
• "并行加速比是软件好坏的标志"：通信模式（分区质量）与问题规模同样重要；小算例的加速比差不能定罪软件。

专题深入：隐式求解器——"大时间步"的代价与回报

隐式时间推进（Beam-Warming 近似因式分解）的账本是单步贵、步数少：每步解块三对角方程组（$O(n)$ 大矩阵的带状求解），但 CFL 可到 10-100（甚至定常问题可到 1000）。回报的条件：定常或慢瞬态问题（时间精度要求低）——隐式的自由度在长步长。代价：每步内存与算力是显式的数倍、并行效率差（三角求解的依赖链）、非定常精度受限。显式 vs 隐式的分工：激波捕捉与湍流 LES 用显式（时间精度刚性）、定常 RANS 用隐式（大步长）——时间尺度决定算法。

多重网格与隐式/局部时间步的叠加是定常求解的"加速三件套"：局部时间步（每格用自己的最大步长）、隐式残差光顺、多重网格——收敛速度的数量级提升来自三件套的叠加而非单一改进。

数学细节展开：边界条件的"特征变量"实现

边界条件的正确实现走特征变量路线：把守恒变量变换到特征空间（左特征向量对角化），按特征方向指定或外推，再变回守恒变量。亚声速出流的完整实现：一条特征（$u - a$）指向域外（来自域内信息），两条指向域内（来自边界条件）——指定压力 + 外推速度与密度，特征变换后各方向的"指定/外推"分配严格按特征符号。用中心差分硬写边界条件（不区分特征方向）是边界反射波污染全场的根源——边界条件的特征实现是双曲问题数值方法的最高级细节。

工程案例细节：升阻力系数验收的行业标准

机翼 CFD 的升阻力系数验收有行业标准：网格收敛性（GCI 指数——三套网格的外推）、基准对照（风洞数据或文献数据）、残差与监测量（升力线平台 + 力矩不漂）、物理合理性（压力分布光滑、无虚假分离）。验收标准是行业审计的法定语言——没有 GCI 的 CFD 报告在工程审核中没有可信度。