IcoFoam

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:

1. J.H. Ferziger and M. Peric Computational Methods for Fluid Dynamics 3rd ed. Springer 2002.
2. Hrvoje Jasak, PhD 1996, PDF of thesis posted at: Error analysis and estimation in the Finite Volume method with applications to fluid flows.
3. 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);
}
 
 
// ************************************************************************* //