icoFoam solves the incompressible laminar Navier-Stokes equations using the PISO algorithm. The code is inherently transient, requiring an initial condition (such as zero velocity) and boundary conditions. The icoFOAM code can take mesh non-orthogonality into account with successive non-orthogonality iterations. The number of PISO corrections and non-orthogonality corrections are controlled through user input.

A highly annotated version of the code is shown below. Ferziger and Peric, below, give a fairly standard overview of PISO, but OpenFOAM uses a notation and form closer to Jasak and Rusche's theses. References to PISO documentation are:

- J.H. Ferziger and M. Peric
*Computational Methods for Fluid Dynamics*3rd ed. Springer 2002. - Hrvoje Jasak, PhD 1996, PDF of thesis posted at: Error analysis and estimation in the Finite Volume method with applications to fluid flows.
- Henrik Rusche, PhD 2002 Computational fluid dynamics of dispersed two-phase flows at high phase fractions.

int main(int argc, char *argv[]) { // these header files contain source code for common tasks that are used in numerous applications. # include "setRootCase.H" # include "createTime.H" # include "createMesh.H" # include "createFields.H" # include "initContinuityErrs.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; // use the runTime object to control time stepping for (runTime++; !runTime.end(); runTime++) { Info<< "Time = " << runTime.timeName() << nl << endl; # include "readPISOControls.H" # include "CourantNo.H" //set up the linear algebra for the momentum equation. The flux of U, phi, is treated explicity //using the last known value of U. fvVectorMatrix UEqn ( fvm::ddt(U) + fvm::div(phi, U) - fvm::laplacian(nu, U) ); // solve using the last known value of p on the RHS. This gives us a velocity field that is // not divergence free, but approximately satisfies momentum. See Eqn. 7.31 of Ferziger & Peric solve(UEqn == -fvc::grad(p)); // --- PISO loop---- take nCorr corrector steps for (int corr=0; corr<nCorr; corr++) { // from the last solution of velocity, extract the diag. term from the matrix and store the reciprocal // note that the matrix coefficients are functions of U due to the non-linearity of convection. volScalarField rUA = 1.0/UEqn.A(); // take a Jacobi pass and update U. See Hrv Jasak's thesis eqn. 3.137 and Henrik Rusche's thesis, eqn. 2.43 // UEqn.H is the right-hand side of the UEqn minus the product of (the off-diagonal terms and U). // Note that since the pressure gradient is not included in the UEqn. above, this gives us U without // the pressure gradient. Also note that UEqn.H() is a function of U. U = rUA*UEqn.H(); // calculate the fluxes by dotting the interpolated velocity (to cell faces) with face normals // The ddtPhiCorr term accounts for the divergence of the face velocity field by taking out the // difference between the interpolated velocity and the flux. phi = (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, U, phi); // adjusts the inlet and outlet fluxes to obey continuity, which is necessary for creating a well-posed // problem where a solution for pressure exists. adjustPhi(phi, U, p); // iteratively correct for non-orthogonality. The non-orthogonal part of the Laplacian is calculated from the most recent // solution for pressure, using a deferred-correction approach. for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { // set up the pressure equation fvScalarMatrix pEqn ( fvm::laplacian(rUA, p) == fvc::div(phi) ); // in incompressible flow, only relative pressure matters. Unless there is a pressure BC present, // one cell's pressure can be set arbitrarily to produce a unique pressure solution pEqn.setReference(pRefCell, pRefValue); pEqn.solve(); // on the last non-orthogonality correction, correct the flux using the most up-to-date pressure if (nonOrth == nNonOrthCorr) { //The .flux method includes contributions from all implicit terms of the pEqn (the Laplacian) phi -= pEqn.flux(); } } // end of non-orthogonality looping # include "continuityErrs.H" // add pressure gradient to interior velocity and BC's. Note that this pressure is not just a small // correction to a previous pressure, but is the entire pressure field. Contrast this to the use of p' // in Ferziger & Peric, Eqn. 7.37. U -= rUA*fvc::grad(p); U.correctBoundaryConditions(); } // end of the PISO loop runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } // end of the time step loop Info<< "End\n" << endl; return(0); } // ************************************************************************* //