Week 11 - Simulation of Flow through a pipe in OpenFoam
Aim:-Simulation of Flow through a pipe in OpenFoam
Abstract
The flow through a pipe is the case of internal flow. The surface is the duct in the case of air, it is a pipe in the case of water. The flow is driven because of the pressure difference created by the pumps, air handling units, fans, etc. They find their major application refrigeration and air conditioning industries.
The fluid velocity changes from zero to wall and becomes maximum at the center. In fluid flow, it is preferred to take average velocity because it remains constant when the flow is incompressible and the cross-sectional area is constant.
Assumptions
- The steady laminar incompressible flow of a fluid is taken.
- There is no velocity in the radial direction.
- There is no acceleration since the flow is steady.
- Reynolds number for the laminar flow is 2100.
- Only the wedge shape of the pipe is considered because the flow is axis- symmetric.
- The equation governing the fluid flow is described by the Hagen-Poiseuille equation.
Hagen-Poiseuille equation

The assumptions of the equations are that the fluid is incompressible and Newtonian, the flow is laminar through a pipe of the constant cross-section of length that is longer than the diameter, and there is no acceleration of the fluid in the pipe.
The velocity profile remains unchanged in the hydrodynamically fully developed region, the shear stress at the wall also remains unchanged.
Hydrodynamically fully developed region `(∂u(r,x))/(∂x)=0` which means u =u(r)

Matlab coding and script for blockMeshDict file
% Simulation of Flow through a pipe
clear all
close all
clc
% Providing intial values for pipe simulation
%Laminar Flow
Re = 2100;
% Density
rho = 1000;
%Diameter of Pipe
D = 0.02;
%Radious of pipe
R = D/2;
%Wedge angle
theta = 3 ;
%Length of the pipe
L = 0.06*Re*D;
% Total Length of the pipe
L_p = L + 0.3;
% Dynamic viscosity
mu = 0.00089;
nu = mu/rho; %kinematic viscosity
% Hagen Poiseuille flow equations
v_avg = (Re*mu)/(rho*D);
v_max = 2*(v_avg);
del_p = (8*mu*L_p*v_avg)/(R^2);
kinematic_p = del_p/rho;
%vertices co-ordinate
v0 = [0 0 0];
v1 = [0 R*cosd(theta/2) -R*sind(theta/2)]; %point 1
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)];
grading = 0.2;
b1=50;
b2=1;
b3=200;
% creating blockmesh dict file for openfoam
h1 ='/*--------------------------------*- C++ -*----------------------------------*\';
h2 =' ========= |';
h3 =' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
h4 =' \\ / O peration | Website: https://openfoam.org';
h5 =' \\ / A nd | Version: 9';
h6 =' \\/ M anipulation |';
h7 ='\*---------------------------------------------------------------------------*/';
f1 = fopen('BLOCKMESHDict.txt','w');
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,'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','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'%2s\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 hex (0 1 2 0 3 4 5) (%d %d %d)',50,1,200);
fprintf(f1,' SimpleGrading (1 0.1 1)\n\n');
fprintf(f1,');\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','(0 1 4 3)',');','}');
fprintf(f1,'mergePatchPairs \n(\n);\n');
fprintf(f1,'%2s\n\n','// ************************************************************************* //');
fclose(f1);
Creating geometry of wedge profile
/*--------------------------------*- 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.996573e-03 -2.617695e-04)
(0 9.996573e-03 2.617695e-04)
(2.820000e+00 0 0)
(2.820000e+00 9.996573e-03 -2.617695e-04)
(2.820000e+00 9.996573e-03 2.617695e-04)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (200 50 1) SimpleGrading (1 2.000000e-01 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (2.820000e+00 1.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
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
(
);
//************************************************************************* //
Control Dict
/*--------------------------------*- 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 50;
deltaT 0.01;
writeControl timeStep;
writeInterval 200;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Transport Properties
/*--------------------------------*- 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.9e-7;
// ************************************************************************* //
Pressure
/*--------------------------------*- 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;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.018;
}
wall
{
type zeroGradient;
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //
Velocity
/*--------------------------------*- 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;
}
inlet
{
type fixedValue;
value uniform (0.0937312 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //
Geometry Created by blockMeshDict file

Velocity distribution of fully developed profile


Pressure distribution of fully developed profile


Simulation graph
Velocity profile at inlet

Velocity at distance x=0.5m

Velocity at distance x=1 m

Velocity at distance x=2 m

Pressure Profile

Shear stress:
- When the fluid is in motion, shear stresses are developed due to the particles in the fluid moving relative to one another.
- For a fluid flowing in a pipe Fluid velocity will be zero at the pipe wall. Velocity will increase while moving towards the center of the pipe. Shear Force is normally present because the adjacent layer of the fluid move with different velocities compared to each other.
- Fluids Obeying Newton’s law where the values of “mu are constant are known as Newtonian fluids. If ‘mu’ is constant, the shear stress is linearly dependent on the velocity gradient.
By considering the velocity of this relative motion, the Shear rate can be calculated shear rate is defined as a measurement of the extent or rate of relative motion between the adjacent layers of the moving fluid. The shear rate for fluid flowing between two parallel plates one moving at a constant speed and one is stationary.

Where.`μ/y` Is velocity gradient.
Shear stress is proportional to the shear rate where viscosity is the proportionality coefficient. This is known as Newton’s law of viscosity

Calculating Shear Stress
% Simulation of Flow through a pipe
clear all
close all
clc
% Providing intial values for pipe simulation
%Laminar Flow
Re = 2100;
% Density
rho = 1000;
%Diameter of Pipe
D = 0.02;
%Radious of pipe
R = D/2;
%Wedge angle
theta = 3 ;
%Length of the pipe
L = 0.06*Re*D;
% Total Length of the pipe
L_p = L + 0.3;
% Dynamic viscosity
mu = 0.00089;
nu = mu/rho; %kinematic viscosity
% Hagen Poiseuille flow equations
v_avg = (Re*mu)/(rho*D);
v_max = 2*(v_avg);
del_p = (8*mu*L_p*v_avg)/(R^2);
kinematic_p = del_p/rho;
% calculating shear stress
r = linspace(-0.01,0.01,100);
for i=1:length(r)
tau(i)=(2*mu*v_max*abs(r(i)))/R^2;
end
plot(r,tau)
grid on
xlabel('Radius of pipe in meter')
ylabel('Shear stress in pascal')
title('Shear stress along the radius of pipe')
legend('shear stress of the profile')
Shear stress graph:-

not equal to zero because it is calculated at 1m length in which the flow is not fully developed so it is not equal to zero.