buoyantPimpleFoam解析

引言

buoyantPimpleFoam为一个瞬态的、附加浮力项的针对暖通领域的求解器。在查阅本资料前,首先请阅读buoyantSimpleFoam解析。buoyantSimpleFoam为对应的稳态版本。buoyantPimpleFoam在buoyantSimpleFoam基础上增加了瞬态功能。

控制方程与算法

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

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

  1. 略去临点的影响;

  2. 假定\(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中的的迭代过程可以表示为下面几个步骤:

  1. 通过方程(6)求解获得\(\rho^*\)

  2. 通过\(\rho^*\),组建速度方程,获得\(\bfU^*\)

  3. 组建能量方程,更新温度等变量,更新可压缩性\(\psi^{corr}\),通过状态方程更新\(\rho^{corr1}\)

  4. 通过第二步的\(\bfU^*\)组建\(\bfHbyA^*\),构建压力柏松方程(39),求解获得\(p_{\mathrm{rgh}}^*\)

  5. 通过\(p_{\mathrm{rgh}}^*\)更新速度\(\bfU^{**}\),求解方程(42),获得\(\rho^{**}\)

  6. 通过方程(41),获得\(\rho^{corr2}\)

  7. 通过方程(43)判断连续性误差;

  8. 回到第二步进行迭代,直至通过密度方程求解的\(\rho^{**}\)与状态方程更新的密度\(\rho^{corr2}\)趋向于一致为止;

  9. 迭代为止之后,令通过密度方程求解的\(\rho\)与状态方程更新的密度\(\rho^{corrn}\)完全相等;

关键代码

针对上述求解流程,流程第一步,对应于代码中最开始的密度方程求解:

#include "rhoEqn.H"

流程第二步,需要组建速度方程。buoyantPimpleFoam中的速度方程左侧并无特殊。但方程右侧,采用reconstruct()函数保证无振荡:

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

上述代码的源项对应方程(15)(16)。在速度方程更性之后,求解能量方程,进入流程第三步,并通过下述代码进行可压缩性的更新:

thermo.correct();

以及密度的更新:

rho = thermo.rho();

流程第四步,需要组建压力柏松方程:

fvScalarMatrix p_rghDDtEqn
(  
    psi*correction(fvm::ddt(p_rgh))
  + fvc::ddt(rho) 
  + fvc::div(phi)
  - fvm::laplacian(rhorAUf, p_rgh)
);

上述代码对应方程(39)。流程第五步,在求解压力之后,更新速度,更新密度:

phi += p_rghEqn.flux();
...
#include "rhoEqn.H"

流程第六步,更新\(\rho^{corr2}\)

thermo.rho() -= psi*p_rgh;
...
thermo.rho() += psi*p_rgh;

流程第七步,判断连续性误差:

#include "compressibleContinuityErrs.H"

流程第八步,继续进行迭代。流程第九步,令通过密度方程求解的\(\rho\)与状态方程更新的密度\(\rho^{corrn}\)完全相等:

rho = thermo.rho();