Pipe Flow Simulation for various wedge angles OPENFOAM part1
[Part 1/2]:
[1]INTRODUCTION
In this project we aim to simulate a section of pipe of Wedge shaped geometry to understand various aspects of flow through a pipe.The fluid chosen is water. The flow is laminar and incompressible. The specified Reynolds No. is 2100
Internal flow occurs when the fluid is filled entirely in a conduit and is driven by a difference in pressure.
Fluid velocity in the pipe changes from zero at the surface because of the no-slip boundary condition to a maximum at the pipe centre.
The average velocity `Vavg` remains constant in an incompressible flow when the cross-section of the pipe is constant.
Reynolds Number Re: The transition from Laminar to turbulent flow occurs when Laminar(Re`<=`2100) to Transition(`2100<=Re<=4000`) to turbulent(`Re>=4000`) for a smooth pipe and it has no viscous dissipation.It depends on the geometry, face roughness, flow velocity, surface temperature etc. Reynolds No. is the ratio of Inertial Force to viscous force.
`Re = (Vavg*rho*D)/(mu`
Hydrodynamic Entry length: The region from the pipe inlet to the point at which the velocity boundary layer merges at the centreline is called the hydrodynamic entrance region. The region beyond the entrance region in which the velocity profile is fully developed is the fully debveloped region.
L(laminar) = `0.06*Re*D`
Velocity profile for laminar flow: the velocity profile for a fully developed laminar flow in a pipe is parabolic with a maximum at the centreline and minimum at R=0
`u(r)= (R^2)/(4mu)*((dp)/(dx))*(1-(r^2)/(R^2))`.
The average velocity = `(R^2)/(8mu)*((dp)/(dx))`.Hence the final velocity profile can be expressed as `u(r)=2*Vavg*(1-(r^2)/(R^2))`
Hagen poiseuille Flow: It states that \" for a specified Flow rate, the pressure drop and the pumping power is inversely proportional to the length of the pipe and the viscosity of the fluid, but is inversely proportional to the fourth power of the radius(or diameter) of the pipe.
Volumetric Flow Rate, `Q =Vavg*Ac = (DeltaP*pi*D^4)/(128*mu*L)`
Thus the average velocity is Vavg =`(DeltaP*D^2)/(32*mu*L)`
[II]SPECIFICATIONS OF THE PROBLEM:
1. Given reynolds Number = 2100
2. The working fluid is water at 25`@C`
3. The entry length, L(entrance)= 0.06*Re*D =2.52 m
4. Total length of the pipe taken for the simulation is 2.82 m
5. Density of the fluid, water = `rho = 997.0 kgm^(-3)`
6. The dynamic viscosity of water is `mu = `0.89*10^-3`Pa.S`
7. The kinematic viscosity of water is `8.927*10^-7` `m^2s^-1`
CALCULATED VALUES: (Analytical Solution)
1. Average velocity Vavg= `(Re)/(pho)*(mu)/(D)` = 0.0937 `ms^-1`
2. Maximum velocity, Vmax = 2. Vavg = `0.187 ms^-1`
3. Pressure Drop `DeltaP = (32*(mu)*L*Vavg)/(D^2)` = 18.8197 Pa
4. Kinematic Pressure Drop `DeltaP_(kinematic)` =`(DeltaP)/rho`= 0.0189 m
[III]Matlab Program:
In this project we have used Matlab for the generation of BlockMesh Dict file. For different wedge angles Matlab generates the required file.
Wedge angle 2, 3 and 5 are used here.
clear all
close all
clc
%% Input parameters
Re = 2100 ; %given Reynold\'s No.
d =0.02; % Diameter of the pipe
r = d/2 ; % Radius of the pipe
L_entrance = 0.06*Re*d; % Length of the pipe
fprintf(\'Enter 2,3,4,5 for wedge boundary conditions nEnter 10, 25, 45 for symmetry boundary conditions\');
theta = input(\'enter the angle in angles:n\');% angle of wedge
mu = 0.89e-3; %Viscosity of water at 25 degree C
rho = 997.0; % Density of water at 25 degree C
nu = mu/rho; % Kinematic viscosity of water at 25 degree C
nx = 400;
ny = 40;
nz = 1;
%% Calculation of Pressure and velocity
L= L_entrance + 0.3; % length of the wedge
v_avg = Re*mu/(d*rho); % Average Velocity
v_max = 2*v_avg; % max velocity
del_P = (32*mu*v_avg*L)/(d^2); % Pressure Drop
kin_P = del_P/rho; % Kinematic Pressure Drop
%% calculating the 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)];
%% String Format specifying the BlockMesh file
head1 = \"/*--------------------------------*- C++ -*----------------------------------*\";
head2 = \"========= | \";
head3 = \"\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \";
head4 = \" \\\\ / O peration | Website: https://openfoam.org \";
head5 = \" \\\\ / A nd | Version: 6 \";
head6 = \" \\\\/ M anipulation | \";
head7 = \"\\*---------------------------------------------------------------------------*/\";
b4 = blanks(4);
b8 = blanks(8);
b12 = blanks(12);
b16 = blanks(16);
dict1 = \"version 2.0;\";
dict2 = \"format ascii;\";
dict3 = \"class dictionary;\";
dict4 = \"object blockMeshDict;\";
dict5 = \"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //;\";
conv = \"convertToMeters 1;\";
vertex1 = \"vertices\";
%blocks1 = b4 + \"hex (0 1 2 0 3 4 5 3)(50 1 200) simpleGrading (0.2 1 1)\";
blocks1 = b4 + \"hex (0 3 5 2 0 3 4 1)(400 40 1) simpleGrading (1 0.1 1)\";
footer = \"// ************************************************************************* //\";
%% creating BlockMeshdict file for openFoam
f1 = fopen(\'blockMeshdict_2.txt\',\'wt\');
fprintf(f1,\'%s\\n\', head1);
fprintf(f1,\'%s\\n\', head2);
fprintf(f1,\'%s\\n\', head3);
fprintf(f1,\'%s\\n\', head4);
fprintf(f1,\'%s\\n\', head5);
fprintf(f1,\'%s\\n\', head6);
fprintf(f1,\'%s\\n\', head7);
fprintf(f1,\'FoamFile\\n{\\n\');
fprintf(f1,\'%s\\n\',dict1);
fprintf(f1,\'%s\\n\',dict2);
fprintf(f1,\'%s\\n\',dict3);
fprintf(f1,\'%s\\n\',dict4);
fprintf(f1,\'%s\\n\',dict5);
fprintf(f1,\'%s\\n\',conv);
%printing the vertices
fprintf(f1,\'vertices\\n\');
fprintf(f1,\'(\\n\');
fprintf(f1,\'%s\\t(%d %d %d)\\n\',b4,v0(1),v0(2),v0(3)); %point 0
fprintf(f1,\'%s\\t(%d %f %f)\\n\',b4,v1(1),v1(2),v1(3));%point 1
fprintf(f1,\'%s\\t(%d %f %f)\\n\',b4,v2(1),v2(2),v2(3));%point 2
fprintf(f1,\'%s\\t(%f %d %d)\\n\',b4,v3(1),v3(2),v3(3));%point 3
fprintf(f1,\'%s\\t(%f %f %f)\\n\',b4,v4(1),v4(2),v4(3));%point 4
fprintf(f1,\'%s\\t(%f %f %f)\\n\',b4,v5(1),v5(2),v5(3));%point 5
fprintf(f1,\');\\n\\n\');
% printing the blocks
fprintf(f1,\'blocks\\n\');
fprintf(f1,\'(\\n\');
%fprintf(f1, \'hex (0 1 2 0 3 4 5 3) (%d %d %d)\', nx,ny,nz);
fprintf(f1, \'hex (0 3 5 2 0 3 4 1) (%d %d %d)\', nx,ny,nz);
fprintf(f1, \' simpleGrading (1 0.1 1)\\n\\n\');
fprintf(f1, \');\\n\\n\');
% printing the edges
fprintf(f1,\'edges\\n(\\n\');
fprintf(f1,\'%sarc 1 2 (0 %0.3f 0)\\n\',b4,r);
fprintf(f1,\'%sarc 4 5 (%0.4f %0.3f 0)\\n\',b4,L,r);
% fprintf(f1,\'%sarc 1 2 (%d %d %d)\\n\',b4,r);
% fprintf(f1,\'%sarc 4 5 (%d %d %d)\\n\',b4,L,r);
fprintf(f1, \');\\n\\n\');
% printing the boundaries
fprintf(f1, \'\\n\\n boundary\\n(\\n\');
%fprintf(f1, \'%swall\\n%s{\',b4,b4);
%fprintf(f1, \'%stype wall;\\n%sfaces\\n%s(\\n%s(2 5 4 1)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
fprintf(f1, \'%stop\\n%s{\',b4,b4);
fprintf(f1, \'%stype patch;\\n%sfaces\\n%s(\\n%s(1 2 5 4)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
% fprintf(f1, \'\\t pipewall\\n\');
% fprintf(f1, \'\\t{\\n\');
% fprintf(f1, \'\\t\\t type wall;\\n\');
% fprintf(f1, \'\\t\\t faces\\n\');
% fprintf(f1, \'\\t\\t(\\n\');
% fprintf(f1,\'\\t\\t\\t (%d %d %d %d)\\n\',2,5,4,1);
% fprintf(f1,\'\\t\\t\\t (%d %d %d %d)\\n\',1,2,5,4);
% fprintf(f1, \'\\t\\t);\\n\');
% fprintf(f1, \'\\t}\\n\');
fprintf(f1, \'%sinlet\\n%s{\',b4,b4);
%fprintf(f1, \'%stype patch;\\n%sfaces\\n%s(\\n%s(0 0 2 1)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
fprintf(f1, \'%stype patch;\\n%sfaces\\n%s(\\n%s(0 1 2 0)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
% fprintf(f1, \'\\t inlet\\n\');
% fprintf(f1, \'\\t{\\n\');
% fprintf(f1, \'\\t\\t type patch;\\n\');
% fprintf(f1, \'\\t\\t faces\\n\');
% fprintf(f1, \'\\t\\t(\\n\');
% fprintf(f1,\'\\t\\t\\t (%d %d %d %d)\\n\',0,0,2,1);
% fprintf(f1,\'\\t\\t\\t (%d %d %d %d)\\n\',0,1,2,0);
% fprintf(f1, \'\\t\\t);\\n\');
% fprintf(f1, \'\\t}\\n\');
fprintf(f1, \'%soutlet\\n%s{\',b4,b4);
fprintf(f1, \'%stype patch;\\n%sfaces\\n%s(\\n%s(3 4 5 3)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
% fprintf(f1, \'\\t outlet\\n\');
% fprintf(f1, \'\\t{\\n\');
% fprintf(f1, \'\\t\\t type patch;\\n\');
% fprintf(f1, \'\\t\\t faces\\n\');
% fprintf(f1, \'\\t\\t(\\n\');
% fprintf(f1,\'\\t\\t\\t (%d %d %d %d)\\n\',3,4,5,3);
% fprintf(f1, \'\\t\\t);\\n\');
% fprintf(f1, \'\\t}\\n\');
fprintf(f1, \'%swedgeFront\\n%s{\',b4,b4);
%fprintf(f1, \'%stype wedge;\\n%sfaces\\n%s(\\n%s(0 3 5 2)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
fprintf(f1, \'%stype wedge;\\n%sfaces\\n%s(\\n%s(0 2 5 3)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
% fprintf(f1, \'\\t front\\n\');
% fprintf(f1, \'\\t{\\n\');
% fprintf(f1, \'\\t\\t type wedge;\\n\');
% fprintf(f1, \'\\t\\t faces\\n\');
% fprintf(f1, \'\\t\\t(\\n\');
% fprintf(f1,\'\\t\\t\\t (%d %d %d %d)\\n\',0,3,5,2);
% fprintf(f1, \'\\t\\t);\\n\');
% fprintf(f1, \'\\t}\\n\');
fprintf(f1, \'%swedgeBack\\n%s{\',b4,b4);
fprintf(f1, \'%stype wedge;\\n%sfaces\\n%s(\\n%s(0 1 4 3)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
% fprintf(f1, \'\\t back\\n\');
% fprintf(f1, \'\\t{\\n\');
% fprintf(f1, \'\\t\\t type wedge;\\n\');
% fprintf(f1, \'\\t\\t faces\\n\');
% fprintf(f1, \'\\t\\t(\\n\');
% fprintf(f1,\'\\t\\t\\t (%d %d %d %d)\\n\',0,1,4,3);
% fprintf(f1, \'\\t\\t);\\n\');
% fprintf(f1, \'\\t}\\n\');
fprintf(f1, \'%saxis\\n%s{\',b4,b4);
fprintf(f1, \'%stype empty;\\n%sfaces\\n%s(\\n%s(0 3 3 0)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
% fprintf(f1, \'\\t axis\\n\');
% fprintf(f1, \'\\t{\\n\');
% fprintf(f1, \'\\t\\t type empty;\\n\');
% fprintf(f1, \'\\t\\t faces\\n\');
% fprintf(f1, \'\\t\\t(\\n\');
% fprintf(f1,\'\\t\\t\\t (%d %d %d %d)\\n\',0,3,3,0);
% fprintf(f1, \'\\t\\t);\\n\');
% fprintf(f1, \'\\t}\\n\');
fprintf(f1, \');\\n\\n\');
% fprintf(f1, \');\\nmergePatchpairs();\\n\\n\');
% fprintf(f1, \'%s\', footer);
fprintf(f1, \'mergePatchPairs\\n(\\n\');
fprintf(f1, \');\\n\');
fclose(f1)
The BlockMeshDict file for the wedge angle 3 is shown below.
/*--------------------------------*- 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)
(0 0.009998 0.000175)
(0 0.009998 -0.000175)
(2.820000 0 0)
(2.820000 0.009998 0.000175)
(2.820000 0.009998 -0.000175)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (400 40 1) simpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 0.010 0)
arc 4 5 (2.8200 0.010 0)
);
boundary
(
top
{ type patch;
faces
(
(1 2 5 4)
);
}
inlet
{ type patch;
faces
(
(0 1 2 0)
);
}
outlet
{ type patch;
faces
(
(3 4 5 3)
);
}
wedgeFront
{ type wedge;
faces
(
(0 2 5 3)
);
}
wedgeBack
{ type wedge;
faces
(
(0 1 4 3)
);
}
axis
{ type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //ControlDict 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 simpleFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 12;
deltaT 0.000008;
writeControl timeStep;
writeInterval 50;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Pressure Boundary Condition:
/*--------------------------------*- 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
{
axis
{
type empty;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
top
{
type zeroGradient;
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //Velocity Boundary condition:
/*--------------------------------*- 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.09375 0 0);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (0.09375 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //[IV] Mesh Profile
[V] Simulation Results
1. Wedge Angle: 3
Contour plots
A. Velocity Profile near the entrance[Hydrodynamic Boundary Layer Formation]
B. Velocity Profile near the exit[Fully Developed Flow]
C. Velocity Profile curves
i) Velocity Plot at x=0.02 m from the entrance
Here in this region the wall has the highest viscosity and there is a formation of boundary layer which merges slowly with the increase in the axial distance which can be seen later from the contour plots i.e when the fluid reaches the hydrodynamic entry length.
The velocity profile near the entrance shows 3 regions:
i)Viscous sublayer: This is the region of formation of velocity Boundary layer
ii)Transition Buffer layer: This layer gradually decreases with the viscous sublayer and merges with the turbulent core with increasing distance from the entrance
iii) Turbulent layer: The velocity profile is almost flat in this region which appears with the increase in length
ii) Velocity profile at x = 0.1 m from the entrance
This plot shows the Velocity Boundary Layer increases with size at the cost of shrinking turbulent core and the Transitional Buffer Layer
iii) Velocity Profile at x= 0.25 m from the entrance
iv) Velocity Profile at x= 0.5 m from the entrance
The velocity has increased as it moves further and the hydrodynamic Boundary layer is still present at the top and slowly trying to form the curve as it approaches towards the hydrodynamic entrance length.
v) Velocity Profile at x =1.0 m from the entrance
Here some buffers are still observed at the top wall region
vi)Velocity profile for fully developed laminar flow(at x = L_entrance = 2.52 m)
Shear Stress Distribution for fully developed flow:
It is zero at the centre line and maximum at the walls and there is slight non-linearity near the edges.
2. Velocity Profiles for other wedge angles:
WEDGE ANGLE 2
A. Mesh Profile:
B. Velocity profile near entrance(Hydrodynamic Boundary Layer Formation)
C.Velocity profile near exit(fully developed flow)
WEDGE ANGLE 4
A. Mesh Profile
wedge
inlet(Hydrodynamic Boundary layer formation region-viscous sublayer)
exit(fully developed region)
B. Velocity profile near the entrance(Hydrodynamic Boundary Layer Formation)
C. Velocity profile at the exit(Fully developed Flow)
[VI]CONCLUSIONS
Comparison of maximum velocities `V_max`
`V_(avg) = V_(INLET)`or the velocity provided at the entry length.
We have compared the maximum velocities for three different wedge angles compared the datas with the data obtained from Hagen Poisueille Equation.
We have also compared the mean and max Courant No. and the Simulation(Execution time) taken for the 3 different wedge angles.
The computed velocities has errors and does not reach the velocities that is given by Hagen Poisueille equation but definitely tends to do so.
The execution time increases in an ascending order for wedge angles 2-3-4.