Week 11 - Simulation of Flow through a pipe in OpenFoam
AIM: To write a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes. For a baseline mesh, show that the velocity profile matches with the Hagen poiseuille's equation.
Hagen Poiseuille's Law:
It 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.
The assumptions of the equation are that the fluid is incompressible and Newtonian; the flow is laminar through a pipe of constant circular cross-section that is substantially longer than its diameter; and there is no acceleration of fluid in the pipe. For velocities and pipe diameters above a threshold, actual fluid flow is not laminar but turbulent, leading to larger pressure drops than calculated by the Hagen–Poiseuille equation.
Considering a laminar fluid flow inside a pipe. The surface contact between the water and the wall of the pipe will cause the water molecules to slow down and come to a complete stop near the walls. This layer of water molecules will also cause the adjacent layer of water flow to slow down due to the result of friction. To compensate for this velocity, water drops in the core of the pipe needs to maintain a constant velocity to keep the volumetric flow rate constant. As a result, the velocity boundary layer develops along the boundary walls of the pipe and starts to develop in the flow direction until it hits the center of the pipe.
To run the simulation the geometry is divided into patches and applied boundary conditions.
Calculation:
Reynold Number Re = 1400
Diameter D = 0.025
Kinematic Viscosity = 1.0218*e^(-6)*m^(2)/s
Dynamic Viscosity(mu) = 0.00102Pa.s
Density(rho) = 998kg/m^3
According to Hagen Poiseullie's Derivation:
Velocity Vavg = Re*mu/(rho*D)=0.058m/s
Maximum Velocity Umax = 2*V = 0.116m/s
Length at entry Le = 0.05*Re*D = 1.75m
Length L = 1.75m
Pressure Drop = deltaP = 8*mu*V*L/r^2 = 5.38Pa
Kinematic Pressure Drop = deltaPk = deltaP/rho = 0.005Pa
Darcy Friction Factor = f = 64/Re = 0.045Matlab Code to create blockMeshDict:
Matlab is used to generate the blockMeshDict file for the simulation.
clear all
close all
clc
% Creating blockMeshDict file for wedge creation;
D = 0.025; % Diamater of Pipe
R = 1400; % Reynolds Number
L = (0.06*R*D); % Length of Pipe
theta = 2; % Wedge Angle
% File parsing for BlockMeshDict file;
f1 = fopen('blockMeshDict.txt','w');
fprintf(f1,'%s','/*--------------------------------*- C++ -*----------------------------------*');
fprintf(f1,'\n');
fprintf(f1,'%s',' ========= |');
fprintf(f1,'\n');
fprintf(f1,'%s',' \ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'\n');
fprintf(f1,'%s',' \ / O peration | Website: https://openfoam.org');
fprintf(f1,'\n');
fprintf(f1,'%s',' \ / A nd | Version: 6');
fprintf(f1,'\n');
fprintf(f1,'%s',' \/ M anipulation |');
fprintf(f1,'\n');
fprintf(f1,'%s','*---------------------------------------------------------------------------*/');
fprintf(f1,'\n');
fprintf(f1,'%s','FoamFile');
fprintf(f1,'\n');
fprintf(f1,'%s','{');
fprintf(f1,'\n');
fprintf(f1,'%s',' version 2.0;');
fprintf(f1,'\n');
fprintf(f1,'%s',' format ascii;');
fprintf(f1,'\n');
fprintf(f1,'%s',' class dictionary;');
fprintf(f1,'\n');
fprintf(f1,'%s',' object blockMeshDict;');
fprintf(f1,'\n');
fprintf(f1,'%s','}');
fprintf(f1,'\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s','convertToMeters 1;');
fprintf(f1,'\n');
fprintf(f1,'\n');
% vertice creation;
fprintf(f1,'%s','vertices');
fprintf(f1,'\n');
fprintf(f1,'(\n');
fprintf(f1,' (%d %d %d)\n',0,0,0);
fprintf(f1,' (%d %d %d)\n',0,(D/2)*cosd(theta/2),(D/2)*sind(theta/2));
fprintf(f1,' (%d %d %d)\n',0,(D/2)*cosd(theta/2),-(D/2)*sind(theta/2));
fprintf(f1,' (%d %d %d)\n',L,0,0);
fprintf(f1,' (%d %d %d)\n',L,(D/2)*cosd(theta/2),(D/2)*sind(theta/2));
fprintf(f1,' (%d %d %d)\n',L,(D/2)*cosd(theta/2),-(D/2)*sind(theta/2));
fprintf(f1,');\n');
% Block creation;
fprintf(f1,'blocks\n');
fprintf(f1,'(\n');
fprintf(f1,' hex (0 1 2 0 3 4 5 3) (50 1 100) simpleGrading (%d %d %d)\n',0.1,1,10);
fprintf(f1,');\n');
% Edge creation;
fprintf(f1,'edges\n');
fprintf(f1,'(\n');
fprintf(f1,' arc 1 2 (0 %d 0)\n',D/2);
fprintf(f1,' arc 4 5 (%d %d 0)\n',L,D/2);
fprintf(f1,');\n');
% Boundary creation;
fprintf(f1,'boundary\n');
fprintf(f1,'(\n');
%inlet face;
fprintf(f1,' inletwall\n');
fprintf(f1,' {\n');
fprintf(f1,' type patch;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (0 1 2 0)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
%Outlet face;
fprintf(f1,' outletwall\n');
fprintf(f1,' {\n');
fprintf(f1,' type patch;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (3 4 5 3)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
%Top face;
fprintf(f1,' Noslipwall\n');
fprintf(f1,' {\n');
fprintf(f1,' type wall;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (1 4 5 2)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
%wedge front;
fprintf(f1,' front\n');
fprintf(f1,' {\n');
fprintf(f1,' type wedge;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (0 3 5 2)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
%wedge Back;
fprintf(f1,' back\n');
fprintf(f1,' {\n');
fprintf(f1,' type wedge;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (0 1 4 3)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
%wedge axis;
fprintf(f1,' axis\n');
fprintf(f1,' {\n');
fprintf(f1,' type empty;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (0 3 3 0)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
fprintf(f1,');\n');
%Patchpairs;
fprintf(f1,'mergePatchPairs\n');
fprintf(f1,'(\n');
fprintf(f1,');\n')
fprintf(f1,'// ************************************************************************* //');
fclose(f1);Result of matlab file:
the above matlab code generates the following blockMeshDict file
blockMeshDict 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.0124981 -0.000218155)
(0 0.0124981 0.000218155)
(2.1 0 0)
(2.1 0.0124981 -0.000218155)
(2.1 0.0124981 0.000218155)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 100) simpleGrading (0.1 1 10)
);
edges
(
arc 1 2 (0 0.0125 0)
arc 4 5 (2.1 0.0125 0)
);
boundary
(
inletWall
{
type patch;
faces
(
(0 1 2 0)
);
}
outletWall
{
type patch;
faces
(
(3 4 5 3)
);
}
Noslipwall
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
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: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 0.1;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Velocity File:
/*--------------------------------*- 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.05 0 0);
boundaryField
{
inletWall
{
type fixedValue;
value uniform (1 0 0);
}
outletWall
{
type zeroGradient;
}
Noslipwall
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Pressure File:
/*--------------------------------*- 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
{
inletWall
{
type zeroGradient;
}
outletWall
{
type fixedValue;
value uniform 1;
}
Noslipwall
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Results & Conclusion:
The results generated by using OpenFoam are:
at angle 2:
Pressure Variation:
Velocity Variation:
Velocity at initial
Velocity at final
It can be observed that the inlet velocity of 0.05m/s and we got final velocity 0.1m/s. The velocity at walls is zero and increasing at the centre of pipe.
From the calculations we have hydrodynamic length of 1.75m.