WRITE YOUR OWN SOLVER

New boundary conditions, utility and solver programming in OpenFOAM

OpenFOAM: Programming

Disclaimer:

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 is a C++ library used primarily to create executable, known as applications. The applications are further categorized into two groups: 'solvers' which are each designed to solve a specific problem in continuum mechanics and 'utilities' to perform tasks that involve data manipulation. Wikki Ltd is provider and major contributor to foam-extend, the community version of the OpenFOAM software. A solver needs to be compiled into a binary file whereas boundary conditions need to be compiled into a shared obeject library. The shared object library can be dynamically linked to any solver without re-compiling.

OpenFOAM code for the N-S equation:

OpenFOAM like any other numerical solver, converts the Partial Differential Equations (PDEs) into a set of linear algebraic equations, [A] {x} = {b}, where {x} and {b} are volFields (geometricField). [A] is a fvMatrix which is created by the discretization of a geometricField and inherits the algebra of its corresponding field. It supports many of the standard algebraic matrix operations.

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

• fvc::ddt(A) = ∂A/∂t
• fvc::ddt(rho, A) = ∂(ρA)/∂t
• fvc::d2dt2(rho, A) = ∂(ρ.∂A/∂t)/∂t
• fvc::laplacian(s, A) = ∇(s.∇A)
• fvc: explicit evaluation of derivatives (i.e. based on known values). fvm: implicit evaluation of derivatives

  solve(fvm::ddt(phi) == kappa*fvm::laplacian(phi))

  solve(fvm::ddt(phi) == kappa*fvc::laplacian(phi))

  solve (fvm::ddt(phi) == kappa*0.5*(fvm::laplacian(phi) + fvc::laplacian(phi)))

solve (
   fvm::ddt(k) + fvm::div(phi, k) - fvm::laplacian(nu() + nut, k) ==
   nut*magSqr(symm(fvc)::grad(U)) - fvm::Sp(epsilon/k, k)
); 

Pressure-based Equation

Mass conservation:

solve (
    fvm::ddt(rho)
  + fvc::div(phi)
); 

Momentum equation:

fvVectorMatrix UEqn  (
    fvm::ddt(rho, U)
  + fvm::div(phi, U)
  - fvm::laplacian(mu, U)
);
solve(UEqn == -fvc::grad(p)); 

Energy equation:

solve (
    fvm::ddt(rho, e)
  + fvm::div(phi, e)
  - fvm::laplacian(mu, e)
    ==
  - p*fvc::div(phi/fvc::interpolate(rho))
  + mu*magSqr(symm(fvc::grad(U)))
);
T = e/Cv; 

Density-based compressible flow solver

solve (
    fvm::ddt(rho)
  + fvm::div(phiv, rho)
);

solve (
    fvm::ddt(rhoU)
  + fvm::div(phiv, rhoU)
    ==
  - fvc::grad(p)
);

solve (
    fvm::ddt(rhoE)
  + fvm::div(phiv, rhoE)
    ==
  - fvc::div(phiv2, p)
);

fvVectorMatrix UEqn (	
  rho * fvm::d2dt2(U) == fvm::laplacian(2*mu + lambda, U)
   + fvc::div (
     - (mu + lambda) * gradU  + lambda * tr(gradU)*I
     + (sigma & gradU)
  )
);

fvScalarMatrix TEqn (
  fvm::ddt(T) + fvm::div(phi,T) - fvm::laplacian(DT,T)  // = 0
);
Teqn.solve();

Can you write a solver which supercedes all those defined separately for single phase incompressible flows?

Syntax and Conventions

A short and detailed presentation on C++ programming methods and its implementation in OpenFOAM can be found here.

C++ is a strongly type language: checks for types of all data and functions and gives compiler errors if they do not fit. C++ is a statically typed [that error checking is done only during compilation], compiled [a separate application called compiler is needed and executable needs to be created], case-sensitive, free-form [unlike FORTRAN or Python where indentation specifies a block of code it does not recognize end of the line as a terminator . Thus, it does not matter where a statement in a line is put] programming language that supports procedural, object-oriented programming. An interpreted language such as Python may also be object-oriented but error checking is done during run-time and no compilation is required. C++ has backward compatibility to C, that is the programs written in C can be compiled in C++ environment without any modification. The Java Virtual Machine and Python interpretor are written in C++

github.com/ UnnamedMoose/ BasicOpenFOAMProgrammingTutorials provides a good starting point for those looking to develop apps in OpenFOAM. Some codes from this repository has also been included in later section of this page.

Ubuntu on WSL (Windows Sub-system for Linux) uses g++ to compile C++ codes. Use command `whereis g++` to check if g++ is installed or not. If not, use command `sudo apt-get install g++` to install the compiler. g++ hello.cpp: to complie the code, use command ./a.out to run the compiled code (a.out is the default name of executable). g++ hello.cpp -o hello: to name the executable 'hello'. Use ./hello to run the executable created.

C++ is a family of programming language designated as "Object Oriented". Object orientation focuses on the 'objects' instead of the functions. An 'object' belongs to a class of objects with the same features or attributes or properties. For example, 'student' can be a 'class' of objects and roll number, student name, grade ... would be the attributes. The class of objects defines:

• The construction of the object
• The destruction of the object
• Attributes of the object (member data)
• Methods which manipulate the object (member functions)

class foamEqn {
  public:  //declaration of public members
    int m;
  private: //declaration of hidden members
}

#include <iostream>
void fileNames(string fname) {
  cout << fname << "*.txt\n";
}

int main() {
  fileNames("mat");
  fileNames("bc");
  return 0;
} 

The :: operator - It is called the scope-resolution operator and by prefixing a type-name with this, it tells compiler to look in the global namespace for the type. In other words, :: links either of a variable, function, class, typedef... to a namespace or to a class. If there is no left hand side before :: then it refers to the global namespace. Thus, in 'fvm::laplacian', fvm is a namespace and laplacian is a function (it can also be another namespace of a variable or an onject). Foam: namespace for OpenFOAM, constant: namespace for collection of constants, fv: namespace for finite-volume, fvc (Finite Volume Calculus): namespace of functions to calculate explicit derivatives - given a field creats a new field, fvm (Finite Volume Method): namespace of functions to calculate implicit derivatives returning a matrix.

• Statements end with a semi-colon, blocks of code such as those in loops are enclosed in curly braces { ... }
• There are essentially three uses of the preprocessor in C++: directives, constants, and macros.
• Directives are commands that tell the preprocessor to skip part of a file, include another file or define a constant or macro. Directives always begin with a hash sign (#) and for readability they are placed flush to the left of the sentence.
• All other uses of the preprocessor involve processing (#)defined constants or macros. Typically, constants and macros are written in ALL CAPITAL LETTERS to indicate they are special.
• Conditional Compilation: options that tell preprocessor to remove lines of code before compiler take over the files. They include #if, #elif, #else, #ifdef, and #ifndef. An #if or #if/#elif/#else block or a #ifdef or #ifndef block must be terminated with a closing #endif.
• #include: this directive tells the preprocessor to grab the text of a file and place it directly into the current file.
• #ifndef: this directive stands for IF Not DEFined
• #define: define a constant. e.g. #define PI (3.14156)
• OpenFOAM screen output is very similar to rudimentary C++ with its std::cout, std::nl and std::endl replaced by Foam::Info, Foam::nl, and Foam::endl
• << and >> are output and input operators and 'endl' is a manipulator that generates a new line (equivalent to \n). Instead of using 'cin' and 'cout', OpenFOAM defined a new output stream 'Info' and it is recommended to use this one instead since it takes care of write-outs for parallel simulations.
• namespace is a method to distinguish two variables or functions having same names. It is analogous to roll numbers to differentiate two students having exactly same name.
• << can also be used to combine text and variable though '+' operator is used to concatenate two strings.
• OpenFOAM has implemented (inherited) few variables to suit the programming requirement.
  • label = integer (int in C++): the size of variable types in bytes are not explicitly defined in C++
  • scalar = float or double
  • word = single word (inherited rom C++ string object)
  • string = quoted text
  • list = array bounded by () e.g. (1 2 4 8)
  • prefixList = array of type 'prefix'e.g. labelList, wordList, scalarList
  • 'Info': similar to 'cout', but also works in parallel applications. Example: Info << "Some message.... " << endl;
• using namespace Foam::constant → The 'using' directive here is used to refer to a particular item within a namespace.
• Foam::scalar Foam::noiseFFT::p0 = 1e-5 → Here the variable p0 in scopes 'Foam::scalar' and 'Foam::noiseFFT' are assigned a value of 1e-5.
• typedef: create a new name for an existing type. e.g. typedef int single; → single is another name for int, typedef float real → real is another name for float. Note that the term 'void' represents the absence of type. In case of a 'function' it tells that the function does not return anything.
• Constants refer to fixed values that the program cannot alter, they are also called literals. Constants can be defined either using #define preprocessor or 'const' keyword.
• x = &myObj defines x as pointer to object myObj. & implies object is a reference i.e. it refers to an object of the same class. x->read() accesses the member function read(). e.g. turbulence->correct() in $FOAM_SOLVERS/ incompressible/ pimpleFoam/ pimpleFoam.C
• -> is used to call a member data (or function) through a pointer. OpenFOAM uses 'pointer' functionality to make run-time choices (such as selection of turbulence model) possible.
• A '*' after an object name implies that a pointer should be constructed.
• 'forAll' defined in UList.H is an alias to 'for-loop' defined to march through all entries of a list of objects of any class.
• If an argument variable should be changed inside a function, the type of the argument must be a reference i.e. void change(double& x1). The reference parameter x1 here will be a reference to the argument to the function instead of a local variable in the function.

Operators: The usual precedence-rules apply to all of the following.

• &: Inner product for vectors and tensors, ^: Cross product of two vectors
• %: Modulo operator
• && ||: The logical and and or operators, - : Unary negation, ! : Logical negation
• < > >= <=: Relational comparisons, == != : : Equality and inequality-operators
• ? and : Ternary operator for cond ? a : b. The condition must evaluate to bool, the fields a and b must be of the same type.
• The keyword bool can be used to force a boolean conversion of scalar values. A threshold of -/+ 0.5 is used to define true/false.
• pos(x): if x > 0: 1.0 else 0.0, pos0(x) : if x ≥ 0: 1.0 else 0.0
• neg(x): if x < 0: 1.0 else 0.0
• neg0(x): if x ≤ 0: 1.0 else 0.0
• sign(x): if x is positive: 1.0 else if x is negative: -1.0
• mag(x): Absolute value |x| Implemented for all data types. Yields a scalar. Can be also be used to convert a bool to a scalar.
• magSqr(x): Square of the magnitude |x|^2 Implemented for all non-logical data types. Yields a scalar
• pow(x,y): Power x^y, exp(x): Exponential function e^x, log(x): Natural logarithm, log10(x): Logarithm base 10, sqr(x): Square x^2, sqrt(x): Square root, cbrt(x) : Cubic root

• 'polyMesh' class holds mesh definition objects whereas polyBoundaryMesh is a list of polyPatches
• To access fields defined at cell centers of the mesh you need to use the class volField which can be accessed by adding the header volFields.H

volScalarField p (   //Create the scalar volume field
  IOobject (  //Assign and initialize volField to mesh
    "p", runTime.timeName(), mesh,
    IOobject::MUST_READ, IOobject::AUTO_WRITE
  ),
  mesh
); 

• volScalarField field of scalars (eg. pressure)
• volVectorField field of vectors (eg. velocity)
• volTensorField field of 2nd rank tensors (eg. stress)

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. One can add his/her own aliases such as: alias ubc='cd $FOAM_SOURCE/finiteVolume/fields/fvPatchFields/derived' - to move to folder containing all boundary conditions.
3. bin: The bin directory contains shell scripts, such as paraFoam,foamNew, foamLog...
4. Other userful aliases can be: alias fv='cd $FOAM_SRC/finiteVolume', alias usrSrc='cd $WM_PROJECT_USER_DIR/src', alias fo='cd $FOAM_SRC/OpenFOAM/db/functionObjects'. Note there should be no whitespace character between 'useSrc' and '=' and make sure that the folder exists.
5. $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'
6. $FOAM_SRC/finiteVolume/finiteVolume: this folder contains discretization schemes namely convectionSchemes, ddtSchemes, fv, fvm, d2dt2Schemes, divSchemes, fvc, fvSchemes, fvSolution, gradSchemes, laplacianSchemes, snGradSchemes. Use alias 'fv' to move into $FOAM_SRC/finiteVolume directory.
7. $FOAM_SRC/OpenFOAM: This core library includes the definitions of the containers used for the operations, the field definitions, the declaration of the mesh and of all the mesh features such as zones and sets
8. $FOAM_SRC/TurbulenceModel: contains the source code of turbulence models
9. $FOAM_APPBIN: The path to the folder containing the executables. Note that the $FOAM_APP refers to the folder containing the source codes.
10. $FOAM_TUTORIALS: sample cases (simulation data) to demonstrate OpenFOAM functionalities, tutorial folder by typing 'tut' on command line
11. test: cases that test and validate individual models or aspects of OpenFOAM functionality. Sample codes to help understand the usage of OpenFOAM libraries
12. etc: The etc directory contains environment set-up files, global controlDict and default thermoData. scripts -foamyHexMeshDict, mesh, meshQualityDict, postProcessing, setConstraintTypes, solvers, surface
13. doc: documentation such as user guides in PDF format, Doxygen
14. 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.
15. $WM_PROJECT_DIR: Full path to locate binary executables of OpenFOAM release
16. 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'. Similarly, command "find $FOAM_SRC -iname \*fixedValue\*.[CH]" lists all the occurrences of fixedValue.C and fixedValue.H files in folder $FOAM_SRC.
17. $WM_PROJECT_DIR/src/OpenFOAM/lnInclude: list of *.H and *.c files linked to specific utility.
18. $WM_PROJECT_USER_DIR is the user's directory for running OpenFOAM cases.
19. All the source codes related to boundary conditions can be found in the src files of respective boundary condition in finiteVolume/ fields/ fvPatchField/ derived under $FOAM_SRC. In this folder, you can also check which all boundary conditions types are available with the installed version of OpenFOAM.
20. Environment variables pointing to local (user's) workspace applications and libraries: $FOAM_PROJECT_USER_DIR, $FOAM_USER_APPBIN and $FOAM_USER_LIBBIN

export WM_OPTIONS=$WM_ARCH$WM_COMPILER$WM_PRECISION_OPTION$WM_LABEL_OPTION$WM_COMPILE_OPTION

• export WM_COMPILER=Gcc
• export WM_PRECISION_OPTION=DP
• export WM_COMPILE_OPTION=Opt

• export WM_ARCH=$(uname -s)
• export WM_LABEL_OPTION=Int$WM_LABEL_SIZE

File Structure OpenFOAM

• Codes separation into files: *.C and *H
• Comment and indentation style such 4 characters for indentation
• Common approach to common needs e.g. I/O, construction of objects, stream support, handling function parameters, const and non-const access
• Blank lines, no trailing whitespace, no spaces around brackets
• Porting and multiple platform support handled through environment variables
• All tools, compiler versions and paths can be controlled with environment variables

• OpenFOAM follows object oriented programming where the types (int, double) can be seen as classes and the variables assigned to a type are objects of that class (int a). Object orientation focuses on the objects instead of the functions. An object belongs to a class of objects with the same attributes. The class defines:
  • The construction of the object
  • The destruction of the object
  • Attributes of the object (member data)
  • Methods which manipulate the object (member functions)
• 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. Make/files: names all the source files (.C), it specifies the boundary condition library name and location of the output file.
• 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.
• OpenFOAM usage 'wmake' utility for compilation - a reference copy is available in folder '$FOAM_SRC'.
• OpenFOAM represents scalars, vectors and matrices as tensor fields. A zero rank tensor is a scalar, a first rank tensor is a vector and a second rank tensor is a (3 × 3) matrix. OpenFOAM contains a C++ class library named primitive ($FOAM_SRC/OpenFOAM/primitives/) which contains classes for the tensor mathematics.
• If a tensor is defined as "tensor T(1, 2, 3, 4, 5, 6, 7, 8, 9)" then T.yz() ≡ T(2,3) = 6.

  XiFoam.C

  EXE = $(FOAM_APPBIN)/XiFoam → instructs what the new executable name will be

  LIB = $(FOAM_USER_LIBBIN)/libAcousticAnalogy → instructs what the new library name will be

The compiler searches for the included header files in the following order, specified with the -I option in 'wmake':

1. $WM_PROJECT_DIR/src/OpenFOAM/lnInclude directory
2. local lnInclude directory
3. local directory
4. platform dependent paths set in files in the $WM_PROJECT_DIR/wmake/rules/$WM_ARCH/ directory. e.g. /usr/X11/include and $(MPICH_ARCH_PATH)/include
5. other directories specified explicitly in the Make/options file with the -I option

The Make/options file contains the full directory paths and library names using the syntax: EXE_LIBS = -L<libraryPath> -l<libraryName> The directory paths are preceeded by the -L flag, the library names are preceeded by the -l (small L) flag.

Make/options: specifies directories to search for include files and libraries to link the solver against.

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/sampling/lnInclude

EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools \
    -lsampling

• EXE: Location for the executable. Use $FOAM_APPBIN or $FOAM_USER_APPBIN
• LIB: location for the library. Use $FOAM_LIBBIN or $FOAM_USER_LIBBIN
• EXE_INC: location of include paths. Use $LIB_SRC. Each library soft-likts all files into lnInclude directory for easy inclusion of search paths
• EXE_LIBS, LIB_LIBS: Link libraries for executables or other libraries

Note:

V3.x / V5.x / V7.x

EXE_INC = \
  -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
  -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
  -I$(LIB_SRC)/transportModels \
  -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
  -I$(LIB_SRC)/fvOptions/lnInclude \
  -I$(LIB_SRC)/finiteVolume/lnInclude \
  -I$(LIB_SRC)/meshTools/lnInclude \
  -I$(LIB_SRC)/sampling/lnInclude

EXE_LIBS = \
  -lturbulenceModels \
  -lincompressibleTurbulenceModels \
  -lincompressibleTransportModels \
  -lfiniteVolume -lmeshTools -lfvOptions -lsampling 

V8.0

EXE_INC = \
  -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
  -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
  -I$(LIB_SRC)/transportModels/lnInclude \
  -I$(LIB_SRC)/finiteVolume/lnInclude \
  -I$(LIB_SRC)/meshTools/lnInclude \
  -I$(LIB_SRC)/sampling/lnInclude

EXE_LIBS = \
  -lmomentumTransportModels \
  -lincompressibleMomentumTransportModels \
  -lincompressibleTransportModels \
  -lfiniteVolume -lmeshTools -lfvOptions -lsampling 

Following <optional Arguments> may be specified for building libraries. Argument Type of compilation are:

• lib: Build a statically-linked library
• libso: Build a dynamically-linked library
• libo: Build a statically-linked object file library
• jar: Build a JAVA archive
• exe: Build an application independent of the specified project library

wclean

On execution, wmake builds a dependency list file with a .dep file extension and a list of files in a Make/$WM_OPTIONS directory. To remove these files the user can run the wclean script. This script located in folder $WM_DIR cleans up the wmake control directory Make/$WM_OPTIONS and removes the 'lnInclude' directories generated for libraries.

By design, OpenFOAM is intended to be used inside Linux OS and a sound knowledge of this operating system will help one manoeuvre the folders easily. E.g. cp -r is a command required to copy a folder. -r option stands for recursive copy. Similarly, $(!!) can be used to access the output from last command. 'tree' command list the content of current folder in a hierarchical tree structure: `tree -d -L 1 $WM_PROJECT_DIR`

FVM solution method

 volVectorField cen = mesh.C();
 scalar sumRT = 0.0;
 forAll(cen, k) { 
   sumRT = sumRT + (mesh.V()[k] * rho[k] * T[k]));
 }

fvVectorMatrix UEqn (
 fvm::ddt(U) + 0.5*fvm::div(phi, U) - fvm::laplacian(nu, U)
); UEqn.solve();

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   {               //ddtScheme is a 'class' under namespace 'fv'
  default         Euler;   
}

gradSchemes  {               //gradScheme is a 'class' under namespace 'fv'
  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
}

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

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 )
} 

To access fields defined at cell centers of the mesh you need to use the class volField. Information about the mesh: Functions that provide information about the mesh and are used without any arguments.

Utility statements

• runTime.write(): write to stadard output
• Centre of faces: Cf[faceI].x() / patchName.Cf()[i][0], Cf[faceI].y() / patchName.Cf()[i][1], Cf[faceI].z() / patchName.Cf()[i][2]. E.g. volVectorField centres = mesh.C() creates a vector field names 'centres' which holds X, Y and Z coordinates of the centroids of the mesh cells.
• Access time: this->db().time().value()
• Loop over all faces of a patch: forAll(Cf, faceI) { field[faceI] = ... } Note that forAll is a OpenFOAM macro and can be replaced with standard for loop in C++ with new iterator instead of iterator 'patchI' or 'faceI'. forAll(Cf, faceI) = for (int faceI=0; Cf.size()<faceI; faceI++)
• forAll(U, i) is equivalent to "for (int i=0; patch.size()<i; i++)". zpatch.Cf()[i][1] - indices [0], [1] and [2] are used to access the X-, Y- and Z-coordinates respectively.
• volVectorField cen = mesh.C() - create a variable of class 'volVectorField' named 'cen' which holds the centroids of all cells in the mesh.
• surfaceVectorField srf = mesh.Sf() - create a variable of class 'surfaceVectorField' named 'srf' which holds the surface area vector of all faces in the mesh.
• surfaceScalarField smg = mesh.magSf() - create a variable of class 'surfaceScalarField' named 'smg' which holds the magnitude of surface areas of all faces in the mesh.
• surfaceVectorField scn = mesh.Cf() - create a variable of class 'surfaceVectorField' named 'scn' which holds the centroids of all faces in the mesh.
• surfaceScalarField ffx = Phi() - create a variable of class 'surfaceScalarField' named 'ffx' which holds the flux through all the faces the mesh. In other words, phi = mesh.Sf() & Uf where & denotes the dot product.

forAll(Cf, faceI) {
 if((Cf[faceI].z() > Z0) && (Cf[faceI].z() < Z1) &&
  (Cf[faceI].y() > Y0) && (Cf[faceI].y() < Y1)  ) {
  if ( t < 10.0) {
    field[faceI] = 1.0;
  }
  else {
    field[faceI] = 5.0;
  }
 }
}

Code with forAll loop

forAll(mesh.boundaryMesh(), patchI)
Info << "Patch " << patchI << ": " << mesh.boundary()[patchI].name() << " with "
<< mesh.boundary()[patchI].Cf().size() << " faces. Starts at total face "
<< mesh.boundary()[patchI].start() << endl; 

Code with standard C++ loop

for (int i = 0; i < mesh.boundaryMesh().size(); i++)
Info << "Patch " << i << ": " << mesh.boundary()[i].name() << " with "
<< mesh.boundary()[i].Cf().size() << " faces. Starts at total face "
<< mesh.boundary()[i].start() << endl; 

Header Files

createFields.H: declares all the field variables and initializes the solution. Contains field and material property definition used by the equation.

volScalarField T (
  IOobject (
  "T", runTime.timeName(), mesh, 
  IOobject::MUST_READ, 
  IOobject::AUTO_WRITE
  ), mesh
 ); 

In OOP like C++ used in OpenFOAM, a dot operator is used to access methods and attributes of an object. For example, piple.loop() refers to loop() method defined in object 'pimple'. Similarly, runTime.x() defined below where 'x' stands for attributes of object runTime.

• When the object is created or asked to read: MUST_READ, READ_IF_PRESENT, NO_READ
• When the object is destroyed or asked to write: AUTO_WRITE, NO_WRITE
• runTime is a variable of OpenFOAM class Time: used to control the time-stepping through the code

IOdictionary transportProperties (  
 IOobject (
   "transportProperties",
   runTime.constant(), mesh,
   IOobject::MUST_READ,
   IOobject::NO_WRITE
 )
);
dimensionedScalar DT( transportProperties.lookup("DT") ); 

In the following example NO_READ implies data does not need to be read from any dictionary: rather – construct using ::mag() function specified inside the code. Then write it out to the terminal windows.

volumeScalarField magU (
IOobject (
  "magU", runTime.timeName ( ),  mesh,
  IOobject::NO_READ,
  IO::AUTO_WRITE
  ),
  ::mag(U)
) ;
magU.write(); 

Source code: pimpleFoam - description as in V6 is Transient solver for incompressible, turbulent flow of Newtonian fluids, with optional mesh motion and mesh topology changes. Turbulence modelling is generic i.e. laminar, RAS or LES may be selected The lines starting with .. are the content of 'included' header files. Some changes in code structure hss been made to reduce number of lines. The source code can be found in the folder $FOAM_APP/ solver/ incompressible / pimpleFoam

• A( ): returns the central coefficients of a fvVectorMatrix
• H( ): returns the H-operation source of a fvVectorMatrix.

  

• Sf(): returns cell face area vector of a fvMesh
• flux(): returns face flux field from a fvScalarMatrix
• correctBoundaryConditions(): correct the boundary field of volVectorField.
• adjustPhi(): Adjust the inlet and outlet fluxes to ensure continuity (mass balance)
• constrainPressure(): Update pressure boundary conditions to ensure flux consistency

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

/*fvCFD.H will add all the class declarations needed to access mesh, fields, 
tensor algebra, fvm/fvc operators, time, parallel communication and so on...*/
#include "fvCFD.H"

#include "dynamicFvMesh.H"
#include "singlePhaseTransportModel.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "CorrectPhi.H"
#include "fvOptions.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main part of the program, argc and argv prodives command line options
int main(int argc, char *argv[])   {
  
  // Check and set root case folder: $FOAM_SRC/OpenFOAM/include -Checks the 
  // basic folder structure, verifies there is a control dict present ...
  // also deals with parsing command line arguments and options.
  #include "setRootCase.H"
  
  #include "listOptions.H"
  Foam::argList args(argc, argv);
    if (!args.checkRootCase()) {
      Foam::FatalError.exit();  
	//exit() is a member function of class 'FatalError' in namespace 'Foam'
    }
  #include "listOutput.H"

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

  // Read mesh: $FOAM_SRC/OpenFOAM/include
  // Create fvMesh (instances of objects (or classes) called mesh). mesh.C() 
  // and .Cf() return vector fields denoting centres of each cell and internal
  // face. Calling the mesh.C().size() method yields total size of the mesh.
  #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
          )
        );

  //Initialize field: read 0/U, 0/p, 0/k ... transportPropties... dictionaries
  //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",                   //name of the dictionary, field
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,   //must read 'p' field every timestep
        IOobject::AUTO_WRITE
      ),
      mesh
    );

    Info<< "Reading field U\n" << endl;
    volVectorField U  (        //define U volume vector field
      IOobject (
        "U",                   //name of the dictionary
        runTime.timeName(),    // look for dictionary "U"
        mesh,
        IOobject::MUST_READ,   //must read 'U' field every timestep
        IOobject::AUTO_WRITE
      ),
     mesh
    );

    #include "createPhi.H"    // Flux field phi of velocity U at cell faces

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

    singlePhaseTransportModel laminarTransport(U, phi);

    autoPtr<incompressible::turbulenceModel> 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"  // Read adjustTimeStep, maxCo, maxDeltaT
      
      #include "CourantNo.H"         // Calculate Mean and Max Courant Numbers 

      /*$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"         // Get the deltaT defined in controlDict
      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())   {    
        
		//Set-up the linear algebra for momentum equation and solve. UEqn.H 
		//corresponds to the momentum predictor, and pEqn.H corresponds to the
		//corrector step(s).
		#include "UEqn.H"

        // --- 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;
}
// ************************************************************************* //

fvVectorMatrix UEqn (
  fvm::ddt(U) + fvm::div(phi, U) - fvm::laplacian(nu, U)
);
if (piso.momentumPredictor()) {
  solve(UEqn == -fvc::grad(p));
} 

while (piso.correct()) {
  volScalarField rAU(1.0/UEqn.A());
  volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
  surfaceScalarField phiHbyA (
    "phiHbyA",
    fvc::flux(HbyA) + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
  );
  adjustPhi(phiHbyA, U, p);
  // Update the pressure BCs to ensure flux consistency
  constrainPressure(p, U, phiHbyA, rAU);
  while (piso.correctNonOrthogonal()) {
    fvScalarMatrix pEqn (
      fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
    );
    pEqn.setReference(pRefCell, pRefValue);
    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
    if (piso.finalNonOrthogonalIter()) {
      phi = phiHbyA - pEqn.flux();
    }
  }
  #include "continuityErrs.H"
  U = HbyA - rAU*fvc::grad(p);
  U.correctBoundaryConditions();
} 

Create time

Create mesh for time = 0

PIMPLE: No convergence criteria found

PIMPLE: No corrector convergence criteria found
        Calculations will do 5 corrections

PIMPLE: Operating solver in transient mode with 5 outer correctors

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type LES
Selecting LES turbulence model SpalartAllmarasDDES
Selecting LES delta type vanDriest
Selecting LES delta type cubeRootVol
--> FOAM Warning :
    From function void Foam::LESModels::cubeRootVolDelta::calcDelta()
    in file LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C at line 55
    Case is 2D, LES is not strictly applicable

Selecting patchDistMethod meshWave
LES
{
    LESModel        SpalartAllmarasDDES;
    turbulence      on;
    printCoeffs     on;
    delta           vanDriest;
    vanDriestCoeffs
    {
        delta           cubeRootVol;
        cubeRootVolCoeffs
        {
            deltaCoeff      2;
        }
    }
    Ce              1.048;
    sigmaNut        0.66666;
    kappa           0.41;
    Cb1             0.1355;
    Cb2             0.622;
    Cw2             0.3;
    Cw3             2;
    Cv1             7.1;
    Cs              0.3;
    CDES            0.65;
    ck              0.07;
}

No MRF models present

No finite volume options present
Courant Number mean: 0.2884196 max: 16.346364

Starting time loop

IOdictionary transportProperties (
 IOobject (
  "transportProperties",   //Name of the dictionary
  runTime.constant(),     //Constant over run-time
  mesh ,
  IOobject::MUST_READ,     //Read at creation
  IOobject::NO_WRITE       //Do not write at each time-step
 )
) ;

dimensionedScalar nu (     //Create variable 'nu' and search in transportProperties
 transportProperties.lookup ("nu")
);

volVec torField U (        //volVectorField read in from disk
 IOobject (
  "U" ,
   runTime.timeName() ,    //Associated with runTime database, directory name
   mesh,
   IOobject::MUST_READ
 ),
 mesh
);

Develop New Solver

Steps to create a new solver: Following content is based on numerous documents available online. Some of the references are: OpenFOAM Programming – the basic classes by Prof Gavin Tabor, Introductory OpenFOAM Course University of Genoa, DICCA Dipartimento di Ingegneria Civile, Chimica e Ambientale, The Open Source CFD Toolbox - Programmer’s Guide, Basic OpenFOAM Programming Tutorial: Adding Passive Scalar Transport Equation to icoFoam by Vuko Vukcevic on youTube, OpenFOAM programming tutorial by Tommaso Lucchini ...

Step-01

• solvers: folder for the source code for the distributed solvers.
• test: folder for the source code of several test cases that show the usage of few OpenFOAM classes.
• utilities: folder for the source code for the distributed utilities.
• Allwmake: script file which compiles all the content of solvers and utilities. To compile the test cases in 'test' folder, move to the identified 'test' case directory and compile it by typing 'wmake'.
• foamNewBC and foamNewApp utilities can be used to create a directory with source and compilation files for a new application. e.g. foamNewBC -f -v expInlet will create a folder 'expInlet' in the working directory (directory from which command was issued) and create expInletFvPatchVectorField.C, expInletFvPatchVectorField.H and make directory.

Alternatively, if you want to start from ground zero, use command 'foamNewApp mySolver'. It will create a folder mySolver having subfolder Make. It will also create a file called mySolver.C. Note that it does not create any header file like mySolver.H or createFields.H. The source file mySolver.C will have empty 'Description' field. However, this is the hardest way to develop a solver as almost all controls (e.g. PISO) needs to be typed in. This is apparent from the source code of pimpleFoam that has been added in previous paragraphs.

Step-02

Step-03

Existing 'make' file         |  New 'make' file
XiFoam.C                     |  userXiFoam.C

EXE = $(FOAM_APPBIN)/XiFoam  |  EXE = $(FOAM_USER_APPBIN)/userXiFoam 

$(FOAM_USER_APPBIN)/userXiFoam instructs the compiler to name the application userXiFoam and to copy the executable in the directory $FOAM_USER_APPBIN.

Step-04

Step-05

Info<< "Reading field phi (passive scalar)\n" << endl;
surfaceScalarField phi (
 IOobject (  //phi is an input/ouput object type
   "phi",
   runTime.timeName(),  //The variable is bound to runtime and mesh
   mesh,
   IOobject::NO_READ,
   IOobject::AUTO_WRITE
 ),
 fvc::flux(U)
); 

dimensionedScalar nu (
  "nu",
  dimViscosity,
  transportProperties.lookup("nu")
); 

Step-06

Step-07

Step-08

divSchemes {
  default          none;
  div(phi,U)       bounded Gauss linearUpwind grad(U);
  div(phi,k)       bounded Gauss upwind;
  ...

U {
  solver          smoothSolver;
  smoother        symGaussSeidel;
  tolerance       1e-05;
  relTol          0;
} 

Implementation of a new user-defined B.C.

Step-00:

Step-01:

foamNewBC utility can be used to create a directory with source and compilation files for a new application. e.g. foamNewBC -f -v expInlet will create a folder 'expInlet' in the working directory (directory from which command was issued) along with additional files expInletFvPatchVectorField.C and expInletFvPatchVectorField.H as well as make directory. expInletFvPatchVectorField.C: is the actual source code of the application. This file contains the definition of the classes and functions. expInletFvPatchVectorField.H: header files required to compile the application. This file contains variables, functions and classes declarations.

Step-02:

$(derivedFvPatchFields)/totalTemperature/totalTemperatureFvPatchScalarField.C

LIB = $(FOAM_LIBBIN)/libfiniteVolume

avgControlled/avgControlledFvPatchScalarField.C

LIB = $(FOAM_USER_LIBBIN)/libfiniteVolume

Official release file

  
EXE_INC = \
  -I$(LIB_SRC)/triSurface/lnInclude \
  -I$(LIB_SRC)/meshTools/lnInclude  \

LIB_LIBS =     \
  -lOpenFOAM   \
  -ltriSurface \
  -lmeshTools

Modified file

EXE_INC = \
  -I$(LIB_SRC)/triSurface/lnInclude \
  -I$(LIB_SRC)/meshTools/lnInclude  \
  -I$(LIB_SRC)/finiteVolume/lnInclude

LIB_LIBS =     \
  -lOpenFOAM   \
  -ltriSurface \
  -lmeshTools  \
  -lfiniteVolume

Step-03:

Notice the following sections in the expInletFvPatchVectorField.H file. Different types of private data are declared. These are the variables required for the implementation of the new BC. In implementation of exponential velocity profile, we need to use vectors and scalars, therefore we can keep the lines scalar scalarData_; and vector data_; Other lines can be deleted as those datatypes are not needed.

Also, as we will use two vectors in our implementation, we need to duplicate lines containing "vector data_" (total 6 occurrences). Note that all data members should be of type 'private'. The OpenFOAM convention is to have trailing underscore in the names.

// Private Data
 //- Single valued scalar quantity, e.g. a coefficient
   scalar scalarData_;
 //- Single valued Type quantity, e.g. reference pressure pRefValue_
 // Other options include vector, tensor
   vector data_;

 //- Field of Types, typically defined across patch faces
 // e.g. total pressure p0_. Other options include vectorField
   vectorField fieldData_;

 //- Type specified as a function of time for time-varying BCs
   autoPtr<Function1<vector>> timeVsData_;

 //- Word entry, e.g. pName_ for name of the pressure field on database
   word wordData_;
  
 //- Label, e.g. patch index, current time index
   label labelData_;

 //- Boolean for true/false, e.g. specify if flow rate is volumetric_
   bool boolData_; 

Step-04:

Original section of code

    fixedValueFvPatchVectorField(p, iF),
    meanVelocity_(0.0),
    n_(Zero),
    R_(Zero)

Modified section of code

    fixedValueFvPatchVectorField(p, iF),
    meanVelocity_(0.0),
    n_(1, 0, 0),
    R_(0, 1, 0)

Step-05:

Foam::expInletFvPatchVectorField:: 
expInletFvPatchVectorField (
  const fvPatch& p,
  const DimensionedField<vector, volMesh>& iF
)
:
  fixedValueFvPatchVectorField(p, iF),
  scalarData_(0.0),
  data_(Zero),
  fieldData_(p.size(), Zero), 
  timeVsData_(),
  wordData_("wordDefault"),
  labelData_(-1),
  boolData_(false) {
  }

Foam::scalar Foam::expInletFvPatchVectorField::t() const {
    return db().time().timeOutputValue();
}

Another option is to search for the words fieldData, wordData, labelData and boolData to delete all these lines. In editor 'vim', use :g/fieldData/d to delete all lines containing "profile" (remove /d i.e. use :g/fieldData/ to check or view the lines that will get deleted).

Step-06:

dict.lookup will look for these keywords in the input dictionary

fixedValueFvPatchVectorField(p, iF),
    meanVelocity_(dict.lookup<scalar>("meanVelocity")),
    n_(pTraits<vector>(dict.lookup("n"))),
    R_(pTraits<vector>(dict.lookup("R"))) 

Step-07:

(
   data_
 + fieldData_
 + scalarData_*timeVsData_->value(t())
); 

Step-08:

//Get dimensions of the bounding box: MIN and MAX points
  boundBox bb(patch().patch().localPoints(), true);

//Calculate coordinates of centre point of the patch
  vector ctr = 0.5*(bb.max() + bb.min());
  vector Dh = (bb.max() - bb.min());

//Access centre of the faces in the patch
  const vectorField& c = patch().Cf();

//Computes scalar field to be used for defining the parabolic profile
//& denotes dot product of two vectors in OpenFOAM
  scalarField coord = 2*((c - ctr) & R_)/(Dh & R_); 

Step-09:

n_*meanVelocity_*(1.0 - sqr(coord))

Step-10:

Step-11:

if (mag(n_) < SMALL || mag(R_) < SMALL) {
  FatalErrorIn("expInletFvPatchVectorField(dict)")
  << "n and/or R given with zero  or close to zero."
  << abort(FatalError);   //Abort exeuction - solver will exit
}
n_ /= mag(n_); //This is equivalent to n_ = n_/mag(n_)
R_ /= mag(R_); //This is equivalent to R_ = R_/(mag(y_)

Step-12:

vIn  {
  type         expInlet;  //Name of the boundary condition
  meanVelcoity 2.5;       //User defined values
  n            (1 0 0);
  y            (0 1 0);
} 

Add the name of the library in the controlDict.

Step-13:

Step-14:

These steps are based on user instruction document from WolfDynamics: "Introductory OpenFOAM Course From 16th to 20th July, 2018", University of Genoa, Dipartimento di Ingegneria Civile, Chimica e Ambientale. Different names have been used to avoid following the code verbatim as mentioned in this document. Using different variable names necessitated or enforced user to carefully look into the codes.

Porous Media

The viscous and inertial source terms are defined in the member function addViscousInertialResistance and the power law from is defined in the addPowerLawResistance member function.

Cavitation Models

Reference: Implementing a Zwart-Gerber-Belamri cavitation model, Marcus Jansson, Applied Thermal and Fluid Sciences, Linkoping University

• Solves for two isothermal, incompressible fluids using a Volume of Fluids (VoF) approach where momentum equations are only solved for the mixture.
• Mass transfer through cavitation with the models Merkle, Kunz and Schnerr-Sauer.
• Solves for the liquid fractiona and includes effect of gravity.
• interPhaseChangeFoam uses alphaEqnSubCycle to correct the phase interfaces.
• alphaEqnSubCycle (and alphaEqn) includes vDotAlphal which is part of the phaseChangeTwoPhaseMixtures library.
• vDotAlphal includes mDotAlphal which is different for every cavitation model
• Member functions: rRb, alphaNuc, pCoeff, mDotAlphal, mDotP - as described in the Rayleigh-Plesset and Schnerr-Sauer model equations above.
• The volume source terms vDotvAlphal (vaporization) and vDotcAlphal (condensation) are components of the pair called vDotAlphal. This is found in the phaseChangeTwoPhaseMixtures class, which is a part of the interPhaseChageFoam solver.

• cp -r $FOAM_APP/ ... /phaseChangeTwoPhaseMixtures $WM_PROJECT_USER_DIR/ ... /phaseChangeTwoPhaseMixtures
• cd $WM_PROJECT_USER_DIR/src/phaseChangeTwoPhaseMixtures
• cp -r SchnerrSauer Zwart
• Rename Schnerr-Sauer files: cd Zwart, mv SchnerrSauer.C Zwart.C, mv SchnerrSauer.H Zwart.H

Step-2: Replace each 'SchnerrSauer' with 'Zwart' in Zwart.C and Zwart.H. Use Linux command: sed -i s/"SchnerrSauer"/"Zwart"/g Zwart.*

Step-3: Update wmake 'files' such that it ends with the following two lines:

...
Zwart/Zwart.C
LIB = $(FOAM_USER_LIBBIN)/libphaseChangeTwoPhaseMixtures

Step-4: Zwart.H - Remove the private member data for n_, dNuc_, rRb_ and alphaNuc in Zwart.H Add the new member data between Cv_ and p0_.

//- Nucleation site volume
  dimensionedScalar rNuc_;
//- Nucleation site radius
  dimensionedScalar Rb_; 

Step-5: Zwart.C - Remove the constructors for n_ and dNuc_. Add in the constructor, between Cv_ and p0_: rNuc_("rNuc", dimless, phaseChangeTwoPhaseMixtureCoeffs_), Rb_("Rb", dimLength, phaseChangeTwoPhaseMixtureCoeffs_),

Step-6: Zwart.C - Remove the entire member functions rRb and alphaNuc. Replace the member functions pCoeff and mDotAlpha

pCoeff: Zwart Model

Foam::tmp<Foam::volScalarField>
Foam::phaseChangeTwoPhaseMixtures::Zwart::pCoeff (
  const volScalarField& p
) 
const {
  return (3*rho2())*sqrt(2/(3*rho1())) /(Rb_*sqrt(mag(p - pSat()) + 0.01*pSat()));
} 

mDotAlphal (Schnerr-Sauer Model)

Foam::Pair<Foam::tmp<Foam::volScalarField>>
Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotAlphal() const  {
  const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
  volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
  volScalarField pCoeff(this->pCoeff(p));
  return Pair<tmp<volScalarField>>
  (
   Cc_*limitedAlpha1*pCoeff*max(p - pSat(), p0_),
   Cv_*(1.0 + alphaNuc() - limitedAlpha1)*pCoeff*min(p - pSat(), p0_)
  );
} 

mDotAlphal (Zwart Model)

Foam::Pair<Foam::tmp<Foam::volScalarField>>
Foam::phaseChangeTwoPhaseMixtures::Zwart::mDotAlphal() const {
  const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
  volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
  volScalarField pCoeff(this->pCoeff(p));
  return Pair<tmp<volScalarField>>
  (
    Cc_*pCoeff*max(p - pSat(), p0_),
    Cv_*rNuc_*pCoeff*min(p - pSat(), p0_)
  );
} 

mDotP (Zwart Model)

Foam::Pair<Foam::tmp<Foam::volScalarField>>
Foam::phaseChangeTwoPhaseMixtures::Zwart::mDotP() const {
  const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
  volScalarField pCoeff(this->pCoeff(p));
  volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
  return Pair<tmp<volScalarField>>
  (
    Cc_*(1.0 - limitedAlpha1)*pos(p - pSat())*pCoeff,
    (-Cv_)*rNuc_*limitedAlpha1*neg(p - pSat())*pCoeff
  );
} 

Step-7: remove n_ and dNuc_ from the read function and add rNuc_ and Rb_.

bool Foam::phaseChangeTwoPhaseMixtures::Zwart::read() {
  if (phaseChangeTwoPhaseMixture::read())  {
    phaseChangeTwoPhaseMixtureCoeffs_ = optionalSubDict(type() + "Coeffs");
    phaseChangeTwoPhaseMixtureCoeffs_.lookup("Cc") >> Cc_;
    phaseChangeTwoPhaseMixtureCoeffs_.lookup("Cv") >> Cv_;
    phaseChangeTwoPhaseMixtureCoeffs_.lookup("rNuc") >> rNuc_;
    phaseChangeTwoPhaseMixtureCoeffs_.lookup("Rb") >> Rb_;
    return true;
  } else {
    return false;
  }
} 

Step-8: Compile the class again using wmake. Go to the run directory and copy the throttleSchnerr tutorial. Remove the old time directories. In constant/ transportProperties, change phaseChangeTwoPhaseMixture from SchnerrSauer to Zwart to use the new cavitation model. Add the new model coeffcients after the section with SchnerrSauerCoeffs. Note the addition of rNuc and rRb and the values for Cc and Cv.

ZwartCoeffs {
  Cc 	0.01;
  Cv 	50;
  rNuc 	5.0e-04;
  Rb 	1.0e-06;
} 

Step-9: Solve the case, now with Zwart-Gerber-Belamri cavitation model by typing 'interPhaseChangeFoam' on command prompt.

Kunz mass transfer model implementation can be found here.

Implement Non-Newtonian Viscosity Model

Step-01:

Step-02:

Step-03:

Info << " Reading field nu" << nl << endl;

volScalarField nu (
  IOobject (
   "nu", runTime.timeName (),
   mesh,
   IOobject::MUST_READ, IOobject::AUTO_WRITE
  ), mesh
);

Step-04:

volScalarField J2 = 0.5* magSqr(symm (fvc::grad (U ))) + dimensionedScalar ("tiny", dimensionSet (0,0,-2,0,0,0,0), 0.0001);

nu = sqr(pow(sqr(eta)*J2, 0.25) + pow(tau_y/2, 0.5)) / (rho*sqrt(J2));

Step-05:

Step-06:

Step-07:

transportModel Casson;
CassonCoeffs {
  eta   eta   [1 -1 -1 0 0 0 0 ]  5.00;
  tau_y tau_y [1 -1 -2 0 0 0 0 ]  12.0;
  rho   rho   [1 -3  0 0 0 0 0]   1100;
}

Step-08:

Dynamic Mesh

A good report on Implementation of dynamicMesh in OpenFOAM can be found here. The implementation of point-wise deformation of meshes can be found in this tutorial document.

Example implementation in rhoPimpleFoam and explanations are described below.

#include "fvCFD.H"
#include "dynamicFvMesh.H"  //Inclusion of classes to handle the dynamic mesh
#include "psiThermo.H"
#include "turbulenceModel.H"
#include "bound.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"

//Creates mesh ojbect as a dynamicFvMesh and not as fvMesh
#include "createDynamicFvMesh.H" 
#include "initContinuityErrs.H"

pimpleControl pimple(mesh);

#include "readControls.H"
#include "createFields.H"
#include "createFvOptions.H"
#include "createPcorrTypes.H"
#include "createRhoUf.H"
 Includes:
  surfaceVectorField rhoUf (
   IOobject (
    "rhoUf",
    runTime.timeName(),
    mesh,
    IOobject::READ_IF_PRESENT,
    IOobject::AUTO_WRITE
   ),
  fvc::interpolate(rho*U)
  );

#include "CourantNo.H"
#include "setInitialDeltaT.H"

//-----------------------------------------------------------------------------
while (runTime.run()) {
 {
  //Absolute convective flux due to the flow motion
  volScalarField divrhoU (
   "divrhoU",
    fvc::div(fvc::absolute(phi, rho, U))
  );
  runTime++;
  volVectorField rhoU("rhoU", rho*U);
  
  // Do any mesh changes: Triggers the mesh motion updating the new point
  // positions and calculating the mesh flux
  mesh.update();
  
  if (mesh.changing() && correctPhi) {
   phi = mesh.Sf() & rhoUf;
   #include "correctPhi.H"
   // Make the fluxes relative to the mesh-motion
   fvc::makeRelative(phi, rho, U);
  }
 }
} 

Electromagnetics

Variables in electromagnetism:

• Q: Electric charge [C]
• A: magnetic potential vector
• H: magnetic field strength (A/m or N/Wb or Oersted)
• J: current density [A/m²]
• D: electric flux density or electric displacement [C/m²]
• E: electric field strength or electric field intensity [V/m]
• B: magnetic flux density [T] or [Wb/m²]
• ε: absolute electrical permittivity, ε₀ = absolute permittivity of vacuum = 8.854E-12 [F/m]
• ε_r: relative electrical permittivity or Dielectric Constant (defined as 'k' in electrostaticFoam)
• μ: absolute magentic permeability (subscript '0' indicates permeability of vacuum or air = 4π×10^-7 H/m)
• ρ: volume density of charge (∇D = ρ) or electrical resistivity (electrical resistance, R = ρ × Length / Area) - [Ω-m] or [Ohm-m]
• σ: electrical conductivity = 1/ρ - [1/Ω-m]

Magnetic Potential Equation

B = fvc::curl(A); 
Je = -sigma*(fvc::grad(ElPot));
solve ( 
  fvm::laplacian(A) == muMag * sigma * (fvc::grad(ElPot))
); 

Electric Potential Equation

solve (
  fvm::laplacian(sigma, ElPot)
); 

1. magneticFoam: solver for the magnetic field generated by permanent magnets
2. electrostaticFoam: solver for electrostatics - the static current (current does not change with time). Note that current is flow of charge. The equations which are solved are:

   

   solve (
      fvm::laplacian(phi) + rho/epsilon0
   ); 

   

   rhoFlux = -k*mesh.magSf()*fvc::snGrad(phi);
   solve (
     fvm::ddt(rho) + fvm::div(rhoFlux, rho)
   );

3. mhdFoam: Solver for magnetohydrodynamics (MHD) - incompressible, laminar flow of a conducting fluid under the influence of a magnetic field. An applied magnetic field H acts as a driving force, at present boundary conditions cannot be set via the electric field E or current density J. The fluid viscosity nu, conductivity sigma and permeability mu are read in as uniform constants.

   A fictitious magnetic flux pressure pB is introduced in order to compensate for discretisation errors and create a magnetic face flux field which is divergence free as required by Maxwell's equations.

   However, in this formulation discretisation error prevents the normal stresses in UB from cancelling with those from BU, but it is unknown whether this is a serious error. A correction could be introduced whereby the normal stresses in the discretised BU term are replaced by those from the UB term, but this would violate the boundedness constraint presently observed in the present numerics which guarantees div(U) and div(H) are zero.

ELECTROMAGNETIC FOAM

Step-01: Create a folder say 'eMagFoam' in your $WM_PROJECT_USER_DIR directory: cd $WM_PROJECT_USER_DIR ❯ mkdir eMagFoam ❯ cd eMagFoam. Note that foamNewSource shall create app (solver) name based on the folder

Step-02: Generate a new solver template files using command - foamNewSource App emagFoam. This will create emagFoam.C file and Make folder. The Make/files will contain following two lines:

eMagFoam.C

EXE = $(FOAM_APPBIN)/eMagFam

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude

EXE_LIBS = \
    -lfiniteVolume

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude

EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools

Step-03: Update the eMagFoam.C file by adding following two include statments. Note that the createMesh.H should be included before 'createField.H' else compilation step (wmake) will throw error.

#include "createMesh.H"

#include "createFields.H"

Add the solve statements to solve PDE and described earlier.

solve (
  fvm::laplacian(sigma, ElPot)
);

solve (
  fvm::laplacian(A) == muMag * sigma * (fvc::grad(ElPot))
);

Step-04: Add statements to write the results into a new 'time' directory.

Step-05: construct and initialize muMag, sigma, Elpot, A, B in createFields.H which will contain the statements and expression as per this file.

Step-06: Compile the new solver using 'wmake'. If the process is successful, you will get following statement as last line of the output: -lm -o /home/X/ OpenFOAM/X-8/ platforms/linux64GccDPInt32Opt/ bin/eMagFoam where 'X' stands for the user name.

Step-07: Create mesh, A, ElPot and sigma dictionaries in 0/ folder. Update the controlDict with solver eMagFoam. Run the case. Note the units of A (Magnetic Potential) is [1 1 -2 0 0 -1 0], ElPot (Electric Potential, φ) has unit of [1 2 -3 0 0 -1 0] and unit of sigma (electrical conductivity) is [-1 -3 3 0 2 0 0].

Step-08: Check the result in ParaView.

Use postProcess functionality to extract components such as (a) postProcess -func 'components(A)' -time 1 and (b)postProcess -func 'components(B)' -time 1.

The case file can be downloaded from here. Use the utilities blockMesh, changeDictionary, setFields and then solver eMagFoam in sequence.

Acoustics: Develop acoustic / noise models.

acoustic = density0 * fvc::div(fvc::div(U*U));  //quadrupole distribution 
monoPole = -fvc::div(psurfaceSource);
acousticSources =  monoPole + acoustic;
fvScalarMatrixpaEqn (
  1.0/c/c * fvm::d2dt2(pa) - fvm::laplacian(pa); pa = p - pRef
); 
solve(paEqn == acoustic);  // sound pressure as acoustic sources 

Disclaimer:

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