Simulation of Flow through a pipe in OpenFoam
JORGE LOZANO
Introduction to CFD using Matlab and OpenFOAM
Simulation of Flow through a pipe in OpenFoam
OBJECTIVE
- Write a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes
- Reynolds number based on pipe diameter and inlet velocity should be 2100
- Working fluid - water
- Calculate the length of the pipe
- Calculate length of the pipe using the entry length formula for laminar flow through a pipe
- Show that entry length is enough to produce a fully developed flow. You will have to compare velocity profiles at different positions along the length of the pipe.
- Solver – icoFoam.
Theory
Hagen-Poiseuille equation
The Poiseuille’s law states that the flow of liquid depends on following factors like the pressure gradient (∆P), the length of the narrow tube (L) of radius (r) and the viscosity of the fluid (η) along with relationship among them.
Hagen—Poiseuille equation Continuum mechanics. To calculate the flow through each lamina, we multiply the velocity from above and the area Poiseuille the lamina. Finally, put this expression in the form of a differential equation dropping the term quadratic in dr. 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.
- The Pressure Gradient (∆P) : Shows the difference in the pressure between the two ends of the tube, determined by the fact any fluid will always flow from high pressure (p1) to low pressure region(p2) and the flow rate is determined by the pressure gradient (P1 – P2)
- Radius of tube: The liquid flow varies directly with the radius to the power 4.
- Viscosity (η): The flow of the fluid varies inversely with the viscosity of the fluid and as the viscosity of the fluid increases, the flow decreases vice versa.
- Length of the Tube (L): The liquid flow is inversely proportional to the length of the tube, therefore longer the tube, greater is the resistance to the flow.
- Resistance(R): The resistance is described by 8Ln / πr4 and therefore the Poiseuille’s law becomes
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.
The analytical calculations to obtain the optimum initial conditions as well as the geometry are:
- Fluid: Water @60 C
- Reynolds number: 2100
- Density:983.2 kg/m^3
- Kinematic viscosity:0.475e^-6
- Diameter: 0.127m
- Solver:icoFoam
Calculations and previous results:
- Le=0.06*Re*D=16.002m
- L= Le+0.3=16.302m
- Vavg= (Re*mu)/(rho*D) = 7.84559 m/s
- Vmax=2*Vavg=15.6912 m/s
- P1-P2=1183.75118 Pa
- Tao=(2*mu*Vmax)/(R^2) = 0.0017636 N/m^2
The flow is usually expressed at outlet pressure. It follows that the resistance R is proportional to the length L of the resistor, which is true. The Hagen—Poiseuille equation is useful in determining the flow rate of intravenous fluids that may be achieved using various sizes of peripheral and central cannulas.
As commented before, the velocity along the pipe will obtain a parabola like graph when having a fully develop water flow.
Procedure
To apply the wedge pipe-flow equation several ico-Foam simulations were processed to see the effect of the shear stress, velocity variation, and the efficiency of the angle applied to the pre-processor.
By obtaining from analytical equations the initial conditions as well as the geometry requirements, by using MATLAB a blockMeshDict code was created being able to change several values obtain a pipe-section of a desired angle and length.
The following code shows the blockMeshDict file, created using MATLAB and the result obtained, knowing that the first simulation will be applying a 1 deg wedge pipe:
Matlab blockMeshDict code:
clear all
close all
clc
%%Input data for analysis of flow
%Diameter of the pipe in m
D=0.127;
%Radius of the pipe in m
R=D/2;
%Reynolds number
Re=2100;
%Entrance Length in m
L_e=0.06*Re*D;
%Total length of the pipe in m
L=0.3+L_e;
%wedge angle in degree
theta=1;
%Dynamic viscosity of water at 60 deg. Cel. in Ns/n^2
mu=0.0004660;
%Density of water at 60 deg. Cel. in kg/m^3
rho=983.2;
%kinematic viscosity of water at 60 deg. Cel. in m^2/s
nu=mu/rho;
%%calculating the solution according to Hagen-Poiseuille\'s flow equation
%Average velocity
V=(Re*mu)/(rho*D)
%Maximum velocity
u_max=2*V
%pressure gradient
del_pdel_x=(8*mu*V)/(R^2)
%Pressure drop
delta_p=(8*mu*V*L)/(R^2)
%kinematic pressure
kin_pressure=delta_p/rho
%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)];
f1=fopen(\'BlockMeshDict.txt\',\'w\');
fprintf(f1,\'%s\\n\',\'/*------------------------------*- C++ -*-----------------------*\\\');
fprintf(f1,\'%s\\n\',\'============ | \');
fprintf(f1,\'%s\\n\',\'\\\\ / | \');
fprintf(f1,\'%s\\n\',\' \\\\ / F ield | OpenFOAM:The Open Source CFD Toolbox\');
fprintf(f1,\'%s\\n\',\' \\\\ / O peration | Website: https://openfoam.org \');
fprintf(f1,\'%s\\n\',\' \\\\ / A nd | Version:7 \');
fprintf(f1,\'%s\\n\',\' \\\\/ M anipulation | \');
fprintf(f1,\'%s\\n\',\'\\*--------------------------------------------------------------*/\');
fprintf(f1,\'%s\\n\',\'FoamFile\');
fprintf(f1,\'%s\\n\',\'{\');
fprintf(f1,\'%s\\n\',\' version 2.0;\');
fprintf(f1,\'%s\\n\',\' format ascii;\');
fprintf(f1,\'%s\\n\',\' class dictionary;\');
fprintf(f1,\'%s\\n\',\' object blockMeshDict;\');
fprintf(f1,\'%s\\n\',\'}\');
fprintf(f1,\'%s\\n\\n\',\'// ********************************************************** //\');
fprintf(f1,\'%s\\n\\n\',\'convertToMeters 1;\');
%Defining the vertices of the geometry
fprintf(f1,\'%s\\n\',\'vertices\');
fprintf(f1,\'%s\\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\',v5(1),v5(2),v5(3));
fprintf(f1,\'%s\\n\\n\',\');\');
%Defining the blocks
fprintf(f1,\'%s\\n\',\'blocks\');
fprintf(f1,\'%s\\n\',\'(\');
fprintf(f1,\'\\thex (0 1 2 0 3 4 5 3)(%d %d %d) simpleGrading (0.1 1 1)\\n\',50,1,400);
fprintf(f1,\'%s\\n\\n\',\');\');
%Defining the edges
fprintf(f1,\'%s\\n\',\'edges\');
fprintf(f1,\'%s\\n\',\'(\');
fprintf(f1,\'\\tarc %d %d (%d %d %d)\\n\',1,2,0,R,0);
fprintf(f1,\'\\tarc %d %d (%d %d %d)\\n\',4,5,L,R,0);
fprintf(f1,\'%s\\n\',\');\');
%Defining the boundaries
fprintf(f1,\'%s\\n\',\'boundary\');
fprintf(f1,\'%s\\n\',\'(\');
fprintf(f1,\'\\tinlet\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype patch;\\n\');
fprintf(f1,\'\\t\\tfaces \\n\\t\\t(\\n\');
fprintf(f1,\'\\t\\t\\t(0 2 1 0)\\n\');
fprintf(f1,\'\\t\\t);\\n\\t}\\n\');
fprintf(f1,\'\\toutlet\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype patch;\\n\');
fprintf(f1,\'\\t\\tfaces \\n\\t\\t(\\n\');
fprintf(f1,\'\\t\\t\\t(3 4 5 3)\\n\');
fprintf(f1,\'\\t\\t);\\n\\t}\\n\');
fprintf(f1,\'\\tfront\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype wedge;\\n\');
fprintf(f1,\'\\t\\tfaces \\n\\t\\t(\\n\');
fprintf(f1,\'\\t\\t\\t(0 3 5 2)\\n\');
fprintf(f1,\'\\t\\t);\\n\\t}\\n\');
fprintf(f1,\'\\tback\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype wedge;\\n\');
fprintf(f1,\'\\t\\tfaces \\n\\t\\t(\\n\');
fprintf(f1,\'\\t\\t\\t(0 1 4 3)\\n\');
fprintf(f1,\'\\t\\t);\\n\\t}\\n\');
fprintf(f1,\'\\ttop\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype wall;\\n\');
fprintf(f1,\'\\t\\tfaces \\n\\t\\t(\\n\');
fprintf(f1,\'\\t\\t\\t(1 2 5 4)\\n\');
fprintf(f1,\'\\t\\t);\\n\\t}\\n\');
fprintf(f1,\'\\taxis\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype empty;\\n\');
fprintf(f1,\'\\t\\tfaces \\n\\t\\t(\\n\');
fprintf(f1,\'\\t\\t\\t(0 3 3 0)\\n\');
fprintf(f1,\'\\t\\t);\\n\\t}\\n\');
fprintf(f1,\');\\n\\n\');
fprintf(f1,\'mergePatchPairs\\n(\\n);\\n\\n\');
fprintf(f1,\'// ******************************************************************* //\');
fclose(f1)blockMeshDict code generated:
/*------------------------------*- 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;
object blockMeshDict;
}
// ********************************************************** //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.0634976 -0.000554135)
(0 0.0634976 0.000554135)
(16.302 0 0)
(16.302 0.0634976 -0.000554135)
(16.302 0.0634976 0.000554135)
);
blocks
(
hex (0 1 2 0 3 4 5 3)(50 1 400) simpleGrading (0.1 1 1)
);
edges
(
arc 1 2 (0 0.0635 0)
arc 4 5 (16.302 0.0635 0)
);
boundary
(
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)
);
}
top
{
type wall;
faces
(
(1 2 5 4)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ******************************************************************* //The following codes are some of the important files used to select the initial conditions to perform an incompressible simulation.
Pressure initial condition “0”:
/*--------------------------------*- 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
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 1183.75118;
}
front
{
type wedge;
}
back
{
type wedge;
}
top
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Velocity initial condition “0”:
/*--------------------------------*- 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 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (7.8456 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
top
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //
ControlDict initial simulation setup:
/*--------------------------------*- 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 latestTime;
latestTime 5;
stopAt endTime;
endTime 15;
deltaT 0.002;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
The following code is to perfrom the simulation while using four processors:
/*--------------------------------*- 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 decomposeParDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
numberOfSubdomains 4;
method hierarchical;
hierarchicalCoeffs
{
n (4 1 1);
delta 0.001;
order xyz;
}
distributed yes;
// ************************************************************************* //
Results
By performing a transient simulation of 5 seconds where the flow fully develops, obtaining a courant average number of 0.194, which depends on the time step selected, the grid size and the mean velocity of the flow.
The following picture represents the last iterations performed on the simulation.
Last simulation iterations performed:
Partial mesh wedge pipe:
Pressure Contour:
Velocity Contour:
Reynolds number through the overall pipe length:
Shear Stress through the radial length:
The following pictures present the flow development over the length of the wedge-pipe, maintaining the mass flow:
Velocity development 0 m from inlet (RESULT):
Velocity development 0.1 m from inlet (RESULT):
Velocity development 1 m from inlet (RESULT):
Velocity development 16.302 m from inlet (RESULT):
The following cases represents the same simulation setup as before although it’s modified on the next conditions:
Total time: 15 seconds
Wedge pipe angle to test: 10deg., 25deg. and 45deg.
The following figures represent the block mesh for each angle: 10, 25 and 45 respectively:
The following figures represent the final code lines with their respective computational time for each: 10, 25 and 45 respectively:
Having as expected, by having a greater angle it results in more computational time, and greater inconsistencies in the convergence of the results.
By comparing the results from the velocity over the radial length of the 10, 25 and 45 degrees according to the 1-degree simulation which result in a more accurate solution, the next two graphs are obtain, being the second one a shear stress error precision graph
Error precision velocity vs radial length graph
Error precision shear stress vs radial length graph
Conclusion
Several CFD simulations were performed by using the icoFoam solver from OpenFOAM,, applying the Hagen-Poiseuille flow equations to determine several initial conditional and important system values, being able to see through the development of the pipe the flow development reaching the maximum velocity (maintaining the mass flow through the wedge pipe) as well as the shear stress values at the end of the wedge pipe generated being able to code a mesh by using MATLAB and comparing the error efficiency for the velocity and shear stresses within a radial distance a the final distance along the other angles analyzed: 10, 25 and 45 degrees. The shear stress obtain at the surface with the closer cells are as calculated initially with a certain amount of error precision as well as the velocities (average and maximum) due to the non- slip condition (at the top wall) and the conservation of mass and momentum according to the Navier-stokes equations applied. According to the efficiency of the simulations it isn’t as precise as expected although it results in close values as a practical application due to the fast-computational time taken for all cases.