CFD: 可压 + 稳态
可压缩求解器相对于不可压缩求解器还可以继续进行区分：
是否需要考虑体积力（比如温度变化引起的浮力）；
弱可压缩还是强可压缩（前者典型的为人居环境的温度变化，后者主要是超音速流动）；
对于超音速流动，也即可压缩性特别强的情况，得益于双曲方程的数学特征，这些方法通常采用瞬态计算且使用密度作为主要变量，然后通过状态方程求解压力，即密度基求解器。但密度基求解器在应用于弱可压缩流的情况下效率低下。一方面的解释是在不可压缩领域，压力与密度的耦合非常弱，另一方面的解释是不可压缩假定下音速趋向于无穷大，导致时间步长过小。
针对弱可压缩领域，fluid模块为一个压力基、可压缩、适用于全流速、也包含体积力的求解器，主要用于求解普适性可压缩流动，也可用于求解低音速/超音速流动。其可以处理瞬态问题，也可以处理稳态问题。本文讨论稳态算法。
控制方程
首先有稳态质量方程、动量方程、以及状态方程（在这里不考虑能量方程，因为其对压力速度耦合求解不起决定作用）：
(1)\[
\nabla\cdot\rho\bfU=0,
\]
(2)\[
\nabla \cdot (\rho\mathbf{U} \mathbf{U})=-\nabla p+\nabla \cdot\tau +\rho\bfg,
\label{mom}
\]
(3)\[
p=\rho RT.
\]
方程(1)没有时间项，并无明显的变量，更像是一种限定性条件。因此稳态求解器中不存在对质量方程的求解。方程(1)中的\(\rho\bfU\)可用于组建质量通量。上述三个方程，可用于迭代求解三个未知量：速度、压力、密度。在讨论体积力之前，首先忽略\(\rho\bfg\)的影响。
离散
首先对方程(2)通过高斯定理进行对速度\(\bfU\)的离散，组建速度方程有：
(4)\[
\int {\nabla \cdot \left( {\rho\mathbf{U}\mathbf{U}} \right)\mathrm{d}V =
\int {
\rho\mathbf{U}\mathbf{U}\cdot\rd\bfS} = } \sum_f {{{\left( {\rho^n{\mathbf{U}^n}{\mathbf{U}^{*}}} \right)}_f}} \cdot\bfS_f = \sum_f F_f^n \mathbf{U}_f^{*} ,
\]
(5)\[
\int \nabla p \mathrm{d}V
=\int p \mathrm{d}\bfS
=\sum_f p_f^n\bfS_f,
\]
其中上标\(^n\)表示为当前迭代步（已知），上标\(^{*}\)表示下一个迭代步（待求），下标\(_f\)表示网格单元面上的值，\(\bfS_f\)表示网格单元的各个面的面矢量，\(F_f\)为质量通量，\(\mu\)为运动粘度（假定粘度为常数）。此处略去粘性项\(\tau\)的离散且假定粘度为\(0\)。需要注意的是，目前并没有处理\(\nabla\cdot\tau\)这一项，因为其对求解流程并不影响。
注意：
OpenFOAM中可压缩求解器的通量phi定义为\(\rho_f\bfU_f\cdot\bfS_f\)，不可压缩求解器的通量phi定义为\(\bfU_f\cdot\bfS_f\)。
整理方程(4)、(5)有：
(6)\[
\sum_f F_f^n \mathbf{U}_f^{*}
= -\sum_f p_f^n\bfS_f.
\]
需要注意的是，方程(6)左侧的对流项采用了线性化处理。方程(6)中的\(\bfU_f\)需要从体心速度进行插值来获得，在此步可以引入各种插值格式。假设使用中心线性格式（均一网格）：
(7)\[
\mathbf{U}_f^{*} = \frac{{\mathbf{U}_\mathrm{P}^{*} + \mathbf{U}_\mathrm{N}^{*}}}{2}.
\]
(8)\[
p_f^n = \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}.
\]
将方程(7)，(8)代入到方程(6)有
(9)\[
\sum_f {F_f^n \frac{{\mathbf{U}_\mathrm{P}^{*} + \mathbf{U}_\mathrm{N}^{*}}}{2}} = -\sum_f \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}\bfS_f,
\]
(10)\[
\sum_f { {{\frac{{F_f^n}}{2}} } } \mathbf{U}_\mathrm{P}^{*} = - \sum_f { { {\frac{{F_f^n}}{2} } \mathbf{U}_\mathrm{N}^{*}} } -\sum \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}\bfS_f.
\]
将上式左右两边同时除以网格单元体积\(V_\rP\)，同时简写为
(11)\[
{A_\mathrm{P}^n}\mathbf{U}_\mathrm{P}^{*}{\rm{ + }}\sum_f {A_\mathrm{N}^n\mathbf{U}_\mathrm{N}^{*}} = -\frac{1}{V_\rP}\sum_f \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}\bfS_f,
\]
其中\(A_\rP^n\)，\(A_\rN^n\)分别表示当前网格点与相邻网格点的离散系数：
(12)\[
A_\mathrm{P}^n= \frac{1}{V_\rP}\sum_f \frac{F_f^n}{2} ,
\]
(13)\[
A_\mathrm{N}^n= \frac{1}{V_\rP}{\frac{{F_f^n}}{2}} ,
\]
可见，\(A_\mathrm{P}^n\)与\(A_\mathrm{N}^n\)为关于\(\rho^n\bfU^n\)的代数式。方程(11)构成一个稀疏线性系统。求解方程(11)即可获得速度\(\mathbf{U}^{*}\)。对方程(11)进行转换有：
(14)\[
\mathbf{U}_\mathrm{P}^{*} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A^n_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^n\bfS_f.
\]
其中\(\bfHbyA\)定义为：
(15)\[
\bfHbyA_\rP^{*}=\frac{-\sum_NA_\rN^n\bfU_\rN^{*}}{A_\rP^n}.
\]
注意：
方程(15)中如果未进行动量预测，则右边取\(\bfU_\rN^n\)。
同时，面上插值有：
(16)\[
\mathbf{U}_f^{*} = \mathbf{HbyA}^{*}_f - \frac{1}{{{A_f^n}}} \left(\frac{1}{V_\rP}\sum_f p_f^n\bfS_f\right)_f,
\]
(17)\[
\bfHbyA_f^{*}=\left(\frac{-\sum_NA_\rN^n\bfU_\rN^{*}}{A_\rP^n}\right)_f.
\]
注意：
在求解速度方程后，求解器求解能量方程，并进行thermo.correct()函数，对\(\psi=1/RT\)进行更新。因此，在下文中，用\(\psi^{corr}\)表示更新后的\(1/RT\)。
方程(14)在收敛的情况下可以写为：
(18)\[
\mathbf{U}_\mathrm{P}^{n+1} = \mathbf{HbyA}^{n+1}_\mathrm{P} - \frac{1}{{{A^{n+1}_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{n+1}\bfS_f.
\]
方程(18)与方程(14)相减有（需将\(A\)线性化，类似\(A^{n+1}=A^n\)）：
(19)\[
\mathbf{U}_\mathrm{P}^{'} = \mathbf{HbyA}^{'}_\mathrm{P}
-\frac{1}{{{A^{n}_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{'}\bfS_f.
\]
方程(19)提供了速度修正量与压力修正量的关系。在这里略去临点的影响，方程(19)可以写为：
(20)\[
\mathbf{U}_\mathrm{P}^{'} = - \frac{1}{{{A^{n}_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{'}\bfS_f.
\]
方程(20)提供了速度修正与压力修正的一对一关系。将其与方程(14)加和有：
(21)\[
\mathbf{U}_\mathrm{P}^{**} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A^n_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{*}\bfS_f.
\]
(22)\[
\mathbf{U}_f^{**} = \mathbf{HbyA}^{*}_f - \frac{1}{{{A^n_f}}} \frac{1}{V_\rP}\left(\sum_f p_f^{*}\bfS_f\right)_f.
\]
现在考虑连续性方程：
(23)\[
\sum_f (\rho_f^*+\rho_f')\mathbf{U}_f^{**} \cdot \bfS_f=0.
\]
其中\(\rho_f'\)表示修正密度，\(\bfU_f'\)表示修正速度。可见可压缩与不可压缩最重要的区别在于连续性方程出现了密度修正。
低速流压力更新
在弱可压缩的情况下，求解完动量方程之后，需要求解能量方程，并更新密度中间量：
(24)\[
\rho^{corr1}=p^t/\psi^{corr1}
\]
方程(23)中的密度修正\(\rho_f^*+\rho_f'\)直接用密度中间量替代：
(25)\[
\sum_f \rho_f^{corr1}\mathbf{U}_f^{**} \cdot \bfS_f=0.
\]
将方程(22)代入到方程(25)中有：
(26)\[
\sum_f \rho_f^{corr1} \mathbf{HbyA}^{*}_f \cdot \bfS_f= \sum_f\frac{1}{{{A^n_f}}} \frac{1}{V_\rP}\left(\sum_f p_f^{*}\bfS_f\right)_f \cdot \bfS_f.
\]
可见，方程(26)可用于求解获得压力\(p^*\)。
高速流transonic压力更新
在强可压缩的情况下，要考虑压力的对流效应，而不是密度中间量。结合状态方程：
(27)\[
\rho_f'=\psi^{corr}p_f'
\]
也即
\[
\rho^{corr}=p^n/\psi^{corr}
\]
方程(23)可以写为
(28)\[
\sum_f \left(\rho_f^n+\psi^{corr}p_f'\right)(\mathbf{U}_f^{*}+\bfU_f') \cdot \bfS_f=0.
\]
展开为
(29)\[
\sum_f \left(\left(\rho_f^n+\psi^{corr}p_f'\right)\bfU_f^{*}
+
\rho_f^n\bfU_f'
\right)\cdot\bfS_f=0.
\]
(30)\[
\sum_f \left(\left(\rho_f^n+\psi^{corr}p_f'\right)\bfU_f^{*}
\right)\cdot\bfS_f
=
-\sum_f \left(
\rho_f^n\bfU_f'
\right)\cdot\bfS_f
\]
注意：
方程(30)中的\(p'\bfU'\)为两个修正量的乘积，被省略。
因此在这里引入略去邻点影响的假定：
(31)\[
\mathbf{U}_\mathrm{P}^{'} = - \frac{1}{V_\rP}\frac{1}{{{A_\mathrm{P}^n}}} \sum_f p_f^{'}\bfS_f.
\]
(32)\[
\bfU_f'=-\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f’\bfS_f\right)_f
\]
将方程(32)、(16)代入到(30)有：
\[\begin{split}
\sum_f \left(\rho_f^n+\psi^{corr}p_f'\right)
\mathbf{HbyA}^{*}_f \cdot\bfS_f
-
\sum_f \left(\rho_f^n+\psi^{corr}p_f'\right)
\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\cdot\bfS_f
\hspace{6em}
\\
=
\sum_f \rho_f^n
\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{'}\bfS_f\right)_f \cdot\bfS_f\hspace{3em}
\end{split}\]
\[\begin{split}
\sum_f \left(\rho_f^n+\psi^{corr}p_f'\right)
\mathbf{HbyA}^{*}_f \cdot\bfS_f
-
\sum_f \left(\psi^{corr}p_f'\right)
\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\cdot\bfS_f
\hspace{6em}
\\
=
\sum_f \rho_f^n
\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{*}\bfS_f\right)_f \cdot\bfS_f
\hspace{2em}
\end{split}\]
\[\begin{split}
\sum_f \rho_f^n
\mathbf{HbyA}^{*}_f \cdot\bfS_f
+
\sum_f \psi^{corr}p_f'
\left(\mathbf{HbyA}^{*}_f-\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\right)\cdot\bfS_f
\hspace{6em}
\\
=
\sum_f \rho_f^n
\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{*}\bfS_f\right)_f \cdot\bfS_f
\hspace{2em}
\end{split}\]
\[\begin{split}
\sum_f \rho_f^n
\mathbf{HbyA}^{*}_f \cdot\bfS_f
+
\sum_f \psi^{corr}(p_f^*-p_f^n)
\left(\mathbf{HbyA}^{*}_f-\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\right)\cdot\bfS_f
\hspace{6em}
\\
=
\sum_f \rho_f^n
\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{*}\bfS_f\right)_f \cdot\bfS_f
\hspace{2em}
\end{split}\]
\[\begin{split}
\sum_f \rho_f^n
\mathbf{HbyA}^{*}_f \cdot\bfS_f
-
\sum_f \psi^{corr}p_f^n
\left(\mathbf{HbyA}^{*}_f-\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\right)\cdot\bfS_f
\hspace{9em}
\\
+
\sum_f \psi^{corr}p_f^*
\left(\mathbf{HbyA}^{*}_f-\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\right)\cdot\bfS_f
=
\sum_f \rho_f^n
\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{*}\bfS_f\right)_f \cdot\bfS_f
\hspace{3em}
\end{split}\]
方程最后的形式为：
(33)\[\begin{split}
\sum_f \left(\rho_f^n-\psi^{corr}p_f^n\right)
\mathbf{HbyA}^{*}_f \cdot\bfS_f
+
\sum_f \psi^{corr}p_f^*\mathbf{HbyA}^{*}_f\cdot\bfS_f
\hspace{17em}
\\
+
\sum_f \psi^{corr}(p_f^n - p_f^*)
\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\cdot\bfS_f
=
\sum_f \rho_f^n
\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{*}\bfS_f\right)_f \cdot\bfS_f
\hspace{3em}
\end{split}\]
方程(33)左边的第三项可以进行省略，演变为：
(34)\[\begin{split}
\sum_f \left(\rho_f^n-\psi^{corr}p_f^n\right)
\mathbf{HbyA}^{*}_f \cdot\bfS_f
+
\sum_f \psi^{corr}p_f^*\mathbf{HbyA}^{*}_f\cdot\bfS_f
\hspace{17em}
\\
=
\sum_f \rho_f^n
\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{*}\bfS_f\right)_f \cdot\bfS_f
\hspace{6em}
\end{split}\]
方程(34)可以用来求解\(p^{*}\)。
至此，低速流与高速流的\(p^*\)更新方法简述完毕。现将方程(14)与方程(31)加和有：
(35)\[
\mathbf{U}_\mathrm{P}^{**} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{V_\rP}\frac{1}{{{A_\mathrm{P}^{n}}}} \sum_f p_f^{*}\bfS_f.
\]
将\(p^{*}\)代入到方程(35)有\(\mathbf{U}_\mathrm{P}^{**}\)。随后，密度通过状态方程进行更新：
(36)\[
\rho^*=\psi^{corr}p^*
\]
总之，可压缩算法中的的迭代过程可以表示为下面几个步骤：
通过方程(14)求解获得\(\bfU^*\)，并组建\(\bfHbyA^*\)；
通过能量方程更新可压缩性\(\psi^{corr}\)；
低速流通过方程(26)求解获得\(p^*\)，高速流通过方程(34)求解获得\(p^*\)；
通过方程(35)求解获得\(\bfU^{**}\)；
通过方程(36)更新获得\(\rho^{*}\)；
回到第一步，继续进行求解；
稳态求解器需要对压力进行场松弛，速度、能量等方程进行矩阵松弛增加稳定性。同样的，fluid模块也可以调用SIMPLEC算法。具体内容请参考不可压稳态SIMPLE算法。
同时需要注意的是，高速流与低速流的压力方程特征发生了变化，对比低速流的方程(26)与高速流的方程(34)，低速流为压力的拉普拉斯方程。高速流的压力方程含有压力的对流项。
体积力
可压缩求解器在用于模拟暖通环境的时候，要附加浮力的影响。在跨音速的情况下也通常忽略体积力。因此在有体积力的时候，主要考虑弱可压缩的情况。当然了，也可以顺着思路进行跨音速体积力的扩展。
在存在体积力的时候，动量方程中的重力项以及压力项可以进行数值处理。定义
(37)\[
p_\mathrm{rgh}=p-\rho \bfg \cdot\bfh
\]
其中\(\bfh\)表示网格单元体心的位置矢量。对方程(37)进行梯度操作：
(38)\[
\nabla p_\mathrm{rgh}=\nabla p-\bfg\cdot \bfh \nabla \rho - \rho \bfg
\]
即为
(39)\[
-\nabla p+\rho \bfg=-\bfg\cdot \bfh \nabla \rho -\nabla p_\mathrm{rgh}
\]
将方程(39)代入到方程(2)中有：
(40)\[
\nabla \cdot (\rho\mathbf{U} \mathbf{U})=-\bfg\cdot \bfh \nabla \rho -\nabla p_\mathrm{rgh}+\nabla\cdot\tau
\]
接下来看压力方程的处理。针对方程(31)，在考虑方程(40)后演变为：
(41)\[
\mathbf{U}_\mathrm{P}^{'} = - \frac{1}{{{A^{n}_\mathrm{P}}}} \frac{1}{V_\rP}
\left(
\sum_f p_{\mathrm{rgh},f}^{'}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
+\sum_f \bfg_f\cdot \bfh_f \rho_f^{'}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right)
\]
方程(41)提供了速度修正与压力修正的一对一关系。但发现，其中还存在密度修正。如果顺着这个思路进行，压力方程将进一步出现未知变量密度。
为了处理这个问题，我们可以将方程(41)中的密度修正忽略掉。这是一种方法，但处理起来比较激进。另外一种方法是在能量方程更新可压缩性\(\psi^{corr}\)后，在其基础上，通过状态方程更新密度中间量：
(42)\[
\rho^{corr}=p^n/\psi^{corr}
\]
方程(42)在前面也已经出现过。在这个情况下，认为\(\rho^n\)与密度修正量\(\rho^{'}\)的关系，满足：
(43)\[
\rho^{corr}=\rho^n+\rho^{'}
\]
而非\(\rho^{n+1}=\rho^n+\rho^{'}\)。这样，由于\(\rho^{corr}\)已经通过方程(42)求出，因此压力方程仅仅存在未知量压力。依据此思想，结合修正方程(41)有：
(44)\[
\mathbf{U}_\mathrm{P}^{**} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A^n_\mathrm{P}}}} \frac{1}{V_\rP}
\left(
\sum_f p_{\mathrm{rgh},f}^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
+\sum_f \bfg_f\cdot \bfh_f \rho_f^{corr}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right).
\]
(45)\[
\mathbf{U}_f^{**} = \mathbf{HbyA}^{*}_f - \frac{1}{{{A^n_f}}} \frac{1}{V_\rP}
\left(
\sum_f p_{\mathrm{rgh},f}^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
+\sum_f \bfg_f\cdot \bfh_f \rho_f^{corr}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right)_f.
\]
(46)\[
\mathbf{U}^{**} =\mathbf{U}^{*}+\bfU^{'}
\]
将方程(45)代入到方程(25)中有压力泊松方程：
(47)\[\begin{split}
\sum_f \rho_f^{corr} \mathbf{HbyA}^{*}_f \cdot \bfS_f
-
\sum_f
\frac{\rho^{corr}_f}{{{A^n_f}}} \frac{1}{V_\rP}
\left(
\sum_f \bfg_f\cdot \bfh_f \rho_f^{corr}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right)_f
\cdot \bfS_f
\hspace{5em}
\\
= \sum_f\frac{\rho^{corr}_f}{{{A^n_f}}} \frac{1}{V_\rP}
\left(
\sum_f p_{\mathrm{rgh},f}^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right)_f
\cdot \bfS_f.
\end{split}\]
方程(47)的连续形式为：
(48)\[
\nabla\cdot\left(\rho^{corr} \left(\mathbf{HbyA} - \frac{1}{A}\bfg\cdot\bfh\nabla\rho \right) \right)
=
\nabla\cdot\left(\rho^{corr}\frac{1}{A} \nabla p_{\mathrm{rgh}}^*\right)
\]
方程(48)可用于求解获得压力\(p^*\)。
同时，在考虑体积力的时候，更新速度的时候要考虑体积力的影响：
(49)\[
\mathbf{U}_\mathrm{P}^{**} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A^n_\mathrm{P}}}} \frac{1}{V_\rP}
\left(
\sum_f p_{\mathrm{rgh},f}^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
+\sum_f \bfg_f\cdot \bfh_f \rho_f^{corr}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right)
\]
关键代码
fluid模块中的速度方程并无特殊。在速度方程更性之后，fluid模块中会求解能量方程，并通过下述代码进行可压缩性的更新：
高速流transonic中的压力方程通过下面的代码组建：
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
+ fvm::div(phid, p)
- fvm::laplacian(rhorAAtUf, p)
==
fvModels.source(psi, p, rho.name())
);
其中很明显具有压力的对流项。低速流中的压力方程通过下面的代码组建：
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
- fvm::laplacian(rhorAAtUf, p)
==
fvModels.source(psi, p, rho.name())
);
其为一个拉普拉斯方程。在求解后，需要对密度进行更新：
若考虑体积力，方程(44)中速度方程的更新采用了reconstruct()函数：
U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
其中HbyA不需要进行解释。phig的代码为：
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
其可以翻译为：
\[
-\left(\frac{\rho^{corr}_f}{A^n_\rP}\right)_f \bfg_f\cdot\bfh_f (\nabla\rho^{corr})_f\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\]
其中有关于面法向梯度fvc::snGrad()的离散，请参考《无痛苦NS方程笔记》。p_rghEqn.flux()中的p_rhoEqn代码如下：
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
);
其中的p_rghEqn.flux()可以翻译为：
\[
\left(\frac{\rho^{corr}_f}{A^n_\rP}\right)_f (\nabla p_{\mathrm{rgh}})_f\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\]
更详细的关于p_rghEqn.flux()的解释，请参考《无痛苦NS方程笔记》。因此对于下述代码：
+rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
可以理解为：
\[
+
\frac{1}{A^n_\rP}
\mathrm{recon}\left(
\left(
-\left(\frac{\rho^{corr}}{A^n_\rP}\right)_f \bfg_f\cdot\bfh_f (\nabla\rho^{corr})_f\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
- \left(\frac{\rho^{corr}}{A^n_\rP}\right)_f (\nabla p_{\mathrm{rgh}})_f \frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right)\frac{(A^n_\rP)_f}{\rho^{corr}_f}
\right)
\]
约去相关项并展开reconstruct()函数：
\[
+
\frac{1}{A^n_\rP}
\frac{\sum \bfn_f \cdot\left(
- \bfg_f\cdot\bfh_f (\nabla\rho^{corr})_f\frac{\bfS_f}{|\bfS_f|}
- (\nabla p_{\mathrm{rgh}})_f \frac{\bfS_f}{|\bfS_f|}
\right) |\bfS_f|}
{
\sum \bfn_f\cdot\bfS_f
}
\]
其连续性形式为：
\[
-\frac{1}{A^n}\bfg\cdot\bfh\nabla\rho^{corr}
-\frac{1}{A^n} \nabla p_{\mathrm{rgh}}
\]
更详细的关于reconstruct()的解释，请参考《无痛苦NS方程笔记》。