CFD: 可压 + 稳态 + 体积力

本文讨论的算法对应OpenFOAM中的fluid模块以及isothermalFluid模块。同样对应OpenFOAM-10之前的rhoSimpleFoam、buoyantFoam（稳态部分）、buoyantSimpleFoam。

这些代码主要用于求解可压缩流体以及浮力驱动流。fluid模块主要基于isothermalFluid模块，仅仅在其基础上添加能量方程。在OpenFOAM模块化之前，有无体积力采用不同的求解器来处理，如buoyantSimpleFoam（有体积力），rhoSimpleFoam（无体积力）。在OpenFOAM模块化后都被处理到fluid模块中来减少代码复用。

控制方程

首先有稳态质量方程、动量方程、以及状态方程。能量方程请参考《无痛苦NS方程笔记》，因为其对压力速度耦合求解不起决定作用，在后续不会对其进行离散。有：

(1)\[ \nabla\cdot\rho\bfU=0, \]

(2)\[ \nabla \cdot (\rho\mathbf{U} \mathbf{U})=-\nabla p+\nabla \cdot\tau +\rho\bfg, \label{mom} \]

(3)\[ \rho=f_{\mathrm{EOS}}(p,T). \]

方程(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}^{*}\)。

Warning

可压缩求解器中的系数\(A\)已经乘以密度。

能量方程

求解完动量方程之后，需要求解能量方程。能量方程的离散与速度方程的离散思路相同，在这里不做介绍。求解能量方程后，进入能量方程中的thermo.correct()函数函数。该方程会更新\(T,C_p, C_v\)等一些列热力学参数。有关如何从能量方程反推\(T\)，请参考《无痛苦NS方程笔记》。

Warning

在OpenFOAM模块化求解器中，能量方程可以不进行求解。这种情况下其构成isothermalFluid模块，也即等温可压求解器。

假定更新后的温度为\(T^*\)（对于isothermalFluid模块\(T\)不变）。对于rhoThermo模型类，有流体的可压缩性\(\psi^{*}\)以及不同的状态方程中的密度：

(14)\[ \psi^{*}=f_{\mathrm{EOS}} (T^*) \]

(15)\[ \varrho^{*}=f_{\mathrm{EOS}}(\psi^{*},p^n) \]

Warning

\(\varrho^{*}\)存储在rhoThermo中，而不是求解器的\(\rho\)。注意在迭代过程中，rho对应的\(\rho\)不同于thermo.rho()的\(\varrho^{*}\)。

对于不同的状态方程，方程(14)，(15)存在不同的形式。但总体原则是通过状态方程来更新可压缩性以及密度。

对于psiThermo模型类，其认为密度与压力存在线性关系，因此不会显性的更新密度。需要用密度的时候，用\(\varrho^{*}=p^n \psi^{*}\)计算即可。

弱可压缩压力方程

低速流动可以认为是可压缩性比较小的。因此压力方程的推导与不可压缩基本一致。对方程(11)进行转换有：

(16)\[ \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\)定义为：

(17)\[ \bfHbyA_\rP^{*}=\frac{-\sum_NA_\rN^n\bfU_\rN^{*}}{A_\rP^n}. \]

注意：

方程(17)中如果未进行动量预测，则右边取\(\bfU_\rN^n\)。

同时，面上插值有：

(18)\[ \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, \]

(19)\[ \bfHbyA_f^{*}=\left(\frac{-\sum_NA_\rN^n\bfU_\rN^{*}}{A_\rP^n}\right)_f. \]

方程(16)在收敛的情况下可以写为：

(20)\[ \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. \]

方程(20)与方程(16)相减有（需将\(A\)线性化，类似\(A^{n+1}=A^n\)）：

(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. \]

方程(21)提供了速度修正量与压力修正量的关系。在这里略去临点的影响，方程(21)可以写为：

(22)\[ \mathbf{U}_\mathrm{P}^{'} = - \frac{1}{{{A^{n}_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{'}\bfS_f. \]

方程(22)提供了速度修正与压力修正的一对一关系。将其与方程(16)加和有：

(23)\[ \mathbf{U}_\mathrm{P}^{**} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A^n_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{*}\bfS_f. \]

(24)\[ \mathbf{U}_f^{**} = \mathbf{HbyA}^{*}_f - \frac{1}{{{A^n_f}}} \frac{1}{V_\rP}\left(\sum_f p_f^{*}\bfS_f\right)_f. \]

现在考虑连续性方程，我们希望考虑速度修正以及密度修正后的方程，满足：

(25)\[ \sum_f (\varrho_f^*+\rho_f') \mathbf{U}_f^{**} \cdot \bfS_f = 0. \]

其中\(\rho_f'\)表示修正密度。

注意：

可压缩求解器与不可压缩求解器最重要的区别在于连续性方程出现了密度修正。

方程(25)将用于组建压力方程。为了将连续性方程中的密度变量与压力变量结合起来，不管什么类型的状态方程，均有：

(26)\[ \rho_f'=\psi^{*}p_f' \]

那么收敛后的连续性方程可以写为

(27)\[ \sum_f \left(\varrho_f^*+\psi^{*}p_f'\right)(\mathbf{U}_f^{**}) \cdot \bfS_f=0. \]

在一些状态方程中，\(\psi=0\)，这可以严格的演变成弱可压缩压力方程。一些状态方程中\(\psi\neq 0\)，但如果按照低速流的假定\(\psi\approx 0\)，其也可以演变为弱可压缩压力方程。总之，只要认为方程(25)中的\(\psi\approx 0\)，就可以推导出弱可压缩压力方程。其可以写为：

(28)\[ \sum_f \varrho_f^* \mathbf{U}_f^{**} \cdot \bfS_f = 0. \]

将方程(24)代入到方程(28)中有：

(29)\[ \sum_f \varrho_f^* \mathbf{HbyA}^{*}_f \cdot \bfS_f= \sum_f\frac{\varrho_f^*}{{{A^n_f}}} \frac{1}{V_\rP}\left(\sum_f p_f^{*}\bfS_f\right)_f \cdot \bfS_f. \]

可见，方程(29)可用于求解获得压力\(p^*\)。

最后可以看出，弱可压缩流动的压力方程与不可压缩的压力方程基本一致，在不可压缩流中，连续性方程写为\(\sum_f\mathbf{U}_f^{**} \cdot \bfS_f = 0\)，其与方程(28)的区别就是略去\(\varrho_f^* \)而已。

注意：

\(\frac{\varrho_f^*}{{{A^n_f}}}\)在代码中一般被写为rhorAUf，其中的密度是求解能量方程后更新的密度（可以加适量松弛）。

强可压缩压力方程

在强可压缩的情况下，\(\psi\neq 0\)，密度修正不可忽略。方程(27)可以展开为

(30)\[ \sum_f \left(\left(\varrho_f^*+\psi^{*}p_f'\right)\bfU_f^{*} + \varrho_f^*\bfU_f' \right)\cdot\bfS_f=0. \]

(31)\[ \sum_f \left(\left(\varrho_f^*+\psi^{*}p_f'\right)\bfU_f^{*} \right)\cdot\bfS_f = -\sum_f \left( \varrho_f^*\bfU_f' \right)\cdot\bfS_f \]

注意：

方程(31)中的\(p'\bfU'\)为两个修正量的乘积，被省略。

因此在这里引入略去邻点影响的假定：

(32)\[ \mathbf{U}_\mathrm{P}^{'} = - \frac{1}{V_\rP}\frac{1}{{{A_\mathrm{P}^n}}} \sum_f p_f^{'}\bfS_f. \]

(33)\[ \bfU_f'=-\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f’\bfS_f\right)_f \]

将方程(33)、(18)代入到(31)有：

\[\begin{split} \sum_f \left(\varrho_f^*+\psi^{*}p_f'\right) \mathbf{HbyA}^{*}_f \cdot\bfS_f - \sum_f \left(\varrho_f^*+\psi^{*}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 \varrho_f^* \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(\varrho_f^*+\psi^{*}p_f'\right) \mathbf{HbyA}^{*}_f \cdot\bfS_f - \sum_f \left(\psi^{*}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 \varrho_f^* \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 \varrho_f^* \mathbf{HbyA}^{*}_f \cdot\bfS_f + \sum_f \psi^{*}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 \varrho_f^* \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 \varrho_f^* \mathbf{HbyA}^{*}_f \cdot\bfS_f + \sum_f \psi^{*}(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 \varrho_f^* \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 \varrho_f^* \mathbf{HbyA}^{*}_f \cdot\bfS_f - \sum_f \psi^{*}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^{*}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 \varrho_f^* \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}\]

方程最后的形式为：

(34)\[\begin{split} \sum_f \left(\varrho_f^*-\psi^{*}p_f^n\right) \mathbf{HbyA}^{*}_f \cdot\bfS_f + \sum_f \psi^{*}p_f^*\mathbf{HbyA}^{*}_f\cdot\bfS_f \hspace{17em} \\ + \sum_f \psi^{*}(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 \varrho_f^* \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}\]

方程(34)左边的第三项，形成类似\(\nabla\cdot(p^*\nabla p^n)\)的方程项，其不可离散，在这里将其略去。因此方程(34)可以写为：

(35)\[\begin{multline} \sum_f \varrho_f^* \mathbf{HbyA}^{*}_f \cdot\bfS_f - \sum_f \psi^{*}p_f^n \mathbf{HbyA}^{*}_f \cdot\bfS_f + \sum_f \psi^{*}p_f^*\mathbf{HbyA}^{*}_f\cdot\bfS_f \\ = \sum_f \varrho_f^* \frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{*}\bfS_f\right)_f \cdot\bfS_f \end{multline}\]

其可以用来求解\(p^{*}\)。可见，强可压缩相比较于弱可压缩的区别主要是增加了第二项以及第三项。同时，在\(\psi\approx 0\)的时候，强可压缩流压力方程退化为弱可压缩流压力方程。

速度修正

求解压力后，有：

(36)\[ \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^{*}\)代入到方程(36)有\(\mathbf{U}_\mathrm{P}^{**}\)。同时可以继续更新密度：

(37)\[ \rho^*=f_{\mathrm{EOS}}(\psi^{*},p^*) \]

至此一个迭代步求解完毕，随后进入到下一个迭代步骤。

总之，可压缩算法中的的迭代过程可以表示为下面几个步骤：

通过方程(11)求解获得\(\bfU^*\)，并组建\(\bfHbyA^*\)；
通过能量方程更新可压缩性\(\psi^{*}\)；
低速流通过方程(29)求解获得\(p^*\)，高速流通过方程35求解获得\(p^*\)；
通过方程(37)更新获得\(\rho^{*}\)；
通过方程(36)求解获得\(\bfU^{**}\)；
回到第一步，继续进行求解；

稳态求解器需要对压力进行场松弛，速度、能量等方程进行矩阵松弛增加稳定性。同样的，fluid模块也可以调用SIMPLEC算法。具体内容请参考CFD: 不可压 + 稳态算法。

最后要注意，稳态可压缩算法的连续性误差的计算方法与不可压缩算法相同。在这里要区分于CFD: 可压 + 瞬态算法的连续性误差的计算。

体积力

可压缩求解器在用于模拟暖通环境的时候，要附加浮力的影响。在跨音速的情况下也通常忽略体积力。因此在有体积力的时候，主要考虑弱可压缩的情况。当然了，也可以顺着思路进行跨音速体积力的扩展。

在存在体积力的时候，动量方程中的重力项以及压力项可以进行数值处理。定义

(38)\[ p_\mathrm{rgh}=p-\rho \bfg \cdot\bfh \]

其中\(\bfh\)表示网格单元体心的位置矢量。对方程(38)进行梯度操作：

(39)\[ \nabla p_\mathrm{rgh}=\nabla p-\bfg\cdot \bfh \nabla \rho - \rho \bfg \]

即为

(40)\[ -\nabla p+\rho \bfg=-\bfg\cdot \bfh \nabla \rho -\nabla p_\mathrm{rgh} \]

将方程(40)代入到方程(2)中有：

(41)\[ \nabla \cdot (\rho\mathbf{U} \mathbf{U})=-\bfg\cdot \bfh \nabla \rho -\nabla p_\mathrm{rgh}+\nabla\cdot\tau \]

接下来看压力方程的处理。针对方程(32)，在考虑方程(41)后演变为：

(42)\[ \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) \]

方程(42)提供了速度修正与压力修正的一对一关系。但发现，其中还存在密度修正。如果顺着这个思路进行，压力方程将进一步出现未知变量密度。为了处理这个问题，我们可以将方程(42)中的密度修正忽略掉。这是一种方法，但处理起来比较激进。

另外一种方法是在能量方程更新可压缩性\(\psi^{*}\)后，在其基础上，通过状态方程更新密度中间量：

(43)\[ \rho^{*}=f_{\mathrm{EOS}}(p^n, \psi^{*}) \]

方程(43)在前面也已经出现过。在这个情况下，认为\(\rho^n\)与密度修正量\(\rho^{'}\)的关系，满足：

(44)\[ \rho^{*}=\rho^n+\rho^{'} \]

这样，由于\(\rho^{*}\)已经通过方程(43)求出，因此压力方程仅仅存在未知量压力。依据此思想，结合修正方程(42)有：

(45)\[ \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^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f| \right). \]

(46)\[ \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^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f| \right)_f. \]

(47)\[ \mathbf{U}^{**} =\mathbf{U}^{*}+\bfU^{'} \]

将方程(46)代入到方程(28)中有压力泊松方程：

(48)\[\begin{split} \sum_f \rho_f^{*} \mathbf{HbyA}^{*}_f \cdot \bfS_f - \sum_f \frac{\rho^{*}_f}{{{A^n_f}}} \frac{1}{V_\rP} \left( \sum_f \bfg_f\cdot \bfh_f \rho_f^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f| \right)_f \cdot \bfS_f \hspace{5em} \\ = \sum_f\frac{\rho^{*}_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}\]

方程(48)的连续形式为：

(49)\[ \nabla\cdot\left(\rho^{*} \left(\mathbf{HbyA} - \frac{1}{A}\bfg\cdot\bfh\nabla\rho^* \right) \right) = \nabla\cdot\left(\rho^{*}\frac{1}{A} \nabla p_{\mathrm{rgh}}^*\right) \]

方程(49)可用于求解并获得存在体积力的压力\(p^*\)。

同时，在考虑体积力的时候，更新速度的时候要考虑体积力的影响：

(50)\[ \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^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f| \right) \]

关键代码

fluid模块中的速度方程并无特殊。在速度方程更性之后，fluid模块中会求解能量方程，并通过下述代码进行\(\psi,C_p,C_v\)等变量的更新：

thermo.correct();

强可压缩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())
);

其为一个拉普拉斯方程。在求解后，需要对密度进行更新：

rho = thermo.rho();

若考虑体积力，方程(45)中速度方程的更新采用了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^{*}_f}{A^n_\rP}\right)_f \bfg_f\cdot\bfh_f (\nabla\rho^{*})_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^{*}_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^{*}}{A^n_\rP}\right)_f \bfg_f\cdot\bfh_f (\nabla\rho^{*})_f\frac{\bfS_f}{|\bfS_f|}|\bfS_f| - \left(\frac{\rho^{*}}{A^n_\rP}\right)_f (\nabla p_{\mathrm{rgh}})_f \frac{\bfS_f}{|\bfS_f|}|\bfS_f| \right)\frac{(A^n_\rP)_f}{\rho^{*}_f} \right) \]

约去相关项并展开reconstruct()函数：

\[ + \frac{1}{A^n_\rP} \frac{\sum \bfn_f \cdot\left( - \bfg_f\cdot\bfh_f (\nabla\rho^{*})_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^{*} -\frac{1}{A^n} \nabla p_{\mathrm{rgh}} \]

更详细的关于reconstruct()的解释，请参考《无痛苦NS方程笔记》。