Project 1- Creating your own openFOAM solver

Solving scalar transport equation in Openfoam:

A scalar transport equation (given below) is to be solved in a computational domain along with the momentum and continuity equation for unsteady incompressible flow. Icofoam solver is customized to include the solution for scalar transport equation for this purpose.

`(dels)/(delt) + grad.(vecUs) = 0`

s -> scalar to be transported in the domain

`vecU` -> velocity vector in the domain.

Time derivative of 's' is dependent on convection of 's' using the momentum in the domain. Diffusion and source terms are not present in the scalar transport equation.

To prepare the solver:

Icofoam folder inside $FOAM_SOLVERS/incompressible is copied to run directory. The file, icofoam.C is renamed to scalarFoam.cpp. Make/files(shown below) is modified to compile this renamed file and to place the executable in $FOAM_USER_APPBIN with appropriate name

scalarFoam.cpp

EXE = $(FOAM_USER_APPBIN)/scalarFoam

Then, createFields.H file updated to include the below snippet to create a volume scalar field 's' which is stored as cell center values in the mesh.

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

scalarFoam.cpp is updated to create the scalar transport equation that includes the time derivative and convective term. This equation is solved once every time-step outside the piso loop so that 's' is updated based on the final values of `vecU` and pressure, 'p'. The following code snippet is added after the piso correct while loop ends and before runTime.write() is called, so that, the calculated values of 's' at cell centers are written out at prescried time-steps.

        fvScalarMatrix SEqn
        (
            fvm::ddt(s)
          + fvm::div(phi, s)
        );

        SEqn.relax();
        SEqn.solve();

After the files are saved, wmake command is called inside the scalarFoam folder which compiles the updated files and creates an executable in $FOAM_USER_APPBIN.

To prepare the case:

cavity tutorial is copied from $FOAM_TUTORIALS/incompressible/icoFoam folder to $FOAM_RUN folder. Boundary conditions and solver settings have to be updated to include the scalar that is to be solved. First the boundary conditions for 's' is set by creating a new file named 's' inside the '0' folder as shown below. 's' is initialized to 0 in the internal domain which will be modified using setFieldDict. Boundary conditions for s is set as zeroGradient for all the 2D boundaries.

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2112                                  |
|   \\  /    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
{
    movingWall
    {
        type            zeroGradient;
    }
    fixedWalls
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
}

In the file, system/controlDict, application is set to scalarFoam and endTime is set to 1.

In the file, system/fvSchemes, grad(s) and div(phi,s) are set inside gradSchemes and divSchemes respectively with Gauss linear.

In the file, system/fvSolution following code snippet is added to specify the solution procedure for SEqn created in scalarFoam solver. "PCG" solver used for 'p' is not suitable for asymmetric matrix formed by SEqn. Hence the following options are used.

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

In the current problem, we assume the scalar is sedimented in the lower half of the cavity domain. So 's' is set to 1 in the lower half of domain and 0 in the upper half. When the top wall moves in positve x direction, a vortex will be created in the domain. Due to this, the scalar is expected to mix in the domian following the vectors of velocity in the domain. To set 's' as prescribed above, setFieldsDict is created in system directory with the contents as shown below.

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

defaultFieldValues
(
    volScalarFieldValue s 0
);

regions
(
    // Set cell values
    // (does zerogradient on boundaries)
    boxToCell
    {
        box (0 0 0) (0.1 0.05 0.01);

        fieldValues
        (
            volScalarFieldValue s 1
        );
    }
);

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

Solution procudure:

blockMesh command is called inside the cavity folder to create the mesh as per the settings in system/blockMeshDict. polyMesh is created in constant folder which will be used by the solver.

setFields command is called, which will modify the 's' file in 0 folder as per the setup in the file, setFieldsDict.

scalarFoam command is called, which will solve for `vecU`, 'p' and 's' in the domain from 0 to 1 seconds with `Deltat`=0.005 writing out the fields every 0.1 seconds as set in controlDict file. Solution is seemed to converge well in each time-step by judging the residuals at the time-steps.

The solution is visualized in paraview using the command paraFoam. Node values of 's' is shown as contour on domain and an animation is exported from paraview as shown below.