Week 10 - Simulation of Flow through a pipe in OpenFoam
Hagen- Poiseuille's equations :
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. It can be successfully applied to air flow in lung alveoli, or the flow through a drinking straw or through a hypodermic needle.
The Hagen–Poiseuille equation describes the relationship between pressure, fluidic resistance and flow rate, analogous to voltage, resistance, and current, respectively, in Ohm's law for electrical circuits ( V = R I ). Both electrical resistance and fluidic resistance are proportional to the length of the device.
`Deltap=(8muLQ)/(piR^4)`
`Deltap`=pressure difference between the two ends
`mu`= dynamic viscosity
L = length of pipe
Q = Volumetric flow rate
R = pipe radius
PROGRAM CODES AND FILES :
First, we are going to write a MatLab code to automate the generation of the blockMeshDict file, which takes our respective inputs and outputs a fully functional blockMeshDict file.
% SIMULATION OF FLOW THROUGH A PIPE
% Code to automate generation of blockMeshDict file
clear all
close all
clc
% INPUTS
% Assumptions
% Diameter of Pipe (m)
D = 0.02;
% Radius of Pipe (m)
R = D/2;
% Reynold's Number
Re = 2100;
% Wedge Angle (degrees)
theta = 4;
% Total Length of the Pipe (m)
L = 2.8;
% Assuming working fluid is water at 25° C
% Density of water (kg/m³)
rho = 997;
% Dynamic Viscosity of Water (Pa.s)
mu = 8.90 * 10^-4;
% Kinematic Viscosity of Water (m^2/s)
nu = mu/rho;
% Hydrodynamic Entrance Length (m)
L_h = 0.06 * Re * D;
% HAGEN - POISEUILLE'S FLOW EQUATIONS
% Average Velocity (m/s)
V_avg = (mu * Re)/(rho * D);
% Maximum Velocity (m/s)
V_max = 2 * V_avg;
% Pressure Drop (Pa)
delta_p = (32*mu*V_avg*L)/(D^2);
% Kinematic Pressure Drop (Pa)
kinematic_delta_p = delta_p/rho;
% DEFINING THE GEOMETRY
% Vertices
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)];The blockMeshDict file has been written automatically by the MatLab program as shown above, which we use to create the geometry of the wedge and run the simulation for the same.
% Creating BlockMesh File
f1 = fopen('blockMeshDict.txt');
% Header
fprintf(f1,'%s\n %s\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n%s\n\n',h1,h2,h3,h4,h5,h6,h7,h8,h9,h10,h11,h12,h13,h14);
% Conversion Factor
fprintf(f1,'%s\n\n','convertToMeters 1;');
% Vertices
fprintf(f1,'%s\n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t(%d %d %d)\n',V0);
fprintf(f1,'\t(%d %f %f)\n',V1);
fprintf(f1,'\t(%d %f %f)\n',V2);
fprintf(f1,'\t(%.2f %d %d)\n',V3);
fprintf(f1,'\t(%.2f %f %f)\n',V4);
fprintf(f1,'\t(%.2f %f %f)\n',V5);
fprintf(f1,'%s\n\n',');');
% Blocks
fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'\thex (0 3 4 1 0 3 5 2) (50 10 1) simpleGrading (1 1 1)\n');
fprintf(f1,'%s\n\n',');');
% Edges
fprintf(f1,'%s\n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'\tarc 2 1 (%d %.3f %d)\n',0,R,0);
fprintf(f1,'\tarc 4 5 (%.3f %.3f %d)\n',L,R,0);
fprintf(f1,'%s\n\n',');');
% Boundary
fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
% Axis
fprintf(f1,'\taxis\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype empty;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 3 3 0)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Inlet
fprintf(f1,'\tinlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 2 1 0)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Outlet
fprintf(f1,'\toutlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(3 4 5 3)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Front Face
fprintf(f1,'\tfront\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wedge;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 3 5 2)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Back Face
fprintf(f1,'\tback\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wedge;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 1 4 3)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Top Wall
fprintf(f1,'\tpipeWall\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wall;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(1 2 5 4)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','// ************************************************************************* //');After running the program the the blockMesh file have been automatically created which can paste it instead of the previous blockMesh file as shown below:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.009994 -0.000349)
(0 0.009994 0.000349)
(2.80 0 0)
(2.80 0.009994 -0.000349)
(2.80 0.009994 0.000349)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (50 10 1) simpleGrading (1 1 1)
);
edges
(
arc 2 1 (0 0.010 0)
arc 4 5 (2.800 0.010 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
inlet
{
type patch;
faces
(
(0 2 1 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)
);
}
pipeWall
{
type wall;
faces
(
(1 2 5 4)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //ControlDict File:
I'd like for you to change the values of the directory according to the requirements.
p File:
Changing the values of the pressure file.
p File:
Changing the values of the pressure file.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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 1;
}
front
{
type wedge;
}
back
{
type wedge;
}
pipeWall
{
type zeroGradient;
}
}
// ************************************************************************* //
U File:
Changing the values of the velocity file.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (1 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
pipeWall
{
type noSlip;
}
}
// ************************************************************************* //physical properties:
Changing the values of the physical properties.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object physicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9268e-07;
// ************************************************************************* //
After creating this all files as mentioned above, check the blockMesh.
If blockMesh is okay then check with checkMesh command.
If checkMesh is okay then go on by typing the icoFoam which will create the multiple files upto the time interval 0 to 5 which will take some time for creating it.
If the icoFoam is good then by typing the paraFoam we can can analyze our geometry.
Geometry with Meshing:
As you can see below, there is a Geometry with the meshing of length at the x-axis '0 to 2.8' at the y-axis '0 to 0.009994' and at the z-axis '-0.000349 to 0.000349'.
Angle and diameter as '4 rad and 0.002'.
Velocity U:
This is the velocity contour of the geometry from the inlet to the outlet where the velocity of the geometry is less near the walls and more at farther to the walls.
The velocity and pressure along the profile of the geometry from the inlet to the outlet are shown below:
Post-process Shear stress and validate wall shear stress for fully developed flow:
Shear stress is a deformation of material by slippage along a plane.
Mathematical expression: `tau=mu(du)/(dx)`
Shear stress graph of fully developed flow:
Analytical shear stress graph obtained by MATLAB:
Conclusion:
The fully developed flow velocity profile is obtained at length 1.9 from the inlet of the geometry.
The pressure is decreasing from the inlet to the outlet in the geometry.
By comparing the pressure drop and maximum velocity for fully developed flow we can say that the oulet pressure decreses rapidly.