Challenge

Project 1- Creating your own openFOAM solver

KKavendra Thakare16 May 2022

Aim: To implement a customised solver in OpenFOAM.

Project Description:
Use the icoFoam solver in the applications/solvers/incompressible/icoFoam
directory to create a customised solver - scalarFoam

The scalar transport equation is defined as:

`frac{dels}{delt} + nabla*(Us) = 0`

where `s` is the scalar quantity being transported by the velocity `U`.

Steps to make customised solver and its compilation:

To copy the respective solver source codes from the $FOAM_SOLVERS.
The incompressible icoFoam solver source code is copied from the
$FOAM_SOLVERS/incompressible/icoFoam to the my_solvers directory.
The copied solver source code directory is renamed to scalarFoam from icoFoam .
The command mv -r icoFoam scalarFoam has been used to rename the directory.
The solver code directory includes three main components as follows
1. Make directory which further contains files and options files.
  The file files contains the source code file and location of the executable where its get stored.
  The following changes are made into the files.
  
  The file options contains the list of the included executables and linked shared libraries to be compiled
  with the source code. The option file is kept as it is.
2. creatFields.H , The header file which includes the IOdictionaries of the parameters used in the code.
  The IOobject for the scalar `s` has been added to the existing creatFields.H .
  The edited code is as follows:
3. scalarFoam.C ,
  The source code file which has been renamed to the scalarFoam.C as required by using appropriate
  renaming command. Since the name of the source code file should be same as the name of the solver.
  The following modifications are made into the source code to include the scalar `s` transport solution
  as mentioned earlier.The required scalar transport equation for scalar `s` has been added
  towards end of the file after correction of the velocity filed.
After making all the necessary modifications into the solver source code directory.
The solver has been compiled successfully using the commands,
wclean - to clear the previously build files and
wmake - to compile the solver.
The compilation results are as shown in the figure.

The solver executable scalarFoam* is saved into the $FOAM_USER_APPBIN.
The following screenshot shows the same.

The scalarFoam solver directory is made up of following files and their respective codes are as follows:

Make/files:

scalarFoam.C

EXE = $(FOAM_USER_APPBIN)/scalarFoam

Make/options:

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

EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools

creatFields.H

Info<< "Reading transportProperties\n" << endl;

IOdictionary transportProperties
(
    IOobject
    (
        "transportProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

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

Info<< "Reading field p\n" << endl;
volScalarField p
(
    IOobject
    (
        "p",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Reading field U\n" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Reading field s\n" << endl;
volScalarField s
(
    IOobject
    (
        "s",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);


#include "createPhi.H"

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

scalarFoam.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    scalarFoam

Group
    grpIncompressibleSolvers

Description
    Transient solver for scalar transport

    \heading Solver details
    The solver uses the PISO algorithm to solve the continuity equation:

        \f[
            \div \vec{U} = 0
        \f]

    and momentum equation:

        \f[
            \ddt{\vec{U}}
          + \div \left( \vec{U} \vec{U} \right)
          - \div \left(\nu \grad \vec{U} \right)
          = - \grad p
        \f]

    Where:
    \vartable
        \vec{U} | Velocity
        p       | Pressure
    \endvartable

    \heading Required fields
    \plaintable
        U       | Velocity [m/s]
        p       | Kinematic pressure, p/rho [m2/s2]
    \endplaintable

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "pisoControl.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Transient solver for incompressible, laminar flow"
        " of Newtonian fluids."
    );

    #include "postProcess.H"

    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMesh.H"

    pisoControl piso(mesh);

    #include "createFields.H"
    #include "initContinuityErrs.H"

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;

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

        #include "CourantNo.H"

        // Momentum predictor

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

        if (piso.momentumPredictor())
        {
            solve(UEqn == -fvc::grad(p));
        }

        // --- PISO loop
       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);

            // Non-orthogonal pressure corrector loop
           while (piso.correctNonOrthogonal())
           {
                // Pressure corrector

              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();
        }
       
       fvScalarMatrix sEqn
        (
            fvm::ddt(s)
          + fvm::div(phi,s)
        );

         sEqn.solve();
         s.correctBoundaryConditions();


        runTime.write();

        runTime.printExecutionTime(Info);
    }

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //

Steps to make tutorial case for customised solver and get the simulation results:

To copy the respective tutorial case from the $FOAM_TUTORIALS.
The cavity case from incompressible icoFoam solver tutorial cases, is copied from the
$FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity to the run directory.
The copied tutorial case directory is renamed to scalarcavity from cavity .
The command mv -r cavity scalarcavity has been used to rename the directory.
The tutorial consists of three main components as explained below
- 0 directory , It contains initial and boundary conditions for all variables used for the simulation.
  Here 0 directory contains files for velocity `U` , dynamic pressure `p` and scalar `s` .
  The different time directories as written alongside 0 directory into the case directory
  as per the write interval settings into the controlDict.
  - file `U`:
  - file `p`:
  - file `s`: This file is modified from the copied version of the file `p`.
- constant directory:
  It contains the constant properties dictionaries and the polyMesh directory which is
  generated after the Mesh has been generated (here after blockMesh command).
  Here the dictionary for transportProperties is included.
- system directory:
  It contains the blockMeshDict, controlDict, fvSchemes, fvSolutions, and setFieldsDict.
  Except setFieldsDict file all other files have been already part of the existing tutorial case system directory.
  The setFieldsDict file has been copied from the damBreak tutorial case from the $FOAM_TUTORIALS/multiphase/interFoam/laminar/damBreak/damBreak/system/setFieldsDict
  - The blockMeshDict file has not been modified.
  - The controlDict file has been modified for the solver to be used
    (scalarFoam instead of icoFoam) as follows,
  - The file fvSchemes enlists all the schemes used for various types of numerical calculations.
    The same is modified according to the scalar transport equation we added.
    Specifically the Gauss upwind - divSchemes for div(phi ,s) has been added to the existing file.
  - The file fvSolution contains the solvers for all of the variables used in the code.
    Here the solver for the scalar `s` has been added as shown in the figure.
    Specifically, the PBiCG solver with DILU preconditioner has been used for the solver of scalar `s`.
  - The setFieldsDict file is used to set the internal fields during runtime.
    The same is used to set field of scalar `s`.
    The circular scalar field has been set using sphereToCell command with specified centre and radius.
    The file contents are as follows,
After all of the above modifications into the tutorial case scalarcavity for the use of solver scalarFoam.
The tutorial case has been simulated with the following sequence of the command.
1. blockMesh - to generate mesh. The directory polyMesh has been generated into the constant subdirectory.
  The results are as follows.
  
  polyMesh directory is generated after blockMesh command execution,
2. checkMesh - to check the mesh quality. The results are as follows.
3. setFields - to set the internal scalar field. The results are as follows.
4. scalarFoam - to solve the case using the scalarFoam solver, which has been customised as explained above.
  The final time step results are as follows.
5. paraFoam - The alias command to generate simulation.foam file and open it into the Windows native ParaView.
  Where the path to open the Windows native Paraview is supplied into the alias command along
  with making of the simulation.foam file.
  The results from the paraview for initial and final time steps with intermediate step are as follows.
  - Initial scalar field as set by setFields command.
  - Intermediate scalarfield at time 0.5 sec as solution progresses.
  - Final scalar field at the end of 1 sec.

The code of all the files from the scalarcavity case are as follows.

0/U file

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    movingWall
    {
        type            fixedValue;
        value           uniform (1 0 0);
    }

    fixedWalls
    {
        type            noSlip;
    }

    frontAndBack
    {
        type            empty;
    }
}

// ************************************************************************* //

0/p file

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    movingWall
    {
        type            zeroGradient;
    }

    fixedWalls
    {
        type            zeroGradient;
    }

    frontAndBack
    {
        type            empty;
    }
}

// ************************************************************************* //

0/s file

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      s;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 0 0 0 0];

internalField   uniform 300;

boundaryField
{
    movingWall
    {
        type            zeroGradient;
    }

    fixedWalls
    {
        type            zeroGradient;
    }

    frontAndBack
    {
        type            empty;
    }
}

// ************************************************************************* //

constant/transportProperties

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

nu              0.01;


// ************************************************************************* //

system/blockMeshDict

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

scale   0.1;

vertices
(
    (0 0 0)
    (1 0 0)
    (1 1 0)
    (0 1 0)
    (0 0 0.1)
    (1 0 0.1)
    (1 1 0.1)
    (0 1 0.1)
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (20 20 1) simpleGrading (1 1 1)
);

edges
(
);

boundary
(
    movingWall
    {
        type wall;
        faces
        (
            (3 7 6 2)
        );
    }
    fixedWalls
    {
        type wall;
        faces
        (
            (0 4 7 3)
            (2 6 5 1)
            (1 5 4 0)
        );
    }
    frontAndBack
    {
        type empty;
        faces
        (
            (0 3 2 1)
            (4 5 6 7)
        );
    }
);

mergePatchPairs
(
);

// ************************************************************************* //

system/controlDict

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     scalarFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         1;

deltaT          0.005;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


// ************************************************************************* //

system/fvSchemes

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
   // grad(s)         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
    div(phi,s)      Gauss upwind; //linear;
}

laplacianSchemes
{
    default         Gauss linear orthogonal;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         orthogonal;
}


// ************************************************************************* //

system/fvSolution

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

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

    s
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-06;
        relTol          0;
    }

}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}


// ************************************************************************* //

system/setFieldsDict

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      setFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

defaultFieldValues
(
    volScalarFieldValue s 300
);

regions
(
    sphereToCell
    {
        centre (0.05 0.05 0); 
        radius 0.03;
        fieldValues
        (
            volScalarFieldValue s 500
        );
    }
);


// ************************************************************************* //

As the desired results are obtained from the simulation of the test case scalarcavity using the customised solver scalarFoam.
Thus the solver scalarFoam has been successfully customised for scalar transport and tested successfully on the test case scalarcavity.

Project 1- Creating your own openFOAM solver

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.