版本对应：
OpenFOAM-11中的fluid模块，对应OpenFOAM-10之前的rhoSimpleFoam、rhoPimpleFoam、buoyantFoam。其可以处理可压缩流体以及浮力驱动流。fluid模块主要基于isothermalFluid模块，仅仅在其基础上添加能量方程。在OpenFOAM-11之前，有无体积力采用不同的求解器来处理，buoyantSimpleFoam（有体积力），rhoSimpleFoam（无体积力）。在OpenFOAM-11都被处理到fluid模块中来减少代码复用。
因此本篇文章讨论的算法适用于OpenFOAM新版的fluid模块。以及rhoSimpleFoam、buoyantFoam（稳态部分）、buoyangSimpleFoam。
CFD: 可压 + 瞬态
可压缩求解器相对于不可压缩求解器还可以继续进行区分：
是否需要考虑体积力（比如温度变化引起的浮力）；
弱可压缩还是强可压缩（前者典型的为人居环境的温度变化，后者主要是超音速流动）；
对于超音速流动，也即可压缩性特别强的情况，得益于双曲方程的数学特征，这些方法通常采用瞬态计算且使用密度作为主要变量，然后通过状态方程求解压力，即密度基求解器。但密度基求解器在应用于弱可压缩流的情况下效率低下。一方面的解释是在不可压缩领域，压力与密度的耦合非常弱，另一方面的解释是不可压缩假定下音速趋向于无穷大，导致时间步长过小。
针对弱可压缩领域，fluid模块为一个压力基、可压缩、适用于全流速、也包含体积力的求解器，主要用于求解普适性可压缩流动，也可用于求解低音速/超音速流动。其可以处理瞬态问题，也可以处理稳态问题。本文讨论瞬态算法。
控制方程
fluid模块（buoyantFoam，甚至更老的buoyantPimpleFoam）中求解的为下述瞬态连续性方程、动量方程、以及状态方程：
其实还需要附加能量方程。但能量方程不涉及到速度-压力-密度耦合问题，此处略。上述三个方程，可用于迭代求解三个未知量：速度、压力、密度。首先，方程(2)中的重力项以及压力项可以进行数值处理。定义
其中\(\bfh\)表示网格单元体心的位置矢量。对方程(4)进行梯度操作：
即为
方程(7)为求解的动量方程。
离散
不同与稳态连续性方程，方程(1)可以用于更新密度。因此方程(1)可以离散为：
如果对时间项采用Euler格式进行离散，空间项隐性格式有：
方程(9)存在两个未知量：密度\(\rho^{t+\dt}\)与速度\(\bfU^{t+\dt}\)。不可解。因此，在对对流项离散的过程中，需要使用当前时间步的速度\(\bfU^{t}\)。在这种情况下，离散方程可以写为：
因为对流项的离散并没有采用\(\bfU^{t+\dt}\)，方程(10)求解的密度并非当前时间步最终收敛的密度\(\rho^{t+\dt}\)。
注意：
方程(1)若需要去掉时间项，就变成了稳态的连续型方程。稳态的连续性方程并无明显的变量，是一种限定性条件。
对方程(7)通过高斯定理进行对速度\(\bfU\)的离散，组建速度方程有：
Warning
方程(11)右侧之所以\(\rho\)带星号，是因为密度方程已经求结果。在这里是一个已知量。
Warning
方程(12)\(\rho\)带星号与之前一样的结果。第一次迭代的时候\(\bfU\)采用当前\(t\)的值。因此\(F_f^t=\rho_f^*\bfU_f^t\cdot\bfS_f\)。
其中上标\(^t\)表示为当前时间步（已知），出来密度外的上标\(^{*}\)表示当前时间步下的迭代步（待求），下标\(_f\)表示网格单元面上的值，\(\bfS_f\)表示网格单元的各个面的面矢量，\(F_f\)为质量通量。此处略去粘性项\(\tau\)的离散因为其不影响压力速度耦合。
为了防止方程(13)与方程(14)的离散产生数值振荡。其一般采取面法向梯度的形式来进行：
需要注意的是，方程(17)左侧的对流项采用了线性化处理。方程(17)中的\(\bfU_f\)需要从体心速度进行插值来获得，在此步可以引入各种插值格式。假设使用中心线性格式（均一网格）：
其可简写为
其中\(A_\rP^{t}\)，\(A_\rN^{t}\)分别表示当前网格点与相邻网格点的离散系数：
可见，\(A_\mathrm{P}\)与\(A_\mathrm{N}\)为关于\(\rho^{*}\)的代数式（且当前为已知量）。方程(20)构成一个稀疏线性系统。求解方程(20)即可获得速度\(\mathbf{U}^{*}\)。
对方程(20)进行转换有：
其中\(\bfHbyA\)定义为：
同时，面插值有：
注意：
在求解速度方程后，求解器求解能量方程，并进行thermo.correct()函数，对\(\psi=1/RT\)进行更新。因此，在下文中，用\(\psi^{corr1}\)表示更新后的\(1/RT\)。
方程(23)在收敛的情况下可以写为：
在这里需要做一系列的假定：
略去临点的影响；
\(A^{t+\dt}_\mathrm{P}=A^{t}_\mathrm{P}\)，否则其为一个非线性方程组不可解，也可以认为\(\dt\)间隔比较小的时候忽略其影响；
方程(28)需要提供了速度修正与压力修正的一对一关系。但发现针对浮力驱动流，其中还存在密度修正。如果顺着这个思路进行，压力方程将进一步出现未知变量密度导致无法封闭。
Warning
这也是可压缩求解器与不可压缩求解器的主要区别。
为了处理这个问题，我们可以将方程(28)中的密度修正忽略掉。这样，即存在\(\mathbf{U}_f^{**}\)与\(p^{*}\)的一对一关系，结合连续性方程，可以用来求解压力。这是一种方法，但处理起来比较激进，可能会导致变量增量过大，引起发散。
如果想考虑密度修正的话，则在压力方程中，会存在一个未知的密度变量导致不封闭。但同时需要注意的是，密度变量也可以从状态方程中得出，目前状态方程我们并没有使用。只是在求解能量方程后，更新了可压缩性得到了\(\psi^{corr1}\)。在此基础上，其实可以通过状态方程更新密度中间量：
在这里，可以认为\(\rho^{corr1}\)与密度修正量\(\rho^{'}\)的关系满足：
由于\(\rho^{corr1}\)已经通过方程(29)求出，因此压力方程仅仅存在未知量压力。依据此思想，将方程(28)与方程(23)加和有：
注意：
这里\(\mathbf{U}^{**} =\mathbf{U}^{*}+\bfU^{'}\)，而不是\(\mathbf{U}^{t+\dt} =\mathbf{U}^{*}+\bfU^{'}\)，这是因为已经做了一系列的假定。
同时，方程(32)中的密度已经被求出，因此可以在后续中用来组建压力方程。
压力方程
首先考虑低速流的情况，低速流认为\(\mathbf{U}_f^{**}\)需要满足连续性方程：
现在我们通过方程(33)来构造压力柏松方程。首先，方程(33)要对当先的速度\(\bfU^{**}\)进行实施，来构成一种限定性条件，即：
再次，方程(34)中的对流项存在未知变量密度\(\rho^{t+\dt}\)，因此要做线性化：
同样，方程(35)的时间项存在未知的变量\(\rho^{t+\dt}\)。OpenFOAM采用了一种巧妙的思想来进行设计，在这里假定读者是熟悉PISO算法需要若干次迭代的：
在进行1次迭代后：
在进行2次迭代后：
同时要注意，\(\rho^{corr}\)就是需要进行多次修正，然后趋向于\(\rho^{t+\dt}\)。假定PISO算法做了2次迭代认为\(\rho^{corr3}\)趋向于\(\rho^{t+\dt}\)。则方程方程(39)就变为：
参考这个思想，结合状态方程，将方程(39)代入到方程(35)中有：
注意：
这里\(\frac{p^{*}_\mathrm{rgh}-p^t_\mathrm{rgh}}{\dt}\)对应correction(fvc::ddt(p_rgh))函数。correction()函数类似于fvm::ddt()函数，不同的是fvm::ddt()在同一个时间步之内，一直表示的为\((p^{*}-p^t)/\dt\)，其中\(p^{*}\)表示待求变量。若进行一次压力求解，fvm::ddt()变为\((p^{**}-p^t)/\dt\)。继续进行压力求解，则变为\((p^{***}-p^t)/\dt\)。correction(fvm::ddt())则在最初表示的为\((p^{*}-p^t)/\dt\)。进行一次压力求解，变成\((p^{**}-p^{*})/\dt\)。再进行一次压力求解，变成\((p^{***}-p^{**})/\dt\)。同时也可以看出，在最终一个时间步收敛的情况下，\(p^{***}\)趋向于\(p^{**}\)，因此correction(fvm::ddt())作用趋向于0。有关correction()函数的详细内容请参考无痛苦N-S方程笔记。
可见，方程(40)可用于求解获得压力\(p_\mathrm{rgh}^*\)。将\(p_\mathrm{rgh}^*\)代入到方程(32)有\(\mathbf{U}_\mathrm{P}^{**}\)。通过更新后的压力\(p_{\mathrm{rgh},f}^{*}\)以及速度\(\mathbf{U}_\mathrm{P}^{**}\)，依据密度更新的两种方法，进行下述两步：
依据状态方程更新密度。有压力修正量为：
(41)\[ p_{\mathrm{rgh},f}^{'}=p_{\mathrm{rgh},f}^{*}-p_{\mathrm{rgh},f}^{t} \]基于压力修正量，可以进一步更新密度：
(42)\[ \rho^{corr2}=\rho^{corr1}+\psi^{*}(p_{\mathrm{rgh},f}^{*}-p_{\mathrm{rgh},f}^{t}) \]依据连续性方程更新密度。通过速度\(\mathbf{U}_\mathrm{P}^{**}\)组建新的通量，求解连续性方程。
(43)\[ \frac{\rho^{**}-\rho^t}{\Delta t}+\sum_f \rho^{corr1}_f\bfU_f^{**}\cdot\bfS_f=0, \]
上述两种更新密度的方法，在一个时间步下收敛后，必然一致。因此，可以通过\(\rho^{**}\)与\(\rho^{corr2}\)的差来判断当前时间步下是否收敛：
至此，我们有了满足连续性方程的速度场\(\mathbf{U}_\mathrm{P}^{**}\)。但由于在上述推导过程中，引入了大量的假设，\(\mathbf{U}_\mathrm{P}^{**}\)与\(p^{*}\)并不是严格满足动量方程。同时，方程(44)也未必满足。在这里可以把求得的这些场当作初始场，在一个时间步内再进行迭代求解。直到方程(44)满足为止。
总之，buoyantPimpleFoam中的的迭代过程可以表示为下面几个步骤：
通过方程(10)求解获得\(\rho^*\)；
通过\(\rho^*\)，组建速度方程，获得\(\bfU^*\)；
组建能量方程，更新温度等变量，更新可压缩性\(\psi^{corr}\)，通过状态方程更新\(\rho^{corr1}\)；
通过第二步的\(\bfU^*\)组建\(\bfHbyA^*\)，构建压力柏松方程(40)，求解获得\(p_{\mathrm{rgh}}^*\)；
通过\(p_{\mathrm{rgh}}^*\)更新速度\(\bfU^{**}\)，求解方程(43)，获得\(\rho^{**}\)；
通过方程(42)，获得\(\rho^{corr2}\)；
通过方程(44)判断连续性误差；
回到第二步进行迭代，直至通过密度方程求解的\(\rho^{**}\)与状态方程更新的密度\(\rho^{corr2}\)趋向于一致为止；
迭代为止之后，令通过密度方程求解的\(\rho\)与状态方程更新的密度\(\rho^{corrn}\)完全相等；
关键代码
针对上述求解流程，流程第一步，对应于代码中最开始的密度方程求解：
#include "rhoEqn.H"
流程第二步，需要组建速度方程。buoyantPimpleFoam中的速度方程左侧并无特殊。但方程右侧，采用reconstruct()函数保证无振荡：
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
上述代码的源项对应方程(15)、(16)。在速度方程更性之后，求解能量方程，进入流程第三步，并通过下述代码进行可压缩性的更新：
thermo.correct();
以及密度的更新：
rho = thermo.rho();
流程第四步，需要组建压力柏松方程：
fvScalarMatrix p_rghDDtEqn
(
psi*correction(fvm::ddt(p_rgh))
+ fvc::ddt(rho)
+ fvc::div(phi)
- fvm::laplacian(rhorAUf, p_rgh)
);
上述代码对应方程(40)。流程第五步，在求解压力之后，更新速度，更新密度：
phi += p_rghEqn.flux();
...
#include "rhoEqn.H"
流程第六步，更新\(\rho^{corr2}\)：
thermo.rho() -= psi*p_rgh;
...
thermo.rho() += psi*p_rgh;
流程第七步，判断连续性误差：
#include "compressibleContinuityErrs.H"
流程第八步，继续进行迭代。流程第九步，令通过密度方程求解的\(\rho\)与状态方程更新的密度\(\rho^{corrn}\)完全相等：
rho = thermo.rho();