Challenge

Project 1- Creating your own OpenFOAM solver

PPavan Vittanala3 Nov 2023

AIM: Use the icoFoam solver to create your own solver - scalarFoam.

Description:

Solver Name: icoFoam

Directory: applications/solvers/incompressible/icoFoam

The icoFoam is transient solver for incompressible and laminar flows. It solves the Navier-Stokes equation using the PISO(Pressure Iimplicit with Spliting of Operators) algorithm.

The scalar transport equation is defined as:

`(dels)/(delt) + nabla.(us) = 0`

where s is the scalar quantity being transported by the velocity u.

Procedure:

First, I copied the icoFoam solver from applications/solvers/incompressible/icoFoam and renamed it as scalarFoam.
Here, the copied the icofoam tutorial from the $FOAM_TUTORIALS/incompressible/icoFoam/cavity as a test case, and moved it into the scalarFoam solver file.

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
    icoFoam

Group
    grpIncompressibleSolvers

Description
    Transient solver for incompressible, laminar flow of Newtonian fluids.

    \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
    (
        "Scalar solver for scalar quantity being transported by velocity U"
        
    );

    #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

        fvScalarMatrix sEqn
        (
            fvm::ddt(s)
          + fvm::div(phi, s)
          //- fvm::laplacian(nu, U)
        );
        
         solve(sEqn);
         
         

        runTime.write();

        runTime.printExecutionTime(Info);
    }

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

    return 0;
}


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

Here, the scalar quantity 's' is defined as sEqn and stored it in fvScalarMatrix. In the sEqn the LHS side the first term is the time derivative of s i.e.., `(dels)/(delt)` and the second term is a convection term i.e.., divergence of velocity cector field (`nabla.(us)`).

createFields.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
);

dimensionedScalar su
(
    "su",
    dimViscosity,
    transportProperties
);*/

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 sRefCell = 0;
scalar sRefValue = 0.0;
setRefCell(s,  sRefCell, sRefValue);
mesh.setFluxRequired(s.name());
*/

Here, created IOobject fields of U and s in this file.

Make/files
```
scalarFoam.C

EXE = scalarFoam
```

Make/options

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

EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools

Then run the solver by using the command 'wmake' and results as follows:

Here, no error is found, and our solver is working fine. Then, Now setup the test case files i.e.., scalarCavity files as follows

Now inside the sclarCavity there are '0', 'constant', 'system' folders.

system/ blockMeshDict

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

scale   1;

vertices
(
    (0 0 0)
    (5 0 0)
    (5 1 0)
    (0 1 0)
    (0 0 1)
    (5 0 1)
    (5 1 1)
    (0 1 1)
);

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

edges
(
);

boundary
(
    inlet
    {
        type patch;
        faces
        (
            (0 4 7 3)
        );
    }
    outlet
    {
        type patch;
        faces
        (
        
            (2 6 5 1)
            
        );
    }
    frontAndBack
    {
        type empty;
        faces
        (
          (3 7 6 2) 
          (1 5 4 0)
          
        );
    }
    topAndbottom
    {
    	type wall;
        faces
        (
        
   	(0 3 2 1)
        (4 5 6 7)
        
    	);
    }
);

mergePatchPairs
(
);


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

Here, the test case is a rectangle of dimensions 5x1x1 meter. And the boundries are named as inlet, outlet, frontAndback, topAndbottom. And defined the boundary type of each boundary as shown above.

system/ controlDict

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

application     scalarFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         2;

deltaT          0.005;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Here, the solver is renamed as scalarFoam and also defined the start time and endtime and deltaT to our case setup as shown above.

fvSchemes

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

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    
}

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

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default 		linear;
}

snGradSchemes
{
    default         orthogonal;
}


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

fvSolutions

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

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


}




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

setFieldsDict

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

defaultFieldValues
(
    volScalarFieldValue s 0
    
);

regions
(
    boxToCell
    {
        box (1 0.2 0.2) (2 0.8 0.8);
        		
        fieldValues
        (
     
            volScalarFieldValue s 1
        );
    
    }
);


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

Here, a setFieldsDict is setup as shown above. Here two scalarfields are defined at s=0 and s=1 fieldvalues. One phase is at setfield of 0 and another phase is made at setfield of 1.

0/s file

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

dimensions      [0 0 0 0 0 0 0];

internalField   uniform 0;



boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
    topAndbottom
    {
        type            zeroGradient;
    }
}


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

0/ U file

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2306                                 |
|   \\  /    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 (2 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           $internalField;
    }

    outlet
    {
        type            zeroGradient;
        
    }

    frontAndBack
    {
        type            empty;
    }
    topAndbottom
    {
        type            noSlip;
    }
}


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

And here the velocity at inlet is 2m/sec.

Post-Processing:

Run the solver on the test case by using the following commands:

1. blockMesh: Used to create the mesh as the specifications provided inside the blockMeshDict file

2. setFields: Used to setup the phases of the case file as per instructions given inside the setFieldsDict

3. scalarFoam: Finally run the scalarFoam solver to start the stimulation.

4. paraFoam: Used to visualize the results.

Results:

Here, the courant number is resulted as 0.5

And the residual of 's' as intial residual = 0.00307136 and final residual = 2.85487e-17

At timestep = 0

At the final timestep = 2

Here the animated results are attached in the given link: https://drive.google.com/file/d/13bwnoTbHvILuHNaVAKnk0sMyDUubU1Xl/view?usp=drive_link

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.