两类pEqn.H求解可压缩流动计算时稳定性问题

[yfclark]

有对pimple算法求解pEqn.H有研究的老哥吗，关于这个算法有一些困惑，(1)https://github.com/OpenFOAM/OpenFOAM-2.3.x/blob/master/applications/solvers/combustion/reactingFoam/reactingFoam.C
中密度方程为什么要求解两次，一次出现在runtime++后面，一次包含在pEqn.H里面。
（2）reactingFoam和rhoreactingFoam中的pEqn.H求解的压力方程并不相同，而且reactingFoam的pEqn.H和rhoPimpleFoam里面的pEqn.H基本上是一致的。有使用过rhoreactingFoam的大佬吗？rhoreactingFoam在求解可压缩流动是是否比reactingFoam有优势？我现在使用的是reactingFoam这一类的pEqn.H,计算过程中一个时间步内压力变化很大，导致温度出现了震荡现象，导致物性也发生了震荡，无法收敛。
下面贴出这两中pEqn.H
rhoPimpleFoam/reactingFoam

if (pimple.transonic())
{
    surfaceScalarField phid
    (
        "phid",
        fvc::interpolate(psi)
       *(
            (fvc::interpolate(HbyA) & mesh.Sf())
          + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
        )
    );

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

    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvm::ddt(psi, p)
          + fvm::div(phid, p)
          - fvm::laplacian(rhorAUf, p)
          ==
            fvOptions(psi, p, rho.name())
        );

        fvOptions.constrain(pEqn);

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

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

rhoreactingFoam

 thermo.rho() -= psi*p;

    volScalarField rAU(1.0/UEqn.A());
    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();

    if (pimple.transonic())
    {
        surfaceScalarField phiHbyA
        (
            "phiHbyA",
            (
                (fvc::interpolate(rho*HbyA) & mesh.Sf())
              + rhorAUf*fvc::ddtCorr(rho, U, phi)
            )/fvc::interpolate(rho)
        );

        fvOptions.makeRelative(phiHbyA);

        surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);

        phiHbyA *= fvc::interpolate(rho);

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

        while (pimple.correctNonOrthogonal())
        {
            fvScalarMatrix pEqn
            (
                pDDtEqn
              - fvm::laplacian(rhorAUf, p)
             ==
                fvOptions(psi, p, rho.name())
            );

            fvOptions.constrain(pEqn);

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

            if (pimple.finalNonOrthogonalIter())
            {
                phi = phiHbyA + pEqn.flux();
            }
        }
    }
    thermo.rho() += psi*p;

只贴了transonic部分的

[cccrrryyy]

我最近也在思考这个问题，不知道你这边有没有什么进展。我看了OpenFOAM-4和OpenFOAM-6还有一些其他版本的，分享一下我目前的一些看法。

就你发的这个2.3.x而言，主程序里面runtime++后面的#include "rhoEqn.H"是一个根据速度分布求解密度的方程，应该只是为了有一个初始的密度场，不然后续的计算无法进行。在真正开始pimple循环后，EEqn.H末尾有一个thermo.correct()。这一步根据计算后的能量场，更新了温度、粘度、可压缩比（psi）等等物性，所以在pEqn.H的开头，出现了一步rho = thermo.rho()，实际上进行了更新密度的操作（rho = psi * p，注意这个时候psi已经变了），之后再进行压力的求解。注意，在pEqn.H的79行和80行，再次调用"rhoEqn.H"（计算根据速度场应该得到的密度场，或者说是没有更新之前的密度场），以及 "compressibleContinuityErrs.H"（计算前面所说的密度场和由热物性所得到的密度场之间的区别，即密度场更新前和更新后的差异）。不知道我以上的理解对不对。

我还没弄明白的问题和你一样，reactingFoam和rhoReactingFoam的压力方程为什么代码上有差异。看看有没有别人回答吧。