Challenge

Week 8 - Simulation of a backward facing step in OpenFOAM

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 solver based on the described physics of the flow in problem statement.
  • Explain 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.
  • Next, preform two case study, firstly, form simulation of the flow without using any grading factor (i.e., GF = 1) . Secondly, simulate the flow with grading factor of 0.2. The cells should be finer near the walls (including the the step wall).

Objective:-

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

Problem Statement:-

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

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 (includig the step wall).

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 key board or right click -> Open Terminal. After starting terminal yo will see following window.

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.

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 incopressible folders.

In above picture, we can see the plethora of solver folders. For this particular solver we need solver for laminar-viscous flow. So here we choose ‘icoFoam’

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.)

Therefore, we will select icoFoam folder by typing ‘cd icoFoam’ and by typing ‘ls’ we explore this folder in terminal.

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.

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’.

And then 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 specification for our block mesh in system folder.

As we can see in above picture, 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’.

Now, make following changes in 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 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 identifier by which openFoam reads to understand type of file its reading in this case it is ‘blockMeshDict’. it contain 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, its 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 meter, 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 we 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 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 following sketch.


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

Next, we are going to define block for meshing. If we look at he mesh format first letter, in block parenthesis, is ‘hex’ which define 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 grid 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

Lastly, we give simple grading which is cell expansion ratios

Grading factor or expansion ratios:-

The expansion ratios enables the mesh to be graded, or refined, in specified directions. The ratio is that of the width of the end cell δe">δeδealong one edge of a block to the width of the start cell δs">δsδ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 x1,x2andx3">x1,x2andx3x1,x2andx3directions respectively with only 3 expansion ratio. For 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 20 1) simpleGrading (1 1 1) // block1
hex (3 2 4 5 14 13 15 16) (200 20 1) simpleGrading (1 1 1) // block2
hex (2 7 6 4 13 18 17 15) (200 20 1) simpleGrading (1 1 1) // block3
hex (1 8 7 2 12 19 18 13) (200 20 1) simpleGrading (1 1 1) // block4
hex (10 9 8 1 21 20 19 12) (200 20 1) simpleGrading (1 1 1) // block5
);

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

edges
(
);

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

mergePatchPairs
(
);

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

After, editing of blockMeshDic file save it and close it.

Furthermore, we are going to edit ‘controlDict’ file to provides time limits - start time and end time. Also to represent deltat. While selection of deltat CFL value shoud be less than 1. otherwise, the solution could stopped forcibly or it may diverge at the time of converging. So make following changes in ‘controlDict’ file.

Edit codes following way

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


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

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

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

Initial setting for Velocity

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

fixedWalls
{
type noSlip;
value uniform (0 0 0);
}

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

fixedWalls
{
type zeroGradient;
}

frontAndBack
{
type empty;
}
}

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

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

blockMesh

CheckMesh

After, checking mesh we get following confirmation

At this stage, we can check our mesh is as per requirement or not by using paraView typing paraFoam command.

Take a look at above mesh. One can notice it is not uniform all over the surface. So let us just adjust the mesh grid following way.

hex (0 1 2 3 11 12 13 14) (80 5 1) simpleGrading (1 1 1) // block1
hex (3 2 4 5 14 13 15 16) (80 5 1) simpleGrading (1 1 1) // block2
hex (2 7 6 4 13 18 17 15) (120 5 1) simpleGrading (1 1 1) // block3
hex (1 8 7 2 12 19 18 13) (120 5 1) simpleGrading (1 1 1) // block4
hex (10 9 8 1 21 20 19 12) (120 10 1) simpleGrading (1 1 1) // block5

After making 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 simulation will start.

Now, let us view the results in paraView by using ‘paraFoam’ command

Mesh grid:-

Velocity Distribution along the domain area:-

Pressure Distribution along the domain area:-

To Calculate velocity at distance 0.085 from inlet modify following parameters.

Velocity Plot:-

With Grading Factor 0.2 (cell near fixed walls)

Now, lets solve same problem for grading factor 0.2

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

Output mesh from above codes:-

Velocity Distribution along the domain area:-

Pressure Distribution along the domain area:-

Velocity Plot:-

With Grading Factor 0.2 (cell near All walls)

Next,use grading factor 0.2 for all walls in the system.

blocks
(
hex (0 1 2 3 11 12 13 14) (80 5 1) // block1
simpleGrading
(
(
(0.5 0.5 5) // 50% y-dir, 50% cells, expansion = 5
(0.5 0.5 0.2) // 50% y-dir, 50% cells, expansion = 0.2 (1/5)
)
(
(0.8 0.9 5) // 80% y-dir, 90% cells, expansion = 5
(0.2 0.1 1) // 20% y-dir, 10% cells, expansion = 1
)
1 // z-direction expansion ratio
)
hex (3 2 4 5 14 13 15 16) (80 5 1) // block2
simpleGrading
(
(
(0.5 0.5 5) // 50% y-dir, 50% cells, expansion = 5
(0.5 0.5 0.2) // 50% y-dir, 50% cells, expansion = 0.2 (1/5)
)
(
(0.2 0.1 1) // 20% y-dir, 10% cells, expansion = 1
(0.8 0.9 0.2) // 80% y-dir, 90% cells, expansion = 5
)
1 // z-direction expansion ratio
)
hex (2 7 6 4 13 18 17 15) (120 5 1) // block3
simpleGrading
(
(
(0.5 0.5 5) // 50% y-dir, 50% cells, expansion = 5
(0.5 0.5 0.2) // 50% y-dir, 50% cells, expansion = 0.2 (1/5)
)
(
(0.2 0.1 1) // 20% y-dir, 10% cells, expansion = 1
(0.8 0.9 0.2) // 80% y-dir, 90% cells, expansion = 5
)
1 // z-direction expansion ratio
)
hex (1 8 7 2 12 19 18 13) (120 5 1) // block4
simpleGrading
(
(
(0.5 0.5 5) // 50% y-dir, 50% cells, expansion = 5
(0.5 0.5 0.2) // 50% y-dir, 50% cells, expansion = 0.2 (1/5)
)
(
(0.8 0.9 5) // 80% y-dir, 90% cells, expansion = 5
(0.2 0.1 1) // 20% y-dir, 10% cells, expansion = 1
)
1 // z-direction expansion ratio
)
hex (10 9 8 1 21 20 19 12) (120 10 1) // block5
simpleGrading
(
(
(0.5 0.5 5) // 50% y-dir, 50% cells, expansion = 5
(0.5 0.5 0.2) // 50% y-dir, 50% cells, expansion = 0.2 (1/5)
)
(
(0.9 0.9 5) // 90% y-dir, 90% cells, expansion = 5
(0.1 0.1 1) // 10% y-dir, 10% cells, expansion = 1
)
1 // z-direction expansion ratio
)
);

Output mesh from above codes:-

Velocity Distribution along the domain area:-

Pressure Distribution along the domain area:-

Velocity Plot:-

Results

- Up to this point, we have done plotting velocity plot for both, without grading factor and with 0.2 grading factor, cases. If we take a look at all the three line graph, we can see, the plot for velocity without grading factor seems more smooth and accurate than other two plot. So one can say, grading can disrupt out result. But it is incorrect. Because, if we take good look at it, we can notice, the mesh is fine at step wall along the X-axis but very coarse along Y-axis. Therefore, we are not getting smooth plot for velocity at 0.085 point at Y-axis.

- To understand this thing clearly, we can see below figure, where we can see- on the left side there is coarse mesh and on the right side fine mesh. Above picture, there respective graphs are presented, which can really convince anyone how at any given point mesh produce a plot. The one on the right side mesh is giving very smooth plot except one small pick near point 0.01 and on the other hand, for coarse structured mesh, the curve of graph the curve is not smooth, there are sharp corner forming at the top of the curve.

- After learning how mesh actually works, we should use fine mesh near the point where we expect our results. In CFD fluid always shows large changes near walls and obstacles along its way and to capture those changes effectively we should always use grading or finer mesh near a wall.

- To get effective and reliable results, one should know where to use finer grading factor. In a bid to obtain more effective results, suppose, if we use large number of fine elements, it could use our computation time and also require large computational memory(RAM). Therefore, use of grading is important, so that, we can use fine mesh to capture the flow at boundary wall and obstacles, course mesh at more open area.

- We can also see importance of fine mesh by observing following picture. As grading factor increasing near step, the resolution of image in improving. Means, as the number of element is increasing, the flow result is also increasing.

- We can see more fine result at grading factor 0.15. Hence, always use fine mesh in the area where maximum chnges happen in velocity, pressure or other property of fluid, to capture more detailed and accurate results.

Complete blockMeshDict file for Grading factor = 0.2

/*--------------------------------*- 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.01;

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

blocks
(
hex (0 1 2 3 11 12 13 14) (80 5 1) // block1
simpleGrading
(
(
(0.5 0.5 5) // 50% y-dir, 50% cells, expansion = 5
(0.5 0.5 0.2) // 50% y-dir, 50% cells, expansion = 0.2 (1/5)
)
(
(0.8 0.9 5) // 80% y-dir, 90% cells, expansion = 5
(0.2 0.1 1) // 20% y-dir, 10% cells, expansion = 1
)
1 // z-direction expansion ratio
)
hex (3 2 4 5 14 13 15 16) (80 5 1) // block2
simpleGrading
(
(
(0.5 0.5 5) // 50% y-dir, 50% cells, expansion = 5
(0.5 0.5 0.2) // 50% y-dir, 50% cells, expansion = 0.2 (1/5)
)
(
(0.2 0.1 1) // 20% y-dir, 10% cells, expansion = 1
(0.8 0.9 0.2) // 80% y-dir, 90% cells, expansion = 5
)
1 // z-direction expansion ratio
)
hex (2 7 6 4 13 18 17 15) (120 5 1) // block3
simpleGrading
(
(
(0.5 0.5 5) // 50% y-dir, 50% cells, expansion = 5
(0.5 0.5 0.2) // 50% y-dir, 50% cells, expansion = 0.2 (1/5)
)
(
(0.2 0.1 1) // 20% y-dir, 10% cells, expansion = 1
(0.8 0.9 0.2) // 80% y-dir, 90% cells, expansion = 5
)
1 // z-direction expansion ratio
)
hex (1 8 7 2 12 19 18 13) (120 5 1) // block4
simpleGrading
(
(
(0.5 0.5 5) // 50% y-dir, 50% cells, expansion = 5
(0.5 0.5 0.2) // 50% y-dir, 50% cells, expansion = 0.2 (1/5)
)
(
(0.8 0.9 5) // 80% y-dir, 90% cells, expansion = 5
(0.2 0.1 1) // 20% y-dir, 10% cells, expansion = 1
)
1 // z-direction expansion ratio
)
hex (10 9 8 1 21 20 19 12) (120 10 1) // block5
simpleGrading
(
(
(0.5 0.5 5) // 50% y-dir, 50% cells, expansion = 5
(0.5 0.5 0.2) // 50% y-dir, 50% cells, expansion = 0.2 (1/5)
)
(
(0.9 0.9 5) // 90% y-dir, 90% cells, expansion = 5
(0.1 0.1 1) // 10% y-dir, 10% cells, expansion = 1
)
1 // z-direction expansion ratio
)
);

edges
(
);

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

mergePatchPairs
(
);

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

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.