pEqn.H

Go to the documentation of this file.

{
    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
 
    volScalarField rAU(1.0/UEqn.A());
    surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU));
    volVectorField HbyA("HbyA", U);
    HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
    tUEqn.clear();
 
    bool closedVolume = false;
 
    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
 
    if (simple.transonic())
    {
        surfaceScalarField phid
        (
            "phid",
            (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
        );
 
        phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
        while (simple.correctNonOrthogonal())
        {
            fvScalarMatrix pEqn
            (
                fvc::div(phiHbyA)
              + fvm::div(phid, p)
              - fvm::laplacian(rhorAUf, p)
            ==
                fvOptions(psi, p, rho.name())
            );
 
            // Relax the pressure equation to ensure diagonal-dominance
            pEqn.relax();
 
            pEqn.setReference
            (
                pressureControl.refCell(),
                pressureControl.refValue()
            );
 
            pEqn.solve();
 
            if (simple.finalNonOrthogonalIter())
            {
                phi = phiHbyA + pEqn.flux();
            }
        }
    }
    else
    {
        closedVolume = adjustPhi(phiHbyA, U, p);
 
        if (adjustFringe)
        {
            oversetAdjustPhi(phiHbyA, U);
        }
 
        while (simple.correctNonOrthogonal())
        {
            fvScalarMatrix pEqn
            (
                fvc::div(phiHbyA)
              - fvm::laplacian(rhorAUf, p)
            ==
                fvOptions(psi, p, rho.name())
            );
 
            pEqn.setReference
            (
                pressureControl.refCell(),
                pressureControl.refValue()
            );
 
            pEqn.solve();
 
            if (simple.finalNonOrthogonalIter())
            {
                phi = phiHbyA + pEqn.flux();
            }
        }
    }
 
    #include "incompressible/continuityErrs.H"
 
    // Explicitly relax pressure for momentum corrector
    p.relax();
 
    volVectorField gradP(fvc::grad(p));
 
    U = HbyA - rAU*cellMask*gradP;
    U.correctBoundaryConditions();
    fvOptions.correct(U);
 
 
    bool pLimited = pressureControl.limit(p);
 
    // For closed-volume cases adjust the pressure and density levels
    // to obey overall mass continuity
    if (closedVolume)
    {
        p += (initialMass - fvc::domainIntegrate(psi*p))
            /fvc::domainIntegrate(psi);
    }
 
    if (pLimited || closedVolume)
    {
        p.correctBoundaryConditions();
    }
 
    rho = thermo.rho();
 
    if (!simple.transonic())
    {
        rho.relax();
    }
}