Simulation Of Flow Through a Pipe - OpenFoam
PROJECT
SIMULATION OF FLOW THROUGH A PIPE
INTRODUCTION
The purpose of this project is to simulate the Laminar Incompressible flow through a pipe. The project is divided into two part , in the first part the objective is to make a simulation. With the given Reynold’s number and Pipe Diameter calculation of inner velocity and the entrance length of the pipe has to be done .
We are considering symmetry condition for this problem else it would be computationally expensive and minimum wedge angle for proper simulation of flow .
Assume that the liquid exhibits laminar flow in a pipe. In the pipe there is bunch of circular layers (lamina) of liquid each having a velocity determined by the radial distance from the center of the tube .
Fluid velocity in a pipe changes from zero at the surface because of the no slip boundary condition to a maximum at the center.
For Laminar Flows the Re <= 2100 .
Entrance Length
The entrance region refers to that portion of a pipe untill the velocity profile is fully developed.
Entrance Number : `El_e=l_e/d`
`l_e` = Length of fully developed velocity profile
`d` = tube or duct diameter
`El_(laminar) = 0.06*R_e`
`l_e/d=0.06*R_e`
`l_e = 0.06*R_e*d`
Hagen Poiseuille Flow :
It applies to Laminar Flow of an incompressible fluid of viscosity through a tube of length and radius .
Flow Rate is directly proportional to the pressure difference , and inversely proportional to the length of the tube and viscosity of the fluid.
Flow Rate : `Q = V_(avg)A_c = (DeltaPpiD^4)/(128muL)`
Average Velocity : `V_(avg) = (DeltaPD^2)/(32muL)`
Assumed Properties :
Given Reynold number : 2100
Diameter : 0.02 m
Fluid temperature : 25 degree Celsius
Entry Length: `l_e = 0.06*R_e*d` = 2.52m
Density of Water : `rho = 997.0 kgm^-3`
Dynamic viscosity : mu = 0.89*10^-3 Pa.S
Kinematic Viscosity : `nu = 8.927*10^-7 m^2s^-1`
Analytical Solution:
Average velocity :`V_(avg) = (Re*mu)/(rho*D) = 0.09373 ms^-1`
Maximum Velocity : `V_(max) = 2*V_(avg) = 0.1875 ms^-1`
Pressure Drop : `DeltaP = (32muLV_(avg))/D^2 = 16.817 Pa`
Kinematic Pressure Drop : `DeltaP_(kinematic) = (DeltaP)/rho = 0.01687 m`
GEOMETRY AND MESHING :
This Project requires the use of MATLAB to code and generate the blockMeshDict file for OpenFoam.
MATLAB Script:
MatLab Code to automatically generate blockMeshDict file for any given wedge angle:
Wedge angle : 4
clear all
close all
clc
% Input
theta = input('Enter the angle in angles '); % Angle of wedge
Re=2100;
D=0.02;
L = 0.06*Re*D; % Length of the pipe
R=D/2;
mu = 0.89e-3; % Viscosity of water at 25 deg C
rho = 997.0; % Density of water at 25 deg C
nu = mu/rho; % Kinematic Viscosity of water at 25 deg C
% Calculation of pressure and velocity:
v_avg = Re*mu/(D*rho); % Average Velocity
v_max = 2*v_avg; % Max Velocity
del_P = (32*mu*v_avg*L)/(D^2); % Pressure Drop
kin_P = del_P/rho; % KinematicPressure Drop
n1 = 50; %Number of cells along the pipe
n2 = 1; %Number of cells radially
n3 = 300; %Number of cells on the curve
gfactor=0.2;
%Vertices of the wedge
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)];
% Printing the file
text1='/*--------------------------------*- C++ -*----------------------------------*\';
text2='| ========= | |';
text3='| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |';
text4='| \\ / O peration | Version: 4.1 |';
text5='| \\ / A nd | Web: www.OpenFOAM.org |';
text6='| \\/ M anipulation | |';
text7='\*---------------------------------------------------------------------------*/';
text8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
%Creating file
f1 = fopen('blockMeshDict.txt','wt');
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\n',text7);
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',text8);
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)',n1,n2,n3);
fprintf(f1,' simpleGrading (%0.1f 1 1)\n);\n\n',gfactor);
fprintf(f1,'edges\n(\n\tarc %d %d (%d %0.2f %d)\n\tarc %d %d (%0.2f %0.2f %d)\n);\n\n', 1, 2, 0, R, 0, 4, 5, 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','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','inlet','patch',0,1,2,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','pipe','wall',1,4,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','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','front','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','back','wedge',0,1,4,3);
fprintf(f1,');\n\n');
fprintf(f1,'mergePatchPairs\n(\n);\n\n');
fprintf(f1,'%s',text8);
fclose(f1);
Geometry:
Fig: Axi-symmetric geometry using the wedge patch type
Main Geometric types :
• patch: generic type containing no geometric or topological information about the mesh, e.g. used for an inlet or an outlet.
• wall: for patch that coincides with a solid wall, required for some physical modelling, e.g. wall functions in turbulence modelling.
• empty: dimension for which no solution is required.
• wedge: for 2 dimensional axi-symmetric cases, e.g. a cylinder, the geometry is specified as a wedge of small angle (e.g. 1◦) and 1 cell thick, running along the centre line, straddling one of the coordinate planes, as shown in figure the axi-symmetric wedge planes must be specified as separate patches of wedge type.
BlockMeshDict File generated for Wedge angle 4:
/*--------------------------------*- 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 9.993908e-03 -3.489950e-04)
(0 9.993908e-03 3.489950e-04)
(2.520000e+00 0 0)
(2.520000e+00 9.993908e-03 -3.489950e-04)
(2.520000e+00 9.993908e-03 3.489950e-04)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 300) simpleGrading (0.2 1 1)
);
edges
(
arc 1 2 (0 0.01 0)
arc 4 5 (2.52 0.01 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Mesh Obtained:
Boundary and Initial Conditions :
Every patch includes a type entry that specifies the type of boundary condition.
Following is the conditions stated for Pressure and Velocity:
| axis | pipe | inlet | oulet | front | back | |
| Pressure (p) | empty | zeroGradient | zeroGradient | fixedValue | wedge | wedge |
| uniform(0.0168682) | ||||||
| Velocity | empty | fixedValue | fixedValue | |||
| uniform (0 0 0) | uniform (0.0937312 0 0) | zeroGradient | wedge | wedge |
Pressure Boundary Condition
/*--------------------------------*- 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 volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
axis
{
type empty;
}
pipe
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0168682;
}
front
{
type wedge;
}
backVelocity Boundary Condition
/*--------------------------------*- 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 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);
}
pipe
{
type fixedValue;
value uniform (0 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //Control
/*--------------------------------*- 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 "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10;
deltaT 0.01;
writeControl timeStep;
writeInterval 10;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Post Processing:
Velocity Contour:
l Profile at INLET
l Fully Developed Profile
Velocity Profile:
l At INLET
l At MIDDLE
l At OUTLET (Fully Developed)
As it can be seen that the velocity profile is a branch of a Parabola. The analytical velocity is 0.1875 and we have obtained 0.1647 with some relative errors.
Pressure Contour:
Pressure Drop:
From the graph it can be seen that that the pressure is non linear at the the beginning and linear onwards .This is because at the beginning the flow is non- developed while later it is developed . The pressure drop analytical is 0.1682 and obtained is 0.021 with some relative error.
Wall Shear stress:
To obtain this graph we have to select the patch pipe .
From the graph it can be seen that the wall shear stress is constant at the latter part of the length of the pipe . This is because of the fact that the velocity becomes constant.