控制方程与算法
buoyantSimpleFoam中求解的为下述稳态质量方程、动量方程、以及状态方程：
(1)\[
\nabla\cdot\rho\bfU=0,
\]
(2)\[
\nabla \cdot (\rho\mathbf{U} \mathbf{U})=-\nabla p+\nabla \cdot\tau +\rho \bfg,
\]
(3)\[
p=\rho RT.
\]
其中方程(1)没有时间项，并无明显的变量，更像是一种限定性条件。因此稳态求解器中不存在对质量方程的求解。同时，方程(1)中的\(\rho\bfU\)可用于组建质量通量。
同时，buoyantSimpleFoam中能量方程在求解能量变量后，对其直接求解即可，不涉及到速度-压力-密度耦合问题，此处略。上述三个方程，可用于迭代求解三个未知量：速度、压力、密度。
首先，方程(7)中的重力项以及压力项可以进行数值处理。定义
(4)\[
p_\mathrm{rgh}=p-\rho \bfg \cdot\bfh
\]
其中\(\bfh\)表示网格单元体心的位置矢量。对方程(4)进行梯度操作：
(5)\[
\nabla p_\mathrm{rgh}=\nabla p-\bfg\cdot \bfh \nabla \rho - \rho \bfg
\]
即为
(6)\[
-\nabla p+\rho \bfg=-\bfg\cdot \bfh \nabla \rho -\nabla p_\mathrm{rgh}
\]
将方程(4)代入到方程(2)中有：
(7)\[
\nabla \cdot (\rho\mathbf{U} \mathbf{U})=-\bfg\cdot \bfh \nabla \rho -\nabla p_\mathrm{rgh}+\nabla\cdot\tau.
\label{mom}
\]
对方程(7)通过高斯定理进行对速度\(\bfU\)的离散，组建速度方程有：
(8)\[
\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^{*} ,
\]
(9)\[
\int \nabla p_\mathrm{rgh} \mathrm{d}V
=\int p_\mathrm{rgh} \mathrm{d}\bfS
=\sum_f p_{\mathrm{rgh},f}^n\bfS_f,
\]
(10)\[
\int \bfg\cdot \bfh \nabla \rho \mathrm{d}V
=\int \bfg_f\cdot \bfh_f \rho_f \mathrm{d}\bfS
=\sum_f \bfg_f\cdot \bfh_f \rho_f^n\bfS_f,
\]
其中上标\(^n\)表示为当前迭代步（已知），上标\(^{*}\)表示下一个迭代步（待求），下标\(_f\)表示网格单元面上的值，\(\bfS_f\)表示网格单元的各个面的面矢量，\(F_f\)为质量通量，\(\mu\)为运动粘度（假定粘度为常数）。此处略去粘性项\(\tau\)的离散且假定粘度为\(0\)。需要注意的是，目前并没有处理\(\nabla\cdot\tau\)这一项，因为其对求解流程并不影响。为了防止方程(9)与方程(10)的离散产生数值振荡。其一般采取面法向梯度的形式来进行：
(11)\[
\int \nabla p_\mathrm{rgh} \mathrm{d}V
=\sum_f p_{\mathrm{rgh},f}^n\bfS_f
=\sum_f p_{\mathrm{rgh},f}^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|,
\]
(12)\[
\int \bfg\cdot \bfh \nabla \rho \mathrm{d}V
=\sum_f \bfg_f\cdot \bfh_f \rho_f^n\bfS_f
=\sum_f \bfg_f\cdot \bfh_f \rho_f^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|,
\]
注意：
以\(\int \nabla p_\mathrm{rgh} \mathrm{d}V\)这一项举例。其常规离散方法，如方程(9)，可以写为fvc::grad(p_rgh)。若参用面法向梯度的形式，如方程(11)，可以写为fvc::reconstruct(fvc::snGrad(p_rgh)*mesh.magSf())。
整理方程(11)、(12)、(8)有：
(13)\[
\sum_f F_f^n \mathbf{U}_f^{*}
= -\sum_f p_{\mathrm{rgh},f}^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
-\sum_f \bfg_f\cdot \bfh_f \rho_f^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|.
\]
需要注意的是，方程(13)左侧的对流项采用了线性化处理。方程(13)中的\(\bfU_f\)需要从体心速度进行插值来获得，在此步可以引入各种插值格式。假设使用中心线性格式（均一网格）：
(14)\[
\mathbf{U}_f^{*} = \frac{{\mathbf{U}_\mathrm{P}^{*} + \mathbf{U}_\mathrm{N}^{*}}}{2}.
\]
将方程(14)代入到方程(13)有
(15)\[
\sum_f {F_f^n \frac{{\mathbf{U}_\mathrm{P}^{*} + \mathbf{U}_\mathrm{N}^{*}}}{2}} = -\sum_f p_{\mathrm{rgh},f}^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
-\sum_f \bfg_f\cdot \bfh_f \rho_f^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|.
\]
(16)\[
\sum_f { {{\frac{{F_f^n}}{2}} } } \mathbf{U}_\mathrm{P}^{*}
+
\sum_f { { {\frac{{F_f^n}}{2} } \mathbf{U}_\mathrm{N}^{*}} }
=
-\sum_f p_{\mathrm{rgh},f}^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
-\sum_f \bfg_f\cdot \bfh_f \rho_f^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|.
\]
将上式左右两边同时除以网格单元体积\(V_\rP\)，同时简写为
(17)\[
{A_\mathrm{P}^n}\mathbf{U}_\mathrm{P}^{*}{\rm{ + }}\sum_f {A_\mathrm{N}^n\mathbf{U}_\mathrm{N}^{*}}
=
-\frac{1}{V_\rP}
\left(
\sum_f p_{\mathrm{rgh},f}^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
+\sum_f \bfg_f\cdot \bfh_f \rho_f^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right),
\]
其中\(A_\rP^n\)，\(A_\rN^n\)分别表示当前网格点与相邻网格点的离散系数：
(18)\[
A_\mathrm{P}^n= \frac{1}{V_\rP}\sum_f \frac{F_f^n}{2} ,
\]
(19)\[
A_\mathrm{N}^n= \frac{1}{V_\rP}{\frac{{F_f^n}}{2}} ,
\]
可见，\(A_\mathrm{P}\)与\(A_\mathrm{N}\)为关于\(\rho^n\)的代数式。方程(17)构成一个稀疏线性系统。求解方程(17)即可获得速度\(\mathbf{U}^{*}\)。对方程(17)进行转换有：
(20)\[
\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}^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
+\sum_f \bfg_f\cdot \bfh_f \rho_f^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right).
\]
其中\(\bfHbyA\)定义为：
(21)\[
\bfHbyA_\rP^{*}=\frac{-\sum_NA_\rN^n\bfU_\rN^{*}}{A_\rP^n}.
\]
同时，面上插值有：
(22)\[
\mathbf{U}_f^{*} = \mathbf{HbyA}^{*}_f - \frac{1}{{{A_f^n}}}
\left(
\frac{1}{V_\rP}
\left(
\sum_f p_{\mathrm{rgh},f}^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
+\sum_f \bfg_f\cdot \bfh_f \rho_f^n\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right)
\right)_f,
\]
(23)\[
\bfHbyA_f^{*}=\left(\frac{-\sum_NA_\rN^n\bfU_\rN^{*}}{A_\rP^n}\right)_f.
\]
注意：
在求解速度方程后，求解器求解能量方程，并进行thermo.correct()函数，对\(\psi=1/RT\)进行更新。因此，在下文中，用\(\psi^{corr}\)表示更新后的\(1/RT\)。
在低速流中，可以认为密度的变化并不是特别的大。方程(20)在收敛的情况下可以写为：
(24)\[
\mathbf{U}_\mathrm{P}^{n+1} = \mathbf{HbyA}^{n+1}_\mathrm{P} - \frac{1}{{{A^{n+1}_\mathrm{P}}}} \frac{1}{V_\rP}
\left(
\sum_f p_{\mathrm{rgh},f}^{n+1}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
+\sum_f \bfg_f\cdot \bfh_f \rho_f^{n+1}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right).
\]
在这里需要做一系列的假定：
略去临点的影响；
假定\(A_\rP\)的变化较小（\(A^{n+1}_\mathrm{P}=A^{n}_\mathrm{P}\)）；
方程(24)与方程(20)相减有：
(25)\[
\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).
\]
方程(25)需要提供了速度修正与压力修正的一对一关系。但发现，其中还存在密度修正。如果顺着这个思路进行，压力方程将进一步出现未知变量密度。为了处理这个问题，我们可以将方程(25)中的密度修正忽略掉。这是一种方法，但处理起来比较激进。另外一种方法是，参考rhoSimpleFoam解析，由于在速度方程之后要进行能量方程的计算，能量方程会更新可压缩性\(\psi^{corr}\)，在更新可压缩性后，在\(\psi^{corr}\)的基础上，可以通过状态方程更新密度中间量：
(26)\[
\rho^{corr}=p^n/\psi^{corr}
\]
从而，认为\(\rho^n\)与密度修正量\(\rho^{'}\)的关系，满足：
(27)\[
\rho^{corr}=\rho^n+\rho^{'}
\]
而非\(\rho^{n+1}=\rho^n+\rho^{'}\)。这样，由于\(\rho^{corr}\)已经通过方程(10)求出，因此压力方程仅仅存在未知量压力。依据此思想，将方程(25)与方程(20)加和有：
(28)\[
\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).
\]
(29)\[
\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.
\]
(30)\[
\mathbf{U}^{**} =\mathbf{U}^{*}+\bfU^{'}
\]
在这里，低速流认为\(\mathbf{U}_f^{**}\)需要满足连续性方程：
(31)\[
\sum_f \rho_f^{corr}\mathbf{U}_f^{**} \cdot \bfS_f=0.
\]
注意：
高速流则同样满足
\[
\sum_f (\rho_f^{n}+\rho_f')(\mathbf{U}_f^{*}+\bfU^{'} )\cdot \bfS_f=0.
\]
但需要将\(\rho_f^{corr}\)分解为\(\rho_f^{n}+\rho_f'\)处理。有关高速流的推导请参考rhoSimpleFoam解析。
同样的，方程(15)中的密度采用已知的密度（低速流密度变动充分小）。将方程(12)代入到方程(15)中有：
(32)\[\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}\]
可见，方程(16)可用于求解获得压力\(p^*\)。将\(p^{*}\)代入到方程(12)有\(\mathbf{U}_\mathrm{P}^{**}\)。随后，密度再次通过状态方程进行更新：
(33)\[
\rho^*=\psi^{corr}p^*
\]
总之，buoyantSimpleFoam中的的迭代过程可以表示为下面几个步骤：
通过方程(20)求解获得\(\bfU^*\)，并组建\(\bfHbyA^*\)；
通过能量方程更新可压缩性\(\psi^{corr}\)；
通过状态方程更新中间变量密度\(\rho^{corr}\)；
通过方程(16)求解获得\(p_{\mathrm{rgh}}^*\)；
通过方程(12)求解获得\(\bfU^{**}\)；
通过方程(18)更新获得\(\rho^{*}\)；
回到第一步，继续进行求解；
在这里需要注意的是，buoyantSimpleFoam为一个稳态求解器。同样的需要对压力进行场松弛，速度、能量等方程进行矩阵松弛增加稳定性。同样的，buoyantSimpleFoam也可以调用SIMPLEC算法，但求解器并未进行植入。不同于rhoSimpleFoam、simpleFoam，buoyantSimpleFoam主要考虑了浮力源项的作用，这主要带来了两点变化：
附加浮力的求解器压力变量变为了\(p_{\mathrm{rgh}}\)；
附加除了压力项之外的求解器，速度修正与压力修正的关系会变得更加复杂，如方程(25)；
为了防止修正量过大，密度修正、压力修正都可以进行场松弛；
需要注意的是。buoyantSimpleFoam在整个迭代过程中，为一个迭代求解器，因此对于各种变量的更新，某些过程可以忽略。一些可以忽略的过程比如有：
速度方程可以不进行求解；
密度中间变量可以不进行更新；
那么具体的，是否需要进行这些必要操作，是问题依赖的。但这并不影响最终向收敛过程的推进。