Week 8 - BlockMesh Drill down challenge
The objective of this project was to simulate the flow in backward facing step in openfoam using blockmesh. The geomtery of the step is as follows:
The boundary conditions are also given in the above image, with the left-most face as inlet, right-most as outlet, and no slip boundaries on top and bottom walls. To implement this geometry using blockMesh, the blocks are made as follows:
The vertices are numbered in the manner mentioned in the above image.
The general structure of blockMesh dict is 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
vertices
(
(0 0 0)
(0.08 0 0)
(0.08 0.01 0)
(0 0.01 0)
(0.2 0.01 0)
(0.2 0 0)
(0.2 -0.01 0)
(0.08 -0.01 0)
(0 0 0.00005)
(0.08 0 0.00005)
(0.08 0.01 0.00005)
(0 0.01 0.00005)
(0.2 0.01 0.00005)
(0.2 0 0.00005)
(0.2 -0.01 0.00005)
(0.08 -0.01 0.00005)
);
blocks
(
hex (0 1 2 3 8 9 10 11) (100 20 1)
simpleGrading
(
1
(
(0.5 0.5 1) //divide the block in 2 equal divisions along y,
(0.5 0.5 1) //50% cells in each division with grading of 1
)
1
)
hex (1 5 4 2 9 13 12 10) (100 20 1)
simpleGrading
(
1
(
(0.5 0.5 1)
(0.5 0.5 1)
)
1
)
hex (7 6 5 1 15 14 13 9) (100 20 1)
simpleGrading
(
1
1
1
)
);
edges
(
);
boundary
(
noslipWall
{
type wall;
faces
(
(11 10 2 3)
(10 12 4 2)
(9 8 0 1)
(14 15 7 6)
(15 9 1 7)
);
}
inlet
{
type patch;
faces
(
(8 11 0 3)
);
}
outlet
{
type patch;
faces
(
(5 4 12 13)
(6 5 13 14)
);
}
frontAndBack
{
type empty;
faces
(
(0 3 2 1)
(1 2 4 5)
(7 1 5 6)
(8 9 10 11)
(9 13 12 10)
(15 14 13 9)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
To capture the boundary layer graded mesh is used. The blocks are multi-graded with cells stretched towards top and bottom walls. The cells in x-direction in each block are 100, in y-direction 20 and in z-direction there is only 1 cell, rendering it as a 2D simulation.
The pressure boundary file is as follows:
*Note that the it is the pressure divided by density.
/*--------------------------------*- 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
{
noslipWall
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value $internalField;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //
As seen in above file the outlet is given as fixedValue with value as zero. The boundary at walls and inlet is given at zeroGradient. zeroGradient boundary condition simply extrapolates the quantity to the outlet patch from the nearest cell value. The frontAndBack are given as empty faces, as it is 2D simulation.
The U boundary file is 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 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;
}
noslipWall
{
type noSlip;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //
As seen in above file, the value of velocity is given as 1m/s in x-direction. The boundary at walls is given at noSlip which implies that the velocity at wall is zero. The outlet is given as zeroGradient. zeroGradient boundary condition simply extrapolates the quantity to the outlet patch from the nearest cell value. The frontAndBack are given as empty faces, as it is 2D simulation.
The transportProperties are 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 "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.01;
// ************************************************************************* //
In this file the kinematic viscoisty nu is given with value of 0.01 `m^2/s`
The controlDict file is 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.05;
deltaT 5e-6;
writeControl adjustableRunTime;
writeInterval 0.01;
purgeWrite 1;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable false;
// ************************************************************************* //
The solver used is icoFoam. The solver uses PISO(Pressure Implicit with Splitting Operators) algorithm to solve continuity equation `Delta*u=0` and momentum equation `(del(u))/(delt)+Delta*(uoxu)-Delta*(vDeltau)=Deltap`
The startTime is given as 0 and endTime as 0.05. The value of deltaT is chosen so that the Courant number does not blow up during the course of simulation. To check whether convergence is achieved, pressure residuals are plotted simultaneously while running simulation using gnu plot.
The fvSchemes file is 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 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;
}
// ************************************************************************* //
In this file, the various schemes to be used to discretize each term of equations are mentioned.
The fvSolution file is 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 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;
}
// ************************************************************************* //
In this file, the solver to be used to solve the linear system of equation for u and p is mentioned with the tolearnces. Also, the number of corrector loops for the PISO algorithm is specified here.
The simulation is run for grading factors of 1,0.8,0.5,0.2. The blockMesh files for different cases are as follows:
For 0.2 grading factor
/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
vertices
(
(0 0 0)
(0.08 0 0)
(0.08 0.01 0)
(0 0.01 0)
(0.2 0.01 0)
(0.2 0 0)
(0.2 -0.01 0)
(0.08 -0.01 0)
(0 0 0.00005)
(0.08 0 0.00005)
(0.08 0.01 0.00005)
(0 0.01 0.00005)
(0.2 0.01 0.00005)
(0.2 0 0.00005)
(0.2 -0.01 0.00005)
(0.08 -0.01 0.00005)
);
blocks
(
hex (0 1 2 3 8 9 10 11) (100 20 1)
simpleGrading
(
0.2
(
(0.5 0.5 5)
(0.5 0.5 0.2)
)
1
)
hex (1 5 4 2 9 13 12 10) (100 20 1)
simpleGrading
(
5
(
(0.5 0.5 5)
(0.5 0.5 0.2)
)
1
)
hex (7 6 5 1 15 14 13 9) (100 20 1)
simpleGrading
(
5
5
1
)
);
edges
(
);
boundary
(
noslipWall
{
type wall;
faces
(
(11 10 2 3)
(10 12 4 2)
(9 8 0 1)
(14 15 7 6)
(15 9 1 7)
);
}
inlet
{
type patch;
faces
(
(8 11 0 3)
);
}
outlet
{
type patch;
faces
(
(5 4 12 13)
(6 5 13 14)
);
}
frontAndBack
{
type empty;
faces
(
(0 3 2 1)
(1 2 4 5)
(7 1 5 6)
(8 9 10 11)
(9 13 12 10)
(15 14 13 9)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
For 0.5 grading factor:
/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
vertices
(
(0 0 0)
(0.08 0 0)
(0.08 0.01 0)
(0 0.01 0)
(0.2 0.01 0)
(0.2 0 0)
(0.2 -0.01 0)
(0.08 -0.01 0)
(0 0 0.00005)
(0.08 0 0.00005)
(0.08 0.01 0.00005)
(0 0.01 0.00005)
(0.2 0.01 0.00005)
(0.2 0 0.00005)
(0.2 -0.01 0.00005)
(0.08 -0.01 0.00005)
);
blocks
(
hex (0 1 2 3 8 9 10 11) (100 20 1)
simpleGrading
(
0.5
(
(0.5 0.5 2)
(0.5 0.5 0.5)
)
1
)
hex (1 5 4 2 9 13 12 10) (100 20 1)
simpleGrading
(
2
(
(0.5 0.5 2)
(0.5 0.5 0.5)
)
1
)
hex (7 6 5 1 15 14 13 9) (100 20 1)
simpleGrading
(
2
2
1
)
);
edges
(
);
boundary
(
noslipWall
{
type wall;
faces
(
(11 10 2 3)
(10 12 4 2)
(9 8 0 1)
(14 15 7 6)
(15 9 1 7)
);
}
inlet
{
type patch;
faces
(
(8 11 0 3)
);
}
outlet
{
type patch;
faces
(
(5 4 12 13)
(6 5 13 14)
);
}
frontAndBack
{
type empty;
faces
(
(0 3 2 1)
(1 2 4 5)
(7 1 5 6)
(8 9 10 11)
(9 13 12 10)
(15 14 13 9)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
For 0.8 grading factor:
/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
vertices
(
(0 0 0)
(0.08 0 0)
(0.08 0.01 0)
(0 0.01 0)
(0.2 0.01 0)
(0.2 0 0)
(0.2 -0.01 0)
(0.08 -0.01 0)
(0 0 0.00005)
(0.08 0 0.00005)
(0.08 0.01 0.00005)
(0 0.01 0.00005)
(0.2 0.01 0.00005)
(0.2 0 0.00005)
(0.2 -0.01 0.00005)
(0.08 -0.01 0.00005)
);
blocks
(
hex (0 1 2 3 8 9 10 11) (100 20 1)
simpleGrading
(
0.8
(
(0.5 0.5 1.25)
(0.5 0.5 0.8)
)
1
)
hex (1 5 4 2 9 13 12 10) (100 20 1)
simpleGrading
(
1.25
(
(0.5 0.5 1.25)
(0.5 0.5 0.8)
)
1
)
hex (7 6 5 1 15 14 13 9) (100 20 1)
simpleGrading
(
1.25
1.25
1
)
);
edges
(
);
boundary
(
noslipWall
{
type wall;
faces
(
(11 10 2 3)
(10 12 4 2)
(9 8 0 1)
(14 15 7 6)
(15 9 1 7)
);
}
inlet
{
type patch;
faces
(
(8 11 0 3)
);
}
outlet
{
type patch;
faces
(
(5 4 12 13)
(6 5 13 14)
);
}
frontAndBack
{
type empty;
faces
(
(0 3 2 1)
(1 2 4 5)
(7 1 5 6)
(8 9 10 11)
(9 13 12 10)
(15 14 13 9)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
The grading factor 1 file is as shown previously.
The mesh generated are as follows:
For grading factor 1 For grading factor 0.8
For grading factor 0.5 For grading factor 0.2
The pressure residuals as solution progresses are as follows:
For grading factor 1 For grading factor 0.8
For grading factor 0.5 For grading factor 0.2
The residuals seem to be converging after about 2000 iterations. The residuals remain constant thereafter till 10000 iterations or 0.05s which is endTime. The solution could be stopped right after the graph becomes flat though here it is made to run till endTime for consistency while comparing the profiles.
The velocity contours are as follows:
For grading factor 1 For grading factor 0.8
For grading factor 0.5 For grading factor 0.2
On closer look at the velocity profile, we can observe that the "red" region becomes smaller in width as grading factor decreases. There is also a noticeable difference in the "white" region's width as grading factor decreases, indicating that the velocity is almost zero near wall and slowly goes on increasing in case of 0.2 grading factor, while for grading factor 1, it the boundary layer is not captured properly. This is obvious due to the fact that lower the grading factor, smaller and tightly packed the cells will be towards the walls, more will be the accuracy in capturing the boundary layer.
The velocity profile at 0.085 from inlet:
For grading factor 1 For grading factor 0.8
For grading factor 0.5 For grading factor 0.2
There is not any noticeable difference in velocity profile just after step, but upon a closer look at max velocities:
| Factor=1 | Factor=0.8 | Factor=0.5 | Factor=0.2 | |
| Max Velocity | 1.10333m/s | 1.10593m/s | 1.10832m/s | 1.11035m/s |
It is seen that Grading factor 0.2 has highest velocity followed by grading factor 1. This suggest that grading factor 0.2 captures the velocity profile the best and grading factors 0.8 and 0.5 actually perform worse than factor 1.
To investigate further the streamlines are plotted with velocity:
For grading factor 1 For grading factor 0.8
For grading factor 0.5 For grading factor 0.2
As it is seen, there is a flow reversal at the bottom left corner of the step, which is best captured by grading factor 0.2 and the worse by grading factor 0.8 and 1. It reinforces the fact that decreasing grading factor better captures physics of the problem.
Courant number comparison:
| Grading 1 | Grading 0.8 | Grading 0.5 | Grading 0.2 | |
| meanCo | 0.00330228 | 0.00329862 | 0.00329088 | 0.00327635 |
| maxCo | 0.011004 | 0.0119319 | 0.0145055 | 0.0233568 |
As can be seen, the courant number decreases as grading factor increases. This is due to fact that the 'dx' decreases as grading factor is decreased, and for constant timestep 5e-6 and velocity 1m/s, causing increase in Courant number as 'dx' is in denominator. Thus we must reduce timestep when we decreasing grading factor. For present case the Courant number is well below 1 and thus timestep is not reduced. But for complicated geometries with complex flow physics, timestep will have to reduced as grading factor decreases to keep the solution stable.
In conclusion, decreasing grading factor captures the physics better, but it comes at the cost of reduced timestep. So the grading factor must be chosen optimally to get closer to real physics of the problem and also keeping in mind the computational time required to reach the convergence.