版本对应:

OpenFOAM-11中的fluid模块,对应OpenFOAM-10之前的rhoSimpleFoam、rhoPimpleFoam、buoyantFoam。其可以处理可压缩流体以及浮力驱动流。fluid模块主要基于isothermalFluid模块,仅仅在其基础上添加能量方程。在buoyancy模块开启后,可以处理浮力驱动流。

CFD:可压+体积力+稳态算法

可压稳态浮力驱动SIMPLE算法相对于可压稳态SIMPLE算法最重要的区别是增加了浮力项(也即重力\(\bfg\))。其他算法的处理大同小异。因此本文不会从基本的控制方程开始,而是仅仅分析浮力项的数值处理。在学习本文算法之前,请首先学习可压稳态SIMPLE算法

控制方程与算法

本文讨论的控制方程,与可压稳态SIMPLE算法中罗列的控制方程最重要的区别在于重力项:

(1)\[ \nabla\cdot\rho\bfU=0, \]
(2)\[ \nabla \cdot (\rho\mathbf{U} \mathbf{U})=-\nabla p+\nabla \cdot\tau +\rho \bfg \label{mom} \]
(3)\[ p=\rho RT. \]

方程(2)中的重力项以及压力项可以进行数值处理。定义

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

将方程(6)代入到方程(2)中有:

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

低速压力方程

可压稳态SIMPLE算法中,速度修正量与压力修正量的关系可以写为:

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

在浮力驱动流中,结合方程(6),方程(8)演变为:

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

方程(9)提供了速度修正与压力修正的一对一关系。但发现,其中还存在密度修正。如果顺着这个思路进行,压力方程将进一步出现未知变量密度。

为了处理这个问题,我们可以将方程(9)中的密度修正忽略掉。这是一种方法,但处理起来比较激进。另外一种方法是,参考可压稳态SIMPLE算法,由于能量方程会更新可压缩性\(\psi^{corr}\),在更新可压缩性后,在\(\psi^{corr}\)的基础上,可以通过状态方程更新密度中间量:

(10)\[ \rho^{corr}=p^n/\psi^{corr} \]

从而,认为\(\rho^n\)与密度修正量\(\rho^{'}\)的关系,满足:

(11)\[ \rho^{corr}=\rho^n+\rho^{'} \]

而非\(\rho^{n+1}=\rho^n+\rho^{'}\)。这样,由于\(\rho^{corr}\)已经通过方程(10)求出,因此压力方程仅仅存在未知量压力。依据此思想,结合修正方程(9)有:

(12)\[ \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). \]
(13)\[ \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. \]
(14)\[ \mathbf{U}^{**} =\mathbf{U}^{*}+\bfU^{'} \]

在这里,低速流认为\(\mathbf{U}_f^{**}\)需要满足连续性方程:

(15)\[ \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'\)处理。

同样的,方程(15)中的密度采用已知的密度(低速流密度变动充分小)。将方程(13)代入到方程(15)中有:

(16)\[\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)的连续形式为:

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

注意:

可以看出其中的未修正的速度,方程的左侧,可以继续写为\(\mathbf{HbyA}-\frac{\bfg\cdot\bfh}{A}\nabla\rho\)

可见,方程(16)可用于求解获得压力\(p^*\)。将\(p^{*}\)代入到方程(12)\(\mathbf{U}_\mathrm{P}^{**}\)。随后,密度再次通过状态方程进行更新:

(18)\[ \rho^*=\psi^{corr}p^* \]

在这里需要注意的是,参考可压稳态SIMPLE算法,其中的速度更新方法为:

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

其中并没有考虑体积力浮力的作用。在浮力驱动流中,要考虑浮力体积力的影响,也即需要通过方程(12)的方式来更新\(\mathbf{U}_\mathrm{P}^{**}\)。在这里将其进行重写:

(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}^{*}\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) \]

可压稳态浮力驱动SIMPLE算法为一个稳态求解器。同样的需要对压力进行场松弛,速度、能量等方程进行矩阵松弛增加稳定性。同样的,可压稳态浮力驱动流也可以调用SIMPLEC算法。由于算法考虑了浮力源项的作用,这主要带来了两点变化:

  • 附加浮力的求解器压力变量变为了\(p_{\mathrm{rgh}}\)

  • 附加除了压力项之外的求解器,速度修正与压力修正的关系会变得更加复杂,如方程(9)

  • 为了防止修正量过大,密度修正、压力修正都可以进行场松弛;

需要注意的是。可压稳态浮力驱动流SIMPLE算法在整个迭代过程中,为一个迭代求解器,因此对于各种变量的更新,某些过程可以忽略。一些可以忽略的过程比如有:

  • 速度方程可以不进行求解;

  • 密度中间变量可以不进行更新;

那么具体的,是否需要进行这些必要操作,是问题依赖的。但这并不影响最终向收敛过程的推进。

关键代码

可压稳态浮力驱动SIMPLE算法的关键代码与可压稳态SIMPLE算法无主要差异。唯一的区别在与压力方程。另外,方程(12)中速度方程的更新采用了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^{corr}_f}{A^n_\rP}\right)_f \bfg_f\cdot\bfh_f (\nabla\rho^{corr})_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^{corr}_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^{corr}}{A^n_\rP}\right)_f \bfg_f\cdot\bfh_f (\nabla\rho^{corr})_f\frac{\bfS_f}{|\bfS_f|}|\bfS_f| - \left(\frac{\rho^{corr}}{A^n_\rP}\right)_f (\nabla p_{\mathrm{rgh}})_f \frac{\bfS_f}{|\bfS_f|}|\bfS_f| \right)\frac{(A^n_\rP)_f}{\rho^{corr}_f} \right) \]

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

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

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