Week 8 - Simulation of a backward facing step in OpenFOAM
Aim: To simulate an incompressible, laminar, viscous flow through a backward-facing step in OpenFOAM and post-processing in Paraview.
Introduction :
Backward Facing Step (BFS) is widely known for its application in the studies of turbulence in internal flows. 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.
To start a simulation in OpenFOAM, the first step is to choose a tutorial that is similar to a case. The tutorial might not be exactly similar, so, we make necessary changes, might be in the geometry, mesh, solver settings, etc.
In our case, we are interested in simulating the flow through a backward-facing step similar to the 'Cavity' case in the OpenFOAM tutorials. First, we copy the Cavity folder(Incompressible-->icoFoam) from tutorials to the run folder, where the simulations are run. After copying, the first thing that needs to be changed is the geometry and mesh settings in the blockMeshDict file(system->blockMeshDict).
Geometry modeling:
In order to create the above geometry, we need to change the vertices in the blockMesh file. We make an assumption that the geometry is made up of 5 blocks as shown below.
The geometry is made up of 21 points. Then the blocks are to be defined using the 'hex' cmd. For example,
hex(0 1 2 3 11 12 13 14) (1 1 1) simplegrading(1 1 1)
While defining the box, we orient the box in such a way that the first two vertices(0 1) in the hex represent the direction of the x-axis. Similarly, the second and third vertices(1 2) represent the y-axis, and the first and fifth vertices(0 11) represent the z-axis. This forms the local coordinate system. Similarly, other boxes can be defined and oriented on the top or side of the boxes.
As the physics that we want to capture takes place near the change in the crossectional area(at the step), meshing has to be finer in that region. That region comprises the right side of the 0 and 1 blocks, and, the left side of the 2,3, and 4 blocks.
Meshing:
The mesh refinement can be changed by changing the values in simple grading cmd. The first value represents the size of the last cell to the size of the first cell in the x-direction. The second and third represent the same in y and z-directions. The number of elements in x,y, and z-directions is represented by the second element, next to simple grading.
Boundary conditions:
The left side of the geometry is assigned input bc whereas the right side is assigned as the outlet bc.
inlet - faces -
(3 0 11 14)
(5 3 14 16)
outlet - faces -
(7 6 17 18)
(8 7 18 19)
(9 8 19 20)
The inlet and outlet are given patch type bc which indicates that mass transfer can take place.
Whenever specifying the face notations, the order should be in such a way that the direction of the points should specify the normal using the right-handed thumb rule. The normal of each and every face should be pointing outwards.
Top - faces -
(4 5 16 15)
(6 4 15 17)
Bottom-faces-
0 1 12 11)
(10 9 20 21)
step-side - faces -
(12 1 10 21)
The top, bottom, and step-side faces are given walls type bc which indicates no mass transfer happens through them.
The remaining walls i.e, the front and back faces of the geometry are given empty. Empty indicates that the flow field is not calculated in that face or simply the flow field varies in the XY plane.
step-1 - editiing the BlockMeshDict file
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 0.01;
vertices
(
(0 0 0) // vertex 0
(8 0 0) // vertex 1
(8 0.5 0) // vertex 2
(0 0.5 0) // vertex 3
(8 1 0) // vertex 4
(0 1 0) // vertex 5
(20 1 0) // vertex 6
(20 0.5 0)// vertex 7
(20 0 0) // vertex 8
(20 -1 0) // vertex 9
(8 -1 0) // vertex 10
(0 0 1) // vertex 11
(8 0 1) // vertex 12
(8 0.5 1) // vertex 13
(0 0.5 1) // vertex 14
(8 1 1) // vertex 15
(0 1 1) // vertex 16
(20 1 1) // vertex 17
(20 0.5 1)// vertex 18
(20 0 1) // vertex 19
(20 -1 1) // vertex 20
(8 -1 1) // vertex 21
);
blocks
(
hex (0 1 2 3 11 12 13 14) (50 5 1) simpleGrading (1 1 1)
hex (3 2 4 5 14 13 15 16) (50 5 1) simpleGrading (1 1 1)
hex (2 7 6 4 13 18 17 15) (150 5 1) simpleGrading (1 1 1)
hex (1 8 7 2 12 19 18 13) (150 5 1) simpleGrading (1 1 1)
hex (10 9 8 1 21 20 19 12) (150 10 1) simpleGrading (1 1 1)
);
boundary
(
inlet
{
type patch;
faces
(
(3 0 11 14)
(5 3 14 16)
);
}
outlet
{
type patch;
faces
(
(7 6 17 18)
(8 7 18 19)
(9 8 19 20)
);
}
frontAndBack
{
type empty;
faces
(
(12 13 14 11)
(13 15 16 14)
(18 17 15 13)
(19 18 13 12)
(20 19 12 21)
(0 3 2 1)
(3 5 4 2)
(2 4 6 7)
(1 2 7 8)
(10 1 8 9)
);
}
topandbottom
{
type wall;
faces
(
(4 5 16 15)
(6 4 15 17)
(0 1 12 11)
(10 9 20 21)
);
}
stepside
{
type wall;
faces
(
(1 10 21 12)
);
}
);
// ************************************************************************* //After changing the vertices, creating blocks, meshing, and assigning the boundary faces in the blockMeshDist file. Paraview is opened, which is a 3D data visualization software, which shows the geometry, the mesh created and also plots the contours of the pressure, velocity, temperature fields, etc at different time steps(based on the solving).
Geometry :
backward facing step 3D view
meshing with no grading factor(case1)
mesh refined case (gf=0.2)
BlockmeshDict file for case-2:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 0.01;
vertices
(
(0 0 0) // vertex 0
(8 0 0) // vertex 1
(8 0.5 0) // vertex 2
(0 0.5 0) // vertex 3
(8 1 0) // vertex 4
(0 1 0) // vertex 5
(20 1 0) // vertex 6
(20 0.5 0)// vertex 7
(20 0 0) // vertex 8
(20 -1 0) // vertex 9
(8 -1 0) // vertex 10
(0 0 1) // vertex 11
(8 0 1) // vertex 12
(8 0.5 1) // vertex 13
(0 0.5 1) // vertex 14
(8 1 1) // vertex 15
(0 1 1) // vertex 16
(20 1 1) // vertex 17
(20 0.5 1)// vertex 18
(20 0 1) // vertex 19
(20 -1 1) // vertex 20
(8 -1 1) // vertex 21
);
blocks
(
hex (0 1 2 3 11 12 13 14) (50 5 1) simpleGrading (0.2 0.2 1)
hex (3 2 4 5 14 13 15 16) (50 5 1) simpleGrading (0.2 5 1)
hex (2 7 6 4 13 18 17 15) (150 5 1) simpleGrading (5 5 1)
hex (1 8 7 2 12 19 18 13) (150 5 1) simpleGrading (5 0.2 1)
hex (10 9 8 1 21 20 19 12) (150 10 1) simpleGrading (5 0.2 1)
);
boundary
(
inlet
{
type patch;
faces
(
(3 0 11 14)
(5 3 14 16)
);
}
outlet
{
type patch;
faces
(
(7 6 17 18)
(8 7 18 19)
(9 8 19 20)
);
}
frontAndBack
{
type empty;
faces
(
(12 13 14 11)
(13 15 16 14)
(18 17 15 13)
(19 18 13 12)
(20 19 12 21)
(0 3 2 1)
(3 5 4 2)
(2 4 6 7)
(1 2 7 8)
(10 1 8 9)
);
}
topandbottom
{
type wall;
faces
(
(4 5 16 15)
(6 4 15 17)
(0 1 12 11)
(10 9 20 21)
);
}
stepside
{
type wall;
faces
(
(1 10 21 12)
);
}
);
// ************************************************************************* //As the initial pressure and velocity fields are defined for the previous geometry. We need to change the faces assigned for different boundaries and the boundary names(according to the case). In the pressure and velocity files, the units are also mentioned, OpenFOAM takes pressure values input as the kinematic pressure i.e, static pressure/density. The initial pressure values at the inner volumes are given as 0 shown by internalField uniform 0. Then to calculate the pressure on the inlet face, the zeroGradient type method, which extrapolates the pressure value from the inner volumes, is used. The outlet pressure value is given as zero gauge pressure. The walls which include top, bottom, and stepside faces are given a zeroGradient type boundary condition.
Pressure :
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ 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;
}
topandbottom
{
type zeroGradient;
}
stepside
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //Similarly, in the velocity file, the faces and the boundary names are changed. Then the values at the inner volumes are given as 0 indicated by internalField uniform(000) (as pressure is a vector). The inlet velocity is calculated based on the flow type i.e, laminar flow(low Re flow,10)
`Re = (v .L)/nu`
where v- inlet velocity to be calculated = 10m/s
L - characteristic length = 0.01m
nu - kinematic viscosity = 0.01m^2/s
The outlet velocity values are calculated using the zeroGradient method. The walls which include the top, bottom, and stepside faces are given noSlip condition, which means the velocity value at the adjacent fluid layer to the wall is zero.
velocity:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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 (10 0 0);
}
outlet
{
type zeroGradient;
}
topandbottom
{
type noSlip;
}
stepside
{
type noSlip;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //Then the transport Properties which is present in the constant folder are changed.The value of kinematic viscosity is changed to 0.01 m^2/s.
Transport properties file:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object physicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 1e-2;
// ************************************************************************* //After changing the parameters which simulate the flow, the simulation run time, time steps, etc are changed in the controlDict file.
For case 1, the time-step size is given as 1e-3.
For case 2, as the mesh size decreases, the time step size needs to be decreased as well to maintain the CFL criterion. So, the time-step size is given out to be 5e-5.
controlDict file:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 2;
deltaT 1e-3;
writeControl runTime;
writeInterval 0.1;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Post-processing:
Initial velocity
Initial pressure:
velocity contours :
case1
case 2
pressure contours:
case1
case2
Velocity profiles before and after the step are plotted below:
case1
case2
Observations and Conclusion:
1. From the velocity profiles, we can observe that the velocity near the wall is almost 0(due to the no-slip condition), and the velocity gradually increases i.e, it follows a parabolic curve.
2. When the transition happens near the step, the velocity decreases in the upper section of the duct(as the area increases), and the profile followed in the lower section is also a parabola, but it doesn't reach the maximum velocity as it reaches before the step.
3. When the mesh size was decreased, the curve becomes much smoother/accuracy increases.