SIMPLE scheme (2/3): Openfoam implementation-SimpleFoam

A walk through SimpleFoam wiki page ¹. This is my study note of simpleFoam wiki page with links that helped me understand.

Before simple loop

createFields.H

The included file createFields.H creates volumeScalarField p and volVectorField U from the 0 directory. Then it includes createPhi.H which creates surfaceScalarField phi with value fvc::flux(U) :

Info<< "Reading/calculating face flux field phi\n" << endl;

surfaceScalarField phi
(
    IOobject
    (
        "phi",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    fvc::flux(U)
);

It calls the fvc::flux() function defined in src/finiteVolume/finiteVolume/fvc/fvcFlux.C as

 Foam::tmp<Foam::surfaceScalarField> Foam::fvc::flux
 (
     const volVectorField& vvf
 )
 {
     return scheme<vector>
     (
         vvf.mesh(),
         "flux(" + vvf.name() + ')'
     )().dotInterpolate(vvf.mesh().Sf(), vvf);
 }

This returns a surfaceScalarField with value = Sf . (face interpolated value of volumeVectorField U).

Reference p and turbulence instantiation

label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, simple.dict(), pRefCell, pRefValue);
mesh.setFluxRequired(p.name());


singlePhaseTransportModel laminarTransport(U, phi);

autoPtr<incompressible::turbulenceModel> turbulence
(
    incompressible::turbulenceModel::New(U, phi, laminarTransport)
);

laminarTransport is an object of the singlePhaseTransportModel class. This has a pointer to viscosity model.

 Foam::singlePhaseTransportModel::singlePhaseTransportModel
 (
     const volVectorField& U,
     const surfaceScalarField& phi
 )
 :
     IOdictionary
     (
         IOobject
         (
             "transportProperties",
             U.time().constant(),
             U.db(),
             IOobject::MUST_READ_IF_MODIFIED,
             IOobject::NO_WRITE
         )
     ),
     viscosityModelPtr_(viscosityModel::New("nu", *this, U, phi))
 {}

turbulence is a pointer to the new turbulence model created by the statement incompressible::turbulenceModel::New(U, phi, laminarTransport)

Checking MRF options:

#include “createMRF.H” has only one line IOMRFZoneList MRF(mesh);

createIOobject (src/finiteVolume/cfdTools/general/MRF/IOMRFZoneList.C)

createIOobject member function of IOMRFZoneList check for the MRFProperties dictionary. If not present, it will display “No MRF models present”.

Checking fvOptions:

#include “createFvOptions.H” (finiteVolume/cfdTools/general/include/createFvOptions.H) has the following content.

fv::options&amp; fvOptions(fv::options::New(mesh));

if (!fvOptions.optionList::size())
{
    Info << "No finite volume options present" << endl;
}

options is a class derived from IOdictionary and optionList classes. optionList is derive from fvList.

SIMPLE loop

Now, we are entering the simple loop. See src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C line 138. Loop exits if criteria is satisfied.

loop()

bool Foam::simpleControl::loop()
{
    solutionControl::setFirstIterFlag(true, true);

    read();

    Time& runTime = const_cast<Time&>(mesh_.time());

    if (initialised_ && criteriaSatisfied())
    {
        Info<< nl
            << algorithmName_
            << " solution converged in "
            << runTime.timeName() << " iterations" << nl << endl;

        // Set to finalise calculation
        runTime.writeAndEnd();
    }
    else
    {
        initialised_ = true;
        storePrevIterFields();
    }

    return runTime.loop();
}

criteriaSatisfied()

bool Foam::simpleControl::criteriaSatisfied()
{
    if (residualControl_.empty())
    {
        return false;
    }

    bool achieved = true;
    bool checked = false;    // safety that some checks were indeed performed

    const dictionary& solverDict = mesh_.solverPerformanceDict();
    forAllConstIters(solverDict, iter)
    {
        const entry& solverPerfDictEntry = *iter;

        const word& fieldName = solverPerfDictEntry.keyword();
        const label fieldi = applyToField(fieldName);

        if (fieldi != -1)
        {
            Pair<scalar> residuals = maxResidual(solverPerfDictEntry);
            checked = true;

            const bool absCheck =
                (residuals.first() < residualControl_[fieldi].absTol);

            achieved = achieved && absCheck;

            if (debug)
            {
                Info<< algorithmName_ << " solution statistics:" << endl;

                Info<< "    " << fieldName << ": tolerance = "
                    << residuals.first()
                    << " (" << residualControl_[fieldi].absTol << ")"
                    << endl;
            }
        }
    }

    return checked && achieved;
}

maxResidual

maxResidual function extract the maximum residual for the specified field. It returns a pair with initial residual as first member, the final residual as second (last) member. Hence, absCheck checks if initial residual is less than the tolerance specified in the residualControl sub-dictionary of SIMPLE dictionary. See a sample SIMPLE dictionary below.

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent      yes;

    residualControl
    {
        p               1e-5;
        U               1e-6;
    }
}

absCheck pass if residual is lesser than corresponding value given in sub-dictionary for all the entries (Here, both for p and U). The for loop forAllConstIters(solverDict, iter) shows the looping.

UEqn.H

Correct Boundary Velocity

MRF.correctBoundaryVelocity(U);

MRFZoneList is a base class which has the above function.

 void Foam::MRFZoneList::correctBoundaryVelocity(volVectorField& U) const
 {
     forAll(*this, i)
     {
         operator[](i).correctBoundaryVelocity(U);
     }
 }

It invokes the same function correctBoundaryVelocity for all the MRFZones. For each MRF zone, it sets the rotating solid body velocity $\overrightarrow{ \Omega} \times \overrightarrow{r}$ on included boundary faces² ³.

void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
{
    if (!active_)
    {
        return;
    }

    const vector Omega = this->Omega();

    // Included patches
    volVectorField::Boundary& Ubf = U.boundaryFieldRef();

    forAll(includedFaces_, patchi)
    {
        const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];

        vectorField pfld(Ubf[patchi]);

        forAll(includedFaces_[patchi], i)
        {
            label patchFacei = includedFaces_[patchi][i];

            pfld[patchFacei] = (Omega ^ (patchC[patchFacei] - origin_));
        }

        Ubf[patchi] == pfld;
    }
}

After this, acceleration and the calculation of the Coriolis force is done. See the simpleFoam wiki page. I skip this as we are not interested in MRF zone right now to keep it simple.

assemble UEqn

tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      + turbulence->divDevReff(U)
     ==
        fvOptions(U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();

apply under relaxation

UEqn.relax();

Apply under relaxation. See page # 554 in Moukalled and simpleFoam wiki page.

constrain the source

fvOptions.constrain(UEqn);

momentum Predictor

If momentum predictor is enabled (it is not necessary):

if (simple.momentumPredictor())
    {
        Info<< "Inside if(simple.momentumPredictor())" << endl;
        solve(UEqn == -fvc::grad(p));

        fvOptions.correct(U);
    }

If momentum predictor is disabled, no implicit calculation of U will be done but U is solved explicitly after the pressure equation. It depends on the importance of the laplacian in the momentum equation. if we have important diffusivity in the momentum, we should solve the momentum equation implicitly.

the momentum predictor is a first approximation of the velocity field obtained from the solution of a momentum equation assembled with the previous time-step pressure field. This velocity doesn’t fulfill the restriction given by the continuity equation which is achieved laterly by the PISO method. Sometimes you can avoid the solution of the momPred saving time, but it can lead to a bad convergence in the PISO loop. Is a user’s choice based in the problem behavior.⁴

Hi, in the PISO method, U is predicted and then corrected in the pressure-velocity coupling loop, so that, if you don’t predict the velocity, the velocity used to start the coupling loop is a derivative of the velocity in the last time step as is done in interFoam and other solvers. In the code you’ve cited the matrix for the momentum equation is assembled so that you have the A and H methods available for the use in the PISO loop, the line: HbyA = rAU*UEqn.H(); calculates a first approximation for U using the H operator applied to the velocity of the last-timestep or which was predicted when you solved the UEqn (if you did it) and the A operator. ⁵

pEqn.H

Calculate rAU

Extract the semi-central coeffs from the UEqn, reverse them and call rAU. This is the famous “operator splitting” ⁶

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

Calculate HbyA

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

HbyA is calculated while satisfying the boundary conditions (this is what constrainHbyA does). See simpleFoam wiki page

flux of HbyA

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

Here the flux function is same as that employed in createFields.H. Input HbyA is a volVectorField and output is a surfaceScalarField.

makeRelative

MRF.makeRelative(phiHbyA);

In a grid which is not deforming and moves with a moving frame of reference the continuity equation has to be satisfied within this grid. This means that the sum of the relative fluxes over the control volumes have to be zero:

void Foam::MRFZone::makeRelative(surfaceScalarField& phi) const
{
    makeRelativeRhoFlux(geometricOneField(), phi);
}
 
template<class RhoFieldType>
void Foam::MRFZone::makeRelativeRhoFlux
(
    const RhoFieldType& rho,
    surfaceScalarField& phi
) const
{
...
    scalarField& phii = phi.primitiveFieldRef();
 
    // Internal faces
    forAll(internalFaces_, i)
    {
        label facei = internalFaces_[i];
        phii[facei] -= rho[facei]*(Omega ^ (Cfi[facei] - origin_)) & Sfi[facei];
    }
 
    makeRelativeRhoFlux(rho.boundaryField(), phi.boundaryFieldRef());
}

adjustPhi

adjustPhi(phiHbyA,U,p);

adjustPhi looks at all the boundary patches of the mesh, adding the mass flux. Some boundaries are adjustable, such as a pressure boundary condition. If the mass balance doesn’t equal zero, it will modify the flux at the adjustable boundary conditions, and fail if it can’t achieve a mass balance of zero. If setPressure is true, it will do all this, if it is false, it doesn’t bother attempting this.It doesn’t touch the internal field. Just boundary patches that can be adjusted, such as zero gradient.⁷,⁸

pEqn in OpenFOAM

The following equation is equivalent to the eqn(21) in first post⁹. However, here H includes source term also.

\begin {equation} u_{P} = \frac{H[u^{*}]}{a_{P}^{*}} - \frac{1}{a_{P}^{*}}\nabla p_{P} \end {equation}

consider this already include the diagonal modification due to under relaxation.

u_{P} = u_{P}^{*}+ u_{P}^{'} \\ p_{P} = p_{P}^{n-1}+ p_{P}^{'}

\begin {equation} u_{P}^{'} = \frac{H[u^{'}]}{a_{P}^{*}} - \frac{\nabla p_{P}^{'}}{a_{P}^{*}} \end {equation}

Note that here it is assumed that the diagonal coefficients of the corrector equation are the same as in the momentum equation previously solved. Taking the divergence of the above equation and setting it to zero (we want to obtain a velocity field $u_{P}$ which satisfies the continuity equation) we get an equation for the pressure:

\begin {equation} \sum_{f} \frac{H[u^{*}]}{a_{P}^{*}}|_{f}\sdot S_{f}+\sum_{f} \frac{H[u^{'}]}{a_{P}^{*}}|_{f}\sdot S_{f} - \sum_{f} \frac{\nabla p_{P}}{a_{P}^{*}}|_{f}\sdot S_{f} = 0 \end {equation}

In SIMPLE, the velocity correction of the neighbouring points are not considered, hence the second term is omitted, to get:

\begin {equation} \sum_{f} \frac{H[u^{*}]}{a_{P}^{*}}|_{f}\sdot S_{f} - \sum_{f} \frac{\nabla p_{P}}{a_{P}^{*}}|_{f}\sdot S_{f} = 0 \end {equation}

\begin {equation} \underbrace{\frac{1}{a_{P}^{*}}}_{\text{rAtU}} \nabla p_{P}|_{f}\sdot S_{f} = \sum_{f} \underbrace{\frac{H[u^{*}]}{a_{P}^{*}}|_{f} \sdot S_{f} }_{\text{phiHbyA}} \end {equation}

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

Here, fvc::div will call fvc::surfaceIntegrate which sums up the phiHbyA in all faces and add to or sunbtract from cell centre according to if the cell is owner or neighbour. For laplacian, see gaussLaplacianScheme(1/3): Theory and implementation in OpenFOAM

constrainPressure

constrainPressure(p, U, phiHbyA, rAtU(), MRF);

This function updates the pressure gradient for the patches where a fixedFluxPressure boundary condition is chosen.

[From simpleFoam wiki page]

For the simple method the velocity at point P is given by (1):

If we dot product with the surface are vector $S_{f}$ (change $u_{P} \text{ by } u_{P} -\Omega \times R$ for MRF case in both of the below two equations ):

\begin {equation} u_{P}\sdot S_{f} =\underbrace{ \frac{H[u^{*}]}{a_{P}^{*}}\sdot S_{f}}_{\text{phiHbyA}}- \underbrace{ \frac{1}{a_{P}^{*}}}_{\text{rAtU}}\nabla p_{P}\sdot S_{f} \end {equation}

Finally we rearrange for the pressure gradient, divide by the magnitude of the surface are vector $S_{f}$ and evaluate at the boundary we get

\begin {equation} \nabla p_{P}\sdot n_{bf} = \frac{1}{|S_{f}| \sdot rAtU} \Big(phiHbyA - u_{P}\sdot S_{f} \Big) \end {equation}

The source code can be found in $FOAM_SRC/finiteVolume/cfdTools/general/constrainPressure/constrainPressure.C:

template<class RAUType, class MRFType>
void Foam::constrainPressure
(
    volScalarField& p,
    const volVectorField& U,
    const surfaceScalarField& phiHbyA,
    const RAUType& rAU,
    const MRFType& MRF
)
{
    constrainPressure(p, geometricOneField(), U, phiHbyA, rAU, MRF);
}
 
template<class RhoType, class RAUType, class MRFType>
void Foam::constrainPressure
(
    volScalarField& p,
    const RhoType& rho,
    const volVectorField& U,
    const surfaceScalarField& phiHbyA,
    const RAUType& rhorAU,
    const MRFType& MRF
)
{
    const fvMesh& mesh = p.mesh();
 
    volScalarField::Boundary& pBf = p.boundaryFieldRef();
 
    const volVectorField::Boundary& UBf = U.boundaryField();
    const surfaceScalarField::Boundary& phiHbyABf =
        phiHbyA.boundaryField();
    const typename RAUType::Boundary& rhorAUBf =
        rhorAU.boundaryField();
    const surfaceVectorField::Boundary& SfBf =
        mesh.Sf().boundaryField();
    const surfaceScalarField::Boundary& magSfBf =
        mesh.magSf().boundaryField();
 
    forAll(pBf, patchi)
    {
        if (isA<fixedFluxPressureFvPatchScalarField>(pBf[patchi]))
        {
            refCast<fixedFluxPressureFvPatchScalarField>
            (
                pBf[patchi]
            ).updateSnGrad
            (
                (
                    phiHbyABf[patchi]
                  - rho.boundaryField()[patchi]
                   *MRF.relative(SfBf[patchi] & UBf[patchi], patchi)
                )
               /(magSfBf[patchi]*rhorAUBf[patchi])
            );
        }
    }
}

We see that fixedFluxPressure boundary condition sets the pressure gradient to the provided value such that the flux on the boundary is that specified by the velocity boundary condition. The pressure gradient is not set directly by the boundaryPatch but is set by the function constrainPressure which is called within the pressure equation.

Non-orthogonal pressure corrector loop

// Non-orthogonal pressure corrector loop
    while (simple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
        );

        pEqn.setReference(pRefCell, pRefValue);

        pEqn.solve();

        if (simple.finalNonOrthogonalIter())
        {
            phi = phiHbyA - pEqn.flux();
        }
    }

setting a reference pressure

 template<class Type>
 void Foam::fvMatrix<Type>::setReference
 (
     const label celli,
     const Type& value,
     const bool forceReference
 )
 {
     if ((forceReference || psi_.needReference()) && celli >= 0)
     {
         source()[celli] += diag()[celli]*value;
         diag()[celli] += diag()[celli];
     }
 }

Used to enforce pressure value in a point. Explained here. Could not understand fully, esp diag()[celli] += diag()[celli] part. Why doubling the diag() at that point ? coeff of p added in matrix and a value is added in source so that p[celli] = value as the result of the solution??

solve the pEqn.

pEqn.solve();

Getting divergence free flux of phi

phi -= pEqn.flux();

Now, we solved for pressure with non-div-free phi. But the literature (Jasak’s thesis, Issa et al.) tells us, that we can correct the fluxes using the pressure field and this way ensure div-free condition. This is the famous phi -= pEqn.flux(); Explained here

Please read Alberto Passalacqua’s post here. Pressure equation is div(phi) = div(U_f) = 0. Must read: fumiya’s explanation of flux()

Explicit relaxation of pressure

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

Momentum corrector

    // Momentum corrector
    U = HbyA - rAtU()*fvc::grad(p);
    U.correctBoundaryConditions();
    fvOptions.correct(U);

Rhie Chaw Interpolation in OpenFOAM

The introduction of neighbour pressure terms in the momentum and pressure equation is done in the spirit of Rhie and Chow. It is achieved by two contributions. Please see ¹⁰:

The first one is the way the Laplacian term is disretized
- Laplacian-term of p uses the value of the gradient of p on the cell face
- The gradient is calculated using neighbouring cells, and not neighbouring faces
- The gradient-term of p is calculated from the cell face values of p
The second contribution is achieved by the flux() operator after the solution of the pressure equation
- phi does not include any effect of pressure when solving the continuity equation
- rUA does not include any effect of pressure when solving for continuity and final velocity corrector

Please read Rhie-Chow Interpolation, Mustafa Bhotvawala ¹¹

Rhie-Chow propose an intermediate velocity where the effect of pressure is removed, which is then interpolated from the cell nodes to cell faces. HbyA is the component of velocity that does not contain pressure. The flux phiHbyA on the RHS is just the flux of the HbyA term. There is no pressure term in the flux. This uncouples the laplacian term in the LHS from any pressure term on the RHS

With this background, read simpleFoam wiki page