Week 8 - Simulation of a backward facing step in OpenFOAM

Simulation of Backward Facing Step in OpenFoam:

In this Challenge we are going to simulate an laminar-incompressible-viscous flow through a backward facing step geometry. For our Simulation we will use the tutorial problems given with the OpenFoam that uses laminar-incompressible-viscous flow. We can find this example by going to the folder as shown: openfoam8/tutorial/incompressible/icoFoam/cavity/cavity.

this example tutorial use all the same governing equation needed for the simulation. We will be changing the geometry as per our problem in this challenge and then simulate. For changing the geometry we will be going inside the system folder with our cavity folder. There we can see the blockMeshDict. we will be editing this blockMesh file to change our geometry and mesh as required.

Geometry:

First we will see our geometry for our backward facing step as given in the problem

In our case for a backward facing step so we will first create a geometry of a backward facing step. we will be using 5 blocks to create our geometry as shown below in 2D from:

we have got an idea fo how our block will be after seeing this ablove pictures now we will have to create points that will help use create this blocks in a 3D form. as we know the blockMeshDict is in C++ format the indexing for this points will start from  0. Now see the following figure to understand the indexing used to create this blocks:

Now we will be entering the coordinates of each point in our blockMeshDict

convertToMeters 1;

vertices
(
    (0 0 0)
    (0.08 0 0)
    (0.08 0.005 0)
    (0 0.005 0)
    (0.08 0.01 0)
    (0 0.01 0)
    (0.2 0.01 0)
    (0.2 0.005 0)
    (0.2 0 0)
    (0.2 -0.01 0)
    (0.08 -0.01 0)
    
    (0 0 0.01)
    (0.08 0 0.01)
    (0.08 0.005 0.01)
    (0 0.005 0.01)
    (0.08 0.01 0.01)
    (0 0.01 0.01)
    (0.2 0.01 0.01)
    (0.2 0.005 0.01)
    (0.2 0 0.01)
    (0.2 -0.01 0.01)
    (0.08 -0.01 0.01)
    
);

after creating the points now we will have to enter the code create the blocks form this points/vertices.

also first we will be using simple grading with grading factor 1 so the code for creating the blocks and meshing with simple grading is given below

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

Boundary conditions:

now we have our mesh ready, so we will be giving boundary conditions as required for our problem. Following is the code for our boundary conditions

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

inlet and outlet are of the type patch , noSlipWalls are of type wall and the frontAndBack are of type empty since our focus is in 2D from of the backward facing step.

Now we have our blockMeshDict ready. we can use blockMesh command in out terminal to check whether there are any errors in our file.

Initial Conditions:

After this we will be giving inital conditions for our problems that are stored in the folder '0' . we have two initial conditions U and p. so we will edit the velocity file first and then the pressure file.

for our velocity initial condition of uniform value 1m/s in positive x direction was given and the type of inlet was set to be fixed value. for the outlet the type was set to zeroGradient and for noSlipWalls the type was noSlip.

following is the code for the velocity initial conditions:

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

    outlet
    {
        type            zeroGradient;
    }
    
    noSlipWalls
    {
        type            noSlip;
    }

    frontAndBack
    {
        type            empty;
    }
}

similary for the pressure initial conditions were set to type zeroGradient for the inlet and noSlipWalls, uniform fixed value 0  for the outlet.

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    
    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    noSlipWalls
    {
        type            zeroGradient;
    }

    frontAndBack
    {
        type            empty;
    }
}

Time step:

now lastly we will have to change the timestep so that our simulation run fine and the Courant number does not exceed the limiting value or else the solution will blow up . We will go to the controlDict file and edit it for our required time step values. So the timestep was set to deltaT = 0.001 and the endtime for the simulation was set to 0.5  , we will be writing our file after every 20 timesteps so we have the datafile saved. 

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         0.5;

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

now we will move back to the constant folder and edit the value of kinematic viscoucity to 1e-4 . this can be done in the transportPorperties file.

nu              [0 2 -1 0 0 0 0] 0.0001;

As all the setup for our simulation is now completed we will now be simulating our problem using the icoFoam solver so we will go to the terminal at the location where we are inside our main folder for our program there we will be giving the icoFoam command in the terminal and our simulation will be running.

after the simulation is completed we can see our results in paraview , to open paraview we will use paraFoam command. Now inside paraview we will first see our mesh

Mesh

here we can see that near the backward step the elements are of same size as we have used the grading factor of 1.

Results (GF = 1) :

Velocity Contour

we will see the velocity contour of the backward step

Pressure Contour

Plot Over Line

we have to also make a plot over line at a distance of 0.085 m from the inlet for that we can see the location of our plot over line as below

there forms a zone of low pressure at near the backward step that can be clearly seen in the pressure contour.

now at this line we will see the variation of velocity and pressure along the Yposition

U-vs-Y

the x-velocity and the total vleocity are differen in the lower region of the y position indicating that the recirculations have been formed in this region.

P-vs-Y

Grading Factor = 0.2

we have seen our results with a mesh that is uniform all over the geometry now we will be changing the grading factor so that boundary regions are captured properly

Now changing the grading factor to 0.2 . accordingly we need to change our code for our blockMeshDict

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

all teh other thing will be kept same as earlier and we will run our simulation using the icoFoam command.

we can see the mesh and results in the paraview usig the paraFoam command

Mesh

from the above images we can that there is a grading in our mesh and the region near our bakward step can be captured properly.

Results (GF =0.2):

Velocity Contour:

transient 

Pressure Contour

Plot Over Line:

the location for a plot over line is agaon the same at x = 0.085 m from the inlet

V-vs-Y

the velocity distribution for the plot over line can be seen below

P-vs-Y

the pressure distributino for teh plot ove line can be seen below

From this simulation it is clear that when flow passes through a backward facing step there forms low pressure region near the step and in the region downstream the step there forms recirculations due to lower pressure and high velocity flowing from above. 

Importance of Grading/Growth Factor:

• From the pressure vs Y-position plots fro both the grading factors we can see that the plot for grading factor 0.2 is smoother and it has higher pressure range.
• If we obeserve the velocity contours of both cases we fiind that the boundary layer of the later case is smaller and crealy defined this is due to the small size of elements at the boundary.