版本对应:
OpenFOAM-11中的fluid模块,对应OpenFOAM-10之前的rhoSimpleFoam、rhoPimpleFoam、buoyantFoam。其可以处理可压缩流体以及浮力驱动流。fluid模块主要基于isothermalFluid模块,仅仅在其基础上添加能量方程。在OpenFOAM-11之前,有无体积力采用不同的求解器来处理,buoyantSimpleFoam(有体积力),rhoSimpleFoam(无体积力)。在OpenFOAM-11都被处理到fluid模块中来减少代码复用。
因此本篇文章讨论的算法适用于OpenFOAM新版的fluid模块。以及rhoPimpleFoam、buoyantFoam(瞬态部分)、buoyantPimpleFoam。
CFD: 可压 + 瞬态
可压缩求解器相对于不可压缩求解器还可以继续进行区分:
是否需要考虑体积力(比如温度变化引起的浮力);
弱可压缩还是强可压缩(前者典型的为人居环境的温度变化,后者主要是超音速流动);
对于超音速流动,也即可压缩性特别强的情况,得益于双曲方程的数学特征,这些方法通常采用瞬态计算且使用密度作为主要变量,然后通过状态方程求解压力,即密度基求解器。但密度基求解器在应用于弱可压缩流的情况下效率低下。一方面的解释是在不可压缩领域,压力与密度的耦合非常弱,另一方面的解释是不可压缩假定下音速趋向于无穷大,导致时间步长过小。
针对弱可压缩领域,fluid模块为一个压力基、可压缩、适用于全流速、也包含体积力的求解器,主要用于求解普适性可压缩流动,也可用于求解低音速/超音速流动。其可以处理瞬态问题,也可以处理稳态问题。本文讨论瞬态算法。
控制方程
fluid模块(buoyantFoam,甚至更老的buoyantPimpleFoam)中求解的为下述瞬态连续性方程、动量方程、以及状态方程:
其实还需要附加能量方程。但能量方程不涉及到速度-压力-密度耦合问题,此处略。上述三个方程,可用于迭代求解三个未知量:速度、压力、密度。首先,方程(2)中的重力项以及压力项可以进行数值处理。定义
其中\(\bfh\)表示网格单元体心的位置矢量。对方程(4)进行梯度操作:
即为
方程(7)为求解的动量方程。
密度方程
不同与稳态连续性方程,方程(1)可以用于更新密度。因此方程(1)可以离散为:
如果对时间项采用Euler格式进行离散,空间项隐性格式有:
方程(9)存在两个未知量:密度\(\rho^{t+\dt}\)与速度\(\bfU^{t+\dt}\)。不可解。因此,在对对流项离散的过程中,需要使用当前时间步的速度\(\bfU^{t}\)。在这种情况下,离散方程可以写为:
因为对流项的离散并没有采用\(\bfU^{t+\dt}\)。
方程(10)可以用于求解密度。但求解的密度并非当前时间步最终收敛的密度\(\rho^{t+\dt}\)。
注意:
方程(1)若需要去掉时间项,就变成了稳态的连续型方程。稳态的连续性方程并无明显的变量,是一种限定性条件。因此稳态求解器不存在密度方程。
速度方程
对方程(7)通过高斯定理进行对速度\(\bfU\)的离散,组建速度方程有:
Warning
方程(11)右侧之所以\(\rho\)带星号,是因为密度方程已经求结果。在这里是一个已知量。
Warning
方程(12)\(\rho\)带星号与之前一样的结果。第一次迭代的时候\(\bfU\)采用当前\(t\)的值。因此\(F_f^t=\rho_f^*\bfU_f^t\cdot\bfS_f\)。
其中上标\(^t\)表示为当前时间步(已知),出来密度外的上标\(^{*}\)表示当前时间步下的迭代步(待求),下标\(_f\)表示网格单元面上的值,\(\bfS_f\)表示网格单元的各个面的面矢量,\(F_f\)为质量通量。此处略去粘性项\(\tau\)的离散因为其不影响压力速度耦合。
为了防止方程(13)与方程(14)的离散产生数值振荡。其一般采取面法向梯度的形式来进行:
需要注意的是,方程(17)左侧的对流项采用了线性化处理。方程(17)中的\(\bfU_f\)需要从体心速度进行插值来获得,在此步可以引入各种插值格式。假设使用中心线性格式(均一网格):
其可简写为
其中\(A_\rP^{t}\),\(A_\rN^{t}\)分别表示当前网格点与相邻网格点的离散系数:
可见,\(A_\mathrm{P}\)与\(A_\mathrm{N}\)为关于\(\rho^{*}\)的代数式(且当前为已知量)。方程(20)构成一个稀疏线性系统。求解方程(20)即可获得速度\(\mathbf{U}^{*}\)。
对方程(20)进行转换有:
其中\(\bfHbyA\)定义为:
同时,面插值有:
能量方程
能量方程的离散与速度方程的离散思路相同,在这里不做介绍。求解完动量方程之后,需要求解能量方程获得\(T^*\)。
对于基于可压缩性\(\psi^{*}\)的算法(例如hePsiThermo),在求解能量后,通过thermo.correct()
函数更新流体的可压缩性\(\psi^{*}\):
对于基于可压缩性\( \varrho ^{*}\)的算法(例如heRhoThermo),还需要更新热物理模型库中的密度,假定理想气体:
注意,这里的密度\( \varrho \)存储在thermo
中,而不是求解器的rho
。因此在迭代过程中,rho
不同于thermo.rho()
,\(\rho \neq \varrho\)。同时,不同的状态方程有不同的\( \varrho ^{*}\)更新方法。
压力方程
方程(23)在收敛的情况下可以写为:
在这里需要做一系列的假定:
略去临点的影响;
\(A^{t+\dt}_\mathrm{P}=A^{t}_\mathrm{P}\),否则其为一个非线性方程组不可解,也可以认为\(\dt\)间隔比较小的时候忽略其影响;
方程(30)需要提供了速度修正与压力修正的一对一关系。但发现针对浮力驱动流,其中还存在密度修正。如果顺着这个思路进行,压力方程将进一步出现未知变量密度导致无法封闭。
Warning
这也是可压缩求解器与不可压缩求解器的主要区别。可压缩求解器的速度压力耦合还需要考虑密度。
为了处理这个问题,我们可以将方程(30)中的密度修正忽略掉。这样,即存在\(\mathbf{U}_f^{**}\)与\(p^{*}\)的一对一关系,结合连续性方程,可以用来求解压力。这是一种方法,但处理起来比较激进,可能会导致变量增量过大,引起发散。
如果想考虑密度修正的话,则在压力方程中,会存在一个未知的密度变量导致不封闭。但同时需要注意的是,密度变量之前从状态方程中已经得出一个更新值\(\varrho^*\),目前并没有使用。因此在这里,可以认为\(\varrho^*\)与密度修正量\(\rho^{'}\)的关系满足(heRhoThermo类方法):
由于\(\varrho^*\)已求出,因此压力方程仅仅存在未知量压力。依据此思想,将方程(30)与方程(23)加和有:
注意:
这里\(\mathbf{U}^{**} =\mathbf{U}^{*}+\bfU^{'}\),而不是\(\mathbf{U}^{t+\dt} =\mathbf{U}^{*}+\bfU^{'}\),这是因为已经做了一系列的假定。同时,方程(33)中的密度已经被求出,因此方程(33)可以在后续中用来组建压力方程。
另外,如果采用hePsiThermo类方法,方程(33)中的\(\varrho_f^{*}\)直接采用状态方程求出。
首先考虑低速流的情况,低速流认为\(\mathbf{U}_f^{**}\)需要满足连续性方程:
现在我们通过方程(34)来构造压力柏松方程。首先,方程(34)存在多个未知变量如未知的密度,以及未知的速度,且不存在压力。因此要对其进行改动。首先要对当先的速度\(\bfU^{**}\)进行实施,来构成一种限定性条件,即:
再次,方程(35)中的对流项存在未知变量密度\(\rho^{t+\dt}\),因此要做线性化:
同样,方程(36)的时间项存在未知的变量\(\rho^{t+\dt}\)。问题是,用方程(36)代替方程方程(34)是否过于可行?
这是可以的,首先,PISO算法是需要进行若干次迭代的,也即\(\varrho^{*}\)就是需要进行多次内迭代,然后趋向于\(\rho^{t+\dt}\),\(\bfU^{**}\)经过多次内迭代,趋向于\(\bfU^{t+\dt}\)。假定PISO算法做了n次迭代,每次迭代更新到的密度为\(\tilde{\varrho}\):
方程(37)在第一次内迭代后变成
第二次内迭代后变成
第三次内迭代后变成
如果第三次内迭代后认为收敛,有
方程(37)实际就变为:
在最终,\(\rho^{t+\Delta t}-\varrho^{***} \)实际失去作用。因此,方程(36)可以改写为:
其中\(\rho^{**}\)表示获得\(\varrho^{*}\)当前内迭代步后下一个待求的内迭代步的密度。方程(39)的第一个时间项可以转换为压力的隐性离散项,第二个时间项为密度的显性离散,且每一次内迭代都会进行变化。因此其可以继续写为:
注意:
这里\(\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方程笔记。
可见,方程(41)可用于求解获得压力\(p_\mathrm{rgh}^*\)。将\(p_\mathrm{rgh}^*\)代入到方程(33)有\(\mathbf{U}_\mathrm{P}^{**}\)。
另外,对于\(\frac{\rho^{t+\dt}-\rho^t }{\dt}\),是否可以直接写成\(\psi^{*}\frac{ p^{*}_\mathrm{rgh}- p^t_\mathrm{rgh}}{\dt}\)(注意其中\(\psi^{*}\)已知,\( p^{*}_\mathrm{rgh}\)未知)。这样写也是可以的。在一些求解器中,确实这样操作。在这种情况下, 算例层面需要指定hePsiThermo。
连续性误差与迭代
在前文中已经讨论过,因为是瞬态,因此存在密度的传输方程,因此有3种方法来更新密度:
依据密度/压力修正来推进密度(heRhoThermo类算法):\(\varrho^{**}=\varrho^*+\psi^* p_\mathrm{rgh}^{'}\);
依据修正的可压缩性来推进密度(hePsiThermo类算法):\(\varrho^{**}=\psi^* p_\mathrm{rgh}^{*}\);
依据方程(10)更新的密度\(\rho^{**}\);
第1种与第2种的区别是heRhoThermo以及hePsiThermo的区别。除此之外,不管采用heRhoThermo还是hePsiThermo,这两种更新密度的方法,在一个时间步下收敛后,必然与第三种方法一致。因此,可以通过\(\varrho^{**}\)与\(\rho^{**}\)的差来判断当前时间步下是否收敛:
至此,我们有了满足连续性方程的速度场\(\mathbf{U}_\mathrm{P}^{**}\)。但由于在上述推导过程中,引入了大量的假设,\(\mathbf{U}_\mathrm{P}^{**}\)与\(p^{*}\)并不是严格满足动量方程。同时,方程(42)也未必满足。在这里可以把求得的这些场当作初始场,在一个时间步内再进行迭代求解。直到方程(42)满足为止。
总之,可压缩瞬态的迭代过程可以表示为下面几个步骤:
通过方程(10)求解获得\(\rho^*\);
通过\(\rho^*\),组建速度方程,获得\(\bfU^*\);
组建能量方程,更新温度等变量,更新可压缩性\(\psi^{*}\),通过状态方程更新\(\varrho^{*}\);
通过第二步的\(\bfU^*\)组建\(\bfHbyA^*\),构建压力柏松方程(41),求解获得\(p_{\mathrm{rgh}}^*\);
通过\(p_{\mathrm{rgh}}^*\)更新速度\(\bfU^{**}\),求解方程(10),获得\(\rho^{**}\);
通过方程(42)判断连续性误差;
回到第二步进行迭代,直至通过密度方程求解的\(\rho\)与状态方程更新的密度\(\varrho\)趋向于一致为止;
关键代码
针对上述求解流程,流程第一步,对应于代码中最开始的密度方程求解,获得\(\rho^*\):
#include "rhoEqn.H"
流程第二步,需要组建速度方程。b方程右侧,采用reconstruct()
函数保证无振荡:
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
上述代码的源项对应方程(15)、(16)。在速度方程更性之后,求解能量方程,进入流程第三步,并通过下述代码进行温度相关变量的更新:
thermo.correct();
以及密度的更新\(\varrho^*\):
rho = thermo.rho();
流程第四步,需要组建压力柏松方程:
fvScalarMatrix p_rghDDtEqn
(
psi*correction(fvm::ddt(p_rgh))
+ fvc::ddt(rho)
+ fvc::div(phi)
- fvm::laplacian(rhorAUf, p_rgh)
);
上述代码对应方程(41)。流程第五步,在求解压力之后,更新速度,更新密度:
phi += p_rghEqn.flux();
...
#include "rhoEqn.H"
流程第六步,更新\(\varrho^{**}\):
const volScalarField psip0(psi*p);
...
thermo.correctRho(psi*p - psip0);
流程第七步,判断连续性误差:
#include "compressibleContinuityErrs.H"
流程第八步,继续进行迭代。流程第九步,令通过密度方程求解的\(\rho\)与状态方程更新的密度\(\varrho\)完全相等:
rho = thermo.rho();