Week 8 - Simulation of a backward facing step in OpenFOAM

AIM : To simulate an incompressible-laminar-viscous flow through the backward facing step geometry.

 We should perform a transient simulation. The solver can be chosen based on the described physics of the flow. 

Here based on the decribed physics of the flow we have chosen icoFoam solver.

Case 1 - Simulate the flow without using any grading factor (i.e., GF = 1)

Case 2 - Simulate the flow with grading factor of 0.2. The cells should be finer near the walls (includig the step wall).

You can use the following mesh parameters but you can feel free to modify these values.

Mesh specification

1. Number of cells along the x direction (longer dimension) = 200
2. Number of cells along the longer y direction = 20;

The typical approach adopted to performing numerical simulation in OpenFOAM involves the following steps:

1. Identifying the tutorial in the 'FOAM_TUTORIALS' folder that contains the appropriate governing equations relevant to the physics of the flow for the given problem. A typical file structure of an OpenFOAM tutorial consists of the below mentioned directories and files:
   1. A 'system' directory with the 'controlDict', 'fvSchemes' and 'fvSolution' files used for specifying the numerical schemes and solution procedure along with control on the simulation run properties. Additionally, for the given numerical problem, a file by the name 'blockMeshDict' mentions the creation of different blocks representing the geometry as well as its mesh properties.
   2. A 'constant' directory that contains a 'polyMesh' directory outlining the various aspects of the mesh such as 'boundary', 'points', 'faces', etc., and a file that defines the physical properties, in this case, 'transportProperties' specifying the value of kinematic viscosity.
   3. A time directory, usually called the '0' directory that consists of files representing the different fields or fluid flow properties with values as required at the start of simulation run or at time $t=0$.
2. Copying the now identified tutorial folder from the 'FOAM_TUTORIALS' folder to the 'FOAM_RUN' folder, specifically created to allow simulation run of OpenFOAM numerical simulation cases, with either the same name or preferably a name relevant to the geometry to be discretized.
3. Replacing the geometry and boundary list of the pre-existing tutorial with the geometry and boundary list pertaining to the given problem and defining mesh parameters with regards to the flow physics, in this case by updating the 'blockMeshDict' file (Section 2.1 Modelling).
4. Updating the fields' files present inside the '0' directory viz., 'U' and 'p', with the boundary list as specified in the 'blockMeshDict' file and appropriate values at the relevant boundary to ensure the given numerical problem has the right boundary conditions and that the solution would be correctly initialised (Section 2.2 Boundary Conditions).
5. Updating the 'controlDict' file, in this case the 'deltaT' and 'endTime' parameters with values that take into consideration domain specific attributes such as length of the domain in the direction of flow, fluid flow velocity and avoiding violation of Courant-Friedrichs-Lewy Condition or the CFL number by limiting its mean and maximum value during simulation run below 1 (Section 2.4 Determination of transient parameters).

STEPS INVOLVED IN SIMULATION:

1. In order to start the simulation process, First start terminal in Linux Operating system using ctrl+alt+T buttons on key board or right click -> Open Terminal. After starting terminal yo will see following window.
2. Now, we have to simulate our problem in openFOAM and we are going to solve our problem using tutorials present in openFOAM, therefore, just type ‘tut’ to open OpenFOAM tutorial directory. And type ‘ls’ to identify what kind of tutorials present inside this directory.
3. We have to solve incopressible flow simulation for this problem so select incompressible directory by typing ‘cd incompressible/’ and then
again type ‘ls’ to explore folders in incompressible folders.
4. we will select icoFoam folder by typing ‘cd icoFoam/’ and by typing ‘ls’ we explore this folder in terminal.
icoFoam is the solver in openFoam which is use as transient solver for incopressible, laminar flow of Newtonian fluids. (Newtonian fluid is a fluid in which the viscous stress arising from its flow, at every point, are linearly correlated to the local strain rate - the rate of change of its
deformation over time. Newtonian fluids are the simplest mathematical models of fluids that account for viscosity.)
5. Here, we have, again to option, cavity and elbow. They both are part of the icoFoam solver and cavity section use to create required domain and black mesh so select ‘cavity’ folder.
6. In cavity, we are, again, going to see cavity folder, where we have to define our fluid domain and generate block mesh. But we are not going to make any changes here in tutorial section. First of all, OpenFoam does no let anyone make any changes in their official tutorial and secondly, any changes with their official tutorials my corrupt openFoam. An in a bid to correct it, we may need to reinstall open foam all over again. So, we are going to copy this cavity folder using ‘cp -r cavity’ command and paste it in openFoam run folder by typing ‘cp -r cavity
$FOAM_RUN/cavity_test’.
7. And then open run folder by typing ‘cd $FOAM_RUN’. where we can see cavity_test file.
8. Now, we can make any changes to the cavity_test file, which will never affect the openFoam official files. Open cavity_test folder. Here, we can provide specification for our block mesh in system folder.
9. There are 4 script file system folder and by opening blockMeshDict we can create domain and provide required specification for block mesh. So, open this file by typing ‘gedit blockMeshDict’.
10. Now, make following changes in blockMeshDict script file. The blockMeshDict file is a script file that contains information regarding all
aspects of the geometry.

Methodology:
Step 1: Selection of suitable solver
OpenFOAM has a wide variety of solvers for specific use-cases. Based on the requirement, one has to be selected. In this report,
the physics being solved is the flow of incompressible viscous laminar flow. The suitable solver for such a case is icoFOAM.
Step 2: Setting up the case
OpenFOAM, being an open-source software lacks a graphical user interface. Therefore, the case has to be set up using the folder
structure inherent to it. As a beginner, the best approach is to copy a tutorial case that uses the required solver and change the
settings as required.
The cavity tutorial, which uses the icoFOAM solver was copied for this report.
Step 3: Creating a mesh
Once the folder structure is set up, the mesh for the geometry has to be generated. For this purpose, blockMesh, the inherent
feature of OpenFOAM was used. The blockMesh requires the definition of rectangular blocks within the geometry and hence the
geometry was split into 3 blocks as shown below. First, the vertices should be defined and then the blocks should be created from
the vertices such that the local coordinated system of the block is the same.

OpenFOAM is a 3D solver and therefore blocks should also be 3D. Hence the blocks were created based on the configuration
shown above. A 2D case can be simulated by specifying a single cell in the z-direction. In addition to the blocks, the number of
cells and the grading should be defined for each block. 

Two variations of the same case were simulated. In the first case, no grading was used in the mesh. In the second one, a grading of
0.2 was used and the results were compared.
Grading is a term used to defined the relative size of the last cell in a block with respect to the first cell.
The above parameters should be defined in the file, blockMeshDict inside the settings folder. The blockMeshDict file used in this
simulation is shown below:

Boundary condition specification

The domain specifications are provided in the following figure.

To create the geometry as shown in the figure, we need edit the blockMeshDict file that is available in the system folder of the case setup.

Geometry Creation:

Let us consider creating the geometry using 5 blocks as shown in the figure below,

NOTE: We have demonstrated the geometry creation using 5 blocks and have indexed the vertices this way. This can be varied, according to your comfort.

The Isometric View, 

Now let’s, start indexing the Vertices on the Back side of the Blocks, Since openFoam uses C++ programming, we start the index from “0”. (The Numbering of the vertices can be varied according to your comfort)

Similarly, the indexing for the frontside of the Blocks are,

To View the Bigger Picture,

Finally, to define the blocks, we need to specify the right order of vertex indices

For example:

to create the 1st Block for this case,

hex (0 1 2 3 11 12 13 14) (1 1 1) simpleGrading (1 1 1)

Similarly proceed by defining the vertices of the other blocks in the Geometry.

To View geometry that you have created, clear out all the entries inside the boundary section of the blockMeshDict file and execute the blockMesh command in the terminal.

Now open paraview and view the geometry that you have created.

Meshing:

Now to mesh the geometry fill in the number of cells in x, y and z direction inside the second parenthesis of the Hex command,

i.e., hex(0 1 2 3 11 12 13 14) (10 10 1) simplegrading(1 1 1).

since this is a 2D problem, increasing the number of cells in the z direction, doesn’t make sense. To view the mesh in paraFoam, select the surface with edges option and view the mesh.

Similarly to refine the mesh near the walls of the geometry, you need to play around with the numbers in the third parenthesis. We leave that for you to figure out, how It can be done.

Boundary Conditions:

Before we run the simulations, we need to specify the initial and boundary conditions.

The given boundary conditions are,

1.Inlet: (Type patch)

The faces that contribute to inlet boundary are shown above. If we are looking at these faces from the outside, the definition of inlet boundary will be in counter-clockwise direction,

(0 11 14 3)

(3 14 16 5)

2.Outlet: (Type patch)

Similarly, we have 3 faces contributing to the outlet boundary condition. Since we are viewing the faces from outside,  hence the definition of the outlet boundary will be in counter-clockwise direction,

(9 8 19 20)

(8 7 18 19)

(7 6 17 18)

Now for the other boundary look at the images below and repeat the steps to define the boundaries.

3.Front & Back: (Type empty)

4.No Slip walls: (Type wall)

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


convertToMeters 0.001;

vertices
(
    (0 0 0) //0
    (80 0 0) //1
    (80 5 0) //2
    (0 5 0) //3
    (0 0 1) //4
    (80 0 1) //5
    (80 5 1) //6
    (0 5 1) //7
  
       (0 10 0) //8
   (80 10 0) //9
   (0 10 1) //10
   (80 10 1) //11
    (200 5 0) //12
    (200 10 0) //13
    (200 10 1) //14
    (200 5 1) //15
   (200 0 1) //16
    (200 0 0) //17
    (200 -10 0) //18
    (80 -10 1) //19
   (200 -10 1) //20
    (80 -10 0) //21
    );
    
    


blocks
(
    hex (0 1 2 3 4 5 6 7) (200 20 1) simpleGrading (1 1 1)
     hex (3 2 9 8 7 6 11 10) (200 20 1) simpleGrading (1 1 1)
    hex (2 12 13 9 6 15 14 11) (200 20 1) simpleGrading (1 1 1)
    hex (1 17 12 2 5 16 15 6) (200 20 1) simpleGrading (1 1 1)
     hex (21 18 17 1 19 20 16 5) (200 20 1) simpleGrading (1 1 1)
     
);

edges
(
);

boundary
( 
    outlet
   {
      type patch;
      faces
      (
        
       (12 13 14 15)
       (17 12 15 16)
       (18 17 16 20)
        );
        }
        frontAndback
        {
           type empty;
           faces
           (
           (0 3 2 1)
           (3 8 9 2)
           (2 9 13 12 )
           (1 2 12 17)
           (21 1 17 18)
           (4 5 6 7)
           (7 6 11 10)
           (6 15 14 11)
           (5 16 15 6)
           (19 20 16 5)
           );
           }
           noslipwall
           {
              type wall;
              faces
              (
                 (10 11 9 8)
                 (11 14 13 9)
                 (4 0 1 5)
                  (21 19 5 1)
                 (19 21 18 20)
                 );
                 }
                 inlet
                 {
                    type patch;
                    faces
                    (
                      
                      
                      (0 4 7 3)
                      (3 7 10 8)
                      );
                      }
                 
       
   
  );  

mergePatchPairs
(
);

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

first we create blocks by giving the vertixes of the block using the dimension of the blocks in the x y and z direction. second we have called 'hexa' defines type of meshig element used and next we have number of grid points is specified along x,y z and we have simplegrading factor which is ratio of the size of the end cell to the size of the starting cell here simplegrading (1,1,1) the size of the cell are equally spaced in all three direction.

next we specify the boundary condtion at inlet ,outlet,frontandback and at no slip the inlet and outlet are specifed as patches as they are considerd as the faces where fluid enters and leaves and front and bottom faces are considered as empty  and top and bottom faces are considered as slip walls as the velocity will be zero.

in the icoFoam slover it writes out the results for each field u,p in to time directories.for this case the entries in the controlDict are shown below 

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

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         0.2;

deltaT          0.00001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

 next we set boundary condtion i.e p and u for our faces we have p file as shown below

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

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

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    noslipwall
    {
        type            zeroGradient;
    }
   
      frontAndBack
      {
        type              empty;
        }
        
}
// ************************************************************************* //

 the first line specifies the dimensions of the field here kinematic pressure i,e m^2s^-2. the boundary consists of the inlet patch for which its given as zero gradient meaning normal gradient of pressure is zero for inlet face and for the outlet the unifrom value is given and for the frontandback faces its given empty as resprent the the front and back planes of the 2d cases therefore it ste has empty.

we have velocity file has

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     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
{
    inlet
    {
        type            fixedValue;
        value           uniform (1 0 0);
    }

    outlet
    {
        type           zeroGradient;
    }

    noslipwall
    {
        type            noSlip;
    }
   
       frontAndBack
       {
         type             empty;
         }
         
           
}

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

the first line resprents the dimensions for velocity which is m/sec hence m^1 and sec^-1 the internal field is initililased as unifrom zero which is expressed in three compoments i.e unifrom (0 0 0) the inlet wall moves at spped of 1 m/sec in the x direction hence fixed value condtion also but with unifrom(1 0 0) the other faces considred as walls hence unifrom value of ( 0 0 0).

 after applying boundary condtio now we run our slover i.e icoFoam and i can be visulaised in the paraFoam

without grading factor 

 with 0.2 grading the blockMeshDict file will be

blocks
(
    hex (0 1 2 3 4 5 6 7) (200 20 1) simpleGrading (0.2 0.2 1)
     hex (3 2 9 8 7 6 11 10) (200 20 1) simpleGrading (0.2 0.2 1)
    hex (2 12 13 9 6 15 14 11) (200 20 1) simpleGrading (7 0.2 1)
    hex (1 17 12 2 5 16 15 6) (200 20 1) simpleGrading (7 0.2  1)
     hex (21 18 17 1 19 20 16 5) (200 20 1) simpleGrading (7 0.2 1)
     
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (200 20 1) simpleGrading (0.5 0.5  1)
     hex (3 2 9 8 7 6 11 10) (200 20 1) simpleGrading (0.5 0.5 1)
    hex (2 12 13 9 6 15 14 11) (200 20 1) simpleGrading (7 0.5 1)
    hex (1 17 12 2 5 16 15 6) (200 20 1) simpleGrading (7 0.5  1)
     hex (21 18 17 1 19 20 16 5) (200 20 1) simpleGrading (7 0.5 1)
     
);

 conclusion 

               as from above plot there is small amont of change in the velocity profiles but as grading size decreases we get fine mesh hence we have fine meshing at grading factor of 0.2.

ert