CFD Online Discussion Forums - Update of face fluxes in simpleFoam

Hello to all,

I am new o CFD in general and just started studying the SIMPLE algorithm (in simpleFoam).

The continuity equation for incompressible flows is:

$\nabla \cdot \vec{v} = 0$

And the momentum equation in a semi-discritezed manner is:

$a_P \vec{v}_P + \sum_N a_N \vec{v}_N = - \nabla P$

or alternatively

$a_P \vec{v}_P = H(\vec{v}) - \nabla P$

We can now write an expression for $\vec{v}_P$ as:

$\vec{v}_P = a_P^{-1} H(\vec{v}) - a_P^{-1}\nabla P$

To derive the pressure correction equation we must replace the $\vec{v}_P$ into the continuity equation. Which results in:

$\nabla \cdot \left [a_P^{-1} H(\vec{v}) \right ] = \nabla \cdot \left [a_P^{-1}\nabla P \right ]$

We must also compute the face fluxes. In this case:

$F = \vec{v}_f \cdot \vec{S_f} = \left (a_P^{-1} H(\vec{v}) \right )_f- \left (a_P^{-1}\nabla P \right )_f$

In simpleFOAM:

$a_P^{-1}$ is defined as:

Code:

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

$a_P^{-1} H(\vec{v})$

Code:

volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));

$\left (a_P^{-1} H(\vec{v}) \right )_f$

Code:

surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));

Pressure Equation as:

Code:

       

 fvScalarMatrix pEqn

        (

            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)

        );

rAtU() is equal to rAU for the SIMPLE algorithm

The velocity is reconstructed as:

Code:

U = HbyA - rAtU()*fvc::grad(p);

My doubt is in the update of the fluxes.

In the code, we have :

Code:

phi = phiHbyA - pEqn.flux();

Is this equivalent to having (?):
$F = \left (a_p^{-1} H(\vec{v}) \right )_f - \left ( a_p^{-1} H(\vec{v}) - a_p^{-1} \nabla p \right )\cdot \vec{S_f}$

In addition, why does the equation change in a compressible case (rhoSimpleFoam) to:

Code:

phi = phiHbyA + pEqn.flux();

Sorry in advance if this is too basic. I am still "breaking stone". I would also appreciate recommendations for literature showing how these coupling algorithms are implemented in OpenFoam