Simulation of a backward facing step in OpenFOAM
In this project, I will be simulating an incompressible-laminar-viscous flow through a backward facing step. This will be solved using the software OpenFOAM running on Debian.
The solver for this simulation was chosen to be icoFoam, I decided what solver to use from this guide on OpenFOAM's website
This is from their guide, for icoFoam:
We can notice that in our case, the flow is laminar and incompressible and the solution is transient as well.
I will be solving for 2 cases, one without Grading Factor and the other with a Grading Factor.
Grading Factor and its importance in Meshing
When solving a problem numerically, it is essential to do meshing. When less elements are generated in the mesh, we trade off some accuracy of the true solution for faster runtimes and more stability. A Graded mesh essentially reduces the size of the elements in the mesh in one part of the domain while increasing the size in another part. If done correctly, a graded mesh can improve the accuracy of a solution. A graded mesh allows us to increase the accuracy of the solution in the part of the domain we are interested in (near the step walls, in this case), while sacrificing some accuracy in other parts we are not interested in. The advantage is that, since the number of elements remains constant, the stability (courant number) of the solution is not affected.
Case 1: No Grading Factor
First I start off by navigating to the tutorials folder in OpenFOAM, then I find the icoFoam solver under incompressible flow. I decide to use the cavity solver, and I proceed to copy the cavity solver from the tutorials folder to the run folder under the name week8_case1
95b82f586eb4: ~>> tut
95b82f586eb4: tutorials>> cd incompressible/
95b82f586eb4: incompressible>> cd icoFoam/
95b82f586eb4: icoFoam>> cd cavity/
95b82f586eb4: cavity>> ls
Allclean Allrun cavity cavityClipped cavityGrade
95b82f586eb4: cavity>> cp -r cavity/ $FOAM_RUN/week8_case1
95b82f586eb4: cavity>> run
95b82f586eb4: run>> cd week8_case1/
95b82f586eb4: week8_case1>> ls
0 constant system
We can notice that there are 3 folders inside. The 0 folder is for the initial conditions, constant folder contains BlockMeshDict, fvSchemes, fvSolutions and controlDict. The system folder contains transportProperties.
These are the files that I changed from their default contents which came from the tutorial.
system Folder
blockMeshDict - This file contains the geometry of the domain, the faces, edges, etc. It also contains the mesh properties. The mesh used in this case is a uniformly sized mesh containing 200 elements in the longest X direction and 20 in the longest Y direction.
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 9
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
// vertices from point 0 to 10
(0 0 0)
(0.08 0 0)
(0.08 0.005 0)
(0 0.005 0)
(0.08 0.01 0)
(0 0.01 0)
(0.2 0.01 0)
(0.2 0.005 0)
(0.2 0 0)
(0.2 -0.01 0)
(0.08 -0.01 0)
// vertices from point 11 to 21
(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.1)
(0 0.01 0.1)
(0.2 0.01 0.1)
(0.2 0.005 0.1)
(0.2 0 0.1)
(0.2 -0.01 0.1)
(0.08 -0.01 0.1)
);
blocks
(
// block 1
hex (0 1 2 3 11 12 13 14) (80 5 1) simpleGrading (1 1 1)
//block 2
hex (3 2 4 5 14 13 15 16) (80 5 1) simpleGrading (1 1 1)
//block 3
hex (2 7 6 4 13 18 17 15) (120 5 1) simpleGrading (1 1 1)
//block 4
hex (1 8 7 2 12 19 18 13) (120 5 1) simpleGrading (1 1 1)
//block 5
hex (10 9 8 1 21 20 19 12) (120 10 1) simpleGrading (1 1 1)
);
edges
(
);
boundary
(
inlet
{
type patch;
faces
(
(11 14 3 0)
(14 16 5 3)
);
}
outlet
{
type patch;
faces
(
(6 17 18 7)
(8 7 18 19)
(9 8 19 20)
);
}
otherwall //no slip walls
{ type wall;
faces
(
//top
(15 4 5 16)
(15 17 6 4)
//step
(21 12 1 10)
//bottom left
(11 0 1 12)
//bottom right
(21 10 9 20)
);
}
frontAndBack
{
type empty;
faces
( //back
(0 3 2 1)
(3 5 4 2)
(2 4 6 7)
(1 2 7 8)
(10 1 8 9)
//front
(13 15 16 14)
(12 13 14 11)
(20 19 12 21)
(19 18 13 12)
(18 17 15 13)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //controlDict - This folder contains some essential values which controls the simulation (like timeperiod, timestep, precision, etc.) In this case, the simulation was run for 5 seconds (5 cycles in total) and with a timestep of 0.001s.
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 9
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 5;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
//adaptive timestep
/*
adjustTimeStep true;
maxCo 1;*/
// ************************************************************************* //folder fvSolutions and fvSchemes were not changed for this simulation
0 folder
U (Velocity) - This contains the initial and boundary velocity conditions for the domain. The initial conditions are 0.1 m/s across the domain. The boundary condition at the inlet is 0.2 m/s.
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 9
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.1 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.2 0 0);
}
outlet
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
otherwall
{
type noSlip; // zero at walls
}
}
// ************************************************************************* //p (pressure) - This contains the initial and boundary pressure conditions for the domain. The initial conditions are 0 Pa across the domain. The boundary condition at the outlet is 0 Pa.
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 9
\/ 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;
}
frontAndBack
{
type empty;
}
otherwall
{
type zeroGradient;
}
}
// ************************************************************************* //constant Folder
Upon running the blockMesh command, that script also generates a polyMesh folder in this folder, which contains details of the mesh.
transportProperties - This contains other values required to run the simulation but remain constant through the solution (such as Dynamic Viscosity `nu`)
Dynamic viscosity was chosen to be 1e-5.
transportProperties
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 9
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 1e-5;
// ************************************************************************* //Now, I begin the simulation process.
1. Running the blockMesh command
95b82f586eb4: week8_case1>> blockMesh
/*---------------------------------------------------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 9
\/ M anipulation |
*---------------------------------------------------------------------------*/
Build : 9-a0f1846504f2
Exec : blockMesh
Date : Nov 11 2021
Time : 09:13:19
Host : "95b82f586eb4"
PID : 678
I/O : uncollated
Case : /home/openfoam/run/week8_case1
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time
Reading "blockMeshDict"
Creating block mesh from
"system/blockMeshDict"
Creating block edges
No non-planar block faces defined
Creating topology blocks
Creating topology patches
Creating block mesh topology
Check topology
Basic statistics
Number of internal faces : 5
Number of boundary faces : 20
Number of defined boundary faces : 20
Number of undefined boundary faces : 0
Checking patch -> block consistency
Creating block offsets
Creating merge list .
Creating polyMesh from blockMesh
Creating patches
Creating cells
Creating points with scale 1
Block 0 cell size :
i : 0.001 .. 0.001
j : 0.001
k : 0.1
Block 1 cell size :
i : 0.001 .. 0.001
j : 0.001 .. 0.001
k : 0.1
Block 2 cell size :
i : 0.001 .. 0.001
j : 0.001 .. 0.001
k : 0.1
Block 3 cell size :
i : 0.001 .. 0.001
j : 0.001
k : 0.1
Block 4 cell size :
i : 0.001 .. 0.001
j : 0.001 .. 0.001
k : 0.1
Writing polyMesh
----------------
Mesh Information
----------------
boundingBox: (0 -0.01 0) (0.2 0.01 0.1)
nPoints: 6842
nCells: 3200
nFaces: 13020
nInternalFaces: 6180
----------------
Patches
----------------
patch 0 (start: 6180 size: 10) name: inlet
patch 1 (start: 6190 size: 20) name: outlet
patch 2 (start: 6210 size: 410) name: otherwall
patch 3 (start: 6620 size: 6400) name: frontAndBack
End
Generated Mesh
The following commands are executed
icoFoam
paraFoam
This means the solution is now ready, paraFoam opens the paraView window which contains our solution
This is the solution in gif format
direct gif link: https://gfycat.com/gratefulalarmedindianpangolin
direct youtube link: https://youtu.be/pPJlCernYNk
Observations: We can notice that the flow through the inlet section is laminar until the step wall is encountered. The velocity of the flow in the inlet region is high at the centre and progressively decreases going towards the walls. When the flow encounters the step wall for the first time, there are minor vortices formed near the step wall and at the top of the wall which dissipates as the flow approaches steady state. Slowly the flow approaches steady state, at around 125th timestep (2.5s) into the simulation. There is no significant change in flow after it attains steady state.
Case 2 (Grading factor of 0.2) :
First I start off by navigating to the tutorials folder in OpenFOAM, then I find the icoFoam solver under incompressible flow. I decide to use the cavity solver, and I proceed to copy the cavity solver from the tutorials folder to the run folder under the name week8_case2
95b82f586eb4: ~>> tut
95b82f586eb4: tutorials>> cd incompressible/
95b82f586eb4: incompressible>> cd icoFoam/
95b82f586eb4: icoFoam>> cd cavity/
95b82f586eb4: cavity>> ls
Allclean Allrun cavity cavityClipped cavityGrade
95b82f586eb4: cavity>> cp -r cavity/ $FOAM_RUN/week8_case2
95b82f586eb4: cavity>> run
95b82f586eb4: run>> cd week8_case2/
95b82f586eb4: week8_case2>> ls
0 constant system
We can notice that there are 3 folders inside. The 0 folder is for the initial conditions, constant folder contains BlockMeshDict, fvSchemes, fvSolutions and controlDict. The system folder contains transportProperties.
These are the files that I changed from their default contents which came from the tutorial.
system Folder
blockMeshDict - This file contains the geometry of the domain, the faces, edges, etc. It also contains the mesh properties. The mesh used in this case is a graded mesh containing 200 elements in the longest X direction and 20 in the longest Y direction. The grading factor is 0.2 and in such a way that elements near the step wall is smaller than elements on the other side of the domain.
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 9
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
// vertices from point 0 to 10
(0 0 0)
(0.08 0 0)
(0.08 0.005 0)
(0 0.005 0)
(0.08 0.01 0)
(0 0.01 0)
(0.2 0.01 0)
(0.2 0.005 0)
(0.2 0 0)
(0.2 -0.01 0)
(0.08 -0.01 0)
// vertices from point 11 to 21
(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.1)
(0 0.01 0.1)
(0.2 0.01 0.1)
(0.2 0.005 0.1)
(0.2 0 0.1)
(0.2 -0.01 0.1)
(0.08 -0.01 0.1)
);
blocks
(
// block 1
hex (0 1 2 3 11 12 13 14) (80 5 1) simpleGrading (0.2 5 1)
//block 2
hex (3 2 4 5 14 13 15 16) (80 5 1) simpleGrading (0.2 0.2 1)
//block 3
hex (2 7 6 4 13 18 17 15) (120 5 1) simpleGrading (5 0.2 1)
//block 4
hex (1 8 7 2 12 19 18 13) (120 5 1) simpleGrading (5 5 1)
//block 5
hex (10 9 8 1 21 20 19 12) (120 10 1) simpleGrading (5 1 1)
);
edges
(
);
boundary
(
inlet
{
type patch;
faces
(
(11 14 3 0)
(14 16 5 3)
);
}
outlet
{
type patch;
faces
(
(6 17 18 7)
(8 7 18 19)
(9 8 19 20)
);
}
otherwall //no slip walls
{ type wall;
faces
(
//top
(15 4 5 16)
(15 17 6 4)
//step
(21 12 1 10)
//bottom left
(11 0 1 12)
//bottom right
(21 10 9 20)
);
}
frontAndBack
{
type empty;
faces
( //back
(0 3 2 1)
(3 5 4 2)
(2 4 6 7)
(1 2 7 8)
(10 1 8 9)
//front
(13 15 16 14)
(12 13 14 11)
(20 19 12 21)
(19 18 13 12)
(18 17 15 13)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //controlDict - This folder contains some essential values which controls the simulation (like timeperiod, timestep, precision, etc.) In this case, the simulation was run for 5 seconds (5 cycles in total) and with a timestep of 0.001s.
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 9
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 5;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
//adaptive timestep
/*
adjustTimeStep true;
maxCo 1;*/
// ************************************************************************* //folder fvSolutions and fvSchemes were not changed for this simulation
0 folder
U (Velocity) - This contains the initial and boundary velocity conditions for the domain. The initial conditions are 0.1 m/s across the domain. The boundary condition at the inlet is 0.2 m/s.
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 9
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.1 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.2 0 0);
}
outlet
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
otherwall
{
type noSlip; // zero at walls
}
}
// ************************************************************************* //p (pressure) - This contains the initial and boundary pressure conditions for the domain. The initial conditions are 0 Pa across the domain. The boundary condition at the outlet is 0 Pa.
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 9
\/ 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;
}
frontAndBack
{
type empty;
}
otherwall
{
type zeroGradient;
}
}
// ************************************************************************* // constant Folder
Upon running the blockMesh command, that script also generates a polyMesh folder in this folder, which contains details of the mesh.
transportProperties - This contains other values required to run the simulation but remain constant through the solution (such as Dynamic Viscosity `nu`)
Dynamic viscosity was chosen to be 1e-5.
transportProperties
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 9
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 1e-5;
// ************************************************************************* //Now, I begin the simulation process.
1. Running command blockMesh
54ae5cb21649: week8_case2>> blockMesh
/*---------------------------------------------------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 9
\/ M anipulation |
*---------------------------------------------------------------------------*/
Build : 9-a0f1846504f2
Exec : blockMesh
Date : Nov 12 2021
Time : 08:52:33
Host : "54ae5cb21649"
PID : 152
I/O : uncollated
Case : /home/openfoam/run/week8_case2
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time
Reading "blockMeshDict"
Creating block mesh from
"system/blockMeshDict"
Creating block edges
No non-planar block faces defined
Creating topology blocks
Creating topology patches
Creating block mesh topology
Check topology
Basic statistics
Number of internal faces : 5
Number of boundary faces : 20
Number of defined boundary faces : 20
Number of undefined boundary faces : 0
Checking patch -> block consistency
Creating block offsets
Creating merge list .
Creating polyMesh from blockMesh
Creating patches
Creating cells
Creating points with scale 1
Block 0 cell size :
i : 0.00200653 .. 0.000401307
j : 0.001
k : 0.1
Block 1 cell size :
i : 0.00200653 .. 0.000401307
j : 0.001 .. 0.001
k : 0.1
Block 2 cell size :
i : 0.00040166 .. 0.0020083
j : 0.001 .. 0.001
k : 0.1
Block 3 cell size :
i : 0.00040166 .. 0.0020083
j : 0.001
k : 0.1
Block 4 cell size :
i : 0.00040166 .. 0.0020083
j : 0.001 .. 0.001
k : 0.1
Writing polyMesh
----------------
Mesh Information
----------------
boundingBox: (0 -0.01 0) (0.2 0.01 0.1)
nPoints: 6842
nCells: 3200
nFaces: 13020
nInternalFaces: 6180
----------------
Patches
----------------
patch 0 (start: 6180 size: 10) name: inlet
patch 1 (start: 6190 size: 20) name: outlet
patch 2 (start: 6210 size: 410) name: otherwall
patch 3 (start: 6620 size: 6400) name: frontAndBack
End
Generated Mesh
The following commands are executed
icoFoam
paraFoam
This means the solution is now ready, paraFoam opens the paraView window which contains our solution
This is the solution in gif format
Direct gif link: https://gfycat.com/downrightsecretchinesecrocodilelizard
Direct youtube link: https://youtu.be/RqeeqkB5O-0
view at 4x speed for optimal experience
Observations: We can notice that the flow through the inlet section is laminar until the step wall is encountered. The velocity of the flow in the inlet region is high at the centre and progressively decreases going towards the walls. When the flow encounters the step wall for the first time, there are minor vortices formed near the step wall and at the top of the wall which dissipates as the flow approaches steady state. Slowly the flow approaches steady state, at around 125th timestep (2.5s) into the simulation. There is no significant change in flow after it attains steady state.
Velocity Profile at x= 0.085 for both Cases
Plot over line tool was used to generate a csv file for both cases, the file contains the distribution of velocity for both cases at x=0.085m. A simple plot was generated from those values
Case 1
Case 2
Observations: We can notice that case 1's curve is quite smoother than case 2's curve at all points. Additionally, the maximum value for Case 1: 0.285561 is slightly higher than that of Case 2: 0.279766
- Sarath Kumar Uhi brother my name sarath and have a problem need ur support could please contact my watzup 7902298078