可压缩求解器大多使用非稳态算法。得益于双曲方程的数学特征,这些方法使用密度作为主要变量,然后通过状态方程求解压力,即密度基求解器。但密度基求解器在应用于不可压缩流的情况下效率低下。一方面的解释是在不可压缩领域,压力与密度的耦合非常弱,另一方面的解释是不可压缩假定下音速趋向于无穷大,导致时间步长过小。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算法中的的迭代过程可以表示为下面几个步骤:
通过方程(14)求解获得\(\bfU^*\),并组建\(\bfHbyA^*\);
通过能量方程更新可压缩性\(\psi^{corr}\);
低速流通过方程(24)求解获得\(p^*\),高速流通过方程(33)求解获得\(p^*\);
通过方程(34)求解获得\(\bfU^{**}\);
通过方程(35)更新获得\(\rho^{*}\);
回到第一步,继续进行求解;
稳态求解器需要对压力进行场松弛,速度、能量等方程进行矩阵松弛增加稳定性。同样的,isoThermalFluid模块也可以调用SIMPLEC算法。具体内容请参考不可压稳态SIMPLE算法。
同时需要注意的是,高速流与低速流的压力方程特征发生了变化,对比低速流的方程(24)与高速流的方程(33),低速流为压力的拉普拉斯方程。高速流的压力方程含有压力的对流项。