第 16 章：可压缩流体流动的计算（Fluid Flow Computation: Compressible Flows）
16.1 历史（Historical）
CFD 方法传统上分成两大族：density-based（密度基）和 pressure-based（压力基）。density-based 方法长期主导航空工业中跨声速与超声速流动的模拟；这一格局在 SIMPLE 算法首次提出时已经基本成型。SIMPLE 是一种 pressure-based 方法，最初正是为了弥补密度基方法在低速区域上的不足而被开发的，并且在解决不可压缩流动与低 Mach 数流动上表现得相当高效。从历史上看，这一分工并非偶然，而是与两类流动的数学性质直接相关：跨声速与超声速流动中密度变化显著，密度基方法通过把密度作为主变量能够自然地捕捉激波；而低速不可压缩流动中密度几乎不变，密度作为主变量反而会因为方程组的病态（特征速度趋于零）而难以收敛。
早期的工作集中于把任一族方法扩展到对方传统擅长的流型。Harlow 与 Amsden 是最早尝试模拟全速度域流体流动的研究者之一，其文献 [1, 2] 即本章参考文献列表中的相关条目。在其工作中，将压力作为主变量（相对于密度）的做法被提出来作为优势，理由是压力的变化幅度不论 Mach 数大小都能保持在有限范围内——这避免了密度在低 Mach 数时变得对扰动极敏感的问题。然而真正提供清晰解决方案的，是 Patankar 的工作。Patankar 的研究使基于 SIMPLE 的方法得以真正发展为能够处理全速度域流体流动的方法。最关键的进展是对压力方程进行了重新构造，将密度修正与速度修正纳入其中，使方程的类型从不可压缩流动的纯椭圆型，转变为跨声速与超声速可压缩流动的双曲型。这一变化使 SIMPLE 系列方法能够横跨整个 Mach 数谱进行求解，其中压力扮演着双重角色——在 Mach 数较高（高度可压缩）的极限下影响密度，在 Mach 数较低（不可压缩）的极限下影响速度——以实现质量守恒。这一段历史回顾为后续各节奠定基调：本章所介绍的可压缩 SIMPLE 算法正是这一历史线索的具体落地，它把 Patankar 当年的压力方程改造思路形式化为今天所见的离散压力修正方程。
16.2 引言（Introduction）
pressure-based 方法的一个重要优点是，它能在不同的 Mach 数流型下求解流体流动而无需任何为促进收敛或稳定计算而添加的人为处理。这种能力源自压力在可压缩流动中所起的双重作用。这种"无需人为处理"的特性是 pressure-based 方法相对 density-based 人工压缩性技术的重要优势——后者需要精心设计预条件子才能在低速下保持数值稳定，而 pressure-based 方法只需要在算法结构上自然地把压力与密度的耦合纳入方程即可。这一双重作用可以通过考察以下两个极端情形来最佳地描述。
第一种情形是 Mach 数非常低时。此时通过动量守恒建立流场所需的压力梯度非常小，以至于它不会显著地影响密度，流动可以被视为不可压缩的。因此密度-压力之间以及密度-速度之间的关联都极弱，意味着密度的变化对速度变化并不敏感。在这种情况下，连续性方程不再能被视作密度的方程，而是起到约束速度场的作用——换言之，连续性方程不再独立地决定密度，而是作为对速度场的一个约束（不可压缩条件 \(\nabla \cdot \mathbf{v} = 0\)）。这一观察与第 15 章中不可压缩 SIMPLE 算法的设计完全一致：压力仅通过修正速度场来间接修正密度。
第二种情形是高超声速。此时速度的变化相对于速度本身的量值来说变得较小，意味着压力的变化会显著地影响密度。于是压力通过状态方程直接作用于密度以满足质量守恒，而连续性方程可被视为密度的方程。此时连续性方程的物理角色发生反转：从"约束速度场"变为"求解密度场"。这种角色转换正是 Mach 数从低到高跨越声速时压力方程从椭圆型变为双曲型的物理根源。
上述两种极限情形揭示了压力在可压缩流动中扮演的双重角色。它清楚地表明，压力通过状态方程作用于密度场、通过动量方程中的梯度作用于速度场，从而实现质量守恒。这种双重角色解释了为何 pressure-based 方法能够成功预测全速度域的流体流动——核心在于：只要在压力修正方程中正确地同时纳入这两种作用（通过密度修正 \(C_q p'\) 和速度修正 \(-\mathbf{D}_v \nabla p'\)），算法就能在不同的 Mach 数极限下自动地呈现出相应的物理行为，而不需要为每个流型单独调整算法结构。
然而这一事实并未阻止 density-based 阵营的研究者采用人工压缩性技术来开发可求解全速度域流体流动的方法。为了克服由此产生的刚度矩阵带来的性能退化，研究者引入了对所得刚度矩阵的预条件处理，并陆续有若干使用该技术的方法出现在文献中。类似地，也有多种 pressure-based 方法被开发出来用于预测全速度域流体流动，分别采用交错网格或同位（collocated）变量格式。在这些方法中，有些采用原始变量作为因变量，有些则以动量分量作为因变量。还有研究者采用流向密度阻滞（stream-wise directed density-retardation）的概念，由依赖于 Mach 数的监控函数控制，以便在跨声速和超声速流型中反映守恒律的双曲特性——这一思路的核心是用 Mach 数作为权重，在低 Mach 区抑制密度变化、在高 Mach 区恢复密度变化。另一些技术则在高 Mach 数时采用一阶迎风格式来计算控制体面上的密度、在低 Mach 数时采用中心差分格式——这相当于根据 Mach 数自动切换密度插值格式。
本章将在前一章发展的同位压力基技术的基础上进行扩展，使其能够模拟所有 Mach 数值的流体流动。所采用的方法易于实现、精度高，且不需要任何额外显式阻尼项来增强稳健性或精确地分辨激波。这一承诺的具体兑现方式将在后文的推导中逐一展现：在压力修正方程中通过 \(C_q p'\) 项引入密度修正、在动量方程中加入体粘度项、在能量方程中完整地处理所有可压缩特有的项，以及在边界条件中针对 Mach 数分别设计 inlet/outlet 的处理方式。
16.3 守恒方程（The Conservation Equations）
求解可压缩流体问题的守恒方程包括连续性方程、动量方程与能量方程，三者必须同时求解，这与第 15 章中不可压缩流动只解连续性与动量、能量方程通常被忽略掉的做法有本质区别。对于表现为理想气体的牛顿流体，这些方程可以写成如下形式。
连续性方程（方程 16.1）为
动量方程（方程 16.2）为
其中体粘度相关项 \(- \nabla[(2/3)\mu(\nabla \cdot \mathbf{v})]\) 仅在可压缩情形下保留下来（不可压缩时该散度为零）。
此处所采用的能量方程（方程 16.3）形式是用温度表示的版本：
方程右端除扩散项 \(\nabla \cdot (k \nabla T)\) 外，还包括：与比热变化率 \(D c_p/D t\) 相关的项、与压力物质导数 \(D p/D t\) 相关的项（带 \(2/3\) 系数，对应体粘度的可压缩贡献）、体粘度耗散项 \(\mu W\)（其中 \(W = \nabla \cdot \mathbf{v}\)）、粘性耗散项 \(\mu U\)（\(U\) 为应变率张量的二阶不变量），以及单位体积热源/汇 \(\dot{q}_V\)。这些项在第 15 章的不可压缩能量方程中没有出现，是可压缩情形所特有的。
上述方程组应补充一个联系密度与压力、温度的状态方程 \(\rho = \rho(p, T)\)，对于理想气体可写成 \(\rho = p/(R T)\)（方程 16.4），其中 \(R\) 是气体常数。该状态方程把流体动力学量（速度、压力）与热力学量（温度）耦合起来，是流体动力学与热力学之间桥梁的形式化表达。在后续的离散推导中，这一耦合将通过密度修正场 \(\rho' = C_q p'\) 与压力修正 \(p'\) 直接关联起来。
在后文的推导中，上标 \(n\) 表示迭代开始时的取值，上标 \(*\) 表示迭代过程中已更新一次的取值，上标 \(**\) 表示同一次迭代中已更新两次的取值。这套标记约定与第 15 章中的不可压缩情形保持一致，但内涵因可压缩效应而变得更丰富——\(*\) 与 \(**\) 之间的差异不再仅由速度修正贡献，还包括密度修正与压力修正的耦合。
16.4 动量方程的离散（Discretization of the Momentum Equation）
对方程 (16.2)（动量方程）在图 16.1 所示的控制体 \(C\) 上进行离散，所得到的形式与第 15 章中给出的不可压缩形式基本相同。两者之间只有两处差异，分别涉及密度的插值到控制体面上的方式以及新增的 \(-\nabla[(2/3)\mu(\nabla \cdot \mathbf{v})]\) 项——该项含体粘度。这两处差异虽然"局部"，但却折射出可压缩流与不可压缩流在动量方程层面上的本质差别。
先看第一处差异：在可压缩流动中密度不再是常量，由于它被存储在控制体质心处，要在控制体面上计算质量流率就需要对它进行插值。采用线性插值剖面（即中心差分）会在高速情况下引起振荡——具体地说，在激波附近密度梯度很大，线性插值会在激波前后产生非物理的"过冲"（数值振荡）。因此应该使用有界（bounded）的迎风偏置方案来完成这一任务。第 11 章和第 12 章中介绍的任一种有界对流格式都可以采用，包括 MINMOD、SMART、STOIC、MUSCL、SUPERBEE 等 TVD 类格式以及 NVF 类的 Q（universal limiter）格式等。在 OpenFOAM® 实现中，对应的处理方式是用 linearInterpolate(rho*U) 或在 Rhie-Chow 步骤中对密度采用迎风格式 upwind<scalar>(mesh,mDot).interpolate(rho) 来得到面上的密度。
第二处差异是涉及 \(\nabla(\mu \nabla \cdot \mathbf{v})\) 的新增项。该项此前尚未被离散化，它的离散形式是利用方程 (2.85) 得到的：根据该方程，一个标量梯度的体积分被转化为面积分，进而转化为控制体面上通量的求和
面上的速度向量散度按下式计算
其中 \(\phi = u\)、\(v\) 或 \(w\) 的梯度按下式线性插值到面上
上式中的权重 \(g_C\)、\(g_F\) 来自标准的线性插值系数。应当注意，\(\nabla \cdot \mathbf{v}\) 在每个面上都被独立计算（即对每个分量单独插值再求和），这避免了把散度作为一个整体插值可能引入的不一致性。
动量方程的最终离散形式由方程 (15.70) 给出，其系数由方程 (15.71) 给出，再加上一项
进入源项。这相当于把体粘度贡献直接作为源项加到动量方程的代数形式中，而不需要修改主对角元或邻接系数——这一处理方式对应于 Stokes 假设下粘度对称张量的常用近似。
与不可压缩流动的情形类似，代数方程被欠松弛（under relaxed），并写成方程 (15.78) 的形式，该形式适合用于压力修正方程的推导。欠松弛在可压缩情形下通常更为关键：高速流动中 \(\dot{m}^*_f\) 的数值往往较大，源项更新幅度也较大，欠松弛可以避免迭代初期的过冲与发散。
16.5 压力修正方程（The Pressure Correction Equation）
可压缩流动的压力修正方程可以通过对不可压缩情形下相应方程的简单扩展来得到。两者的差异与密度的变化有关：可压缩情形下需要定义一个密度修正场 \(\rho'\)，并通过一个压力-密度关系把它与压力修正场联系起来。然而，由此带来的实质差异会体现在边界条件的处理上，这些将留到本章后面再解释。
对于理想气体，压力与密度之间的关系写成 \(\rho R T = p\)。利用这一关系，可以对密度修正与压力修正之间的关系作 Taylor 展开，得到
修正后的压力、密度、速度与质量流率场定义为
而以修正场表示的半离散连续性方程可以写成
其中
二阶修正项 \(\rho'_f \mathbf{v}'_f \cdot \mathbf{S}_f\) 通常被忽略，因为与其他项相比它要小得多。这一近似对收敛速率几乎没有影响，仅在求解过程最初的几次迭代中起作用。此外，最终的解并不受影响，因为收敛时修正场为零。
利用 Rhie-Chow 插值对 \(\mathbf{v}^*_f\) 与 \(\mathbf{v}'_f\)，\(\dot{m}^*_f\) 与 \(\dot{m}'_f\) 分别表示为
以及
此处同样略去二阶项。注意已将 \(\rho'_f\) 用 \(C_{q,f} p'_f\) 替代。方程 (16.14) 中带下划线的项与不可压缩算法中遇到的项一样困难，通常被略去。略去该项后，质量流率的修正变为
该式右端的第一项类似于不可压缩情形下出现的项，而第二项是新增的密度修正贡献。这一项很重要，它把压力修正方程从椭圆型方程转化为能够分辨超声速与高超声速情况下可能出现的激波的双曲型方程。这使得可压缩 SIMPLE 算法能够在不需要任何特殊预条件处理的情况下用于预测全速度域的流体流动。
通过对方程 (16.15) 作简单归一化处理可获得更深入的洞见：用 \(\dot{m}^*_f \cdot \mathbf{S}_f \cdot C_{p,f} \rho^*_f\) 去除方程两端，可使 \(p'_f\) 项前的权重为 1，并使 \(\nabla p'_f\) 项前的权重正比于 \(1/(M^2)\)（其中 \(M\) 为当地 Mach 数），即
对 Mach 数较低的流动，\(\nabla p'\) 修正项占主导地位，使方程回到椭圆形式，与不可压缩情形相同。而在 Mach 数很高的流动中，\(p'_f\) 修正项不能再被忽略，从而赋予修正方程以双曲特性。这两种行为的结合使该算法能够预测全速度域流体流动。
将方程 (16.14) 代入连续性方程 (16.11)，得到可压缩形式的压力修正方程：
对带下划线项的处理方式不同，便得到 SIMPLE 系列算法的不同变体。略去带下划线项，便得到 SIMPLE 算法的压力修正方程：
该方程可分解为 transient-like term、convection-like term、diffusion-like term 三部分与一个 source-like term 之和的形式 (16.18)，其中 \(a^0_{F,p}\) 与 \(a^0_{C,p}\) 的完整表达由方程 (16.19) 给出，源项 \(b^0_{C,p}\) 同样在 (16.19) 中显式给出，并包含一个通常被忽略的非正交项。
值得强调的是，convection-like 项是在压力修正方程的推导过程中自然出现的，它的存在对于算法分辨全速度域流动的能力至关重要。对于 Mach 数较高的流动，密度修正以对流形式传播（即呈双曲行为），描述这一现象的数学算子是一阶散度算子。因此，与不可压缩流动中只出现 diffusion-like 项从而使压力修正方程表现为椭圆特性不同，形式为 \(p' + C\) 的压力修正解不再满足方程。这一点说明：对不可压缩流动而言，无论在边界上指定何压力值都不会影响解；但对可压缩流动而言，必须在边界上明确给定压力的具体值，因为所选的值会影响最终的解。
还应注意，由于收敛时修正场为零，用于离散 convection-like 项的格式阶次对最终结果的精度没有影响。但 \(\dot{m}^*_f\) 的情况不是这样：在其计算中使用高阶格式的确能改善算法的激波捕捉特性。因此为了增强稳健性，对 convection-like 项的离散采用迎风格式是有益的。此外，按照第 15 章所述忽略 diffusion-like 项的非正交贡献，压力修正方程与其系数化为 (16.19) 所列的形式。
求得压力修正场 \(p'\) 后，按下列各式对速度、压力、密度与质量流率场进行修正：
其中 \(k_\rho\) 是密度的欠松弛因子。
16.6 能量方程的离散（Discretization of The Energy Equation）
能量方程 (16.3) 中的非稳态项、对流项与扩散项的离散遵循前面章节给出的一般过程，此处不再重复。
16.6.1 额外项的离散（Discretization of the Extra Terms）
本节重点讨论能量方程右端所出现的、在一般标量方程离散中未曾处理的新项，这些项专属于能量方程。它们之中有许多作为源项处理，在对其积分时取单元体质心处的取值，以保证二阶精度。
16.6.1.1 比热项（The Specific Heat Term）
包含比热项的离散化按下式进行：
16.6.1.2 物质导数项（The Substantial Derivative Term）
压力的物质导数项的离散化为
16.6.1.3 体粘度耗散项（The Dissipation Term）
含体粘度的耗散项的离散形式由下式给出：
16.6.1.4 粘性耗散项（The Viscous Dissipation Term）
粘性耗散项的离散化方式与体粘度项相似，其表达为
其中 \(U_C^{**}\) 包含
各项均在控制体质心 \(C\) 处取值，并乘以 \(V_C\)。
16.6.1.5 源/汇项（The Source/Sink Term）
单位体积内的热源/汇项被离散化为
上述各项的离散形式被代入能量方程，得到能量方程的代数形式，下一节将作介绍。
16.6.2 能量方程的代数形式（The Algebraic Form of the Energy Equation）
对非稳态项采用一阶 Euler 格式、对对流项在 deferred correction 框架下采用高分辨率格式，可以将能量方程的最终代数形式写成
其中系数为
与其他变量类似，能量方程的求解通常也需要欠松弛。
16.7 可压缩 SIMPLE 算法（The Compressible SIMPLE Algorithm）
同位可压缩 SIMPLE 算法的各要素如图 16.2 所示，可归纳为以下步骤。整个流程在结构上与第 15 章的不可压缩 SIMPLE 算法一致，主要差异在于：加入了能量方程的求解（第七步），以及在压力修正方程中处理密度修正贡献（第四到第六步）。
第一步，为计算 \(t + \Delta t\) 时刻的解，以 \(t\) 时刻的收敛值 \(p^{(n)}\)、\(\mathbf{v}^{(n)}\)、\(\rho^{(n)}\)、\(T^{(n)}\)、\(\dot{m}^{(n)}\) 作为初始猜测。这一步继承了第 15 章的做法，但因为可压缩情形下还需要温度与密度的初值，所以初猜的字段更多。
第二步，求解由方程 (16.2) 给出的动量方程以得到新的速度场 \(\mathbf{v}^*\)。动量方程的离散形式与第 15 章相同，但其中密度的插值与体粘度项需要按 16.4 节的处理方式给出。
第三步，利用状态方程 \(\rho = p/(RT)\) 计算新的密度场 \(\rho^*\)。这一步是新增的——不可压缩情形下密度为常量、不需要从压力反推密度；而可压缩情形下压力的每次更新都会立即改变密度。
第四步，用 Rhie-Chow 插值技术（方程 16.13）更新控制体面上的质量流率，得到满足动量守恒的质量流率场 \(\dot{m}^*\)。这一步在不可压缩与可压缩两种情形下形式相同，但插值中涉及的密度取值不同。
第五步，利用新的质量流率计算压力修正方程 (16.19) 的系数并求解，得到压力修正场 \(p'\)。系数中包括 convection-like 项（来自密度修正）和 diffusion-like 项（来自速度修正），这是与不可压缩情形最关键的差异。
第六步，利用方程 (16.20)–(16.23) 更新控制体质心处的压力、密度、速度场以及控制体面上的质量流率，得到满足连续性的场。这组修正公式体现了三个变量（压力、密度、速度）的协同更新——压力被加上 \(k_p p'\)，密度被加上 \(k_\rho C_q p'\)，速度被减去 \(\mathbf{D}_v \nabla p'\)，三者通过 \(p'\) 共同收敛到满足质量守恒的状态。
第七步，求解能量方程 (16.29)–(16.30) 得到新的温度场 \(T^*\)。能量方程在每一外迭代中也必须求解，以便温度场能够反过来影响下一轮动量方程中的粘度等热物性参数。
第八步，将 \(\mathbf{v}^{**}\)、\(\dot{m}^{**}\)、\(\rho^{**}\)、\(T^*\)、\(p^*\) 分别设为下一迭代中速度、质量流率、密度、温度与压力的初始猜测。第九步，回到第二步并重复，直到收敛。第十步，将 \(t + \Delta t\) 时刻的解设为收敛解。第十一步，把当前时间推进到下一时刻 \(t + \Delta t\)。第十二步，回到第一步并重复，直到达到最后一个时间步。
整个流程的关键在于第二步到第七步构成一个"内迭代"，每完成一次该循环，五个场都被更新一次；连续执行多次直到所有修正场为零，即得到当前时刻的稳态（或瞬态）解。需要注意的是，密度、温度、压力之间通过状态方程和能量方程形成耦合：压力更新 → 密度更新 → 状态方程一致性检查 → 能量方程更新温度 → 温度反过来影响粘度等物性，这构成了一个完整的可压缩流求解回路。
16.8 边界条件（Boundary Conditions）
一般而言，对动量方程的边界条件处理在不可压缩与可压缩流动之间没有差别。因此动量方程方面所需的修改即第 15 章所讨论的（壁面滑移、no-slip、symmetry、inlet/outlet 的速度指定等），此处不再重复。然而，压力修正方程方面则出现实质性差异，这也构成本节的主要内容。能量方程的边界条件遵循一般标量变量 \(\phi\) 的处理方式（包括 inlet、outlet、Dirichlet、von Neumann、对称条件），同样不再重复。本节专注于压力修正方程在 inlet 与 outlet 处相对于不可压缩情形的增量。
对于边界控制体（如图 16.3 所示），连续性方程写成
其中边界面的贡献被单独标出，\(\dot{m}^*_b\) 表示边界质量通量、\(\dot{m}'_b\) 表示其修正。在离散推导中，把内表面与边界面分开处理有助于更清晰地追踪每一步贡献的来源。
采用与不可压缩情形相同的 Rhie-Chow 插值在边界面上的表达，边界面的速度、质量流率与质量流率修正分别表示为
这些表达式与第 15 章所给出的对应表达式之间的唯一区别是质量流率修正方程 (16.34) 多了一项
，它源于密度修正与压力修正的耦合——这一项在不可压缩情形下不存在（因为密度不变、\(C_{q,b} = 0\)），是可压缩情形下边界修正的"标志性"增量。
此外，在壁面与对称平面上，不可压缩与可压缩流动的边界条件实现没有差别——因为壁面与对称平面上 \(\dot{m}_b = 0\)（无论密度是否变化），密度修正项 \(\rho'_b = C_{q,b} p'_b\) 仍然为 0 或被自然抑制。第 15 章中关于壁面与对称边界条件对压力修正方程的边界元素所做修改在此同样适用，无须再述。
剩下要讨论的是 inlet 与 outlet 的边界条件。对可压缩流动，所需施加的条件由 Mach 数值决定。对无粘流动，当流动从亚声速变为超声速时，方程的数学类型从椭圆型变为双曲型。这意味着在入口/出口处，所需的边界条件数量也随之变化——亚声速边界需要指定较少的变量（方程是椭圆型，只有部分特征线指向域内），而超声速边界需要指定全部变量（方程是双曲型，全部特征线都从边界发出）。其实现细节将在下文给出。
16.8.1 入口边界条件（Inlet Boundary Conditions）
在域的入口处流动可能是亚声速或超声速，需要分别处理，因为流动方程可能属于椭圆型或双曲型。这一区分的本质是特征线方向在亚声速与超声速情形下的差异：亚声速时部分特征线指向域内、部分指向域外，因此只有部分变量需要在入口处被指定；超声速时全部特征线都从边界出发（即全部从域外向域内传播），因此所有变量都需要在入口处被指定。下文按这两个子情形分别给出处理方式。
16.8.1.1 入口亚声速流动（Subsonic Flow at Inlet）
在亚声速情况下，入口处可以施加多种条件，包括指定速度、指定静压与速度方向，或指定总压与速度方向。当域内可能过渡到超声速时，应使用最后一种类型。
第一种，指定速度（\(p_b = ?\)、\(\dot{m}_b = ?\)、\(\mathbf{v}_b\) 指定）。与不可压缩流动不同，对于可压缩流动，由于密度依赖于压力，即使在入口指定了速度，质量通量 \(\dot{m}'_b = \rho'_b \mathbf{v}^*_b \cdot \mathbf{S}_b \neq 0\) 仍然未知。在入口边界，乘压力修正 \(p'_b\) 的系数由 \(a^0_{p_b} = C_{q,b} \dot{m}^*_b / \rho^*_b\) 给出。为了实现压力修正方程，\(p'_b\) 用内部节点表示，相应地修改这些节点的系数。对于常值剖面的情形（即 \(p_b = p_C\)），\(a^0_{C,p}\) 系数由方程 (16.36) 给出，其中包括内部面贡献与边界面贡献之和。\(p_b\) 的取值仍然按照不可压缩流动一章所述由内部外推得到。
第二种，指定静压与速度方向（\(p_b = p_{specified}\)、\(\hat{\mathbf{e}}_v\) 指定、\(\dot{m}_b = ?\)、\(\mathbf{v}_b = ?\)）。在指定静压的情形下，\(p_b\) 已知，因此 \(p'_b\) 被置零，从而 \(\rho'_b\) 也为零。因此其实现方式与不可压缩情形相似，入口被处理为 Dirichlet 边界条件。在已知速度方向的情况下，其大小按不可压缩情形下的方法用方程 (16.33) 计算，得到与方程 (15.137) 相似的表达式。压力修正方程的系数由方程 (16.37) 给出。
第三种，指定总压与速度方向（\(p_{0,b} = p_{0,specified}\)、\(\hat{\mathbf{e}}_v\) 指定、\(\dot{m}_b = ?\)、\(\mathbf{v}_b = ?\)）。对这种情形，边界处速度的大小与压力均未知，但通过总压方程相关联：
其中下标 \(b\) 表示边界，\(p_{0,b}\) 为总压，\(p_b\) 为静压，\(\gamma\) 为比热比，\(M_b\) 为 Mach 数
方程 (16.37) 可改写为用总压表示的静压表达式
其中 \(\hat{\mathbf{e}}_v\) 为速度方向单位向量。对该式关于 \(\dot{m}^*_b\) 求导，得到
将其代入方程 (15.163)，可得到一个压力修正作为质量通量修正函数的方程：
把方程 (16.42) 给出的等价表达式替换方程 (16.34) 中的 \(p'_b\)，质量通量修正变为
将方程 (16.43) 给出的 \(\dot{m}'_b\) 代入展开的连续性方程，便可得到修正后的边界控制体系数：
最后还应提到，亚声速入口的能量方程边界条件通常是指定静温 \(T_b\) 或总温 \(T_{0,b}\)。如果指定静温，则与 Dirichlet 条件类似；如果指定总温，则每次迭代中由总温方程
解出静温并作为已知量处理。因此施加的也是 Dirichlet 型边界条件。
16.8.1.2 入口超声速流动（Supersonic Flow at Inlet）
指定静压、速度与温度（\(p_b = p_{specified}\)、\(\mathbf{v}_b = \mathbf{v}_{specified}\)、\(T = T_{specified}\)）。在超声速入口处，必须指定所有变量（压力、速度与温度）。这等价于 Dirichlet 型条件，意味着 \(\dot{m}'_b = p'_b = 0\)。因此边界控制体的 \(a^0_{C,p}\) 系数化为方程 (16.46) 给出的形式，仅含内部面贡献而无边界面贡献。
16.8.2.1 出口亚声速流动（Subsonic Flow at Outlet）
指定压力（\(p_b = p_{specified}\)、\(\dot{m}_b = ?\)、\(\mathbf{v}_b = ?\)）。在亚声速出口处通常指定压力。因此压力修正 \(p'_b\) 置零，而质量流率修正由下式计算
由于 \(\mathbf{v}^*_b\) 未知，习惯上假设其方向与上游速度 \(\mathbf{v}^*_C\) 相同。压力修正方程中 \(a_C\) 系数的表达写成方程 (16.48) 给出的形式。对于能量方程则采用零通量 Neumann 型边界条件。
指定质量流率（\(\dot{m}_b = \dot{m}_{specified}\)、\(p_b = ?\)、\(\mathbf{v}_b = ?\)）。对于出口处指定质量流率的情形，\(\dot{m}'_b\) 为零，直接从压力修正方程中略去，边界元素的系数无需修改。在方程 (16.34) 中令 \(\dot{m}'_b\) 为零，可得到边界处压力修正作为边界控制体质心压力修正函数的表达式
据此可计算边界的压力与密度。对于能量方程同样采用零通量 Neumann 型边界条件。
16.8.2.2 出口超声速流动（Supersonic Flow at Outlet）
在超声速出口处，所有变量都不应被指定，压力、速度、密度与温度的值由域的内部外推得到。因此 \(\dot{m}_b\) 与 \(p_b\) 都从内部单元外推。这等价于在压力修正上施加 Neumann 边界条件，由此得到修正后的 \(a_C\) 系数形式 (16.50)。
16.9 计算指针（Computational Pointers）
16.9.1 uFVM
对可压缩流动，算法的主要修改出现在压力方程的构造过程中——需要加入 convection-like 项。这部分添加在 cfdAssembleMdotTerm 中，由 Listing 16.1 给出的脚本实现。脚本计算各面 local_mdot_f，并据此更新 local_FLUXCf1、local_FLUXCf2、local_FLUXVf 各项，把 \((m_f / \rho_f)(\partial \rho / \partial p) P'\) 形式的项加到压力修正方程的组装中。
另一处重要的修改在于边界条件的处理——现在必须考虑一个可变的密度场。例如超声速入口条件在压力修正方程中的实现如 Listing 16.2 所示，其形式与 Listing 16.1 相同，只是具体应用到边界面处理流程中。
16.9.2 OpenFOAM®
本节将 simpleFoam 扩展为可处理全速度域的可压缩流体流动。这需要以下修改：(i) 加入与连续性和动量方程同时求解的能量方程；(ii) 使用一个联系密度与温度、压力的状态方程；(iii) 对压力修正方程以及若干边界条件作必要的修改。所得代码记为 simpleFoamCompressible，许多扩展以补充 include 文件的形式加入，如 Listing 16.3 所示。其中 upwind.H 与 gaussConvectionScheme.H 用于强制对压力修正方程的对流项采用迎风离散；bound.H 类用于把变量限制在一定范围内。
createFields.H 现在包含密度场以及其他与可压缩流物理相关的变量与常数的定义。psiTermo 类（Listing 16.4）提供对 OpenFOAM® 库中热物理关系的访问，包括方程 (16.4) 所描述的理想气体定律。压力与焓场在 thermo 类（Listing 16.5）中定义，在 createFields.H 中以引用方式访问。
速度的定义与不可压缩版本相同，但质量通量在定义中显式包含密度（Listing 16.6）。在可压缩流动中，对某些变量（如密度与压力）设置物理极限可以增强稳健性，特别是在最初几次迭代中，以避免变量取到非物理值（如负密度或负压力）。因此可以在 case 定义中设定上下限（Listing 16.7），并在 createFields.H 中读入。
动量方程以稍作修改的语法定义，使其能反映密度与热物理性质关系。其线性化公式的语法在 Listing 16.8 中给出。第一个指令定义动量方程的有限体积离散（尽管采用向量形式实现，速度的三个分量仍按分离方式求解）。系统被隐式松弛后由迭代求解器求解。求得动量方程后得到速度场的一个新猜测，但该速度场不一定满足连续性方程。
为强制质量守恒，需要构造压力修正方程来修正速度。按照方程 (16.19)，所用的语法如 Listing 16.9 所示。为避免棋盘效应，mDot 质量通量场用 Rhie-Chow 插值计算，但此时还考虑基于热物性模型求出的密度场（Listing 16.10）。值得一提的是，密度用迎风格式插值到面，以模拟可压缩流动的双曲行为。压力修正方程被完整组装并用 Listing 16.11 所示的语法求解。
求得压力修正方程后，依赖于压力修正的变量被更新。对质量通量场按 Listing 16.12 更新，这与不可压缩情形下的通量修正相似；其中 flux() 函数（Listing 16.13）通过对矩阵的上下系数与单元值相乘来计算修正通量。最后，控制体质心处的速度、密度、压力按方程 (16.20)、(16.21)、(16.22) 更新，如 Listing 16.14 所示，其中 alphaP 是压力与密度更新的显式松弛因子 \(k_p\)，对 SIMPLE 求解器的稳定性是必要的。
为了反映可压缩效应，能量方程被引入，相关的温度被求解。在 OpenFOAM® 中，能量方程以比静态焓 \(h = C_p T\) 的形式给出（由方程 (3.61)），其求解过程如 Listing 16.15 所示。一旦能量方程被求解，新的焓被用来更新温度与气体属性（如比热）。
除主求解器外，新的总压与总温边界条件被实现用于亚声速入口补丁，这些是模拟可压缩流动时常用的边界条件。边界条件在目录 derivedFvPatchFields 中定义。对于亚声速入口的总压条件（totalPressureComp），updateCoeffs() 函数被修改为 Listing 16.16 给出的形式，其中根据方程 (16.40) 计算边界静压并存入 newp 变量。对于总入口压力条件，还需要 totalPressureCorrectorComp 来处理入口处的质量流率修正影响边界元素对角系数的效应（方程 (16.43) 与 (16.44)）。其实现如 Listing 16.17 所示，把对系数的贡献重置为零，仅保留 gradientInternalCoeffs() 用于正确返回值，保证对散度算子和 Laplacian 算子都施加一个边界条件。修正后的压力修正方程边界元素对角系数在 Listing 16.18 中给出，其中 deltaM() 函数按方程 (16.44) 原样实现，Laplacian 算子贡献被去除。
对于能量方程在亚声速入口的总温条件（totalTemp），按方程 (16.45) 把总温转化为静温并以 Dirichlet 边界条件形式施加。updateCoeffs() 函数被修改为 Listing 16.19 给出的形式：每次迭代时根据入口总温与边界速度更新温度值。
totalVelocity 函数（Listing 16.20）实现了在亚声速入口应用总条件时对速度场所需的更新。为此使用 Dirichlet 边界条件，速度值根据边界上的质量通量 mDot 迭代求出。算法从通量场 mDot 中取出值（该值在调用 flux() 函数求解压力修正方程之后被更新），再除以面面积与对应密度。引入 clip 变量用于防止入口出现回流。neg 函数对负值返回 1、对正值返回 0，从而不允许外向速度被接受。inletDir 变量给出用户定义的入口速度方向。
16.10 小结（Closure）
本章把前一章所发展的不可压缩分块压力基方法扩展到能够处理全速度域的可压缩流体流动。这需要修改压力修正方程使其包含一个 convection-like 项，从而把方程的类型从椭圆型变为双曲型。还需要修改动量方程、求解能量方程，并加入状态方程。同样关键的是可压缩流动模拟所需的特殊边界条件。本章给出了若干边界条件以及一些实现细节。
对基础不可压缩代码所做的修改相对而言是局部的，但却能带来其求解能力的实质性扩展。下一章将介绍处理时间平均 Navier-Stokes 方程所需的额外技术，这些方程是求解湍流问题所必需的。
本章个人批注
读完本章，我对 Moukalled 等教材把"压力基方法从不可压缩扩展到可压缩"这一技术脉络的整体思路有了更清楚的认识。本章的核心信息并不是某个具体公式，而是 SIMPLE 系列方法能"横跨整个 Mach 数谱"这一事实背后的两个机制：(1) 压力在低 Mach 数时通过梯度作用于速度、在高 Mach 数时通过状态方程作用于密度；(2) 压力修正方程因为新增了 convection-like 项，方程类型从椭圆变为双曲。
这两条机制合在一起就解释了为什么 SIMPLE 可以在不作人为调整的情况下同时处理不可压缩流与高超声速流。第 16.5 节关于压力修正方程为何必须严格指定边界压力值（而不是像不可压缩那样任意指定）的论述（第 380–382 行）是我之前没有清晰理解过的：椭圆型方程的常数解 \(p' + C\) 总是满足方程，因此不可压缩情形下边界压力的相对值不影响解；而可压缩情形下由于方程变为双曲，\(p' + C\) 不再是解，边界压力必须显式给定。
关于 uFVM 与 OpenFOAM® 实现：本章在 16.9.1 与 16.9.2 中给出了清晰的"代码 ↔ 方程"对照表，特别是 Listing 16.1 (cfdAssembleMdotTerm) 与 Listing 16.9 (pressure correction equation assembly) 都对应 (16.19) 中的 \(a^0_{F,p}\) 与 \(a^0_{C,p}\)。这种"方程式 ↔ 代码行"的并排展示对实现层面的学习非常有帮助。但同时也应注意，Listing 16.17 中显式把 valueInternalCoeffs/gradientBoundaryCoeffs/valueBoundaryCoeffs 都返回 0，仅保留 gradientInternalCoeffs 的非平凡表达——这是 OpenFOAM® 实现"修正边界元素对角系数"的典型做法，对照 (16.43)/(16.44) 看代码逻辑会更顺畅。
未覆盖的源文习题（16.11 共 7 道）按 skill 规则跳过，不在内概述中展开。其中 Exercise 1 关于气体供应管网的迭代求解与 Exercise 2/3 关于收扩喷管中压力修正方程的建立，是对 16.5 节理论的延伸应用。
与上下章的衔接（一段话）
第 15 章为不可压缩情形下的 SIMPLE/SIMPLEC/PISO/PRIME 体系建立了完整的理论推导与 OpenFOAM® 实现，并把 Rhie-Chow 插值作为在同位网格上避免棋盘问题的关键工具。本章作为第 15 章的直接续篇，把同样的压力基框架扩展到可压缩流——保留 SIMPLE 的整体结构（动量方程 → 压力修正方程 → 场修正 → 能量方程 → 收敛判断），但在压力修正方程中加入密度修正所引起的 convection-like 项、把方程从椭圆型转化为双曲型，从而能够处理跨声速、超声速与高超声速流动中的激波；并在动量方程中加入体粘度项 \(- (2/3)\nabla(\mu \nabla \cdot \mathbf{v})\)、在能量方程中显式加入比热、压力物质导数、体粘度耗散、粘性耗散与热源/汇等额外项；边界条件也针对 Mach 数区分亚声速/超声速入口与出口、区分指定速度/静压/总压等多种类型。作者用这一章把"压力基 SIMPLE 体系能通吃全速度域"的论点完整地展现出来，同时也清楚地指出了下一章的方向：把同一框架叠加到时间平均 Navier-Stokes 方程上以处理湍流问题。