Week 8 - Simulation of a backward facing step in OpenFOAM
OBJECTIVE:
- To simulate an incompressible laminar viscous flow through the backwards facing step geometry and perform a transient simulation.
- Case 1-To simulate the flow without using any grading factor (GF =1)
- Case 2-To simulate the flow with a grading factor (GF=0.2). The cells should be finer near the walls(including the step walls).
Incompressible Laminar viscous flow.
Incompressible flow:
The flow in whch there is no change in the density is called the Incompressible flow.
Laminar flow:
A Laminar flow is one where the fluid motion is regular is no mixing.It can be visualised by the pathlines that do not intersect each other.
In fluid dynamics, laminar flow is characterized by fluid particles following smooth paths in layers, with each layer moving smoothly past the adjacent layers with little or no mixing.
At low velocities, the fluid tends to flow without lateral mixing and adjacent layers slide past one another. There are no cross-currents perpendicular to the direction of flow nor eddies or swrils of fluids.
Mathematically, Reynold number is the ratio of inertial force to the viscous force.
Viscous Flow:
A type of fluid flow in which there is a continuous steady motion of the particles, the motion at a fixed point always remains constant. Also called streamline flow, laminar flow, Steady flow.
Mesh grading
Grading Factor`GF=(size of the last cell)/(size of the first cell)`
The domain specifications are given in question
In order to start the simulation process, First start terminal in linux OS using right click > open Terminal
Now, we have to simulate our problem in openFoam in tutorials, type 'tut' to open openFoam tutorial directory. And typle 'ls' to identify files present inside the 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 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 incompressible, laminar flow of Newtoninan fluids.
Therefore, we will select icoFoam folder by typing 'cd icoFoam' and by typing 'ls' we explore this folder in terminal.
Here, We have options cavity and elbow. They both are part of the icoFoam solver and cavity section use to create required domain and blockMesh so select 'Cavity'.
In cavity,We are again going to see cavity folder, Where we have to define our fluid domain and generate blockMesh.
But we are not going to make any changes here in tutorial section. openFoam doesnot let anyone to make changes in their official tutorial and secondly any changes with their offical tutorials may corrupt openFoam.
So, we are going to copy this cavity folder using 'cp -r cavity' command and paste it in openFoam run folder ny typing ;cp - cavity $FOAM_RUN/cavity_test' .
Then open run folder by typing 'run', We can cavity_test file in it.
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 blockMesh in system folder.
As we can see in above picture, there are 4 scrpit file system folder and by opening blockMeshDict we can domain and provide required specification for blockMesh. 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 geometry.
When we open this file, first we are going to see header information in the form of comment which is not executable
Then, we will see file identifier by which openFoam reads to understand type of file its reading in this case it is 'blockMeshDict'. It contains information in the dicitonary format, means key and the value format. As we can see version is the key and '2.0' is its value thats why it is called key value pair, we dont make any changes here
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 blockMeshDict file go to the lines where the word hex can be found. That is the line where we specify the faces of the bocks and also the number of cells in X,Y and Z direction.
Vertices:
Afterwards in codes we provide list of vertices by providing x,y and z co ordinates 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)
(8 0 0)
(8 0.5 0)
(0 0.5 0)
(8 1 0)
(0 1 0)
(20 1 0)
(20 0.5 0)
(20 0 0)
(20 -1 0)
(8 -1 0)
(0 0 0.1)
(8 0 0.1)
(8 0.5 0.1)
(0 0.5 0.1)
(8 1 0.1)
(0 1 0.1)
(20 1 0.1)
(20 0.5 0.1)
(20 0 0.1)
(20 -1 0.1)
(8 -1 0.1)
);Next we are going to define block for meshing. If we look at the mesh format first letter in block parenthesis in "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.
Given mesh specification:
1.Number of cells along the x direction =200
2.Number of cells along the longer y direction =20
Grading Factor(GF):
The grading factor enables the mesh to be graded or refined in specified directions.
The simple description specifies unifrom expansion in the local x1,x2 and x3 direction respectively, In first case we not going to provide any GF 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) (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)
);Next, define boundary of the problem, in this we have to define faces of the our 3rd 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 blockMeshDict file save it and close it.
Furthermore, we are going to edit 'controlDict' file to provide time limits -start time and end time. Also to represent deltat. While selection of deltat CFL value should be less than 1. Otherwise, the solution could stopped forcibly or it may diverge at the time of coverging, So make followig changes in'controlDict' file.
Edit 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 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 stands for initial values 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 follwing confirmation
At this stage, we can check our mesh is as per requirement or not by using paraview by typing paraFoam command
If mesh is not unifrom then mesh is made uniform all over the surface by adjusting the grid
Now the mesh will be as shown in the figure below
After deploying mesh, now we are ready to run the simulation. Simulation that we selected was icoFoam. Hence 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 the inlet modify following parameters.
Velocity plot:-
With grading Factor (GF=0.2) [cell near fixed walls]
Now let us 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)
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 5 1)
);Output mesh from above code
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)
simpleGrading
(
(
(0.5 0.5 5)
(0.5 0.5 0.2)
)
(
(0.8 0.9 5)
(0.2 0.1 1)
)
1
)
hex (3 2 4 5 14 13 15 16) (80 5 1)
simpleGrading
(
(
(0.5 0.5 5)
(0.5 0.5 0.2)
)
(
(0.2 0.1 1)
(0.8 0.9 0.2)
)
1
)
hex (2 7 6 4 13 18 17 15) (120 5 1)
simpleGrading
(
(
(0.5 0.5 5)
(0.5 0.5 0.2)
)
(
(0.2 0.1 1)
(0.8 0.9 0.2)
)
1
)
hex (1 8 7 2 12 19 18 13) (120 5 1)
simpleGrading
(
(
(0.5 0.5 5)
(0.5 0.5 0.2)
)
(
(0.8 0.9 5)
(0.2 0.1 1)
)
1
)
hex (10 9 8 1 21 20 19 12) (120 10 1)
simpleGrading
(
(
(0.5 0.5 5)
(0.5 0.5 0.2)
)
(
(0.9 0.9 5)
(0.1 0.1 1)
)
1
)
);Output mesh from above code
Velocity Distribution along the domain area:-
Pressure Distribution along the domain area:-
Velocity Plot:
RESULTS:
1. We have done velocity plotting for both without grading factor and with (GF=0.2). If we take a look at all the three line graph, we can see the plot for velocity without G.F seems more smooth and accurate than other two plot. Therefor grading can disrupt output result. 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.085point at Y-axis.
2. If the mesh id not fine it gives average results because of capturing less elements at each node, To capture accurate results its required for finer mesh to capture flow properties of fluid. By comparing both the plots G.F=0.2 gives better reults than simple grading 1
3.Velocity magnitude profiles gives similar results for both cases in general at 0.085m there is an decrease in velocity but with a mesh grading factor finer results are obtained with accurate results at step height and boundary walls.
4.A graded mesh determines the element density in each region of the model based on the complexity of the geometry.More complex geometry with more angles and cuts will have a higher density of elements to better capture the topographical features. The mesh density grades away from the complex regions.