Another steady-state SIMPLE solver for compressible flow

I developed another steady-state SIMPLE solver for compressible flow which is in the well-known form of SIMPLE method.

In simpleFoam, the equation of the corrector step is designated to solve the pressure.

But the well-known SIMPLE method solves the pressure equation for p'.

I created a new rhoSimpleFoam which solves the equation for p' by applying the following changes to the original rhoSimpleFoam:

1. Adds new volScalarField pPrime (= p')
2. Uses U* instead of HbyA in the pressure equation.

1.
In createFields.H, new volScalarField pPrime is read.
The boundary conditions are identilal to those for p, except that the fixed value usually on the outlet boundary is zero.

Info<< "Reading field pPrime\n" << endl;
volScalarField pPrime
(
    IOobject
    (
        "pPrime",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

2.

In pEqn.H, U* is used instead of HbyA.

volScalarField rAU(1.0/UEqn().A());

//volVectorField HbyA("HbyA", U);

//HbyA = rAU*UEqn().H();

// At this stage, U = U*

surfaceScalarField phiUstar

(

    "phiUstar",

    fvc::interpolate(rho*U) & mesh.Sf()

);

fvOptions.makeRelative(fvc::interpolate(rho), phiUstar);

closedVolume = adjustPhi(phiUstar, U, p);

while (simple.correctNonOrthogonal())

{

    fvScalarMatrix pPrimeEqn

    (

        fvc::div(phiUstar)

      - fvm::laplacian(rho*rAU, pPrime)

      ==

        fvOptions(psi, pPrime, rho.name())

    );

    pPrimeEqn.setReference(pRefCell, pRefValue);

    fvOptions.constrain(pPrimeEqn);

    pPrimeEqn.solve();

    if (simple.finalNonOrthogonalIter())

    {

        phi = phiUstar + pPrimeEqn.flux();

    }

}

And correct p and U.

p += pPrime;

U -= rAU*fvc::grad(pPrime);

p.correctBoundaryConditions();

U.correctBoundaryConditions();

fvOptions.correct(U);

Finally, p' is reset to zero and go to the next timestep.

pPrime = min

         (

             max

             (

                 pPrime, dimensionedScalar("0", dimPressure, scalar(0))

             ),

             dimensionedScalar("0", dimPressure, scalar(0))

         );