rhoPimpleFoam解析


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

1. 引言

rhoPimpleFoam为一个主要用于暖通类计算的可压缩求解器,其也可以用于亚音速流动。另外一个类似的求解器为buoyantPimpleFoam,二者的主要区别在于后者考虑了浮力的作用。buoyantPimpleFoam不能处理接近音速流动,其主要用于暖通类应用,如温度引起的浮力驱动流。rhoPimpleFoam则倾向为一个无重力普适性可压缩求解器。本文对rhoPimpleFoam的解析在icoFoam解析的基础上进行适配,建议首先阅读icoFoam解析获取更基础的算法信息。

2. 控制方程与算法

首先有可压缩流体的N-S方程:

\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}}$表示有效导热,$\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}即为求解的transonic压力方程。其中的$\psi \mathbf{HbyA}_\mathrm{P}$即为代码中的phid。 对比所有的可压缩求解器和不可压缩求解器,会发现不可压缩求解器的压力方程为一个纯粹的泊松方程。在可压缩求解器中,压力方程包含了对流项和时间项,其实际上是对流波动方程。可压缩求解器中压力方程中的对流项和扩散项的相对重要性和速度有关,在速度较低的情况下,拉普拉斯项占主要作用,可压缩压力方程近似于不可压缩的压力泊松方程。在速度较高的情况下,对流项占主导地位,其反映了流场本身的双曲特性。这样,可压缩求解器中的压力方程等于求解连续性方程(密度)。

3. 代码分析
下面进入代码分析,速度方程以及能量方程请参考buoyantPimpleFoam解析,以及CFD中的能量方程。我们主要看压力方程pEqn.H
if (pimple.transonic())
{
    surfaceScalarField phid //transonic中的phid
    (
        "phid",
        fvc::interpolate(psi)
       *(
            (fvc::interpolate(HbyA) & mesh.Sf())
          + fvc::ddtPhiCorr(rAU, rho, U, phi)
        )
    );

    fvOptions.relativeFlux(fvc::interpolate(psi), phid);

    volScalarField Dp("Dp", rho*rAU);

    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn //方程($\ref{NS5}$)
        (
            fvm::ddt(psi, p)
          + fvm::div(phid, p)
          - fvm::laplacian(Dp, p)
          ==
            fvOptions(psi, p, rho.name())
        );

        fvOptions.constrain(pEqn);

        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
            phi == pEqn.flux();
        }
    }
}
else
{
    surfaceScalarField phiHbyA 
    (
        "phiHbyA",
        fvc::interpolate(rho)
       *(
            (fvc::interpolate(HbyA) & mesh.Sf())
          + fvc::ddtPhiCorr(rAU, rho, U, phi)
        )
    );

    fvOptions.relativeFlux(fvc::interpolate(rho), phiHbyA);

    volScalarField Dp("Dp", rho*rAU);

    while (pimple.correctNonOrthogonal())
    {
        // Pressure corrector
        fvScalarMatrix pEqn //方程($\ref{p2}$)
        (
            fvm::ddt(psi, p)
          + fvc::div(phiHbyA)
          - fvm::laplacian(Dp, p)
          ==
            fvOptions(psi, p, rho.name())
        );

        fvOptions.constrain(pEqn);

        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
            phi = phiHbyA + pEqn.flux();
        }
    }
}

东岳流体dyfluid.com
勘误、讨论、补充内容请前往CFD中国