buoyantSimpleFoam解析

引言

在查阅本资料前,首先请阅读rhoSimpleFoam解析。这是因为rhoSimpleFoam与buoyantSimpleFoam的求解思路大体一致。区别在于buoyantSimpleFoam考虑了浮力的作用,主要用于模拟温度引起的浮力驱动流。rhoSimpleFoam并没有包含浮力项,其认为密度差异不能带来流体的流动,因此不能处理热对流问题。

控制方程与算法

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). \]

在这里需要做一系列的假定:

  1. 略去临点的影响;

  2. 假定\(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}\)已经通过方程(26)求出,因此压力方程仅仅存在未知量压力。依据此思想,将方程(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解析

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

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

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

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

总之,buoyantSimpleFoam中的的迭代过程可以表示为下面几个步骤:

  1. 通过方程(20)求解获得\(\bfU^*\),并组建\(\bfHbyA^*\)

  2. 通过能量方程更新可压缩性\(\psi^{corr}\)

  3. 通过状态方程更新中间变量密度\(\rho^{corr}\)

  4. 通过方程(32)求解获得\(p_{\mathrm{rgh}}^*\)

  5. 通过方程(29)求解获得\(\bfU^{**}\)

  6. 通过方程(33)更新获得\(\rho^{*}\)

  7. 回到第一步,继续进行求解;

在这里需要注意的是,buoyantSimpleFoam为一个稳态求解器。同样的需要对压力进行场松弛,速度、能量等方程进行矩阵松弛增加稳定性。同样的,buoyantSimpleFoam也可以调用SIMPLEC算法,但求解器并未进行植入。不同于rhoSimpleFoam、simpleFoam,buoyantSimpleFoam主要考虑了浮力源项的作用,这主要带来了两点变化:

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

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

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

需要注意的是。buoyantSimpleFoam在整个迭代过程中,为一个迭代求解器,因此对于各种变量的更新,某些过程可以忽略。一些可以忽略的过程比如有:

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

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

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

关键代码

buoyantSimpleFoam中的速度方程左侧并无特殊。但方程右侧,采用reconstruct()函数保证无振荡:

    UEqn
 ==
    fvc::reconstruct
    (
        (
          - ghf*fvc::snGrad(rho)
          - fvc::snGrad(p_rgh)
        )*mesh.magSf()
    )

其对应方程(11)(12)。在速度方程更性之后,求解能量方程,并通过下述代码进行可压缩性的更新:

thermo.correct();

随后,为了获得更新的\(\rho^{corr}\),在组建压力方程之前,更新密度:

rho = thermo.rho();

随后,将速度通量通过phiHbyA来定义,即方程(32)左边第一项:

surfaceScalarField phiHbyA
(
    "phiHbyA",
    fvc::interpolate(rho)*fvc::flux(HbyA)
);

方程(32)左边第二项表示浮力项,通过phig来表示:

surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());

随后,组建压力方程更新压力后,再次更新密度获得\(\rho^*\)

rho = thermo.rho();