控制方程与算法
buoyantPimpleFoam中求解的为下述瞬态连续性方程、动量方程、以及状态方程：
(1)\[
\frac{\p\rho}{\p t}+\nabla\cdot\rho\bfU=0,
\]
(2)\[
\frac{\p\rho\bfU}{\p t}+\nabla \cdot (\rho\mathbf{U} \mathbf{U})=-\nabla p+\nabla \cdot\tau +\rho \bfg,
\]
(3)\[
p=\rho RT.
\]
不同与稳态连续性方程，方程(1)看起来可以用于更新密度。方程(1)可以离散为：
(4)\[
\int\int\left(\frac{\p\rho}{\p t}+\nabla\cdot\rho\bfU\right)\rd V\rd t=0,
\]
(5)\[
\frac{\rho^{t+\dt}-\rho^t}{\Delta t}+\sum_f \rho^{t+\dt}_f\bfU_f^{t+\dt}\cdot\bfS_f=0,
\]
方程(5)存在两个未知量：密度\(\rho^{t+\dt}\)与速度\(\bfU^{t+\dt}\)。不可解。因此，在对对流项离散的过程中，需要使用当前时间步的速度\(\bfU^{t}\)。在这种情况下，离散方程可以写为：
(6)\[
\frac{\rho^{*}-\rho^t}{\Delta t}+\sum_f \rho^t_f\bfU_f^t\cdot\bfS_f=0,
\]
因为对流项的离散并没有采用\(\bfU^{t+\dt}\)，方程(6)求解的密度并非当前时间步最终收敛的密度\(\rho^{t+\dt}\)。
注意：
方程(1)若需要去掉时间项，就变成了稳态的连续型方程。稳态的连续性方程并无明显的变量，是一种限定性条件。
buoyantPimpleFoam中的能量方程在求解能量变量后，对其直接求解即可，不涉及到速度-压力-密度耦合问题，此处略。上述三个方程，可用于迭代求解三个未知量：速度、压力、密度。首先，方程(2)中的重力项以及压力项可以进行数值处理。定义
(7)\[
p_\mathrm{rgh}=p-\rho \bfg \cdot\bfh
\]
其中\(\bfh\)表示网格单元体心的位置矢量。对方程(7)进行梯度操作：
(8)\[
\nabla p_\mathrm{rgh}=\nabla p-\bfg\cdot \bfh \nabla \rho - \rho \bfg
\]
即为
(9)\[
-\nabla p+\rho \bfg=-\bfg\cdot \bfh \nabla \rho -\nabla p_\mathrm{rgh}
\]
将方程(7)代入到方程(2)中有：
(10)\[
\frac{\p\rho\bfU}{\p t}+\nabla \cdot (\rho\mathbf{U} \mathbf{U})=-\bfg\cdot \bfh \nabla \rho -\nabla p_\mathrm{rgh}+\nabla\cdot\tau.
\label{mom}
\]
对方程(10)通过高斯定理进行对速度\(\bfU\)的离散，组建速度方程有：
(11)\[
\int\int \frac{\p\rho\bfU}{\p t} \rd V \rd t =
V_\rP \rho_\rP^{*}(\bfU^{*}_\rP - \bfU^t_\rP),
\]
(12)\[
\int\int {\nabla \cdot \left( {\rho\mathbf{U}\mathbf{U}} \right)\mathrm{d}V\rd t = } \Delta t\sum_f {{{\left( {\rho^{t}{\mathbf{U}^t}{\mathbf{U}^{*}}} \right)}_f}} \cdot\bfS_f = \Delta t\sum_f F_f^{t} \mathbf{U}_f^{*} ,
\]
(13)\[
\int\int \nabla p_\mathrm{rgh} \mathrm{d}V\rd t
=\Delta t\int p_\mathrm{rgh} \mathrm{d}\bfS
=\Delta t\sum_f p_{\mathrm{rgh},f}^t\bfS_f,
\]
(14)\[
\int\int \bfg\cdot \bfh \nabla \rho \mathrm{d}V\rd t
=\Delta t\int \bfg_f\cdot \bfh_f \rho_f \mathrm{d}\bfS
=\Delta t\sum_f \bfg_f\cdot \bfh_f \rho_f^{*}\bfS_f,
\]
其中上标\(^t\)表示为当前时间步（已知），上标\(^{*}\)表示当前时间步下的迭代步（待求），下标\(_f\)表示网格单元面上的值，\(\bfS_f\)表示网格单元的各个面的面矢量，\(F_f\)为质量通量。此处略去粘性项\(\tau\)的离散且假定粘度为\(0\)。需要注意的是，目前并没有处理\(\nabla\cdot\tau\)这一项，因为其对求解流程并不影响。为了防止方程(13)与方程(14)的离散产生数值振荡。其一般采取面法向梯度的形式来进行：
(15)\[
\Delta t\sum_f p_{\mathrm{rgh},f}^t\bfS_f
=\Delta t\sum_f p_{\mathrm{rgh},f}^t\frac{\bfS_f}{|\bfS_f|}|\bfS_f|,
\]
(16)\[
\Delta t\sum_f \bfg_f\cdot \bfh_f \rho_f^{*}\bfS_f
=\Delta t\sum_f \bfg_f\cdot \bfh_f \rho_f^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|,
\]
整理方程(15)、(16)、(12)、(11)有：
(17)\[\begin{split}
\frac{\rho_\rP^{*}V_\mathrm{P}}{\Delta t} (\bfU^{*}_\rP - \bfU^t_\rP)+\sum_f F_f^{t} \mathbf{U}_f^{*}
\hspace{20em}
\\
= -\sum_f p_{\mathrm{rgh},f}^t\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
-\sum_f \bfg_f\cdot \bfh_f \rho_f^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|.
\end{split}\]
需要注意的是，方程(17)左侧的对流项采用了线性化处理。方程(17)中的\(\bfU_f\)需要从体心速度进行插值来获得，在此步可以引入各种插值格式。假设使用中心线性格式（均一网格）：
(18)\[
\mathbf{U}_f^{*} = \frac{{\mathbf{U}_\mathrm{P}^{*} + \mathbf{U}_\mathrm{N}^{*}}}{2}.
\]
将方程(18)代入到方程(17)有
(19)\[\begin{split}
\sum_f \left(\frac{F_f^{t}}{2V_\rP} + \frac{\rho_\rP^{*} }{\Delta t} \right) \mathbf{U}_\mathrm{P}^{*}
+
\sum_N { { {\frac{{F_f^{t}}}{2V_
\rP} } \mathbf{U}_\mathrm{N}^{*}} }
\hspace{25em}\\
=
-\frac{1}{V_
P}\sum_f p_{\mathrm{rgh},f}^t\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
-\frac{1}{V_
P}\sum_f \bfg_f\cdot \bfh_f \rho_f^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
+\frac{\rho_\rP^{*} }{\Delta t} \bfU^t_\rP.
\hspace{7em}
\end{split}\]
其可简写为
(20)\[\begin{split}
{A_\mathrm{P}^{t}}\mathbf{U}_\mathrm{P}^{*}{\rm{ + }}\sum_N {A_\mathrm{N}^{t}\mathbf{U}_\mathrm{N}^{*}}
\hspace{27em}
\\
=
-\frac{1}{V_\rP}
\left(
\sum_f p_{\mathrm{rgh},f}^t\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
+\sum_f \bfg_f\cdot \bfh_f \rho_f^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right)+\frac{1}{\Delta t} \bfU^t_\rP,
\hspace{2em}
\end{split}\]
其中\(A_\rP^{t}\)，\(A_\rN^{t}\)分别表示当前网格点与相邻网格点的离散系数：
(21)\[
A_\mathrm{P}^{t}= \frac{\rho_\rP^{*} }{\Delta t}+\sum_f \frac{F_f^{t}}{2V_
P} ,
\]
(22)\[
A_\mathrm{N}^{t}= \frac{1}{V_\rP}{\frac{{F_f^{t}}}{2}} ,
\]
可见，\(A_\mathrm{P}\)与\(A_\mathrm{N}\)为关于\(\rho^{t}\)的代数式。方程(20)构成一个稀疏线性系统。求解方程(20)即可获得速度\(\mathbf{U}^{*}\)。对方程(20)进行转换有：
(23)\[
\mathbf{U}_\mathrm{P}^{*} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A^{t}_\mathrm{P}}}} \frac{1}{V_\rP}
\left(
\sum_f p_{\mathrm{rgh},f}^t\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
+\sum_f \bfg_f\cdot \bfh_f \rho_f^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right).
\]
其中\(\bfHbyA\)定义为：
(24)\[
\bfHbyA_\rP^{*}=\frac{1}{A_\rP^{t}}\left(-\sum_NA_\rN^{t}\bfU_\rN^{*}+\frac{1}{\Delta t} \bfU^t_\rP\right).
\]
同时，面插值有：
(25)\[
\mathbf{U}_f^{*} = \mathbf{HbyA}^{*}_f - \frac{1}{{{A_f^{t}}}}
\left(
\frac{1}{V_\rP}
\left(
\sum_f p_{\mathrm{rgh},f}^t\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
+\sum_f \bfg_f\cdot \bfh_f \rho_f^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right)
\right)_f,
\]
(26)\[
\bfHbyA_f^{*}=\frac{1}{A_f^{t}}\left(-\sum_NA_\rN^{t}\bfU_\rN^{*}+\frac{1}{\Delta t} \bfU^t_\rP\right)_f.
\]
注意：
在求解速度方程后，求解器求解能量方程，并进行thermo.correct()函数，对\(\psi=1/RT\)进行更新。因此，在下文中，用\(\psi^{corr1}\)表示更新后的\(1/RT\)。
方程(23)在收敛的情况下可以写为：
(27)\[\begin{split}
\mathbf{U}_\mathrm{P}^{t+\dt}
\hspace{40em}
\\
= \mathbf{HbyA}^{t+\dt}_\mathrm{P} - \frac{1}{{{A^{t+\dt}_\mathrm{P}}}} \frac{1}{V_\rP}
\left(
\sum_f p_{\mathrm{rgh},f}^{t+\dt}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
+\sum_f \bfg_f\cdot \bfh_f \rho_f^{t+\dt}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right).
\hspace{8em}
\end{split}\]
在这里需要做一系列的假定：
略去临点的影响；
假定\(A_\rP\)的变化较小（\(A^{t+\dt}_\mathrm{P}=A^{t}_\mathrm{P}\)）；
方程(27)与方程(23)相减有：
(28)\[
\mathbf{U}_\mathrm{P}^{'} = - \frac{1}{{{A^{t}_\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).
\]
方程(28)需要提供了速度修正与压力修正的一对一关系。但发现，其中还存在密度修正。如果顺着这个思路进行，压力方程将进一步出现未知变量密度。为了处理这个问题，我们可以将方程(28)中的密度修正忽略掉。这样，即存在\(\mathbf{U}_f^{**}\)与\(p^{*}\)的一对一关系，结合连续性方程，可以用来求解压力。这是一种方法，但处理起来比较激进，可能会导致变量增量过大，引起发散。
如果考虑密度修正的话，则在压力方程中，会存在一个未知的密度变量。同时需要注意的是，密度变量也可以从状态方程中得出。例如，在求解能量方程后，会更新可压缩性\(\psi^{*}\)，在更新可压缩性后，在\(\psi^{corr1}\)的基础上，可以通过状态方程更新密度中间量：
(29)\[
\rho^{corr1}=p^t/\psi^{corr1}
\]
在这里，可以认为\(\rho^{corr}\)与密度修正量\(\rho^{'}\)的关系满足：
(30)\[
\rho^{corr1}=\rho^{*}+\rho^{'}
\]
由于\(\rho^{corr1}\)已经通过方程(29)求出，因此压力方程仅仅存在未知量压力。依据此思想，将方程(28)与方程(23)加和有：
(31)\[
\mathbf{U}_\mathrm{P}^{**} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A^t_\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^{corr1}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right).
\]
(32)\[
\mathbf{U}_f^{**} = \mathbf{HbyA}^{*}_f - \frac{1}{{{A^t_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^{corr1}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right)_f.
\]
注意：
这里\(\mathbf{U}^{**} =\mathbf{U}^{*}+\bfU^{'}\)，而不是\(\mathbf{U}^{t+\dt} =\mathbf{U}^{*}+\bfU^{'}\)，这是因为已经做了一系列的假定。
在这里，低速流认为\(\mathbf{U}_f^{**}\)需要满足连续性方程：
(33)\[
\frac{\rho^{t+\dt}-\rho^t}{\dt}+\sum_f \rho_f^{t+\dt}\mathbf{U}_f^{t+\dt} \cdot \bfS_f=0.
\]
现在我们通过方程(33)来构造压力柏松方程。首先，方程(33)要对当先的速度\(\bfU^{**}\)进行实施，来构成一种限定性条件，即：
(34)\[
\frac{\rho^{t+\dt}-\rho^t}{\dt}+\sum_f \rho_f^{t+\dt}\mathbf{U}_f^{**} \cdot \bfS_f=0.
\]
再次，方程(34)中的对流项存在未知变量密度\(\rho^{t+\dt}\)，因此要做线性化：
(35)\[
\frac{\rho^{t+\dt}-\rho^t}{\dt}+\sum_f \rho_f^{corr1}\mathbf{U}_f^{**} \cdot \bfS_f=0.
\]
同样，方程(35)的时间项存在未知的变量\(\rho^{t+\dt}\)。我们需要将\(\rho^{t+\dt}\)通过待求变量\(p^{*}\)进行表示（调用状态方程）：
(36)\[
\frac{\rho^{t+\dt}-\rho^t}{\dt}
=\frac{\psi^{corr1}p^{*}_\mathrm{rgh}-\rho^{corr1}_\mathrm{rgh}}{\dt}
+\frac{\rho^{corr1}-\rho^t}{\dt}
=\psi^{corr1}\frac{p^{*}_\mathrm{rgh}-p^{corr1}_\mathrm{rgh}}{\dt}
+\frac{\rho^{corr1}-\rho^t}{\dt}
\]
其中的\(p^{corr1}_\mathrm{rgh}\)表示当前已知的压力，在压力方程求解之前，可以用当前时间步已知的来替代：
(37)\[
\psi^{corr1}\frac{p^{*}_\mathrm{rgh}-p^{corr1}_\mathrm{rgh}}{\dt}
+\frac{\rho^{corr1}-\rho^t}{\dt}
=
\psi^{corr1}\frac{p^{*}_\mathrm{rgh}-p^{t}_\mathrm{rgh}}{\dt}
+\frac{\rho^{corr1}-\rho^t}{\dt}
\]
因此，方程(33)可以写为：
(38)\[
\psi^{corr1}\frac{p^{*}_\mathrm{rgh}-p^t_\mathrm{rgh}}{\dt}
+\frac{\rho^{corr1}-\rho^t}{\dt}
+\sum_f \rho_f^{corr1}\mathbf{U}_f^{**} \cdot \bfS_f=0.
\]
将方程(32)代入到方程(38)中有：
(39)\[\begin{split}
\psi^{corr1}\frac{p^{*}_\mathrm{rgh}-p^t_\mathrm{rgh}}{\dt}
+\frac{\rho^{corr1}-\rho^t}{\dt}
+\sum_f \rho_f^{corr} \mathbf{HbyA}^{*}_f \cdot \bfS_f
\hspace{15em}
\\
-\sum_f
\frac{\rho^{corr1}_f}{{{A^t_f}}} \frac{1}{V_\rP}
\left(
\sum_f \bfg_f\cdot \bfh_f \rho_f^{corr1}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right)_f
\cdot \bfS_f
\hspace{7em}
\\
= \sum_f\frac{\rho^{corr1}_f}{{{A^t_f}}} \frac{1}{V_\rP}
\left(
\sum_f p_{\mathrm{rgh},f}^{*}\frac{\bfS_f}{|\bfS_f|}|\bfS_f|
\right)_f
\cdot \bfS_f.
\hspace{3em}
\end{split}\]
注意：
这里\(\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方程笔记。
可见，方程(39)可用于求解获得压力\(p_\mathrm{rgh}^*\)。将\(p_\mathrm{rgh}^*\)代入到方程(32)有\(\mathbf{U}_\mathrm{P}^{**}\)。通过更新后的压力\(p_{\mathrm{rgh},f}^{*}\)以及速度\(\mathbf{U}_\mathrm{P}^{**}\)，依据密度更新的两种方法，进行下述两步：
依据状态方程更新密度。有压力修正量为：
(40)\[
p_{\mathrm{rgh},f}^{'}=p_{\mathrm{rgh},f}^{*}-p_{\mathrm{rgh},f}^{t}
\]
基于压力修正量，可以进一步更新密度：
(41)\[
\rho^{corr2}=\rho^{corr1}+\psi^{*}(p_{\mathrm{rgh},f}^{*}-p_{\mathrm{rgh},f}^{t})
\]
依据连续性方程更新密度。通过速度\(\mathbf{U}_\mathrm{P}^{**}\)组建新的通量，求解连续性方程。
(42)\[
\frac{\rho^{**}-\rho^t}{\Delta t}+\sum_f \rho^{corr1}_f\bfU_f^{**}\cdot\bfS_f=0,
\]
上述两种更新密度的方法，在一个时间步下收敛后，必然一致。因此，可以通过\(\rho^{**}\)与\(\rho^{corr2}\)的差来判断当前时间步下是否收敛：
(43)\[
\sum (\rho^{**}_\rP-\rho^{corr2}_\rP)V_\rP =? 0
\]
至此，我们有了满足连续性方程的速度场\(\mathbf{U}_\mathrm{P}^{**}\)。但由于在上述推导过程中，引入了大量的假设，\(\mathbf{U}_\mathrm{P}^{**}\)与\(p^{*}\)并不是严格满足动量方程。同时，方程(43)也未必满足。在这里可以把求得的这些场当作初始场，在一个时间步内再进行迭代求解。直到方程(43)满足为止。
总之，buoyantPimpleFoam中的的迭代过程可以表示为下面几个步骤：
通过方程(6)求解获得\(\rho^*\)；
通过\(\rho^*\)，组建速度方程，获得\(\bfU^*\)；
组建能量方程，更新温度等变量，更新可压缩性\(\psi^{corr}\)，通过状态方程更新\(\rho^{corr1}\)；
通过第二步的\(\bfU^*\)组建\(\bfHbyA^*\)，构建压力柏松方程(39)，求解获得\(p_{\mathrm{rgh}}^*\)；
通过\(p_{\mathrm{rgh}}^*\)更新速度\(\bfU^{**}\)，求解方程(42)，获得\(\rho^{**}\)；
通过方程(41)，获得\(\rho^{corr2}\)；
通过方程(43)判断连续性误差；
回到第二步进行迭代，直至通过密度方程求解的\(\rho^{**}\)与状态方程更新的密度\(\rho^{corr2}\)趋向于一致为止；
迭代为止之后，令通过密度方程求解的\(\rho\)与状态方程更新的密度\(\rho^{corrn}\)完全相等；