Simulation of Flow through a pipe in OpenFOAM
Known:
Fluid: Water
Reynolds number: 2100
Density: 997 kg/m3
Kinematic viscosity: 1.004e-6 m2/s
OpenFOAM solver to be used: IcoFoam
HAGEN-POISEUILLE FLOW:
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.
Shear stress distribution:
`tau=(-(delp)/(delx))*r/2`
Velocity distribution:
`u=(1/(4mu))*(-(delp)/(delx))*(R^2-r^2)`
Pressure drop:
`delta P=(32*mu*V*L)/(D^2) = 105.041 Pa
HYDRODYNAMIC ENTRANCE LENGTH:
It is the length required by the viscous fluid to develop a fully developed flow in the pipe. Hydrodynamic entrance length depends on the viscosity of the flowing fluid, diameter of the pipe and Velocity of flow
Assumption:
diameter: 10mm
Flow is incompressible, irrotational and laminar
Calculation of Length of pipe:
L=Le+0.3
Le= 0.06*Re*D=1.26m
L=1.56m
Calculation of Average and Maximum velocity:
We know that:`Re=(V*D)/nu`
Hence: V= 2100*0.010/1.004e-6 = 0.21084 m/s (average velocity)
Maximum velocity: 2*Vavg= 0.42168 m/s
Pressure drop: `delta P=(32*mu*V*L)/(D^2) = 105.041 Pa
Maximum shear stress= 2*mu*V/r= 0.084 N/m2
Matlab code for automatic creation of blockmesDict file for any wedge angle and mesh grading factor:
%creating blockmeshdict file
close all
clear all
clc
L=1.56; %Length of pipe
theta=45; % angle of wedge
D=0.01; %Diameter of pipe
r=0.01/2; %radius of pipe
nx=30; %elements in x-direction
ny=1; %elements in y-direction
nz=200; %elements in z-direction
block=fopen(\'blockMeshDict\',\'w\');
fprintf(block,\'FoamFile\\n{\\n\');
fprintf(block,\'%12s \\t 2.0; \\n\',\'version\');
fprintf(block,\'%11s \\t ascii; \\n\',\'format\');
fprintf(block,\'%10s \\t dictionary; \\n\',\'class\');
fprintf(block,\'%11s \\t blockMeshDict; \\n}\\n\\n\',\'object\');
fprintf(block,\'convertToMeters 1;\\n\\n\');
fprintf(block,\'vertices\\n(\\n\');
fprintf(block,\'\\t(%f %f %f)\\n\',0,0,0); %vertex 0
fprintf(block,\'\\t(%f %f %f)\\n\',0,r*cosd(theta/2),r*sind(theta/2)); %vertex 1
fprintf(block,\'\\t(%f %f %f)\\n\',0,r*cosd(theta/2),-r*sind(theta/2)); %vertex 2
fprintf(block,\'\\t(%f %f %f)\\n\',L,0,0); %vertex 3
fprintf(block,\'\\t(%f %f %f)\\n\',L,r*cosd(theta/2),r*sind(theta/2)); %vertex 4
fprintf(block,\'\\t(%f %f %f)\\n);\\n\\n\',L,r*cosd(theta/2),-r*sind(theta/2)); %vertex 5
fprintf(block,\'blocks\\n(\\n\');
fprintf(block,\'\\thex (0 2 1 0 3 5 4 3) (%d %d %d) SimpleGrading (0.1 1 1)\\n);\\n\\n\', nx,ny,nz);
fprintf(block,\'edges\\n(\\n\');
fprintf(block,\'\\tarc 1 2(%f %f %f) \\n\', 0,r,0);
fprintf(block,\'\\tarc 4 5(%f %f %f) \\n);\\n\\n\', L,r,0);
fprintf(block,\'boundary\\n(\\n\');
fprintf(block,\'\\t inlet\\n\');
fprintf(block,\'\\t{\\n\');
fprintf(block,\'\\t\\ttype patch;\\n\');
fprintf(block,\'\\t\\tfaces\\n\');
fprintf(block,\'\\t\\t(\\n\');
fprintf(block,\'\\t\\t\\t(0 1 2 0)\\n\');
fprintf(block,\'\\t\\t);\\n\');
fprintf(block,\'\\t}\\n\');
fprintf(block,\'\\t outlet\\n\');
fprintf(block,\'\\t{\\n\');
fprintf(block,\'\\t\\ttype patch;\\n\');
fprintf(block,\'\\t\\tfaces\\n\');
fprintf(block,\'\\t\\t(\\n\');
fprintf(block,\'\\t\\t\\t(3 5 4 3)\\n\');
fprintf(block,\'\\t\\t);\\n\');
fprintf(block,\'\\t}\\n\');
fprintf(block,\'\\t top\\n\');
fprintf(block,\'\\t{\\n\');
fprintf(block,\'\\t\\ttype wall;\\n\');
fprintf(block,\'\\t\\tfaces\\n\');
fprintf(block,\'\\t\\t(\\n\');
fprintf(block,\'\\t\\t\\t(1 4 5 2)\\n\');
fprintf(block,\'\\t\\t);\\n\');
fprintf(block,\'\\t}\\n\');
fprintf(block,\'\\t back\\n\');
fprintf(block,\'\\t{\\n\');
fprintf(block,\'\\t\\ttype wedge;\\n\');
fprintf(block,\'\\t\\tfaces\\n\');
fprintf(block,\'\\t\\t(\\n\');
fprintf(block,\'\\t\\t\\t(0 2 5 3)\\n\');
fprintf(block,\'\\t\\t);\\n\');
fprintf(block,\'\\t}\\n\');
fprintf(block,\'\\t front\\n\');
fprintf(block,\'\\t{\\n\');
fprintf(block,\'\\t\\ttype wedge;\\n\');
fprintf(block,\'\\t\\tfaces\\n\');
fprintf(block,\'\\t\\t(\\n\');
fprintf(block,\'\\t\\t\\t(0 3 4 1)\\n\');
fprintf(block,\'\\t\\t);\\n\');
fprintf(block,\'\\t}\\n\');
fprintf(block,\'\\t axis\\n\');
fprintf(block,\'\\t{\\n\');
fprintf(block,\'\\t\\ttype empty;\\n\');
fprintf(block,\'\\t\\tfaces\\n\');
fprintf(block,\'\\t\\t(\\n\');
fprintf(block,\'\\t\\t\\t(0 3 3 0)\\n\');
fprintf(block,\'\\t\\t);\\n\');
fprintf(block,\'\\t}\\n\');
fprintf(block,\');\\n\\n\');
fprintf(block,\'mergePatchPairs\\n(\\n);\\n\');
fclose(block);BlockMeshDict file created from MATLAB code:
The vertices were created, after that from \'hex\' command the wedge block was formed. Two circular sectors was formed with the help of \'arc\' command. The boundary names were specified.
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.004999 0.000087)
(0.000000 0.004999 -0.000087)
(1.560000 0.000000 0.000000)
(1.560000 0.004999 0.000087)
(1.560000 0.004999 -0.000087)
);
blocks
(
hex (0 2 1 0 3 5 4 3) (30 1 200) SimpleGrading (0.1 1 1)
);
edges
(
arc 1 2(0.000000 0.005000 0.000000)
arc 4 5(1.560000 0.005000 0.000000)
);
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)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);Initial Velocity conditions i.e. U file:
Inlet: 0.21 m/s
Outlet: zerogradient boundary condition
top wall: no-slip boundary condition
Front wedge: wedge boundary condition
Back wedge: wedge boundary condition
Axis: empty
/*--------------------------------*- 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.210840 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //Initial Pressure conditions i.e. p file:
Inlet: zerogradient boundary condition
Outlet: 0 value
top wall: zerogradient boundary condition
Front wedge: wedge boundary condition
Back wedge: wedge boundary condition
Axis: empty
/*--------------------------------*- 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;
}
top
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //Control Dict file:
Time step size: 0.005 s
End time for simulation: 25 s
/*--------------------------------*- 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 25;
deltaT 0.005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Transport properties 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 \"constant\";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 1.004e-6;
// ************************************************************************* //Geometry and Meshing:
Below shown is the wedge geometry of the pipe with mesh grading of 0.1;
Number of elements in x-direction: 200
Number of elements in y-direction: 30
Number of elements in z-direction: 1
POST PROCESSING OF VELOCITY AND SHEAR STRESS:
Velocity profile at start & end of the pipe:
In the figure, the flow at the pipe has not been developed, we can see that as the flow proceeds through the length of the pipe, it gets distributed in a parabolic profile where flow at the wall stands still with zero velocity due to no-slip conditon, while the flow at the centre of the pipe reaching its maximum value
Flow at Start:
Flow at End:
Velocity Plots at various locations inside the pipe:
x = 0.1m
x = 0.5m
x = 1m
x = 1.26m (entrance length)
x = 1.5m
From the above plots, it can be concluded that, the flow at the pipe has not been developed,at x=0.1 m from the start, the flow is seen to be pretty undeveloped with a trapezoidal velocity plot, but we can see that as the flow proceeds through the length of the pipe, it gets distributed in a parabolic profile where flow at the wall stands still with zero velocity due to no-slip conditon, while the flow at the centre of the pipe reaching its maximum value. The value of maximum velocity for the fully developed steady state flow is 0.418255 m/s, while theoretically it should be 0.42168 m/s. Hence, the simulatio has an error of 0.812 %
This error can be reduced by many ways as follows:
- Refining the mesh of the wedge block
- Reducing the time step size
Shear stress plot:
The shear stress plot shows a good agreement with the theoretical concept of the laminar viscous pipe flow, which states that the shear stress is zero at the centre of the pipe and is maximum at the pipe walls (due to no-slip condition). From the plot, it can be seen that, the shear stress varies linearly along the Y-axis, with zero value at Y=0 .i.e. at the centre of the pipe and maximum at Y = radius i.e. at the walls