Simulation of Flow through a pipe in OpenFoam Part 1
AIM
The aim of this project is to write a matlab program to automatically generate blockmeshDict file for any wedge angle and mesh grading factor
OBJECTIVES
The objectives of this project are
- Simulation of the laminar flow over a pipe
- post processing of velocity and shear stress profiles
LAMINAR FLOW THROUGH PIPE
Laminar flow is defined by fluid particles following smooth paths in layers, with each layer moving smoothly past the adjacent layers with little or no mixing. In laminar flow, the motion of the particles of the fluid is very orderly with particles close to a solid surface moving in straight lines parallel to that surface. The characteristics of laminar flow are
- 'No slip at the boundary.
- Due to viscosity, there is a shear stress between the fluid layers.
- The flow is rotational.
- Due to viscous shear, there is continous dissipation of energy.
- No mixing between the fluid layers.
The hydrodynamic entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces propagating from the interior wall of a pipe. This region is characterized by a non-uniform flow. The fluid enters a pipe at a uniform velocity, then fluid particles in the layer in contact with the surface of the pipe come to a complete stop due to the no slip condition. Due to viscous forces within the fluid, the layer in contact with the pipe surface resists the motion of adjacent layers and slows adjacent layers of fluid down gradually, forming a velocity profile. For the conservation of mass to hold true, the velocity of layers of the fluid in the center of the pipe increases to compensate for the reduced velocities of the layers of fluid near the pipe surface. This develops a velocity gradient across the cross section of the pipe.Reynolds Number is defined as the ratio of inertia force to the viscous force.
`R_e=frac[rhoV_(avg)D][mu]`OR `R_e=frac[V_(avg)D][nu]`
The following calculations are done for the analytical solution of the Hagen Poiseuille law
- Reynolds Number = 2100
- Kinematic viscosity `(nu)`= 1e-6 `m^2/s`
- Dynamic Viscosity `(mu) = 0.00089 `(N-s)/m^2`
- Density of Water = 1000 `(kg)/m^3`
- Diameter of pipe = 20mm or 0.02m
- Hydrodynamic Length `L_eapprox0.06R_eD` `=2.52m`
- Average velocity `V_(avg)=frac[R_enu}[D]``= 0.09345 m/s`
- Maximum velocity `V_(max)=2V_(avg)``=0.1869m/s`
- Pressure `Deltap=frac[8V_(avg)muL][R^2]=18.763N/m^2`
Now in OpenFoam we have to provide the pressure boundary conditions in terms of kinematic pressure
`P_k=frac[DeltaP][rho]=0.018m^2/s^2`
MATLAB CODE FOR BLOCKMESH FILE
clear all;
close all;
clc;
% Laminar flow simulation parameters
Re=2100;
rho=1000;
D=0.02;
R=D/2;
L_ent=0.06*Re*D;
mu=0.00089;
nu=mu/rho;
thetha=2;
%boundary conditions
v_avg= Re*nu/D;
v_max=2*v_avg;
L=L_ent+0.3;
deltap= (8*v_avg*mu*L)/R^2;
% Vertices co-ordinates in the local cordinate system
v0=[0 0 0];
v1=[0 R*cosd(thetha/2) -R*sind(thetha/2)];
v2=[0 R*cosd(thetha/2) R*sind(thetha/2)];
v3=[L 0 0];
v4=[L R*cosd(thetha/2) -R*sind(thetha/2)];
v5=[L R*cosd(thetha/2) R*sind(thetha/2)];
grading=0.2;
b1=50;
b2=1;
b3=200;
%creating blockmesh dict file for mesh generation of the wedge
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 %10s\n','version',' 2.0;');
fprintf(f1,'%16s %12s\n','format',' ascii;');
fprintf(f1,'%16s %17s\n','class','dictionary;');
fprintf(f1,'%16s %20s\n','object','blockMeshDict;');
fprintf(f1,'%2s\n','}');
fprintf(f1,'%2s\n\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'%2s\n\n','converToMeters 1;');
fprintf(f1,'vertices \n (\n');
fprintf(f1, '\t(%d %d %d)\n',v0);
fprintf(f1, '\t(%d %d %d)\n',v1);
fprintf(f1, '\t(%d %d %d)\n',v2);
fprintf(f1, '\t(%d %d %d)\n',v3);
fprintf(f1, '\t(%d %d %d)\n',v4);
fprintf(f1, '\t(%d %d %d)\n',v5);
fprintf(f1,' %2s\n\n', ');')
fprintf(f1,'blocks \n (\n');
fprintf(f1,' \t %s (%d %d %d) %s (%d 1 1) \n','hex (0 1 2 0 3 4 5 3)',b1,b2,b3,'simpleGrading',grading);
fprintf(f1,' %2s\n\n', ');');
fprintf(f1,'edges \n (\n');
fprintf(f1,'\t %s (0 %d 0)\n','arc 1 2',R);
fprintf(f1,'\t %s (%d %d 0)\n','arc 4 5',L,R);
fprintf(f1,' %2s\n\n', ');')
fprintf(f1,'boundary \n (\n');
fprintf(f1,' \t %5s \n','axis','{');
fprintf(f1,'\t %13s %s\n',' type','empty;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 3 3 0)',');','}');
fprintf(f1,' \t %5s \n','wall','{');
fprintf(f1,'\t %13s %s\n',' type','wall;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(2 5 4 1)',');','}');
fprintf(f1,' \t %5s \n','inlet','{');
fprintf(f1,'\t %13s %s\n',' type','patch;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 0 2 1)',');','}');
fprintf(f1,' \t %5s \n','outlet','{');
fprintf(f1,'\t %13s %s\n',' type','patch;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(3 4 5 3)',');','}');
fprintf(f1,' \t %5s \n','wedgeFront','{');
fprintf(f1,'\t %13s %s\n',' type','wedge;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 3 5 2)',');','}');
fprintf(f1,' \t %5s \n','wedgeBack','{');
fprintf(f1,'\t %13s %s\n',' type','wedge;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n %s\n\n','(0 1 4 3)',');','}',');');
fprintf(f1,'mergePatchPairs \n (\n );\n');
fprintf(f1,'%2s\n\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fclose(f1);Now the main objective of the matlab program is to generate the blockMeshDict file for creating the mesh for the given wedge angle and mesh grading factor. The first few lines of the code defines the parameters and the equations which are being used for the simulation. Next we define the 5 vertices to create the wegde for simulating the flow inside the wedge. After that we are using a technique known as file parsing to create the text file for the blockmesh file and then it is transfered to the systems folder. First we define the file name by usinf fopen command and we apply the fprintf statement to print all the text and define the boundaries in the fprinf command to define these wedge and the arc. Then we use fclose statement to close the file and to end the file parsing technique.
BLOCKMESHDICT FILE
/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
converToMeters 1;
vertices
(
(0 0 0)
(0 0.00999848 -0.000174524)
(0 0.00999848 0.000174524)
(2.82 0 0)
(2.82 0.00999848 -0.000174524)
(2.82 0.00999848 0.000174524)
);
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
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 0 2 1)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
wedgeFront
{
type wedge;
faces
(
(0 3 5 2)
);
}
wedgeBack
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
The blockmeshdict file was created for the wedge angle of 2 degrees from the matlab program. Now we can see that the blockmeshDict file is automatically created by using the file parsing technique. When we run the matlab program the text file is automatically generates the blockmesh file for the given wedge angle and mesh grading factor. Here we can see that for the wedge is formed by 5 vertices as given in the matlab code and the boundary field has been generated by the matlab code
CONTROLDICT 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 dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10;
deltaT 0.001;
writeControl timeStep;
writeInterval 200;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
The controldict file shows the that the start time is 0 seconds and the simulation will run for 10 seconds with the time step of 0.001 seconds.
BOUNDARY CONDITIONS
P 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
{
axis
{
type empty;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.018;
}
wall
{
type zeroGradient;
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
The p file gives us the pressure conditions for the geometry where we can see that at the inlet the pressure is zero gradient. The zero gradient indicates that it is Neumann boundary condition at the outlet. The Neumann BC computes the information from the inner nodes using various interpolation schemes from the outer nodes. Here we are supplying 0.018 metre as the boundary condition at the outlet as Openfoam deals with Kinematic pressure.
U 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.09345 0 0);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (0.09345 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
The U file gives the velocity conditions for the simulation. at the inlet we can see that the flow is defined which is the Dirichlet Boundary condition and at the outlet we are providing Neumann boundary condition. The no slip wall conditions are applied to top and bottom walls with front and back as empty walls. At the inlet we are saying that the flow is coming at the rate of 0.09345 m/s.
TRANSPORT PROPERTIES 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 dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 1e-6;
This file shows the kinematic viscosity information for the simulation. The kinematic viscosity is given as 1e-6 m^2/s
RESULTS
Velocity Distribution across the pipe
Velocity distribution shows the development of the hydrodynamic boundary layer as we travel along the length of the pipe. At the inlet there is step formation of boundary layers which shows that the flow is developing and then towards the outlet of the pipe we see there is a constant thickness in the boundary layer. This implies that the flow has developed fully after the hydrodynamic length of the pipe and it also gives us how the velocity is distributed along the length and radius of the cylinder.
VELOCITY PROFILE FOR DIFFERENT LENGTHS OF PIPE
The velocity profile is shown in half section as we are only considered the wedge as the entire cylinder. As for the comparison between the analytical and the simulation results, there seem to be a good match for the velocity profiles which are shown in the plots above. At the inlet there is formation of boundary layer, so the first plot captures the effect and as we go along the length of the pipe the velocity profile forms the non linear variation as it can be observed in the analytical solution. When the flow fully develops we can see that the curve exactly fits the profile shown in the theory of laminar flow across the pipe. it can be inferred that after a certain distance the parabolic curve is formed during the flow development and remains unchanged.
PRESSURE PROFILE
SHEAR STRESS
Note that when we are at the entrance of the pipe the shear stress profile there is no stress in the lower layers of the fluid but as we climb up towards the end of the pipe at the mid section we can see that there is a sudden non linear increase in the value of the shear stress which is shown in the first plot of the shear stress. In the second plot towards the length of boundary layer the shear stress almost has a linear variation across the radius of the wedge. But when the flow fully develops the shear stress plot looks like step formation in the plot which agrees with the Hagen Poiseuille equation profile. the shear stress gradient suddenly drops during the exit flow and it sharply increases again at the exit of the pipe.
COMPARISION BETWEEN ANALYTICAL AND NUMERICAL SOLUTION
Velocity Comparision
| Wedge Angle | maximum velocity in simulation(m/s) | maximum velocity in HP equation(m/s) | Percentage Error |
| 2 | 0.18808 | 0.1896 | 0.152 |
| 3 | 0.18911 | 0.1896 | 0.049 |
Now in the table there is a good match in the results between the analytical v/s numerical simulation for the velocity and the same can be seen for the kinematic pressure as well as the error is very small.
Kinematic Pressure Comparision
| Wedge Angle | kinematic pressure in simulation(m^2/s^2) | Kinematic pressure in HP equation(m^2/s^2) | Percentage Error |
| 2 | 0.01875 | 0.018 | 0.075 |
| 3 | 0.01884 | 0.018 | 0.0864 |
CONCLUSION
- Numerical solution agrees with Hagen Poiseuille law
- The maximum velocity has good match with the analytical and the numerical simulation as the error is very small
- For various wedge angles and mesh grading factor the error gradually increases with increase in wedge length but there is also decrease in error as we give finer mesh.