Week 11 - Simulation of Flow through a pipe in OpenFoam
#Purpose of the Project:
1.Creation of Matlab Code for automatically generating computational mesh for any degree of wedge angle.
2.Velocity profile of outcomes from simulation matches with Hagen Poiseuille's equation.(icoFoam is used as solver)
3.Representaion of fully developed velocity profile.
4.Post processing of Velocity and Shear stress of flow through pipe in Paraview.
Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law ,is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section.
Hydrodynamic Entrance Length:
The entrance length is the length in a tube or duct after an obstruction - until the flow velocity profile is fully developed.
Problem values and variables considering water as a medium:
1.Diameter of Pipe- 0.02 (metres)
2.Radius of pipe - 0.01 m.
3.Density -997 kg/m^3.
4.Dynamic Viscocity-8.9*10^-4 N/m^2*s
5.Reynold Number -2100
6.Hydrodynamic Entrance Length (Le)
Le=0.06*Re*D
.Hydrodynamic Entrance Length (Le)
Le=0.06*Re*D
Le=2.52 m
7.Inlet Velocity (v) can be calculated by:
`Re=(rho*v*D)/mu`
v=0.093 `m*s^-1`
Hagen Poiseuille's equation Analytical Solution:
1.Average Velocity=`V^(avg)`=`(Re*mu)/(rho*D)`=0.093 `m*s^-1`
2.Max. Velocity ==`V^(max)=2*V^(avg)`=0.186 `m*s^-1`
3.Pressure Drop==`deltaP`= `(8*mu*Le*V^(avg))/R^2`=16.686`N*m^-2`
4.Normalised Pressure Drop=`(deltaP)/rho`=0.0167 `m^2*s^-2`
5.Shear Stress=`tau=mu*(du)/(dy)`
%%Matlab code for Generating Mesh automatically for any wedge angle
clear all;
close all;
clc
D=0.02%Diameter of pipe in m
R=D/2 %Radius of pipe
theta=input('Enter the value of theta in degrees---'); %%wedge angle in degrees
L=3; %%Length of pipe in m
%Wedge dimension at 4 degrees
Point1=D/2*cosd(theta/2);
Point2=D/2*sind(theta/2);
%%File for writing blockMesh code
f2=fopen('blockMeshDict.txt','w')
%%Vertex and other parameters defined
Vertex_0=[0 0 0];
Vertex_1=[0 0 L];
Vertex_2=[Point1 -Point2 0];
Vertex_3=[Point1 Point2 0];
Vertex_4=[Point1 -Point2 L];
Vertex_5=[Point1 Point2 L];
Block=[0 2 3 0 1 4 5 1]
Cells=[50 1 300]
Grading=[0.1 1 1];
Inlet_face=[0 2 3 0];
Outlet_face=[1 4 5 1];
wall=[2 3 5 4];
wedge1=[0 1 4 2];
wedge2=[0 1 5 3];
axis=[0 1 1 0];
fprintf(f2,['/*------------------------------*- C++ -*----------------------------------',...
'*\\',...
'\n| ========= | |',...
'\n| \\\\ / F ield | OpenFOAM: TheOpen Source CFD Toolbox |',...
'\n| \\\\ / O peration | Version: plus |',...
'\n| \\\\ / A nd | Web: www.OpenFOAM.com |'...
'\n| \\\\ / M anipulation | |',...
'\n\\*-------------------------------------------------------------------------*/',...
'\nFoamFile',...
'\n{',...
'\n version 2.0;',...
'\n format ascii;',...
'\n class dictionary;',...
'\n object blockMeshDict;',...
'\n}',...
'\n// ******************************************************************* //'])
fprintf(f2,'\n%s %.3f','convertToMeters', 1,';');
fprintf(f2,'\n')
fprintf(f2,'\nvertices');
fprintf(f2,'\n(');
fprintf(f2,'\n(');
fprintf(f2,'%12.6f',Vertex_0);
fprintf(f2,')');
fprintf(f2,'\n(');
fprintf(f2,'%12.6f',Vertex_1);
fprintf(f2,')');
fprintf(f2,'\n(');
fprintf(f2,'%12.6f',Vertex_2);
fprintf(f2,')');
fprintf(f2,'\n(');
fprintf(f2,'%12.6f',Vertex_3);
fprintf(f2,')');
fprintf(f2,'\n(');
fprintf(f2,'%12.6f',Vertex_4);
fprintf(f2,')');
fprintf(f2,'\n(');
fprintf(f2,'%12.6f',Vertex_5);
fprintf(f2,')');
fprintf(f2,'\n);');
fprintf(f2,'\n');
fprintf(f2,'\nblocks');
fprintf(f2,'\n(');
fprintf(f2,'\n');
fprintf(f2,' hex (');
fprintf(f2,'%2.0f',Block);
fprintf(f2,') ');
fprintf(f2,'(');
fprintf(f2,'%4.0f',Cells);
fprintf(f2,') ');
fprintf(f2,'simpleGrading')
fprintf(f2,' (');
fprintf(f2,'%5.1f',Grading);
fprintf(f2,')');
fprintf(f2,'\n);');
fprintf(f2,'\n');
fprintf(f2,'\nedge')
fprintf(f2,'\n(')
fprintf(f2,'\n);')
fprintf(f2,'\n')
fprintf(f2,'\nboundary')
fprintf(f2,['\n(',...
'\ninlet',...
'\n{',...
'\n type patch;',...
'\n faces',...
'\n(',...
'\n('])
fprintf(f2,'%2.0f',Inlet_face)
fprintf(f2,[')',...
'\n);',...
'\n}',])
fprintf(f2,['\noutlet',...
'\n{',...
'\n type patch;',...
'\n faces',...
'\n(',...
'\n('])
fprintf(f2,'%2.0f',Outlet_face)
fprintf(f2,[')',...
'\n);',...
'\n}',])
fprintf(f2,['\nwall',...
'\n{',...
'\n type wall;',...
'\n faces',...
'\n(',...
'\n('])
fprintf(f2,'%2.0f',wall)
fprintf(f2,[')',...
'\n);',...
'\n}',])
fprintf(f2,['\nfront',...
'\n{',...
'\n type wedge;',...
'\n faces',...
'\n(',...
'\n('])
fprintf(f2,'%2.0f',wedge1)
fprintf(f2,[')',...
'\n);',...
'\n}',])
fprintf(f2,['\nback',...
'\n{',...
'\n type wedge;',...
'\n faces',...
'\n(',...
'\n('])
fprintf(f2,'%2.0f',wedge2)
fprintf(f2,[')',...
'\n);',...
'\n}',])
fprintf(f2,['\naxis',...
'\n{',...
'\n type empty;',...
'\n faces',...
'\n(',...
'\n('])
fprintf(f2,'%2.0f',axis)
fprintf(f2,[')',...
'\n);',...
'\n}',])
fprintf(f2,'\n);')
fprintf(f2,'\nmergePatchPairs')
fprintf(f2,'\n(')
fprintf(f2,'\n);')
fprintf(f2,'\n// ********************************************************* //')
fclose(f2); %%File close after writing
Openfoam Case setup details:
For simulating flow in a pipe,axisymmetric condition is used.Wedge boundary condition is applied for only simulating for some degree of wedge angle like 2,3,4 degrees.
Solver used is icoFoam considering flow is incompressible,laminar.
BlockMesh file for creating Mesh at wedge angle =`2^o` :
%%blockMeshDict file for generating Mesh for Wedge angle equal 2 degree%%
/*------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: TheOpen Source CFD Toolbox |
| \\ / O peration | Version: plus |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\ / M anipulation | |
\*-------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// ******************************************************************* //
convertToMeters 1.000
;
vertices
(
( 0.000000 0.000000 0.000000)
( 0.000000 0.000000 3.000000)
( 0.009998 -0.000175 0.000000)
( 0.009998 0.000175 0.000000)
( 0.009998 -0.000175 3.000000)
( 0.009998 0.000175 3.000000)
);
blocks
(
hex ( 0 2 3 0 1 4 5 1) ( 50 1 300) simpleGrading ( 0.1 1.0 1.0)
);
edge
(
);
boundary
(
inlet
{
type patch;
faces
(
( 0 2 3 0)
);
}
outlet
{
type patch;
faces
(
( 1 4 5 1)
);
}
wall
{
type wall;
faces
(
( 2 3 5 4)
);
}
front
{
type wedge;
faces
(
( 0 1 4 2)
);
}
back
{
type wedge;
faces
(
( 0 1 5 3)
);
}
axis
{
type empty;
faces
(
( 0 1 1 0)
);
}
);
mergePatchPairs
(
);
// ********************************************************* //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 latestTime;
startTime 0;
stopAt endTime;
endTime 30;
deltaT 5e-3;
writeControl timeStep;
writeInterval 40;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Pressure Boundary Condition:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0167;
}
wall
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //Velocity Boundary Condition:
/*--------------------------------*- 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 (0 0 0.093);
}
outlet
{
type zeroGradient;
}
wall
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //For `2^o` wedge angle:
Velocity representation at Inlet
Velocity representation at Outlet:
Velocity profile at Inlet:
Velocity profile at 2 metre from inlet of pipe:
Velocity profile after reaching hydrodynamic length of pipe 2.52 m :
Hydrodynamic entrance length of the pipe dimensions used for this simulation is 2.52 metres.So after this region ,velocity profile is fully developed and will not change further.
Shear Stress graph at 2.6 m from inlet of the pipe:
Above shear stress plot shows that it is maximum at the wall and minimum at the centre of the pipe.
For `3^o` degree:
For`4^o`degree:
Outcomes:
1.Maximum velocity obtained from Hagen poiseuille's equation by Analytical solution is 0.186m/sec.
Inspite of using different wedge angle for simulation that is 2,3 and 4,Maximum velocity obtained is 0.18452.
Results obtained from simulation poses error of 0.00148 m/s.
2.Velocity profile is fully developed after hydrodynamic entrance length.
3.Shear stress profile across radial length matches with theoritical profile.