Simulation of Flow over Backward Facing Step using OpenFOAM
Objective
The objective of the project is to simulate the flow over a Backward Facing Step using OpenFOAM, and to determine the effect of mesh grading factor on the magnitude of velocity of the flow profile.
The above figure shows the geometry of the backward facing step through which the flow is to be simulated. The figure also shows the boundary conditions specifications to be applied to the geometry. The mesh specifications for meshing the above geometry are given below.
Mesh specification
- Number of cells along the x direction (longer dimension) = 200
- Number of cells along the y direction = 10
- Use a graded mesh. Grading factor near all walls should be 0.2
Meshing the Geometry
In order to mesh the geometry according to the given specifications, the geometry should be divided into several parts to make the discretization easier.
The given geometry is separated into 5 blocks as shown below.
The isometric view, and the front and back views of the geometry are given below.
The blocks are indexed accordingly as seen in the figure above.
Thus, the backward facing step geometry is divided into the following blocks with their respective indices, as shown below.
The geometry is created and meshed by modifying the BlockMeshDict file in OpenFOAM, according to the Mesh Specifications given above.
The BlockMeshDict file with the code for modifying the geometry and creating the mesh is given below.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0) // vertex 0
(0.08 0 0) // vertex 1
(0.08 0.005 0) // vertex 2
(0 0.005 0) // vertex 3
(0.08 0.01 0) // vertex 4
(0 0.01 0) // vertex 5
(0.2 0.01 0) // vertex 6
(0.2 0.005 0) // vertex 7
(0.2 0 0) // vertex 8
(0.2 -0.01 0) // vertex 9
(0.08 -0.01 0) // vertex 10
(0 0 0.001) // vertex 11
(0.08 0 0.001) // vertex 12
(0.08 0.005 0.001) // vertex 13
(0 0.005 0.001) // vertex 14
(0.08 0.01 0.001) // vertex 15
(0 0.01 0.001) // vertex 16
(0.2 0.01 0.001) // vertex 17
(0.2 0.005 0.001) // vertex 18
(0.2 0 0.001) // vertex 19
(0.2 -0.01 0.001) // vertex 20
(0.08 -0.01 0.001) // vertex 21
);
blocks
(
hex (0 1 2 3 11 12 13 14) (200 10 1) simpleGrading (1 1 1)
hex (3 2 4 5 14 13 15 16) (200 10 1) simpleGrading (1 1 1)
hex (2 7 6 4 13 18 17 15) (200 10 1) simpleGrading (1 1 1)
hex (1 8 7 2 12 19 18 13) (200 10 1) simpleGrading (1 1 1)
hex (10 9 8 1 21 20 19 12) (200 10 1) simpleGrading (1 1 1)
);
edges
(
);
boundary
(
inlet
{
type patch;
faces
(
(0 11 14 3)
(3 14 16 5)
);
}
outlet
{
type patch;
faces
(
(9 8 19 20)
(8 7 18 19)
(7 6 17 18)
);
}
frontAndBack
{
type empty;
faces
(
(0 3 2 1)
(3 5 4 2)
(4 6 7 2)
(2 7 8 1)
(1 8 9 10)
(11 12 13 14)
(14 13 15 16)
(13 18 17 15)
(12 19 18 13)
(21 20 19 12)
);
}
noSlipWalls
{
type wall;
faces
(
(16 15 4 5)
(15 17 6 4)
(10 21 12 1)
(11 0 1 12)
(10 9 20 21)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //The above BlockMeshDict file specifies the vertices, coordinate axes, number of mesh points, and the mesh grading factor for the given geometry. There are 200 cells in the x direction, 10 cells in the y direction and 1 cell in the z direction.
Mesh grading factor represents the amount by which the mesh cells are stretched in a particular direction (mesh stretching). Mesh graind factor is given by the following ratio.
Mesh Grading Factor (for a particular axis) = Size of end cell/Size of starting cell
In the above BlockMesh code, the mesh grading factor is given as 1, which indicates that the sizes of end cell and starting cell are the same. i.e. no mesh stretching. The mesh grading factor is modified accordingly to obtain finer or coarser grids near specific regions (in this case, the boundary layer) to better capture the change in flow parameters at the region, and to determine the effect of mesh grading factor on the results. The geometry is meshed with hexahedral elements with the given specifications.
Boundary Conditions Specification
The above code also provides information about the boundaries of the geometry. The inlet and outlet faces are defined. The walls of the geometry are identified and No Slip boundary condition is given to them to satisy the no slip boundary layer condition. Since the given problem involves a 2D Simulation, the front and back walls of the geometry are specified as empty.
Initial Conditions Specification
Initial Conditions are the conditions given at the starting time step from which the solver computes the values at later timesteps. The initial conditions for the Velocity (U) and the Pressure (p) are to be specified to setup the simulation. These values are specified in the '0' folder in the modified tutorial file in OpenFOAM.
Initial Velocity Conditions
The initial conditions for the flow velocity (U) is given below.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ 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;
}
noSlipWalls
{
type noSlip;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //The above code indicates that a uniform flow velocity value of 1 m/s is specified to the inlet of the geometry. The frontAndBack walls of the geometry are specified empty and the walls are given the noSlip condition (indicating 0 velocity) as seen before.
Initial Pressure Conditions
The initial conditions for the Pressure (p) is given below.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ 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;
}
noSlipWalls
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //The above code indicates that a uniform pressure value of 0 is specified at the outlet.
ControlDict Specification
The ControlDict dictionary is used to specify the main case controls. This includes timing information, write format, optional libraries etc. that can be loaded at runtime.
The time step value (dt) is specified according to the CFL condition, in order to obtain a stable solution.
The ControlDict file for the present simulation is given as follows.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 0.1;
deltaT 1e-5;
writeControl timeStep;
writeInterval 100;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //The ControlDict file is modified to specify the starting time as 0 and ending time as 1 sec. The time step value is given as 1e-5, to satisfy the CFL condition.
Solver Specification
The solver used to solve the given problem is the icoFoam solver. icoFoam solver solves the incompressible Navier-Stokes Equations using the PISO algorithm. The code is inherently transient and requires an initial condition and boundary conditions.
Simulation of Results
Once the mesh specifications, initial and boundary conditions are specified appropriately, the simulation setup is complete.
The geometry is meshed using the blockMesh command. icoFoam command is used to start the simulation is using the icoFoam solver. The simulation runs from the start time specified in the controlDict file to the end time (in this case, 0 to 1 sec). Once the simulation is complete, the results are written for each time step as per the criteria specified in the controlDict file, and can be post-processed to view the fluid flow.
Post-processing of Results
The results obtained after the simulation are post-processed using the ParaView post-processor. ParaView is invoked using the paraFoam command. The results for the flow velocity (U) and pressure (p) are viewed for each time step, till the steady state solution is achieved.
Mesh grading Factor = 1
A mesh grading factor of 1 implies that the starting and ending cells in a particular direction are of the same size, and that the mesh is uniform throughout the domain. Thus, there is no mesh stretching involved.
Mesh
The above figure indicates the meshed geometry with a mesh grading factor of 1.
Velocity Contour Plot
Pressure Contour Plot
Mesh grading factor = 0.8
A mesh grading factor of 0.8 implies that the size of the end cell is 0.8 times the size of the starting cell. Thus, there is mesh stretching involved here.
Mesh
The above figure indicates the mesh geometry with a mesh grading factor of 0.8.
Velocity Contour Plot
Pressure Contour Plot
Mesh grading factor = 0.5
A mesh grading factor of 0.5 implies that the size of the end cell is half the size of the starting cell. Thus, there is mesh stretching involved here.
Mesh
The above figure shows the meshed geometry with a mesh grading factor of 0.5.
Velocity Contour Plot
Pressure Contour Plot
Mesh grading factor = 0.2
A mesh grading factor of 0.2 implies that the size of ending cell is 0.2 times the size of the starting cell.
Mesh
Velocity Contour Plot
Pressure Contour Plot
Line Plots
In order to capture the effects of mesh grading near the boundary walls, a line plot is plotted at a distance of 0.085 m from the inlet, along the y axis. The line plots for different mesh grading factors are shown below.
Mesh grading factor = 1
Mesh grading factor = 0.8
Mesh grading factor = 0.5
Mesh grading factor = 0.2
The above plots show the variation of magnitude of flow velocity along the domain, at a distance of 0.085 m from the inlet. For the mesh grading factors of 1, 0.8, 0.5, the velocity plots are relatively the same, which implies negligible changes in magnitude with changing mesh grading factors. For the mesh grading factor of 0.2, the plot varies and is more defined. This can be attributed to the finer mesh produced by the mesh stretching, which leads to better capturing of the flow phenomena near the boundary walls.
A plot comparing the different velocity plots for different mesh factors is shown below.
Variation of Velocity contours near the step region
To capture the variation of magnitude of velocity near the step region in the geometry, a line plot of velocity across the domain length (0 to 0.2 m) is plotted for the different mesh grading factors.
Mesh grading factor = 1
Mesh grading factor = 0.8
Mesh grading factor = 0.5
Mesh grading factor = 0.2
- The above plots show the variation of magnitude of velocity across the domain length, for different mesh grading factors.
- In all of the above plots, the velocity is initially 0 till 0.08 m from the inlet, when there is a sudden spike in the velocity magnitude.
- This region represents the step region in the geometry, which is 0.08 m from the inlet.
- The plots for grading factors of 1, 0.8, 0.5 look relatively the same with no discernible differences.
- Whereas, the plot with mesh grading factor of 0.2 is different compared to the other three. Firstly, the magnitude of maximum velocity is higher (approximately 0.8 m/s) compared to the other three plots (0.76 m/s).
- In addition to this there is a larger dip in the 0.2 plot compared to the other plots.
- Thus, it can be stated that the mesh with a grading factor of 0.2, captures the flow phenomena near the step region in a more detailed manner when compared to the other grading factors. Hence, a finer mesh near the step region will enable better capturing of flow phenomena and results.
- In order to get further enhanced results and improve accuracy, the number of cells can be increased with this mesh grading factor of 0.2.