sprayFoam计算等压环境下点火延迟时间

香柏树

目前在尝试使用sprayFoam等压环境下的C12H26点火延迟时间，由于sprayFoam中通过泊松方程求解压力（见原始pEqn.H），网格内监测压力显示不能保持等压状态。因此参考chemFoam的源代码，想把泊松方程改为完全气体状态方程P=rhoRT(见修改的pEqn.H)，同时在sprayFoam.C中注释掉压力修正的过程（见修改的sprayFoam.C）。另外能量方程也改为求解Enthalpy的形式。在OpenFOAM6环境下编译没问题，但是计算两个时间步后即浮点溢出了，计算设置见下楼的图1，错误提示下楼的图2。请问各位前辈有没有解决方法啊？
原始pEqn.H

rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();

volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));

if (pimple.nCorrPISO() <= 1)
{
    tUEqn.clear();
}

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

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

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

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

        if (pimple.finalNonOrthogonalIter())
        {
            phi == pEqn.flux();
        }
    }
}
else
{
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        (
            fvc::flux(rho*HbyA)
          + rhorAUf*fvc::ddtCorr(rho, U, phi)
        )
    );

    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);

    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);

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

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

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

#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"

// Explicitly relax pressure for momentum corrector
p.relax();

// Recalculate density from the relaxed pressure
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
Info<< "rho max/min : " << max(rho).value()
    << " " << min(rho).value() << endl;

U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);

if (thermo.dpdt())
{
    dpdt = fvc::ddt(p);
}

修改的pEqn.H

{
    rho = thermo.rho();
    p = rho*Rspecific*thermo.T();  
}

修改的sprayFoam.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    sprayFoam

Description
    Transient solver for compressible, turbulent flow with a spray particle
    cloud.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "turbulentFluidThermoModel.H"
#include "basicSprayCloud.H"
#include "psiReactionThermo.H"
#include "CombustionModel.H"
#include "radiationModel.H"
#include "SLGThermo.H"
#include "pimpleControl.H"
#include "fvOptions.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "postProcess.H"

    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createControl.H"
    #include "createTimeControls.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
    #include "compressibleCourantNo.H"
    #include "setInitialDeltaT.H"
    
    turbulence->validate();

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.run())
    {
        #include "readTimeControls.H"
        #include "compressibleCourantNo.H"
        #include "setDeltaT.H"

        runTime++;

        Info<< "Time = " << runTime.timeName() << nl << endl;

        parcels.evolve();

        #include "rhoEqn.H"
        #include "YEqn.H"
        #include "EEqn.H"
        #include "pEqn.H"

        /*
        if (!pimple.frozenFlow())
        {
            #include "rhoEqn.H"

            // --- Pressure-velocity PIMPLE corrector loop
            while (pimple.loop())
            {
                // #include "UEqn.H"
                #include "YEqn.H"
                #include "EEqn.H"

               
                // --- Pressure corrector loop
                while (pimple.correct())
                {
                    #include "pEqn.H"
                }

                if (pimple.turbCorr())
                {
                    turbulence->correct();
                }
            }


            rho = thermo.rho();

            if (runTime.write())
            {
                combustion->Qdot()().write();
            }
        }
        else
        {
            if (runTime.writeTime())
            {
                parcels.write();
            }
        }
         */

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //

香柏树

68816195-0e44-46d3-9f95-df533596c80d-图片.png
f834dba5-f097-4c55-ac86-66b125a60c38-图片.png

李东岳

@香柏树在 sprayFoam计算等压环境下点火延迟时间中说：

因此参考chemFoam的源代码，想把泊松方程改为完全气体状态方程P=rhoRT(见修改的pEqn.H)

你说的点火之类我不太熟悉。但针对你的这方面描述，chemFoam不涉及到整个场的连续性。sprayFoam要满足质量守恒。你这样做从算法上有点奇怪。

后续的截图看起来发散了，除以0了，当然原因肯定是非常本质的。

香柏树

@李东岳谢谢李老师指导，这么修改确实不太正确。我从算法上看看如何能够合理些。

星星星星晴

@香柏树除以0可能是最大的问题，但是0在哪呢，请问你怎么处理parcel 极小的部分的，会不会是parcel部分，因为你点火了，就出现了这种情况，在clockTime以后就是计算Co 然后solving parcel了，也有可能是你parcel的maxCo=1， parcel在你计算的时候已经飞出了对应的cell 也有可能报错，尝试分析分析这俩地方，看看你的parcel都什么样，看看你的MaxCo减小到0.1或者0.01看看，还是turbulence爆了

星星星星晴

@香柏树另外你直接把整个PIMPLE循环给关掉了，是不是保留一下然后PISO模式运行就行呢, 瞬态

香柏树

@星星星星晴谢谢前辈！我尝试下。

流体中二少年

前辈，你好！请问sprayFoam中的液滴蒸发率怎么获得？