rhoPimpleFoam解析


李东岳


1. 引言

可压缩流动的求解可分为不同的方法。在严格双曲的情况下,瞬态可压缩流动可利用双曲特征的优势采用显性密度基步进求解。在具备混合型特征的可压缩流动中,通常需要进行额外的稳定性适配。一些比较常见的方法有预条件法、压力速度耦合法等。rhoPimpleFoam为一个瞬态的普适性无重力可压缩求解器,主要为求解亚音速流动以及暖通类应用设计,但也可用于求解跨/超音速流动。另外类似的求解器为buoyantPimpleFoam,二者的主要区别在于后者考虑了浮力的作用,并且buoyantPimpleFoam只能处理亚音速流动,其主要用于温度变化引起的浮力驱动流。rhoPimpleFoam中默认存在粘度项,并且植入的能量方程进行了简化,这会导致方程特征的改变因而激波并不会很尖锐。方程的显性离散也并没有受益于高超音速流动的双曲特征。本文对rhoPimpleFoam的解析在icoFoam解析的基础上进行适配,建议首先阅读icoFoam解析获取更基础的算法信息。

2. 控制方程与算法

参考icoFoam解析中讨论的动量方程并添加密度项,以及考虑CFD中的能量方程讨论的能量方程有:

\begin{equation} \frac{\partial \rho}{\partial t}+\nabla \cdot \left(\rho \mathbf{U} \right)=0 \label{NS1} \end{equation} \begin{equation} \frac{\partial \rho \mathbf{U} }{\partial t}+\nabla \cdot \left(\rho \mathbf{U} \mathbf{U} \right)-\nabla \cdot \tau=-\nabla p \label{NS2} \end{equation} \begin{equation} \frac{\partial \rho h}{\partial t}+\nabla \cdot (\rho \mathbf{U} h) + \frac{\partial \rho K}{\partial t}+\nabla \cdot (\rho \mathbf{U} K) - \nabla \cdot (\alpha_\mathrm{eff}\nabla h) =\frac{\partial p}{\partial t} \label{hfinal} \end{equation}
同时附加状态方程 \begin{equation} \rho=\psi p \label{bar} \end{equation}

其中的$\mathbf{U}$为速度矢量,$p$为压力,$\rho$为密度,$\tau$为剪切应力,$h$为比焓,$K$表示机械能,$\alpha_{\mathrm{eff}}$表示有效导热系数,其等于$\frac{\nu_t}{\mathrm{Pr}}$,$\psi=1/RT$表示可压缩性。

2.1. 密度方程与速度方程

首先对方程\eqref{NS2}中的时间项进行对速度$\bfU$关于时间$t$的欧拉全隐离散有

\begin{equation} \int \int {\frac{{\partial\rho \mathbf{U}}}{{\partial t}}\mathrm{d}V\mathrm{d}t = \left( {\rho_\rP^{t+\Delta t}\mathbf{U}_\mathrm{P}^* -\rho_\rP^{t} \mathbf{U}_\mathrm{P}^t} \right)\;} {V_\mathrm{P}}, \label{mom1} \end{equation}

其中$\bfU_\rP^*$为待求速度,但方程\eqref{mom1}中继续存在未知量$\rho_\rP^{t+\dt}$,因此,在离散动量方程之前,需要对密度方程\eqref{NS1}进行更新获得$\rho_\rP^{t+\dt}$。但由于方程\eqref{NS1}中的速度为当前的速度,因此求解的密度为显性解$\rho_\rP^{t+\dt}$。随后,对流项隐性离散

\begin{equation} \int \int {\nabla \cdot \left( {\rho\mathbf{U}\mathbf{U}} \right)\mathrm{d}V\mathrm{d}t = \int \int {\rho\mathbf{U}\mathbf{U}\mathrm{d}\bfS\mathrm{d}t} = } \sum {{{\left( {{\rho^t\mathbf{U}^*}{\mathbf{U}^t}} \right)}_f}} \bfS_f\Delta t = \Delta t\sum {{F_f^t \mathbf{U}_f^*} } , \label{mom2} \end{equation}

其中的$F_f^t$表示质量通量,其等于$\rho_f\bfU_f\cdot\bfS_f$。方程\eqref{NS2}中的$\nabla\cdot\tau$通过湍流模型进行模化,在此不做离散处理。压力项显性离散

\begin{equation} -\int \int \nabla p \mathrm{d}V\mathrm{d}t=-\int\int p \mathrm{d}\bfS\mathrm{d}t=-\Delta t\sum \left(p_f^t\bfS_f\right), \label{gradp} \end{equation}

将方程\eqref{mom1}、\eqref{mom2}以及\eqref{gradp}整理有(此处未考虑粘性项):

\begin{equation} \left( {\rho_\rP^{t+\dt}\mathbf{U}_\mathrm{P}^* -\rho_\rP^t \mathbf{U}_\mathrm{P}^t} \right)\; {V_\mathrm{P}} + \Delta t\sum F_f^t \mathbf{U}_f^* =-\Delta t\sum \left(p_f^t\bfS_f\right) \label{predic} \end{equation}
方程\eqref{predic}中$\mathbf{U}_f^*$被定义为面上的预测速度。面速度需要从体心速度进行插值来获得,在此步可以引入各种插值格式。假设在均一网格上使用中心线性格式

\begin{equation} \mathbf{U}_f^* = \frac{\bfU_\rP^* + \bfU_\rN^*}{2}. \label{linear} \end{equation}
\begin{equation} p_f^t = \frac{p_\mathrm{P}^t + p_\mathrm{N}^t}{2}. \label{linear2} \end{equation}

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

\begin{equation} \left( {\rho_\rP^{t+\dt}\mathbf{U}_\mathrm{P}^* -\rho_\rP^t \mathbf{U}_\mathrm{P}^t} \right)\; {V_\mathrm{P}} + \Delta t\sum F_f^t \frac{\bfU_\rP^* + \bfU_\rN^*}{2} =-\Delta t\sum \left(\frac{p_\mathrm{P}^t + p_\mathrm{N}^t}{2}\bfS_f\right) \end{equation}
整理有:
\begin{equation} \left(\frac{\rho_\rP^{t+\dt}}{\Delta t}+\frac{\sum F_f^t}{2V_\rP} \right)\bfU_\rP^* + \frac{\sum F_f^t}{2V_\rP}\bfU_\rN^* = \frac{\rho_\rP^{t}}{\Delta t}\bfU_\rP^t-\frac{1}{V_\rP}\sum \left(\frac{p_\mathrm{P}^t + p_\mathrm{N}^t}{2}\bfS_f\right) \label{predic2} \end{equation}
将上述方程简写为
\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^*{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^*} = S_\mathrm{P}^t-\frac{1}{V_\rP}\sum \frac{{p_\mathrm{P}^t + p_\mathrm{N}^t}}{2}\bfS_f, \label{apanmom} \end{equation}
其中 \begin{equation}\label{Ap} A_\mathrm{P}=\frac{\rho_\rP^{t+\dt}}{\Delta t}+\frac{\sum F_f^t}{2V_\rP}, \end{equation} \begin{equation} A_\mathrm{N}= \frac{{F_f^t}}{2V_\rP} , \end{equation} \begin{equation} S_\mathrm{P}^t= \frac{\rho_\rP^{t}}{{\Delta t}}\mathbf{U}_\mathrm{P}^t. \end{equation}

在这里在此需要强调的是,上述方程中$\rho_\rP^{t+\dt}$为已知量,因为其已经通过密度方程获得显性解。求解方程\eqref{apanmom}即可获得速度$\bfU^*$。此处的速度并不满足质量守恒。

2.2. 压力泊松方程

压力基求解器需要利用连续性方程构建压力泊松方程,进而修正速度。在获得速度$\bfU^*$后,利用$A_\rP$与$A_\rN$可以组建速度$\bfHbyA^*$收敛量(这里并未考虑压力的作用):

\begin{equation} \bfHbyA_\rP^{t+\dt} = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{t+\dt}} + S_\rP^t} \right) \label{hbya2} \end{equation}
同时有:
\begin{equation} \bfU_\rP^{t+\dt}=\bfHbyA_\rP^{t+\dt} - \frac{1}{{{A_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{t+\dt}\bfS_f \label{Up2} \end{equation}
考虑连续性方程,其最终的离散形式为:
\begin{equation} \frac{\rho^{t+\dt}-\rho^t}{\dt}+\sum_f \rho_f^{t+\dt}\mathbf{U}_f^{t+\dt} \cdot \bfS_f=0. \label{cont} \end{equation}
将方程\eqref{Up2}进行面插值,同时代入至\eqref{cont}有:
\begin{equation} \frac{\rho^{t+\dt}-\rho^t}{\dt}+\sum_f \rho_f^{t+\dt}\left(\bfHbyA_\rP^{t+\dt} - \frac{1}{{{A_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{t+\dt}\bfS_f\right)_f \cdot \bfS_f=0. \label{cont2} \end{equation}
整理有:
\begin{equation} \frac{\rho^{t+\dt}-\rho^t}{\dt} + \sum_f \rho_f^{t+\dt}\bfHbyA_f^{t+\dt} \cdot \bfS_f = \sum_f \rho_f^{t+\dt}\left(\frac{1}{{{A_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{t+\dt}\bfS_f\right)_f \cdot \bfS_f. \label{pres} \end{equation}
即为亚音速情况下连续形式的压力泊松方程:
\begin{equation} \frac{\rho^{t+\dt}-\rho^t}{\dt}+\nabla \cdot \left(\rho\mathbf{HbyA}^{t+\dt} \right)=\nabla \cdot \left( \frac{\rho^{t+\dt}}{{{A}}}\nabla p^{t+\dt} \right) \label{p2pre} \end{equation}
2.3. 迭代算法

由于当前$\bfHbyA^{t+\dt}$以及$p^{t+\dt}$均为未知的,因此压力方程不能求解。下面就讨论一种迭代算法来处理这种问题。方程\eqref{p2pre}中出现多个未知量的原因主要在于方程\eqref{Up2}的代入。如果忽略邻点的影响,方程\eqref{Up2}中的速度可以改写为:

\begin{equation} \bfU_\rP^{**}=\bfHbyA_\rP^* - \frac{1}{{{A_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^*\bfS_f \label{Up21} \end{equation} 在这种情况下,将\eqref{Up21}代入到方程\eqref{cont},可组建下述压力泊松方程:
\begin{equation} \frac{\rho^{t+\dt}-\rho^t}{\dt}+\nabla \cdot \left(\rho\mathbf{HbyA}^{*} \right)=\nabla \cdot \left( \frac{\rho^{t+\dt}}{{{A}}}\nabla p^{*} \right) \label{p2pre2} \end{equation}

另外值得一提的是,在本文的算法中,可压缩求解器中的密度一个是通过方程\eqref{NS1}进行更新,一个是通过压力以及状态方程进行求解。俩种方式最终的密度应该相同。因此,在rhoPimpleFoam中,通过对比二者的密度差进行误差测试。具体的,在每一次压力迭代求解过程中,更新速度以及通量,求解密度方程获得通过求解连续性方程的密度解。另一方面,通过状态方程获得密度二者的差,即为连续性误差。总而言之,可压缩求解器中的迭代过程可以表示为下面几个步骤:

  1. 依据初始速度组建通量,显性离散方程\eqref{NS1},获得下一个时间步的密度$\rho^{t+\dt}$;
  2. 通过方程\eqref{apanmom},获得$\bfU^*$以及$\mathbf{HbyA}^*$;
  3. 求解方程\eqref{p2pre2},获得压力$p^*$;
  4. 通过方程\eqref{Up21}获得速度$\bfU^{**}$以及通量,在此通过方程\eqref{NS1}更新密度$\rho^{t+\dt}$;
  5. 依据状态方程,以及压力$p^*$,更新$\rho^*$,
  6. 通过$\rho^*$与$\rho^{t+\dt}$,求得连续性误差;
  7. 回到第一步迭代求解几次,在时间步内收敛的情况下,$\rho^* \rightarrow \rho^{t+\dt}$;

对比所有的可压缩求解器和不可压缩求解器,会发现不可压缩求解器的压力方程为一个纯粹的泊松方程。在可压缩求解器中,压力方程包含了对流项和时间项,其实际上是对流波动方程。可压缩求解器中压力方程中的对流项和扩散项的相对重要性和速度有关,在速度较低的情况下,拉普拉斯项占主要作用,可压缩压力方程近似于不可压缩的压力泊松方程。在速度较高的情况下,对流项占主导地位,其反映了流场本身的双曲特性。这样,可压缩求解器中的压力方程等于求解连续性方程(密度)。

3. 验证算例

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

上述算例均可在OpenFOAM自带的tutorials里面找到,在此不做过多介绍。本算例为东岳流体为采用NACA0012网格,测试rhoPimpleFoam捕获超音速激波的能力,并测试不同离散格式对结果的影响。

3.1. 模型设置

算例网格为采用blockMesh生成的2D纯六面体结构网格。右侧为出口,除出口外均为自由流边界条件。算例默认为无粘度层流。初始速度为600 $\mathrm{m}/\mathrm{s}$。

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

3.2. 算例运行

相关算例运行首先需要执行网格转换命令:

./Allrun

即可运行。用户可打开Allrun文件检查所运行的命令。

运算结果可参考下图(左图)。很明显,rhoPimpleFoam并没有捕获到尖锐的尖端。算例结果可以通过调整离散格式进一步的提高,用户可以尝试调整离散格式,来获得右图的更好的间断解。然而,在右图中,rhoPimpleFoam捕获的间断虽有提高但仍不明晰。用户可以尝试对能量方程进行适配、或采用密度基求解器进行求解来获得更好的间断解。


采用rhoPimpleFoam,并结合不同的离散格式对激波进行捕获


采用密度基求解器对激波进行捕获

点击下载

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