Week 11 - Simulation of Flow through a pipe in OpenFoam
Our objective is to simulate a flow through a pipe. The governing equation for this will be Hagen-Poiseullie equation, and the assumptions made are ,
- Steady Incompressible Laminar Flow
- 2-Dimentinal flow
- Fluid properties are constant along the length of the pipe
- Axisymetric flow
Then we computing the results of flow through pipe and comparing with thereotical value.
We assume the length of the pipe to be 2 m and the diameter to be 15mm. Instead of simulating for the entire pipe, we simulate of a wedge of paticular angle as it is axisymetric flow. We choose water as fluid and assume its renyolds number Re as 2100 and density p as 1000 kg/m^3. The hydrodynamic boundary layed occurs at the walls of the pipe due to no slip condition,that is the fluid near the walls will have zero velocity. And the velocith of fluid along the axis will be high at any given point on the axis.And velocity both above and below the axis will gradually decrease, this is due to the viscous property of the fluid which produces a friction force between the adjucent fluid layers that resist its flow. Therefore a velocity profile is created as below.

In the above image we can clearly visuvalise the entrance length(Le), it is the distance between the inlet and the point wehre fully developed flow occurs. And the average velocity remains conatant along the axis of the pipe. for laminar flow the entrance length is given by,
`Le approx 0.06 * Re*D`
And similarly we compute ,
Maximum velocity ,`Vmax =`2 * `Vavg`
Vmax = 0.187 `ms^-1`
Pressure drop , `Delta P =( 32 *mu *V*L)/ (D^2)`
`Delta P =`34.53 pa
Kinematic pressure drop = `(Delta P)/ (1000)`
= 0.03453
Procedure:
- We first choose the wedge angle as 3 .
- Then we write a code in matkab for auto generating the blockmeshDic file ,to create a mesh
- We first initialise the constants and then the formula to calculate the velocity,pressure ,entery length.
- With the mentioned theta we write a formula to determine the remaining points for a vertices.
- We define the edges,boundary and faces for generating the mesh on a text file blockmesh.
- Once we finish creating the blockmesh we choose the icoFoam solver to solve our problem.\
- Then with the solved values for the constants we apply this inlet velocity and pressure in the initial conditions file.
- We then simulate for time of 10 sec and we obtain results for every 0.1 seconds.
- Also we have to mention the Kinematic viscocity at 25 degree celcius and write it into Transport Properties.
- After generating mesh we open paraview and visuvalise the mesh.
- Then we simulate for both velocity and pressure and obtain the plot at different points includind for the entrance length (Le) and find the maximum velocity to compare then we plot for the shear stress graph.
clear all
close all
clc
Re = 2100;
p= 1000;
m=0.00089;
D=0.015;
R=D/2;
L_e=(0.06*Re*D);
L= L_e+0.3;
n=m/p;
v_avg=Re*(n/D);
vmax=2*v_avg;
theta=3;
Pdiff= (32*m*L*v_avg)/(D^2)
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)];
h1='/*--------------------------------*- C++ -*----------------------------------*\';
h2='| ========= | |';
h3='| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |';
h4='| \\ / O peration | Version: 2.2.1 |';
h5='| \\ / A nd | Web: www.OpenFOAM.com |';
h6='| \\/ M anipulation | |';
h7='\*---------------------------------------------------------------------------*/';
h8='FoamFile ';
h9=' {';
h10=' version 2.0;';
h11=' format ascii;';
h12=' class dictionary;';
h13=' object blockMeshDict;';
h14=' }';
h15=' // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
h16='convertToMeters 1;';
h17='mergePatchPairs';
h18='(';
h19=');';
h20='// ************************************************************************* //';
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',h7);
fprintf(f1,'%s \n',h8);
fprintf(f1,'%s \n',h9);
fprintf(f1,'%s \n',h10);
fprintf(f1,'%s \n',h11);
fprintf(f1,'%s \n',h12);
fprintf(f1,'%s \n',h13);
fprintf(f1,'%s \n',h14);
fprintf(f1,'%s \n',h15);
fprintf(f1,'%s \n',h16);
fprintf(f1,'vertices \n ( \n ');
fprintf(f1,' (%d %d %d)\n ',v0);
fprintf(f1,' (%d %d %d)\n ',v1);
fprintf(f1,' (%d %d %d)\n ',v2);
fprintf(f1,' (%d %d %d)\n ',v3);
fprintf(f1,' (%d %d %d)\n ',v4);
fprintf(f1,' (%d %d %d)\n ',v5);
fprintf(f1,'\n %2s\n',');');
fprintf(f1,'\n\n');
fprintf(f1,' blocks \n ( ');
fprintf(f1,'\n\t %s (%d %d %d) %s (%d 1 1)','hex(0 1 2 0 3 4 5 3',50,1,200,'Simplegrading',0.2);
fprintf(f1,'\n %2s \n' ,');');
fprintf(f1,'\n\n');
fprintf(f1,'edge \n ( \n ');
fprintf(f1,'\n %3s (0 %d 0) \n','arc 1 2',R);
fprintf(f1,'\n %3s (%d %d 0) \n','arc 4 5',L,R);
fprintf(f1,'\n %2s \n' ,');');
fprintf(f1,'\n\n');
fprintf(f1,'boundary \n ( \n ');
fprintf(f1,' \n %5s \n','axis');
fprintf(f1,' \n %5s \n','{' );
fprintf(f1,'\n\t %5s %5s \n','type','empty;');
fprintf(f1,' \n %18s \n','faces');
fprintf(f1,' \t %14s','(');
fprintf(f1,'\n\t\t %20s \n','(0 3 3 0)');
fprintf(f1,'\t %14s \n' ,');');
fprintf(f1,' \n %5s \n','}' );
fprintf(f1,' \n %5s \n','Inlet');
fprintf(f1,' \n %5s \n','{' );
fprintf(f1,'\n\t %5s %5s \n','type','patch;');
fprintf(f1,' \n %18s \n','faces');
fprintf(f1,' \t %14s' ,'(');
fprintf(f1,'\n\t\t %20s \n','(0 1 2 0)');
fprintf(f1,'\t %14s \n' ,');');
fprintf(f1,' \n %5s \n','}' );
fprintf(f1,' \n %5s \n','Wall');
fprintf(f1,' \n %5s \n','{' );
fprintf(f1,'\n\t %5s %5s \n','type','wall;');
fprintf(f1,' \n %18s \n','faces');
fprintf(f1,' \t %14s' ,'(');
fprintf(f1,'\n\t\t %20s \n','(2 1 4 5)');
fprintf(f1,'\t %14s \n' ,');');
fprintf(f1,' \n %5s \n','}' );
fprintf(f1,' \n %5s \n','Back');
fprintf(f1,' \n %5s \n','{' );
fprintf(f1,'\n\t %5s %5s \n','type','wedge;');
fprintf(f1,' \n %18s \n','faces');
fprintf(f1,' \t %14s' ,'(');
fprintf(f1,'\n\t\t %20s \n','(0 2 5 3)');
fprintf(f1,'\t %14s \n' ,');');
fprintf(f1,' \n %5s \n','}' );
fprintf(f1,' \n %5s \n','Front');
fprintf(f1,' \n %5s \n','{' );
fprintf(f1,'\n\t %5s %5s \n','type','wedge;');
fprintf(f1,' \n %18s \n','faces');
fprintf(f1,' \t %14s' ,'(');
fprintf(f1,'\n\t\t %20s \n','(0 3 4 1)');
fprintf(f1,'\t %14s \n' ,');');
fprintf(f1,' \n %5s \n','}' );
fprintf(f1,' \n %5s \n','Outlet');
fprintf(f1,' \n %5s \n','{' );
fprintf(f1,'\n\t %5s %5s \n','type','patch;');
fprintf(f1,' \n %18s \n','faces');
fprintf(f1,' \t %14s' ,'(');
fprintf(f1,'\n\t\t %20s \n','(3 5 4 3)');
fprintf(f1,'\t %14s \n' ,');');
fprintf(f1,' \n %5s \n','}' );
fprintf(f1,'\n %2s \n' ,');');
fprintf(f1,'\n\n');
fprintf(f1,'%s \n',h17);
fprintf(f1,'%s \n',h18);
fprintf(f1,'%s \n',h19);
fprintf(f1,'%s \n',h20);
fclose(f1);
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.2.1 |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 4.998287e-02 -1.308847e-03)
(0 4.998287e-02 1.308847e-03)
(3.003000e-01 0 0)
(3.003000e-01 4.998287e-02 -1.308847e-03)
(3.003000e-01 4.998287e-02 1.308847e-03)
);
blocks
(
hex(0 1 2 0 3 4 5 3 (50 1 200) Simplegrading (2.000000e-01 1 1)
);
edge
(
arc 1 2 (0 5.000000e-02 0)
arc 4 5 (3.003000e-01 5.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
Inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
Wall
{
type wall;
faces
(
(2 1 4 5)
);
}
Back
{
type wedge;
faces
(
(0 2 5 3)
);
}
Front
{
type wedge;
faces
(
(0 3 4 1)
);
}
Outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10;
deltaT 0.001;
writeControl timeStep;
writeInterval 100;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.0937312 0 0);
boundaryField
{
axis
{
type empty;
}
fixedWalls
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type wedge;
}
Back
{
type wedge;
}
Inlet
{
type fixedValue;
value uniform (0.0937312 0 0);
}
Outlet
{
type zeroGradient;
}
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
axis
{
type empty;
}
fixedWalls
{
type zeroGradient;
}
front
{
type wedge;
}
Back
{
type wedge;
}
Inlet
{
type zeroGradient;
}
Outlet
{
type fixedValue;
value uniform 0.018;
}
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.91e-7;
// ************************************************************************* //
Results:
Mesh:


Velocity Distribution:
Velocity at inlet and at the exit


Pressure:

Plots:
At x=0.02

At x= 0.05

At x=0.5

At x= 1.89
At the entery length Le

At x= 2.19
At the end point

Pressure distribution along X:

Shear Stress:

Conclusion:
As we expected the velocity graph is parabollic and the maximum velocity throuch simulation is 0.162 `ms^-1` . The pressure plot is linear as we gave a small amout of inlet pressure. the poin 1.89 in where th fully developed profile is obtained. The analytical pressure drop is 0.03453 Nm^2/s and the simulation pressure is 0.031.