版本对应:

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

CFD:可压+稳态算法

可压缩求解器大多使用非稳态算法。得益于双曲方程的数学特征,这些方法使用密度作为主要变量,然后通过状态方程求解压力,即密度基求解器。但密度基求解器在应用于不可压缩流的情况下效率低下。一方面的解释是在不可压缩领域,压力与密度的耦合非常弱,另一方面的解释是不可压缩假定下音速趋向于无穷大,导致时间步长过小。rhoSimpleFoam为一个压力基、稳态的、无重力、可压缩、适用于全流速的求解器,主要用于求解普适性可压缩流动,也可用于求解低音速/超音速流动。在处理若可压缩流的情况下,类似的求解器为buoyantSimpleFoam,二者的主要区别在于后者考虑了浮力的作用,主要用于温度引起的浮力驱动流,并且后者只能处理亚音速流动。在应用于超音速的情况下,类似的求解器为rhoCentralFoam,后者使用中心迎风格式可以更尖锐的捕获激波不连续。相对于simpleFoam,rhoSimpleFoam需要处理压力-速度-密度三者的耦合问题。由于存在三变量耦合,因此,可压缩求解器对初始条件的选取更为苛刻。

控制方程与算法

首先有稳态质量方程、动量方程、以及状态方程:

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

方程(1)没有时间项,并无明显的变量,更像是一种限定性条件。因此稳态求解器中不存在对质量方程的求解。同时,方程(1)中的\(\rho\bfU\)可用于组建质量通量。上述三个方程,可用于迭代求解三个未知量:速度、压力、密度。首先对方程(2)通过高斯定理进行对速度\(\bfU\)的离散,组建速度方程有:

(4)\[ \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^{*} , \]
(5)\[ \int \nabla p \mathrm{d}V =\int p \mathrm{d}\bfS =\sum_f p_f^n\bfS_f, \]

其中上标\(^n\)表示为当前迭代步(已知),上标\(^{*}\)表示下一个迭代步(待求),下标\(_f\)表示网格单元面上的值,\(\bfS_f\)表示网格单元的各个面的面矢量,\(F_f\)为质量通量,\(\mu\)为运动粘度(假定粘度为常数)。此处略去粘性项\(\tau\)的离散且假定粘度为\(0\)。需要注意的是,目前并没有处理\(\nabla\cdot\tau\)这一项,因为其对求解流程并不影响。

注意:

OpenFOAM中可压缩求解器的通量phi定义为\(\rho_f\bfU_f\cdot\bfS_f\),不可压缩求解器的通量phi定义为\(\bfU_f\cdot\bfS_f\)

整理方程(4)(5)有:

(6)\[ \sum_f F_f^n \mathbf{U}_f^{*} = -\sum_f p_f^n\bfS_f. \]

需要注意的是,方程(6)左侧的对流项采用了线性化处理。方程(6)中的\(\bfU_f\)需要从体心速度进行插值来获得,在此步可以引入各种插值格式。假设使用中心线性格式(均一网格):

(7)\[ \mathbf{U}_f^{*} = \frac{{\mathbf{U}_\mathrm{P}^{*} + \mathbf{U}_\mathrm{N}^{*}}}{2}. \]
(8)\[ p_f^n = \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}. \]

将方程(7)(8)代入到方程(6)

(9)\[ \sum_f {F_f^n \frac{{\mathbf{U}_\mathrm{P}^{*} + \mathbf{U}_\mathrm{N}^{*}}}{2}} = -\sum_f \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}\bfS_f, \]
(10)\[ \sum_f { {{\frac{{F_f^n}}{2}} } } \mathbf{U}_\mathrm{P}^{*} = - \sum_f { { {\frac{{F_f^n}}{2} } \mathbf{U}_\mathrm{N}^{*}} } -\sum \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}\bfS_f. \]

将上式左右两边同时除以网格单元体积\(V_\rP\),同时简写为

(11)\[ {A_\mathrm{P}^n}\mathbf{U}_\mathrm{P}^{*}{\rm{ + }}\sum_f {A_\mathrm{N}^n\mathbf{U}_\mathrm{N}^{*}} = -\frac{1}{V_\rP}\sum_f \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}\bfS_f, \]

其中\(A_\rP^n\)\(A_\rN^n\)分别表示当前网格点与相邻网格点的离散系数:

(12)\[ A_\mathrm{P}^n= \frac{1}{V_\rP}\sum_f \frac{F_f^n}{2} , \]
(13)\[ A_\mathrm{N}^n= \frac{1}{V_\rP}{\frac{{F_f^n}}{2}} , \]

可见,\(A_\mathrm{P}^n\)\(A_\mathrm{N}^n\)为关于\(\rho^n\bfU^n\)的代数式。方程(11)构成一个稀疏线性系统。求解方程(11)即可获得速度\(\mathbf{U}^{*}\)。对方程(11)进行转换有:

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

其中\(\bfHbyA\)定义为:

(15)\[ \bfHbyA_\rP^{*}=\frac{-\sum_NA_\rN^n\bfU_\rN^{*}}{A_\rP^n}. \]

注意:

方程(15)中如果未进行动量预测,则右边取\(\bfU_\rN^n\)

同时,面上插值有:

(16)\[ \mathbf{U}_f^{*} = \mathbf{HbyA}^{*}_f - \frac{1}{{{A_f^n}}} \left(\frac{1}{V_\rP}\sum_f p_f^n\bfS_f\right)_f, \]
(17)\[ \bfHbyA_f^{*}=\left(\frac{-\sum_NA_\rN^n\bfU_\rN^{*}}{A_\rP^n}\right)_f. \]

注意:

在求解速度方程后,求解器求解能量方程,并进行thermo.correct()函数,对\(\psi=1/RT\)进行更新。因此,在下文中,用\(\psi^{corr}\)表示更新后的\(1/RT\)

低速流压力更新

本节我们考虑低速流动,即马赫数比较小的流动。在低速流中,可以认为密度的变化并不是特别的大。方程(14)在收敛的情况下可以写为:

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

方程(18)与方程(14)相减有(需将\(A\)线性化,类似\(A^{n+1}=A^n\)):

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

方程(19)提供了速度修正量与压力修正量的关系。在这里略去临点的影响,方程(19)可以写为:

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

方程(20)提供了速度修正与压力修正的一对一关系。将其与方程(14)加和有:

(21)\[ \mathbf{U}_\mathrm{P}^{**} = \mathbf{HbyA}^{*}_\mathrm{P} - \frac{1}{{{A^n_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{*}\bfS_f. \]
(22)\[ \mathbf{U}_f^{**} = \mathbf{HbyA}^{*}_f - \frac{1}{{{A^n_f}}} \frac{1}{V_\rP}\left(\sum_f p_f^{*}\bfS_f\right)_f. \]

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

(23)\[ \sum_f \rho_f^{n}\mathbf{U}_f^{**} \cdot \bfS_f=0. \]

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

(24)\[ \sum_f \rho_f^{n} \mathbf{HbyA}^{*}_f \cdot \bfS_f= \sum_f\frac{1}{{{A^n_f}}} \frac{1}{V_\rP}\left(\sum_f p_f^{*}\bfS_f\right)_f \cdot \bfS_f. \]

可见,方程(24)可用于求解获得压力\(p^*\)

高速流transonic压力更新

本节我们考虑高速流动,即马赫数比较大的流动。在高速流中,可以认为密度的变化比较明显。对于网格单元\(_\rP\)的连续性方程\(\nabla\cdot\rho\bfU=0\)进行离散后的形式为:

(25)\[ \sum_f (\rho_f^n+\rho_f')(\mathbf{U}_f^{*}+\bfU_f') \cdot \bfS_f=0. \]

其中\(\rho_f'\)表示修正密度,\(\bfU_f'\)表示修正速度。结合状态方程:

(26)\[ \rho_f'=\psi^{corr}p_f' \]

方程(25)可以写为

(27)\[ \sum_f \left(\rho_f^n+\psi^{corr}p_f'\right)(\mathbf{U}_f^{*}+\bfU_f') \cdot \bfS_f=0. \]

展开为

(28)\[ \sum_f \left(\left(\rho_f^n+\psi^{corr}p_f'\right)\bfU_f^{*} + \rho_f^n\bfU_f' \right)\cdot\bfS_f=0. \]
(29)\[ \sum_f \left(\left(\rho_f^n+\psi^{corr}p_f'\right)\bfU_f^{*} \right)\cdot\bfS_f = -\sum_f \left( \rho_f^n\bfU_f' \right)\cdot\bfS_f \]

注意:

方程(29)中的\(p'\bfU'\)为两个修正量的乘积,被省略。

因此在这里引入略去邻点影响的假定:

(30)\[ \mathbf{U}_\mathrm{P}^{'} = - \frac{1}{V_\rP}\frac{1}{{{A_\mathrm{P}^n}}} \sum_f p_f^{'}\bfS_f. \]
(31)\[ \bfU_f'=-\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f’\bfS_f\right)_f \]

将方程(31)(16)代入到(29)有:

\[\begin{split} \sum_f \left(\rho_f^n+\psi^{corr}p_f'\right) \mathbf{HbyA}^{*}_f \cdot\bfS_f - \sum_f \left(\rho_f^n+\psi^{corr}p_f'\right) \frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\cdot\bfS_f \hspace{6em} \\ = \sum_f \rho_f^n \frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{'}\bfS_f\right)_f \cdot\bfS_f\hspace{3em} \end{split}\]
\[\begin{split} \sum_f \left(\rho_f^n+\psi^{corr}p_f'\right) \mathbf{HbyA}^{*}_f \cdot\bfS_f - \sum_f \left(\psi^{corr}p_f'\right) \frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\cdot\bfS_f \hspace{6em} \\ = \sum_f \rho_f^n \frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{*}\bfS_f\right)_f \cdot\bfS_f \hspace{2em} \end{split}\]
\[\begin{split} \sum_f \rho_f^n \mathbf{HbyA}^{*}_f \cdot\bfS_f + \sum_f \psi^{corr}p_f' \left(\mathbf{HbyA}^{*}_f-\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\right)\cdot\bfS_f \hspace{6em} \\ = \sum_f \rho_f^n \frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{*}\bfS_f\right)_f \cdot\bfS_f \hspace{2em} \end{split}\]
\[\begin{split} \sum_f \rho_f^n \mathbf{HbyA}^{*}_f \cdot\bfS_f + \sum_f \psi^{corr}(p_f^*-p_f^n) \left(\mathbf{HbyA}^{*}_f-\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\right)\cdot\bfS_f \hspace{6em} \\ = \sum_f \rho_f^n \frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{*}\bfS_f\right)_f \cdot\bfS_f \hspace{2em} \end{split}\]
\[\begin{split} \sum_f \rho_f^n \mathbf{HbyA}^{*}_f \cdot\bfS_f - \sum_f \psi^{corr}p_f^n \left(\mathbf{HbyA}^{*}_f-\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\right)\cdot\bfS_f \hspace{9em} \\ + \sum_f \psi^{corr}p_f^* \left(\mathbf{HbyA}^{*}_f-\frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\right)\cdot\bfS_f = \sum_f \rho_f^n \frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{*}\bfS_f\right)_f \cdot\bfS_f \hspace{3em} \end{split}\]

方程最后的形式为:

(32)\[\begin{split} \sum_f \left(\rho_f^n-\psi^{corr}p_f^n\right) \mathbf{HbyA}^{*}_f \cdot\bfS_f + \sum_f \psi^{corr}p_f^*\mathbf{HbyA}^{*}_f\cdot\bfS_f \hspace{17em} \\ + \sum_f \psi^{corr}(p_f^n - p_f^*) \frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\cdot\bfS_f = \sum_f \rho_f^n \frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{*}\bfS_f\right)_f \cdot\bfS_f \hspace{3em} \end{split}\]

方程(32)左边的第三项,在rhoSimpleFoam中进行了省略,演变为:

(33)\[\begin{split} \sum_f \left(\rho_f^n-\psi^{corr}p_f^n\right) \mathbf{HbyA}^{*}_f \cdot\bfS_f + \sum_f \psi^{corr}p_f^*\mathbf{HbyA}^{*}_f\cdot\bfS_f \hspace{17em} \\ = \sum_f \rho_f^n \frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^{*}\bfS_f\right)_f \cdot\bfS_f \hspace{6em} \end{split}\]

方程(33)可以用来求解\(p^{*}\)

至此,低速流与高速流的\(p^*\)更新方法简述完毕。现将方程(14)与方程(30)加和有:

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

\(p^{*}\)代入到方程(34)\(\mathbf{U}_\mathrm{P}^{**}\)。随后,密度通过状态方程进行更新:

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

总之,可压缩SIMPLE算法中的的迭代过程可以表示为下面几个步骤:

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

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

  3. 低速流通过方程(24)求解获得\(p^*\),高速流通过方程(33)求解获得\(p^*\)

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

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

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

稳态求解器需要对压力进行场松弛,速度、能量等方程进行矩阵松弛增加稳定性。同样的,isoThermalFluid模块也可以调用SIMPLEC算法。具体内容请参考不可压稳态SIMPLE算法

同时需要注意的是,高速流与低速流的压力方程特征发生了变化,对比低速流的方程(24)与高速流的方程(33),低速流为压力的拉普拉斯方程。高速流的压力方程含有压力的对流项。

关键代码

isoThermalFluid模块中的速度方程并无特殊。在速度方程更性之后,fluid模块中会求解能量方程,并通过下述代码进行可压缩性的更新:

thermo.correct();

高速流transonic中的压力方程通过下面的代码组建:

fvScalarMatrix pEqn
(
    fvc::div(phiHbyA)
  + fvm::div(phid, p)
  - fvm::laplacian(rhorAAtUf, p)
  ==
    fvModels.source(psi, p, rho.name())
);

其中很明显具有压力的对流项。低速流中的压力方程通过下面的代码组建:

fvScalarMatrix pEqn
(
    fvc::div(phiHbyA)
  - fvm::laplacian(rhorAAtUf, p)
  ==
    fvModels.source(psi, p, rho.name())
);

其为一个拉普拉斯方程。在求解后,需要对密度进行更新:

rho = thermo.rho();