Week 11 - Simulation of Flow through a pipe in OpenFoam
OBJECTIVE :
Simulation of laminar flow through pipe in openfoam
INTRODUCTION :
Pipe flow regime laminar and turbulent is determined by Reynolds Number
`R_e=(ρvD)/μ`
for flow through pipe
Re < 2300 laminar flow
2300< Re < 4000 Transitional flow
Re > 4000 turbulent flow
HYDRODYNAMIC ENTRANCE LENGTH :
The flow has uniform velocity at the pipe entrance as shown in below fig. Because of the no-slip condition at the wall, we know that the velocity at the wall must be zero along the entire length of
the pipe. A boundary layer develops along the walls of the channel. The solid surface exerts a retarding shear force on the flow; thus the speed of the fluid in the neighborhood of the surface is reduced.
Sufficiently far from the pipe entrance, the boundary layer developing on the pipe wall reaches the pipe centerline and the flow becomes entirely viscous. The velocity profile shape then changes slightly after the inviscid core disappears. When the profile shape no longer changes with increasing distance x, the flow is called fully developed.
The distance downstream from the entrance to the location at which fully developed flow begins is called the entrance length.(Le)
For laminar flow `L_e=0.06 ReD`
For turbulent flow `L_e=4.4 (Re)^(1/6)D`
Hagen- Poiseuille's Equations :
In nonideal fluid dynamics, the 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.
Velocity Profile
`u(r)=-1/(4µ) (dP)/(dx) (R^2-r^2 )`
Average Velocity
`v_(aνg)=-R^2/8 (dP)/(dx)`
Maximum Veleocity
`v_(max)=2*v_(aνg)`
Flow Rate
`Q=-(ᴨR^4)/(8µ) (dP)/dx`
Pressure Drop
`(dP)/(dx)=(P_1-P_2)/L`
`P_1-P_2=(8*µ*L*v_(aνg))/R^4`
Head Loss
`H_(loss)=(P_1-P_2)/(ρg)`
Matlab Code For reating BlockMesh :
clear all
close all
clc
% Reynolds number for laminar flow
Re=2100;
% Piope Diameter
D=0.004;
%Pipe Radius
r=D/2;
% Hydrodynamic entry length of pipe
L_h=0.06*Re*D;
% Total Length of pipe
L=L_h+0.7;
% wedge angle
theta=4;
% Dynamic Viscosity
mu=0.89e-3;
% Density
rho=997;
%Equations for Hagen-Poiseuille Flow
v_avg=(Re*mu)/(rho*D);
v_max=2*v_avg;
del_p=(8*mu*L*v_avg)/(r^2);
kinematic_p=del_p/rho;
R=linspace(0,r,1000);
% shear stess
for i=1:length(R)
tau(i)=2*mu*v_max*(R(i)/r^2 );
end
plot(R,tau);
xlabel('Radius');
ylabel('Shear Stress');
legend ('Shear Stress Variation');
% Creating blockMeshDict file
header1='/*--------------------------------*- C++ -*----------------------------------*\'
header2=' ========= |'
header3=' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox'
header4=' \\ / O peration | Website: https://openfoam.org'
header5=' \\ / A nd | Version: 9'
header6=' \\/ M anipulation |'
header7='\*---------------------------------------------------------------------------*/'
header8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
f1=fopen('blockMeshDict','w')
fprintf(f1,'%s\n',header1)
fprintf(f1,'%s\n',header2)
fprintf(f1,'%s\n',header3)
fprintf(f1,'%s\n',header4)
fprintf(f1,'%s\n',header5)
fprintf(f1,'%s\n',header6)
fprintf(f1,'%s\n',header7)
fprintf(f1,'FoamFile\n{\n')
fprintf(f1,'%12s \t ascii;\n','format')
fprintf(f1,'%12s \t dictionary;\n','class')
fprintf(f1,'%12s \t blockMeshDict;\n','object')
fprintf(f1,'}\n')
fprintf(f1,'%s\n\n',header8)
fprintf(f1,'convertToMeters 1;\n\n')
fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1,'\t (%d %d %d)\n',0,0,0)
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1,'\t (%d %d %d)\n',L,0,0)
fprintf(f1,'\t (%d %d %d)\n',L,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1,'\t (%d %d %d)\n',L,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1,');\n\n')
fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t hex (0 3 5 2 0 3 4 1) (%d %d %d) simpleGrading (1 0.1 1)\n',300,30,1)
fprintf(f1,');\n\n')
fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1,'\t arc %d %d (%d %d %d)\n',1,2,0,r,0)
fprintf(f1,'\t arc %d %d (%d %d %d)\n',4,5,L,r,0)
fprintf(f1,');\n\n')
fprintf(f1,'boundary\n')
fprintf(f1,'(\n')
fprintf(f1,'\t inlet \n')
fprintf(f1,'\t {\n')
fprintf(f1,'\t\t type patch;\n')
fprintf(f1,'\t\t faces\n')
fprintf(f1,'\t\t (\n')
fprintf(f1,'\t\t\t (0 1 2 0)\n')
fprintf(f1,'\t\t );\n')
fprintf(f1,'\t }\n')
fprintf(f1,'\t outlet \n')
fprintf(f1,'\t {\n')
fprintf(f1,'\t\t type patch;\n')
fprintf(f1,'\t\t faces\n')
fprintf(f1,'\t\t (\n')
fprintf(f1,'\t\t\t (3 5 4 3)\n')
fprintf(f1,'\t\t );\n')
fprintf(f1,'\t }\n')
fprintf(f1,'\t top \n')
fprintf(f1,'\t {\n')
fprintf(f1,'\t\t type wall;\n')
fprintf(f1,'\t\t faces\n')
fprintf(f1,'\t\t (\n')
fprintf(f1,'\t\t\t (1 4 5 2)\n')
fprintf(f1,'\t\t );\n')
fprintf(f1,'\t }\n')
fprintf(f1,'\t front \n')
fprintf(f1,'\t {\n')
fprintf(f1,'\t\t type wedge;\n')
fprintf(f1,'\t\t faces\n')
fprintf(f1,'\t\t (\n')
fprintf(f1,'\t\t\t (0 3 4 1)\n')
fprintf(f1,'\t\t );\n')
fprintf(f1,'\t }\n')
fprintf(f1,'\t back \n')
fprintf(f1,'\t {\n')
fprintf(f1,'\t\t type wedge;\n')
fprintf(f1,'\t\t faces\n')
fprintf(f1,'\t\t (\n')
fprintf(f1,'\t\t\t (0 2 5 3)\n')
fprintf(f1,'\t\t );\n')
fprintf(f1,'\t }\n')
fprintf(f1,'\t axis \n')
fprintf(f1,'\t {\n')
fprintf(f1,'\t\t type empty;\n')
fprintf(f1,'\t\t faces\n')
fprintf(f1,'\t\t (\n')
fprintf(f1,'\t\t\t (0 3 3 0)\n')
fprintf(f1,'\t\t );\n')
fprintf(f1,'\t }\n')
fprintf(f1,');\n\n')
fprintf(f1,'mergePatchPairs\n')
fprintf(f1,'(\n')
fprintf(f1,');\n')
header9='// ************************************************************************* //'
fprintf(f1,'%s\n',header9)
fclose(f1)
Created BlockMesh File :
/*--------------------------------*- 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
(
(0 0 0)
(0 1.998782e-03 6.979899e-05)
(0 1.998782e-03 -6.979899e-05)
(1.204000e+00 0 0)
(1.204000e+00 1.998782e-03 6.979899e-05)
(1.204000e+00 1.998782e-03 -6.979899e-05)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (300 30 1) simpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 2.000000e-03 0)
arc 4 5 (1.204000e+00 2.000000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
Above two figure shows mesh formed by running blockMesh.It can be seen near the wall fine mesh is created to capture gradient in velocity more accurately.
ControlDict :
/*--------------------------------*- 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 10;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Transport Properties :
/*--------------------------------*- 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] 8.9256e-7;
// ************************************************************************* //
Initial And Boundary Condition :
For Velocity
/*--------------------------------*- 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.4686 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.4686 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
For Pressure
/*--------------------------------*- 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;
}
top
{
type zeroGradient;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
RESULTS :
Ab ove figure shows that at inlet velocity is constant except at wall , at wall because of no slip boundary condition it is zero.
As fluid flows further velocity at centre increases and there is core region formed at centre of pipe where velocity is nearly constant as shown in above figure.
Above figure shows boundary layer development in pipe. Region coloured in blue has low velocity is the region where boundary layer is formed, in this region viscus forses are dominant.
Above figure shows that after traveling nearly 0.5m to 0.6m of pipe velocity of does not change with length.This length is Hydrodynamic Entrance Length. After this flow is fully developed.
From equation `L_e=0.06 ReD` for Reynolds Number 2100 and Diametr 0.04m we get Le=0.5m.
Above figure shows velocity profile at fully developed region. At fully developed region velocity profile is parabolically varies withy radius.
Above figure shows pressure variation along length of pipe. Upto the entrance length there is high pressure drop which causes acceleration of flow and act against viscous effect. But in fully developed region pressure drop is linear because it is only used balance shear forces because of viscosity.
Above figure shows variation of shear stess with radius using Hagen- Poiseuille's equations and it can be seen that as radius increase shear stress increses linearly for laminar flow. At the wall shear stress is maximum and at the centre it is zero.
Above figure shows shear stress variation with radius obtained after simulation in openfoam using icoFoam solver.This graph is similar to the graph obtained using equations.