RhoSimpleFoam

rhoSimpleFoam

rhoSimpleFoam is a steady-state solver for compressible, turbulent flow, using the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm. In the newer releases it also includes an option to use the SIMPLEC (Semi-Implicit Method for Pressure Linked Equations Consistent) algorithm. It is a pressure based compressible solver. That means that a pressure equation is solved and the density is related to the pressure via an equation of state.

1 Solution Strategy

The solver follows a segregated solution strategy. This means that the equations for each variable characterizing the system (the velocity $\bold u$, the pressure $p$, the energy (either internal energy or enthalpy depending on the choice of the user) and the variables characterizing turbulence) is solved sequentially and the solution of the preceding equations is inserted in the subsequent equation. The non-linearity appearing in the momentum equation (the face flux $\phi_f$ which is a function of the velocity and density $\rho_f$ ) is resolved by computing it from the velocity and pressure values of the preceding iteration. Note that the density is related to the pressure via an equation of state which can be selected by the user. The dependence from the pressure is introduced to avoid a decoupling between the momentum and pressure equations and hence the appearance of high frequency oscillation in the solution (check board effect). The first equation to be solve is the momentum equation. It delivers a velocity field $\bold u^*$ which is in general not divergence free, i.e. it does not satisfy the continuity equation. After that the energy equation is solved. Then that the momentum and the continuity equations are used to construct an equation for the pressure. The aim is to obtain a pressure field $p^{n}$, which, if inserted in the momentum equation, delivers a divergence free velocity field $\bold u$. After correcting the velocity field, the equations for turbulence are solved. The above iterative solution procedure is repeated until convergence.

A good book which describes the derivation of the pressure equation for compressible cases is the one of .

To summarize the solution algorithm is as follows

• Compute an intermediate velocity $\bold{u^*}$ by means of the momentum equation. This velocity (together with the density $\rho^*$ does not satisfy the continuity equation

• Compute energy equation. The compressibilty $\Psi$ is updated too by the thermodynamic packaged (of course also the temperature the viscosity).

• Compute the pressure $p^{n+1}$ which satisfies the momentum and also the continuity equation.

• Correct the velocity and face flux which satisfy the continuity equation

• Update the density $\rho^{n+1} = \Psi p^{n+1}$

• Update the turbulence quantities

• If not converged go to the beginning of the iteration

The source code can be found in rhoSimpleFoam.C

 
 
int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Steady-state solver for compressible turbulent flow."
    );
 
    #include "postProcess.H"
 
    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createControl.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
    #include "initContinuityErrs.H"
 
    turbulence->validate();
 
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
    Info<< "\nStarting time loop\n" << endl;
 
    while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;
 
        // Pressure-velocity SIMPLE corrector
        #include "UEqn.H"
        #include "EEqn.H"
 
        if (simple.consistent())
        {
            #include "pcEqn.H"
        }
        else
        {
            #include "pEqn.H"
        }
 
        turbulence->correct();
 
        runTime.write();
 
        runTime.printExecutionTime(Info);
    }
 
    Info<< "End\n" << endl;
 
    return 0;
}
 

2 Equations

In the next section the equation to be solved are described. In contrast to incompressible fluids, where the density is constant, in compressible flows the density varies in space and time, i.e., it the density became a thermodynamic property. For a pure mixture the thermodynamic state is fully described by to variables, i.e., $\rho = \rho(p,T)$ the density is a function of the pressure and the temperature. For this reason in compressible flows we need to solve besides the equations for the momentum and the pressure (which is derived from the continuity equations) also an equation which gives us the evolution of the temperature. We will see later that the temperature is obtained via a thermodynamic relation after we solved the energy equation. After that, the density $\rho$ can be obtained from the other two variables via an equation of state.

2.1 Momentum Equation

The equation of motion in rhosimpleFoam are written for a moving frame of reference. They are however formulated for the absolute velocity (the derivation of the equations of motion can be found in https://openfoamwiki.net/index.php/See_the_MRF_development and also in https://diglib.tugraz.at/download.php?id=581303c7c91f9&location=browse. Some additional information can be found in https://pingpong.chalmers.se/public/pp/public_courses/course07056/published/1497955220499/resourceId/3711490/content/UploadedResources/HakanNilssonRotatingMachineryTrainingOFW11-1.pdf):

- $\nabla \cdot ( \rho \bold{u_R} \otimes \bold{u_I}) + \rho \bold{\Omega} \times \bold{u_I} = \nabla p + \nabla \cdot (\bold{R_{vis}} + \bold{R_{tur}} )$

(1)

$\bold{u_R}$ represents the velocity in the relative reference frame, $\bold{u_I}$ is the velocity described in the inertial reference frame, $\bold{\Omega}$ is the angular velocity of the moving reference frame, $\bold{R_{vis}} + \bold{R_{tur}}$ is the sum of the viscous and turbulent stress tensor. It is evident that the above equation can be interpreted as the Favre averaged Navier Stokes equation rather than the Reynolds averaged Navier stokes equation. Favre averaging is a transformation which avoids the appearance of correlation involving the density fluctuations in the continuity and momentum equation. For a more detailed discussion see the thread https://www.cfd-online.com/Forums/openfoam-solving/193246-favre-reynolds-average-openfoam.html or the book of .

The source code can be found in the file UEqn.H:

 
    // Solve the Momentum equation
 
    MRF.correctBoundaryVelocity(U);
 
    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(U)
     ==
        fvOptions(rho, U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();
 
    UEqn.relax();
 
    fvOptions.constrain(UEqn);
 
    solve(UEqn == -fvc::grad(p));
 
    fvOptions.correct(U);

Note that the above equation differs from the equations without considering a moving frame of reference in two terms: $\rho \bold{u_R}$ which if interpolated on the face of the control volume gives the relative face flux and the term $\rho \bold{\Omega} \times \bold{u_I}$ which is simply a explicit source term in the momentum equation. We will see later in the pressure equation how the relative face flux is calculated.

The source code of the acceleration resulting from the description in a moving frame of reference can be found in the following src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C

 
Foam::tmp<Foam::volVectorField> Foam::MRFZoneList::DDt
(
    const volVectorField& U
) const
{
    tmp<volVectorField> tacceleration
    (
        new volVectorField
        (
            IOobject
            (
                "MRFZoneList:acceleration",
                U.mesh().time().timeName(),
                U.mesh()
            ),
            U.mesh(),
            dimensionedVector(U.dimensions()/dimTime, Zero)
        )
    );
    volVectorField& acceleration = tacceleration.ref();
 
    forAll(*this, i)
    {
        operator[](i).addCoriolis(U, acceleration);
    }
 
    return tacceleration;
}
 
 
Foam::tmp<Foam::volVectorField> Foam::MRFZoneList::DDt
(
    const volScalarField& rho,
    const volVectorField& U
) const
{
    return rho*DDt(U);
}

The calculation of the Coriolis force is done in the file src/finiteVolume/cfdTools/general/MRF/MRFZone.C

 
 
void Foam::MRFZone::addCoriolis
(
    const volVectorField& U,
    volVectorField& ddtU
) const
{
    if (cellZoneID_ == -1)
    {
        return;
    }
 
    const labelList& cells = mesh_.cellZones()[cellZoneID_];
    vectorField& ddtUc = ddtU.primitiveFieldRef();
    const vectorField& Uc = U;
 
    const vector Omega = this->Omega();
 
    forAll(cells, i)
    {
        label celli = cells[i];
        ddtUc[celli] += (Omega ^ Uc[celli]);
    }
}

The way the turbulent stresses enter into the momentum equation is well described here: https://www.openfoam.com/documentation/guides/latest/doc/guide-turbulence-ras.html.

2.2 Energy Equation

The derivation of the Favre agraged energy equation for compressible flows can be found in the book of . A description of the terms resulting in the energy equation from the favre averaging can also be found in https://www.cfd-online.com/Wiki/Favre_averaged_Navier-Stokes_equations. The following equation are taken from .

$\frac{\partial}{\partial t} (\overline{\varrho} \tilde{e_0})+\nabla\cdot\left(\overline{\varrho}\tilde{e_0}\widetilde{\mathbf{u}}+\tilde{\mathbf{u}}\overline{p} + \mathbf{q_{vis}} + \mathbf{q_{tur}} - \mathbf{\tilde{R}_{tot} \cdot \tilde{u}} - \overline{\mathbf{R_{vis}} \cdot \mathbf{u''}} -0.5\overline{\rho \mathbf{u''} (\mathbf{u''} \cdot \mathbf{u''})} \right) = 0$

The total energy $\tilde{e_0}$ is defined as follows:

$\tilde{e_0} = \tilde{e} + 0.5 \tilde{\mathbf{u}}\cdot \tilde{\mathbf{u}} + k$. The overline $\overline{ }$ denotes the Reynolds averaging, the tilde the Favre averaging, k the turbulent kinetic energy and the two $''$ denote the fluctuation due to the Favre averaging:

$\phi = \tilde{\phi} + \phi''$.

$\tau$ is the viscous stress tensor.

If we neglect the turbulent kinetic energy in the total energy and the last three terms of the above equation (together with the time derivative since rhoSimpleFoam is a steady solver) we get:

$\nabla\cdot\left(\overline{\varrho}(\tilde{e} + 0.5 \tilde{\mathbf{u}}\cdot \tilde{\mathbf{u}} )\widetilde{\mathbf{u}}+\tilde{\mathbf{u}}\overline{p} + \mathbf{q_{tot}} \right) = 0$

If we insert the definition of the enthalpy h:

$\tilde{h} = \tilde{e} + \frac{ \overline{p} }{\overline{ \rho}}$

in the above equation we get the following equation:

$\nabla\cdot\left(\overline{\varrho}(\tilde{h} + 0.5 \tilde{\mathbf{u}}\cdot \tilde{\mathbf{u}} )\widetilde{\mathbf{u}} + \mathbf{q_{tot}} \right) = 0$

The total heat flux $\mathbf{q_{tot}}$ is the sum of viscous and turbulent heat flux $\mathbf{q_{tot}} = \mathbf{q_{vis}} + \mathbf{q_{tur}}$. It is calculated as follows in case the sensible enthalpy is chosen as quantity solved for:

$\mathbf{q_{tot}} = - \frac{\mu_t}{Pr_t}\nabla h$

and as follows if the sensible internal energy is chosen as quantity solved for:

$\mathbf{q_{tot}} = - \frac{\mu_t}{Pr_t}\frac{c_p}{c_v}\nabla e$

The reason is to have the same equation whether one solves for e or h. The way this is achieved can be found in the following discussion: https://www.cfd-online.com/Forums/openfoam/232183-question-about-fvm-laplacian-turbulence-alphaeff-he-rhosimplefoam.html.

Differences between the two equation arise in the numerical stability. Since for the two equation two different explicit source terms appear, one equation may be more stable compared to the other depending on the magnitude and sign of the explicit source term. Also here a discussion can be found here: https://www.cfd-online.com/Forums/openfoam-solving/228152-when-choose-sensibleinternalenergy-sensibleenthalpy-thermophysicalproperties.html.

The equation solved can be found in the file EEqn.H

 
{
    volScalarField& he = thermo.he();
 
    fvScalarMatrix EEqn
    (
        fvm::div(phi, he)
      + (
            he.name() == "e"
          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
        )
      - fvm::laplacian(turbulence->alphaEff(), he)
     ==
        fvOptions(rho, he)
    );
 
    EEqn.relax();
 
    fvOptions.constrain(EEqn);
 
    EEqn.solve();
 
    fvOptions.correct(he);
 
    thermo.correct();
}

2.2.1 Correction of the thermodynamic properties

After that the energy equation is solved, the thermodynamic properties of the fluid are updated. This is done by the following code line:

 
    thermo.correct();

To elucidate what is done exactly we'll take the hePsiThermo package as example. In this package the density is related to the pressure via the compressibilty $\Psi$.

The function correct looks like this (see the file $FOAM_SCR/thermophysicalModels/basic/psiThermo/hePsiThermo.C):

 
template<class BasicPsiThermo, class MixtureType>
void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::correct()
{
    DebugInFunction << endl;
 
    calculate
    (
        this->p_,
        this->T_,
        this->he_,
        this->psi_,
        this->mu_,
        this->alpha_,
        false           // No need to update old times
    );
 
    DebugInFunction << "Finished" << endl;
}
 
template<class BasicPsiThermo, class MixtureType>
void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate
(
    const volScalarField& p,
    volScalarField& T,
    volScalarField& he,
    volScalarField& psi,
    volScalarField& mu,
    volScalarField& alpha,
    const bool doOldTimes
)
{
    // Note: update oldTimes before current time so that if T.oldTime() is
    // created from T, it starts from the unconverted T
    if (doOldTimes && (p.nOldTimes() || T.nOldTimes()))
    {
        calculate
        (
            p.oldTime(),
            T.oldTime(),
            he.oldTime(),
            psi.oldTime(),
            mu.oldTime(),
            alpha.oldTime(),
            true
        );
    }
 
    const scalarField& hCells = he.primitiveField();
    const scalarField& pCells = p.primitiveField();
 
    scalarField& TCells = T.primitiveFieldRef();
    scalarField& psiCells = psi.primitiveFieldRef();
    scalarField& muCells = mu.primitiveFieldRef();
    scalarField& alphaCells = alpha.primitiveFieldRef();
 
    forAll(TCells, celli)
    {
        const typename MixtureType::thermoType& mixture_ =
            this->cellMixture(celli);
 
        if (this->updateT())
        {
            TCells[celli] = mixture_.THE
            (
                hCells[celli],
                pCells[celli],
                TCells[celli]
            );
        }
 
        psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
 
        muCells[celli] = mixture_.mu(pCells[celli], TCells[celli]);
        alphaCells[celli] = mixture_.alphah(pCells[celli], TCells[celli]);
    }
 
    const volScalarField::Boundary& pBf = p.boundaryField();
    volScalarField::Boundary& TBf = T.boundaryFieldRef();
    volScalarField::Boundary& psiBf = psi.boundaryFieldRef();
    volScalarField::Boundary& heBf = he.boundaryFieldRef();
    volScalarField::Boundary& muBf = mu.boundaryFieldRef();
    volScalarField::Boundary& alphaBf = alpha.boundaryFieldRef();
 
    forAll(pBf, patchi)
    {
        const fvPatchScalarField& pp = pBf[patchi];
        fvPatchScalarField& pT = TBf[patchi];
        fvPatchScalarField& ppsi = psiBf[patchi];
        fvPatchScalarField& phe = heBf[patchi];
        fvPatchScalarField& pmu = muBf[patchi];
        fvPatchScalarField& palpha = alphaBf[patchi];
 
        if (pT.fixesValue())
        {
            forAll(pT, facei)
            {
                const typename MixtureType::thermoType& mixture_ =
                    this->patchFaceMixture(patchi, facei);
 
                phe[facei] = mixture_.HE(pp[facei], pT[facei]);
 
                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
                pmu[facei] = mixture_.mu(pp[facei], pT[facei]);
                palpha[facei] = mixture_.alphah(pp[facei], pT[facei]);
            }
        }
        else
        {
            forAll(pT, facei)
            {
                const typename MixtureType::thermoType& mixture_ =
                    this->patchFaceMixture(patchi, facei);
 
                if (this->updateT())
                {
                    pT[facei] = mixture_.THE(phe[facei], pp[facei], pT[facei]);
                }
 
                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
                pmu[facei] = mixture_.mu(pp[facei], pT[facei]);
                palpha[facei] = mixture_.alphah(pp[facei], pT[facei]);
            }
        }
    }
}

We see from the above code that first the temperature is retrieved from the enthalpy or internal energy. For a good explanation how it's done exactly see https://caefn.com/openfoam/temperature-calculation. After that the temperature dependent variables like the dynamic viscosity $\mu$, the compresibility $\Psi$ and the thermal diffusivity $\alpha$ are calculated.

2.3 Pressure Equation

Note that a easy to understand derivation of the pressure or pressure correction equations for compressible flows can be found in .

The derivation of the pressure pressure equation is very similar to the one used for the incompressbile case (see e.g. https://openfoamwiki.net/index.php/SimpleFoam). Also here after the solution of the momentum equation we get a velocity $\mathbf{u_P^*}$ and a density $\rho_P^*$ (that's the main difference to the incompressible case) at the point P which together do not satisfy the continuity equation. In semi discrete form the mass flux at the point P $m_P^*$ can be written as:

$m_P^* = \rho_P^*\bold{u_P^*} = \rho_P^* \frac{\bold {H[u^*] }}{a_P^* } - \rho_P^*\frac{\nabla p_P^{n-1}}{a_P^* }$

Now we want to add a velocity correction $\mathbf{u_p^'}$, a density correction $\rho_p^'$ and a pressure correction $p_p^'$ in order that the continuity equation is satisfied:

$m_P = (\rho_P^* + \rho_P^')(\bold{u_P^*} + \mathbf{u_p^'}) = (\rho_P^* + \rho_P^')(\frac{\bold {H[u^*] }}{a_P^* } + \frac{\bold {H[u^'] }}{a_P^* } )- (\rho_P^* + \rho_P^') (\frac{\nabla p_P^{n-1}}{a_P^* } - \frac{\nabla p_P^'}{a_P^* })$

After applying the following definition to the above equation $p_P = p_P^{n-1} + p_p^'$, $\rho_P = \Psi p_P = \rho_P^{*} + \rho_p^' = \Psi p_P^{n-1} + \Psi p_p^'$:

$m_P = (\rho_P^* + \underbrace{\Psi p_P - \Psi p_P^{n-1}}_{\rho_P^'} )(\frac{\bold {H[u^*] }}{a_P^* } + \frac{\bold {H[u^'] }}{a_P^* } ) - (\rho_P^* + \rho_P^') \frac{\nabla p_P}{a_P^* }$

If we interpolate the above relation to the cell faces, sum the above equation over all faces of the control volume and neglect the density correction in front of the pressure term we get a equation for the pressure:

$\sum_f m_P|_f = 0 = \sum_f (\rho_P^* + \underbrace{\Psi p_P - \Psi p_P^{n-1}}_{\rho_P^'} )(\frac{\bold {H[u^*] }}{a_P^* } + \rho_P^*\frac{\bold {H[u^'] }}{a_P^* } ) |_f - \sum_f\rho_P^*\frac{\nabla p_P}{a_P^* } |_f$

Note that the contribution of the neighboring cells to the velocity correction $\bold {H[u'] }$ is still unknown and has to be approximated. Different approximation are used in the simple algorithm family. In the following the two available in rhoSimpleFoam are described.

2.3.1 Simple

If we neglect the contribution of the neighbours of the velocity correction (i.e. we set $\bold {H[u'] } = 0$) we get the pressure equation solved in the simple algorithm.

In this way we get following equation:

$\sum_f (\rho_P^* + \underbrace{\Psi p_P - \Psi p_P^{n-1}}_{\rho_P^'} )(\frac{\bold {H[u^*] }}{a_P^* } ) |_f - \sum_f\rho_P^*\frac{\nabla p_P}{a_P^* } |_f = 0$

.

From the above equation we notice that the pressure equation has also a convective term (the finite volume discretisation of the divergence is a sum over the cell volumes). For high Mach number flows this convective term transforms the equation for the pressure from a elliptical equation in a hyperbolic equation .

In case the transonic switch in the fvSolution dictionary is set to false, the contribution of the density correction $\rho^'$ is neglected. For this reason the transonic switch has to be set to false only for low Mach number flow where density fluctuations are small.

The source code can be found in the file pEqn.H:

 
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(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);
 
    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();
 
U = HbyA - rAU*fvc::grad(p);
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();
}
 

2.3.2 Simplec

The starting point for the derivation of the simplec algorithm are the following two equations for the uncorrected velocity $\mathbf{u_P^*}$ and for the velocity correction $\mathbf{u_P^'}$:

$\bold{u_P^*} = \frac{\bold {H[u^*] }}{a_P^* } - \frac{\nabla p_P^{n-1}}{a_P^* }$

$\bold{u_P^'} = \frac{\bold {H[u^'] }}{a_P^* } - \frac{\nabla p_P^{'}}{a_P^* } = \frac{\sum_N \mathbf{a_N u_N^'}}{a_P^* } - \frac{\nabla p_P^{'}}{a_P^* }$

The scope is to add both equation and multiply them by $\rho_p^* + \rho^'$ to get the corrected mass flux $m_P$ as function of the pressure (ore pressure correction) solely in order to obtain a solvable equation. The problem is that the second of the above equations involves also the neighboring contribution to the velocity at the point P. What we want is a explicit relation between $\mathbf{u_P^'}$ and $p^'$. To obtain such a relation we can approximate the velocity correction at the neighboring cells $\mathbf{u_N^'}$ by the velocity at the cell center P $\mathbf{u_P^'}$. For small cell sizes the two velocity should not differ too much and it is a good approximation. In this way we get:

$\sum_N \mathbf{a_N u_N^'} \approx \sum_N \mathbf{a_N u_P^'} = \bold {H_1} \mathbf{ u_P^'}$

$\bold{u_P^'} \approx \frac{\sum_N \mathbf{a_N u_P^'}}{a_P^* } - \frac{\nabla p_P^{'}}{a_P^* } = \frac{\bold {H_1 }}{a_P^* }\mathbf{ u_P^'} - \frac{\nabla p_P^{'}}{a_P^* }$

By doing so we get an explicit relation between the velocity correction $\mathbf{u_P^'}$ and the pressure correction $p_P^'$:

$\bold{u_P^'} = - \frac{\nabla p_P^{'}}{a_P^* - \bold{H_1} }$

Now we can use the above relation, sum it to the equation for the uncorrected velocity $\mathbf{u_P^*}$ and multiply it by the corrected density $\rho_P^* + \rho^'$:

$m_P = (\rho^* + \rho^') ( \mathbf{u_P^*} + \bold{u_P^'}) = (\rho^* + \rho^') ( \frac{\bold {H[u^*] }}{a_P^* } - \frac{\nabla p_P^{n-1}}{a_P^* } - \frac{\nabla p_P^{'}}{a_P^* - \bold{H_1}} )$

If we want to have a relation for the pressure $p_P$ instead of the pressure correction $p^'$ we can use following definition:

$p_P = p_P^{n-1} + p^'$

Rearranged for $p^'$ and inserted above we obtain:

$m_P = (\rho^* + \underbrace{\rho^'}_{\Psi p_P - \Psi p_P^{n-1}}) ( \mathbf{u_P^*} + \bold{u_P^'}) = (\rho^* + \underbrace{\rho^'}_{\Psi p_P - \Psi p_P^{n-1}}) ( \frac{\bold {H[u^*] }}{a_P^* } - \frac{\nabla p_P^{n-1}}{a_P^* } - \frac{1 }{a_P^* - \bold{H_1}} (\nabla p_P- \nabla p_P^{n-1} ) )$

If we neglect the contribution of the terms which involve the product of the density correction times the pressure terms we obtain:

$m_P = (\rho^* + \Psi p_P - \Psi p_P^{n-1}) ( \frac{\bold {H[u^*] }}{a_P^* }) - (\frac{\rho_P^*}{a_P^* } - \frac{\rho_P^*}{a_P^* - \bold{H_1}} ) \nabla p_P^{n-1} - \frac{\rho_P^*}{a_P^* - \bold{H_1}} \nabla p_P$

The final equation is obtained by summing over the faces of the control volume and setting it to zero:

$0 = \sum_f (\rho^* - \Psi p_P^{n-1}) ( \frac{\bold {H[u^*] }}{a_P^* })|_f + \sum_f ( \Psi p_P ) ( \frac{\bold {H[u^*] }}{a_P^* })|_f - \sum_f (\frac{\rho_P^*}{a_P^* } - \frac{\rho_P^*}{a_P^* - \bold{H_1}} ) \nabla p_P^{n-1}|_f - \sum_f\frac{\rho_P^*}{a_P^* - \bold{H_1}} \nabla p_P|_f$

The source code can be found in the file pcEqn.H:

 
rho = thermo.rho();
 
volScalarField rAU(1.0/UEqn.A());
volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
volVectorField HbyA(constrainHbyA(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);
 
volScalarField rhorAtU("rhorAtU", rho*rAtU);
 
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
 
if (simple.transonic())
{
    surfaceScalarField phid
    (
        "phid",
        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
    );
 
    phiHbyA +=
        fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
      - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
    HbyA -= (rAU - rAtU)*fvc::grad(p);
 
    while (simple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvc::div(phiHbyA)
          + fvm::div(phid, p)
          - fvm::laplacian(rhorAtU, p)
         ==
            fvOptions(psi, p, rho.name())
        );
 
        // Relax the pressure equation to maintain 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);
 
    phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
    HbyA -= (rAU - rAtU)*fvc::grad(p);
 
    while (simple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvc::div(phiHbyA)
          - fvm::laplacian(rhorAtU, p)
          ==
            fvOptions(psi, p, rho.name())
        );
 
        pEqn.setReference
        (
            pressureControl.refCell(),
            pressureControl.refValue()
        );
 
        pEqn.solve();
 
        if (simple.finalNonOrthogonalIter())
        {
            phi = phiHbyA + pEqn.flux();
        }
    }
}
 
// The incompressible form of the continuity error check is appropriate for
// steady-state compressible also.
#include "incompressible/continuityErrs.H"
 
// Explicitly relax pressure for momentum corrector
p.relax();
 
U = HbyA - rAtU*fvc::grad(p);
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();
}