• 返回主页
  • rhoPimpleFoam解析

    2017.06.23:修改前言部分


    求解器类型:压力基可压缩瞬态求解器。

    rhoPimpleFoam更像pimpleFoam的可压缩版本,因此可以处理亚音速流动。另外一个类似的求解器为buoyantPimpleFoam,二者的主要区别在于是否考虑浮力,rhoPimpleFoam并没有包含浮力的作用。并且buoyantPimpleFoam不能处理接近音速流动,其主要用于暖通类应用。在可压缩压力基求解器中,压力从压力泊松方程或者压力修正方程中计算而来(依据连续性方程和动量方程联立)。

    首先有可压缩流体的NS方程: \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 \left(\mu \nabla \mathbf{U} \right)=-\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} 具体的有关方程($\ref{NS1}$)以及方程($\ref{NS2}$)中的操作符含义以及其是如何离散的请参考icoFoam解析。下面对其进行简述。

    方程($\ref{NS1}$)在进行半离散化之后有(此处并未包含压力梯度): \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} 将方程($\ref{Upreal}$)代入到方程($\ref{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(\rho \frac{1}{{{A_\mathrm{P}}}}\nabla p \right)=0 \label{p2pre} \end{equation} 对于速度较低的情况,把方程($\ref{bar}$)替换方程($\ref{p2}$)中的时间导数项即为最终的压力方程: \begin{equation} \frac{\partial \psi p}{\partial t}+\nabla \cdot \left(\rho\mathbf{HbyA}_\mathrm{P} \right)-\nabla \cdot \left(\rho \frac{1}{{{A_\mathrm{P}}}}\nabla p \right)=0 \label{p2} \end{equation} 方程($\ref{p2}$)可以继续化简为: \begin{equation} \frac{\partial \psi p}{\partial t}+\nabla \cdot \left( \psi\mathbf{HbyA}_\mathrm{P} \cdot p \right)-\nabla \cdot \left(\rho \frac{1}{{{A_\mathrm{P}}}}\nabla p \right)=0 \label{NS5} \end{equation} 方程($\ref{NS5}$)即为求解的transonic压力方程。其中的$\psi \mathbf{HbyA}_\mathrm{P}$即为代码中的phid

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

    下面进入代码分析,速度方程以及能量方程请参考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中国