Week 11 - Simulation of Flow through a pipe in OpenFoam
Simulation of Laminar Viscous flow through a pipe in OpenFOAM
Aim : To simulate the Laminar Vicous flow through a pipe in OpenFOAM for a sector(wedge angle `theta` < `5^0`) of pipe.
Objectives :
- To generate a blockMeshDict file using MATLAB that can simulate a pipeflow for given wedge angle.
- Compare the CFD results with the velocity and pressure calculated from the Hagen Poiseullie's equation.
- To postprocess the results obtained and plot the pressure and Velocity contours.
Theory :
The Hagen–Poiseuille equation[1], also known as the Hagen–Poiseuille law, Poiseuille law or Poiseuille equation, gives the pressure drop in an incompressible and newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section.
Poiseuille's Equation describes the pressure drop due to the viscosity of the fluid
Fig.1 : A Laminar flow through pipe
Assumptions in the Hagen-Poiseulie Law :
1. The flow is Laminar
2. The fluid is an incompressible and newtonian fluid (a fluid that obeys Newton's law of viscosity `tau = mu*((du)/dx)`)
3. There is no acceleration of the fluid in the pipe flow.
Hydrodynamic entrance length or entrance length : This is distance a fluid flow travels after entering a pipe, before it becomes fully developed.
Entrance length refers to the length of the entry region, the area following the pipe entrance where effects originating from the interior wall of the pipe propagate into the flow as an expanding boundary length. When the boundary layer expands to fill the entire pipe, the developing flow becomes a fully developed flow, where flow characteristics no longer change with increased distance along the pipe.
Mathematically,
Entrance length of a pipe is given by L = `0.06*Re*D`
Here, D is the diameter of the pipe.
It is evident from the above formula that the Hydrodynamic entrance length is not the same in every case. As the Reynolds number (Re) varies the flow characteristics vary and hence the entrance length.
The entire simulation is carried for the following fixed quantities :
Diameter of the pipe = 0.02m.
Properties of the Water at 25 degree centigrade :
Density = 997 kg/m3.
Dynamic viscosity = 8.9 x 10^(-4) N.s/m^2.
Kinematic viscosity = Dynamic viscosity/Density = 8.9 x 10^(-7) m.s^(-2).Calculations :
1. Calculation of the Hydrodynamic entrance length :
1. For a Reynolds number(Re) of 2100, the flow is laminar in case of a pipe flow.
2. Hydrodynamic entrance length (L_e) = 0.06 x Re x D
= 0.06 x 2100 x 0.02
Therefore L_e = 2.52 m.
This means the minimum length of the pipe should be 2.52 m.
The specific reason being the Hagen Poiseullie's law gives the pressure drop values for a fully developed flow.
In our case the fully developed flow exists for a length greater than 2.52 m.
So let us consider Length of pipe as 2.7 m.
2. calculation of the Flow velocity from the Reynolds number fromula
In our case, the flow is laminar. So let the Re be 2100.
Average velocity (V_avg) = Re. μ /(ρ.D)
= 2100.0.89 x 10^(-3) / (1000 x 0.02)
= 0.09345 m/s
Maximum velocity (V_max) = 2 x Average velocity
= 2 x 0.09345
= 0.1869 m/s.
3. Calculation of Pressure drop from the Hagen Poiseullie's equation
Pressure drop as given by the Hagen poiseullie's formula is
ΔP = 32 x μ x L x V_avg / (D^2)
= 32 x 0.89 x 10^(-3) x 2.7 x 0.09345 /(0.02 x 0.02)
= 18.01 Pa.
Kinematic Pressure drop (Δp) = (ΔP)/(ρ) = 0.01807 mGeneration of Mesh using the blockMeshDict :
The following code is written in MATLAB editor to generate the blockMeshDict file
clear; close; clc
d = 0.02; % diameter of the pipe
r = d/2;
l = 2.7;
theta = input('Enter the wedge angle = ');
x = input('Enter the simple grading factor along X-axis:');
y = input('Enter the simple grading factor along Y-axis:');
z = input('Enter the simple grading factor along Z-axis:');
% A program to parse the data from the blockMeshDict file to another
% blockMeshDict file
fid = fopen('blockMeshDict.txt','w');
fprintf(fid,'/*--------------------------------*- \n');
fprintf(fid,'C++ -*----------------------------------*\ \n');
fprintf(fid,' ========= | \n');
fprintf(fid,' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \n');
fprintf(fid,' \\ / O peration | Website: https://openfoam.org \n');
fprintf(fid,' \\ / A nd | Version: 7 \n');
fprintf(fid,' \\/ M anipulation | \n');
fprintf(fid,'\*-------------------------------------------------------------- \n');
fprintf(fid,'-------------*/ \n');
fprintf(fid,'FoamFile \n');
fprintf(fid,'{ \n');
fprintf(fid,' version 2.0; \n');
fprintf(fid,' format ascii; \n');
fprintf(fid,' class dictionary; \n');
fprintf(fid,' object blockMeshDict; \n');
fprintf(fid,'} \n');
fprintf(fid,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * \n');
fprintf(fid,'* * * * * * // \n');
fprintf(fid,'convertToMeters 1; \n');
fprintf(fid,'\n');
fprintf(fid,'vertices \n');
fprintf(fid,'( \n');
fprintf(fid,'(%d %d %d)\n',0,0,0); %index 0
fprintf(fid,'(%d %d %d) \n', 0, r*cosd(theta/2), r*sind(theta/2)); %index 1
fprintf(fid,'(%d %d %d) \n', 0, r*cosd(theta/2), -r*sind(theta/2)); %index 2
fprintf(fid,'(%d %d %d) \n',l, 0 ,0); %index 3
fprintf(fid,'(%d %d %d) \n',l , r*cosd(theta/2), r*sind(theta/2)); % index 4
fprintf(fid,'(%d %d %d) \n', l, r*cosd(theta/2), -r*sind(theta/2)); % index 5
fprintf(fid,'); \n');
fprintf(fid,'blocks \n');
fprintf(fid,'( \n');
fprintf(fid,'hex (0 3 5 2 0 3 4 1) (600 40 1) simpleGrading (%d %d %d)\n',x,y,z');
fprintf(fid,');\n');
fprintf(fid,'edges \n');
fprintf(fid,'( \n');
fprintf(fid,'arc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(fid,'arc %d %d (%d %d %d)\n',4,5,l,r,0);
fprintf(fid,'); \n');
fprintf(fid,'\n');
fprintf(fid,'boundary \n');
fprintf(fid,'( \n');
fprintf(fid,'inlet \n'); %Defining the inlet faces
fprintf(fid,'{ \n');
fprintf(fid,'type patch;\n');
fprintf(fid,'faces \n');
fprintf(fid,'( \n');
fprintf(fid,'(0 1 2 0)\n');
fprintf(fid,'); \n');
fprintf(fid,'}\n');
fprintf(fid,'outlet \n'); % Defining the outlet faces
fprintf(fid,'{ \n');
fprintf(fid,'type patch ;\n');
fprintf(fid,' faces \n');
fprintf(fid,'( \n');
fprintf(fid,'(3 4 5 3)\n');
fprintf(fid,'); \n');
fprintf(fid,'}\n');
fprintf(fid,'top\n'); %Defining the top wall
fprintf(fid,'{\n');
fprintf(fid,'type wall;\n');
fprintf(fid,'faces\n');
fprintf(fid,'( \n');
fprintf(fid,'(1 4 5 2)\n');
fprintf(fid,'); \n');
fprintf(fid,'} \n');
fprintf(fid,'back \n'); % Defining the back wall
fprintf(fid,'{\n');
fprintf(fid,'type wedge; \n');
fprintf(fid,'faces \n');
fprintf(fid,'( \n');
fprintf(fid,'(0 3 5 2)\n');
fprintf(fid,');\n');
fprintf(fid,'}\n');
fprintf(fid,'front \n'); %Defining the front wall
fprintf(fid,'{ \n');
fprintf(fid,'type wedge;\n');
fprintf(fid,'faces \n');
fprintf(fid,'( \n');
fprintf(fid,'(0 3 4 1)\n');
fprintf(fid,');\n');
fprintf(fid,'} \n');
fprintf(fid,'axis\n'); %Defining the axis line (symmetry line)
fprintf(fid,'{ \n');
fprintf(fid,'(');
fprintf(fid,'type empty;\n');
fprintf(fid,'( \n');
fprintf(fid,'(0 3 3 0)\n');
fprintf(fid,');\n');
fprintf(fid,'} \n');
fprintf(fid,'); \n'); % closing the boundary bracket
fprintf(fid,'mergePatchPairs \n');
fprintf(fid,'( \n');
fprintf(fid,'); \n');
fprintf(fid,'// \n');
fprintf(fid,'// ****************************************************************\n');
fprintf(fid,' // ********* // \n');
The blockMeshDict file is located inside the directory - system.
Setting up the initial conditions :
The directory 0 has the files for pressure and the velocity which contains the pressure and velocity values at the initial conditions
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
{
top
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
bottomwall
{
type empty;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.01807;
}
}
// ************************************************************************* //
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 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.1869 0 0);
}
front
{
type wedge;
}
back
{
type wedge;
}
bottomwall
{
type empty;
}
top
{
type noSlip;
}
outlet
{
type zeroGradient;
}
}
// ************************************************************************* //
The transport properties file has the value for the kinematic viscosity which looks 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;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 1e-6;
// ************************************************************************* //
The geometry looks as follows :
The following are the mesh parameters :
number of points : 48681
number of cells : 24000
number of faces : 96040
no. of internal faces : 47360Results :
1. Veocity contour :
It can be observed that a velocity gradient exists at the interface of the fluid and solid. The no slip condition ensures that the fluid velocity is zero at the edge.
As an established fact, in the fully developed flow region the velocity profile is parabolic.
As evident from the above results, the pressure difference(`Deltap`) is around 0.01838 which is very close to the value obtained from the Hagen Poiseulle's equation.
Shear stress :
Conclusion : A laminar viscous pipe flow is simulated in the OpenFOAM and the values of the pressure drop from the numerical analysis and hagen Poiseullie's equation are compared.
References :
1. Hagen Poiseullie's Law : Wikipedia - https://en.wikipedia.org/wiki/Hagen%E2%80%93Poiseuille_equation
2. Hydrodynamic entrance length : Wikipedia - https://en.wikipedia.org/wiki/Entrance_length_(fluid_dynamics)