rhoSimpleFoam解析


李东岳


1. 引言

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

2. 控制方程与算法
2.1. 速度方程
rhoSimpleFoam中求解的为下述稳态质量方程、动量方程、以及状态方程: \begin{equation} \nabla\cdot\rho\bfU=0, \label{C} \end{equation} \begin{equation} \nabla \cdot (\rho\mathbf{U} \mathbf{U})=-\nabla p+\nabla \cdot\tau, \label{mom} \end{equation} \begin{equation} p=\rho RT. \label{eos} \end{equation}

其中方程\eqref{C}没有时间项,并无明显的变量,更像是一种限定性条件。因此稳态求解器中不存在对质量方程的求解。同时,方程\eqref{C}中的$\rho\bfU$可用于组建质量通量。 同时,rhoSimpleFoam中能量方程在求解能量变量后,对其直接求解即可,不涉及到速度-压力-密度耦合问题,此处略。上述三个方程,可用于迭代求解三个未知量:速度、压力、密度。首先对方程\eqref{mom}通过高斯定理进行对速度$\bfU$的离散,组建速度方程有:

\begin{equation} \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}^{n+1}}} \right)}_f}} \cdot\bfS_f = \sum_f F_f^n \mathbf{U}_f^{n+1} , \label{mom2} \end{equation}
\begin{equation} \int \nabla p \mathrm{d}V =\int p \mathrm{d}\bfS =\sum_f p_f^n\bfS_f, \label{gradp} \end{equation}

其中上标$^n$表示为当前迭代步(已知),上标$^{n+1}$表示下一个迭代步(待求),下标$_f$表示网格单元面上的值,$\bfS_f$表示网格单元的各个面的面矢量,$F_f$为质量通量,$\mu$为运动粘度(假定粘度为常数)。此处略去粘性项$\tau$的离散且假定粘度为0。整理方程\eqref{mom2}、\eqref{gradp}有:

\begin{equation} \sum_f F_f^n \mathbf{U}_f^{n+1} = -\sum_f p_f^n\bfS_f. \label{momF} \end{equation}

需要注意的是,方程\eqref{mom}中的左边第二项(对流项)是非线性的。在求解的过程中,要么选择非线性求解器,要么将对流项线性化。由于非线性求解器非常复杂,因此在OpenFOAM均采用线性化处理。具体的,在对方程\eqref{mom2}中的左边第二项对流项离散中,其中的通量$F_f^n$采用当前已知时间步的速度$\mathbf{U}^n$来计算,同时保留速度$\mathbf{U}^{n+1}_\mathrm{P}$为未知量。这种将高阶非线性项降为一阶线性项的过程即为线性化操作。在计算中$\bfU_f$需要从体心速度进行插值来获得,在此步可以引入各种插值格式。假设使用中心线性格式:

\begin{equation} \mathbf{U}_f^{n+1} = \frac{{\mathbf{U}_\mathrm{P}^{n+1} + \mathbf{U}_\mathrm{N}^{n+1}}}{2}. \label{linear} \end{equation}
\begin{equation} p_f^n = \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}. \label{linear2} \end{equation}

将方程\eqref{linear},\eqref{linear2}代入到方程\eqref{momF}有

\begin{equation} \sum_f {F_f^n \frac{{\mathbf{U}_\mathrm{P}^{n+1} + \mathbf{U}_\mathrm{N}^{n+1}}}{2}} = -\sum_f \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}\bfS_f, \end{equation}
\begin{equation} \sum_f { {{\frac{{F_f^n}}{2}} } } \mathbf{U}_\mathrm{P}^{n+1} = - \sum_f { { {\frac{{F_f^n}}{2} } \mathbf{U}_\mathrm{N}^{n+1}} } -\sum \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}\bfS_f. \end{equation}

将上式简写为

\begin{equation} {A_\mathrm{P}^n}\mathbf{U}_\mathrm{P}^{n+1}{\rm{ + }}\sum_f {A_\mathrm{N}^n\mathbf{U}_\mathrm{N}^{n+1}} = -\frac{1}{V_\rP}\sum_f \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}\bfS_f, \label{apanmom} \end{equation}

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

\begin{equation} A_\mathrm{P}^n= \frac{1}{V_\rP}\sum_f \frac{F_f^n}{2} , \end{equation}
\begin{equation} A_\mathrm{N}^n= \frac{1}{V_\rP}{\frac{{F_f^n}}{2}} , \end{equation}

求解方程\eqref{apanmom}即可获得速度$\mathbf{U}^{n+1}$。需要注意的是,当前密度、压力为$\rho^n$与$p^n$,二者并未更新。同时,对方程\eqref{apanmom}进行转换有:

\begin{equation}\label{un+1} \mathbf{U}_\mathrm{P}^{n+1} = \mathbf{HbyA}^{n+1}_\mathrm{P} - \frac{1}{{{A_\mathrm{P}^n}}} \frac{1}{V_\rP}\sum_f p_f^n\bfS_f. \end{equation}

其中$\bfHbyA$定义为:

\begin{equation}\label{hbya} \bfHbyA_\rP^{n+1}=\frac{-\sum_NA_\rN^n\bfU_\rN^{n+1}}{A_\rP^n}. \end{equation}

参考Rhie & Chow插值,有:

\begin{equation}\label{un+1f} \mathbf{U}_f^{n+1} = \mathbf{HbyA}^{n+1}_f - \frac{1}{{{A_f^n}}} \left(\frac{1}{V_\rP}\sum_f p_f^n\bfS_f\right)_f, \end{equation}
\begin{equation}\label{hbyA1f} \bfHbyA_f^{n+1}=\left(\frac{-\sum_NA_\rN^n\bfU_\rN^{n+1}}{A_\rP^n}\right)_f. \end{equation}
2.2. 压力泊松方程

N-S方程求解的关键问题之一在于并没有压力的方程出现。同时,方程\eqref{apanmom}求解的速度$\mathbf{U}^{n+1}_\mathrm{P}$并不满足连续性方程。在下一个迭代步,速度、密度、以及压力需要满足连续性方程,因此相应的量需要进行修正。对于网格单元$_\rP$的连续性方程$\nabla\cdot\rho\bfU=0$进行离散后的形式为 (Demirdžić et al., 1993)

\begin{equation} \sum_f (\rho_f^n+\rho_f')(\mathbf{U}_f^{n+1}+\bfU_f') \cdot \bfS_f=0. \label{cont} \end{equation}

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

\begin{equation} \rho_f'=\frac{1}{RT}p_f' \label{rhotil} \end{equation}

方程\eqref{cont}可以写为

\begin{equation} \sum_f \left(\rho_f^n+\frac{1}{RT}p_f'\right)(\mathbf{U}_f^{n+1}+\bfU_f') \cdot \bfS_f=0. \label{cont2} \end{equation}

展开为

\begin{equation} \sum_f \left(\rho_f^n\bfU_f^{n+1} + \rho_f^n\bfU_f' + \frac{1}{RT}p_f'\bfU_f^{n+1} \right)\cdot\bfS_f=0. \label{cont3} \end{equation}

其中$p'\bfU'$被省略。考虑方程\eqref{un+1},有:

\begin{equation}\label{un+21} \mathbf{U}_\mathrm{P}^{n+2} = \mathbf{HbyA}^{n+2}_\mathrm{P} - \frac{1}{V_\rP}\frac{1}{{{A_\mathrm{P}^{n+2}}}} \sum_f p_f^{n+1}\bfS_f. \end{equation}

方程\eqref{un+21}并不能够用来组建单变量压力方程($\bfHbyA^{n+2}$与$A_\mathrm{P}^{n+2}$未知)。因此在这里引入略去邻点影响的假定(用当前的$\bfHbyA^{n+1}$代替$\bfHbyA^{n+2}$),同时,假定密度依赖变量$A_\rP$不变,有:

\begin{equation}\label{un+2} \mathbf{U}_\mathrm{P}^{n+2} = \mathbf{HbyA}^{n+1}_\mathrm{P} - \frac{1}{V_\rP}\frac{1}{{{A_\mathrm{P}^{n}}}} \sum_f p_f^{n+1}\bfS_f. \end{equation}

方程\eqref{un+2}与\eqref{un+1}相减有:

\begin{equation}\label{uF} \bfU_\rP'=-\frac{1}{V_\rP}\frac{1}{{{A_\mathrm{P}^{n}}}} \sum_f p_f’\bfS_f \end{equation}

参考Rhie & Chow插值,面上的修正速度为

\begin{equation}\label{uFf} \bfU_f'=-\frac{1}{V_\rP}\frac{1}{{{A_f^{n}}}} \left(\sum_f p_f’\bfS_f\right)_f \end{equation}

将方程\eqref{uFf}、\eqref{un+1f}代入到\eqref{cont3}有:

\begin{equation} \sum_f \left(\rho_f^n \left(\mathbf{HbyA}^{n+1}_f - \frac{1}{V_\rP}\frac{1}{{{A_f^n}}} \left(\sum_f p_f^n\bfS_f\right)_f\right) - \rho_f^n \frac{1}{V_\rP}\frac{1}{{{A_f^{n}}}} \left(\sum_f p_f’\bfS_f\right)_f + \frac{1}{RT}p_f'\bfU_f^{n+1} \right)\cdot\bfS_f=0. \label{cont4} \end{equation}

引入修正压力:

\begin{equation} p_f^{n+1}=p_f^n+p_f' \end{equation}

整理方程\eqref{cont4}有

\begin{equation} \sum_f \left(\frac{\rho_f^n}{{{A_f^n}}} \left(\frac{1}{V_\rP}\sum_f p_f^{n+1}\bfS_f\right)_f \right)\cdot\bfS_f = \sum_f \rho_f^n \mathbf{HbyA}^{n+1}_f\cdot\bfS_f + \sum_f \frac{1}{RT}p_f'\bfU_f^{n+1}\cdot\bfS_f \label{press} \end{equation}

其中

\begin{equation} \frac{1}{RT}=\frac{\gamma}{c^2}=\frac{\p\rho}{\p p}. \label{RT} \end{equation}

其中$c$为音速。对于马赫数较小的流动,$c$趋向于无穷大。$1/RT$可以省略,即:

\begin{equation} \sum_f \left(\frac{\rho_f^n}{{{A_f^n}}} \left(\frac{1}{V_\rP}\sum_f p_f^{n+1}\bfS_f\right)_f \right)\cdot\bfS_f = \sum_f \rho_f^n \mathbf{HbyA}^{n+1}_f\cdot\bfS_f. \label{pressIncom} \end{equation}

有连续形式:

\begin{equation} \nabla\cdot\left(\frac{\rho^n}{A^n}\nabla p^{n+1}\right)=\nabla\cdot\rho^n\bfHbyA^{n+1} \label{pressIncom2} \end{equation}

求解方程\eqref{pressIncom2}有马赫数较小的$p^{n+1}$对于马赫数较大的流动,$1/RT$不能省略,继续将方程\eqref{un+1f}代入到\eqref{press}的右侧:

\begin{multline}\label{pressC1} \sum_f \rho_f^n \mathbf{HbyA}^{n+1}_f\cdot\bfS_f + \sum_f \frac{1}{RT}p_f'\bfU_f^{n+1}\cdot\bfS_f \\\\ = \sum_f \rho_f^n \mathbf{HbyA}^{n+1}_f\cdot\bfS_f + \sum_f \frac{p_f'}{RT}\mathbf{HbyA}^{n+1}_f \cdot\bfS_f + \sum_f \frac{p_f'}{RT}\left( - \frac{1}{{{A_f^n}}} \left(\frac{1}{V_\rP}\sum_f p_f^n\bfS_f\right)_f\right)\cdot\bfS_f \\\\ = \sum_f \frac{p^{n+1}}{RT} \mathbf{HbyA}^{n+1}_f\cdot\bfS_f + \sum_f \rho_f'\left( - \frac{1}{{{A_f^n}}} \left(\frac{1}{V_\rP}\sum_f p_f^n\bfS_f\right)_f\right)\cdot\bfS_f \end{multline}

结合方程\eqref{uFf},有

\begin{equation} \sum_f \rho_f'\left( - \frac{1}{{{A_f^n}}} \left(\frac{1}{V_\rP}\sum_f p_f^{n+1}\bfS_f\right)_f\right)\cdot\bfS_f = \sum_f \rho_f' \bfU_f'\cdot\bfS_f \label{presC2} \end{equation}

可见,方程\eqref{presC2}中由于$\rho_f' \bfU_f'$乘积项可以忽略。因此有连续形式

\begin{equation} \nabla\cdot \left(\frac{1}{RT}\bfHbyA^{n+1}p^{n+1}\right) = \nabla\cdot\left(\frac{\rho^n}{A^n}\nabla p^{n+1}\right) \label{pressIncom3} \end{equation}

在求得$p^{n+1}$之后,回代到方程\eqref{un+2}有$\bfU^{n+2}$。这里的$\bfU^{n+2}$未必满足连续性方程,同时由于方程\eqref{un+2}附加假定,需要再次进行速度压力耦合迭代求解。在迭代过程中,由于压力由$p^n$变为$p^{n+1}$,密度需要进行更新来重新计算$A_\rP$以及$A_\rN$。总而言之,可压缩SIMPLE算法中的的迭代过程可以表示为下面几个步骤:

  1. 依据初始条件,构建速度矩阵,依据当前密度,获得$A_\rP$以及$A_\rN$;
  2. 求解方程\eqref{apanmom}获得速度$\bfU^{n+1}$;
  3. 通过方程\eqref{hbya}组建$\mathbf{HbyA}^{n+1}$;
  4. 近音速条件下,求解方程\eqref{pressIncom3}获得压力$p^{n+1}$,低音速情况下,求解方程\eqref{pressIncom2}获得压力$p^{n+1}$;
  5. 通过方程\eqref{un+2}更新速度有$\mathbf{U}^{n+2}$;
  6. 通过状态方程更新密度;
  7. 回到第一步继续迭代几次;
3. 二维叶栅验证算例

rhoSimpleFoam自带若干算例。在OpenFOAM-7中,自带的算例主要为以下几个:

上述算例均可在OpenFOAM自带的tutorials里面找到,在此不做过多介绍。本算例与网格为CFD中文网用户上传的算例,算例针对rhoSimpleFoam做了一定的适配。算例网格为CFD中文网用户自行生成的2D纯六面体结构网格。右侧为出口,左侧为进口。具体如下:



constant文件夹的物性主要设置如下:

详细设置不便于在此处一一说明,可下载算例文件查看。用户可以直接输入

rhoSimpleFoam

进行求解。运算结果可参考下图。参考CFD中文网的相关讨论,本算例捕获到了叶片前缘的局部高速,更详细的结果可自行分析,研究湍流模型、壁面函数以及其他模型参数对结果的影响。

点击下载算例文件


4. 三段翼亚音速流动

本算例文件来自于CFD中文网的相关讨论。算例求解三段翼的亚音速流动($|\bfU|=125$ m/s)。用户进行设置的时候,遇到了严重的收敛问题。算例constant文件夹的设置参考之前的算例。进口以及边界条件设置如下:



由于可压缩求解器温度/压力的耦合作用,可压缩算例的边界条件需要细心的调节。同时:算例system文件夹下的fvOptions提供一定的稳定性限制。本算例对温度$T$进行上下界限制:

limitT
{
    type       limitTemperature;
    min        101;
    max        1000;
    selectionMode all;
}

详细设置不便于在此处一一说明,可下载算例文件查看。用户可以直接输入

rhoSimpleFoam

进行求解。运算结果可参考下图。初步计算表明湍流模型对算例结果的影响较大,用户可尝试使用不同的湍流模型进行计算。

点击下载算例文件



5. 亚音速轴对称射流

本算例文件来自于CFD中文网的相关讨论。算例对标NASA的提供的轴对称亚音速流动。出口流速大概在165m/s。用户进行设置的时候,遇到了1)稳态瞬态结果对不上,2)进口速度震荡问题。算例constant文件夹的设置参考之前的算例。湍流模型采用kOmegaSST。进口以及边界条件设置如下:



详细设置不便于在此处一一说明,可下载算例文件查看。用户可以直接输入

rhoSimpleFoam

进行求解。运算结果可参考下图。

点击下载算例文件



更新历史
2020.12.24Cuber Jeo | 2020.06.29增加算例 | 2020.06.02增加算例 | 2020.03.04创立页面

东岳流体 2014 - 2020
勘误、讨论、补充内容请前往CFD中文网