rhoPimpleFoam解析


李东岳


1. 引言

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

可压缩气体在航空领域非常多见。例如近音速飞行器表面的流动以及高马赫数航天飞行器的外流场。一些叶轮机、涡轮旋转机械内,也会因为高压的存在导致密度产生变化。在这些情况下,都需要用可压缩算法来进行求解。由于压力不仅仅和密度相关,可压缩气体通常被认为是斜压的,即压力的等值面与密度的等值面并非是垂直关系。

可压缩液体的典型工况为密闭容器内注入液体、高压容器内泄压放出液体等,在这些情况下需要考虑液体的可压缩性。虽然在这些工况下密度的变化是很小的,但是却会引起压力波的传递。如果液体流动方向在某一时刻发生骤变,也会引起压力波动以声速进行传递,即水锤以及液体气动活塞(Suponitsky et al., 2017)。水锤是大量工程问题产生的主要原因,如噪音、震动损耗等,严重情况下,这种极高的压力传递会将管道击穿。一些工业器件的使用可以减少水锤的产生,如急泄阀、膨胀箱等。在旧版本的OpenFOAM中,由于负责的继承关系,可压缩液体通过另外一个求解器sonicLiquidFoam来处理。在新版的OpenFOAM中,rhoPimpleFoam包含了sonicLiquidFoam的功能,sonicLiquidFoam求解器也进而被略去。

如果打不开图像,请右键在新标签页打开图像后刷新几次 如果打不开图像,请右键在新标签页打开图像后刷新几次

图1. 左图:Suponitsky等模拟的气动活塞(Suponitsky et al., 2017)。右图:水锤的产生。当流体自由流动的时候,不会有压力波的传递。在阀门突然截止的时候,突然产生压力波的传递,即水锤。
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}\cdot\mathrm{d}\bfS\mathrm{d}t} = } \sum {{{\left( {{\rho^t\mathbf{U}^*}{\mathbf{U}^t}} \right)}_f}} \cdot\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^{*} = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{*}} + S_\rP^t} \right) \label{hbya2} \end{equation}
同时有:
\begin{equation} \bfU_\rP^{*}=\bfHbyA_\rP^{*} - \frac{1}{{{A_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{t}\bfS_f \label{Up2} \end{equation}
\begin{equation} \bfU_f^{*}=\bfHbyA_f^{*} - \frac{1}{{{A_{\mathrm{P},f}}}} \left(\frac{1}{V_\rP}\sum_f p_f^{t}\bfS_f\right)_f \label{Up3} \end{equation}
现考虑连续性方程,其最终的离散形式为:
\begin{equation} \frac{\rho^{t+\dt}-\rho^t}{\dt}+\sum_f (\rho_f^t+\rho_f')(\mathbf{U}_f^{*}+\bfU_f') \cdot \bfS_f=0. \label{cont0} \end{equation}
参考rhoSimpleFoam解析,其最终可以写为
\begin{equation} \frac{\rho^{t+\dt}-\rho^t}{\dt}+\sum_f \left(\rho_f^t\bfU_f^{*} + \rho_f^n\bfU_f' + \frac{1}{RT}p_f'\bfU_f^{*} \right)\cdot\bfS_f=0. \label{cont} \end{equation}
再次参考rhoSimpleFoam解析,略去邻点影响有: \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}
\begin{equation} \bfU_\rP^{t+\dt}=\bfHbyA_\rP^{*} - \frac{1}{{{A_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{t+\dt}\bfS_f \label{Up22} \end{equation}
将方程\eqref{uFf}、\eqref{Up3}同时代入至\eqref{cont}有:
\begin{multline} \frac{\rho^{t+\dt}-\rho^t}{\dt} \\ +\sum_f \left(\rho_f^t\left(\bfHbyA_f^{*} - \frac{1}{{{A_{\mathrm{P},f}}}} \left(\frac{1}{V_\rP}\sum_f p_f^{t}\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^{*} \right)\cdot\bfS_f \\ =0. \label{cont2} \end{multline}
引入修正压力: \begin{equation} p_f^{t+\dt}=p_f^t+p_f' \end{equation} 整理方程\eqref{cont2}有
\begin{equation} \frac{\rho^{t+\dt}-\rho^t}{\dt}+\sum_f \rho_f^n \mathbf{HbyA}^{*}_f\cdot\bfS_f + \sum_f \frac{1}{RT}p_f'\bfU_f^{*}\cdot\bfS_f = \sum_f \left(\frac{\rho_f^n}{{{A_f^n}}} \left(\frac{1}{V_\rP}\sum_f p_f^{t+\dt}\bfS_f\right)_f \right)\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$趋向于无穷大,$RT$项可以省略。有:
\begin{equation} \frac{\rho^{t+\dt}-\rho^t}{\dt}+\sum_f \rho_f^n \mathbf{HbyA}^{*}_f\cdot\bfS_f = \sum_f \left(\frac{\rho_f^n}{{{A_f^n}}} \left(\frac{1}{V_\rP}\sum_f p_f^{t+\dt}\bfS_f\right)_f \right)\cdot\bfS_f \label{pres} \end{equation}
其即为亚音速情况下连续形式的压力泊松方程:
\begin{equation} \frac{\p\rho}{\p t}+\nabla \cdot \left(\rho\mathbf{HbyA} \right)=\nabla \cdot \left( \frac{\rho}{{{A}}}\nabla p \right) \label{p2pre} \end{equation}
需要注意的是,方程\eqref{pres}并不是精准的,因为其中假定邻点的影响可以忽略。因此,其求得的$p^{t+\dt}$并不是当前时间步下收敛的最终结果。同时,在本文的算法中,可压缩求解器中的密度一个是通过方程\eqref{NS1}进行更新,一个是通过压力以及状态方程进行求解。俩种方式最终的密度应该相同。因此,在rhoPimpleFoam中,通过对比二者的密度差进行误差测试。具体的,在每一次压力迭代求解过程中,更新速度以及通量,求解密度方程获得通过求解连续性方程的密度解。另一方面,通过状态方程获得密度二者的差,即为连续性误差。总而言之,可压缩求解器中的迭代过程可以表示为下面几个步骤:

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

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

3. 验证算例
3.1 翼型

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

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

3.1. 模型设置

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

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

3.2. 算例运行

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

./Allrun

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

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


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


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

点击下载
3.2 掺合阀

本算例来自于CFD中文网提供的算例。算例工况如下所示,针对一个T型管,存在两个进口一个出口,不同温度的空气注入后进行混合,求解器需要模拟掺混效果。该算例用户设置后持续发散,经调试后可顺利运行。


工况介绍

计算过程首先需要1)运行foam3DMeshToFoam Tee.mesh来转化网格,2)运行decomposePar来分解网格,3)运行mpirun -np 32 rhoPimpleFoam -parallel来计算求解。初步计算结果如下:


计算结果:温度(上)、密度(中)、速度(下)

本算例网格70万,由于网格数量过多,调试时间过长,因此针对结果提出以下可改进的建议:


高估的湍流粘度 | 小圆管可去掉

点击下载

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