Simulation of Flow through a pipe
Hydrodynamic Entrance length :
The hydrodynamic entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces propagating from the interior wall of a pipe.This region is characterized by a non-uniform flow.The fluid enters a pipe at a uniform velocity, then fluid particles in the layer in contact with the surface of the pipe come to a complete stop due to the no slip condition. Due to viscous forces within the fluid, the layer in contact with the pipe surface resists the motion of adjacent layers and slows adjacent layers of fluid down gradually, forming a velocity profile.For the conservation of mass to hold true, the velocity of layers of the fluid in the center of the pipe increases to compensate for the reduced velocities of the layers of fluid near the pipe surface. This develops a velocity gradient across the cross section of the pipe.
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.
In case of laminar flow, the velocity profile in the fully developed region is parabolic.
Reylonds number = 2100
Fluid -water
The below properties of water are at 25degress
Density of water ρ=997.04kg/m3
Dynamic viscosity ,μ=8.90∗10^−4 pa.sec
kinematic viscosity, ν=μ/ρ=8.90∗10^−4/997.04
=8.97*10^-7 m^2/sec
Diameter of the pipe (d)=0.02 meters;
Hydrodynamic length of the pipe (l)=0.06 * R*d
=0.06*2100*0.02 = 2.52 m
Length of the pipe (L)=0.3+hydrodynamic length
=0.3+2.52 =2.82 m
Average velocity, Vavg= (R_e μ)/ρD =0.09373 m/sec
Maximum velocity, Vmax=2.Vavg =0.1875 m/sec
Pressure drop (p) = (32μLV_avg)/D^2 =18.81 pa
Kinematic pressure drop =∆p/ρ=0.0018862793
The below shown code is written in Matlab
clear all
close all
clc
% input parameters
Re=2100;%reynolds number
D=0.02; %Diameter of the pipe
r=D/2; %radius of the pipe
rho=997; % density of the water at 25 degc
L_entrance=0.06*Re*D; % Entry length of the pipe
L=0.3+L_entrance; %Length of the pipe
mu=8.90e-4; %dynamic viscosity
nu=mu/rho; %Kinematic viscosity
V_avg=(Re*mu)/(rho*D);% Average velocity
v_max=2*V_avg; %Maximum elocity
del_P=(32*mu*L*V_avg)/D^2; % pressure drop
kin_p=del_P/rho;%kinematicPressure drop
fprintf(\'Enter the theta values of 2,3,4 degrees for wedge angle \\n\');
theta=input(\'Enter the wedge angle in degrees\'); %wedge angle
%calculation of the vertices
v1=[0 0 0];
v2=[0 r*cosd(theta/2) -r*sind(theta/2)];
v3=[0 r*cosd(theta/2) r*sind(theta/2)];
v4=[L 0 0];
v5=[L r*cosd(theta/2) -r*sind(theta/2)];
v6=[L r*cosd(theta/2) r*sind(theta/2)];
% header format for the blockmesh dict file
head1= \'/*----------------------------------*- C++ -*-----------------------------------*\\\';
head2= \' ========= |\';
head3= \' \\\\ / F ield |OpenFOAM: The Open Source CFD Toolbox\';
head4= \' \\\\ / O peration |Website: https://openfoam.org\';
head5= \' \\\\ / A nd |Version: 7\';
head6= \' \\\\/ M anipulation |\';
head7= \'\\*----------------------------------------------------------------------------*/\';
dict1=\'version 2.0;\';
dict2=\'format ascii;\';
dict3=\'class dictionary;\';
dict4=\'object blockMeshDict;\';
dict5=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';
conv=\'ConvertToMeters 1;\';
vertices1=\'vertices\';
block1=\'hex (0 3 4 1 0 3 5 2) (350 50 1) simplegrading (1 0.1 1)\';
footer=\'//***********************************//\';
% creating the blockmesh file for OpenFOAM
f1=fopen(\'blockMeshDict.txt\',\'w\');
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,\'%s \\n{\\n\\t\',\'FoamFile\');
fprintf(f1,\'%s \\n\\t\',dict1);
fprintf(f1,\'%s \\n\\t\',dict2);
fprintf(f1,\'%s \\n\\t\',dict3);
fprintf(f1,\'%s \\n}\\n\',dict4);
fprintf(f1,\'%s \\n\\n\',dict5);
fprintf(f1,\'%s \\n\\n\',conv);
fprintf(f1,\'%s \\n(\\n\\t\',vertices1);
fprintf(f1,\'(%f %f %f)\\n\\t\',v1(1),v1(2),v1(3));
fprintf(f1,\'(%f %f %f)\\n\\t\',v2(1),v2(2),v2(3));
fprintf(f1,\'(%f %f %f)\\n\\t\',v3(1),v3(2),v3(3));
fprintf(f1,\'(%f %f %f)\\n\\t\',v4(1),v4(2),v4(3));
fprintf(f1,\'(%f %f %f)\\n\\t\',v5(1),v5(2),v5(3));
fprintf(f1,\'(%f %f %f)\\n);\\n\',v6(1),v6(2),v6(3));
%printing the blocks
fprintf(f1,\'\\nblocks\\n\');
fprintf(f1,\'(\\n\\t\');
fprintf(f1,\'%s\\n);\\n\',block1);
% printing the edges
fprintf(f1,\'\\nedges\\n(\\n\');
fprintf(f1,\'\\tarc 1 2 (0 %f 0)\\n\',r);
fprintf(f1,\'\\tarc 4 5 (0 %f %f )\\n);\\n\',L,r);
% Printing the boundaries
fprintf(f1,\'\\nboundary \\n(\\n\');
fprintf(f1,\'\\taxis\\n\\t{\\n\\t\\ttype empty;\\n\');
fprintf(f1,\'\\t\\tfaces\\n\\t\\t(\\n\\t\\t\\t(0 3 3 0)\\n\\t\\t);\\n\\t}\\n\');
fprintf(f1,\'\\twall\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype wall;\\n\\t\\tfaces\\n\\t\\t(\\n\');
fprintf(f1,\'\\t\\t\\t(2 5 4 1)\\n\\t\\t);\\n\\t}\\n\');
fprintf(f1,\'\\tinlet\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype patch;\\n\\t\\tfaces\\n\');
fprintf(f1,\'\\t\\t(\\n\\t\\t\\t(0 1 2 0)\\n\\t\\t);\\n\\t}\\n\');
fprintf(f1,\'\\toutlet\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype patch;\\n\\t\\tfaces\\n\');
fprintf(f1,\'\\t\\t(\\n\\t\\t\\t(3 4 5 3)\\n\\t\\t);\\n\\t}\\n\');
fprintf(f1,\'\\tfront\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype wedge;\\n\\t\\tfaces\\n\');
fprintf(f1,\'\\t\\t(\\n\\t\\t\\t(0 3 5 2)\\n\\t\\t);\\n\\t}\\n\');
fprintf(f1,\'\\tback\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype wedge;\\n\\t\\tfaces\\n\');
fprintf(f1,\'\\t\\t(\\n\\t\\t\\t(0 1 4 3)\\n\\t\\t);\\n\\t}\\n);\\t\');
fprintf(f1,\'\\nmergePatchpairs \\n(\\n);\\n\');
fprintf(f1,\'\\n%s\',footer);
The blockmeshdict is generated by writing the code in matlab as shown in below for 2 degrees:
/*----------------------------------*- 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)
(0.000000 0.009998 -0.000175)
(0.000000 0.009998 0.000175)
(2.820000 0.000000 0.000000)
(2.820000 0.009998 -0.000175)
(2.820000 0.009998 0.000175)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (350 50 1) simplegrading (1 0.1 1)
);
edges
(
arc 1 2 (0 0.010000 0)
arc 4 5 (0 2.820000 0.010000 )
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchpairs
(
);
//***********************************//After succesfully generation of the blockMeshdictfile.we have to change the pressure and velocity files in the Zero folder(0).Those both files are shown in below.
Initial and Boundary conditions for Pressure 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 volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
axis
{
type empty;
}
wall
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //Initial and boundary Conditions for Velocity 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 volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.09373 0 0);
boundaryField
{
axis
{
type empty;
}
wall
{
type noSlip;
}
inlet
{
type fixedValue;
value uniform (0.09373 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //In this velocity File, we provide the inlet velocity by calculating the average velocity by using hagen poiseuilli\'s equation.The inlet velocity is in X direction
ControlDict 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 \"system\";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10;
deltaT 0.001;
writeControl timeStep;
writeInterval 100;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //In the controlDictfile we have to specify the startTime,endTime and DeltaT.In this we provide the 0 for startTime,10 seconds for endTime and 0.001 for DeltaT.Select the DeltaT by a suitable number not exceeding the Courant number greater than 1.So be careful while selecting the DeltaT.
Transport Properties
/*--------------------------------*- 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] 8e-7;
// ************************************************************************* //After all the values entered in different folders which contains the file.
Now run the icoFoam command,this start the simulation of the project and takes some time to complete the project.After successfully completion of the icoFoam command then type the paraFoam command in Terminal space to view the Mesh genereted file and also the velocity and pressure generated files.By using this files we can generate the graphs and results.
wedge angle for 2 degres:
Mesh File:
Velocity Profiles for the pipe:
At point 0
At point 0.5
At point 1
At point 2.55
At point 2.7
The above graphs are the velocity profiles of the pipe at various locations from the inlet pipe.
At the subsequent intervals of the distances the velocity profile is gradually developing.
At the distance of 2.52m the velocity profile is fully developed from this area onwards the velocity profile looks similiar .It is the length which was calculated by using entrance length formula.
SHEAR STRESS:
It can be calculated by saving the data of all velocity profiles.It can be saved in .CSC file format.After generating the velocity profile at the peak point(i.e fully developed region) we can save the data by the file command.
Thus the generated data can be export to the MS Excel where we can do some operations in order to get the Shear stress.It can be shown below
Firstly we select the Magnitude of X velocity column and the Y points ,then the do the operations as per the formula for velocity gradient(du/dy) and multiple with the Dynamic viscosity (mu *(du/dy)).Thus the shearstress is obtainted and plot the results.
From the data we can select the column A for velocity and column y for the distance between the points for calculating the velocity gadient (du/dy)
From the graph we can see that the shear stress is minimum at the centre point and maximum at the walls.