Simulation of a backward facing step in OpenFOAM

AIM:

• Simulate an incompressible-laminar-viscous flow through the backward facing step geometry. 
• perform transient simulation; choose a solver based on the described physics of the flow in problem statement. 
• Explaining simulation process and way of setting up the problem statement in openFOAM and also, explain every changes that have done inside the reference case files obtained from OpenFOAM tutorials
• Performing two case studies, firstly, form simulation of the flow without using any grading factor (i.e., GF = 1) .

PROCEDURE TO APPROACH THE PROBLEM:

• Open openFoam tutorials and find the solver to solve incompressible-laminar-viscous flow problem.
• Next step is to copy the required file in solver and paste it in the OpenFOAM run folder so that the official file of OpenFoam remains untouched. 
• Find a block mesh script file to create domain space and provide mesh specification.
• Find ‘controlDict’ file to edit start and end time for simulation.
• In the initial folder ‘0’ edit the initial parameter for velocity and pressure according to faces of the given domain. 
• Use blockMesh command to generate the mesh and checkMesh to verify it.
• Open a simulation file in paraView and plot the velocity profile at point 0.085m from inlet. 
• Write a detailed report explaining all the commands and codes used in the project. 

Problem Statement: 

Simulate an incompressible-laminar-viscous flow through the backward facing step geometry, perform case 1.

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

Mesh specification:

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

Boundary condition specification

The domain specifications are provided in the following figure. 

In order to start the simulation process,

• First start terminal in Linux Operating system using ctrl+alt+T buttons on keyboard or right click
• Open Terminal

• 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. 
• We have to solve incompressible flow simulation for this problem so select incompressible directory by typing ‘cd incompressible’ and then again type ls to explore folders in incompressible folders. 

• In the above picture, we can see the plethora of solver folders. For this particular solver we need a solver for laminar-viscous flow.
• So here we choose ‘icoFoam’ icoFoam is the solver in openFoam which is used as a transient solver for incompressible, laminar flow of Newtonian fluids.

• Therefore, we will select icoFoam folder by typing ‘cd icoFoam’ and by typing ‘ls’ we explore this folder in terminal
• We have an 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

• In the cavity again see the 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 the tutorial section. First of all, OpenFoam does not let anyone make any changes in their official tutorial and secondly, any changes with their official tutorials may corrupt openFoam. To correct it, we may need to reinstall open foam all over again. 
• So, we are going to copy this cavity folder using the ‘cp -r cavity’ command and paste it in the openFoam run folder by typing ‘cp -r cavity $FOAM_RUN/cavity_test’. 
• Open run folder by typing ‘cd $FOAM_RUN’, where we can see cavity_test file.

• 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 specifications for our block mesh in system folder.

• As we can see in the above picture, there are 4 script file system folders and by opening blockMeshDict we can create a domain and provide required specification for block mesh.
• So, open this file by typing ‘gedit blockMeshDict’.

• Now, make the following changes in the blockMeshDict script file. 
• The blockMeshDict file is a script file that contains information regarding all aspects of the geometry. 

• When we open this file, first we are going to see header information in the form of a comment which is not executable.

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \    /   O peration     | Website:  https://openfoam.org
    \  /    A nd           | Version:  8
     \/     M anipulation  |
*---------------------------------------------------------------------------*/

• Then, we will see file identifiers by which openFoam reads to understand the type of file it is reading, in this case it is ‘blockMeshDict’. it contains information in the dictionary format, means, key and the value format.
• As we can see, version is the key and ‘2.0’ is its value, that’s why, it's called key value pair. And here we don’t have to make any changes, so we keep it as it is.

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

• Next, we set units in meters, by converting all type of length in to meter. 

convertToMeters 0.01 ;

Mesh generation without any grading factor

• In the same blockMesh file go to the lines where the word hex can be found. That is the line where we specify the faces of the blocks and also the number of cells in X, Y and Z direction.

• Also we specify the faces as the inlet, outlet and walls as shown in the picture below 

• The blockMeshDict file is edited as shown in the above. 
• The next step is to navigate to the controlDict folder and specify the simulation time as 200 seconds.The script file looks as follows

Vertices:

• Afterward, in codes, we provide a list of vertices by providing x,y,and z coordinates locations.
• its always start from vertex (0,0,0). so we have to create following geometry. 

• To create the above geometry, let us consider creating the geometry using 5 blocks as exhibited in below figure. 

• Now, we are going to provide vertices for each block by referring to the following sketch.

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

• Next, we are going to define a block for meshing. If we look at the mesh format, the first letter, in block parenthesis, is ‘hex’ which defines the type of mesh which is hexahedral mesh.
• Then, again in parenthesis, we provide vertices number according to order which follow right hand thumb rule, meaning, thumb is the z axis and order of vertices follow rotational order of our other finger
• Then comes, number of mesh grids that we want to divide along x, y, and z length. We define it as (nx,ny,nz) format.

Given mesh specification:

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

Now, we give simple grading which is cell expansion ratios

Grading factor or expansion ratios:

• The expansion ratios enable the mesh to be graded, or refined, in specified directions. The ratio is that of the width of the end cell along one edge of a block to the width of the start cell along the edge

Expansion ratio = (size of end cell) / (size of starting cell)

• The simple description specifies uniform expansion in the local directions respectively with only 3 expansion ratio. For the first case, we are not going to provide any grading for our problem.

• So here, we have 5 blocks so we have to define hex mesh for 5 blocks.

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

• Next, define the boundary of the problem, in this we have to define faces of our 3d block that should be in order according to the right hand thumb rule. So set faces for this problem the following way.

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

• After editing the blockMeshDict file save it and close it.
• Furthermore, we are going to edit the ‘controlDict’ file to provide time limits - start time and end time.

Also to represent delta t. While selection of delta t CFL value should be less than 1. otherwise, the solution could stop forcibly or it may diverge at the time of converging. So make the following changes in ‘controlDict’ file.

• Edit codes following way

• Then, we have to define the initial value for pressure and velocity. In order to do that, we have to go back to the cavity_copy directory using the ‘cd ..’ command and select a folder named ‘0’, which stands for initial values. In this case initial value for pressure and velocity.

• In this case,we have to set initial boundary conditions for pressure and velocity at each face and also update the name of the boundary same as we changed inlet outlet name in ‘blockMeshDict’ script file as shown in below codes. 

/*--------------------------------*- 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;
    }
    Noslipwalls
    {
        type            noSlip;
    }

    frontAndBack
    {
        type            empty;
    }
}

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

Initial setting for Pressure

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

    Noslipwalls
    {
        type           zeroGradient;

    }

    frontAndBack
    {
        type            empty;
    }
}

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

• At this point, our problem statement is set, we can now proceed to deploy the mesh using the ‘blockMesh’ command and then use the ‘checkMesh’ command to check if there are any errors in our mesh.

BlockMesh

CheckMesh

After, checking mesh we get following confirmation 

• At this stage, we can check if our mesh is as per requirement or not by using paraView typing paraFoam command.
• Adjust the mesh grid the following way. 

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

• After making the above changes, now the mesh is uniformly distributed along the surface. 
• After deploying mesh, now we are ready to run the simulation. Simulation that we selected was icoFoam. Hence, just type ‘icoFoam’ and the simulation will start. 
• Now, let us view the results in paraView by using ‘paraFoam’ command

RESULTS:

Mesh grid:

Velocity Distribution along the domain area:

Pressure Distribution along the domain area:

Velocity Plot:

CONCLUSION:

• Changing the grading factor the mesh becomes finer and the solution becomes very accurate. Overall the change in the size of simple grading factor (mesh size) the results remains the same as approximately. The execution time is lower for a lower meshing grading factor.