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