Flow through a pipe
Objective of the project: 1.The main objective of this project is to simulate the laminar flow through pipe. The wedge portion is taken as the whole circular section (cylinder) may take time to simulate or it may be costly if your computer do not have much resources.
2. The other objectives is to write the Matlab code for blockMeshDict file for different wedge angle , length and radius.
3. The final objective is to post process the results of velocity at different length of the pipe and shear stress at the fully developed portion and to compare it with the analytical solution i.e by hagen poiseuille equation
Assumptions made :
- Steady flow , incompressible and laminar flow
- Flow is two dimenional ( Flow WRT theta is ignored )
- Constant fluid properties (Viscosity , Density etc)
Flow Variables:
The working fluid is water
- Density = 1000 kg.m-3\">kg.m−3
- Dynamic viscosity = 8.90*10-4\">10−4Pa.s\">Pa.s
- `eynolds number = 210`
- Diameter of pipe = 0.01 m\">m
- To compute the flow velocity, we use the formula - Re=ρvDμ\">Re=ρvDμ
- Obtained Inlet velocity = 0.1869 ms-1\">ms−1
Boundary layer
The layer in which the shearing viscous forces are significant, is called the boundary layer.[6] This boundary layer is a hypothetical concept. It divides the flow in pipe into two regions:[6]
- Boundary layer region: The region in which viscous effects and the velocity changes are significant.[6]
- The irrotational (core) flow region: The region in which viscous effects and velocity changes are negligible, also known as the inviscid core.[2]
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. But this is not the fully developed fluid flow until the normalized temperature profile also becomes constant.[6]
In case of laminar flow, the velocity profile in the fully developed region is parabolic but in the case of turbulent flow it gets a little flatter due to vigorous mixing in radial direction and eddy motion.
Hydrodynamic Entrance Length -
It is the distance a flow travels in a pipe before it becomes fully developed. Entrance length refers to the length of the entry region, the area following the pipe entrance where effects originating from the interior wall of the pipe propagate into the flow as an expanding boundary layer. When the boundary layer expands to fill the entire pipe, the developing flow becomes a fully developed flow, where flow characteristics no longer change with increased distance along the pipe.
For laminar flow, the hydrodynamic entrance length is given as -
L=0.05ReD or L=1.05m
Analytical Solution by Hagen-Poiseuille equation:
- Average Velocity= `V^(avg)=(Re.mu)/(rho.D)`=0189ms^-1
- Maximum Velocity=`V^(max)=2.v^(avg)=0.3738ms^-1`
- Pressure Drop=`DeltaP=(32.mu.L.v^(avg))/D^2=63.8749 pa`
- Kinematic pressure drop=`DeltaP^k=(Deltap)/rho=0.063874m^2s^-2`
MATLAB code for automatic generation of blockMeshDict file:
%% MATLAB program for automatically generating the blockMeshDict file for a wedge angle provided by the user
clc
clear all
close all
%% Variable values
diameter=0.02;
r=diameter/2;
L=1.5;
theta=input(\'Enter the value of wedge angle = \');
%% Specifying string for header part
h1=\'/*––––––––––––––––––––––––––––––--*- C++ -*–––––––––––––––––––––––––––––––––-*\\ \';
h2=\'| ========= | | \';
h3=\'| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox | \';
h4=\'| \\\\ / O peration | Version: 4.0 | \';
h5=\'| \\\\ / A nd | Web: www.OpenFOAM.org | \';
h6=\'| \\\\/ M anipulation | | \';
h7=\' \\*–––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––*/ \';
h8=\' FoamFile\';
h9=\'{\';
h10=\' version 2.0;\';
h11=\' format ascii;\';
h12=\' class dictionary;\';
h13=\' object blockMeshDict;\';
h14=\'}\';
h15=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';
conv=\'convertToMeters 1;\';
s1=blanks(5);
s2=blanks(10);
s3=blanks(15);
v0=[0 0 0];
v1=[r*sind(theta/2) r*cosd(theta/2) 0];
v2=[-r*sind(theta/2) r*cosd(theta/2) 0];
v3=[0 0 L];
v4=[r*sind(theta/2) r*cosd(theta/2) L];
v5=[-r*sind(theta/2) r*cosd(theta/2) L];
block=\'hex (0 0 1 2 3 3 4 5) (1 50 200) simpleGrading (1 0.5 1)\';
edge1=\'arc 1 2\';
edge2=\'arc 4 5\';
e1=[0 r 0];
e2=[0 r L];
inlet=\'(0 0 2 1)\';
outlet=\'(3 4 5 3)\';
wedge1=\'(0 1 4 3)\';
wedge2=\'(3 5 2 0)\';
wall=\'(4 1 2 5)\';
axis=\'(0 3 3 0)\';
final=\'// ************************************************************************* //\';
%% Writing the blockMeshDict
f1=fopen(\'blockMeshDict.txt\',\'w\');
fprintf(f1,\'%s \\n\',h1);
fprintf(f1,\'%s \\n\',h2);
fprintf(f1,\'%s \\n\',h3);
fprintf(f1,\'%s \\n\',h4);
fprintf(f1,\'%s \\n\',h5);
fprintf(f1,\'%s \\n\',h6);
fprintf(f1,\'%s \\n\',h7);
fprintf(f1,\'%s \\n\',h8);
fprintf(f1,\'%s \\n\',h9);
fprintf(f1,\'%s \\n\',h10);
fprintf(f1,\'%s \\n\',h11);
fprintf(f1,\'%s \\n\',h12);
fprintf(f1,\'%s \\n\',h13);
fprintf(f1,\'%s \\n\',h14);
fprintf(f1,\'%s \\n\',h15);
fprintf(f1,\'\\n\');
fprintf(f1,\'%s \\n\',conv);
fprintf(f1,\'\\n\');
fprintf(f1,\'%s \\n\',\'vertices\');
fprintf(f1,\'%s \\n\',\'(\');
fprintf(f1,\'%s (%f %f %f) \\n\',s1,v0(1),v0(2),v0(3));
fprintf(f1,\'%s (%f %f %f) \\n\',s1,v1(1),v1(2),v1(3));
fprintf(f1,\'%s (%f %f %f) \\n\',s1,v2(1),v2(2),v2(3));
fprintf(f1,\'%s (%f %f %f) \\n\',s1,v3(1),v3(2),v3(3));
fprintf(f1,\'%s (%f %f %f) \\n\',s1,v4(1),v4(2),v4(3));
fprintf(f1,\'%s (%f %f %f) \\n\',s1,v5(1),v5(2),v5(3));
fprintf(f1,\'%s \\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s \\n\',\'blocks\');
fprintf(f1,\'%s \\n\',\'(\');
fprintf(f1,\'%s %s \\n\',s1,block);
fprintf(f1,\'%s \\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s \\n\',\'edges\');
fprintf(f1,\'%s \\n\',\'(\');
fprintf(f1,\'%s %s (%f %f %f) \\n\',s1,edge1,e1(1),e1(2),e1(3));
fprintf(f1,\'%s %s (%f %f %f) \\n\',s1,edge2,e2(1),e2(2),e2(3));
fprintf(f1,\'%s \\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s \\n\',\'boundary\');
fprintf(f1,\'%s \\n\',\'(\');
fprintf(f1,\'%s %s \\n\',s1,\'inflow\');
fprintf(f1,\'%s %s \\n\',s1,\'{\');
fprintf(f1,\'%s %s \\n\',s2,\'type patch;\');
fprintf(f1,\'%s %s \\n\',s2,\'faces\');
fprintf(f1,\'%s %s \\n\',s2,\'(\');
fprintf(f1,\'%s %s \\n\',s3,inlet);
fprintf(f1,\'%s %s \\n\',s2,\');\');
fprintf(f1,\'%s %s \\n\',s1,\'}\');
fprintf(f1,\'%s %s \\n\',s1,\'outflow\');
fprintf(f1,\'%s %s \\n\',s1,\'{\');
fprintf(f1,\'%s %s \\n\',s2,\'type patch;\');
fprintf(f1,\'%s %s \\n\',s2,\'faces\');
fprintf(f1,\'%s %s \\n\',s2,\'(\');
fprintf(f1,\'%s %s \\n\',s3,outlet);
fprintf(f1,\'%s %s \\n\',s2,\');\');
fprintf(f1,\'%s %s \\n\',s1,\'}\');
fprintf(f1,\'%s %s \\n\',s1,\'front\');
fprintf(f1,\'%s %s \\n\',s1,\'{\');
fprintf(f1,\'%s %s \\n\',s2,\'type wedge;\');
fprintf(f1,\'%s %s \\n\',s2,\'faces\');
fprintf(f1,\'%s %s \\n\',s2,\'(\');
fprintf(f1,\'%s %s \\n\',s3,wedge1);
fprintf(f1,\'%s %s \\n\',s2,\');\');
fprintf(f1,\'%s %s \\n\',s1,\'}\');
fprintf(f1,\'%s %s \\n\',s1,\'back\');
fprintf(f1,\'%s %s \\n\',s1,\'{\');
fprintf(f1,\'%s %s \\n\',s2,\'type wedge;\');
fprintf(f1,\'%s %s \\n\',s2,\'faces\');
fprintf(f1,\'%s %s \\n\',s2,\'(\');
fprintf(f1,\'%s %s \\n\',s3,wedge2);
fprintf(f1,\'%s %s \\n\',s2,\');\');
fprintf(f1,\'%s %s \\n\',s1,\'}\');
fprintf(f1,\'%s %s \\n\',s1,\'wall\');
fprintf(f1,\'%s %s \\n\',s1,\'{\');
fprintf(f1,\'%s %s \\n\',s2,\'type wall;\');
fprintf(f1,\'%s %s \\n\',s2,\'faces\');
fprintf(f1,\'%s %s \\n\',s2,\'(\');
fprintf(f1,\'%s %s \\n\',s3,wall);
fprintf(f1,\'%s %s \\n\',s2,\');\');
fprintf(f1,\'%s %s \\n\',s1,\'}\');
fprintf(f1,\'%s %s \\n\',s1,\'axis\');
fprintf(f1,\'%s %s \\n\',s1,\'{\');
fprintf(f1,\'%s %s \\n\',s2,\'type empty;\');
fprintf(f1,\'%s %s \\n\',s2,\'faces\');
fprintf(f1,\'%s %s \\n\',s2,\'(\');
fprintf(f1,\'%s %s \\n\',s3,axis);
fprintf(f1,\'%s %s \\n\',s2,\');\');
fprintf(f1,\'%s %s \\n\',s1,\'}\');
fprintf(f1,\'%s \\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s \\n\',\'mergePatchPairs\');
fprintf(f1,\'%s \\n\',\'(\');
fprintf(f1,\'%s \\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s \\n\',final);
fclose(f1);OpenFOAM case setup:
A Wedge boundary condition is used in this case as it is a cylinder and flow is axi-symmetric along the axis.
After updating the blockMeshDict and pressure and velocity boundary conditions file, simulation is ran using icoFoam solver which is used to simulate turbulent, incompressible flows. The simulation parameters used are
- `Deltat=1e-5`
- write interval = 10 step
blockMeshDict file generated for `2^(0)`
/*––––––––––––––––––––––––––––––--*- C++ -*–––––––––––––––––––––––––––––––––-*\\
| ========= | |
| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\\\ / O peration | Version: 4.0 |
| \\\\ / A nd | Web: www.OpenFOAM.org |
| \\\\/ M anipulation | |
\\*–––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000175 0.009998 0.000000)
(-0.000175 0.009998 0.000000)
(0.000000 0.000000 1.500000)
(0.000175 0.009998 1.500000)
(-0.000175 0.009998 1.500000)
);
blocks
(
hex (0 0 1 2 3 3 4 5) (1 50 200) simpleGrading (1 0.5 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (0.000000 0.010000 1.500000)
);
boundary
(
inflow
{
type patch;
faces
(
(0 0 2 1)
);
}
outflow
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type wedge;
faces
(
(0 1 4 3)
);
}
back
{
type wedge;
faces
(
(3 5 2 0)
);
}
wall
{
type wall;
faces
(
(4 1 2 5)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
Velocity Boundary Condition
/*--------------------------------*- C++ -*----------------------------------*\\
| ========= | |
| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\\\ / O peration | Version: 4.1 |
| \\\\ / A nd | Web: www.OpenFOAM.org |
| \\\\/ 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
{
inflow
{
type fixedValue;
value uniform (0 0 0.1869);
}
outflow
{
type zeroGradient;
}
wall
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //Pressure Boundary Condition
/*--------------------------------*- C++ -*----------------------------------*\\
| ========= | |
| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\\\ / O peration | Version: 4.1 |
| \\\\ / A nd | Web: www.OpenFOAM.org |
| \\\\/ M anipulation | |
\\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inflow
{
type zeroGradient;
}
outflow
{
type fixedValue;
value uniform 0.063874;
}
wall
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //2 degree
BlockMeshDict file Generation for 3 degree
vertices
(
(0.000000 0.000000 0.000000)
(0.000262 0.009997 0.000000)
(-0.000262 0.009997 0.000000)
(0.000000 0.000000 1.500000)
(0.000262 0.009997 1.500000)
(-0.000262 0.009997 1.500000)
); BlockMeshDict file Generation for 4 degree
vertices
(
(0.000000 0.000000 0.000000)
(0.000349 0.009994 0.000000)
(-0.000349 0.009994 0.000000)
(0.000000 0.000000 1.500000)
(0.000349 0.009994 1.500000)
(-0.000349 0.009994 1.500000)
);
At x=0.05 m
The velocity profile near the entrance is flat at the center of the pipe. It means the flow is still developing.
A t x=0.25 m
At x=0.5 m
At x=0.75 m
At x=1.15 m (Fully developed velocity profile)
A Parabolic velocity profile is obtained after a certain distance in the pipe. This profile remains unchanged thereafter. It concludes that the length of the pipe as calculated theoretically is sufficient to get a fully developed flow at the outlet of the pipe.
At x=1.15 m shear stress
The observed shear stress profile agrees with the theoretical plots found in the literature. Stress is maximum at the walls and zero at the center of the pipe.
Conclusion:
Comparing the maximum velocities obtained by each simulation with the theoretical maximum velocity, the simulation with a wedge angle of 4 degrees gives less error.