Week 8 - BlockMesh Drill down challenge
In this challenge CFD simulation will be performed using OpenFoam. IcoFoam is a solver which solves incompressible navier stokes equation for a flow in consideration. Here a stepped gap between plates will be used as a subject in which flow will be simulated. While performing the simulation scaling factors will be tested to test the changes in the results which are namely; 0.2,0.5 and 0.8.
These scaling factors will be used to make mesh finer at the step because that is the point where a notable change may occur in the flow. First it has to be insured that the flow that develops between the plates must be laminar flow. Because only when the flow becomes laminar, then it can be concluded that the solution is stable and ripe for analysis.
Assume that fluid flows between plates separated by a width. Now since the flow is two-dimensional, it can be assumed that the plates have infinite breadth. But length and the distance between the plates is known, and they are w and L respectively.
Now before simulation through stepped plates is done, it is of paramount importance to find a pressure difference to be applied across the subject to render a laminar flow across the subject in consideration. The pressure difference induces fluid to move and take up a velocity but how laminar flow works is that when pressure induced velocity starts to increase, it settles and sets at a steady state velocity and engenders a steady state velocity profile across the cross section. But it may happen that pressure difference may be at a value which may cause the fluid to flow at a velocity to transcend beyond the operational ceiling required for the flow to be laminar. And hence a formula will be derived for a pressure difference reqiured for laminar flow across a two parallel plates of infinite breadth and finite length and this formula will in turn will be used to decide and calculate the pressure difference to be applied across the step.
Consider a fluid of density `rho` and viscosity `mu`
Let `u_x and u_y` be velocity of fluid at any point in the direction x and y respectively.
Since incompressible navier stokes equation is in equation, continuity and momentum equation will be used to derive the formula.
Let x be the direction of flow and y be the direction transverse to the flow direction.
1)Continuity equation
`(del(rho.u_x))/(delx)+(del(rho.u_y))/(dely)=0`
but since we are considering the fluid to be incompressible, therefore the density is constant everywhere. Hence the G.D.E becomes
`(del(u_x))/(delx)+(del(u_y))/(dely)=0`
now it is known that `u_x and u_y` are function of both x and y.
Now since this governing differential equation for steady state and pressure is applied across x-direction and therefore the fluid will not have any impetus to flow in y-direction and since diffusion effects are neglected, this will imply that
`u_y=0` at all points
this means that
`(del(u_y))/(dely)=0`
when substituted in incompressible continuity equation we get,
`(del(u_x))/(delx)=0`
this means that u_x is constant in x-direction.
2)Momentum equation
`rho((del(u_x))/(delt)+u_x(del(u_x))/(delx)+u_y(del(u_x))/(dely))=-(del(P))/(delx)+mu((del^2(u_x))/(delx^2)+(del^2(u_x))/(dely^2))`
But since steady state is required
`(del(u_x))/(delt)=0`
`(del(u_x))/(delx)=0`....proven by continuity equation
`u_y=0`....as indicated previously
And therefore the momentum equation becomes
`(del(P))/(delx)=mu((del^2(u_x))/(delx^2)+(del^2(u_x))/(dely^2))`
but since `(del(u_x))/(delx)=0` therefore `(del^2(u_x))/(delx^2)=0`
So the equation becomes
`(del(P))/(delx)=mu(del^2(u_x))/(dely^2)`
Consider that pressure gradient is constant across the x-direction,when pressure difference is applied across the subjects length
So the pressure gradient term can be replaced by a constant term
`(del(P))/(delx)=k`
And the equation becomes
`k=mu(del^2(u_x))/(dely^2)`
After integrating above G.D.E equation twices with respect to y, we get
`mu.u_x=(k.y^2)/2+C_1.y+C_2`
now by applying boundary condition
1)`u_x=0` when y=0
we get
`C_2=0`
2)`u_x=0` when y=w (width)
we get
`C_1=(k.w)/2`
And therefore the equation becomes
`u_x=(k.y)/(2.mu) . (y-w)`
Now to find where u_x is maximum, at that point `(du_x)/(dy)=0`
differentiating both the sides,
we get
`(du_x)/(dy)=k/(2.mu) .(2y-w)`
therefore
`k/(2.mu) .(2y-w)=0`
simplifying we get
`y=w/2`
substituting this value in the original equation we get
`u_(xmax)=-(kw^2)/(8mu)`
but it seems that the term means that it is minima instead of maxima, but k which is actually pressure gradient is inherently negative since it is highest at x=0 and lowest at x=L.
So `k=dp/dx=-(deltap)/L`
substituting in the equation, we get
`u_(xmax)=(deltap.w^2)/(8.mu.L)`
it is important to look at volumetric flow rate per unit breadth of plate,
`q=int_0^w u_x.dy`
`q=int_0^w k/(2mu)(y^2-wy).dy`
`q=k/(2.mu)[y^3/(3)-(w.y^2)/2]_0^w`
`q=-(kw^3)/(12.mu)`
substituting `k=-(deltap)/L`
`q=(deltap.w^3)/(12.mu.L)`
but since `u_(xmax)=-(deltap.w^2)/(8.mu.L)`
and further rearranging these terms
`q=u_(xmax)*w*(2/3)`
but in terms of `V_(avg)` (average velocity)
`q=V_(avg)*w`
comparing both the equations, we get
`V_(avg)=2/3*u_(xmax)`
so the formula for `V_(avg)` becomes
`V_(avg)=(deltap.w^2)/(12.mu.L)`
Now its important to establish a relation between Reynold number and pressure difference
`Re=(rho.w.V_(avg))/(mu)`
by rearranging we get
`V_(avg)=(mu.Re)/(rho.w)`
equating both equations for `V_(avg)`
`(mu.Re)/(rho.w)=(deltap.w^2)/(12.mu.L)`
rearranging both sides to pressure differencing
`dp=(12.mu^2Re.L)/(rho.w^3)`
equating V_avg formula of reynolds number and u_(xmax)
`2/3*u_(xmax)=(mu.Re)/(rho.w)`
`u_(xmax)=(3.mu.Re)/(2.rho.w)`
As it can be seen that pressure difference has been made a function of Reynolds number, which inturn has indirectly engendered `u_(xmax)` to be a function of Reynolds number
Since reynolds number determines the type of flow in the fluid, so a appropriate Re will determine the pressure difference value and the velocity maximum for laminar flow.
Now since formula for dp and `u_(xmax)` has been derived, it is important decide on the material that has to be taken into consideration for the value of density and viscosity.
Material selected=Water
viscosity=`mu=1.002e-3`Pa-s
density=`rho=998`kg/m^3
this makes
kinematice viscosity=viscosity/density
`nu=mu/rho`
`nu=1.004e-6`m^2/s
length=L=0.2m from the length of the stepped plates of infinite breadth
the width size has to be ascertained between 0.01m or 0.02m because of the two width sizes due to step.
Now from the formula of dp,
`dp prop (1/w^3)`
smaller pressure difference for whichever width has to be chosen as the width to calculate the pressure difference.
since dp is smaller for 0.02, therefore width =0.02m will be taken as pressure difference for the calculation.
By hit and trial method, it has been found that a steady state laminar velocity profile is acheived for Re=100 or less. So to be on the safer side
Re=50
now knowing what `V_max` it can be decided what value of dt must be for a stable simulation using formula for C.F.L number,
`C=(V_(max)*dt)/dx`
So dt becomes
`dt=(C.dx)/(V_max)`
And the CFL number decided for the safety of the simulation is 0.01
`C=0.01`
Now to avoid human error to creep in, a code was written in octave to facilitate calculation of all the required terms
Since pressure which is available as an input is in the form of kinematice pressure, the formula for kninematic pressure which is used in the code is
`dp_k=(dp)/rho`
`dp=(12.mu^2Re.L)/(rho^2.w^3)`
clear all
close all
clc
viscosity=1.002*10^-3;
density=998;
nu=viscosity/density
L=0.2;
Re=50;
w=0.02;
dp=(12*(viscosity^2)*L*Re)/(density^2*(w^3))
vmax=(3*viscosity*Re)/(2*density*w)
CFL=0.01;
dx=0.001;
dt=(CFL)*dx/vmax;
dtnew=0.001;%selected dt
newCFL=vmax*dtnew/dx%based on selected dtIn this code, first maximum dt is found out using the formula and based on that maximum, a convenient value of dt is chosen and new CFL number is calculated to be compared with the number that appear in the simulation calculation in the command prompt.
The results in the command window are
Now first all the values obtained from these formulas will be tested by running a simulation for a simple fluid flow between two infinitely wide plates of length 0.2m and width size 0.02m.
blockMeshDict file:-
/*--------------------------------*- 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 0.01;
vertices
(
(0 0 0)
(20 0 0)
(20 2 0)
(0 2 0)
(0 0 1)
(20 0 1)
(20 2 1)
(0 2 1)
);
blocks
(
hex (0 1 2 3 4 5 6 7) (200 10 1) simpleGrading (1 1 1)
);
edges
(
);
boundary
(
inletWall
{
type patch;
faces
(
(0 4 7 3)
);
}
outletWall
{
type patch;
faces
(
(1 2 6 5)
);
}
fixedWalls
{
type wall;
faces
(
(2 3 7 6)
(0 1 5 4)
);
}
frontAndBack
{
type empty;
faces
(
(4 5 6 7)
(0 3 2 1)
);
}
);
mergePatchPairs
(
);
Initial Condition:
1)pressure
/*--------------------------------*- 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
{
inletWall
{
type fixedValue;
value uniform 0.00001512;
}
outletWall
{
type fixedValue;
value uniform 0;
}
fixedWalls
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
2) U
/*--------------------------------*- 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
{
inletWall
{
type zeroGradient;
}
outletWall
{
type zeroGradient;
}
fixedWalls
{
type noSlip;
}
frontAndBack
{
type empty;
}
}
Now the properties regarding time step size will be specified in controlDict file
/*--------------------------------*- 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 400;
deltaT 0.001;
writeControl timeStep;
writeInterval 100;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
All the files which were modified in the cavity tutorial of icoFoam tutorial has been exhibited here. Rest all files specifying the kind of solver(fvScheme) and tolerance(fvSolution) has been kept the same.
To run the simulation,first blockMesh command is given which generates mesh in the geometry specified in the blockMeshDict file. The mesh is then viewed in Paraview using ParaFoam command:
Now icoFoam solver will be used to solve the problem. Typing icoFoam command runs the simulation for stipulated time.The file was run for 400s to reach the steady state and verify whether such a state can be reached.
At the calculation of 400th s in the command prompt, maximum CFL number is calculated. This calculated CFL number from simulation will be compared with the CFL number calculated from theoretical maximum velocity.
From the simulation results it is found that CFL number repeated in its value and didn't increase for a number of time steps before 400s which indicates that the system has reached steady state.
Maximum CFL number from simulation=0.00376473
Maximum theoretical CFL number=0.003765
It appears that CFL numbers for both simulation and theoretical calculation match approximately with each other.
So the result at the 100th second is:-
1)pressure
2) Velocity
From the figure of velocity contour in the space which has been displated between the two plates in simulation, it seems that boundary layer has been formed and reached a steady state.
A line plot of velocity across the cross-section at x=0.2 will be displayed
The maximum velocity the fluid reaches in the simulation is 0.00376473 m/s which is quite close to the maximum velocity stipulated theoretically which is equal to 0.0037650 m/s.
This displays conformity of theoretical values with simulation results and conveys the fact that formula derived in theory will hold its validity in simulation till the assumptions held in theory are valid in environment to be simulated. In this case this means that till reynold number is as low as below 100 to enforce a pure laminar flow, simulation will arrive at steady state will conform to laminar flow theory.
Simulation of flow through the stepped plate:-
Now that the entire theoretical validity of formulae has been established. Simulation of fluid will be performed in the stepped space and the results will be compared for different scaling factors applied at the boxes. These scaling factors will be applied in such a manner that mesh is finer near the step corner and the walls. These scaling factors will be applied in both x as well as y direction. The corner at the step ,and the upper and lower walls are a source of disturbance from which any slight nuances are required to be captured.
Though grading factor is different, but the constantDict file,initial condition (0) files are same for each one of them.
1) ControlDict file
/*--------------------------------*- 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 400;
deltaT 0.001;
writeControl timeStep;
writeInterval 1000;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
Initial Condition:-
1)Pressure
/*--------------------------------*- 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
{
inletWall
{
type fixedValue;
value uniform 0.00001512;
}
outletWall
{
type fixedValue;
value uniform 0;
}
fixedWalls
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
2)Velocity
/*--------------------------------*- 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
{
inletWall
{
type zeroGradient;
}
outletWall
{
type zeroGradient;
}
fixedWalls
{
type noSlip;
}
frontAndBack
{
type empty;
}
}
BlockMesh
The entire domain space is divided into 5 cuboids. These cuboids are numbered at the middle and and there location is displayed in the above image. The Front side and back side has been shown to exhibit numbering of points as listed in the blockMeshDict file. Simplegrading facility has been used to make the mesh finer at the upper and lower walls since those are the positions where no-slip condition affects the velocity profile through the corossectional area. The simplegrading is unidirectional grading for control volumes 1,2,3 and 4 but it is multidirectional grading in the y-direction(unidirectional grading in x-direction) for control volume 5 because the mesh needs to be finer at the bottom floor and also near the step corner.
If s is the grading factor in consideration for analysis
1)Cuboid 1
For this cuboid, mesh needs to be finer at plane (1 12 13 2), in the x-direction and at plane(0 1 12 11), in the y-direction.
So x-direction-grading factor is s
y-direction-grading factor is (1/s)
2)Cuboid 2
For this cuboid, mesh needs to be finer at plane (2 3 14 13), in the x-direction and at plane(4 3 14 15), in the y-direction.
So x-direction-grading factor is s
y-direction-grading factor is s
2)Cuboid 3
For this cuboid, mesh needs to be finer at plane (1 12 13 2), in the x-direction and at plane (12 17 6 1), in the y-direction.
So x-direction-grading factor is (1/s)
y-direction-grading factor is (1/s)
2)Cuboid 4
For this cuboid, mesh needs to be finer at plane (2 3 14 13), in the x-direction and at plane (14 19 8 3), in the y-direction.
So x-direction-grading factor is (1/s)
y-direction-grading factor is s
2)Cuboid 5
For this cuboid, mesh needs to be finer at plane (1 12 20 9), in the x-direction and at plane (20 21 10 9) and (12 17 6 1), in the y-direction.
So x-direction-grading factor is (1/s)
But for y-direction the grading is multidirectional
It will be done in this manner(fraction_of_length fraction_of_total_cells grading_factor)
(0.5 0.5 s) for the plane(12 17 6 1)
(0.5 0.5 (1/s)) for the plane (20 21 10 9)
Z-direction is redundant characteristic added to the geometry for presentation and hence does not need any downgrading in the mesh.
Analysis will be done by plotting velocity profiles at 0.08m, 0.081m, 0.085m and 0.2m.
1) Grading factor 0.2
Following the methodology for application of scaling factor as jotted down above, the stepped domain space will be discretized for the scaling grade facter 0.2.
So s=0.2
BlockMeshDict file
/*--------------------------------*- 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)
(0.08 0 0)
(0.08 0.005 0)
(0.08 0.01 0)
(0 0.01 0)
(0 0.005 0)
(0.2 0 0)
(0.2 0.005 0)
(0.2 0.01 0)
(0.08 -0.01 0)
(0.2 -0.01 0)
(0 0 0.01)
(0.08 0 0.01)
(0.08 0.005 0.01)
(0.08 0.01 0.01)
(0 0.01 0.01)
(0 0.005 0.01)
(0.2 0 0.01)
(0.2 0.005 0.01)
(0.2 0.01 0.01)
(0.08 -0.01 0.01)
(0.2 -0.01 0.01)
);
blocks
(
hex (0 1 2 5 11 12 13 16) (80 5 1) simpleGrading (0.2 5 1)
hex (5 2 3 4 16 13 14 15) (80 5 1) simpleGrading (0.2 0.2 1)
hex (1 6 7 2 12 17 18 13) (120 5 1) simpleGrading (5 5 1)
hex (2 7 8 3 13 18 19 14) (120 5 1) simpleGrading (5 0.2 1)
hex (9 10 6 1 20 21 17 12) (120 10 1)
simpleGrading
(
5
(
(0.5 0.5 5)
(0.5 0.5 0.2)
)
1
)
);
edges
(
);
boundary
(
inletWall
{
type patch;
faces
(
(0 11 16 5)
(5 16 15 4)
);
}
outletWall
{
type patch;
faces
(
(6 17 21 10)
(7 18 17 6)
(8 19 18 7)
);
}
fixedWalls
{
type wall;
faces
(
(3 4 15 14)
(3 14 19 8)
(0 1 12 11)
(9 20 12 1)
(9 10 21 20)
);
}
frontAndBack
{
type empty;
faces
(
(0 5 2 1)
(5 4 3 2)
(13 14 15 16)
(11 12 13 16)
(2 3 8 7)
(1 2 7 6)
(9 1 6 10)
(14 13 18 19)
(13 12 17 18)
(12 20 21 17)
);
}
);
mergePatchPairs
(
);
By typing blockMesh in the command window, mesh is generated and stored in the polymesh in the constant folder.
By typing paraFoam command to view the mesh in the paraview,
Now icoFoam solver will be run for 400 seconds and again the result at the end will be viewed in paraview
Pressure contour
Velocity Contour
Now its important to view velocity profiles
i)Velocity profile at x=0.04m
Maximum velocity is 0.00201929m/s
ii)Velocity profile at x=0.081
The velocity profile is from y=-0.01 to y=0.01
Maximum velocity is 0.00191958m/s
iii)Velocity profile at 0.085m
Maximum velocity is 0.00166245m/s
iv)Velocity profile at x=0.2
Maximum velocity is 0.00103938m/s
v)Pressure from x=0 to x=0.2 at y=0.005 at each point.
Pressure at x=0.08 where transition happens in terms of width = 2.29318e-6 `m^2`/`s^2`
vi)Pressure from x=0.08 to x=0.2 at y=-0.005 at each point.
Maximum pressure is 2.15027e-6 `m^2`/`s^2`at x=0.095
Now CFL number after 400th second in the computer simulation:-
Maximum CFL number is 0.00516177
2)Grading factor 0.5
s=0.5
BlockMeshDict file
/*--------------------------------*- 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)
(0.08 0 0)
(0.08 0.005 0)
(0.08 0.01 0)
(0 0.01 0)
(0 0.005 0)
(0.2 0 0)
(0.2 0.005 0)
(0.2 0.01 0)
(0.08 -0.01 0)
(0.2 -0.01 0)
(0 0 0.01)
(0.08 0 0.01)
(0.08 0.005 0.01)
(0.08 0.01 0.01)
(0 0.01 0.01)
(0 0.005 0.01)
(0.2 0 0.01)
(0.2 0.005 0.01)
(0.2 0.01 0.01)
(0.08 -0.01 0.01)
(0.2 -0.01 0.01)
);
blocks
(
hex (0 1 2 5 11 12 13 16) (80 5 1) simpleGrading (0.5 2 1)
hex (5 2 3 4 16 13 14 15) (80 5 1) simpleGrading (0.5 0.5 1)
hex (1 6 7 2 12 17 18 13) (120 5 1) simpleGrading (2 2 1)
hex (2 7 8 3 13 18 19 14) (120 5 1) simpleGrading (2 0.5 1)
hex (9 10 6 1 20 21 17 12) (120 10 1)
simpleGrading
(
2
(
(0.5 0.5 2)
(0.5 0.5 0.5)
)
1
)
);
edges
(
);
boundary
(
inletWall
{
type patch;
faces
(
(0 11 16 5)
(5 16 15 4)
);
}
outletWall
{
type patch;
faces
(
(6 17 21 10)
(7 18 17 6)
(8 19 18 7)
);
}
fixedWalls
{
type wall;
faces
(
(3 4 15 14)
(3 14 19 8)
(0 1 12 11)
(9 20 12 1)
(9 10 21 20)
);
}
frontAndBack
{
type empty;
faces
(
(0 5 2 1)
(5 4 3 2)
(13 14 15 16)
(11 12 13 16)
(2 3 8 7)
(1 2 7 6)
(9 1 6 10)
(14 13 18 19)
(13 12 17 18)
(12 20 21 17)
);
}
);
mergePatchPairs
(
);
So generating the mesh the same way as done for grading factor 0.2 and this mesh is stored in polymesh folder at constant folder. To view this mesh, give command by typing paraFoam.
Now by giving command by typing icoFoam in the command prompt, the solver solves the problem for 400 seconds as directed in controlDict file.
1)Pressure
2) Velocity
Now its important to view velocity profiles
i)Velocity profile at x=0.04m
Maximum velocity is 0.00201875m/s
ii)Velocity profile at x=0.081
The velocity profile is from y=-0.01 to y=0.01
Maximum velocity is 0.0019082m/s
iii)Velocity profile at 0.085m
Maximum velocity is 0.00164059m/s
There seems to be flow with small velocitiy in negative direction close to the step wall
iv)Velocity profile at x=0.2
Maximum velocity is 0.00102752m/s
v)Pressure from x=0 to x=0.2 at y=0.005 at each point.
Pressure at x=0.08 where transition happens in terms of width =2.30054e-6 `m^2`/`s^2`
vi)Pressure from x=0.08 to x=0.2 at y=-0.005 at each point.
Maximum pressure is 2.1399e-6 `m^2`/`s^2`at x=0.095
This graph explains why there is a negative velocity near the step wall in the fluid. This is because fluid at the step wall expriences a pressure gradient in positive y-direction. This is because the fluid is stationary at the step wall but as we go in positive y-direction the fluid flow suddenly graspes a high velocity and therefore lower pressure than that in fluid near step wall due to the velocity. this cuses the flow upwards. This creates a decrease in pressure at the step wall because of fluid flow at the step wall upward but no incoming fluid at the wall. This creates a momentary vacuum which must be filled with fluid flowing from somewhere else at the step wall. So at points in space near x=0.05 at y close to -0.01 a positive gradient of pressure instead of negative pressure gradient occurs as shown in the graph above. This leads to negative fluid velocity flow at x=0.085 near y=-0.01
This phenomena will be observed irrespect of scaling factor used.
Now lets observe CFL number after 400th second in the computer simulation:-
Maximum CFL number is 0.00303659
3)Grading factor 0.8
s=0.8
BlockMeshDict file
/*--------------------------------*- 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)
(0.08 0 0)
(0.08 0.005 0)
(0.08 0.01 0)
(0 0.01 0)
(0 0.005 0)
(0.2 0 0)
(0.2 0.005 0)
(0.2 0.01 0)
(0.08 -0.01 0)
(0.2 -0.01 0)
(0 0 0.01)
(0.08 0 0.01)
(0.08 0.005 0.01)
(0.08 0.01 0.01)
(0 0.01 0.01)
(0 0.005 0.01)
(0.2 0 0.01)
(0.2 0.005 0.01)
(0.2 0.01 0.01)
(0.08 -0.01 0.01)
(0.2 -0.01 0.01)
);
blocks
(
hex (0 1 2 5 11 12 13 16) (80 5 1) simpleGrading (0.8 1.25 1)
hex (5 2 3 4 16 13 14 15) (80 5 1) simpleGrading (0.8 0.8 1)
hex (1 6 7 2 12 17 18 13) (120 5 1) simpleGrading (1.25 1.25 1)
hex (2 7 8 3 13 18 19 14) (120 5 1) simpleGrading (1.25 0.8 1)
hex (9 10 6 1 20 21 17 12) (120 10 1)
simpleGrading
(
1.25
(
(0.5 0.5 1.25)
(0.5 0.5 0.8)
)
1
)
);
edges
(
);
boundary
(
inletWall
{
type patch;
faces
(
(0 11 16 5)
(5 16 15 4)
);
}
outletWall
{
type patch;
faces
(
(6 17 21 10)
(7 18 17 6)
(8 19 18 7)
);
}
fixedWalls
{
type wall;
faces
(
(3 4 15 14)
(3 14 19 8)
(0 1 12 11)
(9 20 12 1)
(9 10 21 20)
);
}
frontAndBack
{
type empty;
faces
(
(0 5 2 1)
(5 4 3 2)
(13 14 15 16)
(11 12 13 16)
(2 3 8 7)
(1 2 7 6)
(9 1 6 10)
(14 13 18 19)
(13 12 17 18)
(12 20 21 17)
);
}
);
mergePatchPairs
(
);
So generating the mesh the same way as done for grading factor 0.2 and this mesh is stored in polymesh folder at constant folder. To view this mesh, give command by typing paraFoam.
Now by giving command by typing icoFoam in the command prompt, the solver solves the problem for 400 seconds as directed in controlDict file.
1)Pressure
2) Velocity
Now its important to view velocity profiles
i)Velocity profile at x=0.04m
Maximum velocity is 0.00201851m/s
ii)Velocity profile at x=0.081
The velocity profile is from y=-0.01 to y=0.01
Maximum velocity is 0.00190397m/s
iii)Velocity profile at 0.085m
Maximum velocity is 0.0016744m/s
iv)Velocity profile at x=0.2
Maximum velocity is 0.00102467m/s
v)Pressure from x=0 to x=0.2 at y=0.005 at each point.
Pressure at x=0.08 where transition happens in terms of width =2.32317e-6 `m^2`/`s^2`
vi)Pressure from x=0.08 to x=0.2 at y=-0.005 at each point.
Maximum pressure is 2.13892e-6 `m^2`/`s^2` at x=0.095
Now CFL number after 400th second in the computer simulation:-
Maximum CFL number is 0.00241154
Now that simulation has been done for all the three scaling factors, it is important to establish conformity of some theoretical basis with the simulation results. It is seen that the pressure graph from x=0 to x=0.2 is not continuous line i.e a straight line. It has one slope till the step at x=0.08 and another slope after 0.08m. Lets establish theoretically what must be the intermediate pressure at the step to be compared with the same obtained from simulation.
By conservation of mass, mass flow rate at the left side of the step must be equal to mass flow rate at the right side. But since density is same for both left and right side since it is incompressible flow, volumetric flow rate is also the same.
So
`Q_1=Q_2`
But it is known that formula for Q is
`Q=(dp.w^3)/(12*mu*L)`for flow between plates of infinite breadth
but for first section
`dp_1=P0-P_i`
and for 2nd section
`dp_2=P_i-P_(end)`
where
P0 is the pressure exerted at the beginning of the flow and Pend is the pressure at the end of the flow. Both pressures are guage pressure and not absolute pressure.
`P_i`or Pint is the intermediate pressure at the step. When the steady state flow is fully developed, an intermediate pressure at the step will occur to conserve the flow rate accross changing cross-section
So in accordance with the formula
For the first section
`Q_1=(dp_1.w_1^3)/(12*mu*L_1)`
So after simplification
`Q_1=((P0-P_i).w_1^3)/(12*mu*L_1)`
Similarly for second section
`Q_2=(dp_2.w_2^3)/(12*mu*L_2)`
after simplification
`Q_2=((P_i-Pend).w_2^3)/(12*mu*L_2)`
Since volume flow rate is conserved,
`Q_1=Q_2`
therefore
`((P0-P_i).w_1^3)/(12*mu*L_1)=((P_i-Pend).w_2^3)/(12*mu*L_2)`
further simplifying
`((P0-P_i).w_1^3)/(L_1)=((P_i-Pend).w_2^3)/(L_2)`
Now since, formula for intermediate pressure was the purpose, therefore after simplifying we get
`P_i=(P0*w_1^3*L_2+Pend*w_2^3*L_1)/(w_1^3*L_2+w_2^3*L_1)`
now since formula for `U_(x max)`was developed for the flow as a function for pressure difference which will be used in place of the formula as a function of Reynolds number
`U_(xmax)=(dp.w^2)/(8.mu.L)`
but since dp that is used in the software as kinematic pressure therefore
`dp_k=(dp)/rho`
substituting in the above formula above
`U_(xmax)=(dp_k*rho*w^2)/(8.mu.L)`
by further simplifying the formula we get
`U_(xmax)=(dp_k*w^2)/(8.nu.L)`
in the formula for `P_i`we got the formula for norm pressure but no matter the form or unit of the P0 and Pend the formula stands true for all and give the value for `P_i`in the same units as P0 and Pend.
Now these formulas which were developed theoretically will be transformed into a code.
Now there are steps followed in the coding which includes from the point dp calculation for the entire length to the prediction of CFL numbers for the three cases of grading factors. Therefore the steps and subsequent formulae that entail the coding is:-
1)Calculation of dp across the whole length of the stepped plate
L=0.2m , w=0.02m and Re=50
`dp=(12*(mu^2)*L*Re)/(rho^2*(w^3))`
2)Calculation of velocity maximum in first section and second section
`dp_1`and `dp_2` are calculated for the two sections
Uxmax formula is calculated for both the section as `U_(xmax1)`and `U_(xmax2)`
3)Calculation of upper limit of dt value for stable simulation
CFL is taken tentatively as 0.01 , dx is taken as calculated without simple grading and velocity maximum is taken as minimum of velocity maximum for both the section
`dt=(CFL)*dx/(vmax_(min))`
4)Prediction of maximum CFL number as it occurs in the simulation using selected value of dt
dt value is selected based on the upper limit of dt value
maximum of the velocity maximum for the two sections is taken and then three new cases are taken in which dx value is altered based on the three scaling factors.
And then based on the dt value selected, altered dx value and maximum of the velocity maximum, maximum CFL number is predicted as it should appear in the calculations in the simulation.
So the code is
clear all
close all
clc
%water is chosen
viscosity=1.002*10^-3;
density=998;
nu=viscosity/density;
%calculation of dp across the length of the stepped plate
L=0.2;%total length
Re=50;%desired Reynolds number
w=0.02; % bigger width for lower dp
dp=(12*(viscosity^2)*L*Re)/(density^2*(w^3))%calculation of dp accross the whole length
%Calculation for stepped plate
%first section
w1=0.01;
L1=0.08;
%second section
w2=0.02;
L2=0.12;
%Calculation of intermediate Pressure
Pend=0;
P0=Pend+dp;
Pint=(P0*L2*w1^3+Pend*L1*w2^3)/(L2*w1^3+L1*w2^3)
%Calculation of velocity maximum in first section and second section
dp1=P0-Pint;
dp2=Pint-Pend;
uxmax1=(dp1*w1^2)/(8*nu*L1)%first section
uxmax2=(dp2*w2^2)/(8*nu*L2)%Second section
%calculation for upper limit of dt using approximate CFL number
vmax_min=min([uxmax1 uxmax2]);
CFL=0.01;%tentative value of CFL
dx=0.001;%tentative value of dx (simple grading not applied)
dt=(CFL)*dx/vmax_min%upper limit of dt
%For prediction of maximum CFL number as it occurs in the simulation using selected value of dt
dtnew=0.001%selected dt
vmax_max=max([uxmax1 uxmax2]);
%scaling factor 0.2
scaling_factor_point_2=0.2
dxnew_point_2=dx*scaling_factor_point_2;%smallest dx when simple grading is done
newCFL_point_2=vmax_max*dtnew/dxnew_point_2%based on selected dt
%scaling factor 0.5
scaling_factor_point_5=0.5
dxnew_point_5=dx*scaling_factor_point_5;%smallest dx when simple grading is done
newCFL_point_5=vmax_max*dtnew/dxnew_point_5%based on selected dt
%scaling factor 0.8
scaling_factor_point_8=0.8
dxnew_point_8=dx*scaling_factor_point_8;%smallest dx when simple grading is done
newCFL_point_8=vmax_max*dtnew/dxnew_point_8%based on selected dtAnd the output in the command window is
Now that theoretical values have been calculated, they can be compared with the simulation values. Lets make a table for better view of the disparity in the values for theory as well as for values for all the cases.
Note:-
Value of maximum velocity is calculated theoretically for section 1 at x=0.08(at the step) but it is to be understood that the maximum velocity in the velocity profile is the same through out the length 0.08m till the step theoretically. So if maximum velocity is captured in the simulation at x=0.04m it is the same maximum velocity that is supposed to be maximum velocity at x=0.08m theoretically. But due to crossectional area change at the step gives rise to diffusion which dials down the velocity profile at x=0.08m and since the purpose of capturing velocity at x=0.08 was to get maximum velocity of the steady state profile developed across the length and then to be compared with theoretical velocity, therefore maximum velocity is taken at x=0.04m from the simulation and compared with the theoretical steady state maximum velocity for section 1(before the step)
These cases were tabulated in excel sheet and line graph was plotted for each to show a trend.
1) Velocity maximum for section 1
Trend for the absolute difference between theoretical and simulation values
2) Velocity maximum for section 2
Trend for the absolute difference between theoretical and simulation values
Analysis :-
Though simulation values have strong conformity with theoretical values but it seems from the trend that as scaling factor increases, the absolute error decreases which implies that greater scaling factor means greater conformity with theoretical values.
Regarding velocity maximum at both the sections, what happens is that maxima occurs at the mid of the section. But the meshing is coarser at the center of the cross_section owing to the simple grading that is done to make the mesh finer at the walls. Since maxima happens at the center of the cross-section, coarser mesh degrades the accuracy of the velocity. Since lesser the grading factor, coarser the mesh at the center of crossection which inturn leads to less accuracy. This is why with greater the scaling factor more is the accuracy observed in velocity maximum.
Now it seems that scaling factor is of no use at all and it may seem that it decreases accuracy. But if velocity profile at x=0.2m are compared for all the three scaling factor, it can be observed that with scaling factor of 0.2 velocity at the border is closest to zero and with scaling factor 0.8 velocity at the border is farthest from zero. So making mesh finer at the points where there may be sharp changes in variables can help capture values with accuracy.
3) Intermediate Pressure
Trend for the absolute difference between theoretical and simulation values
Regarding pressure values, pressure is measured at y=0.005 which is half the width for section 1 and quarter the width of section 2. It may seem that since pressure doesn't have a profile and is unaffected by grading across the cross-section, how come error is the most for longitudinal grading using 0.2 and the least using 0.8. It is because what may seem as error is not error at all. This is because theoretical calculation done to calculate pressure at the step doesn't take into account the effect of diffusion. Since longitudinally mesh is finer at the step for scaling factor 0.2 than that for scaling factor 0.8, so effects of diffusion is captured by using scaling factor 0.2 with greater accuracy than that done by scaling factor 0.8. So in actuality 0.2 scaling factor captures the realistic value of pressure and hence is more accurate than tht done under scaling factor 0.8. So what seems to be an error is basically the real effect due to diffusion captured in the supposed error.
4)CFL value comparison between theoretical and those obtained from simulation
Trend for the absolute difference between theoretical and simulation values
error in the Maximum CFL number seems to decrease as longitudinally as scaling factor increases. This gives the impression that using grading factor reduces accuracy of the whole simulation results. But it is more complicated than that.
First it has to be understood how CFL was calculated theoretically
Theoretically it was assumed that maximum velocity for section 1 is constant through out section 1 longitudinally. But it was explained before that velocity dials down in the simulation which causes the velocity to be less than what was calculated theoretically at the step. But the decrease is not that high to cause such a significant difference between cfl calculated theoretically and in simulation by computer. Therefore issue is not with the dialing down of velocity at the step.
The issue is that it is assumed that maximum cfl number will occur at the step with the maximum velocity in its profile at the step with lowest dx value due to longitudinal grading in the section 1. What happens is that as fluid approaches the step the dx value starts to decrease because of the grading. this should increase CFL number as fluid approaches the step. But due to diffusion at the step,velocity starts to decrease as fuid approaches the step. This would give a overall decreasing effect in CFL number value. So the maximum velocity remains till few millimeters away from the step after which velocity starts to dial down. So the point at which the velocity remains maximum is the point where dx should be found out to calculate maximum cfl number. But this point in space is same for all the three cases but dx value at this point varies a lot from dx value at the step for scaling factor 0.2 and a little less for scaling factor 0.5 and a lot less for scaling factor 0.8. If gradient of dx values is to be considered for simple grading from x=0 to x=0.08, then it is steep for scaling factor 0.2 and close to flat for scaling factor 0.8. So the point till which the velocity maximum is intact with negligible effects due to diffusion is the point where dx value stands to have a value with a significant difference in value for dx at the step for grading factor 0.2 but would have almost negligible difference in value with the dx value at the step for grading factor 0.8.
This is the reason that there is a significant difference between theoretical and simulation CFL value for grading factor 0.2 and much lesser difference between them for grading factor 0.8.
Conclusion:-
It was a right decision to keep mesh finer at the walls and near the step which captured the effects of zero velocity due to no slip condition at the walls and effects of diffusion at the step respectively. But at the same time it has to be kept in mind to keep the mesh finer as well where the variable stands to obtain maximum value which is in this case at the centre of the cross-section.
It also has to be kept in mind that applying simple grading should not make the mesh coarser while making the mesh finer at other points in space. If grading has to be applied then the first thing that has to be ensured is that mesh without simple grading is in itself appropriate enough to carry out simulation with accurate results and then assign simple-grading factor to make mesh finer where gradient in variable value will be steep like it was for velocity near the walls