Week 11 - Simulation of Flow through a pipe in OpenFoam
In this problem we considering an incompressible laminar flow of Reynold's number 2100 through a pipe of diameter 0.02 m and length of 2.3 m. Here water is the flowing fluid with density 1000`kgm^-3` and a dynamic viscosity of 0.00089 Pa.s.
Since the time taken for completing the simulation is much higher for processing the whole pipe, we take a small sector of the pipe for the easiness. Here we taken the wedge angle as `3^0`.
- Hagen Poiseuille's equation
With the help of Hagen Poiseuille's equation we can find out the pressure loss of fully developed flow.
`DeltaP=(32muLV_(avg))/D^2`
Where `DeltaP`=Pressure Drop
`mu`=Dynamic viscosity
`L`=Total length of the pipe-Hydrodynamic entrance length
`V_(avg)`=Average velocity
D=Diameter of the pipe
- Hydrodynamic Entrance Length
The region from the pipe inlet to the point at which the velocity profile is fully developed is called the hydrodynamic entrance region, and the length of this region is called the hydrodynamic entry length `L_h`.

From the formula for Reynold's number, we can find the average velocity,
`Re=(rhov_(avg)D)/mu`
Where,
`rho`=Density of water
So,
- `v_(avg)=(Remu)/(rhoD)=(2100xx0.00089)/(1000xx0.02)=0.09345ms^-1`
Since it is a laminar flow w.k.t,
`L_h/D~~0.05Re`
So,
- Hydrodynamic entrance length, `L_h=2.1m`
- Maximum velocity, `u_max=2v_(avg)=0.1869ms^-1`
Pressure drop can be find out by Hagen Poiseuille's equation.
- `DeltaP=1.3307`Pa
Shear stress, `tau=(2muu_(max))/r_h`
where, `r_h`=distance between centre of the pipe to the wall at full developed flow region a.k.a radius at hydrodynamic entrance region.
here r_h=0.02 m
So,
- `tau=0.03327` Pa
GENERAL STEPS
- Create 'blockMeshDict' file using matlab.
- Make necessary changes in other files too. (eg:- Files like controlDict, fvSchemes, fvSolution,transportProperties,pressure file[p] and velocity file[u] )
- Run the simulation by openFOAM by icoFOAM solver. ( Because the flow is incompressible )
- Plot the graphs.
- Compare the numerical and analytical values
- MATLAB CODE FOR CREATING 'blockMeshDict' FILE
clear all
close all
clc
%FINDING HYDRADYNAMIC ENTRANCE LENGTH,MAXIMUM VELOCITY,PRESSURE DROP and
%SHEAR STRESS
%Assuming diameter of the pipe,
D=0.02;
%Reynolds number
Re=2100;
%Density of water
rho=1000;
%Dynamic Viscosity
mu=0.00089;
%radius of pipe
r=0.01;
%distance from centre of the pipe to the wall
r_h=0.01;
%wedge angle
theta=3;
%Length of the pipe
L=2.3;
%Average velocity,
v_avg = (Re*mu)/(rho*D);
% 1.HYDRODYNAMIC ENTRANCE LENGTH
L_h = 0.05*Re*D; % L = L_h+0.2
%Total length of pipe - hydrodynamic entrance length
L_p = L-L_h;
% 2.MAXIMUM VELOCITY
u_max = 2*v_avg;
% 3.PRESSURE DROP
p_d = (32*mu*L_p*v_avg)/(D^2);
% 4.SHEAR STRESS
t_w = (2*mu*u_max*r)/(r_h^2);
%CREATING blockMeshDict File
%creating file
f1=fopen('blockMeshDict.txt','w');
%opening part
h1='/*--------------------------------*- C++ -*----------------------------------*\';
h2=' ========= |';
h3=' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
h4=' \\ / O peration | Website: https://openfoam.org';
h5=' \\ / A nd | Version: 8';
h6=' \\/ M anipulation |';
h7='\*---------------------------------------------------------------------------*/';
h8='FoamFile';
h9='{';
h10=' version 2.0;';
h11=' format ascii;';
h12=' class dictionary;';
h13=' object blockMeshDict;';
h14='}';
h15='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
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,'\n')
h16='convertToMeters 1;';
fprintf(f1,'%s\n',h16)
fprintf(f1,'\n')
%writing vertices
h17='vertices';
fprintf(f1,'%s\n',h17)
h18='(';
fprintf(f1,'%s\n',h18)
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)];
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)
h19=');';
fprintf(f1,'%s\n',h19)
fprintf(f1,'\n')
%writing blocks
h20='blocks';
h21=' hex (0 3 5 2 0 3 4 1) (200 10 1) simpleGrading (1 0.1 1)';
h22=');';
fprintf(f1,'%s\n',h20)
fprintf(f1,'(\n')
fprintf(f1,'%s\n',h21)
fprintf(f1,'%s\n',h22)
fprintf(f1,'\n')
%writing edges
h23='edges';
h24='(';
h25=' arc 2 1';
h26=' arc 5 4';
h27=');';
fprintf(f1,'%s\n',h23)
fprintf(f1,'%s\n',h24)
fprintf(f1,'%s (%d %d %d)\n',h25,0,r,0)
fprintf(f1,'%s (%d %d %d)\n',h26,L,r,0)
fprintf(f1,'%s\n',h27)
fprintf(f1,'\n')
%writing boundary
h28='boundary';
h29='(';
%inlet
h30=' Inlet';
h31=' {';
h32=' type patch;';
h33=' faces';
h34=' (';
h35=' (0 1 2 0)';
fprintf(f1,'%s\n',h28)
fprintf(f1,'%s\n',h29)
fprintf(f1,'%s\n',h30)
fprintf(f1,'%s\n',h31)
fprintf(f1,'%s\n',h32)
fprintf(f1,'%s\n',h33)
fprintf(f1,'%s\n',h34)
fprintf(f1,'%s\n',h35)
h36=' );';
h37=' }';
fprintf(f1,'%s\n',h36)
fprintf(f1,'%s\n',h37)
%outlet
h38=' Outlet';
h39=' {';
h40=' type patch;';
h41=' faces';
h42=' (';
h43=' (3 5 4 3)';
h44=' );';
h45=' }';
fprintf(f1,'%s\n',h38)
fprintf(f1,'%s\n',h39)
fprintf(f1,'%s\n',h40)
fprintf(f1,'%s\n',h41)
fprintf(f1,'%s\n',h42)
fprintf(f1,'%s\n',h43)
fprintf(f1,'%s\n',h44)
fprintf(f1,'%s\n',h45)
%surface
h46=' surface';
h47=' {';
h48=' type wall;';
h49=' faces';
h50=' (';
h51=' (2 1 4 5)';
h52=' );';
h53=' }';
fprintf(f1,'%s\n',h46)
fprintf(f1,'%s\n',h47)
fprintf(f1,'%s\n',h48)
fprintf(f1,'%s\n',h49)
fprintf(f1,'%s\n',h50)
fprintf(f1,'%s\n',h51)
fprintf(f1,'%s\n',h52)
fprintf(f1,'%s\n',h53)
%axis
h54=' axis';
h55=' {';
h56=' type empty;';
h57=' faces';
h58=' (';
h59=' (0 3 3 0)';
h60=' );';
h61=' }';
fprintf(f1,'%s\n',h54)
fprintf(f1,'%s\n',h55)
fprintf(f1,'%s\n',h56)
fprintf(f1,'%s\n',h57)
fprintf(f1,'%s\n',h58)
fprintf(f1,'%s\n',h59)
fprintf(f1,'%s\n',h60)
fprintf(f1,'%s\n',h61)
%front
h62=' front';
h63=' {';
h64=' type wedge;';
h65=' faces';
h66=' (';
h67=' (0 3 4 1)';
h68=' );';
h69=' }';
fprintf(f1,'%s\n',h62)
fprintf(f1,'%s\n',h63)
fprintf(f1,'%s\n',h64)
fprintf(f1,'%s\n',h65)
fprintf(f1,'%s\n',h66)
fprintf(f1,'%s\n',h67)
fprintf(f1,'%s\n',h68)
fprintf(f1,'%s\n',h69)
%back
h70=' back';
h71=' {';
h72=' type wedge;';
h73=' faces';
h74=' (';
h75=' (0 2 5 3)';
h76=' );';
h77=' }';
fprintf(f1,'%s\n',h70)
fprintf(f1,'%s\n',h71)
fprintf(f1,'%s\n',h72)
fprintf(f1,'%s\n',h73)
fprintf(f1,'%s\n',h74)
fprintf(f1,'%s\n',h75)
fprintf(f1,'%s\n',h76)
fprintf(f1,'%s\n',h77)
fprintf(f1,')\n')
fprintf(f1,'\n')
%Final part
h78='mergePatchPairs';
h79='(';
h80=');';
h81='// ************************************************************************* //';
fprintf(f1,'%s\n',h78)
fprintf(f1,'%s\n',h79)
fprintf(f1,'%s\n',h80)
fprintf(f1,'\n')
fprintf(f1,'%s\n',h81)
- OpenFOAM FILES
1.blockMeshDict
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.996573e-03 2.617695e-04)
(0 9.996573e-03 -2.617695e-04)
(2.300000e+00 0 0)
(2.300000e+00 9.996573e-03 2.617695e-04)
(2.300000e+00 9.996573e-03 -2.617695e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 10 1) simpleGrading (1 0.1 1)
);
edges
(
arc 2 1 (0 1.000000e-02 0)
arc 5 4 (2.300000e+00 1.000000e-02 0)
);
boundary
(
Inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
Outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
surface
{
type wall;
faces
(
(2 1 4 5)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
2.controlDict
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 25;
deltaT 0.01;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
3.fvSolution
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0.05;
}
pFinal
{
$p;
relTol 0;
}
U
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-05;
relTol 0;
}
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
// ************************************************************************* //
4.fvSchemes
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linear;
}
laplacianSchemes
{
default Gauss linear orthogonal;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default orthogonal;
}
// ************************************************************************* //
5.transportProperties
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9e-07;
// ************************************************************************* //
6.p
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ 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 0;
}
surface
{
type zeroGradient;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
7.U
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ 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 (0.0935 0 0);
}
Outlet
{
type zeroGradient;
}
surface
{
type noSlip;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
- OpenFOAM PROCEDURE
- Open the terminal window. (ctrl+alt+T)
- Go to tutorial folder by typing 'tut'
- Go to Incompressible folder by typing 'cd Incompressible/'
- Go to icoFoam folder by typing 'cd icoFoam/'
- Inside that, go to cavity folder by typing 'cd cavity/'
- Copy the cavity folder inside the main folder, to run folder and renamed as 'pipe'. So that the offical file of the OpenFOAM remain untouched.
- Copy the new 'blockMeshDict' file that we made using MATLAB to the system folder inside the 'pipe' folder. Copy and replace it.
- Make necessary changes to other files inside the 0,system,constant folders. (As mentioned above)
- Go to the run folder by typing 'run'
- Open the 'pipe' folder. (cd pipe)
- After completing whole editing use 'blockMesh' command to generate the mesh and 'checkMesh' to verify it.
- Then type 'icoFoam' run the simulation.
- Then open paraView by typing 'paraFoam'.
- When it gets open click 'apply' to show the geometry.
- Plot the velocity profile and Shear stress profile by using 'Plot Over Line' option.
- Edit the point 1 and point 2 values to 0.085, then only we able to visualise the velocity profile at 0.085.
- Take screenshots.
RESULTS
- ISOMETRIC VIEW OF MESHED GEOMETRY

- VELOCITY PROFILE OF FULLY DEVELOPED FLOW

- SHEAR STRESS PROFILE AT X=2.1m FROM INLET

- SHEAR STRESS PROFILE AT X=2.3m FROM INLET

COMPARISON OF RESULTS
| Theoretical Results | Analytical Results | |
| Maximum Velocity (m/s) | 0.1869 | 0.1786 |
| Pressure Drop (Pa) | 1.3307 | 1.3515 |
| Maximum Shear stress (Pa) | 0.0333 | 0.0318 |
CONCLUSION
- The velocity profile we have obtained is parabolic in nature.
- From velocity profile it is clear that the flow velocity is maximum at centre of the pipe and nearly 0 at the walls.
- From the shear stress plots we can observe that the shear stress maximum at walls and less at the centre of the pipe.