引言
bouyantBoussinesqPimpleFoam在OpenFOAM-7之后，就被bouyantPimpleFoam替换掉了。基金会将Boussinesq近似整理到了状态方程中。但在我们大量的侧重过程中，发现bouyantBoussinesqPimpleFoam在模拟大尺度流动的情况下，在计算开始的时候非常容易发散（可以参考这个算例）。因此本解析深入回顾老版本的bouyantBoussinesqPimpleFoam相关算法。
bouyantBoussinesqPimpleFoam求解器认为流场是不可压缩的，其中的浮力项\(\rho\bfg\)通过在动量方程中添加源项来处理。在传热领域内，Boussinesq假定认为在流动中温度的变化是非常小的，因此对应的密度变化也非常小。所以在控制方程中，除了浮力项，其他项的密度被认为是常数。
Boussinesq假定使得方程的非线性特性降低。但由于Boussinesq假定的限定条件，也使得调用Boussinesq假定的求解器存在一定的限制。在工程中，Boussinesq假定主要用于室温下的液体对流、建筑物对流、气相分散等。在温度变化比较大的情况下，不建议使用Boussinesq假定。Boussinesq假定认为流体的密度可以这样计算：
(1)\[
\rho=\rho_{ref}\rho_k, \;\;\rho_k=1-\beta\left(T-T_{ref}\right)
\]
其中\(\rho\)表示流体的密度，\(\rho_{ref}\)表示流体的参考密度，\(T\)表示流体的温度，\(T_{ref}\)表示流体的参考温度，\(\beta\)表示流体的体膨胀系数。
现考虑不可压缩流体，并附加重力源项的动量方程：
(2)\[
\frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot (\mathbf{U} \mathbf{U})=-\nabla p+\nabla \cdot\boldsymbol{\tau}+\frac{\rho}{\rho_{ref}}\bfg
\]
注意其中的压力项\(p\)的单位为压力除掉密度之后的单位。其中的浮力源项也可以写为\(\frac{\rho}{\rho_{ref}}\bfg=\rho_k\bfg\)。参考interFoam解析，其中的\(-\nabla p + \rho_k\bfg\)会进一步被处理为：
(3)\[
-\nabla p + \rho_k\bfg=-\nabla p_\mathrm{rgh}-\bfg\cdot \bfh \nabla \rho_k
\]
最终方程(2)会被整理为：
(4)\[
\frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot (\mathbf{U} \mathbf{U})=-\nabla p_\mathrm{rgh}-\bfg\cdot \bfh \nabla \rho_k +\nabla \cdot\boldsymbol{\tau}
\]
除此之外，温度也需要被传输，其控制方程为：
(5)\[
\frac{\p T}{\p t}+\nabla\cdot(\bfU T)-\nabla\cdot(\alpha_{eff}T)=0
\]
同理，连续性方程为：
(6)\[
\nabla\cdot\bfU=0
\]
方程(15)、(5)、(6)构成bouyantBoussinesqPimpleFoam的控制方程。
控制方程与算法
首先对方程(15)中的时间项进行对速度\(\bfU\)关于时间\(t\)的欧拉全隐离散有 ：
(7)\[
\int \int {\frac{{\partial \mathbf{U}}}{{\partial t}}\mathrm{d}V\mathrm{d}t = \left( {\mathbf{U}_\mathrm{P}^* - \mathbf{U}_\mathrm{P}^t} \right)\;} {V_\mathrm{P}}
\]
对流项隐性离散：
(8)\[
\int \int {\nabla \cdot \left( {\mathbf{U}\mathbf{U}} \right)\mathrm{d}V\mathrm{d}t = \int \int {\mathbf{U}\mathbf{U}\cdot\mathrm{d}\bfS\mathrm{d}t} = } \sum {{{\left( {{\mathbf{U}^*}{\mathbf{U}^t}} \right)}_f}} \cdot\bfS_f\Delta t = \Delta t\sum {{F_f^t \mathbf{U}_f^*} }
\]
其中上标\(^t\)表示为当前的时间步（已知），上标\(^*\)表示预测步（待求），下标\(_f\)表示网格单元面上的值，\(V_\rP\)表示网格单元体积，\(\bfS_f\)表示网格单元的各个面的面矢量，\(\Delta t\)表示时间步长，\(F_f^t\)为通量。\(\nabla \cdot\boldsymbol{\tau}\)项在这里不做分析，可以参考icoFoam解析中的扩散项的离散。
方程(15)中的压力项\(-\nabla p_\mathrm{rgh}\)与浮力项\(-\bfg\cdot \bfh \nabla \rho_k\)的离散有不同的方式。例如压力梯度的离散，可以采用fvc::grad()函数，也可以采用fvc::reconstruct()函数。采用前者可能会引起数值震荡。因此，OpenFOAM中的体积力均采用fvc::reconstruct()函数来进行。在这里请参考无痛苦NS方程笔记有详细的对比说明。
对动量方程离散后，即可组建动量方程矩阵，提取矩阵系数\(A_\rP\)，\(A_\rN\)等。这些过程与不可压缩流的处理方法一致，在这里不重复介绍，可以参考icoFoam解析
在组建压力方程的时候，参考icoFoam解析，首先组建\(\bfHbyA\)：
(9)\[
\mathbf{HbyA}^{t+\dt}_\mathrm{P} = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{t+\dt}} + S_\mathrm{P}^n} \right)
\]
可以看出，组建的\(\bfHbyA\)并没有考虑浮力的贡献。因此在后续组建压力方程的时候，需要考虑体积力（浮力）的附加效应。参考icoFoam解析，没有体积力的压力方程为：
(10)\[
\sum { {\mathbf{HbyA}_f^{t+\dt}} \cdot \bfS_f } = \sum { {\frac{1}{{A_{\mathrm{P},f}}}\left(\frac{1}{V_\rP}\sum p_f^{t+\dt}\bfS_f\right)_f} \cdot \bfS_f }
\]
在bouyantBoussinesqPimpleFoam中，首先其中的\(p\)要改写为\(p_{rgh}\)，其次体积力通量需要考虑进来：
(11)\[
\sum { \left(\mathbf{HbyA}_f^{t+\dt} - \frac{(\bfg\cdot\bfh)_f}{A_\rP}(\nabla\rho_k)_f\right)\cdot \bfS_f } = \sum { {\frac{1}{{A_{\mathrm{P},f}}}\left(\frac{1}{V_\rP}\sum p_{rgh,f}^{t+\dt}\bfS_f\right)_f} \cdot \bfS_f },
\]
再一次的，其中的体积力通常采用面法向梯度的形式来放置震荡，因此方程(11)中的体积力可以改写为
(12)\[
- \frac{(\bfg\cdot\bfh)_f}{A_\rP}(\nabla\rho_k)_f\cdot \bfS_f
=- \frac{(\bfg\cdot\bfh)_f}{A_\rP}\nabla(\rho_k)_f\cdot \frac{\bfS_f}{|\bfS_f|} \bfS_f
\]
那么方程(11)最终形式为：
(13)\[\begin{split}
\sum { \left(\mathbf{HbyA}_f^{t+\dt} \right)\cdot \bfS_f }
-\sum\left(\frac{(\bfg\cdot\bfh)_f}{A_\rP}\nabla(\rho_k)_f\cdot \frac{\bfS_f}{|\bfS_f|} \bfS_f\right)
\\
= \sum { {\frac{1}{{A_{\mathrm{P},f}}}\left(\frac{1}{V_\rP}\sum p_{rgh,f}^{t+\dt}\bfS_f\right)_f} \cdot \bfS_f },
\end{split}\]
方程(13)为bouyantBoussinesqPimpleFoam求解的压力方程。再一次的参考icoFoam解析，其中的速度更新方法为：
(14)\[
\mathbf{U}_\mathrm{P}^{**} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A_\mathrm{P}}}}\frac{1}{V_\rP}\sum p_f^*\bfS_f.
\]
其中并没有考虑体积力浮力的作用。因此，在bouyantBoussinesqPimpleFoam中，其应该写为：
(15)\[
\mathbf{U}_\mathrm{P}^{**} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A_\mathrm{P}}}}\frac{1}{V_\rP}\sum p_f^*\bfS_f
- \mathrm{reconstruct}\left(\frac{(\bfg\cdot\bfh)_f}{A_\rP}\nabla(\rho_k)_f\cdot \frac{\bfS_f}{|\bfS_f|} \bfS_f\right)
\]
除此之外，还应对温度方程进行离散。在更新温度后，需要对\(\rho_k\)进行更新。总体来说，该求解器有以下特点：