Week 11 - Simulation of Flow through a pipe in OpenFoam
Project Aim: Simulation of flow through a pipe in OpenFoam
Project Objectives:
- To generated automated mesh using MATLAB program
- Make sure that number of cells along the "axi symmetric"
- For a baseline mesh, show that the velocity profile matches with the Hagen Poiseuille's equation
- Show that the velocity profile is fully developed
- Post process velocity and shear stress
To simulate the flow of water through a section of a pipe and to validate the CFD results using theoritical values from Hagen Poiseuille's equations.
Theory:
The understanding of the theory would help in assuming the right inputs to the pipe flow simulation and successful validation of the CFD results with the theoritical solutions. Let us try to undertand Hagen Poiseuille's flow.
When the fluid flows through a pipe, its velocity changes as a function of 'x' and 'r'. Mathematically
`u=u(r,x)`
where
`r implies` radial distance from the center line of the pipe
`x implies` distance along the length of the pipe from inlet
The fluid velocity in a pipe changes from zero at the surface because of the no-slip condition to a maximum at the pipe center. In fluid flow, it is convenient to work with an average velocity `V_(avg`, which remains constant in incompressible flow when the cross-sectional area of the pipe is constant.
From the conservation of mass priciple,
when we know the flow rate or the velocity profile, the average velocity can be determined easily. Mathematically, the above equation for the circular pipe reduces to,
Reynold's Number for the pipe flow:
The Entrance Region:
The region from the pipe inlet to the point at which the boundary layer merges at the centerline is called the hydrodynamic entrance region, and the length of this region is called the hydrodynamic entry length `L_h`.
Flow in the entrance region is called hydrodynamically developing flow since this is the region where the velocity profile develops.
The region beyond the entrance region in which the velocity profile is fully developed and remains unchanged is called the hydrodynamically fully developed region. The flow is said to be fully delveloped when the normalized temperature profile remains unchanged as well. Hydrodynamically developed flow is equivalent to fully developed flow when the fluid in the pipe is not heated or cooled since the fluid temperature in this case remains essentially constant throughout.
The velocity profile in the fully developed region is parabolic in laminar flow and somewhat flatter (or fuller) in turbulent flow due to eddy motion and more vigorous mixing in the radial direction. The time-averaged velocity profile remains unchanged when the flowis fully developed, and thus
The shear stress at the pipe wall tw is related to the slope of the velocity profile at the surface. Noting that the velocity profile remains unchanged in the hydrodynamically fully developed region, the wall shear stress also remains constant in that region.
The Entry Length (`L_h`):
`L_(h,laminar) approx 0.05ReD`
Equations to be used for simulating Hagen Poiseuille flow:
1) Velocity at any location in fully developed region at any x:
2) Avg. Velocity:
3) Velocity at centre line:
4) Pressure drop (Laminar):
5) Pressure drop (Laminar/ Turbulent) in terms of friction factor:
6) friction factor:
7) Volume flow rate:
This equation is known as Poiseuille’s law, and this flow is called Hagen–Poiseuille flow in honor of the works of G. Hagen (1797–1884) and J. Poiseuille (1799–1869) on the subject. Note from above equation that for a specified flow rate, the pressure drop and thus the required pumping power is proportional to the length of the pipe and the viscosity of the fluid, but it is inversely proportional to the fourth power of the radius (or diameter) of the pipe. Therefore, the pumping power requirement for a piping system can be reduced by a factor of 16 by doubling the pipe diameter. Of course the benefits of the reduction in the energy costs must be weighed against the increased cost of construction due to using a larger-diameter pipe.
Assumed Parameters for Simulation:
1) Geometry Consideration:
Flow can nit be simulated in the entire cylinder with the limited computer resrources, so only a sector of the circular pipe will be simulated and the geometry will be axi-symmetric.
2) Wedge angle (`theta`) = 4 dregrees
3) Reynold's Number (`Re`) = 2000 (Laminar)
4) Working fluid = water (`rho=997(kg)/m^3`)
5) Temperature of working fluid = `25 ^0C`
6) Dynamic Viscosity of working fluid (`mu`) = `8.9 xx10^-4 Pa.sec`
7) Kinematic viscosity of working fluid (`nu`) = `8.926xx10^-7 m^2/sec`
8) Diameter of the pipe (D) = 0.02m
Analytical Solution:
1) Average Velocity:
`Re=(V_(avg)xxDxxrho)/mu`
`implies V_(avg)=(muxxRe)/(rhoxxD)`
`implies V_(avg)=(8.9xx10^(-4)xx2000)/(997xx0.02)`
`implies V_(avg)=0.089m/sec`
2) Maximum Velocity:
`V_(max)=2xxV_(avg)=2xx0.089=0.178m/sec`
3) Hydrodynamic Entry Length (`L_(h,laminar)`):
`L_h=0.05xxRexxD=0.05xx2000xx0.02=2m`
So, in order to have a fully developed flow `L> >L_h`.
4) So, Assuming `L=3m`
5) Pressure Drop:
`DeltaP=(32muLV_(avg))/D^2=(32xx8.9xx10^-4xx3xx0.089)/0.02^2=19Pa`
6) Kinematic Pressure Drop:
`P_k=(DeltaP)/rho=(19)/997=0.0191m`
7) Shear Stress:
`f=8tau_w/(rhoV_(avg)^2)=(64)/(Re)=64/2000=0.032`
`implies tau_w=(0.032xx997xx0.089^2)/8=0.0316Pa`
MATLAB Code for automatic mesh generation:
% MATLAB script to generate automated mesh for any wedge angle
clear all
close all
clc
% given inputs
% wedge angle
theta = 4; % Wedge angle
Re = 2000; % Reynold's Number
D = 0.02; % Diameter
R = D/2; % Radius
L_h = 0.05*Re*D; % Hydrodynamic entry length
L = 1.5*L_h; % Length of the wedge
mu = 8.9e-4; % Dynamic Viscosity
rho = 997; % Density
nu = mu/rho; % Kinematic Vicosity
% Analytical Solution
v_avg = (mu*Re)/(D*rho); % Average Velocity
v_max = 2*v_avg; % Maximum velocity
del_p = (32*mu*L*v_avg)/(D^2); % Pressure Drop in Pa
kin_p = del_p/rho; % Kinematic Pressure
% Vertices of the wedge
v_0 = [0 0 0];
v_1 = [L 0 0];
v_2 = [L R*cosd(theta/2) R*sind(theta/2)];
v_3 = [0 R*cosd(theta/2) R*sind(theta/2)];
v_4 = [0 R*cosd(theta/2) -R*sind(theta/2)];
v_5 = [L R*cosd(theta/2) -R*sind(theta/2)];
%Printing the file
text1 = '/*--------------------------------*- C++ -*----------------------------------*\';
text2 = '========= |';
text3 = '\\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
text4 = ' \\ / O peration | Website: https://openfoam.org';
text5 = ' \\ / A nd | Version: 7';
text6 = ' \\/ M anipulation |';
text7 = '\*---------------------------------------------------------------------------*/';
text8 = 'FoamFile';
text9 = '{';
text10 = 'version 2.0;';
text11 = 'format ascii;';
text12 = 'class dictionary;';
text13 = 'object blockMeshDict;';
text14 = '}';
text15 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
% Creating file
f1 = fopen('mesh3','w');
fprintf(f1,'%s\n',text1);
fprintf(f1,'%s\n',text2);
fprintf(f1,'%s\n',text3);
fprintf(f1,'%s\n',text4);
fprintf(f1,'%s\n',text5);
fprintf(f1,'%s\n',text6);
fprintf(f1,'%s\n',text7);
fprintf(f1,'%s\n',text8);
fprintf(f1,'%s\n',text9);
fprintf(f1,'%s\n',text10);
fprintf(f1,'%s\n',text11);
fprintf(f1,'%s\n',text12);
fprintf(f1,'%s\n',text13);
fprintf(f1,'%s\n',text14);
fprintf(f1,'%s\n',text15);
fprintf(f1,'%s\n','convertToMeters 1;');
fprintf(f1,'\n');
fprintf(f1,'%s\n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s(%f%s%f%s%f)\n',blanks(5),v_0(1),blanks(1),v_0(2),blanks(1),v_0(3));
fprintf(f1,'%s(%f%s%f%s%f)\n',blanks(5),v_1(1),blanks(1),v_1(2),blanks(1),v_1(3));
fprintf(f1,'%s(%f%s%f%s%f)\n',blanks(5),v_2(1),blanks(1),v_2(2),blanks(1),v_2(3));
fprintf(f1,'%s(%f%s%f%s%f)\n',blanks(5),v_3(1),blanks(1),v_3(2),blanks(1),v_3(3));
fprintf(f1,'%s(%f%s%f%s%f)\n',blanks(5),v_4(1),blanks(1),v_4(2),blanks(1),v_4(3));
fprintf(f1,'%s(%f%s%f%s%f)\n',blanks(5),v_5(1),blanks(1),v_5(2),blanks(1),v_5(3));
fprintf(f1,'%s\n',');');
fprintf(f1,'\n');
fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s %s \n',blanks(5),'hex (0 3 2 1 0 4 5 1) (30 100 1) simpleGrading (1 1 1)');
fprintf(f1,'%s\n',');');
fprintf(f1,'\n');
fprintf(f1,'%s\n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s arc 4 3 (%f %f %f)\n',blanks(5),0,R,0);
fprintf(f1,'%s arc 5 2 (%f %f %f)\n',blanks(5),L,R,0);
fprintf(f1,'%s\n',');');
fprintf(f1,'\n');
fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s %s\n',blanks(5),'inlet');
fprintf(f1,'%s %s\n',blanks(5),'{');
fprintf(f1,'%s %s\n',blanks(10),'type patch;');
fprintf(f1,'%s %s\n',blanks(10),'faces');
fprintf(f1,'%s %s\n',blanks(10),'(');
fprintf(f1,'%s %s\n',blanks(15),'(0 3 4 0)');
fprintf(f1,'%s %s\n',blanks(10),');');
fprintf(f1,'%s %s\n',blanks(5),'}');
fprintf(f1,'%s %s\n',blanks(5),'outlet');
fprintf(f1,'%s %s\n',blanks(5),'{');
fprintf(f1,'%s %s\n',blanks(10),'type patch;');
fprintf(f1,'%s %s\n',blanks(10),'faces');
fprintf(f1,'%s %s\n',blanks(10),'(');
fprintf(f1,'%s %s\n',blanks(15),'(1 2 5 1)');
fprintf(f1,'%s %s\n',blanks(10),');');
fprintf(f1,'%s %s\n',blanks(5),'}');
fprintf(f1,'%s %s\n',blanks(5),'front');
fprintf(f1,'%s %s\n',blanks(5),'{');
fprintf(f1,'%s %s\n',blanks(10),'type wedge;');
fprintf(f1,'%s %s\n',blanks(10),'faces');
fprintf(f1,'%s %s\n',blanks(10),'(');
fprintf(f1,'%s %s\n',blanks(15),'(0 3 2 1)');
fprintf(f1,'%s %s\n',blanks(10),');');
fprintf(f1,'%s %s\n',blanks(5),'}');
fprintf(f1,'%s %s\n',blanks(5),'back');
fprintf(f1,'%s %s\n',blanks(5),'{');
fprintf(f1,'%s %s\n',blanks(10),'type wedge;');
fprintf(f1,'%s %s\n',blanks(10),'faces');
fprintf(f1,'%s %s\n',blanks(10),'(');
fprintf(f1,'%s %s\n',blanks(15),'(0 4 5 1)');
fprintf(f1,'%s %s\n',blanks(10),');');
fprintf(f1,'%s %s\n',blanks(5),'}');
fprintf(f1,'%s %s\n',blanks(5),'pipe');
fprintf(f1,'%s %s\n',blanks(5),'{');
fprintf(f1,'%s %s\n',blanks(10),'type wall;');
fprintf(f1,'%s %s\n',blanks(10),'faces');
fprintf(f1,'%s %s\n',blanks(10),'(');
fprintf(f1,'%s %s\n',blanks(15),'(3 2 5 4)');
fprintf(f1,'%s %s\n',blanks(10),');');
fprintf(f1,'%s %s\n',blanks(5),'}');
fprintf(f1,'%s %s\n',blanks(5),'axis');
fprintf(f1,'%s %s\n',blanks(5),'{');
fprintf(f1,'%s %s\n',blanks(10),'type empty;');
fprintf(f1,'%s %s\n',blanks(10),'faces');
fprintf(f1,'%s %s\n',blanks(10),'(');
fprintf(f1,'%s %s\n',blanks(15),'(0 1 1 0)');
fprintf(f1,'%s %s\n',blanks(10),');');
fprintf(f1,'%s %s\n',blanks(5),'}');
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n',');');
fprintf(f1,'\n')
fprintf(f1,'%s\n','// ************************************************************************* //');
BlockMeshDict File generated:
/*--------------------------------*- 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(3 0.000000 0.000000)
(3 0.009994 0.000349)
(0.000000 0.009994 0.000349)
(0.000000 0.009994 -0.000349)
(3 0.009994 -0.000349)
);
blocks
(
hex (0 3 2 1 0 4 5 1) (30 100 1) simpleGrading (1 1 1)
);
edges
(
arc 4 3 (0.000000 0.010000 0.000000)
arc 5 2 (10.000000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 3 4 0)
);
}
outlet
{
type patch;
faces
(
(1 2 5 1)
);
}
front
{
type wedge;
faces
(
(0 3 2 1)
);
}
back
{
type wedge;
faces
(
(0 4 5 1)
);
}
pipe
{
type wall;
faces
(
(3 2 5 4)
);
}
axis
{
type empty;
faces
(
(0 1 1 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
Initial and Boundary Conditions:
For Velocity:
/*--------------------------------*- 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.089 0 0);
}
outlet
{
type zeroGradient;
}
pipe
{
type noSlip;
}
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: 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.0191;
}
pipe
{
type zeroGradient;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
Kinematic Viscosity:
/*--------------------------------*- 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] 0.0000008926;
// ************************************************************************* //
Control Dict:
/*--------------------------------*- 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 10;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Generated mesh:
30 cells were used along radial direction and 100 cells were used along the length of the wedge. No scaling factor was used to generate the mesh.
3D view of the generated wedge with wedge angle 4 degrees
Results:
Velocity line plots:
For visualizing the velocity profile at different sections of the wedge so as to see whether the flow is fully developed or not, folllowing line plots were plotted
Velocity is plotted on the same graph at distances 0, 0.5, 1, 1.5 and 2m respectively. Theoritically, as shown previously, as per Hagen Poissuille's equations flow is hydrodynamically undeveloped in this region.
- At x = 0m, the velocity is constant through out and its value is same as that of the initial value (0.089 m/s) through out the radius.
- At x = 0.5m, the maximum velocity is around 14.2 m/s, which drops down parabolically to zero at the walls (R = 0.01 m).
- At x = 1 m and x = 1.5 m, flow is still not developed completely, as slight variations in the profiles can be seen.
- At x = 2 , 2.25, 2.5, 2.75 and 3 m, there is hardly any difference in velocity profiles. So, beyond, 2 m , velocity profile is constant with respection to X. This is what was expected as per the theoritical results. Although the maximum velocity at the radial centre of the pipe is not the same as the theoritical value (comparison will be given next).
This was reason to take the wedge length as 3 m, so that the flow analysed must be fully developed to compare with theoritical results.
Shear Stress line plots:
The shear stress plot above is linear (atleast in the region 0.003 - 0.009), and at the boundarys (centre and walls), it has some non-linearity. One of the reasons for this behaviour might be coarser grid and no scaling factor in radial direction. This linear nature shows that at x = 2.5 m the flow is fully developed because mathematically, velocity profile is parabolic and shear stress is propertional to velocity gradient in 'Y' direction. Therefore this proves that the flow is fully developed.
Also, the wall shear stress obtained, shows that the steady value of 0.0297 Pa is obtained near the wall so as to prove that the simulation results are quite satisfying the Hagen Poiseuille's equations. Also the trend shown below is in accordance with the distribution shown while explaining the theory in start of the report.
Validation of the simulation results with Hagen Poiseuille's flow:
For Velocity:
For Wall Shear Stress:
Conclusions:
- The simulation results were in accordance with the Values obtained from the Hagen Poiseuille's Equations with an acceptable degree of error.
- These errors obtained may be reduced further by finer grids but the computation time would increase. This is the reason, many different grids were not simulated.
- The selection of the length of the pipe can be made by identifying the development of the velocity profile also the material selection for the pipe is made by considering the wall shear stress value.