Simulation of a fluid flow through pipe
Simulation of fluid flow through the pipe using MATLAB & OpenFoam -PART 1 of 2
Objective of the project :
- To simulate for internal Laminar flow through the pipe and compute the velocity profiles & Shear stresses at various cross cestions .
- To Set up the sufficient length of the pipe so that the flow is fully developed (i.e) axial component of the velocity remains unchanged
- To correlate the velocity proifle obtained from the simulations with the Hagen poiseuille equation (Equation for u(r))
Scope of the Project:
- Writing a code in MATLAB for the automatic generation of the blockmesh file for OpenFoam through File parsing technique
- To take the advantage of the rotational symmetry of the pipe and also to reduce the computation time , the axisymmetric model is developed (i.e) different wedge angles is modelled and simulations are carried out for those models
- The codes are written in MATLAB for each models by varying the angle in the program that automatically defines the cylindrical co ordinates (r , theta ,z) for the models accordincly
- The simulations are performed in openFoam with appropriate initial and Boundary conditions imposed and all the necessary flow quantities are computed .
Assumptions :
- Steady flow , incompressible and laminar flow
- Flow is two dimenional ( Flow WRT theta is ignored )
- Constant fluid properties (Viscosity , Density etc)
Introduction - Laminar flow through pipes :
As we all aware that the fluid in contact with the surface , it experience a momentum disturbances such fluid particles in the layer in contact with the surface of the pipe come to a complete stop. This layer also causes the fluid particles in the adjacent layers to slow down gradually as a result of friction. To make up for this velocity reduction, the velocity of the fluid at the midsection of the pipe has to increase to keep the mass flow rate through the
pipe constant. As a result, a velocity gradient develops along the pipe.The region of the flow in which the effects of the viscous shearing forcescaused by fluid viscosity are felt is called the velocity boundary layer or just the boundary layer. The hypothetical boundary surface divides the flow in a pipe into two regions: the boundary layer region, in which the viscous effects and the velocity changes are significant, and the irrotational (core) flow region, in which the frictional effects are negligible and the velocity remains essentially constant in the radial direction. The thickness of this boundary layer increases in the flow direction untilthe boundary layer reaches the pipe center and thus fills the entire pipe The region from the pipe inlet to the point at which the boundary layer merges at the centerline is called the hydrodynamic entrance region, and the length of this region is called the hydrodynamic entry length Lh. Flow in the entrance region is called hydrodynamically developing flow since this is the region where the velocity profile develops.The region beyond the entrance region in which the velocity profile is fullydeveloped and remains unchanged is called the hydrodynamically fully developed region.
It is also to be noted that the velpocity proifle for a fully devloped flow i case of a laminar flow will be parabolic as shown above whereas , for turbulent flows the profile will be fuller and also flatter .
Formulation for entry length
The hydrodynamic entry length is usually taken to be the distance from thepipe entrance to where the wall shear stress (and thus the friction factor)reaches within about 2 percent of the fully developed value. In laminar flow, the hydrodynamic entry length is given approximately as
L _Entry = 0.06* RE * D
Problem description :
As far as the flow through the pipes are concerned , it might sound very simple but while simulating these kind of flow problems , there are so many things that we need to consider .Just a few to list : How long the pipe should be so as to develop a fully debeloped flow and what are the fluid initial and boundary conditions . Our challenge here is to find the values of quantities and thus boundary conditions based on the assumptions under which Hagen poiseuille equation is derived.
Modelling Inputs :
- Re = 2100 (essentially laminar)
- Working fluid water : Kinematic viscosity = 1e-06 m^2/s (at 20 deg )
- Dynamic viscosity of the water : `mu = 0.89*10^-3 Pa.S`
- Density of water = 997 Kg/m^3
- Diameter of the pipe D = 0.02 m (Assuming )
- The average velocity for the corresponding values of Re and D is found out from the expressions
`V_(avg) = (Re*mu)/(rho*D) = 0.09373 ms^-1`
- The maximum velocity is computed from the expression
`V_(max) = 2*V_(avg) = 0.1875 ms^-1`
- The pressure drop across the pipe is computed from the Hagen poiseuille equations as follows
`DeltaP = (32muLV_(avg))/D^2 = 18.2 Pa`
- Since openFoam deals with the kinematic quantites , the kinematic pressure drop is calculated as follows
`DeltaP_(kinematic) = (DeltaP)/rho = 0.018 m`
The entry length for the above inputs can be calculated as L_entry = 0.06 * 2100*0.02
= 2.52 m.Thus an extra length of 0.3 m is added for the simulation to setup the fully developed flow .
The velocity profile u(r) for a fully developed laminar flow in a pipe can be obtained through hagen poiseuille derivation as follows
`u(r) = R^2/(4mu)*((dp)/dx)(1-r^2/R^2)`
The average velocity can also be written as `V_(avg) = R^2/(8mu)*((dp)/dx)`
Hence the final expressions for the velocity profile will be
`u(r) = 2V_(avg)(1-r^2/R^2)`
Program for Wedge Boundary conditions :
It is to be noted that , For 2 dimensional axi-symmetric cases , e.g. a cylinder, the geometry is specified as a wedge of small angle (usually less than 5 degrees) .Here the flow is simulated for two cases ( wedge angle 3 and 4 respectively )
Code snippet ::::Wedge angle 3 degrees
% Solving Hagen poiseuille flow problem inside a pipe
% Inputs
clear all
close all
clc
Re = 2100;
rho = 1000 ; % Unit in Kg/m3;
D = 0.02;% Unit in m;
L_ent = 0.06*Re*D;
R =D/2;
mu = 0.00089;%viscosity of water at 25 degrees .unit in Pa.s
nu =mu/rho;%Unit in m2/s
theta = 3; %determines the wedge angle for geometry
% Solving for boundary conditions
v_avg =(Re*nu)/D;
L =L_ent+0.3;
v_max =2*v_avg;
delp = (v_avg*32*mu*L)/(D*D);%Pressure drop according to Hagen poiseuille
Kin_p =delp/rho;% Kinematic pressure drop .Required to impose in openFoam
% Calculating 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)];
grading = 0.2;
b1 = 50;
b2 =1;
b3 =200;
% Formatting openFoam for blockMeshDict
h1=\'/*--------------------------------*- C++ -*----------------------------------*\\\';
h2=\'| ========= | |\';
h3=\'| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |\';
h4=\'| \\\\ / O peration | Version: 4.1 |\';
h5=\'| \\\\ / A nd | Web: www.OpenFOAM.org |\';
h6=\'| \\\\/ M anipulation | |\';
h7=\'\\*---------------------------------------------------------------------------*/\';
h8=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';
f1=fopen(\'blockMeshDict.txt\',\'wt\');
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\\n\',h7);
fprintf(f1,\'FoamFile\\n{\\n\');
fprintf(f1,\'%16s \\t 2.0;\\n\',\'version\');
fprintf(f1,\'%16s \\t ascii;\\n\',\'format\');
fprintf(f1,\'%16s \\t dictionary;\\n\',\'class\');
fprintf(f1,\'%16s \\t blockMeshDict;\\n}\\n\',\'object\');
fprintf(f1,\'%s\\n\\n\',h8);
fprintf(f1,\'convertToMeters 1;\\n\\n\');
fprintf(f1,\'vertices\\n(\\n\');
fprintf(f1,\'\\t(%d %d %d)\\n\',v0(1),v0(2),v0(3));
fprintf(f1,\'\\t(%d %d %d)\\n\',v1(1),v1(2),v1(3));
fprintf(f1,\'\\t(%d %d %d)\\n\',v2(1),v2(2),v2(3));
fprintf(f1,\'\\t(%d %d %d)\\n\',v3(1),v3(2),v3(3));
fprintf(f1,\'\\t(%d %d %d)\\n\',v4(1),v4(2),v4(3));
fprintf(f1,\'\\t(%d %d %d)\\n);\\n\\n\',v5(1),v5(2),v5(3));
fprintf(f1,\'blocks\\n(\\n\');
fprintf(f1,\'\\t hex (0 1 2 0 3 4 5 3) (%d %d %d)\\n\',b1,b2,b3);
fprintf(f1,\'\\t simpleGrading\\n\\t (%d 1 1)\\n);\\n\\n\',grading);
fprintf(f1,\'edges\\n(\\n\');
fprintf(f1,\'\\t arc 1 2 (%d %d %d)\\n\',0,R,0);
fprintf(f1,\'\\t arc 4 5 (%d %d %d)\\n);\\n\\n\',L,R,0);
fprintf(f1,\'boundary\\n(\\n\');
fprintf(f1,\'\\t%s\\n\\t{\\n\\t\\t type %s;\\n\\t\\t faces \\n\\t\\t (\\n\\t\\t\\t (%d %d %d %d) \\n\\t\\t );\\n\\t}\\n\\n\',\'inlet\',\'patch\',0,0,2,1);
fprintf(f1,\'\\t%s\\n\\t{\\n\\t\\t type %s;\\n\\t\\t faces \\n\\t\\t (\\n\\t\\t\\t (%d %d %d %d) \\n\\t\\t );\\n\\t}\\n\\n\',\'outlet\',\'patch\',3,4,5,3);
fprintf(f1,\'\\t%s\\n\\t{\\n\\t\\t type %s;\\n\\t\\t faces \\n\\t\\t (\\n\\t\\t\\t (%d %d %d %d) \\n\\t\\t );\\n\\t}\\n\\n\',\'Wall\',\'wall\',2,5,4,1);
fprintf(f1,\'\\t%s\\n\\t{\\n\\t\\t type %s;\\n\\t\\t faces \\n\\t\\t (\\n\\t\\t\\t (%d %d %d %d) \\n\\t\\t );\\n\\t}\\n\\n\',\'axis\',\'empty\',0,3,3,0);
fprintf(f1,\'\\t%s\\n\\t{\\n\\t\\t type %s;\\n\\t\\t faces \\n\\t\\t (\\n\\t\\t\\t (%d %d %d %d) \\n\\t\\t );\\n\\t}\\n\\n\',\'wedgeFront\',\'wedge\',0,3,5,2);
fprintf(f1,\'\\t%s\\n\\t{\\n\\t\\t type %s;\\n\\t\\t faces \\n\\t\\t (\\n\\t\\t\\t (%d %d %d %d) \\n\\t\\t );\\n\\t}\\n\\n\',\'wedgeBack\',\'wedge\',0,1,4,3);
fprintf(f1,\');\\n\\n\');
fprintf(f1,\'mergePatchPairs\\n(\\n);\\n\\n\');
fprintf(f1,\'%s\',h8);
fclose(f1);
The following file is an automatic generated file for wedge angle 3 which is derived from the above MATLAB code .This file is used in place of the existing blockMeshDict file which is in the incompressible solver file from openFoam software .Note that the solver used here is icoFoam .The same technique is used for wedge angle 4 as well. The only difference is specifying the angle as 4 degrees in matlab code instead of 3 degrees
/*--------------------------------*- C++ -*----------------------------------*\\
| ========= | |
| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\\\ / O peration | Version: 4.1 |
| \\\\ / A nd | Web: www.OpenFOAM.org |
| \\\\/ M anipulation | |
\\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.00999657 -0.000261769)
(0 0.00999657 0.000261769)
(2.82 0 0)
(2.82 0.00999657 -0.000261769)
(2.82 0.00999657 -0.000261769)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading(0.2 1 1)
);
edges
(
arc 1 2 (0 0.01 0)
arc 4 5 (2.82 0.01 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 0 2 1)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
Wall
{
type wall;
faces
(
(2 5 4 1)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wedgeFront
{
type wedge;
faces
(
(0 3 5 2)
);
}
wedgeBack
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Initial conditions ::: Kinematic pressure drop
/*--------------------------------*- 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 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 0.018;
}
wall
{
type zeroGradient;
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //Initial conditions ::: Inlet Average velocity
/*--------------------------------*- 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 volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.0937312 0 0);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (0.0937312 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //Transport Properties :::: Kinematic viscosity
/*--------------------------------*- C++ -*----------------------------------*\\
| ========= | |
| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\\\ / O peration | Version: 4.1 |
| \\\\ / A nd | Web: www.OpenFOAM.org |
| \\\\/ M anipulation | |
\\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location \"constant\";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 1e-06;
// ************************************************************************* //
Time step controls ::
/*--------------------------------*- 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;
location \"system\";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 50;
deltaT 0.01;
writeControl timeStep;
writeInterval 1000;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Results & Discussions ::
The computational domain with respect to above mentioned mesh grading for the wedge angle 3 is shown below .The simulation is ran for 50 seconds.
Velocity contour near the entrance for Wedge angle 3::
Velocity contour near the exit for Wedge angle 3::
A) Velocity profile at location 0.1
B) Velocity profile at location 0.5
C) Velocity profile at location 1
D) Velocity profile at location 2
E) Velocity profile at location 2.5
F) Velocity profile at location 2.8
It is evident from the results that after a certain distance the parabolic velocity profile is formed and this profile remains unchanged in the axial direction which in turn tells us that the entry length we set up with respect to the conditions inside the domain is absolutely sufficient for obtaining the fully developed flow regime .Also shear stress for the fully developed flow is also plotted below.
The simulation is performed for another wedge angle 4 in a similiar way and the results are compared with each other respectively.