Week 11 - Simulation of Flow through a pipe in OpenFoam
TITLE
Simulation of flow through pipe
OBJECTIVES
- Validate Hydro-dynamic length with the numerical result
- Validate the fully developed flow velocity profile with its analytical profile
- Validate maximum velocity and pressured drop for fully developed flow
- Post-process Shear stress and validate wall shear stress for fully developed flow
THEORY
Consider a pipe of diametre 'r' and length 'l'. A fluid enters from the inlet at a velocity 'v' , and pressure 'p'. The friction from the walls of the pipe and other effects originating from the interior wall of the pipe propagate into the flow as an expanding boundary layer. This boundary layer tries to expand throughout the entire length of the pipe and attains an equilibrium point (knows as Hydrodynamic Entrance Length). After which the flow is fully developed. Below image shows us the graphical explanation.
In nonideal fluid dynamics, the Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, Poiseuille law or Poiseuille equation, 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 constant cross section.
Equations are:
In this case as we are assuming that the flow is Laminar, the equation for the hydrodynamic entrance length is given as
L = 0.05 * reynolds number * Diametre(of the pipe).
`L = 0.05*Re*D`
The velocity of the fluid will be
`V = (Re*mu)/(rho*D)`
The maximum velocity ( which will be attained at fully developed flow)
`V_max = 2*V`
Pressure
`P = (32*mu*V*L)/(D^2)`
Kinematic Pressure
`P_k = P/rho` ( density of water assume as 1000)
Maximum Shear stress
`tau_max = (2*mu*V_max)/R`
With the Avove equations we calculate.
L = 0.42m
V = 0.46725 m/s
`V_max` = 0.9345 m/s
P = 349.3161 pa
`P_k` = 0.3493161 pa
`tau_max` = `mu * ((du)/dy)`
CODES
The code below is the Matlab code for file parsing and accessing the data from the file. It also contains code for calculating the equations analytically.
% file parsing for OpenFOAM
clear all
close all
clc
% Diameter of pipe
D = 0.004; % All dimensions in m (Note that arc points to be changed)
% Reynolds Number
Re = 2100;
% Entry length
L = 0.05 * Re * D; % (result = 0.42)
% Length of pipe (For entry length = 0.42)
L_p = L + 0.58;
% Radius of pipe
r = D/2;
% Wedge angle
theta = 2; % In degrees
% points in the Geometry
Y = r*cosd(theta/2);
Z = r*sind(theta/2);
%Dynamic Viscosity
mu = 8.9e-04;
%Density
rho = 1000; %(density of water is 997, assume 1000.)
%HP equation
v_avg = (Re * mu) / (rho/D);
v_max = 2 * (v_avg);
del_p = (8 * mu * L_p * v_avg) / (r^2);
kinematic_p = del_p / rho;
% Calculation of shear stress
R = linspace(0,D/2,1000);
for i = 1:length(R)
tau(i) = 2 *mu * v_max* (abs(R(i)))/(r^2);
end
plot(R,tau)
xlabel('Radius')
ylabel('Shear Stress')
grid on
legend('SHEAR STRESS')
% File parsing
h1= '/*--------------------------------*- C++ -*----------------------------------*\'
h2=' ========= |'
h3=' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox'
h4=' \\ / O peration | Website: https://openfoam.org'
h5=' \\ / A nd | Version: 8'
h6=' \\/ M anipulation |'
h7='\*---------------------------------------------------------------------------*/'
h8='FoamFile'
h9='{'
h10=' version 2.0;'
h11=' format ascii;'
h12=' class dictionary;'
h13=' object blockMeshDict;'
h14='}'
h15='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
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,'%s \n','convertToMeters 1;');
fprintf(f1,'\n vertices');
fprintf(f1,'\n (');
fprintf(f1,'\n\t (%d %d %d)',0,0,0);
fprintf(f1,'\n\t (%d %d %d)',0,Y,Z);
fprintf(f1,'\n\t (%d %d %d)',0,Y,-Z);
fprintf(f1,'\n\t (%d %d %d)',1,0,0);
fprintf(f1,'\n\t (%d %d %d)',1,Y,Z);
fprintf(f1,'\n\t (%d %d %d)',1,Y,-Z);
fprintf(f1,'\n );');
fprintf(f1,'\n blocks');
fprintf(f1,'\n (');
fprintf(f1,'\n hex (0 3 5 2 0 3 4 1) (500 50 1) simpleGrading (1 1 1)');
fprintf(f1,'\n );');
fprintf(f1,'\n');
fprintf(f1,'\n edges');
fprintf(f1,'\n (');
fprintf(f1,'\n arc 1 2 (0 1e-03 0)');
fprintf(f1,'\n arc 4 5 (0 1e-03 0)');
fprintf(f1,'\n );');
fprintf(f1,'\n ');
fprintf(f1,'\n boundary');
fprintf(f1,'\n (');
fprintf(f1,'\n inlet');
fprintf(f1,'\n {');
fprintf(f1,'\n type patch;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (0 1 2 0)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n outlet');
fprintf(f1,'\n {');
fprintf(f1,'\n type patch;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (3 5 4 3)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n front');
fprintf(f1,'\n {');
fprintf(f1,'\n type wedge;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (0 3 4 1)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n back');
fprintf(f1,'\n {');
fprintf(f1,'\n type wedge;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (0 2 5 3)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n top');
fprintf(f1,'\n {');
fprintf(f1,'\n type wall;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (5 2 1 4)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n axis');
fprintf(f1,'\n {');
fprintf(f1,'\n type empty;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (0 3 3 0)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n );');
fprintf(f1,'\n ');
fprintf(f1,'\n mergePatchPairs');
fprintf(f1,'\n (');
fprintf(f1,'\n );');
fprintf(f1,'\n ');
fprintf(f1,'\n // ************************************************************************* //');
fprintf(f1,'\n ');
BLOCKMESHDICT
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 1.999695e-03 3.490481e-05)
(0 1.999695e-03 -3.490481e-05)
(1 0 0)
(1 1.999695e-03 3.490481e-05)
(1 1.999695e-03 -3.490481e-05)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (500 30 1) simpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 1e-03 0)
arc 4 5 (0 1e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
top
{
type wall;
faces
(
(5 2 1 4)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //PRESSURE P
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ 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.3493161;
}
front
{
type wedge;
}
back
{
type wedge;
}
top
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //VELOCITY U
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ 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.46725 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
top
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //CONTROLDICT
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 3;
deltaT 0.00001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //RESULTS
Below image shows the mesh profile and grading used.
Side profile of the wedge created is shown below. So instead of taking a whole circular profile we make use of a wedge geometry to find out all the useful properties. Because the the way these properties vary in the wedge is going to be same throughout its circular profile. This can be called as Finite Volume method.
The following show the results of the properties varying across and through out the wedge with their orientations.
The first one how velocity of the fluid varies. As you can see the velocity is highest at the centre. ( ar radial distance 0)
x-axis : radial distance(m),y-axis : velocity(m/s)
The pressure variation is given below, pressure usually reaches zero as fluid moves towards the outlet. This is according to HP law. As you can see the graph shows the same nature explained in the HP law.
x-axis : radial distance(m),y-axis : pressure (Pa)
Again according to HP law, shear stress is maximum at the walls(because of the no-slip condition) and highest at the centre. From the equation we can understand that shear stress is dependant or directly proportional to the velocity gradient.
x-axis : radial distance(m),y-axis :Shear-stress
Variation of Velocity
Variation of Pressure
Shear stress vs radius in matlab
