• CFD, Fluid Flow, FEA, Heat/Mass Transfer
  • +91 987-11-19-383
  • amod.cfd@gmail.com

New boundary conditions, utility and solver programming in OpenFOAM

OpenFOAM: Programming

Disclaimer:

This content is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OpenFOAM and OpenCFD trademarks.
One of the strengths of OpenFOAM is that new solvers and utilities can be created by its users with some prior knowledge of the underlying method, physics and programming techniques involved. This tutorial by Tommaso Lucchini from Department of Energy Politecnico di Milano serves to generate this understanding.

OpenFOAM code for the N-S equation:

Note that the specifications fvm and fvc are selected by the user from the fvSchemes dictionary in the system dictionary. Here fvm::laplacian means an implicit finite volume discretization for the Laplacian operator, and similarly for fvm::div for the divergence operator. On the other hand, fvc::grad means an explicit finite volume discretization for the gradient operator. The parentheses ( , ) means a product of the enclosed quantities, including tensor products. Similarly, & denotes dot product of two vectors in OpenFOAM. The code to solve Navier-Stokes equation described below is written for demonstration of fvm and fvc class.

Naier-Stokes Equation

Solve (		
  fvm::ddt(rho,U) + fvm::div(U,U) - fvm::laplacian(mu,U) == f - fvc::grad(p)
);		

Similarly, the statement to solve k-equation in OpenFOAM is given by
 solve (
   fvm::ddt(k) +
   fvm::div(phi, k) -
   fvm::laplacian(nu() + nut, k) ==
   nut*magSqr(symm(fvc)::grad(U)) - fvm::Sp(epsilon/k, k)); 
Here k-equation is expressed as

k-equation for turbulence


OpenFOAM can be used to solve equations governed by theory of elasticity (solid mechanics) as well.

Solid mechanics equation

This equation can be solved in OpenFOAM by:
fvVectorMatrix UEqn (	
rho * fvm::d2dt2(U) == fvm::laplacian(2*mu + lambda, U)
    + fvc::div
    (
    - (mu + lambda) * gradU
    + lambda * tr(gradU)*I
    + (sigma & gradU)
    )
);	

OpenFOAM root folder structure

  1. src: the core OpenFOAM libraries with *.C and *H files, including include files. Typing 'foam' on the command line will take you to installation folder. Similarly, app = $WM_PROJECT_DIR/applications = cd $FOAM_APP and util = $WM_PROJECT_DIR/applications/utilities. Command alias | grep "FOAM" retrieves following aliases
    • alias app='cd $FOAM_APP'
    • alias lib='cd $FOAM_LIBBIN'
    • alias run='cd $FOAM_RUN'
    • alias sol='cd $FOAM_SOLVERS'
    • alias src='cd $FOAM_SRC'
    • alias tut='cd $FOAM_TUTORIALS'
    • alias util='cd $FOAM_UTILITIES'
  2. $FOAM_APP/solvers/: solvers (such as icoFoam, pimpleFoam, sonicFoam ...) and utilities (such as blockMesh, setFields, topoSet ...). Solvers are grouped into 12 categoryies: basic, compressible, DNS, financial, incompressible, multiphase, combustion, discreteMethods, electromagnetics heatTransfer, lagrangian and stressAnalysis. Move into the folder solver using command 'sol'
  3. $FOAM_SRC/finiteVolume/finiteVolume: this folder contains discretization schemes namely convectionSchemes, ddtSchemes, fv, fvm, d2dt2Schemes, divSchemes, fvc, fvSchemes, fvSolution, gradSchemes, laplacianSchemes, snGradSchemes. Use alias 'foamfv' to move into $FOAM_SRC/finiteVolume directory.
  4. $FOAM_TUTORIALS: sample cases (simulation data) to demonstrate OpenFOAM functionalities, tutorial folder by typing 'tut' on command line
  5. test: cases that test and validate individual models or aspects of OpenFOAM functionality
  6. etc: scripts -foamyHexMeshDict, mesh, meshQualityDict, postProcessing, setConstraintTypes, solvers, surface
  7. doc: documentation such as user guides in PDF format, Doxygen
  8. In addition to these folders, the root folder also contains a file call 'Allwmake' and a directory called 'wmake'. These are 'make' utility for compilations.
  9. $WM_PROJECT_DIR: Full path to locate binary executables of OpenFOAM release
  10. Use command <find -iname "sampleDict"> to find a file named sampleDict. To search a folder named 'pimpleFoam' in folder '$WM_PROJECT_DIR' use command: find $WM_PROJECT_DIR –type d -name "*pimpleFoam*" where *is known as wildcard character. To find a file replace 'd' with 'f'.
  11. $WM_PROJECT_DIR/src/OpenFOAM/lnInclude: list of *.H and *.c files linked to specific utility

C++ File Structure
  • There needs to be a *.C file with the solver name that calls all of the necessary libraries as well as initializing the time, meshes, and material properties. This file contains the time loop as well as the calls to the individual solvers.
  • The other .H files in the solver directory relate to the initialization and the calls within the time loop.
  • The Make directory includes a file titled options which "contains the full directory paths to locate the header files … and library names" while the file titled "files", lists the entire path and name of the solver once it has been compiled.
  • The difference between compilation of solver and boundary condition lies in the fact that solvers are compiled into executables and boundary conditions are compiled into shared object libraries which can be linked to any solver without need to recompile for every solver.
  • All the source codes related to boundary conditions can be found in the file /src/finiteVolume/fields/fvPatchField/derived. In this folder, you can also check whicl all boundary conditions types are available with the installed version of OpenFOAM.
  • OpenFOAM usage 'wmake' utility for compilation - a reference copy is available in folder '$FOAM_SRC'.
The file structure is shown below.

OpenFOAM solver file heirarchy

FVM solution method

Before further moving into the OpenFOAM program structure, let's look at the governing equations and methods adopted for their numerical solutions.

FVM discretization scheme

The integration over a cell volume need to be converted into surface intergral over the face of the cells using Gauss theorem as explained below.

FVM discretization Gauss Theorem

Let's look at some of the settings in fvScheme dictionary. Note that the variable φ defined above is a generic representation of u, v, w ... and does not represent 'phi' defined below.
ddtSchemes   {
    default         Euler;
}

gradSchemes  {
    default         Gauss linear;
}

divSchemes   {
    default          none;
    div(phi,U)       bounded Gauss linearUpwind grad(U);
    div(phi,k)       bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;       //phi = U.A = Flux field
    div(phi,R)       bounded Gauss upwind;
    div(R)           Gauss linear;
    div(phi,nuTilda) bounded Gauss upwind;       //nuTilda = Spalart-Allmaras model
    div((nuEff*dev2(T(grad(U))))) Gauss linear;  //nuEff = nuLam + nuTurb
}
dev() calculates the deviatoric part of the matrix: ADEV = A - AHYD = A - 1/3 × tr(A)I. OpenFOAMr versions like 4.x onwards, the incompressible solvers will use dev2() [introduced in version 3 for better stability and convergence] otherwise used for compressible flows. The difference between dev2 and dev is the factor two at some position of the code. dev2() calculates the deviatoric part of the matrix: ADEV = A - 2AHYD = A - 2/3 × tr(A)I.
laplacianSchemes {
    default         Gauss linear corrected;
}

interpolationSchemes {
    default         linear;
}
Solvers operate on the lduMatrix class located in folder FOAM_SRC/OpenFOAM/matrices/lduMatrix/ - lduMatrix is a matrix class in which the coefficients of matrix [A] in linear equation [A]{x}={b} are stored as three arrays. One for the lower triangle (l), one for the the diagonal of the matrix (d) and the third array for uppper triangle (u) - hence known as lduMatrix. Inside of lduMatrix.H and lduMatrix.C, the member functions - const scalarField & lower () const, const scalarField & diag () const and const scalarField & upper () const are defined which returns the list of cofficients for the lower, diagonal and upper part of the matrix.

fvSolution: pimpleFoam This dictionary defines the selection of linear solver, smoothing algorithm and convergence tolerance for inner iterative process.

  1. BICCG - Diagonal Incomplete LU preconditioned BiCG solver
  2. diagonalSolver - diagonal solver for both symmetric and asymmetric problems
  3. GAMG - generally known as Geometric-agglomerated Algebraic MultiGrid solver (named Generalised geometric-Algebraic MultiGrid in the user manual)
  4. ICC - Incomplete Cholesky preconditioned Conjugate Gradients solver
  5. PBiCG - Preconditioned bi-conjugate gradient solver for asymmetric lduMatrices using a run-time selectable preconditioner
  6. PCG - Preconditioned conjugate gradient solver for symmetric lduMatrices using a run-time selectable preconditiioner
  7. smoothSolver - Iterative solver using smoother for symmetric and asymmetric matrices which uses a run-time selected smoother
Preconditioner and smoother: Preconditioners improves the convergence of the origina system [A]{x}={b} using a matrix [M] and resulting in a new system [M]-1[A]{x} = [M]-1{b}. Smoother on the other hand is to reduce the dependency of number of iterataions on mesh type.

solvers {
    p     {
        solver           GAMG;
        tolerance        1e-7;
        relTol           0.01;
        smoother         DICGaussSeidel;

    }

    pFinal     {
        $p;
        relTol          0;
    }

    "(U|k|epsilon)"     {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0.1;
    }

    "(U|k|epsilon)Final"     {
        $U;
        relTol          0;
    }
}

PIMPLE {
    nNonOrthogonalCorrectors 0;  
	
	nOuterCorrectors    100; //Outer Loops (Pressure - Momentum Correction )
    nCorrectors         2;   //Inner Loops (Pressure Correction )
} 

Utility statements

  • Centre of faces: Cf[faceI].x() / patchName.Cf()[i][0], Cf[faceI].y() / patchName.Cf()[i][1], Cf[faceI].z() / patchName.Cf()[i][2]
  • Access time: this->db().time().value()
  • Loop over all faces of a patch: forAll(Cf, faceI) { field[faceI] = ... }

Source code: pimpleFoam Large time-step transient solver for incompressible, flow using the PIMPLE (merged PISO-SIMPLE) algorithm. Sub-models include: turbulence modelling, i.e. laminar, RAS or LES, run-time selectable finite volume options, e.g. MRF, explicit porosity. The lines starting with .. are the content of 'included' header files. Some changes in code structure hs been made to reduce number of lines.

\*---------------------------------------------------------------------------*/
/* Make all the FV machinery or templates available to this code: this header
file is located at $FOAM_SRC/finiteVolume/cfdTools/general/include           */
#include "fvCFD.H"

#include "singlePhaseTransportModel.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "IObasicSourceList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])   {
// Check case folder: $FOAM_SRC/OpenFOAM/include
   #include "setRootCase.H"
      ..  #include "listOptions.H"
      ..  Foam::argList args(argc, argv);
      ..   if (!args.checkRootCase()) {
      ..   Foam::FatalError.exit();
      ..  }
      ..  #include "listOutput.H"

//  Create time directories: $FOAM_SRC/OpenFOAM/include
    #include "createTime.H"
      ..  Foam::Info<< "Create time\n" << Foam::endl;
      ..  Foam::Time runTime(Foam::Time::controlDictName, args);

//  Read mesh: $FOAM_SRC/OpenFOAM/include
    #include "createMesh.H"
      ..  Foam::Info
      ..    << "Create mesh for time = "
      ..    << runTime.timeName() << Foam::nl << Foam::endl;
      ..  Foam::fvMesh mesh 
      ..    (Foam::IOobject 
      ..      (Foam::fvMesh::defaultRegion,
      ..       runTime.timeName(),
      ..       runTime,
      ..       Foam::IOobject::MUST_READ  //Must read mesh every timestep
      ..      )
      ..    );

//  Create read U, p, k, epsilon ... dictionaries and initalize fields
//  File is solver specific: located in $FOAM_APP/solvers/incompressible/pimpleFoam	  
    #include "createFields.H"
      ..  Info<< "Reading field p\n" << endl;
      ..  volScalarField p  (      //define p volume scalar field
      ..    IOobject  (
      ..      "p",
      ..      runTime.timeName(),
      ..      mesh,
      ..      IOobject::MUST_READ,
      ..      IOobject::AUTO_WRITE
      ..    ),
      ..   mesh
      ..  );

      ..  Info<< "Reading field U\n" << endl;
      ..  volVectorField U  (        //define U volume vector field
      ..    IOobject (
      ..      "U",
      ..      runTime.timeName(),    // look for dictionary "U"
      ..      mesh,
      ..      IOobject::MUST_READ,
      ..      IOobject::AUTO_WRITE
      ..    ),
      ..   mesh
      ..  );

      ..  #include "createPhi.H"    // create flux field phi at cell faces

      ..  label pRefCell = 0;
      ..  scalar pRefValue = 0.0;
      ..  setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);

      ..  singlePhaseTransportModel laminarTransport(U, phi);

      ..  autoPtr turbulence  (
      ..    incompressible::turbulenceModel::New(U, phi, laminarTransport)
      ..  );
    
    //Declare and initialise the cumulative continuity error.
    //$FOAM_SRC/finiteVolume/cfdTools/general/include	
    #include "initContinuityErrs.H"
      ..  #ifndef initContinuityErrs_H
      ..  #define initContinuityErrs_H
      ..  scalar cumulativeContErr = 0;
      ..  #endif


    pimpleControl pimple(mesh);
	
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    Info<< "\nStarting time loop\n" << endl;  //Time loop as per controlDict
    while (runTime.run())  {                  //While time < endTime
        #include "readTimeControls.H"
        #include "CourantNo.H"

        /*$FOAM_SRC/finiteVolume/cfdTools/general/include:  Reset the timestep 
          to maintain a constant maximum courant Number. Reduction of time-step 
          is immediate, but increase is damped to avoid unstable oscillations */
        #include "setDeltaT.H"
          ..  if (adjustTimeStep) {
          ..    scalar maxDeltaTFact = maxCo/(CoNum + SMALL);
          ..    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
          ..    runTime.setDeltaT (
          ..      min (
          ..        deltaTFact*runTime.deltaTValue(),
          ..        maxDeltaT
          ..      )
          ..    );
          ..  }

        runTime++;

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

        // --- Pressure-velocity PIMPLE corrector loop:  as many outer correctors 
        // as nOuterCorrectors  and as many inner correctors as nCorrectors defined
        // in fvSolution dictionary.
        while (pimple.loop())   {    
            #include "UEqn.H"      //Solution of the momentum equation

            // --- Pressure corrector loop
            while (pimple.correct()) {
                #include "pEqn.H"  //Include the pressure Equation
            }
            if (pimple.turbCorr())  {
                turbulence->correct();
            }
        }
        //Print on the screen information - computational and clock time
        runTime.write();
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }
    Info<< "End\n" << endl;
    return 0;
}
// ************************************************************************* //

Disclaimer:

This content is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OpenFOAM and OpenCFD trademarks.
References
  1. OpenFOAM's basic solvers for linear systems of equations Solvers, preconditioners, smoothers by Tim Behrens, February 18, 2009
  2. OpenFOAM Programmers Guide
  3. Mathematics, Numerics, Derivations and OpenFOAM by Tobias Holzmann at Holzmann CFD
  4. Dictionaries in tutorial folders of OpenFOAM
  5. https://www.youtube.com/watch?v=L94iYGvZr9Q
  6. https://www.youtube.com/watch?v=2FKoMY91SdM&t=3231s
Contact us
Disclaimers and Policies

The content on CFDyna.com is being constantly refined and improvised with on-the-job experience, testing, and training. Examples might be simplified to improve insight into the physics and basic understanding. Linked pages, articles, references, and examples are constantly reviewed to reduce errors, but we cannot warrant full correctness of all content.

>