Week 11 - Simulation of Flow through a pipe in OpenFoam
Simulation of flow through a pipe in OpenFOAM
We consider a cylindrical pipe with the fluid water entering the pipe with the laminar flow. The aim is to simulate the flow through a section of pipe and compare the flow with the Hagen Poisuilles equation.
Hagen Poiseuuiles Equation
In nonideal fluid dynamics, the Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, Poiseuille law or Poiseuille equation, 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.
The assumptions of the equation are that the fluid is incompressible and Newtonian; the flow is laminar through a pipe of constant circular cross-section that is substantially longer than its diameter; and there is no acceleration of fluid in the pipe. For velocities and pipe diameters above a threshold, actual fluid flow is not laminar but turbulent, leading to larger pressure drops than calculated by the Hagen–Poiseuille equation.
In standard fluid-kinetics notation
`Deltap=(8*mu*L*Q)/A^2`
where,
`Deltap` is the pressure difference between the two ends,
`L`is the length of pipe,
`mu`is the dynamic viscosity,
`Q`is the vloumetric flow rate.
`R` is the pipe radius.
`A`is the cross section of pipe.
The equation does not hold close to the pipe entrance.
This is the diagramatic representation of flow through pipe.This is what the fluid profile looks like.
The layer in which the shearing viscous forces are significant, is called the boundry layer. This boundary layer is a hypothetical concept. It divides the flow in pipe into two regions:
- Boundary layer region: The region in which viscous effects and the velocity changes are significant.
- The irrotational (core) flow region: The region in which viscous effects and velocity changes are negligible, also known as the inviscid core.
When the fluid just enters the pipe, the thickness of the boundary layer gradually increases from zero moving in the direction of fluid flow and eventually reaches the pipe center and fills the entire pipe. This region from the entrance of pipe to the point where the boundary layer covers the entire pipe is termed as the hydrodynamic entrance region and length of the pipe in this region is termed as the hydrodynamic entry length. In this region the velocity profile develops and thus the flow is called the hydrodynamically developing flow. After this region, the velocity profile is fully developed and continues unchanged. This region is called the hydrodynamically fully developed region.
Given:
- Reynlods number :2100
- Working fluid :water .
- Fluid temperature:25 `overset(circ)C`
- Diameter of the pipe =0.02m
- Dynamic viscosity(`nu`)=`0.89times10^-3`
- Kinematic Viscosity =`0.89times10^-6``m^2s^-1`
- Also `DeltaP= (32*mu*L*V_("avg"))/(D^2)`
- `theta=4^o`
- L=2.52m
Analytical Solution:
- Average Velocity : `V_("avg") = (R_e*mu)/(rho*D)`= 0.09373`ms^-1`
- Maximum velocity :`V_max = 2*V_("avg")`= 0.18746`ms^-1`
- Pressure Drop = 16.817pa
- Kinematic Pressure Drop : `(DeltaP)/(rho` =0.016817m.
- We use the incompressible folder and we navigate to then use icoFoam for solving the problem.we copy the folder cavity and then bring it torun folder.and make the following changes for different files.
- set of commands:
- cp -r cavity $FOAM_RUN/simulation
- cd $FOAM_RUN/simulation/cavity
- gedit -filename.
- finally ,blockMesh and icoFoam for creating mesh and running simulation respectively.
- Then open paraView for plotting and other purpose of running the simulation.
Matlab code for creating the blockMeshDict file is as follows.
clear all
close all
clc
%wedge angles and inputs.
theta = input('Enter the value of wedge angle=');
Re =2100;
D=0.02;
L=0.06*Re*D;
R=D/2;
mu=0.89e-3;
rho=997.0;
nu=mu/rho;
%pressure and velocity
v_avg=Re*mu/(D*rho);
v_max=2*v_avg;
del_P =(32*mu*v_avg*L)/(D^2);
kin_P=del_P/rho;
n1=50;
n2=1;
n3=300;
gfactor=0.5;
%vertices of the wedge
V0=[0 0 0 ];
V1=[0 R*cosd(theta/2) -R*sind(theta/2)];
V2=[0 R*cosd(theta/2) R*sind(theta/2)];
V3=[L 0 0 ];
V4=[L R*cosd(theta/2) -R*sind(theta/2)];
V5=[L R*cosd(theta/2) R*sind(theta/2)];
%printing the contents of blockmesh file
string1='/*-------------------------------------*- C++ -*---------------------*';
string2=' ====== | ;';
string3=' / Field | ';
string4=' / Operation |OPENFOAM:the open source cfd toolbox ';
string5=' / And |Version :7 ';
string6=' / Mnaipulation |Web :www.OpenFoam.ORG ';
string7='*----------------------------------------------------------------------*/';
string8='FoamFile';
string9='{';
string10=' version 2.0;';
string11=' format acsii;';
string12=' class dictionary;';
string13=' object blockMeshDict;';
string14='}';
string15='//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
string16='//*******************************************************************//';
fpointer=fopen('blockMeshDict.txt','wt');
fprintf(fpointer,'%sn',string1);
fprintf(fpointer,'%sn',string2);
fprintf(fpointer,'%sn',string3);
fprintf(fpointer,'%sn',string4);
fprintf(fpointer,'%sn',string5);
fprintf(fpointer,'%sn',string6);
fprintf(fpointer,'%sn',string7);
fprintf(fpointer,'%sn',string8);
fprintf(fpointer,'%sn',string9);
fprintf(fpointer,'%sn',string10);
fprintf(fpointer,'%sn',string11);
fprintf(fpointer,'%sn',string12);
fprintf(fpointer,'%sn',string13);
fprintf(fpointer,'%sn',string14);
fprintf(fpointer,'%sn',string15);
fprintf(fpointer,'%sn',string16);
fprintf(fpointer,'%sn','convertToMeters 1;');
fprintf(fpointer,'%sn','vertices');
fprintf(fpointer,'%sn','(');
fprintf(fpointer,'%s (%f %f %f) n',blanks(5),V0(1),V0(2),V0(3));
fprintf(fpointer,'%s (%f %f %f) n',blanks(5),V1(1),V1(2),V1(3));
fprintf(fpointer,'%s (%f %f %f) n',blanks(5),V2(1),V2(2),V2(3));
fprintf(fpointer,'%s (%f %f %f) n',blanks(5),V3(1),V3(2),V3(3));
fprintf(fpointer,'%s (%f %f %f) n',blanks(5),V4(1),V4(2),V4(3));
fprintf(fpointer,'%s (%f %f %f) n',blanks(5),V5(1),V5(2),V5(3));
fprintf(fpointer,'%sn',');');
fprintf(fpointer,'n');
fprintf(fpointer,'%sn','blocks');
fprintf(fpointer,'%sn','(');
fprintf(fpointer,'%s %s n',blanks(5),'hex(0 1 2 0 3 4 5 3 ) (50 1 200) simpleGrading (0.5 1 1)');
fprintf(fpointer,'%sn',');');
fprintf(fpointer,'n');
fprintf(fpointer,'%sn','edges');
fprintf(fpointer,'%sn','(');
fprintf(fpointer,'%s arc 1 2 (%f %f %f) n',blanks(5),0,R,0);
fprintf(fpointer,'%s arc 4 5 (%f %f %f) n',blanks(5),0,R,0);
fprintf(fpointer,'%s n',');');
fprintf(fpointer,'n');
fprintf(fpointer,'%sn','boundry');
fprintf(fpointer,'%sn','(');
fprintf(fpointer,'%s %s n',blanks(5),'inlet');
fprintf(fpointer,'%s %s n',blanks(5),'{');
fprintf(fpointer,'%s %s n',blanks(10),'type patch;');
fprintf(fpointer,'%s %s n',blanks(10),'faces');
fprintf(fpointer,'%s %s n',blanks(10),'(');
fprintf(fpointer,'%s %s n',blanks(15),'(0 1 2 0 )');
fprintf(fpointer,'%s %s n',blanks(10),');');
fprintf(fpointer,'%s %s n',blanks(5),'}');
fprintf(fpointer,'%s %s n',blanks(5),'outlet');
fprintf(fpointer,'%s %s n',blanks(5),'{');
fprintf(fpointer,'%s %s n',blanks(10),'type patch;');
fprintf(fpointer,'%s %s n',blanks(10),'faces');
fprintf(fpointer,'%s %s n',blanks(10),'(');
fprintf(fpointer,'%s %s n',blanks(15),'(3 4 5 3)');
fprintf(fpointer,'%s %s n',blanks(10),');');
fprintf(fpointer,'%s %s n',blanks(5),'}');
fprintf(fpointer,'%s %s n',blanks(5),'front');
fprintf(fpointer,'%s %s n',blanks(5),'{');
fprintf(fpointer,'%s %s n',blanks(10),'type wedge;');
fprintf(fpointer,'%s %s n',blanks(10),'faces');
fprintf(fpointer,'%s %s n',blanks(10),'(');
fprintf(fpointer,'%s %s n',blanks(15),'(0 3 5 2)');
fprintf(fpointer,'%s %s n',blanks(10),');');
fprintf(fpointer,'%s %s n',blanks(5),'}');
fprintf(fpointer,'%s %s n',blanks(5),'back');
fprintf(fpointer,'%s %s n',blanks(5),'{');
fprintf(fpointer,'%s %s n',blanks(10),'type wedge;');
fprintf(fpointer,'%s %s n',blanks(10),'faces');
fprintf(fpointer,'%s %s n',blanks(10),'(');
fprintf(fpointer,'%s %s n',blanks(15),'(0 1 4 3)');
fprintf(fpointer,'%s %s n',blanks(10),');');
fprintf(fpointer,'%s %s n',blanks(5),'}');
fprintf(fpointer,'%s %s n',blanks(5),'pipe');
fprintf(fpointer,'%s %s n',blanks(5),'{');
fprintf(fpointer,'%s %s n',blanks(10),'type wall;');
fprintf(fpointer,'%s %s n',blanks(10),'faces');
fprintf(fpointer,'%s %s n',blanks(10),'(');
fprintf(fpointer,'%s %s n',blanks(15),'(1 4 5 2)');
fprintf(fpointer,'%s %s n',blanks(10),');');
fprintf(fpointer,'%s %s n',blanks(5),'}');
fprintf(fpointer,'%s %s n',blanks(5),'axis');
fprintf(fpointer,'%s %s n',blanks(5),'{');
fprintf(fpointer,'%s %s n',blanks(10),'type empty;');
fprintf(fpointer,'%s %s n',blanks(10),'faces');
fprintf(fpointer,'%s %s n',blanks(10),'(');
fprintf(fpointer,'%s %s n',blanks(15),'(0 3 3 0)');
fprintf(fpointer,'%s %s n',blanks(10),');');
fprintf(fpointer,'%s %s n',blanks(5),'}');
fprintf(fpointer,'%s n',');');
fprintf(fpointer,'n');
fprintf(fpointer,'%s n','mergePatchPairs');
fprintf(fpointer,'%s n','(');
fprintf(fpointer,'%s n',');');
fprintf(fpointer,'n');
fprintf(fpointer,'%s n',string16);
fclose(fpointer);
The blockMeshfile is created and then will be replaced by the the file created by the matlab.The blockMesh file can be seen as follows:
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 8
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009994 -0.000349)
(0.000000 0.009994 0.000349)
(2.5200000 0.000000 0.000000)
(2.520000 0.009994 -0.000349)
(2.520000 0.009994 0.000349)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
The boundry conditons can be given as follows .The velocity boundry condition were put in the problem discuused below.
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 8
\/ 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.0937312 0 0);
}
outlet
{
type zeroGradient;
}
pipe
{
type noSlip;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
The pressure boundry conitions were put according to the pressure drop and is as follows:
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 8
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
axis
{
type empty;
}
pipe
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0168682;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
The data related to kinematic viscosity is updated in transportProperties file.
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 8
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.927e-07;
// ************************************************************************* //The data related to simulation time duration etc is updated in controlDict file as shown below:
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 8
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
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;
// ************************************************************************* //The mesh will be created by this file which looks like this:
In different views it can be seen as be seen as follows:
The results of simulation is as follows:
We plot the velocity graph at
x=0.01 from the inlet. and along the y axis.The distance varies along x axis.But the direction of probe is along y.
x=0.5 from the inlet. and along the y axis.The distance varies along x axis.But the direction of probe is along y.
x=1.15 from the inlet. and along the y axis.The distance varies along x axis.But the direction of probe is along y.
x=2.52 from the inlet that is at outlet. and along the y axis.The distance varies along x axis.But the direction of probe is along y.We could say this a profile where the flow is fully develeoped.
Here we could sayy that the plots seen above represents velocity v/s disyance along the radius of the pipe.We could find that due to nosli pconditions along the boundry of the wall the velocity turns out to be zero.Also we would plot the the shear stress as below with the help of
gradient of velocity and then multiply with kinematic viscosity inorder to get suitable results.
All this can be done in paraView .just we need to check out some calculator functions in paraFoam.
x=1.15 from the inlet , and along the y axis.The distance varies along x axis.But the direction of probe is along y.The shear stress can be plotted as follows:
Shear stress profile:
- Generally shear stress in laminar flow is a result of momentum transfer among the randomly moving molecules in a flow.
- Analytically shear stress can be represented by the formulae:
- `tau = mu*(du)/(dy`
- where ,
- `tau` is the shear stress (F/A)
- `mu` is the dynamic viscosity
- `(du)/dy` rate of shear deformation.
Some of the conclusions drawn:
- The plot at x=0.01 the profile is not fully developed we could understand it by the plot itself.It means there is some distance for which we can attain full profile that is called enterace length ,developing flow as explained.
- it can be seen flow starts to develop as the distance x increases and the profile becomes more parabolic as well
- We obtain almost same plots for x=1.15 and x=2.52 hence we could conclude that the fully developed length is somewhere less then or very near to 1.15.
- Here we could notice that the velocity turns out to be zero at the walls its mainly because of the no slip condition.As the flow progresses gradually the boundry layers develop and the adjacent layers tend to slow downBeacuse of this the velocity at the center increases .Hence becuase of the pressure difference obtained from the equations velocity profile matches with the Hagen poiseuilles equation.
- When the fluid moves shear stresses are developed as the particles in the flow move relative to each other .For a fluid flowing in pipe.,the fluid velocity remains zero at the pipe wall .Shear force is generated because the adjacent layer of the fluid that move with differnt velocities.
- If value of `mu` is constant the shear stress is linearly depend on the velocity gradient.As the change in velocity is larger,more shear stress in the walls.
- Hence the computational mesh is generated for wedge angle degree.Also showed that the velocity profile matches the Hagen poiseulles equation .From the above plotted velocity profiles it is proved that the profile is fully developed.
- The plot for shear stress seems almost linear.It is plotted at 1.15 which nearly at the fully developed region.The linearity and smoothness in case of the plot is mainy due to the fine meshing .Which was created above by blockMeshDict file.
- This are the analytical solutions.
- Average Velocity :`V_("avg") = (R_e*mu)/(rho*D)` = 0.09373`ms^-1`.
- Maximum velocity :`V_max = 2*V_("avg")`= 0.18746`ms^-1`
- Pressure Drop = 16.817pa
- Kinematic Pressure Drop : `(DeltaP)/(rho`=0.016817m.
- By observing the graphs we could say that the analytivcal and numerical solutions are very much closer to each othe