Simulation of a backward facing step in OpenFOAM
I. Aim: CFD analysis of 2D transient incompressible laminar flow in a back step type geometry using OpenFoam.
II. Problem Statement: Perform 2D transient incompressible laminar flow in a back step geometry using a suitable solver available in OpenFoam for the following cases.
Case 1 - Simulate the flow without using any grading factor for the grid (mesh) (i.e., GF = 1)
Case 2 - Simulate the flow with a mesh grading factor of 0.2. The cells should be finer near the walls (including the step wall).
Measure the velocity profile at 0.085 m from the inlet of the geometry for both case studies and compare.
III. Geometry:
The geometry domain dimensions and the boundary details are given in the below figure. All dimensions are in 'm'.
IV . Simulation setup in OpenFoam:
1. Mesh Generation:
Since the problem is of transient incompressible laminar flow, a similar case available in the tutorials section of OpenFOAM by the name 'cavity' is used as a reference file for the analysis. The path to the same is:
tutorials/incompressible/icoFoam/cavity/cavity
This case uses the 'icofoam' solver, which is an incompressible laminar flow solver.
The reference case folder 'cavity', is copied to the 'run' folder and renamed as backstep. The backstep folder consists of the initial conditions (0), constant and system folders.
The system folder consists of the blockMeshDict, controlDict, fvSchemes, fvSolution files.
The blockMeshDict contains the information about the mesh. The file is edited to make the geometry of the backstep using vertices and blocks. Also, the mesh is created using hexahedral elements with grading factors under consideration GF = 1 and 0.2. The boundary face and their type are alos defined in this file.
The boundary types are as below:
1. Inlet: Type - Patch
2. Outlet: Type - Patch
3. Walls(including sidestep): Type noSlip
4. Front &Back faces: Type - empty. Since the problem is 2D in nature, the front and back faces of the blocks are defined as empty
The blockMeshDict files for GF=1 and GF=0.2 cases are given below:
blockMeshDict file for GF=1:
/*--------------------------------*- 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.1;
vertices
(
(0 0 0)
(0.08 0 0)
(0.08 0.005 0)
(0 0.005 0)
(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)
(0 0.01 0)
(0 0.01 0.1)
(0.08 0.01 0.1)
(0.2 0 0)
(0.2 0.005 0)
(0.2 0.01 0)
(0.2 0 0.1)
(0.2 0.005 0.1)
(0.2 0.01 0.1)
(0.08 -0.01 0)
(0.2 -0.01 0)
(0.08 -0.01 0.1)
(0.2 -0.01 0.1)
);
blocks
(
hex (0 1 2 3 4 5 6 7) (81 6 1) simpleGrading (1 1 1)
hex (3 2 8 9 7 6 11 10) (81 6 1) simpleGrading (1 1 1)
hex (1 12 13 2 5 15 16 6) (121 6 1) simpleGrading (1 1 1)
hex (2 13 14 8 6 16 17 11) (121 6 1) simpleGrading (1 1 1)
hex (18 19 12 1 20 21 15 5) (121 11 1) simpleGrading (1 1 1)
);
edges
(
);
boundary
(
inlet
{
type patch;
faces
(
(0 4 7 3)
(3 7 10 9)
);
}
outlet
{
type patch;
faces
(
(21 15 12 19)
(15 16 13 12)
(16 17 14 13)
);
}
front&back
{
type empty;
faces
(
(4 5 6 7)
(7 6 11 10)
(20 21 15 5)
(5 15 16 6)
(6 16 17 11)
(1 0 3 2)
(12 1 2 13)
(19 18 1 12)
(2 3 9 8)
(13 2 8 14)
);
}
noslipwall
{
type noSlip;
faces
(
(9 10 11 8)
(8 11 17 14)
(18 19 21 20)
(0 1 5 4)
(18 1 5 20)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
BlockMeshDict file for GF=1:
/*--------------------------------*- 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.1;
vertices
(
(0 0 0)
(0.8 0 0)
(0.8 0.05 0)
(0 0.05 0)
(0 0 1)
(0.8 0 1)
(0.8 0.05 1)
(0 0.05 1)
(0.8 0.1 0)
(0 0.1 0)
(0 0.1 1)
(0.8 0.1 1)
(2 0 0)
(2 0.05 0)
(2 0.1 0)
(2 0 1)
(2 0.05 1)
(2 0.1 1)
(0.8 -0.1 0)
(2 -0.1 0)
(0.8 -0.1 1)
(2 -0.1 1)
);
blocks
(
hex (0 1 2 3 4 5 6 7) (81 6 1) simpleGrading (0.2 20 1)
hex (3 2 8 9 7 6 11 10) (81 6 1) simpleGrading (0.2 0.2 1)
hex (1 12 13 2 5 15 16 6) (121 6 1) simpleGrading (20 20 1)
hex (2 13 14 8 6 16 17 11) (121 6 1) simpleGrading (20 0.2 1)
hex (18 19 12 1 20 21 15 5) (121 11 1) simpleGrading (20 20 1)
);
edges
(
);
boundary
(
inlet
{
type patch;
faces
(
(0 4 7 3)
(3 7 10 9)
);
}
outlet
{
type patch;
faces
(
(21 15 12 19)
(15 16 13 12)
(16 17 14 13)
);
}
front&back
{
type empty;
faces
(
(4 5 6 7)
(7 6 11 10)
(20 21 15 5)
(5 15 16 6)
(6 16 17 11)
(1 0 3 2)
(12 1 2 13)
(19 18 1 12)
(2 3 9 8)
(13 2 8 14)
);
}
noslipwall
{
type noSlip;
faces
(
(9 10 11 8)
(8 11 17 14)
(18 19 21 20)
(0 1 5 4)
(18 1 5 20)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
The above blockMeshDict file is executed using 'blockMesh' command to generate the mesh and boundary types. The geometry with mesh can be viewed in Paraview by executing the 'paraFoam' command and are shown below:
Mesh with GF =1:
The grading factor represents the ratio of last cell size to first cell size. The above figure shows the mesh is uniform (GF = 1) throughout the domain.
Mesh with GF =0.2:
From the above figure, the mesh is graded towards the walls (GF = 0.2), including near the step wall. The cell size increases from the interior to the walls due to the grading factor.
2. Initial and Boundary conditions:
The initial and boundary conditions are next defined by editing the boundary files 'p' and 'U' in the initial conditions folder (0).
Inlet:
Velocity - uniform - 2 m/s
Pressure - zero gradient
Outlet:
Pressure - uniform - 0 Pa.
velocity - zero gradient
Walls:
velocity - noslip
Pressure - zero gradient
Front & Back
velocity - empty
Pressure - empty
The corresponding script files are shown below:
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;
}
frontAndBack
{
type empty;
}
noslipwall
{
type zeroGradient;
}
}
// ************************************************************************* //Velcoity:
/*--------------------------------*- 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 (2 0 0);
}
outlet
{
type zeroGradient;
}
noslipwall
{
type noSlip;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //3. Solver setup:
The solution setup consists of the following steps:
1. Defining the transport properties:
Since it is an incompressible laminar flow problem, the only equation solved is the momentum equation. For this, we have to define the transport property kinematic viscosity. For the present case, a fluid with kinematic viscosity (`nu`= 0.01 m^2/s) is assumed. The transport properties are defined in script file in the folder 'constant' and is given below:
/*--------------------------------*- 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 "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.01;
// ************************************************************************* //2. Solution type and duration:
The solver 'icofoam', which is a transient incompressible laminar solver is used for solving the present case. The solver conditions are defined in the 'controlDict' file in the 'system' folder. The total time duration for the analysis is assumed to be 0.5 sec. The initial time, end time, timestep size are defined in this file. The time step decides the Courant number and stability of the solution and for a convective case it is given as:
`C = U.dt/dx`
where, C = Courant number ;
U = maximum velocity in the flow feild
dt = timetep size
dx = smallest mesh cell size
The courant number should be less than or equal to one to have a stable solution. The Courant number is assumed to be 0.01 to calculate the timestep for the solution to be stable.
The calculated time dt for C =0.01 and U = 2 m/s are:
1. Mesh with GF=1 , dx = 0.001, dt = 5e-6;
2. Mesh with GF=0.2 , dx = 0.0002, dt = 1e-6;
The controlDict files for the two cases are given below:
GF =1:
/*--------------------------------*- 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.5;
deltaT 0.000005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //GF =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;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 0.5;
deltaT 0.000001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //The solver scheme and solution script files present in the system folder used by the 'icoFoam' are given below:
fvSchemes
/*--------------------------------*- 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 fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linear;
}
laplacianSchemes
{
default Gauss linear orthogonal;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default orthogonal;
}
// ************************************************************************* //fvSolution
/*--------------------------------*- 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 fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0.05;
}
pFinal
{
$p;
relTol 0;
}
U
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-05;
relTol 0;
}
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
// ************************************************************************* //After defining the initial time, end time and time step size in the 'controlDict' file, the 'icoFoam' command is executed to run the analysis. After the analysis is carried upto the end time iteratively, 'paraFoam' command is used to post-process the analysis results in Paraview.
4: Post-Processing:
The below plots show the velocity and pressure variation in the domain with respect to time :
1. Velocity:
| GF =1: | GF =0.2 |
| At t = 0 sec | At t = 0 sec |
 |
 |
| At t = 0.1 sec | At t = 0.1 sec |
 |
 |
| At t =0.5 sec | At t =0.5 sec |
 |
 |
The velocity at the inlet is 2 m/s which correspond to the Reynold's number, Re = 6.66 which indicates the flow is laminar. Due to the laminar condition with no-slip at the walls, the velocity profile variation is smooth throughout the cross-section with zero at the walls and maximum at the centreline. Also, the velocity reduced from2 m/s to a lesser value around 1.5 m/s after the back step cross-section because the cross-sectional area is increased leading to a reduction in velocity satisfying continuity equation. The above contour plots are for the initial time, halfway through the time marching and end time. Similar, profiles are observed for the GF = 0.2 case also.
2. Pressure:
| GF =1: | GF =0.2 |
| At t = 0 sec | At t = 0 sec |
 |
|
| At t = 0.1 sec | At t = 0.1 sec |
 |
|
| At t =0.5 sec | At t =0.5 sec |
 |
|
3. Velocity variation at x =0.085 m from inlet:
| GF =1 | GF =0.2 |
 |
 |
The above plots compare the velocity variation at x =0.085 m from inlet, for both GF= 1 and GF =0.2 case. The velocity variation is smoother in the case of GF =1 when compared to the FG=0.2 case. This is because the cells are uniform in the GF=1 case and in case of GF=0.2, there were finer cells at the centerline and due to a gradual increase in cell size away from the centerline, sharp variations in velocity are observed.
Conclusion:
The velocity profiles in both cases are almost similar although different meshes are used. This is because the time duration considered for the analysis and the input velocity is very small. The effect of the grading factor on the results is negligible for the present problem. However, the grading factor influences the gradients and results values at the boundaries. Fine cells should be used at the walls to accurately measure the gradients at the walls.