Week 8 - Simulation of a backward facing step in OpenFOAM

AIM

To Simulate an incompressible-laminar-viscous flow through the backward facing step geometry.

 

OBJECTIVE

1. To perform a transient simulation to simulate an incompressible-laminar-viscous flow through the backward facing step geometry.
2.  Explain the entire simulation procedure (how you set up the case) and also provide explanations to every change you make inside the reference case files obtained from OpenFOAM tutorials.
3.  Perform 2 case studies as described below and compare the results.

                        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 (including the step wall).
     4.  Measure the velocity profile at 0.085 m from the inlet of the geometry for both case studies and compare.
     5.  Explain why Grading/Growth factor is essential for meshing. 

 

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

Desired output

• Velocity profile at 0.085 m from the inlet of the geometry for both case studies and compare.
• The reason why Grading/Growth factor is essential for meshing. 

 

 

INTRODUCTION

Studying the flow behaviour through various structures is an important subject for a CFD engineer.  When fluids flow through different geometries they exhibit different velocity, pressure and flow properties. in this project, simulation of a flow through backward facing step is done by studying the flow of fluids through cemeteries that have sudden-changing shapes. The simulations are refined using multiple meshes - mesh grading factors inside the blockMeshDict file of blockMesh which is also used to create the geometric to be solved in OpenFOAM.

 

STEPS TO BE DONE

1. Choosing a solver
2. Selecting a solver tutorial to create geometry and meshing.
3. Setting up case
• Assigning boundary conditions
• Assigning initial conditions
4. Post processing the result
5. Plotting the result

 

THEORY

Backward Facing Step (BFS) is widely known for its application in the studies on internal flows. Flow separation is a common and interesting phenomenon in fluid mechanics with significant effects in practical application. Thus, it has been the focus of intensive study for many years. A standard set of simple geometric configurations including the flow in a pipe with a sudden expansion, the flow in a pipe with an obstruction such as an orifice, the flow over an obstacle in a channel and the flow in a channel with a sudden expansion, has been proposed as representative test cases.

Backward Facing Step is one of the vital flow geometry for investigating separated flows. In this the flow enters from the inlet passes through a constant area duct and exits into another duct having a cross-sectional area larger than the inlet duct. As this occurs, a recirculation zone occurs at the 'step’ of the duct having a larger cross-sectional area. The flow separation is caused due to the sudden changes in the geometry. This creates a zone of re-circulation and a point of flow reattachment. Even though it is a simple and common problem but it has a wide scope, which makes it a building block for a variety of applications. This recirculation depends on the Reynolds number of the flow, size of the step, expansion ratio, aspect ratio (in 3D case), etc.

 

Schematics of the Flow over a Backward Facing Step Geometry is shown in the diagram below, -

 

Boundary layer Theory

Boundary layer, in fluid mechanics, thin layer of a flowing gas or liquid in contact with a surface such as that of an airplane wing or of the inside of a pipe. The fluid in the boundary layer is subjected to shearing forces. A range of velocities exists across the boundary layer from maximum to zero, provided the fluid is in contact with the surface. Boundary layers are thinner at the leading edge of an aircraft wing and thicker toward the trailing edge. The flow in such boundary layers is generally laminar at the leading or upstream portion and turbulent in the trailing or downstream portion.

Boundary layer separation over a flat plate :

• It is the locus of the point at which the velocity of fluid becomes equal to 99% of the free steam velocity i.e, (u0 or u-infinity).
• The free stream flow which has uniform velocity u0 when passes over a solid surface there is zero velocity (no slip) at the solid surface and gradually increases upto 99% of free stream velocity normal to the flow direction.
• Within the boundary layer fluid is viscous and velocity gradient i.e,(dy/dx) setup thats why the shear stress are present the flow will be rotational in nature.
• beyond the boundary layer the fluid is viscous but there is no velocity gradient that's why shear stress = 0 and it behave like a ideal flow.
• Some distance from the leading edge the boundary layer is laminar in nature and velocity profile is parabolic in character.
• After a certain point there is transition zone in which the flow is becomes unstable and after a particular distance it switch to turbulent flow boundary layer here the velocity profile is parabolic in nature but the turbulent flow boundary layer is logarithmic in character.

.

 

Incompressible flow : In fluid dynamics, incompressible flow refers to a flow in which the density remains constant in any fluid parcel, i.e. any infinitesimal volume of fluid moving in the flow.

 Laminar flow : Fluid (gas or liquid) flow in which the fluid travels smoothly or in regular paths .In laminar flow, sometimes called streamline flow, the velocity, pressure, and other flow properties at each point in the fluid remain constant. Laminar flow over a horizontal surface may be thought of as consisting of thin layers, or laminae, all parallel to each other. The fluid in contact with the horizontal surface is stationary, but all the other layers slide over each other.

 

Viscous Flow : The fluid flow in which frictional effects become signification, are treated as viscous flow. When two fluid layers move relatively to each other, frictional force develops between them which is quantified by the fluid property 'viscosity'. Boundary layer flows are the example of viscous flow.

The flow which contains all these effects is called incompressible -laminar -viscous flow.

 

No-Slip, and Boundary Layer

• No-slip is the boundary condition, which states that the relative velocity of fluid layers at the wall of duct/pipe is zero.
• It is due to the surface roughness of the wall and the adhesion effect between fluid and solid surfaces.
• This makes the velocity of the fluid at the wall zero and in the transverse direction of the flow, velocity varies.
• Velocity variation in the transverse direction is either parabolic or linear for laminar flow
• Velocity variation in the transverse direction is either logarithmic or follows the power law for turbulent flow.
• And due to No-slip, there is momentum disturbance in the transverse direction of the flow, which causes the formation of a boundary layer near the wall.
• Kinematic viscosity is the property of the fluid responsible for the momentum diffusivity, and the effect of viscosity is dominated in the flow.
• In solving viscous dominated flow only Navier-stokes equation is used.

 

Grading Factor:

The expansion ratio enables the mesh to be graded or refined in specified direction. The ratio is that of the width of the end cell.

δ_e – Along one edge of a block to the width of the start cell.

Δ_s – Along the edge as shown in below figure

 

In simple format,

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

The simple description specifies uniform expansion in the local x₁, x₂, and x₃. Directions respectively with only 3 expansion ratio. For first case, we are not going to provide any grading for our problem.

 

Short commands that we are using this project:-

•  ls - list of all files that stored in the directories.
• cd - To change directory
• cd .. - To go back in the previous directory.
• pwd - Present working directory.
• mkdir - To create directory.
• cp - To copy
• tut - tutorials of openfoam
• app - applications
• rm - to remove files 
• rm -r - to remove directories
• gedit -  used to open text editors
• sol -  solvers
• $FOAM_RUN - navigates to openfoam run directory
• paraFoam - opens paraview where the simulation is represented.

 

 

PROCEDURE.

• First, we are creating the run folder in OpenFOAM in which we will be running our simulations.
• And then need to copy the cavity folder to the run folder so that we can make the changes in the blockMeshDict, controlDict, and pressure and velocity files as to run our required simulations, we copying the folder because we are not supposed to make any changes in the original folders in the OpenFOAM.
• Then will make the required changes in the blockMeshDict file that is naming the vertices for and then creating the blocks, in this case, there are 21 vertices and 5 blocks while creating blocks we also need to specify the number of mesh points needed along the x, y, and z-direction and also the simple grading should be done along the boundaries.
• then, we need to assign the boundary conditions here, we need to define all the faces in OpenFOAM refers to patches and giving them type that is the inlet and outlet will be the faces (type patch), for front and back type will be empty faces, for top and bottom which are the boundary conditions in this case type will be wall type and setting up these will constitute a particular boundary, then will save these changes in the blockMeshDict.
•  Then, in the velocity and pressure files, we have to give the value of the velocity at the inlet type will be fixed value here we should assign some value as initial velocity and pressure will be zerogradient and at the outlet velocity type will be zerogradient and pressure will be type fixedvalue that value will be uniform 0, then at front and back velocity type empty and pressure, as well as type empty at the no slip boundaries that is top and bottom velocity, will be type noslip and pressure type empty. And in the transportProperties we can change the nu value in order to obtain the lower courant number.
• And in the controlDict file the simulation time and the time step need be assigned in the way that the courant number should be less than 0.5.
• and then creating the block mesh and then run the simulation using the respective solver after the completion of simulation open the paraview in the paraFoam and plot the velocity profile for the given value

 

Copying a Case folder from the OpenFOAM tutorial

As mentioned earlier, we have decided to use icoFoam solver for this project. Hence the case folder corresponding to this solver is cavity case.

Steps to copy a case folder from OpenFOAM tutorial:

1. Open the Linux terminal via ubuntu.
2. Run the tut command in the terminal. Now the path changes to OpenFOAM Tutorials.

    3. Run the ls command. It lists the variety of flow problems available in the OpenFOAM tutorial.

   4.  Run cd incompressible/ command to change the directory to incompressible flow folder. ls

   5.  Run ls command to display the list of solvers available in OpenFOAM tutorial for incompressible flow.

   6. Run cd icoFoam/ command to change the directory to icoFoam solver folder. 

  7. Run ls command to display the list of projects available for this solver.cd

  8. Run cd cavity/ command, then run ls command to display the list of cases available for this project and then run cd cavity/ again.

  9. Run pwd command to display the current directory path and then copy this path.

 10.  Run cp -r cavity/ $FOAM_RUN/week8_case1. command to copy the case folder from the OpenFOAM tutorial to our run home directory. 

11. Now we can give whatever the name we wish for this case folder. For this project, I have created two copies of cavity case. One for CASE 1 (without grading factor) and another for CASE 2 (with grading factor).

 

2. Geometry Creation

This section describes procedures being followed to create a geometry for backward facing step.

1. From case folder (e.g. week8_case1), go to system --> blockMeshDict file.
2. Open the blockMeshDict file in Notepad or anyother text editor. You can also run gedit command to open the blockMeshDict file in Ubuntu's text editor.

   3.  Dimensions of the domain are

   4.  OpenFOAM allows to create geometry in 3D cartesian coordinates system using vertices.

   5.  Let us consider creating the geometry using 5 blocks as shown below:

  6.  Let us start indexing the vertices for the blocks created. Since OpenFOAM uses C++ programming, we start indexing from '0'. 

  7.  Open the blockMeshDict file. In the vertices dictionary, the following values of the vertices coordinates are specified in the file:

   8.   In the blocks dictionary, 

• CASE 1: Without grading factor

• CASE 2: With Grading factor 0.2

  9.  In the boundary dictionary, the values for inlet, outlet and walls are specified as given below: 

 

   10.   Then mergePatchPairs dictionary is created and left empty inside. This dictionary is responsible for merging the above created blocks and patches (such as inlet, outlet and walls).

 

   11.   Save the file for both the cases (CASE 1 and CASE 2).

 

C++ Code (Grading Factor 1)

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

convertToMeters 1;

vertices
(
    // vertices from point 0 to 10
    (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)
    
    // vertices from point 11 to 21
    (0 0 0.1)
    (0.08 0 0.1)
    (0.08 0.005 0.1)
    (0 0.005 0.1)
    (0.08 0.01 0.1)
    (0 0.01 0.1)
    (0.2 0.01 0.1)
    (0.2 0.005 0.1)
    (0.2 0 0.1)
    (0.2 -0.01 0.1)
    (0.08 -0.01 0.1)
);

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

edges
(
);

boundary
(
    inlet
    {
        type patch;
        faces
        (
            (11 14 3 0)
            (14 16 5 3)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (6 17 18 7)
            (8 7 18 19)
            (9 8 19 20)
        );
    }
    
    otherwall  //no slip walls
    {	type wall;
        faces
        (
        	//top
            (15 4 5 16)
            (15 17 6 4)
            //step
            (21 12 1 10)
            //bottom left
            (11 0 1 12)
            //bottom right
            (21 10 9 20) 
        );
    }
    frontAndBack
    {
        type empty;
        faces
        (	//back
            (0 3 2 1)
            (3 5 4 2)
            (2 4 6 7)
            (1 2 7 8)
            (10 1 8 9)
            //front
            (13 15 16 14)
            (12 13 14 11)
            (20 19 12 21)
            (19 18 13 12)
            (18 17 15 13)
        );
    }
);
mergePatchPairs
(
);
// ************************************************************************* //

 

C++ Code ( Grading Factor = 0.2)

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

convertToMeters 1;

vertices
(
    // vertices from point 0 to 10
    (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)
   
    // vertices from point 11 to 21
    (0 0 0.1)
    (0.08 0 0.1)
    (0.08 0.005 0.1)
    (0 0.005 0.1)
    (0.08 0.01 0.1)
    (0 0.01 0.1)
    (0.2 0.01 0.1)
    (0.2 0.005 0.1)
    (0.2 0 0.1)
    (0.2 -0.01 0.1)
    (0.08 -0.01 0.1)
);

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

edges
(
);

boundary
(
    inlet
    {
        type patch;
        faces
        (
            (11 14 3 0)
            (14 16 5 3)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (6 17 18 7)
            (8 7 18 19)
            (9 8 19 20)
        );
    }
   
    otherwall  //no slip walls
    { type wall;
        faces
        (
        //top
            (15 4 5 16)
            (15 17 6 4)
            //step
            (21 12 1 10)
            //bottom left
            (11 0 1 12)
            //bottom right
            (21 10 9 20)
        );
    }
    frontAndBack
    {
        type empty;
        faces
        ( //back
            (0 3 2 1)
            (3 5 4 2)
            (2 4 6 7)
            (1 2 7 8)
            (10 1 8 9)
            //front
            (13 15 16 14)
            (12 13 14 11)
            (20 19 12 21)
            (19 18 13 12)
            (18 17 15 13)
        );
    }
);
mergePatchPairs
(
);
// ************************************************************************* //

 

3. Mesh Generation

1. Now come to Ubuntu terminal, run blockMesh command to generate the mesh and run checkMesh command to verify the mesh has been created correctly without any error.
2. Run touch week8_case2.foam for CASE 2 (and touch week8_case1.foam for CASE 1) in the terminal to view the generated mesh. The generated mesh for both cases are as follows:

 

4. Initial and Boundary Conditions

Now we need to specify the initial and boundary conditions such as

• Kinematic Pressure (p) and Velocity (U) at inlet and outlet at time t=0--> P and U file
• Transport properties which is the kinematic Viscosity ν --> physicalProperties file
• Start time, End time and Simulation duration --> controlDict file

1. From case folder, go to 0--> Pfile:

 

C++ Code (p File)

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \    /   O peration     | Website:  https://openfoam.org
    \  /    A nd           | Version:  9
     \/     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    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;
    }

    frontAndBack
    {
        type            empty;
    }
    otherwall
    {
    	type 		zeroGradient;
    }
}

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

 

2. Then open the U file:

 

C++ Code (U File)

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

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0.1 0 0);

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

    outlet
    {
        type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
    otherwall
    {
    	type 		noSlip; // zero at walls
    }
}

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

 

3. Then go to constant--> physicalProperties file, kinematic viscosity is taken as constant through out the domain.

 

C++ Code (physical properties)

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

nu              [0 2 -1 0 0 0 0] 1e-5;


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

 

4. And finally go to system--> controlDictfile

 

C++ Code (ControlDict)

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

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         5;

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

//adaptive timestep
/*
adjustTimeStep true;

maxCo 1;*/


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

5. Save all the files.

 

5. Running the icoFoam solver

Now come to Ubuntu Terminal and run the icoFoam solver.

Here we could see the continuity and momentum equation is solved at each time step for pressure and velocity.

In Linux, we can solve this kind of problems and get output in terms of raw data.

Hence we need data visualization tool to view the simulation and plot the result. 

 

6. Post-processing using Paraview

The output we get from running the icoFoam solver in OpenFOAM can be visualized in Paraview. 

1. In OpenFOAM, simply run the command touch week8_case1.foam.
2. Paraview opens. Click Apply button.

 

CASE 1: Without the grading factor - Uniform mesh

Uniform Mesh:

 

Pressure:

 

Velocity:

 

Velocity Plot at length 0.085 from the inlet:

To plot the velocity at length 0.085 from the inlet,

    1.  Go to Filters menu --> Data Analysis --> Plot Over Line.

 

     2.  Set the length at which we need to visualize the velocity plot.

 

 CASE 2: With Grading factor 0.2 - Non uniform mesh

Non uniform Mesh: Grading factor 0.2

 

Pressure:

 

Velocity:

 

Velocity Plot at length 0.085 from the inlet:

Why is grading / growth factor essential for meshing ?

            Mesh Grading basically refers to a technique called mesh stretching, which helps us to stretch the mesh along any particular direction (x , y or z direction). This helps us in increasing the computational efficiency of the simulations.

`"Mesh Grade"="Size of End Cell" /"Size of Starting Cell"`

 

The goal of meshing in CFD is to obtain the best mesh possible without sacrificing reasonable accuracy. A finer mesh at the areas where the fluid shows large changes, leads to higher accuracy. Thereby, we use a finer mesh around the walls of the geometry and also around the point where we expect our results.

 

Coarse and Fine Mesh

A coarse mesh is the one where mesh elements are kept bigger such that they accommodate the geometry with less number of total mesh elements. Coarse mesh is generally what the CFD/FEA simulation is initially run with and then the mesh is refined further until we reach the point of validation where the results are no longer mesh dependent.

A refined mesh is basically a finer mesh with smaller elements. When the geometry has features like sharp angles, edges, etc. a relatively finer mesh is required to capture them.

 

 

CONCLUSION 

• Hence we have simulated an incompressible-laminar-viscous flow through the backward-facing step geometry.
• We can notice that the flow through the inlet section is laminar until the step wall is encountered. The velocity of the flow in the inlet region is high at the centre and progressively decreases going towards the walls. When the flow encounters the step wall for the first time, there are minor vortices formed near the step wall and at the top of the wall which dissipates as the flow approaches steady state. Slowly the flow approaches steady state. There is no significant change in flow after it attains steady state.
• By changing the simple grading factor from the value of 0.2 to 1 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 lower meshing grading factor.