9.5 低马赫数预处理
9.5 低马赫数预处理
在低亚声速马赫数范围，当流速远小于声速时，控制方程（2.19）的对流项出现刚性。以下示例可以说明这一点。三维情况下，五个特征值为
$$ \begin{aligned}(\Lambda_c)_{1,2,3}&=V\quad\text{（对流模态）},\\(\Lambda_c)_{4,5}&=V\pm c\quad\text{（声学模态）}.\end{aligned}\tag{9.39} $$
其中，V 为逆变速度，c 为声速。时间推进时控制方程的刚性由特征条件数决定，其定义为最大与最小特征值之比：
$$ C_{\mathrm{N}}=\frac{\left|(\Lambda_{c})_{\max}\right|}{\left|(\Lambda_{c})_{\min}\right|}=\frac{\left|V\right|+c}{\left|V\right|}=\frac{M+1}{M}.\tag{9.40} $$
允许的局部时间步长受到最快波，即 $ (\Lambda_{c})_{4} $ 的限制。在一个时间步内，最慢波仅传播单元宽度的一小部分，即 $ \Lambda_{\operatorname*{min}}\Delta t\approx(\Lambda_{\operatorname*{min}}/\Lambda_{\operatorname*{max}})h=h/C_{\operatorname{N}} $。因此，较大的条件数 $ C_{\operatorname{N}} $，即 $ M\to0 $ 的情况，会降低波传播效率，减慢定常收敛 [90]。此外，文献 [91,92] 表明，马赫数趋于零时，可压缩流动格式的人工耗散没有正确的尺度关系，因此这类空间离散在低马赫数下精度下降 [93]。
若整个流场速度均较低 $ (M < 0.2) $，则可忽略可压缩性，采用不可压缩方程。不可压缩 Navier–Stokes 方程可用著名的基于压力的格式求解 $ [94-99] $，也可采用人工可压缩性（或伪可压缩性）方法 $ [100-106] $。不过，还存在以下流动情况：
-
高速流动中包含较大的低速区域，例如强收敛喷管上游的亚声速流动。
-
流速较低，但热源引起的密度变化使流动具有可压缩性，例如表面传热或体积加热（燃烧模拟）。
-
不同马赫数下的可压缩和不可压缩流动同时存在，即全速域流动。例如推进系统、高升力构型，以及垂直／短距起降飞行器机动中的流动。
这些情况需要采用可压缩控制方程。为在低马赫数下高效、准确地求解，可以采用预处理。预处理的优势是使一种求解方法适用于所有马赫数。下面推导预处理控制方程。
9.5.1 预处理方程的推导
以一维 Euler 方程说明预处理，其微分形式为
$$ \frac{\partial\vec{W}}{\partial t}+\frac{\partial\vec{F_{c}}}{\partial x}=0,\tag{9.41} $$
其中，$ \vec{W} = [\rho, \rho u, \rho E]^{T} $ 为守恒变量矢量，$ \vec{F}_{c} $ 表示对流通量。为考察预处理对谱半径（计算时间步长所需）及对流通量 Jacobian（迎风耗散所需）的影响，将 Euler 方程（9.41）改写为拟线性形式：
$$ \frac{\partial\vec{W}}{\partial t}+\bar{A}_{c}\frac{\partial\vec{W}}{\partial x}=0\tag{9.42} $$
其中，$ \vec{A}_{c}=\partial\vec{F}_{c}/\partial\vec{W} $ 为对流通量 Jacobian（参见 A.2 节）。
低马赫数预处理的思路是变换控制方程（9.41），使新方程在低马赫数下，即低于 $ M \approx 0.2 $ 时，具有更好的性质。首先，希望使对流特征值趋于一致，限制条件数（见式（9.40）），消除 $ M \rightarrow 0 $ 时的刚性；其次，希望改变数值耗散的尺度关系，提高精度；最后，希望像基于压力的格式一样耦合压力和速度。该变换用另一组流动变量替代守恒变量 $ \vec{W} $。选择有多种 [107]，最常用的是压力、速度分量和温度。将式（9.42）的拟线性形式变换为新变量 $ \vec{W}_{p} $，得到
$$ \bar{P}\frac{\partial\vec{W}_{p}}{\partial t}+\bar{A}_{c}\bar{P}\frac{\partial\vec{W}_{p}}{\partial x}=0.\tag{9.43} $$
其中，$ \vec{P} = \partial\vec{W}/\partial\vec{W}_{p} $ 表示从新变量 $ \vec{W}_{p} $ 到守恒变量 $ \vec{W} $ 的变换矩阵。引入新的通量 Jacobian：
$$ \bar{A}_{c,p}=\bar{A}_{c}\bar{P}=\frac{\partial\vec{F}_{c}}{\partial\vec{W}_{p}},\tag{9.44} $$
即可将式（9.43）写为
$$ \bar{P}\frac{\partial\vec{W}_{p}}{\partial t}+\bar{A}_{c,p}\frac{\partial\vec{W}_{p}}{\partial x}=0.\tag{9.45} $$
若将时间导数前的 $ \bar{P} $ 替换为适当的预处理矩阵 $ \bar{\Gamma} $（下文给出），预处理后的式（9.45）为
$$ \bar{\Gamma}\frac{\partial\vec{W}_{p}}{\partial t}+\bar{A}_{c,p}\frac{\partial\vec{W}_{p}}{\partial x}=0,\tag{9.46} $$
等价地写为
$$ \frac{\partial\vec{W}_{p}}{\partial t}+\bar{\Gamma}^{-1}\bar{A}_{c,p}\frac{\partial\vec{W}_{p}}{\partial x}=0.\tag{9.47} $$
现在，可以采用任意标准空间及时间离散格式，在新变量 $ \vec{W}_{p} $ 中求解式（9.47）。另一种方法是将预处理方程（9.47）变换回守恒变量 $ \vec{W} $，得到
$$ \bar{P}^{-1}\frac{\partial\vec{W}}{\partial t}+\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1}\frac{\partial\vec{W}}{\partial x}=0,\tag{9.48} $$
其中，$ \bar{P}^{-1} = \partial \vec{W}_{p} / \partial \vec{W} $ 是从守恒变量出发的变换矩阵。为方便起见，式（9.48）可写为
$$ \frac{\partial\vec{W}}{\partial t}+\bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1}\frac{\partial\vec{W}}{\partial x}=0,\tag{9.49} $$
或
$$ \frac{\partial\vec{W}}{\partial t}+\bar{P}\bar{\Gamma}^{-1}\bar{A}_{c}\frac{\partial\vec{W}}{\partial x}=0.\tag{9.50} $$
式（9.49）、（9.50）中的项 $ \bar{P}\bar{\Gamma}^{-1} $ 称为守恒变量预处理矩阵。
现在可见，守恒变量形式的预处理方程 $ (9.49) $ 和 $ (9.50) $ 具有以下特点：
-
空间导数乘以 $ \bar{P}\bar{\Gamma}^{-1} $。
-
非定常方程不同于式（9.42）的原始形式，因此解不再具有时间精度。
-
定常解，即 $ \partial\tilde{W}/\partial t=0 $，保持不变。
-
预处理系统的特征值和特征向量为矩阵 $ \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c}=\bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1} $ 的特征值及特征向量，因此不同于式（9.41）的原系统。
主要任务是寻找适当的变量组 $ \vec{W}_{p} $ 及矩阵 $ \bar{\Gamma} $，使预处理系统（9.50）的对流特征值尽可能接近。同时，矩阵在 $ M \rightarrow 0 $ 时也必须有定义，并且希望高马赫数下预处理系统恢复为原始系统。
在 9.5.3 节介绍变换矩阵及预处理矩阵之前，先讨论低马赫数预处理在流动求解器中的实现。
9.5.2 实现
由于采用积分形式求解控制方程，必须将式（9.50）写为
$$ \frac{\partial}{\partial t}\int_{\Omega}\vec{W}\mathrm{d}\Omega+\bar{P}\bar{\Gamma}^{-1}\oint_{\partial\Omega}(\vec{F}_{\mathrm{c}}-\vec{F}_{\mathrm{v}})\mathrm{d}S=\bar{P}\bar{\Gamma}^{-1}\int_{\Omega}\vec{Q}\mathrm{d}\Omega.\tag{9.51} $$
需要注意，式（9.51）对非定常流动不具有守恒性，因此必须采用双时间推进（见 6.3 节），才能得到具有时间精度的解。
采用显式多级格式（6.1 节）求解预处理控制方程（9.51）时，按以下步骤进行：
-
根据矩阵 $ \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c} $ 的谱半径计算新的时间步长。
-
根据预处理通量 Jacobian 的谱半径（如 JST 格式），或特征值与特征向量（如 Roe 迎风格式），计算人工耗散。耗散既可采用守恒变量，也可采用原始变量 $ [107] $ 表示。
-
计算对流通量，既可保持原计算方式，也可采用新变量 $ \vec{W}_{p} $。
-
累加耗散通量和对流通量。根据人工耗散的形式，将整个残差或仅将对流项乘以守恒变量预处理矩阵，即 $ \bar{P}\bar{\Gamma}^{-1} $。
-
将残差乘以 $ \alpha_{k}\Delta t/V $。
-
执行 IRS（可选）。
-
从旧守恒变量 $ \vec{W}^{(0)} $ 中减去残差，得到新守恒变量 $ \vec{W}^{n+1} $。
-
更新边界条件。所有基于特征变量的边界条件，如入口、出口和远场，都需要根据 $ \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c} $ 修改。
隐式格式采用相似过程，但省去步骤 5、6，并以不同方式更新守恒变量（参见 6.2 节）。还必须在隐式算子中计及预处理，因为 Jacobian 为 $ \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c} $。
预处理标量耗散格式（JST）为（见式（4.50））
$$ \begin{aligned}\vec{D}_{I+1/2}=&\ \hat{\Lambda}_{I+1/2}^{\mathrm{S}}(\bar{\Gamma}\bar{P}^{-1})_{I+1/2}\left[\epsilon_{I+1/2}^{(2)}(\vec{W}_{I+1}-\vec{W}_{I})\right.\\&\left.-\epsilon_{I+1/2}^{(4)}(\vec{W}_{I+2}-3\vec{W}_{I+1}+3\vec{W}_{I}-\vec{W}_{I-1})\right].\end{aligned}\tag{9.52} $$
必须包含项 $ \bar{\Gamma}\bar{P}^{-1} $，以便后续将整个残差乘以 $ \bar{P}\bar{\Gamma}^{-1} $，因为谱半径 $ \hat{\Lambda}^{S} $ 已含预处理矩阵，因此不同于式（4.53）。还可以将式（9.52）的 $ \bar{P}^{-1} $ 和 $ \bar{W} $ 合并为 $ W_{p} $，从而用原始变量表示预处理标量耗散格式：
$$ \begin{aligned}\vec{D}_{I+1/2}=&\ \hat{\Lambda}_{I+1/2}^{\mathrm{S}}\left[\epsilon_{I+1/2}^{(2)}\left(\vec{W}_{p,I+1}-\vec{W}_{p,I}\right)\right.\\&\left.-\epsilon_{I+1/2}^{(4)}\left(\vec{W}_{p,I+2}-3\vec{W}_{p,I+1}+3\vec{W}_{p,I}-\vec{W}_{p,I-1}\right)\right].\end{aligned}\tag{9.53} $$
必须认识到，此时仅将对流通量乘以 $ \Gamma^{-1} $，随后将完整残差乘以 $ \bar{P} $。式（9.53）的非守恒形式比关系式（9.52）略简单，但在内部流动或含激波的全速域流动中存在精度问题。
预处理 Roe 迎风格式可写为（参见式（4.91））
$$ \begin{aligned}(\vec{F}_{c})_{I+1/2}&=\frac{1}{2}\left[\vec{F}_{c}(\vec{W}_{R})+\vec{F}_{c}(\vec{W}_{\mathrm{L}})\right.\\&\quad-\left.(\bar{\Gamma}\bar{P}^{-1})_{I+1/2}|\bar{A}_{\mathrm{Roe}}|_{I+1/2}(\vec{W}_{\mathrm{R}}-\vec{W}_{\mathrm{L}})\right],\end{aligned}\tag{9.54} $$
其中，$ |\bar{A}_{Roe}| = \bar{T}_{c,p} |\bar{\Lambda}_{c,p}| \bar{T}_{c,p}^{-1} $。左特征向量 $ (\bar{T}_{c,p}^{-1}) $、右特征向量 $ (\bar{T}_{c,p}) $ 及特征值 $ (\bar{\Lambda}_{c,p}) $ 均对应于矩阵 $ \bar{P} \bar{\Gamma}^{-1} \bar{A}_{c} = \bar{P} \bar{\Gamma}^{-1} \bar{A}_{c,p} \bar{P}^{-1} $。$ (I + 1/2) $ 处的流动变量由式（4.89）的 Roe 平均得到，整个残差乘以 $ \bar{P} \bar{\Gamma}^{-1} $。式（9.54）的预处理 Roe 格式也可以用原始变量 $ \bar{W}_{p} $ 表示，从而得到
$$ \begin{aligned}(\vec{F}_{c})_{I+1/2}&=\frac{1}{2}\left[\vec{F}_{c}(\vec{W}_{p,\mathrm{R}})+\vec{F}_{c}(\vec{W}_{p,\mathrm{L}})\right.\\&\quad-\left.\bar{\Gamma}_{I+1/2}|\bar{A}_{\mathrm{Roe},p}|_{I+1/2}(\vec{W}_{p,\mathrm{R}}-\vec{W}_{p,\mathrm{L}})\right].\end{aligned}\tag{9.55} $$
构成 $ \tilde{A}_{Roe,p} $ 的特征值与特征向量现在由矩阵 $ \Gamma^{-1}\tilde{A}_{c,p} $ 确定（参见式（9.47））。$ \Gamma^{-1}\tilde{A}_{c,p} $ 与 $ \tilde{P}\tilde{\Gamma}^{-1}\tilde{A}_{c} $ 的特征值相同，即为 $ \tilde{\Lambda}_{c,p} $，但特征向量不同。因此，式（9.55）的 Roe 矩阵构造为 $ |\tilde{A}_{Roe,p}| = \tilde{T}_{p} |\tilde{\Lambda}_{c,p}| \tilde{T}_{p}^{-1} $。此时也必须将残差乘以 $ \tilde{P}\tilde{\Gamma}^{-1} $，以变换到守恒变量。低马赫数迎风格式的讨论亦见文献 [108,109]。
9.5.3 矩阵的形式
原始变量 $ \vec{W}_{p} $ 的多种选择中，最常见的形式为
$$ \vec{W}_{\mathrm{p}}=\left[p,\;u,\; v,\;\boldsymbol{w},\;T\right]^{\mathrm{T}}\tag{9.56} $$
因此，后续仅讨论 $ \vec{W}_{p} $ 的这一具体形式。下面给出一般流体及完全气体的变换矩阵、预处理矩阵，以及特征值、左特征向量和右特征向量。
变换矩阵
对于一般流体，从守恒变量到原始变量 $ \bar{P}^{-1} = \partial \vec{W}_{p} / \partial \vec{W} $ 的变换矩阵为 [110]
$$ \bar{P}^{-1}=\left[\begin{array}{ccccc}{{{\frac{\rho h_{\mathrm{T}}+\rho_{\mathrm{T}}(H-q^{2})}{a_{1}}}}}&{{{\frac{\rho_{\mathrm{T}}u}{a_{1}}}}}&{{{\frac{\rho_{\mathrm{T}} v}{a_{1}}}}}&{{{\frac{\rho_{\mathrm{T}}w}{a_{1}}}}}&{{{-\frac{\rho_{\mathrm{T}}}{a_{1}}}}} \\{{{-\frac{u}{\rho}}}}&{{{\frac{1}{\rho}}}}&{{{0}}}&{{{0}}}&{{{0}}} \\{{{-\frac{ v}{\rho}}}}&{{{0}}}&{{{\frac{1}{\rho}}}}&{{{0}}}&{{{0}}} \\{{{-\frac{w}{\rho}}}}&{{{0}}}&{{{0}}}&{{{\frac{1}{\rho}}}}&{{{0}}} \\{{{\frac{1-\rho_{p}(H-q^{2})-\rho h_{p}}{a_{1}}}}}&{{{-\frac{\rho_{p}u}{a_{1}}}}}&{{{-\frac{\rho_{p} v}{a_{1}}}}}&{{{-\frac{\rho_{p}w}{a_{1}}}}}&{{{\frac{\rho_{p}}{a_{1}}}}}\end{array}\right]\tag{9.57} $$
其中，
$$ \begin{array}{l}{q^{2}\;=\;||\vec{ v}||_{2}^{2}=u^{2}+ v^{2}+w^{2},}\\ {a_{1}\;=\;\rho\rho_{p}h_{\mathrm{T}}+\rho_{\mathrm{T}}(1-\rho h_{p}).}\end{array}\tag{9.58} $$
从原始变量到守恒变量的变换矩阵，即 $ \bar{P} = \partial \vec{W} / \partial \vec{W}_{p} $，为 [110]
$$ \bar{P}=\left[\begin{array}{c c c c c}{\rho_{p}}&{0}&{0}&{0}&{\rho_{\mathrm{T}}}\\ {\rho_{p}u}&{\rho}&{0}&{0}&{\rho_{\mathrm{T}}u}\\ {\rho_{p} v}&{0}&{\rho}&{0}&{\rho_{\mathrm{T}} v}\\ {\rho_{p}w}&{0}&{0}&{\rho}&{\rho_{\mathrm{T}}w}\\ {\rho_{p}H-1-\rho h_{p}}&{\rho u}&{\rho v}&{\rho w}&{\rho_{\mathrm{T}}H+\rho h_{\mathrm{T}}}\end{array}\right].\tag{9.59} $$
式（9.57）—（9.59）中密度及焓对压力和温度的导数可写为
$$ \begin{aligned}\rho_{p}&=\rho\alpha_{p},\\\rho_{T}&=-\rho\alpha_{T},\\h_{p}&=\frac{1-\alpha_{T}T}{\rho},\\h_{T}&=c_{p},\end{aligned}\tag{9.60} $$
其中，$ \alpha_{p} $ 和 $ \alpha_{T} $ 分别为定压与定温压缩性系数。声速可由下式计算：
$$ c^{2}=\frac{\rho h_{T}}{a_{1}}.\tag{9.61} $$
对于完全气体（见 2.4.1 节），式（9.60）的压缩性系数变为 $ \alpha_{p}=1/p $ 和 $ \alpha_{T}=1/T $。此时，式（9.57）中从守恒变量到原始变量 $ \bar{P}^{-1}=\partial\vec{W}_{p}/\partial\vec{W} $ 的变换矩阵可写为
$$ \bar{P}^{-1}=\left[\begin{array}{ccccc}{{{(\gamma-1)\frac{q^{2}}{2}}}}&{{{(1-\gamma)u}}}&{{{(1-\gamma) v}}}&{{{(1-\gamma)w}}}&{{{\gamma-1}}} \\{{{-\frac{u}{\rho}}}}&{{{\frac{1}{\rho}}}}&{{{0}}}&{{{0}}}&{{{0}}} \\{{{-\frac{ v}{\rho}}}}&{{{0}}}&{{{\frac{1}{\rho}}}}&{{{0}}}&{{{0}}} \\{{{-\frac{w}{\rho}}}}&{{{0}}}&{{{0}}}&{{{\frac{1}{\rho}}}}&{{{0}}} \\{{{\frac{1}{\rho}\left[\frac{\gamma q^{2}}{2c_{p}}-T\right]}}}&{{{-\frac{\gamma u}{c_{p}\rho}}}}&{{{-\frac{\gamma v}{c_{p}\rho}}}}&{{{-\frac{\gamma w}{c_{p}\rho}}}}&{{{\frac{\gamma}{c_{p}\rho}}}}\end{array}\right].\tag{9.62} $$
对于完全气体，从原始变量到守恒变量的变换矩阵，即 $ \bar{P} = \partial \vec{W} / \partial \vec{W}_{p} $，为
$$ \bar{P}=\left[\begin{array}{ccccc}\frac{\rho}{p}&0&0&0&-\frac{\rho}{T}\\\frac{\rho u}{p}&\rho&0&0&-\frac{\rho u}{T}\\\frac{\rho v}{p}&0&\rho&0&-\frac{\rho v}{T}\\\frac{\rho w}{p}&0&0&\rho&-\frac{\rho w}{T}\\\frac{\rho E}{p}&\rho u&\rho v&\rho w&-\frac{\rho q^{2}}{2T}\end{array}\right]\tag{9.63} $$
其中，$ q^{2} $ 的定义见式（9.58）。
预处理矩阵
对于 Euler 方程，构造预处理矩阵相对容易。Van Leer 等人 $ [111] $ 提出的形式能够达到最低的条件数。文献 $ [91] $ 对该方法作了详细讨论。然而，Navier–Stokes 方程的预处理更为复杂。这是因为黏性项会产生复数波速，使预处理系统难以分析。最广为认可的黏性流动预处理方法
分别由 Choi 和 Merkle [112,113]、Turkel [107,114–116]、Lee 和 van Leer [117,118] 及 Lee [92]、Jorgenson 和 Pletcher [119]，以及 Weiss 和 Smith [48,120] 提出。应用实例见文献 [110,121–135]。
Weiss–Smith 预处理器
对于一般流体，Weiss 和 Smith [48,120] 的预处理矩阵 $ \bar{\Gamma} $ 与式（9.59）中的 $ \bar{P} $ 具有相同形式，只需将 $ \rho_{p} $ 替换为适当的预处理参数 $ \theta $。预处理矩阵的逆矩阵亦如此，即 $ \bar{\Gamma}^{-1} $ 与式（9.57）的 $ \bar{P}^{-1} $ 形式相似。因此，式（9.58）的参数 $ a_{1} $ 改为
$$ a_{1}^{\Gamma}=\rho\theta h_{\mathrm{T}}+\rho_{\mathrm{T}}(1-\rho h_{p}).\tag{9.64} $$
对于完全气体，预处理矩阵 $ \bar{\Gamma} $ 可写为
$$ \bar{\Gamma}=\left[\begin{array}{c c c c c}{\theta}&{0}&{0}&{0}&{-\frac{\rho}{T}}\\ {\theta u}&{\rho}&{0}&{0}&{-\frac{\rho u}{T}}\\ {\theta v}&{0}&{\rho}&{0}&{-\frac{\rho v}{T}}\\ {\theta w}&{0}&{0}&{\rho}&{-\frac{\rho w}{T}}\\ {\theta H-1}&{\rho u}&{\rho v}&{\rho w}&{-\frac{\rho q^{2}}{2T}}\\ \end{array}\right],\tag{9.65} $$
其中，$ \theta $ 仍为预处理参数，其定义将在后文给出。预处理矩阵的逆矩阵为
$$ \bar{\Gamma}^{-1}=\begin{bmatrix}a_{2}\left[\frac{c^{2}}{\gamma-1}-\left(H-q^{2}\right)\right]&&-a_{2}u&-a_{2} v&-a_{2}w& a_{2}\\ -\frac{u}{\rho}&&\frac{1}{\rho}&0&0&0\\ -\frac{ v}{\rho}&&0&\frac{1}{\rho}&0&0\\ -\frac{w}{\rho}&&0&0&\frac{1}{\rho}&0\\ a_{3}\left[1-\theta(H-q^{2})\right]&&-\theta a_{3}u&-\theta a_{3} v&-\theta a_{3}w&\theta a_{3}\end{bmatrix}\tag{9.66} $$
其中采用以下简记：
$$ \begin{gathered}a_{2}=(\gamma-1)\phi,\\a_{3}=\frac{(\gamma-1)\phi T}{\rho},\\\phi=\frac{1}{\theta c^{2}-(\gamma-1)}.\end{gathered}\tag{9.67} $$
式（9.64）—（9.67）的预处理参数 $ \theta $ 有多种选择。文献 [48] 给出一种基于参考速度 $ u_{r} $ 的定义。对于一般流体，其表达式为
$$ \theta=\frac{1}{u_{\mathrm{r}}^{2}}-\frac{\rho_{\mathrm{T}}}{\rho h_{\mathrm{T}}},\tag{9.68} $$
对于完全气体，相应表达式为
$$ \theta=\frac{1}{u_{\mathrm{r}}^{2}}+(\gamma-1)\frac{1}{c^{2}}.\tag{9.69} $$
式（9.68）或（9.69）的参考速度 $ u_{r} $ 由以下关系确定：
$$ u_{\mathrm{r}}=\min\left[\max\left(||\vec{ v}||_{2},\frac{\nu}{\Delta h},\frac{\kappa}{\Delta h},\epsilon\sqrt{\frac{|\Delta p|}{\rho}}\right),c\right],\tag{9.70} $$
其中，$ \Delta h $ 表征控制体的尺寸，$ \Delta p $ 表示相邻控制体之间的压差，$ \epsilon $ 为一个小数（$ \approx 10^{-3} $）。由式（9.70）可见，参考速度受到当地输运速度的约束。当边界层内扩散或热传导占主导时，$ \nu/\Delta h $ 和 $ \kappa/\Delta h $ 两项变得重要。压力项用于防止 $ u_{r} $ 在驻点处趋于零。当 $ u_{r} = c $ 时，$ \bar{\Gamma} $ 与 $ \bar{P} $ 相同，$ \bar{\Gamma}^{-1} $ 则变为 $ \bar{P}^{-1} $。因此，对超声速流动，预处理会按预期关闭。
另一种针对完全气体提出的选择是令 $ [123] $：
$$ \theta=\frac{1}{\beta\gamma R T}=\frac{1}{\beta c^{2}}.\tag{9.71} $$
式（9.71）的参数 $ \beta $ 定义为
$$ \beta=\frac{M_{\mathrm{r}}^{2}}{1+(\gamma-1)M_{\mathrm{r}}^{2}}\tag{9.72} $$
其中，参考马赫数为
$$ M_{\mathrm{r}}^{2}=\max\left[\min(M^{2},1),M_{\min}^{2}\right],\tag{9.73} $$
M 为当地马赫数 $ (M^{2}=||\vec{v}||_{2}^{2}/c^{2}) $。根据文献 $ ^{[123]} $，参数取 $ M_{min}^{2}=KM_{\infty}^{2} $ 和 $ K\approx3 $（也有人取 K=1，甚至 K=0.15；确切取值似乎取决于驻点区域或边界层内的控制体数量）。容易验证，当 $ M\geq1 $ 时，参数 $ \beta $ 等于 $ 1/\gamma $，矩阵 $ \bar{\Gamma} $ 与 $ \bar{P} $ 相同。需要说明，式（9.69）与式（9.71）对 $ \theta $ 的两种定义是等价的。因此，也可借助式（9.70），由 $ M_{r}=u_{r}/c $ 确定参考马赫数。
预处理系统的特征值
预处理系统式（9.49）的特征值矩阵，即 $ \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c}=\bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1} $，为
$$ \bar{\Lambda}_{c,p}=\left[\begin{array}{ccccc}{{{V}}}&{{{0}}}&{{{0}}}&{{{0}}}&{{{0}}} \\{{{0}}}&{{{V}}}&{{{0}}}&{{{0}}}&{{{0}}} \\{{{0}}}&{{{0}}}&{{{V}}}&{{{0}}}&{{{0}}} \\{{{0}}}&{{{0}}}&{{{0}}}&{{{\frac{(a_{4}+1)}{2}V+c^{^{\prime}}}}}&{{{0}}} \\{{{0}}}&{{{0}}}&{{{0}}}&{{{0}}}&{{{\frac{(a_{4}+1)}{2}V-c^{^{\prime}}}}}\end{array}\right],\tag{9.74} $$
其中，$ V = \vec{ v} \cdot \vec{n} $ 表示逆变速度，而
$$ c^{\prime}=\frac{1}{2}\sqrt{V^{2}(a_{4}-1)^{2}+4a_{5}}\tag{9.75} $$
表示修正后的声速。对于一般情形，式（9.74）和（9.75）的参数 $ a_{4} $ 与 $ a_{5} $ 为
$$ \begin{array}{l}a_{4}=\frac{a_{1}}{a_{1}^{\Gamma}},\\a_{5}=\frac{\rho h_{\mathrm{T}}}{a_{1}^{\Gamma}},\end{array}\tag{9.76} $$
其中，$ a_{1} $ 定义于式（9.58），$ a_{1}^{\Gamma} $ 定义于式（9.64）。对于完全气体，若按式（9.69）定义预处理参数 $ \theta $，则上述参数变为 $ a_{4} = \phi $ 和 $ a_{5} = \phi c^{2} $，其中 $ \phi $ 由式（9.67）给出。若采用
式（9.71）对 $ \theta $ 的第二种定义，则式（9.76）的两个参数分别简化为 $ a_{4}=M_{r}^{2} $ 和 $ a_{5}=M_{r}^{2}c^{2} $。可见，$ c^{\prime}\approx(V/2)\sqrt{5} $ 在 $ |M|\rightarrow0 $ 时成立，因此各特征值按预期趋于一致。这样便降低了式（9.40）的条件数所表征的刚性（条件数为 $ C_{N}\approx2.6 $），显著加快时间推进或迭代求解过程的收敛。另一方面，$ a_{4}=1 $ 和 $ c^{\prime}=c $ 在 $ |M|\geq1 $ 时成立，因而恢复 $ \bar{A}_{c} $ 的特征值。
根据式（9.74），预处理系统的谱半径为
$$ \hat{\Lambda}_{c,p}=\left[\frac{(a_{4}+1)}{2}|V|+c^{\prime}\right]\Delta S\tag{9.77} $$
其中，$ \Delta S $ 表示面的面积。该表达式用于计算式（6.14）、（6.18）、（6.20）或（6.22）的时间步长。对于中心格式的人工耗散，它还用于替换式（4.53）的谱半径，另见式（5.33）、（9.52）和（9.53）。
预处理系统的特征向量
预处理会改变对流通量 Jacobian 矩阵的左、右特征向量（参见 A.11 节）。这一点对于基于特征变量的空间离散尤为重要，例如 Roe 迎风格式（4.3.3 节），以及 4.3.4 节介绍的迎风 TVD 格式。
式（9.47）中，预处理后的新通量 Jacobian 矩阵 $ \Gamma^{-1}\bar{A}_{c,p} $ 的左特征向量矩阵可写为 [110]
$$ \bar{T}_{p}^{-1}=\begin{bmatrix}-\frac{\rho_{\mathrm{T}}n_{x}}{a_{1}^{\Gamma}}&0&\frac{a_{6}n_{z}}{c^{\prime}}&-\frac{a_{6}n_y}{c^{\prime}}&-\frac{a_{6}n_{x}}{T}\\ -\frac{\rho_{\mathrm{T}}n_y}{a_{1}^{\Gamma}}&-\frac{a_{6}n_{z}}{c^{\prime}}&0&\frac{a_{6}n_{x}}{c^{\prime}}&-\frac{a_{6}n_y}{T}\\ -\frac{\rho_{\mathrm{T}}n_{z}}{a_{1}^{\Gamma}}&\frac{a_{6}n_y}{c^{\prime}}&-\frac{a_{6}n_{x}}{c^{\prime}}&0&-\frac{a_{6}n_{z}}{T}\\ \frac{1}{2}+a_{7}&\frac{a_{6}n_{x}}{2c^{\prime}}&\frac{a_{6}n_y}{2c^{\prime}}&\frac{a_{6}n_{z}}{2c^{\prime}}&0\\ \frac{1}{2}-a_{7}&-\frac{a_{6}n_{x}}{2c^{\prime}}&-\frac{a_{6}n_y}{2c^{\prime}}&-\frac{a_{6}n_{z}}{2c^{\prime}}&0\end{bmatrix}.\tag{9.78} $$
按照式（9.55），将 $ \Gamma^{-1}\bar{A}_{c,p} $ 的右特征向量矩阵乘以 $ \bar{\Gamma} $，可得 [110]
$$ \bar\Gamma\bar T_p=\frac{a_1^{\Gamma}}{\rho h_T}\begin{bmatrix}-a_{8}n_{x}&-a_{8}n_y&-a_{8}n_{z}&1&1\\-a_{8}u n_{x}&-a_{8}u n_y-c^{\prime}n_{z}&-a_{8}u n_{z}+c^{\prime}n_y&u+\left(a_{9}+c^{\prime}\right)n_{x}&u+\left(a_{9}-c^{\prime}\right)n_{x}\\-a_{8} v n_{x}+c^{\prime}n_{z}&-a_{8} v n_y&-a_{8} v n_{z}-c^{\prime}n_{x}& v+\left(a_{9}+c^{\prime}\right)n_y& v+\left(a_{9}-c^{\prime}\right)n_y\\-a_{8}w n_{x}-c^{\prime}n_y&-a_{8}w n_y+c^{\prime}n_{x}&-a_{8}w n_{z}&w+\left(a_{9}+c^{\prime}\right)n_{z}&w+\left(a_{9}-c^{\prime}\right)n_{z}\\a_{11}-a_{10}n_{x}&a_{12}-a_{10}n_y&a_{13}-a_{10}n_{z}&H+\left(a_{9}+c^{\prime}\right)V&H+\left(a_{9}-c^{\prime}\right)V\end{bmatrix}.\tag{9.79} $$
其中，$ a_{1}^{\Gamma} $ 定义于式（9.64）。式（9.78）和（9.79）的其余参数为
$$ \begin{aligned}a_{6}&=\rho a_{5}\\a_{7}&=\frac{V(a_{4}-1)}{4c^{\prime}}\\a_{8}&=\frac{\rho_{\mathrm{T}}T}{\rho}\\a_{9}&=-\frac{1}{2}V(a_{4}-1)\\a_{10}&=a_{8}H+T h_{\mathrm{T}}\\a_{11}&=c^{\prime}(v n_{z}-w n_y)\\a_{12}&=c^{\prime}(w n_{x}-u n_{z})\\a_{13}&=c^{\prime}(u n_y-v n_{x})\end{aligned}\tag{9.80} $$
其中，$ a_{4} $ 和 $ a_{5} $ 定义于式（9.76），$ c^{\prime} $ 定义于式（9.75）。
特征向量 $ \bar{T}_{c,p} $ 和 $ \bar{T}_{c,p}^{-1} $ 对应于矩阵 $ \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c}=\bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1} $；采用守恒变量形式的预处理 Roe 格式（见式（9.54））需要这些特征向量，它们可通过修正式（9.78）和（9.79）得到。可以证明，$ \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c} $ 的右特征向量矩阵的各列
由 $ \vec{P}\vec{x} $ 构成，其中，$ \vec{x} $ 是 $ \Gamma^{-1}\vec{A}_{c,p} $ 的右特征向量。因此，$ \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c} $ 的右、左特征向量矩阵可表示为
$$ \begin{aligned}\bar{T}_{c,p}&=\bar{P}\bar{T}_{p},\\\bar{T}_{c,p}^{-1}&=\bar{T}_{p}^{-1}\bar{P}^{-1}.\end{aligned}\tag{9.81} $$
将式（9.81）的上述关系代入守恒变量形式的预处理 Roe 格式（9.54），得到
$$ \begin{aligned}(\vec{F}_{c})_{I+1/2}&=\frac{1}{2}\left[\vec{F}_{c}\big(\vec{W}_{\mathrm{R}}\big)+\vec{F}_{c}\big(\vec{W}_{\mathrm{L}}\big)\right.\\&\quad\left.-(\bar{\Gamma}\ \bar{T}_{p}\mid\bar{\Lambda}_{c,p}\mid\bar{T}_{p}^{-1}\bar{P}^{-1})_{I+1/2}\left(\vec{W}_{\mathrm{R}}-\vec{W}_{\mathrm{L}}\right)\right]\end{aligned}\tag{9.82} $$
其中，$ \bar{\Gamma}\bar{T}_{p} $ 由式（9.79）给出，$ \bar{T}_{p}^{-1} $ 由式（9.78）给出。
图 9.11—9.13 给出一个翼型绕流算例：采用式（9.51）的预处理有限体积格式，其中 $ \Gamma $ 按式（9.65）定义。可见，当来流马赫数为 0.01 时，中心格式和 Roe 迎风格式均无法给出正确解。未经预处理的格式不能正确预测压力分布，因而也无法正确预测升力系数（$ C_{L} = 0.323 $ 和 0.324，而正确值为 0.352）。如图 9.13 所示，预处理既能得到正确解（$ C_{L} = 0.353 $），又能显著加快收敛。
需要指出，随附源代码中的两个二维求解器，均对中心格式和 Roe 迎风格式实现了低马赫数预处理。