Simulation of flow through a pipe using Open Foam
OBJECTIVE
- Write a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes.
- To simulate for internal laminar flow through the pipe and compute the velocity profiles at various sections and set up the sufficient length of the pipe so that the flow is fully developed.
- To correlate the velocity profile obtained from the simulations with the Hagen poiseuille\'s equation.
PROBLEM DESCRIPTION
- The wedge shape is formed due to the symmetry of the pipe to reduce the simulation time.
- The length of the pipe using the entry length formula for laminar flow through a pipe.
- The simulations are performed in OpenFOAM with the icoFoam solver.
- The assumptions for the boundary condition are made depending on the Hagen poiseuille\'s equation.
Assumptions
Reynolds number - 2100.
Working fluid - water.
Temperature - 25c.
Density - 997 kg/m^3.
Diameter - 0.01m.
Wedge angle - 3.
Hagen Poiseuille
It 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 the constant cross-section. Equations for analyzing pipe flow, such as the Darcy Weisbach equation for frictional head loss, often apply only to the fully developed flow portion of the pipe flow. If the total pipe length is large compared to the entrance length, then the effect of the entrance length can usually be neglected and the total pipe length can be used in calculations. If the total pipe length is relatively short in comparison with the entrance length, however, then the entrance region may need to be analyzed separately. An estimate of the entrance length is sometimes needed to determine how to proceed with pipe flow calculations. The Reynolds number for pipe flow is needed to calculate the entrance length for turbulent flow or laminar flow.
The diagram at the left illustrates the meanings of \"entrance region\" and \"fully developed flow.\" When fluid enters a pipe its velocity will often be uniform across the pipe cross-section as shown in the diagram. Near the entrance, the fluid in the center of the pipe isn’t affected by the friction between fluid and pipe walls, but as the flow proceeds down the pipe, the effect of the wall friction moves in toward the pipe center, until the pattern of velocity variation across the pipe (called the velocity profile) becomes constant. The entrance portion of the pipe, where the velocity profile is changing is called the entrance region, and the flow after that entrance region is called \"fully developed flow.\" The next two sections will present equations for estimating the length of the entrance region, called the entrance length, for pipe flow.
Flow-through a Circular pipes
Re = D*v*ρ/μ
Now we have a formula for the velocity of liquid moving through the tube as a function of the distance from the center of the tube
![\"{\\displaystyle](\"https://wikimedia.org/api/rest_v1/media/math/render/svg/97ef52092bc6cc3c1450f232d9b6ae0b23fe6a2c\")
At the center of the tube where the liquid is moving fastest (r = 0) with R being the radius of the tube,
![\"{\\displaystyle](\"https://wikimedia.org/api/rest_v1/media/math/render/svg/e7c5335be75ba8e5fec8fe040d9438b04cc4322a\")
D = pipe diameter.
ρ = fluid density.
μ = fluid viscosity.
V = velocity.
Re = Reynold\'s number.
For Laminar Flow, the entrance length, Le, can be estimated from the equation:
Le/D = 0.06 Re.
MAIN PROGRAM
The circular pipe is symmetry so the small wedge angle is taken and the MATLAB is used to create to blockmeshdict file for the Open foam. The code is for the 3-degree wedge angle
clear all
close all
clc
% Input Paramaters
Re = 2100; % given value
% assume values
t = 25;
rho = 997;
nu = 8.9e-4;
d = 0.01;
r = d/2;
theta =3; %wedge angle
l = 0.06*d*Re; % entrance length formula for laminar flow
x = l;
y = r*cosd(theta/2);
z = r*sind(theta/2);
v = Re*nu/(rho*d); % average velocity
v_max = 2*v; % maximum velocity
p = 32*nu*v*l/d^2; % pressure
p_k = p/rho; % pressure value is in terms of kinetic pressure so it is divided by density
t1 = (\'/*------------------------------------ *- C++ -*----------------------------------------*\\\');
t2 = (\' ============ |\');
t3 = (\' \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox\');
t4 = (\' \\\\ / O peration | Website: https://openfoam.org\');
t5 = (\' \\\\ / A nd | Version: 6\');
t6 = (\' \\\\/ M anipulation |\');
% write a file
fid = fopen(\'blockmesh\',\'w+t\')
fprintf(fid,\'%s \\n\',t1)
fprintf(fid,\'%s\\n\',t2)
fprintf(fid,\'%s \\n\',t3)
fprintf(fid,\'%s \\n\',t4)
fprintf(fid,\'%s \\n\',t5)
fprintf(fid,\'%s \\n\',t6)
fprintf(fid,\'/*--------------------------------------------------------------------------------------*/ \\n\')
fprintf(fid,\' FoamFile \\n\')
fprintf(fid,\' { \\n\')
fprintf(fid,\' version 2.0; \\n\')
fprintf(fid,\' format ascii; \\n\')
fprintf(fid,\' class dictionary; \\n\')
fprintf(fid,\' object blockMeshDict; \\n\')
fprintf(fid,\' } \\n\')
fprintf(fid,\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // \\n\\n\')
fprintf(fid,\'convertToMeters 1; \\n\\n\')
fprintf(fid,\'vertices \\n\')
fprintf(fid,\'( \\n\')
fprintf(fid,\'(0 0 0) //point 0 \\n\')
fprintf(fid,\'(0 %f %f) //point 1 \\n\',y,z)
fprintf(fid,\'(0 %f %f) //point 2 \\n\',y,-z)
fprintf(fid,\'(%f 0 0) //point 3 \\n\',x)
fprintf(fid,\'(%f %f %f) //point 4 \\n\',x,y,z)
fprintf(fid,\'(%f %f %f) //point 5 \\n\',x,y,-z)
fprintf(fid,\'); \\n\\n\')
fprintf(fid,\'blocks \\n\')
fprintf(fid,\'( \\n\')
fprintf(fid,\'hex (0 3 5 2 0 3 4 1) (100 50 100) simplegrading (1 1 1) \\n\')
fprintf(fid,\'); \\n\\n\')
fprintf(fid,\'edges \\n\')
fprintf(fid,\'( \\n\')
fprintf(fid,\' arc 1 2 (0 %f 0)\\n\',r)
fprintf(fid,\' arc 4 5 (%f %f 0)\\n\',l,r)
fprintf(fid,\'); \\n\\n\')
fprintf(fid,\'boundary \\n\')
fprintf(fid,\'( \\n inlet \\n \')
fprintf(fid,\' { \\n\')
fprintf(fid,\'\\n type patch; \\n\')
fprintf(fid,\' faces \\n ( \\n\')
fprintf(fid,\' (0 1 2 0)\\n\')
fprintf(fid,\' ); \\n\')
fprintf(fid,\' } \\n \')
fprintf(fid,\'outlet \\n \')
fprintf(fid,\' { \\n\')
fprintf(fid,\'\\n type patch; \\n\')
fprintf(fid,\' faces \\n ( \\n\')
fprintf(fid,\' (3 5 4 3)\\n\')
fprintf(fid,\' ); \\n\')
fprintf(fid,\' } \\n \')
fprintf(fid,\'front\\n \')
fprintf(fid,\' { \\n\')
fprintf(fid,\'\\n type wedge; \\n\')
fprintf(fid,\' faces \\n ( \\n\')
fprintf(fid,\' (1 0 3 4)\\n\')
fprintf(fid,\' ); \\n\')
fprintf(fid,\' } \\n \')
fprintf(fid,\'back\\n \')
fprintf(fid,\' { \\n\')
fprintf(fid,\'\\n type wedge; \\n\')
fprintf(fid,\' faces \\n ( \\n\')
fprintf(fid,\' (2 5 3 0)\\n\')
fprintf(fid,\' ); \\n\')
fprintf(fid,\' } \\n \')
fprintf(fid,\'top \\n \')
fprintf(fid,\' { \\n\')
fprintf(fid,\'\\n type wall; \\n\')
fprintf(fid,\' faces \\n ( \\n\')
fprintf(fid,\' (2 1 4 5)\\n\')
fprintf(fid,\' ); \\n\')
fprintf(fid,\' } \\n \')
fprintf(fid,\'bottom \\n \')
fprintf(fid,\' { \\n\')
fprintf(fid,\'\\n type empty; \\n\')
fprintf(fid,\' faces \\n ( \\n\')
fprintf(fid,\' (0 3 3 0)\\n\')
fprintf(fid,\' ); \\n\')
fprintf(fid,\' } \\n \')
fprintf(fid,\'); \\n\')
fprintf(fid,\'mergePatchPairs \\n\')
fprintf(fid,\'( \\n\')
fprintf(fid,\'); \\n\')
fprintf(fid,\'// **************************************************************************** //\\n\')
% close the file
fclose(fid)Block Mesh File
The block mesh file for the wedge angle 3
/*------------------------------------ *- C++ -*----------------------------------------*\\
============ |
\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\\\ / O peration | Website: https://openfoam.org
\\\\ / A nd | Version: 6
\\\\/ M anipulation |
/*--------------------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0) //point 0
(0 0.004619 0.001913) //point 1
(0 0.004619 -0.001913) //point 2
(1.260000 0 0) //point 3
(1.260000 0.004619 0.001913) //point 4
(1.260000 0.004619 -0.001913) //point 5
);
blocks
(
hex (0 3 5 2 0 3 4 1) (1 1 1) simplegrading (0.2 0.2 1)
);
edges
(
arc 1 2 (0 0.005000 0)
arc 4 5 (1.260000 0.005000 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
front
{
type wedge;
faces
(
(1 0 3 4)
);
}
back
{
type wedge;
faces
(
(2 5 3 0)
);
}
top
{
type wall;
faces
(
(2 1 4 5)
);
}
bottom
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// **************************************************************************** //Control Dict File
/*--------------------------------*- C++ -*----------------------------------*\\
========= |
\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\\\ / O peration | Website: https://openfoam.org
\\\\ / A nd | Version: 6
\\\\/ 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;
// ************************************************************************* //Initial Conditions: Inlet Velocity
/*--------------------------------*- C++ -*----------------------------------*\\
========= |
\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\\\ / O peration | Website: https://openfoam.org
\\\\ / A nd | Version: 6
\\\\/ 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.375 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
bottom
{
type empty;
}
}
// ************************************************************************* //Initial Conditions: Outlet Pressure
/*--------------------------------*- C++ -*----------------------------------*\\
========= |
\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\\\ / O peration | Website: https://openfoam.org
\\\\ / A nd | Version: 6
\\\\/ 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.0675;
}
front
{
type wedge;
}
back
{
type wedge;
}
top
{
type zeroGradient;
}
bottom
{
type empty;
}
}
// ************************************************************************* //RESULTS
The Geometry of the wedge at 3 degrees
Velocity Profiles
Profile value @ 0.1 distance
Profile value @ 0.5 distance
Profile value @ 1 distance
Profile value @ 1.26 distance
INFERENCE
The velocity profile becomes fully developed and follows Hagen Poiseuille law. The entrance length calculated using the laminar flow equation along the x-axis was found to be 1.26. The plot for the velocity less than the entrance length is not uniform it varies along the length which shows the significance of the entrance length in the pipe flow. The thickness of the boundary layer increases in the flow direction until the boundary layer reaches the pipe center and thus fills the entire pipe. The region beyond the entrance region in which the velocity profile is fully developed and remains unchanged. From the plot, it is clear that at the hydrodynamic entry length the value of the velocity is almost constant.
The value of the velocity from the Hagen poiseuille equation is 0.375 m/s and the value of the velocity obtained from the simulations in the Open Foam is 0.3749 m/s. The mean and the maximum courant number from the simulations is 0.0297 & 0.0741.