rhoPimpleFoam解析


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

1. 引言

可压缩流动的求解可分为不同的方法。在严格双曲的情况下,瞬态可压缩流动可利用双曲特征的优势采用显性密度基步进求解。在具备混合型特征的可压缩流动中,通常需要进行额外的稳定性适配。一些比较常见的方法有预条件法、压力速度耦合法等。

rhoPimpleFoam为一个普适性无重力可压缩求解器,主要为求解亚音速流动以及暖通类应用设计,但也可用于求解跨/超音速流动。另外类似的求解器为buoyantPimpleFoam,二者的主要区别在于后者考虑了浮力的作用,并且buoyantPimpleFoam只能处理亚音速流动,其主要用于温度引起的浮力驱动流。同类型求解器还有sonicFoam,其主要针对跨音速流动。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$表示可压缩性。首先对方程\eqref{NS2}进行半离散化后有(此处并未包含压力梯度):

\begin{equation} \mathbf{U}_\mathrm{P} = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}} + E_\mathrm{P}} \right)=\mathbf{HbyA}_\mathrm{P} \label{Up} \end{equation}

附加压力梯度之后有:

\begin{equation} \mathbf{U}_\mathrm{P} = \mathbf{HbyA}_\mathrm{P} -\frac{1}{{{A_\mathrm{P}}}}\nabla p \label{Upreal} \end{equation}

将方程\eqref{Upreal}代入到方程\eqref{NS1}有:

\begin{equation} \frac{\partial \rho}{\partial t}+\nabla \cdot \left(\rho\left(\mathbf{HbyA}_\mathrm{P} -\frac{1}{{{A_\mathrm{P}}}}\nabla p \right)\right)=0 \label{NS3} \end{equation}

对方程($\ref{NS3}$)进行分解:

\begin{equation} \frac{\partial \rho}{\partial t}+\nabla \cdot \left(\rho\mathbf{HbyA}_\mathrm{P} \right)-\nabla \cdot \left( \frac{\rho}{{{A_\mathrm{P}}}}\nabla p \right)=0 \label{p2pre} \end{equation}

对于跨/超音速的情况,将方程\eqref{bar}替换到方程\eqref{p2pre}中为最终的压力方程:

\begin{equation} \frac{\partial \psi p}{\partial t}+\nabla \cdot \left( \psi\mathbf{HbyA}_\mathrm{P} \cdot p \right)-\nabla \cdot \left( \frac{\rho}{{{A_\mathrm{P}}}}\nabla p \right)=0 \label{NS5} \end{equation}

方程\eqref{NS5}即为求解的压力方程。其中的$\psi \mathbf{HbyA}_\mathrm{P}$即为代码中的phid。

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

3. 代码分析

首先,组建速度方程,即方程\eqref{NS2}

tmp tUEqn
(
    fvm::ddt(rho, U) + fvm::div(phi, U)
  + MRF.DDt(rho, U)
  + turbulence->divDevRhoReff(U)
 ==
    fvOptions(rho, U)
);

求解机械能:

K = 0.5*magSqr(U);

求解能量方程,即方程\eqref{hfinal}:

fvScalarMatrix EEqn
    (
        fvm::ddt(rho, he) + fvm::div(phi, he)
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      + (
            he.name() == "e"
          ? fvc::div
            (
                fvc::absolute(phi/fvc::interpolate(rho), U),
                p,
                "div(phiv,p)"
            )
          : -dpdt
        )
      - fvm::laplacian(turbulence->alphaEff(), he)
     ==
        fvOptions(rho, he)
    );

若在跨/超音速的情况下,求解方程\eqref{NS5}

fvScalarMatrix pDDtEqn
    (
        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
      + fvc::div(phiHbyA) + fvm::div(phid, p)
     ==
        fvOptions(psi, p, rho.name())
    );

若亚音速的情况下,求解显性离散对流项的压力方程:

fvScalarMatrix pDDtEqn
    (
        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
      + fvc::div(phiHbyA)
     ==
        fvOptions(psi, p, rho.name())
    );

最后对湍流进行修正。

4. 验证算例

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

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

4.1. 模型设置

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

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

4.2. 算例运行

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

./Allrun

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

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


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


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

点击下载

东岳流体®版权所有
勘误、讨论、补充内容请前往CFD中文网