Simulation of Laminar Flow through a Pipe using Axi Symmetry and Hagen Poiseuille Equation
Objective: To simulate laminar flow through a pipe using wedge angle and axi symmetry
In a pipe flow the profile develops as shown in the figure below(image taken from wikipedia)
For a laminar flow the parameters can be determined as follows:
The diameter of the pipe is assumed to be 0.025m
The Reynolds Number is 2100 which lies in the laminar region.
`Re=(rhoVD)/mu`
From the above equation the inlet velocity is calculated as
`V=(Remu)/(rhoD)=0.0750m/s^2`
From the Hagen Poisuille Equation:
`Q=(DeltaP piD^4)/(128muL)`
`Q=VA=VpiR^2`
Simplifying:
`V=(DeltaPD^2)/(32muL)`
`DeltaP=(32VmuL)/(D^2)=10.78`
The kinematic Pressure is considered in the icoFoam solver
`(DeltaP)/rho=0.0108 m^2/s^2`
the entry Length is calculated using the formula:
`L_e/D=0.06Re`
`L_e=0.06ReD=3.15m`
the length of the pipe is hance taken as 1.89m.
A matlab script is written to generate a blockMesh file
clear all
close all
clc
%inputs
D=input(\'enter Diameter\');
theta=input(\'Enter Wedge angle \');
Re=2100;
rho=997;
mu=0.89e-3;
g=input(\'Mesh Grading\');
s=1e-2;
dr=2e-4;
%pressure drop calculation
V=(Re*mu)/(rho*D);
Vmax=2*V;
A=pi*D^2/4;
Q=A*V;
L=0.06*Re*D;
DP=(32*mu*L*V)/D^2;
dp=DP/rho;
%variables
r=D/2;
x1=r*cosd(theta/2);
x2=r*cosd(theta/2);
y1=-r*sind(theta/2);
y2=r*sind(theta/2);
%generate text 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=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';
f1=fopen(\'blockMeshDict\',\'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,\'FoamFile\\n{\\n\\t version\\t \\t 2.0;\\n\\t format \\t \\t ascii;\\n\\t class \\t \\t dictionary; \\n\\t object \\t \\t blockMeshDict;\\n}\\n \');
fprintf(f1,\'%s\\n\',text8);
fprintf(f1,\'\\nconvertToMeters 1;\\n\');
fprintf(f1,\'vertices\\n(\\n\');
fprintf(f1,\'\\t(0 0 0)\\n\\t(%8.7f %8.7f 0)\\n\\t(%8.7f %8.7f 0)\\n\',x1,y1,x2,y2);
fprintf(f1,\'\\t(0 0 %8.3f)\\n\\t(%8.7f %8.7f %8.3f)\\n\\t(%8.7f %8.7f %8.3f)\\n);\\n\',L,x1,y1,L,x2,y2,L);
fprintf(f1,\'\\nblocks\\n(\\n\\t\\t hex (0 1 2 0 3 4 5 3) (%2.0f 1 %3.0f) simpleGrading ( %2.1f 1 1)\\n);\\n\',r/dr,L/s,g);
fprintf(f1,\'\\nedges\\n(\\n\\tarc 1 2 (0 %8.3f 0)\\n\\tarc 4 5(0 %8.3f %8.3f)\\n);\\n\',r,r,L);
fprintf(f1,\'\\nboundary\\n(\\n\');
fprintf(f1,\'\\taxis\\n\\t{\\n\\t\\ttype empty;\\n\\t\\tfaces\\n\\t\\t(\\n\\t\\t\\t(3 3 0 0)\\n\\t\\t);\\n}\\n\');
fprintf(f1,\'\\twall\\n\\t{\\n\\t\\ttype wall;\\n\\t\\tfaces\\n\\t\\t(\\n\\t\\t\\t(2 5 4 1)\\n\\t\\t);\\n}\\n\');
fprintf(f1,\'\\tinlet\\n\\t{\\n\\t\\ttype patch;\\n\\t\\tfaces\\n\\t\\t(\\n\\t\\t\\t(3 4 5 3)\\n\\t\\t);\\n}\\n\');
fprintf(f1,\'\\toutlet\\n\\t{\\n\\t\\ttype patch;\\n\\t\\tfaces\\n\\t\\t(\\n\\t\\t\\t(0 2 1 0)\\n\\t\\t);\\n}\\n\');
fprintf(f1,\'\\tfront\\n\\t{\\n\\t\\ttype wedge;\\n\\t\\tfaces\\n\\t\\t(\\n\\t\\t\\t(1 4 3 0)\\n\\t\\t);\\n}\\n\');
fprintf(f1,\'\\tback\\n\\t{\\n\\t\\ttype wedge;\\n\\t\\tfaces\\n\\t\\t(\\n\\t\\t\\t(5 2 0 3)\\n\\t\\t);\\n}\\n\');
fprintf(f1,\');\\n\');
fprintf(f1,\'%s\\n\',text8);
The script takes the pipe Diameter , wedge angle and the mesh grading as inputs and generates a blockmesh file in the current case the wedge angle is taken as 4 degrees and a mesh grading of 0.3 is applied. Mesh grading is applied near the pipe wall to capture any boundary layer effects.
The blockmesh file generated is as follows:
/*--------------------------------*-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 0 0)
(0.0124924 -0.0004362 0)
(0.0124924 0.0004362 0)
(0 0 3.150)
(0.0124924 -0.0004362 3.150)
(0.0124924 0.0004362 3.150)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (63 1 315) simpleGrading ( 0.3 1 1)
);
edges
(
arc 1 2 (0 0.013 0)
arc 4 5(0 0.013 3.150)
);
boundary
(
axis
{
type empty;
faces
(
(3 3 0 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(3 4 5 3)
);
}
outlet
{
type patch;
faces
(
(0 2 1 0)
);
}
front
{
type wedge;
faces
(
(1 4 3 0)
);
}
back
{
type wedge;
faces
(
(5 2 0 3)
);
}
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //The boundary conditions applied are as floows:
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
{
axis
{
type empty;
}
wall
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0108;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //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 0 -0.0750);
}
wall
{
type noSlip;
}
axis
{
type empty;
}
outlet
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //The Matlab script used to generate the plot of shear stress and pressure drop is as follows:
clear all
close all
clc
f1=fopen(\'velocity.txt\',\'r\');
fgetl(f1);
A = fscanf(f1,\'%f\');
U=-A(1:2:end);
r=A(2:2:end);
mu=0.89e-3;
for i=2:length(U)
dU(i)=U(i-1)-U(i);
dx(i)=r(i)-r(i-1);
end
tau=mu*(dU./dx);
dU=dU\';
fclose(f1);
f2=fopen(\'pressure.txt\',\'r\');
fgetl(f2);
B=fscanf(f2,\'%f\');
p=B(1:2:end);
x=3.15-B(2:2:end);
% subplot(2,1,1)
% plot(r,U,\'linewidth\',2,\'color\',\'r\')
% hold on
% plot(-r, U,\'linewidth\',2,\'color\',\'r\')
% xlabel(\'Radial Distance(m)\');
% ylabel(\'Velocity(m/s)\');
% Title(\'Velocity Profile\');
% grid on;
% subplot(2,1,2)
figure (1)
plot(r,dU,\'linewidth\',2,\'color\',\'r\')
hold on
plot(-r,dU,\'linewidth\',2,\'color\',\'r\')
xlabel(\'Radial Distance(m)\');
ylabel(\'Wall Shear\');
grid on ;
figure(2)
plot(x,p,\'linewidth\',2,\'color\',\'k\')
xlabel(\'Length of the pipe(m)\');
ylabel(\'Kinematic Pressure(m^2/s^2)\')
title(\'Pressure Drop\');
grid on;
Results:
The mesh obtained for a wedge angle of 4 with a mesh grading of 0.3.
the fully developed velocity profile is in the inage below. the velovity at the top end which is neat the wall is zero due to the no slip condition.At the inlet the development of the flow can be observed.
Given below are velocity profiles at different lengths of the pipe.The first image shows the velocity profile at the inlet. The second omage shows the velocity profile that is still developing and the third image shows the fulley developed velocity profile.
From the fully developed velocity profile it can be seen that the mavimum veocity attained is at radial distance zero which is the centre of the pipe. The maximum velocity attained is 0.13 m/s. the maximum velocity obtained analytically is 0.15m/s.
The pressure drop obtained is shown is the plot below.
From the plot it can be observed that the pressre drop is not linear near the inlet. This is because the flow is still developing.Once the flow has developed then the pressure drop is linear.
The wall shear stress can be seen in the graph below.
It can be observed that the shear stress is minimum near the centre and increases towards the pipe wall.