simpleFOAM

Now, let's look at how the simpleFOAM works.

the simpleFOAM.C file:

#include "fvCFD.H"
#include "singlePhaseTransportModel.H" // transport model with viscosity model
#include "RASModel.H" // turbulent model

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
#   include "setRootCase.H"
#   include "createTime.H"
#   include "createMesh.H"
#   include "createFields.H"
#   include "initContinuityErrs.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

#       include "readSIMPLEControls.H"
#       include "initConvergenceCheck.H"

        p.storePrevIter(); //store the previous value of pressure, necessary for under-relaxation

        // Pressure-velocity SIMPLE corrector
        {
#           include "UEqn.H"
#           include "pEqn.H"
        }

        turbulence->correct();

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;

#       include "convergenceCheck.H"
    }

    Info<< "End\n" << endl;

    return 0;
}

the UEqn.H file:

    // Solve the Momentum equation

    tmp UEqn  // U equation difinition; fvm means finite volume matrix, implicit
    (
        fvm::div(phi, U)    // divergence of rho.U.U, for incompressible fluid, phi = U
      + turbulence->divDevReff(U)  // turbulent momentum, equivalent
    );

    UEqn().relax(); // under-relax the U equation

    eqnResidual = solve    // solving; fvc means finite volume ..., explicit
    (
        UEqn() == -fvc::grad(p)
    ).initialResidual();

    maxResidual = max(eqnResidual, maxResidual);

the eqnResidual and maxResidual are defined respectively as:

• eqnResidual: Initial residual of the equation
• maxResidual: Maximum residual of the equations after one solution step

the pEqn.H file:

    p.boundaryField().updateCoeffs(); // Update the boundary conditions for p

    volScalarField AU = UEqn().A(); // Calculate the Ap coefficient
    U = UEqn().H()/AU; // calculate U
    UEqn.clear();
    phi = fvc::interpolate(U) & mesh.Sf(); // Calculate the flux
    adjustPhi(phi, U, p);

    // Non-orthogonal pressure corrector loop
    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
    {
        fvScalarMatrix pEqn  // define the p equation
        (
            fvm::laplacian(1.0/AU, p) == fvc::div(phi)
        );

        pEqn.setReference(pRefCell, pRefValue);

        // Retain the residual from the first iteration
        if (nonOrth == 0)
        {
            eqnResidual = pEqn.solve().initialResidual();
            maxResidual = max(eqnResidual, maxResidual);
        }
        else
        {
            pEqn.solve(); // solve p equation
        }

        if (nonOrth == nNonOrthCorr)
        {
            phi -= pEqn.flux();  // correct the flux
        }
    }

#   include "continuityErrs.H"  // Calculate continuity errors

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

    // Momentum corrector
    U -= fvc::grad(p)/AU;
    U.correctBoundaryConditions();

References
[1] The SIMPLE algorithm in OpenFOAM. http://openfoamwiki.net/index.php/The_PISO_algorithm_in_OpenFOAM, May 2, 2011

   Send article as PDF